bzoj 4816: [Sdoi2017]数字表格 莫比乌斯反演

题意

定义fibonacci数列。用f[i]表示数列的第i项，那么

f[0]=0

f[1]=1

f[n]=f[n−1]+f[n−2],n≥2

有一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,j的最大公约数。求这n*m个数的乘积是多少。
答案对109+7取模。
对100%的数据，T≤1000,1≤n,m≤106

分析

一开始推公式的时候习惯性地写了∑上去，推了半天才有另外一个人过来告诉我说，这题求的是乘积。。。

我们设n<=m，显然

ans=∏i=1n∏j=1mf(gcd(i,j))

设

d=gcd(i,j),sum(d)=∑i=1n∑j=1m[gcd(i,j)==d]

通过简单的反演可以得到

sum(d)=∑i=1⌊nd⌋⌊ndi⌋⌊mdi⌋μ(i)

显然有

ans=∏d=1nf(d)sum(d)

把上面两条式子合并在一起

ans=∏d=1n∏i=1⌊nd⌋f(d)⌊ndi⌋⌊mdi⌋μ(i)

设T=di

ans=∏T=1n∏d|Tf(d)⌊nT⌋⌊mT⌋μ(Td)

注意到我们可以把两个下取整提出来

ans=∏T=1n(∏d|Tf(d)μ(Td))⌊nT⌋⌊mT⌋

然后就大功告成啦！
注意到对于那两个下去整，可以将[1,n]分为2(n√+m−−√)段，使得每段的指数是相同的。
那么我们只要O(nlogn)预处理出括号内那部分的前缀积，每次快速幂即可。
时间复杂度O(nlogn+Tn√logn)

代码

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>using namespace std;typedef long long LL;const int N=1000005;const int MOD=1000000007;int tot,prime[N],f[N],s[N],ny[N],mu[N];bool not_prime[N];int ksm(int x,int y){    int ans=1;    while (y)    {        if (y&1) ans=(LL)ans*x%MOD;        x=(LL)x*x%MOD;y>>=1;    }    return ans;}void prework(int n){    mu[1]=1;    for (int i=2;i<=n;i++)    {        if (!not_prime[i]) prime[++tot]=i,mu[i]=-1;        for (int j=1;j<=tot&&i*prime[j]<=n;j++)        {            not_prime[i*prime[j]]=1;            if (i%prime[j]==0)            {                mu[i*prime[j]]=0;                break;            }            mu[i*prime[j]]=-mu[i];        }    }    f[0]=0;f[1]=1;    ny[0]=ny[1]=1;    for (int i=2;i<=n;i++) f[i]=(f[i-1]+f[i-2])%MOD,ny[i]=ksm(f[i],MOD-2);    for (int i=1;i<=n;i++) s[i]=1;    for (int i=1;i<=n;i++)        for (int j=i;j<=n;j+=i)            if (mu[j/i]==1) s[j]=(LL)s[j]*f[i]%MOD;            else if (mu[j/i]==-1) s[j]=(LL)s[j]*ny[i]%MOD;    for (int i=2;i<=n;i++) s[i]=(LL)s[i]*s[i-1]%MOD,ny[i]=ksm(s[i],MOD-2);}int solve(int n,int m){    if (n>m) swap(n,m);    int ans=1;    for (int i=1,last;i<=n;i=last+1)    {        last=min(n/(n/i),m/(m/i));        ans=(LL)ans*ksm((LL)s[last]*ny[i-1]%MOD,(LL)(n/i)*(m/i)%(MOD-1))%MOD;    }    return (ans+MOD)%MOD;}int main(){    prework(1000000);    int T;    scanf("%d",&T);    while (T--)    {        int n,m;        scanf("%d%d",&n,&m);        printf("%d\n",solve(n,m));    }}