SDOI2017数字表格

题意：求\(\prod^n_{i=1}\prod^m_{j=1}fib[gcd(i,j)],fib[1]=fib[2]=1\)

\(T\)组数组\(n,m\in[1,1e6],T\in[1,1e3]\)

SOL:

易得\(ans=\prod^n_{g=1}fib(g)^{S(n/g,m/g)}\)

\(S(n,m)=\sum^n_{i=1}\sum^m_{j=1}[gcd(i,j)]==1=\sum^n_{d=1}\mu(d)(n/d)(m/d)\)

默认\(n<=m\)

注意：

1. 结合2个式子一起
2. 新建变量代替乘积，设\(T=gd\)

\(ans=\prod^n_{T=1}(\prod_{g|T}fib(g)^{\mu(\frac{T}{g})})^{(n/T)(m/T)}\)

加上预处理即可

时间复杂度\(O(nlog_n+T\sqrt nlog_n)\)

#include<bits/stdc++.h>
using namespace std;
inline int read(){
    int x=0,f=1;char c=getchar();
    while(!isdigit(c)){if(c=='-')f=-1;c=getchar();}
    while(isdigit(c)){x=(x<<1)+(x<<3)+(c^48);c=getchar();}
    return f==1?x:-x;
}
#define ll long long
const int N=1e6+4,mod=1e9+7;
inline int ksm(int x,ll r){
    int ret=1;
    for(int i=0;(1ll<<i)<=r;i++){
        if((r>>i)&1)ret=(ll)ret*x%mod;
        x=(ll)x*x%mod; 
    }
    return ret;
}
int n,m,mu[N],pri[N],b[N],fib[N],s[N],cnt,ans;
inline void pre(){
    mu[1]=1;
    for(int i=2;i<N;i++){
        if(!b[i]){pri[++cnt]=i;mu[i]=-1;}
        for(int j=1;j<=cnt&&i*pri[j]<N;j++){
            b[i*pri[j]]=1;
            if(i%pri[j])mu[i*pri[j]]=-mu[i];
            else{
                mu[i*pri[j]]=0;
                break;
            }
        }
    }
    for(int i=0;i<N;i++)s[i]=1;
    fib[1]=fib[2]=1;
    for(int i=3;i<N;i++)fib[i]=(fib[i-1]+fib[i-2])%mod;
    for(int i=1,inv;i<N;i++){
        inv=ksm(fib[i],mod-2);
        for(int j=i;j<N;j+=i)
            if(mu[j/i]==1)s[j]=(ll)s[j]*fib[i]%mod;
            else if(mu[j/i]==-1)s[j]=(ll)s[j]*inv%mod; 
    }
    for(int i=2;i<N;i++)s[i]=(ll)s[i]*s[i-1]%mod;
}
int main(){
    pre();
    int T=read();
    while(T--){
        n=read();m=read();
        if(n>m)n^=m^=n^=m;
        ans=1;
        for(int l=1,r;l<=n;l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans=(ll)ans*ksm((ll)s[r]*ksm(s[l-1],mod-2)%mod,(ll)(n/l)*(m/l))%mod;
        }
        cout<<ans<<"\n";
    } 
    return (0-0);
}

标签：mu,数字,表格,int,sum,ans,fib,SDOI2017,prod
来源： https://www.cnblogs.com/aurora2004/p/12505840.html