P3768 的小推导

看了主页的莫反笔记，做到 P3768，看到上杜教筛那里不大理解，所以再来推一次吧。

---

求 \(\sum \limits _{i=1}^n \sum \limits _{j=1}^n ij \gcd(i,j)\)，\(n \leq 10^9\)。

先推式子

\[\newcommand{\frc}[2]{ \lfloor \frac{#1}{#2} \rfloor } \newcommand{\sm}[3]{ \sum \limits _{#1 = #2} ^ {#3} } \begin{aligned} & = \sm{d}{1}{n} d^3 \sum _{i=1}^{\frc{n}{d}} \sum _{j=1}^{\frc{n}{d}} ij [\gcd(i,j)=1] \\ & = \sm{d}{1}{n} d^3 \sum _{i=1}^{\frc{n}{d}} \sum _{j=1}^{\frc{n}{d}} ij \sum \limits _{p|i,p|j} \mu(p) \\ & = \sm{d}{1}{n} d^3 \sum _{p=1}^{\frc{n}{d}} \mu(p) p^2 \sum _{i=1}^{\frc{n}{dp}} i \sum _{j=1}^{\frc{n}{dp}} j \\ & = \sm{d}{1}{n} d^3 \sum_{p=1}^{\frc{n}{d}} \mu(p) p^2 s^2(\frc{n}{dp}) \ \ \ \ where \ s(x) = \frc{x(x+1)}{2} \\ & = \sm{T}{1}{n} s^2(\frc{n}{T}) \sum _{d|T} \mu(\frac{T}{d}) d^3 \left ( \frac{T}{d} \right ) ^2 \\ & = \sm{T}{1}{n} s^2(\frc{n}{T}) T^2 \sum _{d|T} \mu(d) \times \frac{T}{d} \ \ \gets which \ equals \ \varphi(T) \\ & = \sm{T}{1}{n} s^2 (\frc{n}{T}) T^2 \varphi(T) \\ & = \sm{T}{1}{n} s^2(\frc{n}{T}) f(T) \end{aligned} \]

前面那个东西可以整除分块，所以要给后面的东西卷点东西做杜教筛。

我们知道 \((\varphi \times I)(x) =x\)，所以给 \(T^2 \varphi(T)\) 卷个 \(x^2\) 即可化简。

\[\begin{aligned} & \sum _{d|n} d^2 \varphi(d) \times \left (\frac{n}{d} \right ) ^ 2 \\ = & n^2 \sum_{d|n} \varphi(d) \\ = & n^3 \end{aligned} \]

杜教筛：

\[\begin{aligned} g(1)S_f(n) = & S_h(n) - \sm{d}{2}{n} g(d) \times S_f(\frc{n}{d}) \\ 1 \cdot S_f(n) = & s^2(i) - \sm{d}{2}{n} d^2 \times S_f(\frc{n}{d}) \ \ \\ \end{aligned} \]

结束。

#include <cstdio>
#include <map>
#define LL long long
// #define int long long
using namespace std;
const int M = 1000005, m = 1000000;
int mod, inv6, inv2; map<LL, int> t; LL n;
int qpow(int a, int b){
	long long ans = 1ll;
	for(; b; b >>= 1) {if(b & 1) ans = 1ll * ans * a % mod; a = 1ll * a * a % mod;}
	return ans;
}
int inv(int k) {return qpow(k, mod - 2);}
int addd(int a, int b) {a += b; return a > mod ? a-mod : a;}
int mnss(int a, int b) {a -= b; return a < 0 ? a+mod : a;}
void add(int &x, int y) {x += y; if(x > mod) x -= mod;}
void mns(int &a, int b) {a -= b; if(a < 0) a += mod;}
int p[M], tt, phi[M], sf[M]; bool f[M];
void pre(int n = m){
	f[1] = 1; sf[1] = phi[1] = 1;
	for(int i = 2; i <= n; i++){
		if(!f[i]) p[++tt] = i, phi[i] = i-1;
		for(int j = 1; j <= tt; j++){
			int k = p[j] * i; if(k > n) break;
			if(i % p[j] == 0) {f[k] = 1; phi[k] = phi[i] * p[j]; break;}
			f[k] = 1; phi[k] = phi[i] * (p[j]-1);
		}
        sf[i] = addd(sf[i-1], 1ll * i * i % mod * phi[i] % mod) % mod;
	} 
}
int s(LL n) {n %= mod; return 1ll * n * (n+1) % mod * inv2 % mod;}
int ss(LL n) {n %= mod; return 1ll * n * (n+1) % mod * (2*n+1) % mod * inv6 % mod;}
int Sf(LL n){
    if(n <= m) return sf[n];
    if(t[n]) return t[n];
    int ans = 1ll * s(n) * s(n) % mod;
    for(LL l = 2, r; l <= n; l = r+1){
        r = n/(n/l); mns(ans, 1ll * mnss(ss(r), ss(l-1)) * Sf(n/l) % mod);
    }
    return t[n] = ans;
}
int solve(LL n){
    int ans = 0;
    for(LL l = 1, r; l <= n; l = r+1) {
        r = n/(n/l);
        add(ans, 1ll * s(n/l) * s(n/l) % mod * mnss(Sf(r), Sf(l-1)) % mod);
    }
    return ans;
}
signed main(){
    // freopen("temp.out", "w", stdout);
    scanf("%d %lld", &mod, &n);
    inv6 = inv(6); pre(); inv2 = inv(2);
    printf("%d\n", solve(n));
}