[SDOI2017]数字表格

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

$$prod_{d=1}^nprod_{i=1}^nprod_{j=1}^mif(gcd(i,j)==d)f[gcd(i,j)]$$

直接把$f[d]$提出来

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

上面那个东西用莫比乌斯反演+数论分块可以$O(sqrt n)$求

外面套的这一层也可以数论分块求

于是，我们就得到了一个$O(n)$的做法

但是显然还不够

把上面那坨东西拎出来看

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

太熟悉了

$$sum_{i=1}^{n/d}mu(i)[frac{n}{id}][frac{m}{id}]$$

还是老套路，

令$T=id$

直接把$T$在整个式子里面提出来

$$prod_{T=1}^{n}prod_{d|T}f[d]^{[n/T][m/T]mu(T/d)}$$

有一些一样的东西

$$prod_{T=1}^{n}(prod_{d|T}f[d]^{mu(T/d)})^{[n/T][m/T]}$$

然后怎么办。。。。

很明显，已经可以对$[n/T][m/T]$分块了

里面的东西不好处理，那么就暴力预处理F(T)数组

枚举i,j(i*j<=N)那么复杂度：

$frac{n}{1}+frac{n}{2}+.....frac{n}{10^6}$

复杂度O(nlogn+Q√n)

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cstring>
  4 #include<algorithm>
  5 #include<cmath>
  6 using namespace std;
  7 typedef long long lol;
  8 int Q[1001][2];
  9 const int NN=1000000;
 10 int N=0;
 11 int Mod=1e9+7;
 12 int f[NN+5],F[NN+5],inv[NN+5],Inv[NN+5],ans;
 13 int mu[NN+5],tot,prime[NN],n,m;
 14 bool vis[NN+5];
 15 int qpow(int x,int y)
 16 {
 17   int res=1;
 18   while (y)
 19     {
 20       if (y&1) res=1ll*res*x%Mod;
 21       x=1ll*x*x%Mod;
 22       y>>=1;
 23     }
 24   return res;
 25 }
 26 int pow(int x1,int x2,int y)
 27 {
 28   if (y<0)
 29     {
 30       return x2;
 31     }
 32   if (y==0) return 1;
 33   if (y>0)
 34     {
 35       return x1;
 36     }
 37 }
 38 void pre()
 39 {int i,j;
 40   mu[1]=1;
 41   for (i=2;i<=N;i++)
 42     {
 43       if (vis[i]==0)
 44     {
 45       tot++;
 46       prime[tot]=i;
 47       mu[i]=-1;
 48     }
 49       for (j=1;j<=tot;j++)
 50     {
 51       if (1ll*i*prime[j]>N) break;
 52       vis[i*prime[j]]=1;
 53       if (i%prime[j]==0)
 54         {
 55           mu[i*prime[j]]=0;
 56           break;
 57         }
 58       else mu[i*prime[j]]=-mu[i];
 59     }
 60     }
 61   for (i=1;i<=N;i++)
 62     F[i]=1;
 63   for (i=1;i<=N;i++)
 64     {
 65       for (j=1;j<=N;j++)
 66     {
 67       if (1ll*i*j>N) break;
 68       F[i*j]=(1ll*F[i*j]*pow(f[i],inv[i],mu[j]))%Mod;
 69     }
 70     }
 71 }
 72 int main()
 73 {int T;
 74   int i,j;
 75   cin>>T;
 76   for (i=1;i<=T;i++)
 77     {
 78       scanf("%d%d",&Q[i][0],&Q[i][1]);
 79       N=max(N,max(Q[i][0],Q[i][1]));
 80     }
 81   f[1]=1;f[0]=0;
 82   for (i=2;i<=N;i++)
 83     {
 84       f[i]=f[i-1]+f[i-2];
 85       if (f[i]>Mod) f[i]-=Mod;
 86     }
 87   for (i=1;i<=N;i++)
 88     inv[i]=qpow(f[i],Mod-2);
 89   pre();
 90   F[0]=1;
 91   for (i=1;i<=N;i++)
 92     F[i]=(1ll*F[i]*F[i-1])%Mod;
 93   for (j=1;j<=T;j++)
 94     {
 95       n=Q[j][0];m=Q[j][1];
 96       if (n>m) swap(n,m);
 97       int pos=0;
 98       ans=1;
 99       for (i=1;i<=n;i=pos+1)
100     {
101       if (n/(n/i)<m/(m/i)) pos=n/(n/i);
102       else pos=m/(m/i);
103       ans=(1ll*ans*qpow((1ll*F[pos]*qpow(F[i-1],Mod-2))%Mod,1ll*(n/i)*(m/i)%(Mod-1)))%Mod;
104     }
105       printf("%d
",ans);
106     }
107 }