胡小兔的杜教篩學習筆記

很久沒寫數論題，今天在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$的前綴和：

$$\begin{align*} \sum_{i = 1}^{n} (fg)(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(\lfloor\frac{n}{d}\rfloor) \ &= g(1)S(n) + \sum_{d = 2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) \&= S(n) + \sum_{d = 2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) \end{align}$$

那麼移項可得

$$S(n) = \sum_{i = 1}^{n} (f*g)(i) - \sum_{d = 2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor)$$

若是這個$(f*g)(n)$的前綴和很是好求的話，$S(n)$就能夠數論分塊（看！那有個$\lfloor\frac{n}{d}\rfloor$！）而後遞歸求解下去了。（這裏須要用個哈希表之類的東西記憶化。）

栗子

求$\varphi$的前綴和

$\varphi * 1 = Id$，那麼

$$\begin{align*}S(n) &= \sum_{i = 1}^{n} (\varphi1)(i) - \sum_{d = 2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) \ &= \sum_{i = 1}^{n} i - \sum_{d = 2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) \ &= \frac{n (n + 1)}{2}- \sum_{d = 2}^{n}g(d)S(\lfloor\frac{n}{d}\rfloor) \end{align}$$

遞歸求解便可。

個人代碼：

#include <cstdio>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <iostream>
#include <map>
#define space putchar(' ')
#define enter putchar('\n')
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;
}