[莫比乌斯反演] BZOJ 4816 [Sdoi2017]数字表格

如果我们有$fib_n=\prod_{d|n} f_d$
那么 就有

$$\begin{eqnarray} \prod_{i=1}^n\prod_{j=1}^m fib_{gcd(i,j)} &=&\prod_{i=1}^n\prod_{j=1}^m \prod_{d|i,d|j}f_d &=&\prod_{d=1}^{min(n,m)}f_d^{\lfloor {n\over d}\rfloor \lfloor {m\over d}\rfloor} \end{eqnarray}$$
这个分块 指数一样的一起快速幂就好了
怎么求 $f$ 本来以为要什么高超的反演技巧 实际上直接递推
$$f_n={fib_n \over \prod_{d|n,d\ne n} f_d}$$

#include<cstdio>
#include<cstdlib>
#include<algorithm>
using namespace std;
typedef long long ll;

const int N=1000005;
const int P=1e9+7;
const int G=3;

inline ll Pow(ll a,int b){
  ll ret=1;
  for (;b;b>>=1,a=a*a%P)
    if (b&1)
      ret=ret*a%P;
  return ret;
}
inline ll Inv(ll a){
  return Pow(a,P-2);
}

ll f[N],inv[N];

inline void Pre(int n){
  f[0]=0; f[1]=1;
  for (int i=2;i<=n;i++) f[i]=(f[i-1]+f[i-2])%P;
  for (int i=1;i<=n;i++){
    ll inv=Inv(f[i]);
    for (int j=i+i;j<=n;j+=i)
      f[j]=f[j]*inv%P;
  }
  f[0]=inv[0]=1;
  for (int i=1;i<=n;i++)
    f[i]=f[i]*f[i-1]%P,inv[i]=Inv(f[i]);
}

inline int Solve(int n,int m){
  if (n>m) swap(n,m);
  ll ret=1;
  for (int i=1,j;i<=n;i=j+1){
    j=min(n/(n/i),m/(m/i));
    ret=ret*Pow(f[j]*inv[i-1]%P,(ll)(n/i)*(m/i)%(P-1))%P;
  }
  return ret;
}

#define read(x) scanf("%d",&(x))

int main(){
  int T,n,m;
  freopen("t.in","r",stdin);
  freopen("t.out","w",stdout);
  Pre(1e6);
  read(T);
  while (T--){
    read(n); read(m);
    printf("%d\n",Solve(n,m));
  }
  return 0;
}