bzoj 4816: [Sdoi2017]数字表格

Description
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取模。

解题报告:
这题和之前的差不多.设\(g[x]\)为\(gcd==x\)的数对个数，然后就是基本的反演了
\(\sum_{d=1}^nf[d]^{g[d]}\)
\(g[d]=\sum_{i|d}\lfloor m/d\rfloor\lfloor n/d\rfloor\mu(d/i)\)
这个明显可以数论分块，所以我们处理一下\(g[d]^\mu\)的前缀积,对于除法再求个逆即可

#include <algorithm>
#include <iostream>
#include <cstdlib>
#include <cstring>
#include <cstdio>
#include <cmath>
#define RG register
#define il inline
#define iter iterator
#define Max(a,b) ((a)>(b)?(a):(b))
#define Min(a,b) ((a)<(b)?(a):(b))
using namespace std;
typedef long long ll;
const int N=1e6+5,mod=1e9+7;
int q[1005][2],lim,prime[N],num=0,pi[N];bool vis[N];
void solve(){
    pi[1]=1;int to;
    for(int i=2;i<=lim;i++){
        if(!vis[i]){
            prime[++num]=i;
            pi[i]=-1;
        }
        for(int j=1;j<=num && i*prime[j]<=lim;j++){
            to=i*prime[j];vis[to]=true;
            if(i%prime[j])pi[to]=-pi[i];
            else{
                pi[to]=0;
                break;
            }
        }
    }
}
ll f[N],sum[N];
ll qm(ll x,ll k){
    ll sum=1;
    while(k){
        if(k&1)sum*=x,sum%=mod;
        x*=x;x%=mod;k>>=1;
    }
    return sum;
}
ll query(int n,int m){
    int j;ll ret=1;
    for(int i=1;i<=n;i=j+1){
        j=min(n/(n/i),m/(m/i));
        ret=ret*qm(sum[j]*qm(sum[i-1],mod-2)%mod,(ll)(n/i)*(m/i))%mod;
    }
    return ret;
}
ll ni[N],g[N];
void work()
{
    int Q;cin>>Q;
    for(int i=1;i<=Q;i++){
        scanf("%d%d",&q[i][0],&q[i][1]);
        if(q[i][0]>q[i][1])swap(q[i][0],q[i][1]);
        lim=Max(lim,q[i][0]);
    }
    solve();
    f[0]=0;f[1]=1;
    for(int i=2;i<=lim;i++){
        f[i]=f[i-1]+f[i-2];
        if(f[i]>=mod)f[i]-=mod;
    }
    for(int i=0;i<=lim;i++)ni[i]=qm(f[i],mod-2);
    for(int i=1;i<=lim;i++)g[i]=1;
    for(int i=1;i<=lim;i++)
        for(int j=i;j<=lim;j+=i)
            if(pi[j/i]==1)g[j]=g[j]*f[i]%mod;
            else if(pi[j/i]==-1)g[j]=g[j]*ni[i]%mod;
    sum[0]=1;
    for(int i=1;i<=lim;i++){
        sum[i]=sum[i-1]*g[i]%mod;
    }
    for(int i=1;i<=Q;i++){
        printf("%lld\n",query(q[i][0],q[i][1]));
    }
}

int main()
{
    work();
    return 0;
}