【SDOI2017】数字表格

题目大意

给定n,m，求 Πni=1Πmj=1fib[gcd(i,j)] \Pi_{i=1}^{n}\Pi_{j=1}^{m}fib[gcd(i,j)]
10% 1≤n,m≤100 1 ≤ n, m ≤ 100
30% 1≤n,m≤1000 1 ≤ n, m ≤ 1000
30% T≤3 T ≤ 3
100% T≤1000,1≤n,m≤106 T ≤ 1000,1 ≤ n, m ≤ 10^6

题目解析

一看见题，我就信心满满的打了发莫比乌斯反演正解，然后0分滚粗，还没有暴力容斥60分高。。。后来发现是对某些幂指数取了模的原因。

好的，开始正解：
Πni=1Πmj=1fib[gcd(i,j)] \Pi_{i=1}^{n}\Pi_{j=1}^{m}fib[gcd(i,j)]
我们们不妨令 f(x)=Πd|xg(d) f(x)=\Pi_{d|x}g(d)
则有 g(x)=f(x)Πd|x,d≠xg(d) g(x)=\frac{f(x)}{\Pi_{d|x,d\not=x}g(d)}
于是我们可以递推求出 f(x) f(x)和 g(x) g(x)。。。
然后再看原来的方程
Πni=1Πmj=1fib[gcd(i,j)] \Pi_{i=1}^{n}\Pi_{j=1}^{m}fib[gcd(i,j)]
=Πni=1Πmj=1Πd|i,d|jg(d) =\Pi_{i=1}^{n}\Pi_{j=1}^{m}\Pi_{d|i,d|j}g(d)
=Πmin(n,m)d=1g(d)⌊nd⌋×⌊md⌋ =\Pi_{d=1}^{min(n,m)}g(d)^{\lfloor\frac{n}{d}\rfloor\times\lfloor\frac{m}{d}\rfloor}
又因为 ⌊nd⌋×⌊md⌋ \lfloor\frac{n}{d}\rfloor\times\lfloor\frac{m}{d}\rfloor仅有 n−√+m−−√ \sqrt{n}+\sqrt{m}个取值，搞搞就行了。

代码

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<cstdlib>
#include<queue>
#include<stack>
using namespace std;#define MAXN 1000000
#define MAXM
#define INF 0x3f3f3f3f
typedef long long int LL;template<class T>
void Read(T &x){x=0;char c=getchar();bool flag=0;while(c<'0'||'9'<c){if(c=='-')flag=1;c=getchar();}while('0'<=c&&c<='9'){x=x*10+c-'0';c=getchar();}if(flag)x=-x;
}const int MOD = 1000000007;
int n,m;
int fib[MAXN+10],g[MAXN+10];
int sum[MAXN+10],inv[MAXN+10];int ksm(int a,int p,int mod){int ret=1;while(p){if(p&1)ret=1ll*ret*a%mod;a=1ll*a*a%mod;p>>=1;}return ret;
}void init(){g[1]=g[2]=fib[1]=fib[2]=1;for(int i=3;i<=MAXN;++i)g[i]=fib[i]=(fib[i-2]+fib[i-1])%MOD;for(int i=1;i<=MAXN;++i){int inv=ksm(g[i],MOD-2,MOD);for(int j=i*2;j<=MAXN;j+=i)g[j]=1ll*g[j]*inv%MOD;}inv[0]=sum[0]=1;for(int i=1;i<=MAXN;++i)sum[i]=1ll*sum[i-1]*g[i]%MOD,inv[i]=ksm(sum[i],MOD-2,MOD);
}int main(){//freopen("product.in","r",stdin);//freopen("product.out","w",stdout);init();int T;Read(T);while(T--){Read(n),Read(m);if(n>m)swap(n,m);int ans=1,last=0;for(int i=1;i<=n;i=last+1){last=min(n/(n/i),m/(m/i));ans=1ll*ans*ksm(1ll*sum[last]*inv[i-1]%MOD,1ll*(n/i)*(m/i)%(MOD-1),MOD)%MOD;}printf("%d\n",ans);}
}