[bzoj4816][SDOI2017]数字表格

来源：互联网发布：有线网络信号增强器编辑：程序博客网时间：2024/05/21 11:14

题目描述

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么
f[0]=0
f[1]=1
f[n]=f[n-1]+f[n-2],n>=2
Doris用老师的超级计算机生成了一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,
j的最大公约数。Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。答案对10^9+7取模。

推式子

首先令n<=m（套路）
Πni=1Πmj=1fib[(i,j)]
Πnd=1fib[d]∑ni=1∑mj=1[(i,j)=d]
Πnd=1fib[d]∑⌊nd⌋i=1μ(i)∗⌊nd∗i⌋∗⌊md∗i⌋
这样单次可以线性。
但是还不够，考虑套路令T=d*i，然后交换主体。
ΠnT=1(Πd|Tfib[d]μ(Td))⌊nT⌋∗⌊mT⌋
中间只与T有关，设为a[T]，可以n log n预处理。
而这个式子可以单次根号

#include<cstdio>#include<algorithm>#define fo(i,a,b) for(i=a;i<=b;i++)using namespace std;typedef long long ll;const int maxn=1000000+10,mo=1000000007;int a[maxn],pri[maxn],mu[maxn],fib[maxn],inv[maxn],sum[maxn];bool bz[maxn];int i,j,k,l,t,n,m,top,ans,ca;int qsm(int x,int y){    if (!y) return 1;    int t=qsm(x,y/2);    t=(ll)t*t%mo;    if (y%2) t=(ll)t*x%mo;    return t;}void prepare(){    mu[1]=1;    fo(i,2,maxn-10){        if (!bz[i]) pri[++top]=i,mu[i]=-1;        fo(j,1,top){            if ((ll)i*pri[j]>maxn-10) break;            bz[i*pri[j]]=1;            if (i%pri[j]==0){                mu[i*pri[j]]=0;                break;            }            mu[i*pri[j]]=-mu[i];        }    }    fib[1]=1;    fo(i,2,maxn-10) fib[i]=(fib[i-1]+fib[i-2])%mo;    fo(i,1,maxn-10) inv[i]=qsm(fib[i],mo-2);    fo(i,1,maxn-10) a[i]=1;    fo(i,1,maxn-10)        fo(j,1,(maxn-10)/i){            if (!mu[j]) continue;            if (mu[j]==1) a[i*j]=(ll)a[i*j]*fib[i]%mo;            else a[i*j]=(ll)a[i*j]*inv[i]%mo;        }    sum[0]=1;    fo(i,1,maxn-10) sum[i]=(ll)sum[i-1]*a[i]%mo;}int main(){    freopen("product.in","r",stdin);freopen("product.out","w",stdout);    prepare();    scanf("%d",&ca);    while (ca--){        scanf("%d%d",&n,&m);        if (n>m) swap(n,m);        ans=1;        i=1;        while (i<=n){            j=min(n/(n/i),m/(m/i));            t=(ll)sum[j]*qsm(sum[i-1],mo-2)%mo;            ans=(ll)ans*qsm(t,(ll)(n/i)*(m/i)%(mo-1))%mo;            i=j+1;        }        (ans+=mo)%=mo;        printf("%d\n",ans);    }}

0 0