51nod 1575 Gcd and Lcm

题目：https://www.51nod.com/onlineJudge/questionCode.html#!problemId=1575

题意：
有 $T$ 组数据，每组数据给出一个正整数 $n$ ，求 $\sum_{i=1}^{n}{\sum_{j=1}^{i}{\sum_{k=1}^{i}{\text{lcm}(\gcd(i,j),\gcd(i,k))}}} \bmod 2^{32}$ 。
$T \leq 10, n \leq 10^9$ 。

题解：
观察到 $j$ 和 $k$ 是无关的、可轮换的，不妨尝试将 $j$ 和 $k$ 化简成相同的形式。
注意到 $\sum_{j\leq i}{[\gcd(i,j)=d]}=\sum_{\frac{j}{d}\leq\frac{i}{d}}{[\gcd(\frac{i}{d},\frac{j}{d})=1]}=\varphi(\frac{i}{d})$ ，所以我们可以对原式进行初步化简：

$$\sum_{i=1}^{n}{\sum_{j=1}^{i}\sum_{k=1}^{i}{\text{lcm}(\gcd(i,j),\gcd(i,k))}}\\=\sum_{i=1}^{n}{\sum_{j=1}^{i}\sum_{k=1}^{i}{\sum_{d_1}{\sum_{d_2}{[\gcd(i,j)=d_1][\gcd(i,k)=d_2]\text{lcm}(d_1,d_2)}}}}\\=\sum_{i=1}^{n}{\sum_{d_1|i}{\sum_{d_2|i}{\text{lcm}(d_1,d_2)\sum_{j=1}^{i}{[\gcd(i,j)=d_1]}\sum_{k=1}^{i}{[\gcd(i,k)=d_2]}}}}\\=\sum_{i=1}^{n}{\sum_{d_1|i}{\sum_{d_2|i}{\text{lcm}(d_1,d_2)\varphi(\frac{i}{d_1})\varphi(\frac{i}{d_2})}}}$$

为了消除 $d_1$ 和 $d_2$ 相互的限制，我们需要将 $\text{lcm}(d_1,d_2)$ 拆开。
对于一般情况，不难得到 $\text{lcm}(n,m)=\sum_{d}{[\gcd(\frac{n}{d},\frac{m}{d})=1]\frac{n}{d}\frac{m}{d}d}=\sum_{d|n,d|m}{\sum_{d'|\frac{n}{d},d'|\frac{m}{d}}{\mu(d')d'^2\frac{n}{dd'}\frac{m}{dd'}d}}$ ，所以我们可以尝试进行代换：

$$\sum_{i=1}^{n}{\sum_{d_1|i}{\sum_{d_2|i}{\text{lcm}(d_1,d_2)\varphi(\frac{i}{d_1})\varphi(\frac{i}{d_2})}}}\\=\sum_{i=1}^{n}{\sum_{d_1|i}{\sum_{d_2|i}{\sum_{d|d_1,d|d_2}{\sum_{d'|\frac{d_1}{d},d'|\frac{d_2}{d}}{\mu(d')d'^2\frac{d_1}{dd'}\frac{d_2}{dd'}d}}\varphi(\frac{i}{\frac{d_1}{dd'}dd'})\varphi(\frac{i}{\frac{d_2}{dd'}dd'})}}}\\=\sum_{i=1}^{n}{\sum_{d|i}{\sum_{d'|\frac{i}{d}}{d\mu(d')d'^2\sum_{\frac{d_1}{dd'}|\frac{i}{dd'}}{\frac{d_1}{dd'}\varphi(\frac{\frac{i}{dd'}}{\frac{d_1}{dd'}})}\sum_{\frac{d_2}{dd'}|\frac{i}{dd'}}{\frac{d_2}{dd'}\varphi(\frac{\frac{i}{dd'}}{\frac{d_2}{dd'}})}}}}\\=\sum_{i=1}^{n}{\sum_{d|i}{\sum_{d'|\frac{i}{d}}{d\mu(d')d'^2(\sum_{\frac{d''}{dd'}|\frac{i}{dd'}}{\frac{d''}{dd'}\varphi(\frac{\frac{i}{dd'}}{\frac{d''}{dd'}})})^2}}}$$

在上面的最后一步，终于将 $j$ 和 $k$ 对应的 $d_1$ 和 $d_2$ 表示成类似的形式，并且上式可以表示成：

$$\sum_{i=1}^{n}{\sum_{d|i}{\sum_{d'|\frac{i}{d}}{d\mu(d')d'^2(\sum_{\frac{d''}{dd'}|\frac{i}{dd'}}{\frac{d''}{dd'}\varphi(\frac{\frac{i}{dd'}}{\frac{d''}{dd'}})})^2}}}\\=\sum_{i=1}^{n}{\sum_{xyz=i}{x\mu(y)y^2(\sum_{d|z}{d\varphi(\frac{z}{d})})^2}}\\=\sum_{i=1}^{n}{(id\ast id_2\mu\ast(id\ast\varphi)^2)(i)}$$
其中 $\ast$ 表示狄利克雷卷积。

因此所求是一个积性函数的前缀和，不妨设 $f=id\ast id_2\mu\ast(id\ast\varphi)^2$ ，可以发现 $f(p^k)=p^{k-1}((2k+1)\varphi(p^{k+1})+1)(k>0)$ ，或者是 $f(p^k)=p^2f(p^{k-1})+\varphi(p^k)(2p^k-1)(k>0)$ ，可以尝试利用筛法来计算答案。

这里首先给出 $f(p^k)$ 的复杂推导，令 $g=id\ast\varphi,h=id\ast id_2\mu$ ，不难得到 $g(p^k)=(k+1)p^k-kp^{k-1}(k>0)$ 和 $h(1)=1,h(p^k)=p^k(1-p)(k>1)$ ，则有

$$f(p^k)=\sum_{p^i|p^k}{g^2(p^i)h(p^{k-i})}\\=g^2(p^k)+\sum_{i=0}^{k-1}{((i+1)p^i-ip^{i-1})^2)p^{k-i}(1-p)}\\=g^2(p^k)+p^k\sum_{i=0}^{k-1}{-(i+1)^2p^{i+1}+(i+1)(3i+1)p^i-i(3i+2)p^{i-1}+i^2p^{i-2}}\\=g^2(p^k)+p^k(-\sum_{i=0}^{k}{i^2p^i}+\sum_{i=0}^{k-1}{(i+1)(3i+1)p^i}-\sum_{i=0}^{k-2}{((i+1)(3i+5)p^i)}+\sum_{i=0}^{k-3}{(i+2)^2p^i})$$
由于 $-i^2+(i+3)(3i+1)-(i+1)(3i+5)+(i+2)^2=0$ ，故只需考虑 $i=-1,k-2,k-1,k$ 时的系数。现在考虑上和式外面的 $g^2(p^k)$ ，则有
$$f(p^k)[p^k-1]=1\\f(p^k)[p^{2k-2}]=k^2-(k-2)^2+(k-1)(3k-5)-(k-1)(3k-1)=0\\f(p^k)[p^{2k-1}]=-2k(k+1)-(k-1)^2+k(3k-2)=-2k-1\\f(p^k)[p^{2k}]=(k+1)^2-k^2=2k+1$$
故 $f(p^k)=(2k+1)(p^{2k}-p^{2k-1})+p^{k-1}=p^{k-1}((2k+1)\varphi(p^{k+1})+1)$ ，但是这样做非常麻烦，不妨回到原式来做如下推导
$$f(p^k)=\sum_{i=1}^{p^k}\sum_{j=1}^{p^k}{\text{lcm}(\gcd(i,p^k),\gcd(j,p^k))}\\=\sum_{i=0}^{k}\sum_{j=0}^{k}{p^{\max(i,j)}\varphi(p^{k-i})\varphi(p^{k-j})}\\=\sum_{i=0}^{k}{p^i\varphi(p^{k-i})(2\sum_{j=0}^{i-1}{\varphi(p^{k-j})}+\varphi(p^{k-i}))}\\=p^k(2p^k-1)+\sum_{i=0}^{k-1}{\varphi(p^k)(2p^k-(p^{k-i}+p^{k-i-1}))}\\=(2k+2)p^{2k}-2kp^{2k-1}-p^k-\sum_{i=0}^{k-1}{(p^{2k-i}-p^{2k-i-2})}\\=(2k+1)(p^{2k}-p^{2k-1})+p^{k-1}$$
相对简单许多。

由于本题模数固定，且 $N = \max(n)$ 不是很大，所以我采取了分段打表的做法。比如要计算 $f(x)(x\in[L,R])$ 的值，可以先筛出不超过 $\sqrt{R}$ 的质数，再枚举每个质数去计算对 $f(x)(x\in[L,R])$ 的贡献，这样做的复杂度是 $O(\sum_{x=L}^{R}{\Omega(x)}+\frac{\sqrt{R}}{\ln R})=O((R-L+1)\log\log R+\frac{\sqrt{R}}{\ln R})$ ，选取一个恰当的阈值 $K$ ，分别计算 $[1,K],[K+1,2K],\cdots,[\left\lfloor\frac{N-1}{K}\right\rfloor K+1,N]$ 的值，在本地打表装入代码，在处理询问时只需要处理一个长度不超过 $\frac{K}{2}$ 的区间，这样的复杂度是 $O(T(K\log\log N+\frac{\sqrt{N}}{\ln N}))$ 的，本地打表的复杂度和所有数字的质因子个数之和正相关。

显然 $K$ 越小越好，但这样会增加代码的长度，需要具体问题具体分析。对于本题，所存储的值不超过 $mod=2^{32}$ ，用十进制表示需要 10 个字节，用十六进制表示需要 8 个字节，用八十五进制表示需要 5 个字节（ASCII 可见字符集可以支持九十三进制），用 $p$ 进制表示需要 $\left\lceil\log_p{mod}\right\rceil$ 个字节。
假设可以用 $M$ 个字节存表，则 $\left\lceil\log_p{mod}\right\rceil \frac{N}{K} \leq M$ ，取 $K = \left\lfloor\frac{\left\lceil\log_p{mod}\right\rceil N}{M}\right\rfloor$ 即可。

代码：（略去打表相关的代码）

#include <cmath>
#include <stdio.h>
#include <algorithm>
const int maxn = 31623, delta = 100000, maxd = 10001, maxl = delta + 1;
int tot, pr[maxn], d[maxn], tf[maxd];
int solve(int L, int R) // [L, R]
{
    if(L > R)
        return 0;
    static int _val[maxl], _rem[maxl];
    int *val = _val - L, *rem = _rem - L, lim = (int)ceil(sqrt(R));
    for(int i = L; i <= R; ++i)
    {
        val[i] = 1;
        rem[i] = i;
    }
    for(int i = 0; i < tot && pr[i] <= lim; ++i)
        for(int j = R / pr[i] * pr[i]; j >= L; j -= pr[i])
        {
            rem[j] /= pr[i];
            int cnt = 1, pk = pr[i];
            for( ; rem[j] >= maxn && rem[j] % pr[i] == 0; rem[j] /= pr[i], ++cnt, pk *= pr[i], val[j] *= pr[i]);
            if(rem[j] < maxn)
                for( ; d[rem[j]] == pr[i]; rem[j] /= pr[i], ++cnt, pk *= pr[i], val[j] *= pr[i]);
            val[j] *= (cnt << 1 | 1) * (pr[i] - 1) * pk + 1;
        }
    int ret = 0;
    for(int i = L; i <= R; ++i)
    {
        if(rem[i] > 1)
            val[i] *= 3 * (rem[i] - 1) * rem[i] + 1;
        ret += val[i];
    }
    return ret;
}
inline int solve(int n) // [1, n]
{
    int idx = n / delta, low = idx * delta, upp = low + delta;
    return n - low <= upp - n ? tf[idx] + solve(low + 1, n) : tf[idx + 1] - solve(n + 1, upp);
}
int main()
{
    for(int i = 2; i < maxn; ++i)
    {
        if(!d[i])
            pr[tot++] = d[i] = i;
        for(int j = 0, k; (k = i * pr[j]) < maxn; ++j)
        {
            d[k] = pr[j];
            if(d[i] == pr[j])
                break;
        }
    }
    // parse table data from string to array, satisfied that tf[i] = solve(1, i * delta).
    int t;
    scanf("%d", &t);
    while(t--)
    {
        int n;
        scanf("%d", &n);
        printf("%u\n", solve(n));
    }
    return 0;
}