杜教筛

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

[sf sum_{i=1}^n(f*g)(i) = sum_{i=1}^nsum_{xy=i}f(x)g(y) = sum_{y=1}^ng(y)sum_{xyle n}f(x) = sum_{y=1}^ng(y)Sleft(leftlfloorfrac ny ight floor ight) ]

即，

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

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) << '
';
	}
	return 0;
}

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

简单的数学题

[sfsum_{i=1}^nsum_{j=1}^n ijgcd(i,j)\ sum_{i=1}^nsum_{j=1}^n ijcdot id(gcd(i,j))\ sum_{i=1}^nsum_{j=1}^n ijcdot (1*varphi)(gcd(i,j))\ sum_{i=1}^nsum_{j=1}^n ijsum_{dmid gcd(i,j)}varphi(d)\ sum_{i=1}^nsum_{j=1}^n ijsum_{dmid i,dmid j}varphi(d)\ ]

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

[sfsum_{d=1}^nvarphi(d)d^2sum_{i=1}^{lfloor n/d floor}isum_{j=1}^{lfloor n/d floor}j\ sum_{d=1}^nvarphi(d)d^2S(lfloor n/d floor)^2 ]

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

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

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

证明：

[egin{align} &((fcdot g)*(fcdot h))(n)\ &=sum_{dmid n}f(d)g(d)f(frac nd)h(frac nd)\ &= f(n)sum_{nmid d}g(d)h(frac nd)\ &= (fcdot(g*h))(n) end{align} ]

因此，对于 (varphicdot id^2)，定义 (1cdot id^2)，那么 (id^2cdot(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;
}