BZOJ 4816 SDOI2017 数字表格 莫比乌斯反演

题目链接：http://www.lydsy.com/JudgeOnline/problem.php?id=4816

题意概述：

多组询问，给出一个n*m的表格，表格中的格子form[i][j]=f[gcd(i,j)]，f[i]表示fibonacci数列的第i项。

求表格中所有数字的乘积对 10^9+7 取模的结果。、

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

分析：

一眼就可以认出来这是莫比乌斯反演。关键在于式子的推倒。

 说一下推式子的要点，主要是能够灵活地利用反演公式，把式子踢来踢去！

啊这个博客实在是不好打公式......（至少我还没有研究出来），具体过程请看ATP的博客！(https://blog.csdn.net/fromatp/article/details/69944236)

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<cmath>
#include<queue>
#include<set>
#include<map>
#include<vector>
#include<cctype>
using namespace std;
const int mo=1000000007;
const int maxp=80005;
const int maxn=1000005;
typedef long long LL;

int T,N;
int f[maxn],pri[maxp],mu[maxn],pcnt,g[maxn],_g[maxn];
bool ntp[maxn];
struct data{ int n,m; }Q[1005];

void exgcd(LL a,LL b,LL &d,LL &x,LL &y){
    if(!b) d=a,x=1,y=0;
    else exgcd(b,a%b,d,y,x),y-=a/b*x;
}
int inv(int a){
    LL d,x,y;
    exgcd(a,mo,d,x,y);
    return (x%mo+mo)%mo;
}
int qkpow(int x,LL y)
{
    int re=1,tt=x;
    for(int i=0;(1ll<<i)<=y;i++){
        if((1ll<<i)&y) re=1ll*re*tt%mo;
        tt=1ll*tt*tt%mo;
    }
    return re;
}
void ready()
{
    ntp[0]=ntp[1]=1,mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!ntp[i]) pri[++pcnt]=i,mu[i]=-1;
        for(int j=1;j<=pcnt&&1ll*i*pri[j]<=N;j++){
            ntp[i*pri[j]]=1;
            if(i%pri[j]==0){ mu[i*pri[j]]=0; break; }
            mu[i*pri[j]]=-mu[i];
        }
    }
    f[0]=0,f[1]=1,g[0]=_g[0]=g[1]=_g[1]=1;
    for(int i=2;i<=N;i++){
        f[i]=(f[i-1]+f[i-2])%mo;
        g[i]=_g[i]=1;
    }
    for(int i=1;i<=N;i++)
    for(int j=1;1ll*i*j<=N;j++){
        if(!mu[j]) continue;
        if(mu[j]==1) g[i*j]=1ll*g[i*j]*f[i]%mo;
        else _g[i*j]=1ll*_g[i*j]*f[i]%mo;
    }
    for(int i=1;i<=N;i++)
        g[i]=1ll*g[i]*inv(_g[i])%mo;
    for(int i=2;i<=N;i++)
        g[i]=1ll*g[i]*g[i-1]%mo;
}
int solve(int n,int m)
{
    int re=1,last=0;
    for(int t=1;t<=n;t=last+1){
        last=min(n/(n/t),m/(m/t));
        re=1ll*re*qkpow(1ll*g[last]*inv(g[t-1])%mo,1ll*(n/t)*(m/t))%mo;
    }
    return re;
}
int main()
{
    scanf("%d",&T);
    for(int i=1;i<=T;i++){
        scanf("%d%d",&Q[i].n,&Q[i].m);
        if(Q[i].n>Q[i].m) swap(Q[i].n,Q[i].m);
        N=max(N,Q[i].n);
    }
    ready();
    for(int i=1;i<=T;i++)
        printf("%d
",solve(Q[i].n,Q[i].m));
    return 0;
}