PE 439 Sum of sum of divisors | 51nod 1220 约数之和

题目：https://projecteuler.net/problem=439以及http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1220

题意：
令$d(x)$表示$x$的所有因子之和，求$\sum_{i=1}^{N}{\sum_{j=1}^{N}{d(i*j)}}\mod10^9+7$。
$N\le10^9$。

题解：
令$[x]$为布尔表达式$x$的值，即若$x$为真则$[x]=1$，否则$[x]=0$。
那么有$d(x)=\sum_{i=1}^{x}{[i|x]*i}$
于是问题可以写成

$$\sum_{i=1}^{N}{\sum_{j=1}^{N}{d(i*j)}}=\sum_{i=1}^{N}{\sum_{j=1}^{N}{\sum_{d=1}^{i*j}{[d|i*j]*d}}}$$
考虑如何简化掉 $j$，我们可以计算每个$d$对应多少个 $j$。因为$d|i*j$等价于 $\frac{d}{gcd(i,d)}|j$，所以上式可以化为
$$=\sum_{i=1}^{N}{\sum_{d=1}^{i*N}{\lfloor\frac{N*gcd(i,d)}{d}\rfloor}*d}$$
再考虑如何将 $gcd(i,d)$化为互质的情况，令 $d'=gcd(i,d)$，则有
$$=\sum_{d'=1}^{N}{\sum_{i=1}^{N}{\sum_{d=1}^{i*N}{[gcd(i,d)=d']*\lfloor\frac{N*gcd(i,d)}{d}\rfloor*d}}}$$
由于 $\frac{i}{d'}$和 $\frac{d}{d'}$是整数，我们将上式的 $i$和$d$用 $\frac{i}{d'}*d'$和 $\frac{i}{d'}*d'$替换
$$=\sum_{d'=1}^{N}{\sum_{\frac{i}{d'}=1}^{\lfloor\frac{N}{d'}\rfloor}{\sum_{\frac{d}{d'}=1}^{\frac{i}{d'}*N}{[gcd(\frac{i}{d'},\frac{d}{d'})=1]*\lfloor\frac{N*d'}{d}\rfloor*\frac{d}{d'}*d'}}}$$
再进行一个下标变换，将 $\frac{i}{d'}$和 $\frac{d}{d'}$用 $i$和$d$替换
$$=\sum_{d'=1}^{N}{\sum_{i=1}^{\lfloor\frac{N}{d'}\rfloor}{\sum_{d=1}^{i*N}{[gcd(i,d)=1]*\lfloor\frac{N}{d}\rfloor*d*d'}}}\\=\sum_{d'=1}^{N}{d'\sum_{i=1}^{\lfloor\frac{N}{d'}\rfloor}{\sum_{d=1}^{N}{[gcd(i,d)=1]*\lfloor\frac{N}{d}\rfloor*d}}}$$
可以看到 $d$和$d'$之间没有约束，而且 $i$和$d'$的关系比 $d$的关系简单，不妨将$d'$和 $i$交换求和顺序，考虑每个$i$对应哪些 $d'$，尝试将 $d'$约去
$$=\sum_{i=1}^{N}{\sum_{d'=1}^{\lfloor\frac{N}{i}\rfloor}{d'*\sum_{d=1}^{N}{[gcd(i,d)=1]*\lfloor\frac{N}{d}\rfloor*d}}}\\=\sum_{i=1}^{N}{\frac{\lfloor\frac{N}{i}\rfloor*(\lfloor\frac{N}{i}\rfloor+1)}{2}\sum_{d=1}^{N}{[gcd(i,d)=1]*\lfloor\frac{N}{d}\rfloor*d}}$$
而对于处理 $[gcd(i,d)=1]$，想必读到这里的你已经想出了处理它的方法，利用莫比乌斯函数的经典定理 $\sum_{i=1}^{N}{[i|N]*\mu(i)}=[N=1]$可以得到
$$=\sum_{i=1}^{N}{\frac{\lfloor\frac{N}{i}\rfloor*(\lfloor\frac{N}{i}\rfloor+1)}{2}\sum_{d=1}^{N}{\lfloor\frac{N}{d}\rfloor*d\sum_{d'=1}^{N}{[d'|gcd(i,d)]*\mu(d')}}}\\=\sum_{i=1}^{N}{\frac{\lfloor\frac{N}{i}\rfloor*(\lfloor\frac{N}{i}\rfloor+1)}{2}\sum_{d=1}^{N}{\lfloor\frac{N}{d}\rfloor*d\sum_{d'=1}^{N}{[d'|i]*[d'|d]*\mu(d')}}}$$
细心的你一定发现了 $i$和$d$互相无关，它们分别只和 $d'$有关，于是我们交换求和顺序，算每个 $d'$对应的 $i$和$d$，则有
$$=\sum_{i=1}^{N}{\frac{\lfloor\frac{N}{i}\rfloor*(\lfloor\frac{N}{i}\rfloor+1)}{2}\sum_{d=1}^{N}{\lfloor\frac{N}{d}\rfloor*d\sum_{d'=1}^{N}{[d'|i]*[d'|d]*\mu(d')}}}\\=\sum_{d'=1}^{N}{\mu(d')\sum_{\frac{i}{d'}=1}^{\lfloor\frac{N}{d'}\rfloor}{\frac{\lfloor\frac{N}{\frac{i}{d'}*d'}\rfloor*(\lfloor\frac{N}{\frac{i}{d'}*d'}\rfloor+1)}{2}}\sum_{\frac{d}{d'}=1}^{\lfloor\frac{N}{d'}\rfloor}{\lfloor\frac{N}{\frac{d}{d'}*d'}\rfloor*\frac{d}{d'}*d'}}\\=\sum_{d'=1}^{N}{d'*\mu(d')\sum_{i=1}^{\lfloor\frac{N}{d'}\rfloor}{\frac{\lfloor\frac{N}{i*d'}\rfloor*(\lfloor\frac{N}{i*d'}\rfloor+1)}{2}}\sum_{d=1}^{\lfloor\frac{N}{d'}\rfloor}{\lfloor\frac{N}{d*d'}\rfloor*d}}$$
到这里，我们已经得到了一个不错的式子，可以直接进行计算了，无需继续化简。（再化简下去 $\LaTeX$都要炸啦~）
令 $f(x)=x*\mu(x)$， $g(x)=\sum_{i=1}^{x}{\frac{\lfloor\frac{x}{i}\rfloor*(\lfloor\frac{x}{i}\rfloor+1)}{2}}$， $h(x)=\sum_{i=1}^{x}{\lfloor\frac{x}{i}\rfloor*i}$，则所求为 $\sum_{i=1}^{N}{f(i)*g(\lfloor\frac{N}{i}\rfloor)*h(\lfloor\frac{N}{i}\rfloor)}$。
这个式子里用到的 $g(x)$和 $h(x)$共 $O(\sqrt{N})$个，每个的计算复杂度为 $O(\sqrt{x})$，所以把所有所需的 $g$和$h$算出来的复杂度为 $O(\sum_{i=1}^{\sqrt{N}}{\sqrt{\frac{N}{i}}})=O(N^\frac{3}{4})$，但是达不到这个上界。（实际上很容易改成 $O(N^\frac{2}{3})$的，并且 $g(x)=h(x)$）
如果可以快速计算出一段 $f(x)$的值，就可以结合上面的做法解决这个问题，注意到 $f(x)$是两个积性函数的乘积，所以 $f(x)$也是积性函数，我们可以得到类似梅滕斯函数（莫比乌斯函数的前缀和函数）的计算方法，令 $S(x)=\sum_{i=1}^{x}{f(i)}$，则有
$$S(x)=1-\sum_{i=2}^{x}{\sum_{d=1}^{i-1}{[d|i]*i*\mu(d)}}\\=1-\sum_{i=2}^{x}{\sum_{d=1}^{\lfloor\frac{x}{i}\rfloor}{i*d*\mu(d)}}\\=1-\sum_{i=2}^{x}{i*S(\lfloor\frac{x}{i}\rfloor)}$$
于是我们可以利用线性筛预处理出前 $T$个值，那么计算复杂度变为$O(\sum_{i=1}^{\frac{N}{T}}{\sqrt{\frac{N}{i}}})=O(\frac{N}{\sqrt{T}})$，当 $T$取$N^\frac{2}{3}$的时候可以做到 $O(N^\frac{2}{3})$。
于是总体复杂度为 $O(N^\frac{3}{4})$。

代码：

#include <map>
#include <cmath>
#include <cstdio>
using namespace std;
typedef long long LL;
const int maxn = 1000010, mod = 1000000007, inv2 = 500000004;
map<LL, int> Hash;
int n, sqn, tot, prime[maxn], f[maxn], ans;
bool vis[maxn];
inline void inc(int &x, int y)
{
    x += y;
    if(x >= mod)
        x -= mod;
}
inline void dec(int &x, int y)
{
    x -= y;
    if(x < 0)
        x += mod;
}
inline int num1(int x)
{
    return x * (x + 1LL) % mod * inv2 % mod;
}
inline int num1(int L, int R)
{
    int ret = num1(R);
    dec(ret, num1(L - 1));
    return ret;
}
int calc_imu(int x)
{
    if(x <= sqn)
        return f[x];
    if(Hash.count(x))
        return Hash[x];
    int ret = 1;
    for(int i = 2, j; i <= x; i = j + 1)
    {
        j = x / (x / i);
        dec(ret, (LL)num1(i, j) * calc_imu(x / i) % mod);
    }
    return Hash[x] = ret;
}
inline int calc_imu(int L, int R)
{
    int ret = calc_imu(R);
    dec(ret, calc_imu(L - 1));
    return ret;
}
int calc_g(int n)
{
    int ret = 0;
    for(int i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        inc(ret, (j - i + 1LL) * num1(n / i) % mod);
    }
    return ret;
}
int calc_h(int n)
{
    int ret = 0;
    for(int i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        inc(ret, (LL)(n / i) * num1(i, j) % mod);
    }
    return ret;
}
int main()
{
    scanf("%d", &n);
    sqn = (int)ceil(pow(n, 2.0 / 3));
    f[1] = 1;
    for(int i = 2; i <= sqn; ++i)
    {
        if(!vis[i])
        {
            prime[tot++] = i;
            dec(f[i], 1);
        }
        for(int j = 0, k = sqn / i, o; j < tot && prime[j] <= k; ++j)
        {
            vis[o = i * prime[j]] = 1;
            if(i % prime[j] == 0)
            {
                f[o] = 0;
                break;
            }
            else
                dec(f[o], f[i]);
        }
        f[i] = (f[i - 1] + (LL)i * f[i]) % mod;
    }
    for(int i = 1, j; i <= n; i = j + 1)
    {
        j = n / (n / i);
        inc(ans, (LL)calc_imu(i, j) * calc_g(n / i) % mod * calc_g(n / i) % mod);
    }
    printf("%d\n", ans);
    return 0;
}