【LG3768】简单的数学题

【LG3768】简单的数学题

题面

求
\[ (\sum_{i=1}^n\sum_{j=1}^nij\text{gcd}(i,j))\text{mod}p \]

其中\(n\leq 10^{10},5\times 10^8\leq p \leq 1.1*10^9\)。

题解

推柿子：

\[ \sum_{i=1}^n\sum_{j=1}^nij\text{gcd}(i,j)\\ =\sum_{d=1}d\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ijd^2[\gcd(i,j)==1]\\ =\sum_{d=1}d^3\sum_{i=1}^{\lfloor\frac nd\rfloor}\sum_{j=1}^{\lfloor\frac nd\rfloor}ij[\gcd(i,j)==1]\\ \]

令

\[ f(d')=\sum_{i=1}^{n'}\sum_{j=1}^{n'}ij[\gcd(i,j)==d']\\ g(d')=\sum_{d'\mid x}f(x) \]

则

\[ g(d')=\sum_{d'\mid x}\sum_{i=1}^{n'}\sum_{j=1}^{n'}ij[\gcd(i,j)==x]\\ =\sum_{i=1}^{n'}\sum_{j=1}^{n'}ij[d'|\gcd(i,j)]\\ =\sum_{i=1}^{\lfloor\frac {n'}{d'}\rfloor}\sum_{j=1}^{\lfloor\frac {n'}{d'}\rfloor}ijd'^2[1|\gcd(i,j)]\\ =S(\lfloor\frac {n'}{d'}\rfloor)^2d'^2 \]

其中\(S(x)=\sum_{i=1}^x i\)

那么

\[ f(d')=\sum_{d'\mid x}S(\lfloor\frac {n'}x\rfloor)^2x^2\mu (x)\\ f(1)=\sum_{x=1}^{n'}S(\lfloor\frac {n'}x\rfloor)^2x^2\mu (x) \]

代回去

\[ \sum_{d=1}^nd^3\sum_{x=1}^{\lfloor\frac nd\rfloor}S(\lfloor\frac {n}{xd}\rfloor)^2x^2\mu (x) \]

令\(Q=xd\)，

\[ \sum_{d=1}^nd^3\sum_{x=1}^{\lfloor\frac nd\rfloor}S(\lfloor\frac {n}{Q}\rfloor)^2x^2\mu (x)\\ =\sum_{Q=1}^nS(\lfloor\frac {n}{Q}\rfloor)\sum_{d\mid Q}d^3(\frac Qd)^2\mu (\frac Qd)\\ =\sum_{Q=1}^nS(\lfloor\frac {n}{Q}\rfloor)^2\sum_{d\mid Q}dQ^2\mu (\frac Qd)\\ =\sum_{Q=1}^nS(\lfloor\frac {n}{Q}\rfloor)^2Q^2\varphi (Q) \]

现在我们如何求\(Q^2\varphi (Q)\)的前缀和呢？

令\(f=\text{id}^{2}\cdot \varphi\)，令\(g=\text{id}^2\cdot 1\)，

那么\(f*g=\text{id}^3,g=\text{id}^2\)。

就做完了。

#include <iostream> 
#include <cstdio> 
#include <cstdlib> 
#include <cstring> 
#include <cmath> 
#include <algorithm>
#include <map> 
using namespace std; 
const int MAX = 1e7, MAX_N = 1e7 + 5; 
bool nprime[MAX_N]; 
int prime[MAX_N], cnt, phi[MAX_N], f[MAX_N];
long long N; 
int Mod; 
int fpow(int x, int y) { 
    int res = 1;
    while (y) {
        if (y & 1) res = 1ll * res * x % Mod; 
        x = 1ll * x * x % Mod;
        y >>= 1; 
    } 
    return res; 
} 
void sieve() { 
    phi[1] = 1; 
    for (int i = 2; i <= MAX; i++) { 
        if (!nprime[i]) prime[++cnt] = i, phi[i] = i - 1; 
        for (int j = 1; j <= cnt && prime[j] * i <= MAX; j++) { 
            nprime[prime[j] * i] = 1; 
            if (i % prime[j] == 0) { phi[i * prime[j]] = 1ll * phi[i] * prime[j] % Mod; break; } 
            phi[i * prime[j]] = 1ll * phi[i] * phi[prime[j]] % Mod; 
        } 
    } 
    for (int i = 1; i <= MAX; i++) f[i] = (f[i - 1] + 1ll * phi[i] * i % Mod * i % Mod) % Mod; 
}
int inv2, inv6; 
int sqr(int n) { return 1ll * n * n % Mod; } 
int S(long long n) { return n %= Mod, n % Mod * (n + 1) % Mod * inv2 % Mod; } 
int g(long long n) { return n %= Mod, n % Mod * (n + 1) % Mod * (2 * n + 1) % Mod * inv6 % Mod; } 
map<long long, int> mp; 
int get_f(long long n) { 
    if (n <= MAX) return f[n]; 
    if (mp.find(n) != mp.end()) return mp[n]; 
    int res = sqr(S(n)); 
    for (long long l = 2, r; l <= n; l = r + 1) { 
        r = n / (n / l);
        res = (res - 1ll * (g(r) - g(l - 1) + Mod) % Mod * get_f(n / l) % Mod + Mod) % Mod; 
    } 
    return mp[n] = res; 
} 
int main () { 
#ifndef ONLINE_JUDGE 
    freopen("cpp.in", "r", stdin); 
#endif 
    cin >> Mod >> N; 
    sieve(); 
    inv2 = fpow(2, Mod - 2), inv6 = fpow(6, Mod - 2); 
    int ans = 0; 
    for (long long l = 1, r; l <= N; l = r + 1) { 
        r = N / (N / l); 
        ans = (ans + 1ll * sqr(S(N / l)) * ((get_f(r) - get_f(l - 1) + Mod) % Mod) % Mod) % Mod; 
    }
    printf("%d\n", ans); 
    return 0; 
} 

转载于:https://www.cnblogs.com/heyujun/p/10623894.html