P3704 [SDOI2017]数字表格

斯德哥尔摩 · 2018-05-13 19:31:16 · 个人记录

P3704 [SDOI2017]数字表格

题目要求：Ans=\Pi_{i=1}^{n}\Pi_{j=1}^{m}f[gcd(i,j)]

不妨设n<=m

则Ans=\Pi_{d=1}^{n}\Pi_{i=1}^{n}\Pi_{j=1}^{m}f[d]\times [gcd(i,j)==d]

直接把那啥f(d)提取出来就变成了：

显然，那个指数用莫比乌斯+数论分块可以O(\sqrt n)做，外面的一层也可以，所以总复杂度为O(n)

但是显然不够，需要变形。

将指数单独拿出来看：

二话不说莫比乌斯反演：

设T=id，在整个式子里拎出来：

化简一下：

很明显，对那个啥\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor数论分块。

那——里面的东西怎么办？又不能线性筛。。。

等等——不能线性筛就暴力算呀！

数据范围10^6

每个数暴力算到他的倍数里面去

即：\frac{n}{1}+\frac{n}{2}+\frac{n}{3}+...+\frac{n}{10^6}

这玩意差不多15n，所以直接暴力把那个东西的前缀给求出来，就可以做到O(\sqrt n)了。

附代码：

#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 1000010
#define MOD 1000000007
using namespace std;
long long f[MAXN],g[MAXN],sum[MAXN];//直接开long long
int k=0,prime[MAXN],mu[MAXN];
bool np[MAXN];
inline int read(){
    int date=0,w=1;char c=0;
    while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
    while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
    return date*w;
}
long long mexp(long long a,long long b,long long c){//快速幂
    long long s=1;
    while(b){
        if(b&1)s=s*a%c;
        a=a*a%c;
        b>>=1;
    }
    return s;
}
void make(){//莫比乌斯专属线性筛
    int m=MAXN-10;
    f[1]=g[1]=sum[0]=sum[1]=mu[1]=1;
    np[1]=true;
    for(int i=2;i<=m;i++){
        f[i]=(f[i-1]+f[i-2])%MOD;
        g[i]=mexp(f[i],MOD-2,MOD);
        sum[i]=1;
        if(!np[i]){
            prime[++k]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=k&&prime[j]*i<=m;j++){
            np[prime[j]*i]=true;
            if(i%prime[j]==0)break;
            mu[i*prime[j]]=-mu[i];
        }
    }
    for(int i=1;i<=m;i++){//暴力求前缀
        if(!mu[i])continue;
        for(int j=i;j<=m;j+=i)sum[j]=sum[j]*(mu[i]==1?f[j/i]:g[j/i])%MOD;
    }
    for(int i=2;i<=m;i++)sum[i]=sum[i]*sum[i-1]%MOD;
}
void work(){//求解
    int n,m;
    long long inv,ans=1;
    n=read();m=read();
    if(n>m)swap(n,m);
    for(int i=1,last;i<=n;i=last+1){
        last=min(n/(n/i),m/(m/i));
        inv=sum[last]*mexp(sum[i-1],MOD-2,MOD)%MOD;//求逆元
        ans=ans*mexp(inv,1LL*(n/i)*(m/i)%(MOD-1),MOD)%MOD;
    }
    printf("%lld\n",ans);
}
int main(){
    int t=read();
    make();
    while(t--)work();//多组询问
    return 0;
}