洛谷 P5221 Product

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

\[\prod_{i=1}^N\prod_{j=1}^N\frac{lcm(i,j)}{gcd(i,j)}\ (\bmod\ 104857601) \]

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

\(1 \leq N \leq 1000000\)

---

把式子写成：

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

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

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

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

\[\begin{aligned}g(x,n)&=\sum_{i=1}^n\sum_{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)\) ，就有：

\[\begin{aligned} g(1,n)&=\sum_{i=1}^{n}\sum_{j=1}^{n}\sum_{d|gcd(i,j)}\mu(d)\\&=\sum_{d=1}^n\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}\mu(d)\\&=\sum_{d=1}^n\lfloor\frac{n}{d}\rfloor\lfloor\frac{n}{d}\rfloor\mu(d)\end{aligned}\]

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

然后现在我们要求：

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

对 \(\frac{N}{i}\) 也整除分块，这样就可以做到 \(O(N+\sqrt N\log 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;
}