【bzoj4816】[Sdoi2017]数字表格 莫比乌斯反演

题目描述

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取模。

输入

有多组测试数据。

第一个一个数T，表示数据组数。

接下来T行，每行两个数n,m

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

输出

输出T行，第i行的数是第i组数据的结果

样例输入

3
2 3
4 5
6 7

样例输出

1
6
960

---

题解

莫比乌斯反演

$prodlimits_{i=1}^nprodlimits_{j=1}^mf(gcd(i,j))\=prodlimits_{d=1}^{min(n,m)}(prodlimits_{i=1}^nprodlimits_{j=1}^m[gcd(i,j)=d]·f(d))\=prodlimits_{d=1}^{min(n,m)}f(d)^{sumlimits_{i=1}^nsumlimits_{j=1}^m[gcd(i,j)=d]}\=prodlimits_{d=1}^{min(n,m)}f(d)^{sumlimits_{i=1}^{lfloorfrac nd floor}sumlimits_{j=1}^{lfloor frac md floor}[gcd(i,j)=1]}\=prodlimits_{d=1}^{min(n,m)}f(d)^{sumlimits_{i=1}^{lfloorfrac nd floor}sumlimits_{j=1}^{lfloor frac md floor}sumlimits_{k|i&k&j}mu(d)}\=prodlimits_{d=1}^{min(n,m)}f(d)^{sumlimits_{k=1}^{min(lfloorfrac nd floor,lfloorfrac md floor)}mu(k)lfloorfrac n{dk} floorlfloorfrac m{dk} floor}$

到了这一步我们可以选择分块套分块，不过显然时间复杂度不足以应对多组询问。

继续令$D=dk$，可以得到：

$prodlimits_{d=1}^{min(n,m)}f(d)^{sumlimits_{k=1}^{min(lfloorfrac nd floor,lfloorfrac md floor)}mu(k)lfloorfrac n{dk} floorlfloorfrac m{dk} floor}\=prodlimits_{D=1}^{min(n,m)}(prodlimits_{d|D}f(d)^{mu(frac Dd)})^{lfloorfrac nD floorlfloorfrac mD floor}\=prodlimits_{D=1}^{min(n,m)}(t(D))^{lfloorfrac nD floorlfloorfrac mD floor}\(t(D)=prodlimits_{d|D}f(d)^{mu(frac Dd)})$

于是线性筛出$mu$，递推出f，预处理出f的逆元，进而使用$O(nln n)$的时间预处理出t数组。这里有一个小技巧：先枚举$frac Dd$，当它的$mu$值等于0时不作任何处理。亲测可以有效减少时间。

然后预处理出t数组之后分块处理询问即可。

时间复杂度为O(跑得过)$O(nln n+Tsqrt n)$，实际上跑了37s。

#include <cstdio>
#include <cstring>
#include <algorithm>
#define N 1000010
using namespace std;
typedef long long ll;
const int k = 1000000 , mod = 1000000007;
int mu[N] , prime[N] , tot;
ll f[N] , inv[N] , t[N] , sum[N];
bool np[N];
ll pow(ll x , ll y)
{
    ll ans = 1;
    while(y)
    {
        if(y & 1) ans = ans * x % mod;
        x = x * x % mod , y >>= 1;
    }
    return ans;
}
int main()
{
    int i , j , T , n , m;
    ll ans;
    sum[0] = t[1] = f[1] = inv[1] = mu[1] = 1;
    for(i = 2 ; i <= k ; i ++ )
    {
        t[i] = 1 , f[i] = (f[i - 1] + f[i - 2]) % mod , inv[i] = pow(f[i] , mod - 2) % mod;
        if(!np[i]) mu[i] = -1 , prime[++tot] = i;
        for(j = 1 ; j <= tot && i * prime[j] <= k ; j ++ )
        {
            np[i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                mu[i * prime[j]] = 0;
                break;
            }
            else mu[i * prime[j]] = -mu[i];
        }
    }
    for(j = 1 ; j <= k ; j ++ )
    {
        if(!mu[j]) continue;
        for(i = 1 ; i * j <= k ; i ++ )
        {
            if(mu[j] == 1) t[i * j] = t[i * j] * f[i] % mod;
            else t[i * j] = t[i * j] * inv[i] % mod;
        }
    }
    for(i = 1 ; i <= k ; i ++ ) sum[i] = sum[i - 1] * t[i] % mod;
    scanf("%d" , &T);
    while(T -- )
    {
        scanf("%d%d" , &n , &m) , ans = 1;
        for(i = 1 ; i <= n && i <= m ; i = j + 1)
            j = min(n / (n / i) , m / (m / i)) , ans = ans * pow(sum[j] * pow(sum[i - 1] , mod - 2) % mod , (ll)(n / i) * (m / i)) % mod;
        printf("%lld
" , ans);
    }
    return 0;
}