P3704 [SDOI2017]数字表格
斯德哥尔摩
2018-05-13 19:31:16
[P3704 [SDOI2017]数字表格](https://www.luogu.org/problemnew/show/P3704)
题目要求：$Ans=\Pi_{i=1}^{n}\Pi_{j=1}^{m}f[gcd(i,j)]$
不妨设$n<=m$
则$Ans=\Pi_{d=1}^{n}\Pi_{i=1}^{n}\Pi_{j=1}^{m}f[d]\times [gcd(i,j)==d]$
直接把那啥$f(d)$提取出来就变成了：
$Ans=\Pi_{d=1}^{n}f[d]^{\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]}$
显然，那个指数用莫比乌斯+数论分块可以$O(\sqrt n)$做，外面的一层也可以，所以总复杂度为$O(n)$
但是显然不够，需要变形。
将指数单独拿出来看：
$\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]$
二话不说莫比乌斯反演：
$\sum_{i=1}^{n/d} \mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor$
设$T=id$，在整个式子里拎出来：
$\Pi_{T=1}^{n}\Pi_{d|T}f[d]^{\mu(T/d)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$
化简一下：
$\Pi_{T=1}^{n}(\Pi_{d|T}f[d]^{\mu(T/d)})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$
很明显，对那个啥$\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor$数论分块。
那——里面的东西怎么办？又不能线性筛。。。
等等——不能线性筛就暴力算呀！
数据范围$10^6$
每个数暴力算到他的倍数里面去
即：$\frac{n}{1}+\frac{n}{2}+\frac{n}{3}+...+\frac{n}{10^6}$
这玩意差不多$15n$，所以直接暴力把那个东西的前缀给求出来，就可以做到$O(\sqrt n)$了。
附代码：
```cpp
#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 1000010
#define MOD 1000000007
using namespace std;
long long f[MAXN],g[MAXN],sum[MAXN];//直接开long long
int k=0,prime[MAXN],mu[MAXN];
bool np[MAXN];
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
long long mexp(long long a,long long b,long long c){//快速幂
long long s=1;
while(b){
if(b&1)s=s*a%c;
a=a*a%c;
b>>=1;
}
return s;
}
void make(){//莫比乌斯专属线性筛
int m=MAXN-10;
f[1]=g[1]=sum[0]=sum[1]=mu[1]=1;
np[1]=true;
for(int i=2;i<=m;i++){
f[i]=(f[i-1]+f[i-2])%MOD;
g[i]=mexp(f[i],MOD-2,MOD);
sum[i]=1;
if(!np[i]){
prime[++k]=i;
mu[i]=-1;
}
for(int j=1;j<=k&&prime[j]*i<=m;j++){
np[prime[j]*i]=true;
if(i%prime[j]==0)break;
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=m;i++){//暴力求前缀
if(!mu[i])continue;
for(int j=i;j<=m;j+=i)sum[j]=sum[j]*(mu[i]==1?f[j/i]:g[j/i])%MOD;
}
for(int i=2;i<=m;i++)sum[i]=sum[i]*sum[i-1]%MOD;
}
void work(){//求解
int n,m;
long long inv,ans=1;
n=read();m=read();
if(n>m)swap(n,m);
for(int i=1,last;i<=n;i=last+1){
last=min(n/(n/i),m/(m/i));
inv=sum[last]*mexp(sum[i-1],MOD-2,MOD)%MOD;//求逆元
ans=ans*mexp(inv,1LL*(n/i)*(m/i)%(MOD-1),MOD)%MOD;
}
printf("%lld\n",ans);
}
int main(){
int t=read();
make();
while(t--)work();//多组询问
return 0;
}
```