bzoj 4816: [Sdoi2017]数字表格

题目大意:

题解.

\[\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j))\]
\[\prod_{x=1}^nf(x)^{\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j) = x]}\]
\[\prod_{x=1}^nf(x)^{\sum_{i=1}^{\frac{n}{x}}\sum_{j=1}^{\frac{m}{x}}[gcd(i,j) = 1]}\]
\[\prod_{x=1}^nf(x)^{\sum_{d=1}^{\frac{n}{x}}\mu(d)*[\frac{n}{dx}]*[\frac{m}{dx}]}\]
\[\prod_{T=1}^n(\prod_{d|T}f(d)^{\mu(\frac{T}{d})})^{[\frac{n}{T}]*[\frac{m}{T}]}\]

设 : \(g(T) = \prod_{d|T}f(d)^{\mu(\frac{T}{d})}\)

如果可以快速预处理\(g(x)\)就可以做到单次\(O(\sqrt{n}\log n)\)
那么考虑如何快速预处理\(g(T)\)
考虑到\(\mu(x)\)有不少是0,剩下的只有1/-1两种取值.
那么对于单个的\(g(T)\)的计算可以只枚举所有有用的取值.
在T的约数中只选取mu值不为0的值进行统计.
即相同的质因子最多包含一个.所以发现可以找出所有不相同的质因子.
显然个数少到可以状压起来 (35711131719 = 4849845 > 1000000)
所以状压一下每个质因子是否选取.相应的进行函数值的计算.

为了辅助快速找出所有不同的质因子
可以提前线性筛出所有数的最小质因子及除去这个最小因子后成为的数。
单点求出所有的\(g(x)\)后就可以直接分块处理每个询问了.

#include <cstdio>
#include <cstring>
#include <algorithm>
using namespace std;
typedef long long ll;
inline void read(int &x){
x=0;static char ch;static bool flag;flag = false;
while(ch=getchar(),ch<'!');if(ch == '-') ch=getchar(),flag = true;
while(x=(x<<1)+(x<<3)+ch-'0',ch=getchar(),ch>'!');if(flag) x=-x;
}
#define rg register int
#define rep(i,a,b) for(rg i=(a);i<=(b);++i)
#define per(i,a,b) for(rg i=(a);i>=(b);--i)
const int maxn = 1000010;
const int mod = (1e9 + 7);
int pri[maxn],cnt,mu[maxn],minp[maxn],oth[maxn];
bool vis[maxn];
inline ll qpow(ll x,ll p){
ll ret = 1;
for(;p;p>>=1,x=1LL*x*x%mod) if(p&1) ret = 1LL*ret*x % mod;
return ret;
}
inline void liner(int n){
mu[1] = minp[1] = oth[1] = 1;
rep(i,2,n){
if(!vis[i]){
pri[++cnt] = i;
mu[i] = -1;
minp[i] = i;
oth[i] = 1;
}
rep(j,1,cnt){
ll x = 1LL*i*pri[j];
if(x > n) break;
vis[x] = true;
if(i % pri[j] == 0){
minp[x] = minp[i];
oth[x] = oth[i];
mu[x] = 0;
break;
}
mu[x] = -mu[i];
minp[x] = pri[j];
oth[x] = i;
}
}
}
int g[maxn],f[maxn],d[maxn];
#define lowbit(x) (x&-x)
inline int calc(int T){
ll ups = f[T],dns = 1;
static int a[22],cnt;
cnt = 0;
for(int x=T;x!=1;x=oth[x])a[++cnt]=minp[x];
d[0] = 1;
for(int i=1;i<=cnt;++i) d[1<<i-1] = a[i];
for(int s=1,x;s<(1<<cnt);++s){
x = lowbit(s);
d[s] = d[s^x]*d[x];
if(mu[d[s]] == 1) ups = 1LL*ups*f[T/d[s]] % mod;
else dns = 1LL*dns*f[T/d[s]] % mod;
}
return 1LL*ups*qpow(dns,mod-2) % mod;
}
int main(){
liner(1000000);
int n,m;
n = 1000000;
f[0] = 0;f[1] = 1;
rep(i,2,n){
f[i] = f[i-1] + f[i-2];
if(f[i] >= mod) f[i] -= mod;
}
rep(i,1,n) g[i] = calc(i);
rep(i,2,n) g[i] = 1LL*g[i-1]*g[i] % mod;
g[0] = 1;
int T;read(T);
while(T--){
read(n);read(m);
if(n > m) swap(n,m);
int ans = 1;
for(int i=1,j;i <= n;i = j+1){
j = min(n / (n/i),m / (m/i));
ans =1LL*ans*qpow(1LL*g[j]*qpow(g[i-1],mod-2)%mod,1LL*(n/i)*(m/i)%(mod-1)) % mod;
}
printf("%d\n",ans);
}
return 0;
}