原文链接http://www.cnblogs.com/zhouzhendong/p/8666106.html
题目传送门 - BZOJ4816
题意
定义$f(0)=0,f(1)=1,f(i)=f(i-1)+f(i-2)$。
$T$组数据，每组数据两个整数$n,m$，求$prod_{i=1}^nprod_{j=1}^m f(gcd(i,j))$。
$Tleq 1000,1leq n,m leq 10^6$
题解
先推一波式子。
$$prod_{i=1}^nprod_{j=1}^m f(gcd(i,j))\=prod_{d=1}^n f(d)^{sum_{i=1}^{leftlfloorfrac nd ight floor}sum_{j=1}^{leftlfloorfrac md ight floor}[gcd(i,j)=1]}\=prod_{d=1}^n f(d)^{sum_{i=1}^{leftlfloorfrac nd ight floor}sum_{j=1}^{leftlfloorfrac md ight floor}sum_{p|i,p|j}mu(p)}\=prod_{d=1}^n f(d)^{sum_{p=1}^{leftlfloorfrac nd ight floor}mu(p)leftlfloorfrac{n}{pd} ight floorleftlfloorfrac{m}{pd} ight floor}$$
设$D=pd$。
$$prod_{i=1}^nprod_{j=1}^m f(gcd(i,j))\=prod_{d=1}^n f(d)^{sum_{p=1}^{leftlfloorfrac nd ight floor}mu(p)leftlfloorfrac{n}{pd} ight floorleftlfloorfrac{m}{pd} ight floor}\=prod_{D=1}^{n}(prod_{d|D}f(d)^{mu(frac Dd)})^{leftlfloorfrac nD ight floorleftlfloorfrac mD ight floor}$$
先顺手预处理$mu$。
于是我们先预处理$f(x)$以及$f(x)$的逆元，然后再$O(n log n)$预处理出所有$g(D)=prod_{d|D}f(d)^{mu(frac Dd)}$。
然后再预处理出$g$的前缀积以及前缀积的逆元。
这些的复杂度都是$O(n log n)$。
然后回答一个询问的时候再整除分块一下，单次询问复杂度为$O(sqrt n log n)$。
所以总复杂度为$O(n log n+Tsqrt n log n)$。
所有测试点的总时限开了50s。
哈哈应该稳过了吧！
但是！！！！！！
BZOJ毒！瘤！卡！常！！！！！！
于是我们需要常数优化。
考虑处理一个序列的逆元。
我们求出当前序列的前缀积是$O(n)$的。
记$a$为原序列。
记$inv_i$为的$i$项的逆元。
记$presum_i$为序列前$i$项的前缀积。
记$preinv_i$为序列钱$i$项的前缀积的逆元。
则：
$$inv_i=preinv_i imes presum_{i-1}$$
$$preinv_{i-1}=preinv_i imes a_i$$
于是你就可以倒着来求逆元了，复杂度$O(n)$。
但是求$g$的时候复杂度还是$O(n log n)$。
所以总的复杂度还是$O(n log n+Tsqrt n log n)$。
但是你卡常数了！！！
你过了！！！
QAQ
代码
#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const int N=1e6+5,mod=1e9+7;
LL prime[N],u[N],pcnt;
bool isprime[N];
LL f[N],pref[N],invf[N],g[N],preg[N],invg[N];
int T,n,m;
LL Pow(LL a,LL b){
if (!b)
return 1LL;
LL x=Pow(a,b/2);
x=x*x%mod;
if (b&1LL)
x=x*a%mod;
return x;
}
LL Inv(LL a){
return Pow(a,mod-2);
}
void get_prime(int n){
memset(isprime,true,sizeof isprime);
u[1]=1,pcnt=isprime[0]=isprime[1]=0;
for (int i=2;i<=n;i++){
if (isprime[i])
prime[++pcnt]=i,u[i]=-1;
for (int j=1;j<=pcnt&&i*prime[j]<=n;j++){
isprime[i*prime[j]]=0;
if (i%prime[j])
u[i*prime[j]]=-u[i];
else {
u[i*prime[j]]=0;
break;
}
}
}
}
void init(int n){
get_prime(n);
f[0]=0,f[1]=invf[1]=g[0]=1;
for (int i=2;i<=n;i++)
f[i]=(f[i-1]+f[i-2])%mod;
pref[0]=1;
for (int i=1;i<=n;i++)
pref[i]=pref[i-1]*f[i]%mod;
LL v=Inv(pref[n]);
for (int i=n;i>=1;i--)
invf[i]=v*pref[i-1]%mod,v=v*f[i]%mod;
for (int i=1;i<=n;i++)
g[i]=1;
for (int i=1;i<=n;i++)
for (int j=1;i*j<=n;j++){
if (u[j]==-1)
g[i*j]=g[i*j]*invf[i]%mod;
if (u[j]==1)
g[i*j]=g[i*j]*f[i]%mod;
}
preg[0]=1;
for (int i=1;i<=n;i++)
preg[i]=preg[i-1]*g[i]%mod;
invg[n]=Inv(preg[n]);
for (int i=n-1;i>=0;i--)
invg[i]=invg[i+1]*g[i+1]%mod;
}
int main(){
init(1e6);
scanf("%d",&T);
while (T--){
scanf("%d%d",&n,&m);
if (n>m)
swap(n,m);
LL ans=1;
for (int D=1,i;D<=n;D=i+1){
i=min(n/(n/D),m/(m/D));
ans=ans*Pow(preg[i]*invg[D-1]%mod,1LL*(n/D)*(m/D)%(mod-1))%mod;
}
printf("%lld
",ans);
}
return 0;
}