杜教筛

前置知识：莫比乌斯反演

杜教筛

杜教筛可以在非线性时间内求积性函数前缀和。

P4213 【模板】杜教筛（Sum）

我们假设要求前缀和的积性函数为$f(n)$，设前缀和$S(n)=\sum\limits_{i=1}^nf(i)$。

我们找一个合适的函数$g(n)$，我们看看$f,g$两个函数的狄利克雷卷积的前缀和。

$\sum\limits_{i=1}^n(f*g)(i)$

$=\sum\limits_{i=1}^n\sum\limits_{d|i}f(\dfrac{i}{d})g(d)$

先枚举$d$再枚举$i$：

$\sum\limits_{d=1}^n\sum\limits_{d|i\&i\leq n}f(\dfrac{i}{d})g(d)$

把$i$除以$d$

$\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)g({d})$

提出$g(d)$

$\sum\limits_{d=1}^ng(d)\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(i)=\sum\limits_{d=1}^ng(d)S(\lfloor\frac{n}{d}\rfloor)$

也就是说$\sum\limits_{i=1}^n(f*g)(i)=\sum\limits_{i=1}^ng(i)S(\lfloor\frac{n}{i}\rfloor)$

再看$g(1)S(n)$的值：

$g(1)S(n)=\sum\limits_{i=1}^ng(i)S(\lfloor\dfrac{n}{i}\rfloor)-\sum\limits_{i=2}^ng(i)S(\lfloor\dfrac{n}{i}\rfloor)=\sum\limits_{i=1}^n(f*g)(i)-\sum\limits_{i=2}^ng(i)S(\lfloor\dfrac{n}{i}\rfloor)$

只要我们找到$g$使得可以$O(1)$算出$f*g$的前缀和、$g$的前缀和，我们靠这个式子也就可以递归。

递归前，我们可以预处理$n^{\frac{2}{3}}$内的前缀和，再用unordered_map记忆化，就能达到$O(n^{\frac{2}{3}})$的复杂度。

注：unordered_map为C++11的STL

本题，当$f$为$\mu$时，找的$g$为$I$，则$f*g=\epsilon,g(1)=1$

当$f$为$\varphi$时，找的$g$为$I$，则$f*g=id,g(1)=1$

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
const int MAXN=5000005;
bool isp[MAXN+1];
int phi[MAXN+1],mu[MAXN+1],pcnt,n,prime[MAXN+1];
ll ans,smu[MAXN+1],sphi[MAXN+1];
unordered_map<int,ll>Sphi,Smu;
ll Get_smu(unsigned int n) {//mu
	if(n<=MAXN) return smu[n];
	if(Smu[n]) return Smu[n];
	ull ans=1ll;
	for(unsigned int l=2,r; l<=n; l=r+1) {//整除分块
		r=n/(n/l);
		ans-=(r-l+1ll)*Get_smu(n/l);//递归
	}
	return Smu[n]=ans;
} //mu*I=e
ll Get_sphi(unsigned int n) {//phi
	if(n<=MAXN)return sphi[n];
	if(Sphi[n]) return Sphi[n];
	ull ans=n*(n+1ll)>>1;
	for(unsigned int l=2,r; l<=n; l=r+1) {
		r=n/(n/l);
		ans-=(r-l+1ll)*Get_sphi(n/l);
	}
	return Sphi[n]=ans;
}//phi*I=id
int main() {
	sphi[1]=smu[1]=mu[1]=phi[1]=1;
	for(int i=2; i<=MAXN; ++i) {//筛出一些mu,phi并计算前缀和
		if(!isp[i]) {
			prime[++pcnt]=i;
			mu[i]=-1;
			phi[i]=i-1;
		}
		for(int j=1,x; j<=pcnt&&(x=i*prime[j])<=MAXN; ++j) {
			isp[x]=1;
			if(i%prime[j]) {
				mu[x]=-mu[i];
				phi[x]=phi[i]*(prime[j]-1);
			} else {
				mu[x]=0;
				phi[x]=phi[i]*prime[j];
				break;
			}
		}
		sphi[i]=sphi[i-1]+phi[i];
		smu[i]=smu[i-1]+mu[i];//预处理一些前缀和
	}
	int t,n;
	scanf("%d",&t);
	while(t--) {
		scanf("%d",&n);
		printf("%lld %lld\n",Get_sphi(n),Get_smu(n));
	}
	return 0;
}

一道题

P3768 简单的数学题

先用$d$枚举$\gcd(i,j)$

$\sum\limits_{d=1}^n\sum\limits_{i=1}^n\sum\limits_{j=1}^n[\gcd(i,j)=d]ijd$

把$i,j$除以$d$

$\sum\limits_{d=1}^n\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor}[\gcd(i,j)=1] (id)(jd)d$

把三个$d$都提出来

$\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} ij[\gcd(i,j)=1]$

再把$\epsilon=\mu*I$代入

$\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} ij\sum\limits_{k|\gcd(i,j)}\mu(k)$

前面的题也做过，$k|\gcd(i,j)$等价于$k|i\&k|j$

$\sum\limits_{d=1}^nd^3\sum\limits_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{n}{d}\rfloor} ij\sum\limits_{k|i\&k|j}\mu(k)$

把$k$提前

$\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\sum\limits_{k|i\&i\leq\lfloor\frac{n}{d}\rfloor}\sum\limits_{k|j\&j\leq\lfloor\frac{n}{d}\rfloor}ij\mu(k)$

把$i,\mu(k)$提出

$\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum\limits_{k|i\&i\leq\lfloor\frac{n}{d}\rfloor}i\sum\limits_{k|j\&j\leq\lfloor\frac{n}{d}\rfloor}j$

后面的$i,j$都除以$k$

$\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}ik\sum\limits_{j=1}^{\lfloor\frac{n}{dk}\rfloor} jk$

把两个$k$都提出

$\sum\limits_{d=1}^nd^3\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}\mu(k)k^2\sum\limits_{i=1}^{\lfloor\frac{n}{dk}\rfloor}i\sum\limits_{j=1}^{\lfloor\frac{n}{dk}\rfloor} j$

把$d$乘进去，设$S(n)=\sum\limits_{i=1}^ni$

$\sum\limits_{d=1}^n\sum\limits_{k=1}^{\lfloor\frac{n}{d}\rfloor}d^3k^2\mu(k)S^2(\lfloor\dfrac{n}{dk}\rfloor)$

先枚举$T=dk$，再枚举$d|T$，$k$就是$\dfrac{n}{d}$

$\sum\limits_{T=1}^n\sum\limits_{d|T}dT^2\mu(\dfrac{n}{d})S^2(\lfloor\dfrac{n}{T}\rfloor)$

提出$T^2S^2(\lfloor\dfrac{n}{T}\rfloor)$

$\sum\limits_{T=1}^nT^2S^2(\lfloor\dfrac{n}{T}\rfloor)\sum\limits_{d|T}d\mu(\dfrac{n}{d})$

后面的$\sum\limits_{d|T}d\mu(\dfrac{n}{d})$就是$id*\mu$

由$I*\varphi=id$莫比乌斯反演得$id*\mu=\varphi$

$\sum\limits_{T=1}^nT^2S^2(\lfloor\dfrac{n}{T}\rfloor)\varphi(T)=\sum\limits_{T=1}^nS^2(\lfloor\dfrac{n}{T}\rfloor)T^2\varphi(T)$

设$f(n)=T^2\varphi(T)$，若能求出它的前缀和，就可以整出分块了。

考虑用杜教筛。

设$g(n)=n^2$，则$f*g=\sum\limits_{d|n}d^2\varphi(d)\times(\dfrac{n}{d})^2=n^2\sum\limits_{d|n}\varphi(d)$

后面的部分就是$\varphi*I$，是等于$id$的

即$f*g=id^3$，前缀和为$(\dfrac{n(n+1)}{2})^2$。

而且$g$的前缀和为$\dfrac{n(n+1)(2n+1)}{6}$，然后就可以杜教筛了。

注：代码中，我手写了哈希表。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MAXN=5000005;
namespace Hash {//图版哈希表
	const int mo=233333;
	struct hash_table {
		int head[mo],htcnt;
		struct hashtable {
			ll id,val;
			int nxt;
		} ht[mo];
		void add(int u,ll v,ll w) {
			++htcnt;
			ht[htcnt].id=v;
			ht[htcnt].val=w;
			ht[htcnt].nxt=head[u];
			head[u]=htcnt;
		}
		void clr() {
			memset(head,0,sizeof(head));
			memset(&ht,0,sizeof(&ht));
			htcnt=0;
		}
		void insert(ll x,ll i) {
			int k=x%mo;
			add(k,x,i);
		}
		int query(ll x) {
			for(int i=head[x%mo]; i; i=ht[i].nxt)
				if(ht[i].id==x)
					return ht[i].val;
			return -1;
		}
	};
} using namespace Hash;
int prime[MAXN+1],pcnt,phi[MAXN+1],mod;
int sphi[MAXN+1],inv6;//inv6为6的逆元
hash_table Sphi;
ll n;
int Getsphi(ll m) {//查询
	if(m<=MAXN) return sphi[m];
	return Sphi.query(m);
}
int S(int n) {//1~n和
	return (n*(n+1ll)>>1)%mod;
}
int Sp(int n) {//g的前缀和
	return n*(n+1ll)%mod*(n<<1|1ll)%mod*inv6%mod;
}
void Get(ll m,ll d) {//杜教筛
	if(m<=MAXN||~Sphi.query(m))return ;//查询过就不管
	ll ans=S(m%mod);
	ans=ans*ans%mod;//f*g前缀和
	for(ll l=2,r; l<=m; l=r+1) {
		r=m/(m/l);
		Get(m/l,d*l);//先递归
		ans=(ans+mod-(Sp(r%mod)+0ll+mod-Sp((l-1)%mod))%mod*Getsphi(m/l)%mod)%mod;
	}
	Sphi.insert(m,ans);//记忆化
}
int Get_sphi(ll m) {
	Get(m,1);
	return Getsphi(m);//查询
}
int fastpow(int a,int k=mod-2) {//快速幂求逆元
	int base=1;
	while(k) {
		if(k&1)base=base*1ll*a%mod;
		a=a*1ll*a%mod;
		k>>=1;
	}
	return base;
}
int main() {
	scanf("%d%lld",&mod,&n);
	inv6=fastpow(6);
	sphi[1]=phi[1]=1;
	for(int i=2; i<=MAXN; ++i) {//线性筛前缀和
		if(!phi[i]) {
			prime[++pcnt]=i;
			phi[i]=i-1;
		}
		for(int j=1,x; j<=pcnt&&(x=i*prime[j])<=MAXN; ++j) {
			if(i%prime[j]) {
				phi[x]=phi[i]*(prime[j]-1);
			} else {
				phi[x]=phi[i]*prime[j];
				break;
			}
		}
		sphi[i]=(sphi[i-1]+phi[i]*1ll*i%mod*i)%mod;
	}
	ll tmp1,ans=0;
	for(ll l=1,r; l<=n; l=r+1) {//整除分块求答案
		r=n/(n/l);
		tmp1=S((n/l)%mod);
		tmp1=tmp1*tmp1%mod;
		ans=(ans+tmp1*(Get_sphi(r)+0ll+mod-Get_sphi(l-1))%mod)%mod;
	}
	printf("%lld\n",ans);
	return 0;
}