【BZOJ 4816】 4816: [Sdoi2017]数字表格 （莫比乌斯）

4816: [Sdoi2017]数字表格

Time Limit: 50 Sec  Memory Limit: 128 MB
Submit: 666  Solved: 312

Description

Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么

f[0]=0

f[1]=1

f[n]=f[n-1]+f[n-2],n>=2

Doris用老师的超级计算机生成了一个n×m的表格，第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,

j的最大公约数。Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。答案对10^9+7取模。

Input

有多组测试数据。

第一个一个数T，表示数据组数。

接下来T行，每行两个数n,m

T<=1000,1<=n,m<=10^6

Output

输出T行，第i行的数是第i组数据的结果

Sample Input

3
2 3
4 5
6 7

Sample Output

1
6
960

HINT

Source

鸣谢infinityedge上传

【分析】

　　额。。。加法变乘法有时候还是忍不住求和。。

　　推式子。。

　　$$Ans=Pi Pi f(gcd(i,j))$$

　　$$=Pi f(d)^{sum 1 [gcd(i,j)==d]}$$

　　$$=Pi f(d)^{sum_{i=1}^{n/d}sum_{j=1}^{m/d} 1 [gcd(i,j)==1]}$$

　　$$=Pi f(d)^{sum mu(d')*lfloordfrac{n}{d*d'} floor lfloordfrac{m}{d*d'} floor}$$

　　【啊好辛苦

　　设D=d*d'

　　$$Ans=Pi_{D}Pi_{d|D} f(d)^{lfloordfrac{n}{D} floor lfloordfrac{m}{D} floor*mu(dfrac{D}{d})}$$

　　$$=Pi Pi (f(d)^{mu(D/d)})^{lfloordfrac{n}{D} floor lfloordfrac{m}{D} floor}$$

　　设$g(D)=Pi_{d|D}f(d)^{mu(D/d)}$

　　这个$O(nlogn)$预处理，后面的D在询问时$sqrt n $分块即可。

　　因为有快速幂，所以是$O(nlogn)+T*sqrt n*logn$

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<algorithm>
using namespace std;
#define Maxn 1000010
#define Mod 1000000007

int g[Maxn],pri[Maxn],pl,mu[Maxn],f[Maxn];
int inv[Maxn],sm[Maxn];
bool vis[Maxn];

int qpow(int x,int b)
{
    int ans=1;
    while(b)
    {
        if(b&1) ans=1LL*ans*x%Mod;
        x=1LL*x*x%Mod;
        b>>=1;
    }
    return ans;
}

void init()
{
    f[0]=0;f[1]=1;for(int i=2;i<=Maxn-10;i++) f[i]=(f[i-1]+f[i-2])%Mod;
    for(int i=1;i<=Maxn-10;i++) inv[i]=qpow(f[i],Mod-2);
    memset(vis,0,sizeof(vis));
    mu[1]=1;
    for(int i=2;i<=Maxn-10;i++)
    {
        if(!vis[i]) pri[++pl]=i,mu[i]=-1;
        for(int j=1;j<=pl;j++)
        {
            if(i*pri[j]>Maxn-10) break;
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) {mu[i*pri[j]]=0;break;}
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(int i=0;i<=Maxn-10;i++) g[i]=1;
    sm[0]=1;
    int i;
    for(i=1;i<=Maxn-10;i++)
    {
        for(int j=i;j<=Maxn-10;j+=i) if(mu[j/i]!=0)
        {
            if(mu[j/i]==-1) g[j]=1LL*g[j]*inv[i]%Mod;
            else g[j]=1LL*g[j]*f[i]%Mod;
        }
        sm[i]=1LL*sm[i-1]*g[i]%Mod;
    }
}


int main()
{
    init();
    int T;
    scanf("%d",&T);
    while(T--)
    {
        int n,m;
        scanf("%d%d",&n,&m);
        if(n>m) swap(n,m);
        int ans=1;
        for(int i=1;i<=n;)
        {
            int x=n/i,y=m/i,r1=n/x,r2=m/y;
            r1=min(r1,r2);
            ans=1LL*ans*qpow(1LL*sm[r1]*qpow(sm[i-1],Mod-2)%Mod,1LL*x*y%(Mod-1))%Mod;
            i=r1+1;
        }
        printf("%d
",ans);
    }
    return 0;
}

2017-04-26 20:21:03