BZOJ 4176: Lucas的数论 [杜教筛]

4176: Lucas的数论

题意：求\(\sum_{i=1}^n \sum_{j=1}^n \sigma_0(ij)\) \(n \le 10^9\)

---

代入\(\sigma_0(nm)=\sum_{i\mid n}\sum_{j\mid m}[(i,j)=1]\)

反演得到

\[\sum_{d=1}^n \mu(d) (g(\frac{n}{d}))^2 \\
g(n) = \sum_{i=1}^n \sigma_0(i)
\]

杜教筛\(\mu \ \sigma_0\)的前缀和

当然和前面的题一样，\(\sigma_0\)也可以用预处理+分块

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

把上题TLE的杜教筛代码抄上就行了

#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], mu[N], lp[N], si[N];

void sieve(int n) {

	mu[1]=1; si[1]=1;

	for(int i=2; i<=n; i++) {

		if(!notp[i]) p[++p[0]] = i, mu[i] = -1, si[i] = lp[i] = 2;

		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) {

				mu[t] = 0;

				lp[t] = lp[i] + 1;

				si[t] = si[i] / lp[i] * lp[t];

				break;

			}

			mu[t] = -mu[i];

			lp[t] = 2;

			si[t] = si[i] * 2;

		}

		mu[i] += mu[i-1];

		si[i] += si[i-1];

	}

}

namespace ha {

	const int p = 1001001;

	struct ha {

		struct meow{int ne, val, r;} e[10000];

		int cnt, h[p];

		ha() {cnt=0; memset(h, 0, sizeof(h));}

		inline void insert(int x, int val) {

			int u = x%p;

			for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return;

			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;

		}

	} hs, hu;

} using ha::hs; using ha::hu;

int dj_u(int n) {

	if(n <= U) return mu[n];

	if(hu.quer(n) != -1) return hu.quer(n);

	int ans = 1, r;

	for(int i=2; i<=n; i=r+1) {

		r = n/(n/i);

		mod(ans -= (r-i+1) * dj_u(n/i) %mo);

	}

	hu.insert(n, ans);

	return ans;

}

int dj_s(int n) {

	if(n <= U) return si[n];

	if(hs.quer(n) != -1) return hs.quer(n);

	int ans = n, r, now, last = dj_u(1);

	for(int i=2; i<=n; i=r+1, last=now) {

		r = n/(n/i); now = dj_u(r);

		mod(ans -= (ll) (now - last) * dj_s(n/i) %mo);

	}

	hs.insert(n, ans);

	return ans;

}

int solve(int n) {

	dj_u(n); dj_s(n);

	int ans=0, r, last=0, now;

	for(int i=1; i<=n; i=r+1, last=now) {

		r = n/(n/i); now = dj_u(r); ll t = dj_s(n/i); //printf("hi [%d, %d]  %d %d  %lld\n", i, r, now, last, t);

		mod(ans += (ll) (now - last) * t %mo * t %mo);

	}

	return ans;

}

int main() {

	freopen("in", "r", stdin);

	int n=read(); //U = pow(n, 2.0/3);

	sieve(U);

	printf("%d", solve(n));

}