杜教筛

给定函数 f,g， 令 S(n) = Σ_1≤i≤n f(i)， 则有：

\[\sf \sum_{i=1}^n(f*g)(i) = \sum_{i=1}^n\sum_{xy=i}f(x)g(y) = \sum_{y=1}^ng(y)\sum_{xy\le n}f(x) = \sum_{y=1}^ng(y)S\left(\left\lfloor\frac ny\right\rfloor\right) \]

即，

\[\sf g(1)S(n) = \sum_{i=1}^n(f*g)(i) - \sum_{y=2}^ng(y)S(\lfloor n/y\rfloor) \]

int S(int n) {
	int ans = \* \sum_{i=1}^n (f*g)(i) *\
    for(int l=2, r; l<=n; l = r+1) {
      r = min(n, n / (n/l));
      ans -= Sg(l, r) * S(n / l);
    }
  return ans / g1;
}

算法应用的难点是选取 g 使得 g 的前缀和和 f*g 的前缀和都可快速计算

如果用线性筛之类的预处理处前 N^2/3 的答案， 总复杂度就是 O(N^2/3) 的了。

---

洛谷杜教筛模板

对于 φ， 众所周知有 φ * 1 = id， 显然可以快速计算。

对于 u， 众所周知有 u * 1 = ε， 显然可以快速计算。

#include<bits/stdc++.h>
typedef long long LL;
using namespace std;

const int maxn = 2147483647, maxm = 2e6;

int m, p[2000003], v[2000003], phi[2000003], mu[2000003];
LL Smu[2000003], Sphi[2000003];
void euler(int n) {
	phi[1] = mu[1] = 1;
	for(int i = 2; i <= n; ++i) {
		if(!v[i]) p[++m]=v[i]=i, mu[i]=-1, phi[i]=i-1;
		for(int j = 1; j <= m; ++j) {
			if(p[j] > v[i] || p[j] > n/i) break;
			v[i * p[j]] = p[j];
			mu[i * p[j]] = (i%p[j] ? -mu[i] : 0);
			phi[i * p[j]] = phi[i] * (i%p[j] ? p[j]-1 : p[j]);
		}
	}
	Smu[1] = mu[1], Sphi[1] = phi[1];
	for(int i = 2; i <= n; ++i)
		Smu[i] = Smu[i-1] + mu[i],
		Sphi[i] = Sphi[i-1] + phi[i];
}

map<int,LL> ansphi;
map<int,LL> ansmu;

LL SSphi(LL n) {
	if(n <= maxm) return Sphi[n];
	if(ansphi.count(n)) return ansphi[n];
	LL ans = (LL)n * (n+1) / 2;
	for(LL i = 2, j; i <= n; i = j + 1) {
		j = min(n, n / (n/i));
		ans -= (LL)(j - i + 1) * SSphi(n/i);
	}
	return ansphi[n] = ans;
}

LL SSmu(LL n) {
	if(n <= maxm) return Smu[n];
	if(ansmu.count(n)) return ansmu[n];
	LL ans = 1ll;
	for(LL i = 2, j; i <= n; i = j + 1) {
		j = min(n, n / (n/i));
		ans -= (LL)(j - i + 1) * SSmu(n/i);
	}
	return ansmu[n] = ans;
}

int main()
{
	euler(2000000);
	int T; scanf("%d",&T); while(T--) {
		int n; scanf("%d", &n); cout << SSphi(n) << ' ' << SSmu(n) << '\n';
	}
	return 0;
}

去掉 map 的 log 的方法：由于每次递归调用的总筛总是能表示成 \(\lfloor\dfrac Nx\rfloor\)， 而 \(x\ge N^{1/3}\) 的时候都会直接查预处理的表， 所以只需要记录 x， 就可以用数组记忆化了。然而实际表现并不快， 因为每次 N 要变化， 所以每次很多东西都要重新计算。

---

简单的数学题

\[\sf\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j)\\ \sum_{i=1}^n\sum_{j=1}^n ij\cdot id(gcd(i,j))\\ \sum_{i=1}^n\sum_{j=1}^n ij\cdot (1*\varphi)(gcd(i,j))\\ \sum_{i=1}^n\sum_{j=1}^n ij\sum_{d\mid gcd(i,j)}\varphi(d)\\ \sum_{i=1}^n\sum_{j=1}^n ij\sum_{d\mid i,d\mid j}\varphi(d)\\ \]

然后就可以传统艺能——交换求和顺序了

\[\sf\sum_{d=1}^n\varphi(d)d^2\sum_{i=1}^{\lfloor n/d\rfloor}i\sum_{j=1}^{\lfloor n/d\rfloor}j\\ \sum_{d=1}^n\varphi(d)d^2S(\lfloor n/d\rfloor)^2 \]

其中， \(\sf S(n) = \dfrac{n(n+1)}2\)

定义数论函数的点乘：\(\sf (f\cdot g)(n) = f(n)g(n)\)。

有引理：若 f 是完全积性函数， g,h 是数论函数， 则 \(\sf f\cdot(g*h) = (f\cdot g)*(f\cdot h)\)。

证明：

\[\begin{align} &((f\cdot g)*(f\cdot h))(n)\\ &=\sum_{d\mid n}f(d)g(d)f(\frac nd)h(\frac nd)\\ &= f(n)\sum_{n\mid d}g(d)h(\frac nd)\\ &= (f\cdot(g*h))(n) \end{align} \]

因此， 对于 \(\varphi\cdot id^2\)， 定义 \(1\cdot id^2\)， 那么 \(id^2\cdot(\varphi*1) = id^3\)， 就是 f*g 了， 可以杜教筛了。

#include <bits/stdc++.h>
typedef long long LL;
using namespace std;

const int maxm = 8003333;
LL mo;
LL ksm(LL a, LL b) {
	a %= mo;
	LL res = 1ll;
	for (; b; b >>= 1, a = a * a % mo)
		if (b & 1) res = res * a % mo;
	return res % mo;
}

int m, p[maxm], v[maxm], phi[maxm];
LL mS[maxm];
void euler(int n) {
	phi[1] = 1;
	for (int i = 2; i <= n; ++i) {
		if(!v[i]) p[++m] = v[i] = i, phi[i] = i - 1;
		for (int j = 1; j <= m; ++j) {
			if (p[j] > v[i] || p[j] > n / i) break;
			v[i * p[j]] = p[j];
			phi[i * p[j]] = phi[i] * (i % p[j] ? p[j]-1 : p[j]);
		}
	}
	for (int i = 1; i <= n; ++i) {
		mS[i] = (mS[i-1] + (LL)phi[i] * i % mo * i % mo) % mo;	
	}
}


LL n, inv2, inv4, inv6;
map<int,LL> ansS;

LL Sum(LL x) {
	x %= mo;
	return x * (x+1) % mo * inv2 % mo;
}
LL S2(LL x) {
	x %= mo;
	return x * (x+1) % mo * (2*x+1) % mo * inv6 % mo;
}

LL S3(LL x) {
	x %= mo;
	return x * x % mo * (x+1) % mo * (x+1) % mo * inv4 % mo;
}
LL S(LL n) {
	if(n <= 8000000) return mS[n];
	if(ansS.count(n)) return ansS[n];
	LL ans = S3(n);
	for (LL i=2, j; i <= n; i = j + 1) {
		j = min(n, n / (n / i));
		ans -= (S2(j) - S2(i-1) + mo) % mo * S(n / i) % mo;
	}
	ans = (ans%mo + mo) % mo;
	return ansS[n] = ans;
}

int main ()
{ 
	cin >> mo >> n;
	euler(8000000);
	inv2 = ksm(2, mo - 2), inv4 = ksm(4, mo - 2), inv6 = ksm(6, mo - 2);
	LL ans = 0ll;
	for (LL i = 1, j; i <= n; i = j + 1) {
		j = min (n, n / (n / i));
		LL t = Sum (n / i);
		ans = (ans + (S(j) - S(i-1)) * t % mo * t % mo) % mo;
	}
	ans = (ans % mo + mo) % mo;
	cout << ans;
	return 0;
}