\(Description\)

给定 \(n\) 个正整数序列 ,每个序列长度为 \(m\)。

选择至少 \(1\) 个序列,在每个被选择的序列中选择一个元素,求出所有被选择的元素的 \(\gcd\)。

求所有方案的结果之和,答案对 \(1e9+7\) 取模。两种方案不同,当且仅当存在至少一个元素,在一种方案中被选择,在另一种中没有。

\(Input\)

第一行,两个正整数 \(n,m\)。

接下来 \(n\) 行,每行 \(m\) 个正整数,第 \(i\) 行代表序列 。

\(Output\)

第一行,一个整数,代表答案对 \(1e9+7\) 取模的结果。

解析

一道比较难的莫比乌斯反演题,要用到其中的性质,且按套路行事技巧处很多

最后推出的是一个关于欧拉函数的式子,莫比乌斯不见了

好,现在进行套路推导

先设 \(f(x)\) 表示选择至少一个序列,在每个被选择的序列选择一个元素,它们的 \(\gcd = x\) 的方案数。

则易得

\[f(x)=\prod_{i=1}^n(\sum_{j=1}^m [a_{i,j}=x]+1)-1
\]

然后套路 \(F(x)\) 表示同 \(f(x)\) 但涵盖了 \(x\) 的倍数,即

\[F(x)=\sum_{x|d}f(d)=\sum_{x|d}(\prod_{i=1}^n(\sum_{j=1}^m [a_{i,j}=d]+1)-1)
\]

然后我们发现,我们先枚举 \(x\) ,再枚举其倍数 \(d\),而后面 \([a_{i,j}=d]\) 肯定是 \(x\) 的倍数,所以我们可以简化式子

\[\sum_{x|d}(\prod_{i=1}^n(\sum_{j=1}^m [a_{i,j}=d]+1)-1) = \prod_{i=1}(\sum_{j=1}^m [x|a_{i,j}]+1)-1
\]

而此时,为了日后式子的简便即实现,我们设

\(cnt_{i,x}=\sum_{j=1}^m[x|a_{i,j}]\) 表示第 \(i\) 个数列所有是 \(x\) 的倍数的数的个数

再为了枚举得到所有答案,我们设 \(lim\) 表示所有元素的最大值

然后一波推式子,反演

\[\begin{aligned}
\sum_{i=1}^{lim}if(i)
&=\sum_{i=1}^{lim}i\sum_{i|d}F(d)\mu(\frac{d}{i}) \\
&=\sum_{i=1}^{lim}i\sum_{i|d}\prod_{j=1}^n((cnt_{j,d}+1)-1)\mu(\frac{d}{i}) \\
&=\sum_{d=1}^{lim}\sum_{i|d}i\mu(\frac{d}{i})(\prod_{j=1}^n(cnt_{j,d}+1)-1)
\end{aligned}
\]

然后,然后~~~好像没戏了

但,我们有伟大的欧拉!!

\[\begin{aligned}
\varphi(n)
&=\sum_{i=1}^n[\gcd(i,n)=1] \\
&=\sum_{i=1}^n\sum_{d|gcd(i,j)}\mu(d) \\
&=\sum_{d|n}\mu(d)\frac{n}{d}
\end{aligned}
\]

哈哈哈,太棒了!

相同的一部分,代入式子

\[\begin{aligned}
\sum_{d=1}^{lim}\sum_{i|d}i\mu(\frac{d}{i})(\prod_{j=1}^n(cnt_{j,d}+1)-1)
&=\sum_{d=1}^{lim}\varphi(n)(\prod_{j=1}^n(cnt_{j,d}+1)-1)
\end{aligned}
\]

于是这题就这样了。

线性筛 \(\varphi\),预处理 \(cnt\) 数组(根据套路,不要枚举因子而是枚举倍数)。

\(Code\)

#include<cstdio>
#include<iostream>
using namespace std;
typedef long long LL; const int N = 1e5;
const LL mod = 1e9 + 7;
int lim , n , m , a[25][N + 5] , cnt[25][N + 5] , phi[N + 5] , prime[N + 5] , vis[N + 5] , tot;
LL ans; inline void getPhi()
{
phi[1] = 1;
for(register int i = 2; i <= N; i++)
{
if (!vis[i]) phi[prime[++tot] = i] = i - 1;
for(register int j = 1; j <= tot && prime[j] * i <= N; j++)
{
vis[prime[j] * i] = 1;
if (i % prime[j] == 0)
{
phi[prime[j] * i] = phi[i] * prime[j];
break;
}
phi[prime[j] * i] = phi[i] * (prime[j] - 1);
}
}
} inline void getCnt()
{
for(register int i = 1; i <= n; i++)
for(register int j = 1; j <= lim; j++)
for(register int k = 2; k * j <= lim; k++)
cnt[i][j] += cnt[i][j * k];
} int main()
{
freopen("b.in" , "r" , stdin);
freopen("b.out" , "w" , stdout);
scanf("%d%d" , &n , &m);
for(register int i = 1; i <= n; i++)
for(register int j = 1; j <= m; j++)
scanf("%d" , &a[i][j]) , cnt[i][a[i][j]]++ , lim = max(lim , a[i][j]);
getPhi() , getCnt();
for(register int d = 1; d <= lim; d++)
{
LL res = 1;
for(register int j = 1; j <= n; j++) res = res * (LL)(cnt[j][d] + 1) % mod;
ans = (ans + (LL)((LL)phi[d] * (res - 1)) % mod) % mod;
}
printf("%lld" , ans);
}

最新文章

  1. android Intent介绍
  2. JAVA编程思想(第四版)学习笔记----11.4 容器的打印
  3. Spring操作指南-AOP基本示例(基于XML)
  4. UNABLE TO PURGE A RECORD(二)
  5. js闭包问题
  6. 用apktool工具进行apk的编译和反编译
  7. 解决ecshop登陆自动退出的莫名现象
  8. C#导出数据至excel模板
  9. String.equals()
  10. 基于visual Studio2013解决面试题之1309求子集
  11. MulticastSocket绑定端口的问题
  12. SqlSession 同步为注册,因为同步未激活
  13. Fiddler抓包【1】_介绍及界面概述
  14. linux 修改配色
  15. JSP 指令 脚本元素 表达式 声明
  16. Unity应用架构设计(9)——构建统一的 Repository
  17. ppt提取文字
  18. poj3468
  19. javascript将list转换成树状结构
  20. spring学习 十六 spring加载属性文件

热门文章

  1. orcl between and 时间
  2. day25 前端
  3. 【数据库】Postgresql、PG的分区操作:创建、删除指定分区,非分区表转分区表
  4. go-dongle 0.2.0 版本发布了,一个轻量级、语义化的 golang 编码解码、加密解密库
  5. MongoDB - 数据模型的设计模式
  6. Qt操作Json小结
  7. Python怎么引入不同的库?
  8. 第五篇:前端之JQuery
  9. 封装一个python的pymysql操作类
  10. 德摩根定律的证明 De Morgan&#39;s law