洛谷 P5221 Product

({ m CYJian})最近闲的玩起了(gcd)。。他想到了一个非常简单而有意思的式子：

[prod_{i=1}^Nprod_{j=1}^Nfrac{lcm(i,j)}{gcd(i,j)} (mod 104857601) ]

({ m CYJian})已经算出来这个式子的值了。现在请你帮他验算一下吧。({ m CYJian})只给你(0.2s)的时间哦。

(1 leq N leq 1000000)

---

把式子写成：

[=prod_{i=1}^Nprod_{j=1}^Nfrac{ij}{gcd^2(i,j)}=frac{prod_{i=1}^Nprod_{j=1}^Nij}{(prod_{i=1}^Nprod_{j=1}^Ngcd(i,j))^2} ]

上面的部分就是 ((N!)^{2N}) ，对于下面 (prod_{i=1}^Nprod_{j=1}^Ngcd(i,j)) 部分，考虑枚举 (gcd) 是多少，设 (g(x,n)) 表示 (sum_{i=1}^nsum_{j=1}^n[gcd(i,j)==x]) ，那么这个式子就变成了：

[prod_{i=1}^Ni^{g(i,N)} ]

考虑怎么求所有的 (g(i,N)) ，我们写出式子：

[egin{aligned}g(x,n)&=sum_{i=1}^nsum_{j=1}^n[gcd(i,j)==x]\&=sum_{i=1}^{frac{n}{x}}sum_{j=1}^{frac{n}{x}}[gcd(i,j)==1]\&=g(1,frac{n}{x})end{aligned} ]

然后考虑求 (g(1,n)) ，就有：

[egin{aligned} g(1,n)&=sum_{i=1}^{n}sum_{j=1}^{n}sum_{d|gcd(i,j)}mu(d)\&=sum_{d=1}^nsum_{i=1}^{frac{n}{d}}sum_{j=1}^{frac{n}{d}}mu(d)\&=sum_{d=1}^nlfloorfrac{n}{d} floorlfloorfrac{n}{d} floormu(d)end{aligned}]

这个可以整除分块 (sqrt{n}) 完成。

然后现在我们要求：

[prod_{i=1}^Ni^{g(1,frac{N}{i})} ]

对 (frac{N}{i}) 也整除分块，这样就可以做到 (O(N+sqrt Nlog N)) 的复杂度处理出来了。

Code

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cstring>
const int N = 1e6;
const int M = 1e5;
const int p = 104857601;
using namespace std;
int n,mul[N + 5],prime[M + 5],pcnt,ans,s;
bool vis[N + 5];
int mypow(int a,long long x){int s = 1;for (;x;x & 1 ? s = 1ll * a * s % p : 0,a = 1ll * a * a % p,x >>= 1);return s;}
long long calc(int n)
{
    long long ans = 0;
    for (int l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        ans += 1ll * (n / l) * (n / l) * (mul[r] - mul[l - 1]);
    }
    return ans;
}
void prework()
{
    vis[1] = 1;mul[1] = 1;
    for (int i = 2;i <= n;i++)
    {
        if (!vis[i])
        {
            mul[i] = -1;
            prime[++pcnt] = i;
        }
        for (int j = 1;j <= pcnt && prime[j] * i <= n;j++)
        {
            vis[prime[j] * i] = 1;
            mul[prime[j] * i] = -mul[i];
            if (i % prime[j] == 0)
            {
                mul[prime[j] * i] = 0;
                break;
            }
        }
    }
    for (int i = 1;i <= n;i++)
        mul[i] += mul[i - 1];
}
int main()
{
    scanf("%d",&n);
    prework();
    ans = 1;
    for (int i = 1;i <= n;i++)
        ans = 1ll * ans * i % p;
    ans = mypow(ans,n);
    ans = 1ll * ans * ans % p;
    s = 1;
    for (int l = 1,r;l <= n;l = r + 1)
    {
        r = n / (n / l);
        int tmp = 1;
        for (int i = l;i <= r;i++)
            tmp = 1ll * tmp * i % p;
        s = 1ll * s * mypow(tmp,calc(n / l)) % p;
    }
    s = 1ll * s * s % p;
    ans = 1ll * ans * mypow(s,p - 2) % p;
    cout<<(ans + p) % p<<endl;
    return 0;
}