51Nod1965 奇怪的式子

题意：求 $\prod_{i = 1} ^ n \sigma_0(i) ^ {i + \mu(i)}$。$n\leq 10 ^ {11}$，对 $10 ^ {12} + 39$ 取模，$T(T\leq 3)$ 组数据，15s。

前置知识，即使会也可以看看里面的一些定义，下文会用到~

离谱模数 + 奇怪卡常题，不过 Min_25 的应用是非常不错的。

首先我们容易发现 $i + \mu(i)$ 似乎没有任何性质，考虑拆开计算。那么分别计算 $\sigma_0(i) ^ i$ 和 $\sigma_0(i) ^ {\mu(i)}$。

1. 正常的 Min_25 筛

---

第一部分

直接考虑 Min_25 筛，考虑计算 $\prod\sigma_0(i) ^ i$。假设我们已经知道了一些 $x$ 的 $\prod \sigma_0(x) ^ x$，考虑怎样在这些数上所有乘 $p$（$p$ 是一个质数，$s$ 和 $p$ 都互质）。直接推式子：

$$\begin{aligned} &\prod \sigma_0(p\times x) ^ {p\times x}\\ =& \prod \sigma_0(p) ^ {p\times x} \times \prod \sigma_0(x) ^ {p\times x}\\ =& (\sigma_0(p) ^ p) ^ {\sum x} \times (\prod \sigma_0(x) ^ x) ^ p \end{aligned}$$

容易发现我们还需要计算 $\sum x$，这个在 Min_25 递归过程中很容易计算。$p$ 也可以换成 $p ^ k$，和上面的推导是相同的，可以直接使用。

第二部分

下面考虑 $\prod\sigma_0(i) ^ {\mu(i)}$。首先一个比较重要的性质为我们只需要管 $p_1p_2 \dots p_k$ 类似的数，因为 $\mu(i)$ 为 0 的话就没有贡献。而写成这样，我们又可以表示 $\sigma_0(i)$ 为 $2 ^ {d(i)}$，其中 $d(i)$ 表示 $i$ 的质因子个数。那么只需要计算 $\sum d(i)\mu(i)$ 即可。

还是考虑 Min_25 筛的递归过程，我们只需要枚举次幂 $k = 1$ 即可。继续计算时，还是会遇到类似合并两个答案的问题。

假设前面我们已经算好了一些 $\sum d(x) \mu(x)$，现在考虑对这些所有数乘上一个 $p$。还是直接推式子：

$$\begin{aligned} &\sum d(x\times p)\mu(x\times p)\\ =& \sum (d(x) + 1)(-\mu(x))\\ =& -\sum d(x)\mu(x) - \sum \mu(x) \end{aligned}$$

递归的过程需要同时计算 $\mu(x)$，这个是容易的。

---

那么我们就得到了一个使用递归 Min_25 筛的做法。

递归代码

const LL Mod1 = 1000000000039LL, Mod2 = 1000000000038LL;
struct LPow {
    static const int B = 1 << 20;
    LL sk[B], bk[B];
    void init(int x)
    {
        sk[0] = bk[0] = 1;
        for (int i = 1; i < B; ++ i) sk[i] = sk[i - 1] * x % Mod1;
        bk[1] = sk[B - 1] * x % Mod1;
        for (int i = 2; i < B; ++ i) bk[i] = (s128) bk[i - 1] * bk[1] % Mod1;
    }
    
    LL lpow(LL n) { return n %= Mod2, (s128) bk[n >> 20] * sk[n & 1048575] % Mod1; }
} lp[25];
bool mem2;

LL qpow(LL a, LL k = Mod1 - 2)
{
    LL res = 1;
    for (; k; k >>= 1, a = (s128) a * a % Mod1)
        if (k & 1) res = (s128) res * a % Mod1;
    return res;
}

void sieve()
{
    for (int i = 2; i < N; ++ i)
    {
        if (!st[i]) prime[++ cnt] = i;
        for (int j = 1; j <= cnt && i * prime[j] < N; ++ j)
        {
            st[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
}

inline int getid(LL x) { return x <= sq ? id1[x] : id2[n / x]; }

void solveg1()
{
    for (int i = 1; i <= tot; ++ i) g[i] = a[i] - 1;
    for (int i = 1; i <= cnt; ++ i) sp[i] = i;
    for (int i = 1; i <= cnt; ++ i)
    {
        LL le = (LL) prime[i] * prime[i];
        for (int j = 1; j <= tot && a[j] >= le; ++ j)
            g[j] -= g[getid(a[j] / prime[i])] - sp[i - 1];
    }
    // std::cout << g[1] << '\n';
}

void solveg2()
{
    for (int i = 1; i <= tot; ++ i) g[i] = ((s128) a[i] * (a[i] + 1) / 2 - 1) % Mod2;
    for (int i = 1; i <= cnt; ++ i) sp[i] = sp[i - 1] + prime[i];
    for (int i = 1; i <= cnt; ++ i)
    {
        LL le = (LL) prime[i] * prime[i];
        for (int j = 1; j <= tot && a[j] >= le; ++ j)
            g[j] = (g[j] + (s128) (g[getid(a[j] / prime[i])] - sp[i - 1] + Mod2)
                * (Mod2 - prime[i])) % Mod2;
    }
    // std::cout << g[1] << '\n';
}

PLL S1(LL n, int i)
{
    if (prime[i] >= n) return {};
    LL res = -g[getid(n)] + sp[i], tot = res;
    for (int nw = i + 1; nw <= cnt && (LL) prime[nw] * prime[nw] <= n; ++ nw)
    {
        PLL frm = S1(n / prime[nw], nw);
        res -= frm.first + frm.second, tot -= frm.second;
    }
    return {res, tot};
}

PLL S2(LL n, int i)
{
    if (prime[i] >= n) return {1, 0};
    LL res = lp[2].lpow(g[getid(n)] - sp[i] + Mod2), tot = g[getid(n)] - sp[i] + Mod2;
    for (int nw = i + 1; nw <= cnt && (LL) prime[nw] * prime[nw] <= n; ++ nw)
    {
        LL cur = prime[nw];
        for (int k = 1; cur <= n; cur *= prime[nw], ++ k)
        {
            PLL frm = S2(n / cur, nw);
            if (k != 1) frm.second ++;
            LL mul, t = (s128) cur * frm.second % Mod2;
            if (k + 1 < 25) mul = lp[k + 1].lpow(t);
            else mul = qpow(k + 1, t);
            res = (s128) res * mul % Mod1 * qpow(frm.first, cur) % Mod1;
            tot = (tot + (s128) frm.second * cur) % Mod2;
        }
    }
    return {res, tot};
}

void work()
{
    std::cin >> n;
    sq = std::sqrt(n), tot = 0;
    for (LL l = 1, r, t; l <= n; l = r + 1)
    {
        t = n / l, r = n / t;
        a[++ tot] = t;
        if (t <= sq) id1[t] = tot;
        else id2[r] = tot;
    }
    solveg1();
    auto res1 = S1(n, 0);
    // std::cout << (res1.first + Mod2 * 10) % Mod2 << '\n';
    solveg2();
    auto res2 = S2(n, 0);
    LL res = (s128) qpow(2, res1.first + Mod2 * 10) * res2.first % Mod1;
    std::cout << res << '\n';
}

2. 迭代的 Min_25 筛

如果认真观察一些 Min_25 的题目的话，会发现一个事实：如果带递归的 Min_25，一般只能跑 $10 ^ {10}$，但是如果只有第一部分迭代的话，一般可以跑 $10 ^ {11}$，不过时限略大。本题中，上面给出的代码在单组 $10 ^ {11}$ 的情况下需要跑 1min 的时间，很明显的印证了上面所说到的这一点。（递归中还带快速幂，显然无法接受）

考虑使用迭代的方法计算这两部分。可能就和正常理解的 Min_25 不太相同。

---

第二部分

先看第二部分罢。

注意到这个部分递归的时候，$p$ 的次幂都是 $k = 1$，考虑从这里入手。

还是类似于 Min_25 筛的定义，设 $S(n, i)$ 为所有满足 $\mathrm{lpf}(x) > p_i\lor i\in P$ 的数的 $\sum d(x)\mu(x)$ 的和，其中 $\mathrm{lpf}(x)$ 表示 $x$ 的最小质因子，$p_i$ 是第 $i$ 个质数，$P$ 是质数集合。

考虑直接迭代求 $G(n, i)$。和 Min_25 筛第一部分不同的地方在于，我们现在求的是 $G(n, 0)$ 而不是 $G(n, k)$（保证 $p_k > \sqrt n$），所以我们应该倒序枚举 $i$。

注意到我们迭代的时候也用到了 $\sum \mu(i)$，设 $H(n, k)$ 为满足条件的 $\mu(x)$ 的和。

大力推式子，可以得到：

$$\begin{aligned} G(n, i - 1) =& G(n, i) - (G(\dfrac n{p_i}, i) + i) - (H(\dfrac n{p_i}, i) + i)\\ H(n, i - 1) =& H(n, i) - (H(\dfrac n{p_i}, i) + i) \end{aligned}$$

$+i$ 的原因是我们要把 $p_{1\dots i}$ 的贡献除去，每个的 $d(i)\mu(i)$ 和 $\mu(i)$ 都是 -1，所以减掉贡献应该 $+i$。

那么按照这个式子计算，最后 $G(n, 0)$ 就是答案。这个计算过程类似于 Min_25 第一步，并不需要滚动数组。

第一部分

那么再来看第一部分。这个就和上面那一个优化思路就不同了，因为枚举的 $p$ 的次幂不一定是 1 了。

从朴素的方向出发，我们考虑计算每一个质因子的贡献。大力推式子：

$$\begin{aligned} \prod_{i = 1} ^ n \sigma_0(i) ^ i =& \prod_{p\in P} \prod_{k} (k + 1) ^ {\sum_{i = 1} ^ n [p ^ k \mid i\land p ^ {k + 1}\nmid i] i} \end{aligned}$$

指数上的一坨我们展开一下：

$$\begin{aligned} &\sum_{i = 1} ^ n [p ^ k \mid i \land p ^ {k + 1} \nmid i]\\ =& \sum_{i = 1} ^ n [p ^ k | i]i - \sum_{i = 1} ^ n [p ^ {k + 1} | i] i\\ =& p ^ k \sum_{i = 1} ^ {\frac n{p ^ k}}i - p ^ {k + 1} \sum_{i = 1} ^ {\frac n{p ^ {k + 1}}} i\\ =& p ^ k \mathrm{sum}(\dfrac n{p ^ k}) - p ^ {k + 1} \mathrm{sum} (\dfrac n{p ^ {k + 1}}) \end{aligned}$$

其中 $\mathrm{sum}(x)$ 表示 $\sum_{i = 1} ^ x i$，为 $\dfrac {x\times (x + 1)}2$。

这么多 $(p, k)$ 并不好计算，但是注意到 $p > \sqrt n$ 的时候，$k = 1$，似乎是一个不错的办法。

首先计算 $p\leq \sqrt n$。直接暴力计算即可，时间复杂度并不高，反正很快。

下面来看 $p > \sqrt n$ 的情况。注意到 $k = 1$ 的情况是很好表示的，我们可以马上写出式子：

$$\prod_{p\in P\land p > \sqrt n} 2 ^ {p \times \mathrm{sum}(\lfloor \frac np\rfloor )}$$

注意到我们需要对 $\dfrac np$ 进行整除分块，然后需要求出区间素数的和。这个是直接使用 Min_25 第一部分就可以很快得到的（类似 区间素数问题），那么我们也可以只通过迭代计算第一部分的答案了。

---

通过对 Min_25 的改造，对式子的推导，我们最终还是通过简单的 Min_25 第一部分的迭代就可以计算这个答案了。时间复杂度不太会算，据说都是 $O(n ^ {\frac 34})$ 的。新算法能在 3s 左右得到 $10 ^ {11}$ 的答案的，可见还是比递归写法快不少的。

代码实现方面，注意是 $\bmod P$ 还是 $\bmod \varphi(P)$ 是需要注意的。好像不能使用 __int128_t，那就写 long double 乘法罢。

const LL Mod1 = 1000000000039LL, Mod2 = 1000000000038LL;
LL mul1(LL a, LL b)
{
    long double x = (long double) a * b;
    return (a * b - (LL) (x / Mod1) * Mod1 + Mod1);
}

LL mul2(LL a, LL b)
{
    long double x = (long double) a * b;
    return (a * b - (LL) (x / Mod2) * Mod2 + Mod2);
}

LL qpow(LL a, LL k = Mod1 - 2)
{
    LL res = 1;
    for (; k; k >>= 1, a = mul1(a, a))
        if (k & 1) res = mul1(res, a);
    return res;
}

void sieve()
{
    for (int i = 2; i < N; ++ i)
    {
        if (!st[i]) prime[++ cnt] = i;
        for (int j = 1; j <= cnt && i * prime[j] < N; ++ j)
        {
            st[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
}

inline int getid(LL x) { return x <= sq ? id1[x] : id2[n / x]; }
inline LL sum1(LL x) { return x & 1 ? mul2(x, (x + 1) / 2) : mul2(x / 2, x + 1); }
LL& adj(LL &x) { return x += x >> 63 & Mod2; }

void solveg1() // sigma(i) * i
{
    for (int i = 1; i <= tot; ++ i) g[i] = sum1(a[i]) - 1;
    for (int i = 1; i <= tot; ++ i) h[i] = a[i] - 1;
    for (int i = 1; i <= cnt; ++ i) sp[i] = sp[i - 1] + prime[i];
    for (int i = 1; i <= cnt; ++ i) sph[i] = i;
    for (int i = 1; i <= cnt; ++ i)
    {
        LL le = (LL) prime[i] * prime[i], tmp;
        for (int j = 1, id; j <= tot && a[j] >= le; ++ j)
            g[j] = (g[j] + (Mod2 - adj(tmp = g[id = getid(a[j] / prime[i])] - sp[i - 1]))
                 * prime[i]) % Mod2,
            h[j] -= h[id] - sph[i - 1];
    }
    // std::cout << g[1] << '\n';
}

void solveg2()
{
    for (int i = 1; i <= tot; ++ i) h[i] = Mod2 - h[i], g[i] = h[i];
    // std::cout << g[1] << '\n';
    for (int i = cnt; i; -- i)
    {
        LL le = (LL) prime[i] * prime[i];
        for (int j = 1; j <= tot && a[j] >= le; ++ j)
        {
            int nid = getid(a[j] / prime[i]);
            adj(adj(h[j] -= h[nid]) -= i);
            g[j] = (g[j] - h[nid] - i - g[nid] - i + 4 * Mod2) % Mod2;
        }
    }
    // std::cout << g[1] << '\n';
}

void work()
{
    std::cin >> n;
    sq = std::sqrt(n) + 1, tot = 0;
    for (LL l = 1, r, t; l <= n; l = r + 1)
    {
        t = n / l, r = n / t;
        a[++ tot] = t;
        if (t <= sq) id1[t] = tot;
        else id2[r] = tot;
    }
    
    solveg1();
    int lim = 1;
    LL res = 1;
    for (; a[lim + 1] >= sq; ++ lim) ;
    for (int i = 1; i <= cnt; ++ i)
    {
        if (prime[i] > a[lim]) break;
        LL cur = prime[i];
        for (int k = 1; cur <= n; k ++, cur *= prime[i])
        {
            LL mul = (mul2(cur, sum1(n / cur))
                + mul2((Mod2 - cur) * prime[i] % Mod2, sum1(n / cur / prime[i]))) % Mod2;
            res = mul1(res, qpow(k + 1, mul));
        }
    }
    LL mul = 0;
    for (LL l = 1, r, t; l <= n; l = r + 1)
    {
        t = n / l, r = n / t;
        if (getid(r) > lim) continue;
        mul = (mul + mul2(g[getid(r)] - g[std::min(lim, getid(l - 1))] + Mod2, sum1(t))) % Mod2;
    }
    res = mul1(res, qpow(2, mul));

    solveg2();
    res = mul1(res, qpow(2, g[1]));
    // std::cout << g[1] << '\n';
    
    std::cout << res % Mod1 << '\n';
}