还是发出来提醒下自己,免得又给忘了...
杜教筛
- 一种可以在 (O(n^{frac 2 3})) 时间内解决积性函数前缀和的操作.
- 求 (S(n)=sum_{i=1}^n f(i),f) 为积性函数.
- 构造两个积性函数 (g,h) ,使得 (f*g=h).
[egin{align*}
sum_{i=1}^n h(i)&=
sum_{i=1}^n sum_{d|n}f(d)*g(frac i d)\&=
sum_{d=1}^n g(d) cdot sum_{i=1}^{lfloor frac n d
floor} f(i)\&=
g(1)cdot S(n)+sum_{d=2}^n g(d)cdot S(lfloor frac n d
floor)\
herefore g(1)cdot S(n)&=sum_{i=1}^n h(i)-sum_{d=2}^n g(d)cdot S(lfloor frac n d
floor).
end{align*}
]
-
这样,只要我们选取的 (g,h) 的前缀和很好求,计算 (S(n)) 就只有 (O(n^{frac 2 3})) 的复杂度.
-
实现时可以暴力算前 (sqrt n) 项的 (S) ,按照上面的式子递归计算,并用 (hash) 或 (map) 记忆化.
-
常用的几个 (f*g=h) :
-
(1. mu * I=epsilon).
-
(2. varphi * I=id).
-
(3. id * mu=varphi).
-
(4. f *epsilon=f).
-
另外,若 (f(i)=mu(i) cdot i^k,varphi(i) cdot i^k) , 也十分好构造.直接用原来对应的 (gcdot i^k) ,这样得到的 (h) 会多乘一个 (n^k) ,而这个东西的前缀和可以背公式((k) 小)或插值搞出来.((k) 大).
-
例:
简单的数学题
-
题意:求(sum_{i=1}^n sum_{j=1}^n ij*gcd(i,j)) , (nleq 1e10) ,答案对一个质数 (p) 取模.
-
这道题用 (varphi) 反演会简单一些.
[egin{align*}
sum_{i=1}^n sum_{j=1}^n ij*gcd(i,j)&=
sum_{i=1}^n sum_{j=1}^n ij sum_{k|gcd(i,j)} varphi(k)\&=
sum_{k=1}^n varphi(k) sum_{k|i}sum_{k|j}ij\&=
sum_{k=1}^n varphi(k) cdot k^2 cdot (sum_{i=1}^{lfloor frac n k
floor} i)^2.
end{align*}
]
- 现在外层整除分块是 (sqrt n) 的,考虑到 (n) 的大小,需要快速计算出(sum_{k=l}^r varphi(k) cdot k^2) ,需杜教筛解决.
#include"bits/stdc++.h"
using namespace std;
typedef long long LoveLive;
inline LoveLive read()
{
LoveLive out=0,fh=1;
char jp=getchar();
while ((jp>'9'||jp<'0')&&jp!='-')
jp=getchar();
if (jp=='-')
{
fh=-1;
jp=getchar();
}
while (jp>='0'&&jp<='9')
{
out=out*10+jp-'0';
jp=getchar();
}
return out*fh;
}
const int MAXN=4600010;
int p;
LoveLive n;
LoveLive add(LoveLive a,LoveLive b)
{
return (a + b) % p;
}
LoveLive mul(LoveLive a,LoveLive b)
{
return (a * b) % p;
}
LoveLive fpow(LoveLive a,LoveLive b)
{
a%=p;
b%=p-1;
int res=1;
while(b)
{
if(b&1)
res=mul(res,a);
a=mul(a,a);
b>>=1;
}
return res;
}
int prime[MAXN],cnt=0,ism[MAXN];
LoveLive phi[MAXN],sumphi[MAXN];
LoveLive inv6,inv2;
void init(int N)
{
inv6=fpow(6,p-2);
inv2=fpow(2,p-2);
ism[1]=phi[1]=1;
for(int i=2;i<=N;++i)
{
if(!ism[i])
prime[++cnt]=i,phi[i]=i-1;
for(int j=1;j<=cnt && i*prime[j]<=N;++j)
{
ism[i*prime[j]]=1;
if(i%prime[j]==0)
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
else
phi[i*prime[j]]=phi[i]*(prime[j]-1);
}
}
for(int i=1;i<=N;++i)
sumphi[i]=add(sumphi[i-1],mul(mul(i,i),phi[i]));
}
LoveLive s2(LoveLive x)
{
x%=p;
LoveLive res=x+1;
res=mul(res,x);
res=mul(res,2*x+1);
res=mul(res,inv6);
return res;
}
LoveLive s3(LoveLive x)
{
x%=p;
LoveLive tmp=mul(x,x+1);
tmp=mul(tmp,inv2);
return mul(tmp,tmp);
}
map<LoveLive,LoveLive> mp;
LoveLive calc(LoveLive x)
{
if(x<=MAXN-10)
return sumphi[x];
if(mp[x])
return mp[x];
LoveLive res=s3(x);
for(LoveLive l=2,r;l<=x;l=r+1)
{
r=x/(x/l);
LoveLive tmp=add(s2(r),p-s2(l-1));
tmp=mul(tmp,calc(x/l));
res=add(res,p-tmp);
}
return mp[x]=res;
}
int main()
{
p=read();
n=read();
init(MAXN-10);
LoveLive ans=0;
for(LoveLive l=1,r;l<=n;l=r+1)
{
r=n/(n/l);
LoveLive tmp=add(calc(r),p-calc(l-1));
tmp=mul(tmp,s3(n/l));
ans=add(ans,tmp);
}
printf("%lld
",ans);
return 0;
}