斯德哥尔摩 的博客

斯德哥尔摩 的博客

P3704 [SDOI2017]数字表格

posted on 2018-05-13 19:31:16 | under 刷题 |

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)$提取出来就变成了: $Ans=\Pi_{d=1}^{n}f[d]^{\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]}$

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

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

将指数单独拿出来看: $\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]$

二话不说莫比乌斯反演: $\sum_{i=1}^{n/d} \mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor$

设 $T=id$,在整个式子里拎出来: $\Pi_{T=1}^{n}\Pi_{d|T}f[d]^{\mu(T/d)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$

化简一下: $\Pi_{T=1}^{n}(\Pi_{d|T}f[d]^{\mu(T/d)})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}$

很明显,对那个啥 $\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;
}