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

题意

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

$$f[0]=0$$

$$f[1]=1$$
$$f[n]=f[n-1]+f[n-2],n\geq 2$$
有一个 $n×m$的表格，第 $i$行第$j$列的格子中的数是 $f[\gcd(i,j)]$，其中 $\gcd(i,j)$表示 $i,j$的最大公约数。求这n*m个数的乘积是多少。
答案对 $10^9+7$取模。
对 $100\%$的数据， $T\leq1000,1\leq n,m\leq 10^6$

分析

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

我们设n<=m，显然

$$ans=\prod_{i=1}^n\prod_{j=1}^mf(gcd(i,j))$$
设 $$d=gcd(i,j),sum(d)=\sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)==d]$$
通过简单的反演可以得到 $$sum(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{di}\rfloor\lfloor\frac{m}{di}\rfloor\mu(i)$$
显然有 $$ans=\prod_{d=1}^nf(d)^{sum(d)}$$
把上面两条式子合并在一起 $$ans=\prod_{d=1}^n\prod_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(d)^{\lfloor\frac{n}{di}\rfloor\lfloor\frac{m}{di}\rfloor\mu(i)}$$
设 $T=di$ $$ans=\prod_{T=1}^n\prod_{d|T}f(d)^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor\mu(\frac{T}{d})}$$ 注意到我们可以把两个下取整提出来 $$ans=\prod_{T=1}^n(\prod_{d|T}f(d)^{\mu(\frac{T}{d})})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$$
然后就大功告成啦！
注意到对于那两个下去整，可以将 $[1,n]$分为 $2(\sqrt n+\sqrt m)$段，使得每段的指数是相同的。
那么我们只要O(nlogn)预处理出括号内那部分的前缀积，每次快速幂即可。
时间复杂度 $O(nlogn+T\sqrt nlogn)$

代码

#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));
    }
}