51NOD 1220 约数之和 [杜教筛]

1220 约数之和

题意：求\(\sum_{i=1}^n \sum_{j=1}^n \sigma_1(ij)​\)

---

\[\sigma_0(ij) = \sum_{x\mid i}\sum_{y\mid j}[(x,y)=1]\\ \sigma_1(ij) = \sum_{x\mid i}\sum_{y\mid j}x\cdot\frac{j}{y}[(x,y)=1] \\ \]

怎么证明第二个式子？

\[\sigma_1(n) = \prod_i(1 + p_i + p_i^2 + ... + p_i^{a_i}) \]

和第一个式子类似，每个质因子是独立的，x和y枚举了i和j的所有质因子组合，gcd为1保证了对于每个质因子\(p_i\)，x和y有一个指数为0.这样实际上不重不漏的枚举了一个质因子在ij中的所有指数。

反演后得到

\[\sum_{d=1}^n \mu(d) \cdot d \cdot g^2(\frac{n}{d}) \\ g(n) = \sum_{i=1}^n \sigma_1(i)=\sum_{i=1}^n i \lfloor \frac{n}{i} \rfloor \]

杜教筛\(\mu \cdot id\)，卷上\(id\)即可；

\(g\)预处理+分块。

总复杂度\(O(n^{\frac{2}{3}})\)

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
#include <ctime>
using namespace std;
typedef long long ll;
const int N=1664512, mo=1e9+7;
int U=1664510;
inline ll read(){
	char c=getchar(); ll x=0,f=1;
	while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
	while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
	return x*f;
}

inline void mod(int &x) {if(x>=mo) x-=mo; else if(x<0) x+=mo;}
inline void mod(ll &x) {if(x>=mo) x-=mo; else if(x<0) x+=mo;}
bool notp[N]; int p[N/10], lp[N], si[N], mu[N];
void sieve(int n) {
	si[1] = 1; mu[1] = 1;
	for(int i=2; i<=n; i++) {
		if(!notp[i]) p[++p[0]] = i, lp[i] = i, si[i] = 1+i, mu[i] = -1;
		for(int j=1; j <= p[0] && i*p[j] <= n; j++) {
			int t = i*p[j]; notp[t] = 1; 
			if(i % p[j] == 0) {
				lp[t] = lp[i] * p[j];
				if(lp[t] == t) mod(si[t] = si[i] + lp[t]);
				else si[t] = (ll) si[t / lp[t]] * si[lp[t]] %mo;
				mu[t] = 0;
				break;
			}
			lp[t] = p[j];
			si[t] = (ll) si[i] * (1 + p[j]) %mo;
			mu[t] = -mu[i];
		}
	}
	for(int i=1; i<=n; i++) mod(si[i] += si[i-1]), mod(mu[i] = mu[i-1] + mu[i] * i);
}

namespace ha {
	const int p = 1001001;
	struct meow{int ne, val, r;} e[3000];
	int cnt, h[p];
	inline void insert(int x, int val) {
		int u = x%p;
		e[++cnt] = (meow){h[u], val, x}; h[u] = cnt;
	}
	inline int quer(int x) {
		int u = x%p;
		for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return e[i].val;
		return -1;
	}
} using ha::insert; using ha::quer;

inline ll sum(ll n) {return n * (n + 1) / 2 %mo;}
int dj_u(int n) {
	if(n <= U) return mu[n];
	if(quer(n) != -1) return quer(n);
	int ans=1, r;
	for(int i=2; i<=n; i=r+1) {
		r = n/(n/i);
		mod(ans -= (ll) (sum(r) - sum(i-1)) * dj_u(n/i) %mo);
	}
	insert(n, ans);
	return ans;
}

int g(int n) {
	if(n <= U) return si[n];
	int ans=0, r;
	for(int i=1; i<=n; i=r+1) {
		r = n/(n/i);
		mod(ans += (ll) (sum(r) - sum(i-1)) * (n/i) %mo);
	}
	return ans;
}

int solve(int n) {
	int ans=0, r, now, last=0;
	for(int i=1; i<=n; i=r+1, last=now) {
		r = n/(n/i); now = dj_u(r); ll t = g(n/i); //printf("solve [%d, %d] %d  %d   %d\n", i, r, now, last, g(n/i));
		mod(ans += (ll) (now - last) * t %mo * t %mo);
	}
	return ans;
}

int main() {
	freopen("in", "r", stdin);
	sieve(U);
	int n=read();
	printf("%d", solve(n));
}