[莫比乌斯反演] BZOJ4816: [Sdoi2017]数字表格

挺简单的。。。

$$\prod_{i=1}^n \prod_{j=1}^m fib[gcd(i,j)]=\prod_d fib[d]^{\sum_i^n \sum_j^m [gcd(i,j)==d]}$$
上面这个指数是经典的反演的问题，然后就变成：
$$\prod_{d} fib[d]^{\sum_k \mu(k)\lfloor \frac n {dk}\rfloor \lfloor \frac m {dk}\rfloor}$$
套路变形一下
$$\prod_{T=1} (\prod_{d|T} fib[d]^{ \mu(\frac T d)})^{\lfloor \frac n T\rfloor \lfloor \frac m T\rfloor}$$
中间这个 $f(T)=\prod_{d|T} fib[d]^{\mu(\frac T d)}$，直接暴力 $nlogn$ 算。
$$\prod_{T=1} f(T)^{\lfloor \frac n T\rfloor \lfloor \frac m T\rfloor}$$
分块求就好了。

这题我跑的很慢，于是网上看了一下发现一个小trick:(其实是我太菜之前不知道)
求一个数列的每个数的逆元和前缀积的逆元，可以先算 $S_n^{-1}$，然后$O(n)$ 推。
$S_i^{-1}=S_{i+1}^{-1}*a_{i+1}\ \ a_i^{-1}=S_{i}^{-1}*S_{i-1}$

#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=1000005, N=1000000, MOD=1e9+7;
int p[maxn],mu[maxn];
bool vis[maxn];
void get_mu(){
    mu[1]=1;
    for(int i=2;i<=N;i++){
        if(!vis[i]) p[++p[0]]=i, mu[i]=-1;
        for(int j=1;j<=p[0]&&(LL)i*p[j]<=N;j++){
            vis[i*p[j]]=true;
            if(i%p[j]==0){ mu[i*p[j]]=0; break; }
            mu[i*p[j]]=-mu[i];
        }
    }
}
LL Pow(LL a,LL b){
    if(b==-1) return Pow(a,MOD-2);
    LL res=1;
    for(;b;a=a*a%MOD,b>>=1) if(b&1) res=(res*a)%MOD;
    return res; 
}
int _test,n,m,fib[maxn],fib_n[maxn],f[maxn],f_n[maxn];
LL ans;
int f_mul(int L,int R){ return (LL)f[R]*f_n[L-1]%MOD; }
int main(){
    freopen("bzoj4816.in","r",stdin);
    freopen("bzoj4816.out","w",stdout);
    get_mu();
    fib[0]=0; fib[1]=fib_n[1]=1; for(int i=2;i<=N;i++) fib[i]=(fib[i-1]+fib[i-2])%MOD, fib_n[i]=Pow(fib[i],MOD-2);
    for(int i=1;i<=N;i++) f[i]=1;
    for(int i=1;i<=N;i++)
     for(int j=1;(LL)i*j<=N;j++){
        int t=(mu[j]==0?1:(mu[j]==1?fib[i]:fib_n[i]));
        f[i*j]=((LL)f[i*j]*t)%MOD;
     }
    f[0]=f_n[0]=1; for(int i=1;i<=N;i++) f[i]=((LL)f[i]*f[i-1])%MOD, f_n[i]=Pow(f[i],MOD-2);
    scanf("%d",&_test);
    while(_test--){
        scanf("%d%d",&n,&m); int _min=min(n,m);
        ans=1;
        for(int i=1,nxt;i<=_min;i=nxt+1){
            nxt=min(n/(n/i),m/(m/i));
            ans=(ans*Pow(f_mul(i,nxt),(LL)(n/i)*(m/i)))%MOD;
        }
        printf("%lld\n",ans);
    }
    return 0;
}