杜教筛
前置知识：莫比乌斯反演
杜教筛
杜教筛可以在非线性时间内求积性函数前缀和。
我们假设要求前缀和的积性函数为\(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;
}
一道题
先用\(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;
}
