【学术篇】SDOI2017数字表格 莫比乌斯反演

======传======送======门======在======里======面======
去年忘记可以预处理了… 然后就打了10pts的暴力… 现在学了莫比乌斯反演就可以来做了
= = = = = = = = = = = = = = = = = = = = = = = = = = = = = = == = = = = = = =

这个题目看着非常的简单, 就是要求这个式子

i=1Nj=1Mf[gcd(i,j)]

这个式子虽然是 而不是 , 但是还有个gcd, 我们也可以试着化出一个莫比乌斯反演的形式.
我们按照往常的套路来枚举gcd.
ans=i=1Nj=1Md=1min(n,m)f[d]b(gcd(i,j)=d)=d=1Nf[d]i=1Ndj=1Mdb(gcd(i,j)=1)=d=1Nf[d]Ndi=1Mdj=1b(gcd(i,j)=1)

然后指数上这一坨东西我们很 熟悉了, 我们得到过一个结论就是
i=1Nj=1M[gcd(i,j)=1]=x=1min(N,M)μ(x)NxMx

所以就有
ans=d=1Nf[d]min(Nd,Md)x=1μ(x)NxdMxd

然后令 t=xd , 就可以化成这个样纸
ans=t=1Nd|tf[d]min(Nd,Md)td=1μ(td)NtMt=t=1N(d|tf[d]min(Nd,Md)td=1μ(td))NtMt

这样括号外面我们会枚举分块, 只要把前缀和改成前缀积就ok了. 但是括号里面呢?
这个形式没见过啊, 好像筛法也没有筛这种东西的.
但是直觉告诉我们里面的这个东西好像不是很大, 可能大约是在 O(nlog2n)=O(nlogn) 级别的?
然后事实证明确实是差不多这样的(luogu题解里面说大约是15n左右), 所以我们做一波预处理就ok了.

这样总复杂度似乎差不多就是 O(n+nlogn+qn) 的, 不知道算的对不对..大概写O(能过)会更科学一点??

然后就是代码了…

#include <cmath>
#include <cstdio>
const int N=1e6+6;
const int P=1e9+7;
int prime[N],mu[N],f[N],g[N],F[N],tot;
bool notp[N];
inline int gn(int a=0,char c=0){
    for(;c<'0'||c>'9';c=getchar());
    for(;c>47&&c<58;c=getchar())a=a*10+c-48;return a;
}
template<class T>
inline T min(const T&a,const T&b){
    return a<b?a:b;
}
inline int qpow(int a,int b,int s=1){
    for(;b;b>>=1,a=1LL*a*a%P)
        if(b&1) s=1LL*s*a%P;
    return s;
}
void shai(int n){
    f[1]=mu[1]=notp[1]=F[0]=F[1]=1;
    for(int i=2;i<=n;++i){
        f[i]=(f[i-1]+f[i-2])%P; F[i]=1;
        if(!notp[i])prime[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot&&i*prime[j]<n;++j){
            int k=i*prime[j]; notp[k]=1;
            if(i%prime[j]==0){mu[k]=0;break;}
            else mu[k]=-mu[i];
        }
    }
    for(int i=1;i<=n;++i)
        for(int j=i;j<=n;j+=i)
            if(mu[j/i]) F[j]=1LL*F[j]*(mu[j/i]==1?f[i]:qpow(f[i],P-2))%P;
    for(int i=2;i<=n;++i)
        F[i]=1LL*F[i-1]*F[i]%P;
}
int solve(int n,int m){
    int ans=1,last,mn=min(n,m);
    for(int i=1;i<=mn;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ans=1LL*ans*qpow(1LL*F[last]*qpow(F[i-1],P-2)%P,
                1LL*(n/i)*(m/i)%(P-1))%P;
    }
    return (ans%P+P)%P;
}
int main(){ shai(1e6);
    int T=gn(),m,n;
    while(T--) 
        n=gn(),m=gn(),
        printf("%d\n",solve(n,m));
}

哦 对了 还有一件事就是根据费马小定理, G(P1)1(mod p) , 所以外面的指数对 (p1) 取模就好了…

  • 0
    点赞
  • 1
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值