BZOJ:4816: [Sdoi2017]数字表格

4816: [Sdoi2017]数字表格

Time Limit: 50 Sec Memory Limit: 128 MB
Submit: 501 Solved: 222
[Submit][Status][Discuss]

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

算是基础数论吧……

然后我就由于过于自信而gg，又WA又T

推一下式子很容易得到

$$ans=prod_{d=1}^{n}prod_{i|d}F_i^{mu(frac{d}{i})*leftlfloorfrac{n}{d} ight floor*leftlfloorfrac{m}{d} ight floor}$$

一开始我就拿这个直接分段上$O(n^{frac{3}{4}})$然后顺利地T掉了……

把$prod_{i|d}F_i^{mu(frac{d}{i})}$预处理出来就好了。

#include<cstdio>
#include<algorithm>
#define MN 1000001
using namespace std;
 
int read_p,read_ca;
inline int read(){
    read_p=0;read_ca=getchar();
    while(read_ca<'0'||read_ca>'9') read_ca=getchar();
    while(read_ca>='0'&&read_ca<='9') read_p=read_p*10+read_ca-48,read_ca=getchar();
    return read_p;
}
const int MOD=1e9+7;
int T,p[MN],num=0,mu[MN],F[MN],I[MN],n,m,mmh,S[MN];
bool bo[MN];
inline void M(int &x){while(x>=MOD)x-=MOD;while(x<0)x+=MOD;}
inline void _M(int &x){while(x>=MOD-1)x-=MOD-1;while(x<0)x+=MOD-1;}
inline int min(int x,int y){return x<y?x:y;}
inline int mi(int x,int y){
    int mmh=1;
    if (y<0) y+=MOD-1;
    for (;y;x=1LL*x*x%MOD,y>>=1) if (y&1) mmh=1LL*mmh*x%MOD;
    return mmh;
}
int main(){
    register int i,j,k;
    mu[1]=1;
    for (i=2;i<MN;i++){
        if (!bo[i]) p[++num]=i,mu[i]=-1;
        for (j=1;j<=num&&i*p[j]<MN;j++)
        if (bo[i*p[j]]=1,i%p[j])mu[i*p[j]]=-mu[i];else{mu[i*p[j]]=0;break;}
    }
    F[0]=0;F[1]=1;I[0]=I[1]=1;S[1]=1;
    for (i=2;i<MN;i++) M(F[i]=F[i-1]+F[i-2]),I[i]=mi(F[i],MOD-2),S[i]=1;
    for (i=1;i<MN;i++)
    if (mu[i])
    for (j=i,k=1;j<MN;j+=i,k++) S[j]=1LL*S[j]*(mu[i]==1?F[k]:I[k])%MOD;S[0]=1;
    for (i=1;i<MN;i++) S[i]=1LL*S[i]*S[i-1]%MOD;
    T=read();
    while(T--){
        n=read();m=read();
        if (n>m) swap(n,m);
        mmh=1;
        for (i=1;i<=n;i=j+1){
            j=min(n/(n/i),m/(m/i));
            mmh=1LL*mmh*mi(1LL*S[j]*mi(S[i-1],MOD-2)%MOD,1LL*(n/i)*(m/i)%(MOD-1))%MOD;
        }
        printf("%d
",mmh);
    }
}

View Code