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})*\left\lfloor\frac{n}{d}\right\rfloor*\left\lfloor\frac{m}{d}\right\rfloor}$$

一开始我就拿这个直接分段上$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=;read_ca=getchar();

    while(read_ca<''||read_ca>'') read_ca=getchar();

    while(read_ca>=''&&read_ca<='') read_p=read_p*+read_ca-,read_ca=getchar();

    return read_p;

}

const int MOD=1e9+;

int T,p[MN],num=,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<)x+=MOD;}

inline void _M(int &x){while(x>=MOD-)x-=MOD-;while(x<)x+=MOD-;}

inline int min(int x,int y){return x<y?x:y;}

inline int mi(int x,int y){

    int mmh=;

    if (y<) y+=MOD-;

    for (;y;x=1LL*x*x%MOD,y>>=) if (y&) mmh=1LL*mmh*x%MOD;

    return mmh;

}

int main(){

    register int i,j,k;

    mu[]=;

    for (i=;i<MN;i++){

        if (!bo[i]) p[++num]=i,mu[i]=-;

        for (j=;j<=num&&i*p[j]<MN;j++)

        if (bo[i*p[j]]=,i%p[j])mu[i*p[j]]=-mu[i];else{mu[i*p[j]]=;break;}

    }

    F[]=;F[]=;I[]=I[]=;S[]=;

    for (i=;i<MN;i++) M(F[i]=F[i-]+F[i-]),I[i]=mi(F[i],MOD-),S[i]=;

    for (i=;i<MN;i++)

    if (mu[i])

    for (j=i,k=;j<MN;j+=i,k++) S[j]=1LL*S[j]*(mu[i]==?F[k]:I[k])%MOD;S[]=;

    for (i=;i<MN;i++) S[i]=1LL*S[i]*S[i-]%MOD;

    T=read();

    while(T--){

        n=read();m=read();

        if (n>m) swap(n,m);

        mmh=;

        for (i=;i<=n;i=j+){

            j=min(n/(n/i),m/(m/i));

            mmh=1LL*mmh*mi(1LL*S[j]*mi(S[i-],MOD-)%MOD,1LL*(n/i)*(m/i)%(MOD-))%MOD;

        }

        printf("%d\n",mmh);

    }

}

秒客网