[51Nod 1584] 加权约数和

Posted 2021-02-09 fly-in-milkyway

Description

在整理以前的试题时，他发现了这样一道题目：“求 (sumsigma(i))，其中 (1≤i≤N)，(σ(i)) 表示 (i) 的约数之和。”

现在他长大了，题目也变难了，所以麻烦你来帮他解决一道数论题吧。

他需要你求如下表达式的值：

[ sum_{i=1}^Nsum_{j=1}^Nmax(i,j)?σ(i?j) ]
其中 (max(i,j)) 表示 (i) 和 (j) 里的最大值，(sigma(i?j)) 表示 (i?j) 的约数之和。

例如当 (N=2) 的时候，由 (σ(1)=1,σ(2)=1+2=3,σ(4)=1+2+4=7) 可知，答案应为 (1?σ(1?1)+2?σ(1?2)+2?σ(2?1)+2?σ(2?2)=27)。

他发现答案有点大，所以你只需要告诉他答案模 (1000000007) 的值即可。

Input

每个测试点含有多组测试数据。

第一行是一个正整数 (T)，表示接下来有 (T) 组测试数据。

接下来的 (T) 行，每组测试数据占一行。

每行有一个正整数 (N)，含义如描述所示。

Output

共有 (T) 行。对于每组测试数据，输出一行信息"Case #x: y"。

其中 (x) 表示对应的是第几组测试数据，(y) 表示相应的答案模 (1000000007) 的值。

Sample Input

5
1
2
3
4
5

Sample Output

Case #1: 1
Case #2: 27
Case #3: 162
Case #4: 686
Case #5: 1741

HINT

(1≤T≤50000,1≤N≤1000000)

Solution

〖一〗

[ sigma(icdot j)=sumlimits_{pmid i}sumlimits_{qmid j}[left(p,q)=1 ight]dfrac{pcdot j}{q} ]

证明：原式 (=sumlimits_{pmid i}sumlimits_{qmid j}left[left(p,dfrac{j}{q} ight)=1 ight]pcdot q)，设 (i=sum p_i^{a_i},j=sum p_i^{b_i},p=sum p_i^{c_i},q=sum p_i^{d_i})。

• 若 (0<c_ile a_i)，则 (d_i=b_i)，此时可以表示出 (p_i^{(b_i+1) ightarrow(b_i+a_i)})；
• 若 (c_i=0)，则 (0le d_ile b_i)，此时可以表示出 (p_i^{0 ightarrow b_i})。

综上，该式一定可以表示出 (p_i^{0 ightarrow(b_i+a_i)})。

〖二〗

[ sum_{i=1}^nsum_{j=1}^nmax(i,j)?σ(i?j)=2sum_{i=1}^{n}sum_{j=1}^{i}icdotsigma(icdot j)-sum_{i=1}^nicdotsigma(i^2) ]

[ egin{eqnarray} sum_{i=1}^nicdotsigma(i^2)&=&sum_{i=1}^{n}isum_{pmid i}sum_{qmid i}[(p,q)=1]frac{pcdot i}{q}&=&sum_{d=1}^{n}mu(d)sum_{i=1}^{n}isum_{pmid i}sum_{qmid i}[dmid(p,q)]frac{pcdot i}{q}&=&sum_{d=1}^{n}mu(d)sum_{i=1}^{n}isum_{dmid pmid i}sum_{dmid qmid i}frac{pcdot i}{q}&=&sum_{d=1}^{n}mu(d)cdot d^2sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}isum_{pmid i}psum_{qmid i}frac{i}{q}&=&sum_{d=1}^{n}mu(d)cdot d^2sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}ileft(sum_{pmid i}p ight)^2 end{eqnarray} ]

[ egin{eqnarray} sum_{i=1}^{n}sum_{j=1}^{i}icdot sigma(icdot j)&=&sum_{i=1}^{n}sum_{j=1}^{i}isum_{pmid i}sum_{qmid j}[(p,q)=1]frac{pcdot j}{q}&=&sum_{d=1}^{n}mu(d)sum_{i=1}^{n}sum_{j=1}^{i}isum_{pmid i}sum_{qmid j}[dmid(p,q)]frac{pcdot j}{q}&=&sum_{d=1}^{n}mu(d)cdot d^2sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}isum_{pmid i}psum_{j=1}^{i}sum_{qmid j}frac{j}{q}&=&sum_{d=1}^{n}mu(d)cdot d^2sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}isum_{pmid i}psum_{j=1}^{i}sum_{qmid j}q end{eqnarray} ]

而这样只能做到 (O(sqrt n)) 询问，于是进一步化简：

[ egin{eqnarray} sum_{d=1}^{n}mu(d)cdot d^2sum_{i=1}^{leftlfloorfrac{n}{d} ight floor}f(i)&=&sum_{i=1}^nsum_{dmid i}mu(d)cdot d^2cdot fleft(frac{i}{d} ight) end{eqnarray} ]

第一个式子中的 (d) 是 (icdot d) 的约数，第二个式子中的 (i) 表示 (icdot d)，(d) 还是约数，因此两个式子相等。

然后就可以 (O(nln n)) 预处理，(O(1)) 查询了。

Code

#include <cstdio>

const int N = 1000005, mod = 1000000007;
int mu[N], np[N], p[N], tot, a[N], b[N], f[N], g[N], n = 1000000;

int read() {
    int x = 0; char c = getchar();
    while (c < '0' || c > '9') c = getchar();
    while (c >= '0' && c <= '9') x = (x << 3) + (x << 1) + (c ^ 48), c = getchar();
    return x;
}
void euler() {
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if (!np[i]) p[++tot] = i, mu[i] = -1;
        for (int j = 1; j <= tot && i * p[j] <= n; ++j) {
            np[i * p[j]] = 1;
            if (i % p[j] == 0) { mu[i * p[j]] = 0; break; }
            mu[i * p[j]] = -mu[i];
        }
    }
}
int main() {
    euler();
    for (int i = 1; i <= n; ++i)
        for (int j = i; j <= n; j += i)
            if ((a[j] += i) >= mod) a[j] -= mod;
    for (int i = 1; i <= n; ++i) if ((b[i] = b[i - 1] + a[i]) >= mod) b[i] -= mod;
    for (int i = 1; i <= n; ++i)
        for (int j = i; j <= n; j += i) {
            f[j] = (f[j] + 1LL * mu[i] * i * j % mod * a[j / i] % mod * b[j / i]) % mod;
            g[j] = (g[j] + 1LL * mu[i] * i * j % mod * a[j / i] % mod * a[j / i]) % mod;
        }
    for (int i = 2; i <= n; ++i) {
        if ((f[i] += f[i - 1]) >= mod) f[i] -= mod;
        if (f[i] < 0) f[i] += mod; //mu可能是负数, 需要判断正负
        if ((g[i] += g[i - 1]) >= mod) g[i] -= mod;
        if (g[i] < 0) g[i] += mod;
    }
    int T = read();
    for (int i = 1; i <= T; ++i) n = read(), printf("Case #%d: %lld
", i, (2LL * f[n] + mod - g[n]) % mod);
    return 0;
}