BZOJ 4816[SDOI2017]數字表格(莫比烏斯反演)

題目連接php

\(Description\)
\(f_i\)表示\(fibonacci\)數列第\(i\)項,求\(\prod_{i=1}^{n}\prod_{j=1}^{m}f[gcd(i,j)]\)算法

\(T<=10^3,n,m≤10^6\)spa


\(Solution\)
再來推式子(默認\(n<m\))
\[\prod_{i=1}^{n}\prod_{j=1}^mf[gcd(i,j)]\]
按照套路枚舉\(gcd\)
\[\prod_{d=1}^n\prod_{i=1}^{n/d}\prod_{j=1}^{m/d}f[d][gcd(i,j)==1]\]
考慮每一個\(f[d]\)被乘了多少次
\[\prod_{d=1}^nf[d]^{\sum_{i=1}^{n/d}\sum_{j=1}^{m/d}[gcd(i,j)==1]}\]
指數很是熟悉
\[\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\left\lfloor\frac{m}{id}\right\rfloor\]
如今指數能夠數論分塊,對\(f\)作一個前綴積以後外層也能夠數論分塊,這就能夠作到每次詢問\(O(nlogn)\)了。
可是這個式子還能夠繼續推,咱們能夠枚舉\(id\)
\(T=id\)
\[\prod_{T=1}^n\prod_{d|T}f[d]^{\mu(T/d)\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}\]
根據數學運算法則就能夠這樣算
\[\prod_{T=1}^n(\prod_{d|T}f[d]^{\mu(T/d)})^{\lfloor\frac{n}{T}\rfloor\lfloor\frac{m}{T}\rfloor}\]
括號裏面的式子能夠預處理,乘方能夠數論分塊,這樣就能夠單詞詢問\(O(\sqrt n logn)\)了。code

#include<complex>
#include<cstdio>
using namespace std;
const int mod=1e9+7;
const int N=1e6+7;
int Q,n,m,tot;
int prime[N],mu[N],f[N],g[N],pro[N];
bool check[N];
int qread()
{
    int x=0;
    char ch=getchar();
    while(ch<'0' || ch>'9')ch=getchar();
    while(ch>='0' && ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x;
}
int Fpow(long long b,int p)
{
    long long res=1;
    for(;p;p>>=1,b=b*b%mod)
        if(p&1)res=res*b%mod;
    return res;
}
void Init()
{
    check[1]=1;
    mu[1]=f[1]=pro[0]=pro[1]=g[1]=1;
    for(int i=2;i<N;i++)
    {
        if(!check[i])prime[++tot]=i,mu[i]=-1;
        for(int j=1;j<=tot && i*prime[j]<N;j++)
        {
            check[i*prime[j]]=1;
            if(i%prime[j])mu[i*prime[j]]=-mu[i];
            else break;
        }
        f[i]=(f[i-1]+f[i-2])%mod;pro[i]=1;
        g[i]=Fpow(f[i],mod-2);
    }
    for(int i=1;i<N;i++)
    {
        if(!mu[i])continue;
        for(int j=i;j<N;j+=i)
            pro[j]=1ll*pro[j]*(mu[i]==1?f[j/i]:g[j/i])%mod;
    }
    for(int i=2;i<N;i++)
        pro[i]=1ll*pro[i]*pro[i-1]%mod;
}
int main()
{
    Init();
    scanf("%d",&Q);
    while(Q--)
    {
        n=qread();m=qread();
        if(n>m)swap(n,m);
        int ans=1;
        for(int l=1,r;l<=n;l=r+1)
        {
            r=min(n/(n/l),m/(m/l));
            ans=1ll*ans*Fpow(1ll*pro[r]*Fpow(pro[l-1],mod-2)%mod,1ll*(n/l)*(m/l)%(mod-1))%mod;
        }
        printf("%d\n",ans);
    }
    return 0;
}
相關文章
相關標籤/搜索