【刷题记录】SDOI2017数字表格

https://www.luogu.org/problemnew/show/P3704
题目描述
Doris刚刚学习了fibonacci数列。用f[i]表示数列的第i项，那么
f[0]=1,f[1]=1,
f[n]=f[n-1]+f[n-2]
Doris用老师的超级计算机生成了一个n×m的表格，
第i行第j列的格子中的数是f[gcd(i,j)]，其中gcd(i,j)表示i,j的最大公约数。
Doris的表格中共有n×m个数，她想知道这些数的乘积是多少。
答案对10^9+7取模。
莫比乌斯反演

直接上表达式（F(i)表示斐波那契数列第I项）
$ans =\prod_{i=1}^{m}\prod_{j=1}^nF((i,j))$
$=\prod_{i=1}^{min(m,n)}{F(i)}^{\sum_{(j,k)=i}1}$
分子是莫比乌斯反演中常见套路
$=\prod_{i=1}^{min(m,n)}{F(i)}^{\sum_{i|d}μ(\frac d i)[ \frac n d ][\frac m d]}$
$lnans=\sum_{d=1}^{min(m,n)}[\frac n d][\frac m d]\sum_{i|d}lnF(i)μ(\frac d i)$
$ans = \prod_{d=1}^{min(n,m)}(\prod_{i|d}F(i)^{μ(\frac d i )})^{[ \frac n d ][\frac m d]}$
中间括号可以前缀积（理论上也可以线段树）但是注意使用前缀积需要处理逆元
再使用快速乘
这样总时间复杂度
$O(nlogn+Tlogn\sqrt n)$
代码

#include <iostream>
#include <cstdio>
#include <algorithm>
using namespace std;
#define mod 1000000007
#define N 1000000
#define ll long long 
int f[N+6],pre[N+6],prime[N+6],mu[N+6],tot;
bool notprime[N+6];
int fastpow(int x,int n)
{
    if(!n)return 1;
    int ans = fastpow(x,n/2);
    ans = ((ll)ans *ans)%mod;
    if(n&1)
    ans = ((ll)ans *x)%mod;
    return ans;
}
int exgcd(ll &x,ll &y,int a, int b)
{
    if(b==0)
    {
        x = 1;
        y = 0;
        return a;
    }
    int ans = exgcd(x,y,b,a%b);
    ll t = x;
    x = y;
    y = t- (a/b)*y;
    return ans;
}
int inv(int m)
{
    ll x,y;
    exgcd(x,y,m,mod);
    return (int)(((x%mod)+mod)%mod);
}
void init()
{
    pre[0]=pre[1]=f[1]=notprime[1]=mu[1]=1;
    for(int i = 2; i <= N; i ++)
    {
        pre[i]=1;
        f[i]=((ll)f[i-1]+f[i-2])%mod;
        if(!notprime[i])
        {
            prime[++tot]=i;
            mu[i]=-1;
        }
        for(int j = 1 ; j <= tot && i*prime[j] <= N  ; j ++)
        {
            notprime[i*prime[j]]=1;
            if(i%prime[j]!=0)mu[i*prime[j]]=-mu[i];
            else 
            {
                break;
            }
        }
    }
    for(int i = 1; i <= N ; i ++)
    for(int j = i; j <= N ; j += i)
    {   
        int t = f[i];
        if(mu[j/i]==-1)
        {
            t = inv(t);
        }
        else if(mu[j/i]==0)
        {
            t = 1;
        }
        pre[j]=((ll)pre[j]*t)%mod;
    }
    for(int i = 1; i <= N; i ++)
    {
        pre[i]=((ll)pre[i]*pre[i-1])%mod;
    }
}
int main()
{
    int t;
    init();
    scanf("%d",&t);
    while(t--)
    {
        int ans = 1, n,m;
        scanf("%d%d",&n,&m);
        if(n>m)swap(m,n);
        for(int i = 1,last; i <= n; i = last +1 )
        {
            last = min(n/(n/i),m/(m/i));
            ans = ((ll)ans* fastpow((int)(((ll)pre[last]*inv(pre[i-1]))%mod),(n/i)*(m/i)))%mod;
            
        }
        printf("%d\n",ans);
    }
}