bzoj4816 [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

正解：莫比乌斯函数。

水水的一道题，不过卡常数。。推导一波吧。。

$Ans=prod_{i=1}^{n}prod_{j=1}^{m}f(gcd(i,j))$

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

直接跳过吧。。因为中间的都是老套路了。。

$Ans=prod_{d=1}^{min(n,m)}f(d)^{sum_{p=1}^{min(left lfloor frac{n}{d} ight floor,left lfloor frac{m}{d} ight floor)} mu(p)left lfloor frac{n}{dp} ight floorleft lfloor frac{m}{dp} ight floor}$

然后好像没办法往下化简了，不过这个式子用数论分块就能过了。。因为复杂度不是满的，好像是$O(Tn^{frac{3}{4}})$吧。。

反正这就能$AC$了。。

//It is made by wfj_2048~
#include <algorithm>
#include <iostream>
#include <complex>
#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <vector>
#include <cmath>
#include <queue>
#include <stack>
#include <map>
#include <set>
#define rhl (1000000007)
#define inf (1<<30)
#define N (1000010)
#define il inline
#define RG register
#define ll long long
#define File(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)

using namespace std;

int vis[N],inv[N],mu[N],prime[N],n,m,cnt,pos1,pos2;
ll f[N],ans;

il int gi(){
    RG int x=0,q=1; RG char ch=getchar();
    while ((ch<'0' || ch>'9') && ch!='-') ch=getchar();
    if (ch=='-') q=-1,ch=getchar();
    while (ch>='0' && ch<='9') x=x*10+ch-48,ch=getchar();
    return q*x;
}

il ll qpow(RG ll a,RG ll b){
    RG ll ans=1;
    while (b){
    if (b&1) ans=ans*a%rhl;
    a=a*a%rhl,b>>=1;
    }
    return ans;
}

il void pre(){
    f[1]=vis[1]=mu[1]=1;
    for (RG int i=2;i<N;++i){
    if (!vis[i]) prime[++cnt]=i,mu[i]=-1;
    for (RG int j=1,k;j<=cnt;++j){
        k=i*prime[j]; if (k>=N) break; vis[k]=1;
        if (i%prime[j]) mu[k]=-mu[i]; else break;
    }
    f[i]=f[i-1]+f[i-2]; if (f[i]>=rhl) f[i]-=rhl;
    }
    inv[0]=inv[1]=1,f[0]=1;
    for (RG int i=2;i<N;++i)
    mu[i]+=mu[i-1],f[i]*=f[i-1],f[i]%=rhl,inv[i]=qpow(f[i],rhl-2);
    return;
}

il void work(){
    n=gi(),m=gi(),ans=1; if (n>m) swap(n,m);
    for (RG int i=1;i<=n;i=pos1+1){
    pos1=min(n/(n/i),m/(m/i)); RG ll res=0;
    for (RG int j=1;j<=n/i;j=pos2+1){
        pos2=min(n/i/(n/i/j),m/i/(m/i/j));
        res+=(ll)(mu[pos2]-mu[j-1])*(n/i/j)*(m/i/j);
    }
    ans*=qpow(f[pos1]*(ll)inv[i-1]%rhl,res),ans%=rhl;
    }
    printf("%lld
",ans); return;
}

int main(){
    File("product");
    pre(); RG int T=gi();
    while (T--) work();
    return 0;
}