好久没写数论题，今天在51nod抓了一道，发现自己早就把杜教筛忘得一干二净啦~ 所以今天我把杜教筛学习笔记整理一下，防止以后再次忘记 =v=
[Warning] 杜教筛复杂度证明我暂时还不会 >_< 我会抓紧时间学的
前置技能
如果你已经了解了以下某些部分的内容，请跳过该部分。
积性函数
(f)是一个数论函数，若对于任意两个互质的数(a)和(b)有(f(a*b) = f(a)*f(b))，则称(f)是积性函数。
（本文可能用到的）常见积性函数
- (1(n) = 1)
- (Id(n) = n)
- (e(n) = [n == 1])
- (d(n)) 因数个数
- (sigma(n)) 因数之和
- (mu(n)) 莫比乌斯函数
- (varphi(n)) 欧拉函数
狄利克雷卷积
两个数论函数(f)和(g)对于(n)的狄利克雷卷积((f * g)(n))是(sum_{d | n}f(d)g(frac{n}{d}))。
两个积性函数的狄利克雷卷积也是积性函数（此时(n)是传入函数中的参数），即：若(f)和(g)是积性函数，(h(n) = sum_{d | n}f(d)g(frac{n}{d}))，则(h)是积性函数。
常见狄利克雷卷积
- (mu * 1 = e)，即(sum_{d|n}mu(d) = [n == 1])
- (mu * Id = varphi)，即(sum_{d|n}mu(d)frac{n}{d} = varphi(n))
- (varphi * 1 = Id)，即(sum_{d|n}varphi(d) = n)
杜教筛
杜教筛是用来在(O(n^frac{2}{3}))时间内求一些积性函数的前缀和的。
假设我们要求前缀和的函数是(f)，它的前缀和为(S)， 即(S(n) = sum_{i = 1}^{n}f(i))。
我们选择一个积性函数(g)用来和(f)卷一卷。（(g)是哪来的？怎么求？……不知道，因题而定，基本靠猜 = =）
然后求一下(f*g)的前缀和：
[egin{align*} sum_{i = 1}^{n} (f*g)(i) &= sum_{i = 1}^{n}sum_{d|i}g(d)f(frac{i}{d}) \&= sum_{d = 1}^{n}g(d)sum_{d|i} f(frac{i}{d}) \ &= sum_{d = 1}^ng(d)S(lfloorfrac{n}{d}
floor) \ &= g(1)S(n) + sum_{d = 2}^{n}g(d)S(lfloorfrac{n}{d}
floor) \&= S(n) + sum_{d = 2}^{n}g(d)S(lfloorfrac{n}{d}
floor) end{align*}
]
那么移项可得
[S(n) = sum_{i = 1}^{n} (f*g)(i) - sum_{d = 2}^{n}g(d)S(lfloorfrac{n}{d}
floor)
]
如果这个((f*g)(n))的前缀和非常好求的话，(S(n))就可以数论分块（看！那有个(lfloorfrac{n}{d} floor)！）然后递归求解下去了。（这里需要用个哈希表之类的东西记忆化。）
栗子
求(varphi)的前缀和
(varphi * 1 = Id)，那么
[egin{align*}S(n) &= sum_{i = 1}^{n} (varphi*1)(i) - sum_{d = 2}^{n}g(d)S(lfloorfrac{n}{d}
floor) \ &= sum_{i = 1}^{n} i - sum_{d = 2}^{n}g(d)S(lfloorfrac{n}{d}
floor) \ &= frac{n (n + 1)}{2}- sum_{d = 2}^{n}g(d)S(lfloorfrac{n}{d}
floor) end{align*}
]
递归求解即可。
我的代码：
#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <map>
#define space putchar(' ')
#define enter putchar('
')
typedef long long ll;
using namespace std;
template <class T>
void read(T &x){
char c;
bool op = 0;
while(c = getchar(), c < '0' || c > '9')
if(c == '-') op = 1;
x = c - '0';
while(c = getchar(), c >= '0' && c <= '9')
x = x * 10 + c - '0';
if(op) x = -x;
}
template <class T>
void write(T x){
if(x < 0) putchar('-'), x = -x;
if(x >= 10) write(x / 10);
putchar('0' + x % 10);
}
const int N = 5000000, P = 1000000007, inv2 = 500000004, S = 1333333;
ll n, m;
int phi[N], sum[N], prime[N], tot, idx;
int adj[S], nxt[S], val[S];
ll num[S];
bool notprime[N];
void add(ll u, ll v){
num[++idx] = u;
val[idx] = v;
nxt[idx] = adj[u % S];
adj[u % S] = idx;
}
ll find(ll u){
for(int e = adj[u % S]; e; e = nxt[e])
if(num[e] == u) return val[e];
return -1;
}
void init(){
phi[1] = 1;
for(ll i = 2; i <= m; i++){
if(!notprime[i]) prime[++tot] = i, phi[i] = i - 1;
for(int j = 1; j <= tot && i * prime[j] <= m; j++){
notprime[i * prime[j]] = 1;
if(i % prime[j]) phi[i * prime[j]] = phi[i] * (prime[j] - 1);
else{
phi[i * prime[j]] = phi[i] * prime[j];
break;
}
}
}
for(ll i = 1; i <= m; i++)
sum[i] = (sum[i - 1] + phi[i]) % P;
}
ll calc(ll x){
if(x <= m) return sum[x];
ll ret = find(x);
if(ret != -1) return ret;
ret = (x + 1) % P * (x % P) % P * inv2 % P;
for(ll i = 2, last; i <= x; i = last + 1){
last = x / (x / i);
ret = (ret - (last - i + 1) % P * calc(x / i)) % P;
}
ret = (ret + P) % P;
add(x, ret);
return ret;
}
int main(){
read(n), m = pow(n, 2.0 / 3) + 0.5;
init();
write(calc(n)), enter;
return 0;
}