杜教筛用于计算某一类积性函数的前缀和。

思想

构造 Dirchlet 卷积，把一个困难的前缀和转化为两个相对简单的前缀和。

形式化地说，我们有数论函数 f,g，记 F(n) = \sum\limits_{i = 1}^n f(i) ，那么有：

\begin{aligned} \sum\limits_{i = 1}^n f*g(i) &= \sum\limits_{i=1}^n \sum_{x\times y = i}f(x)\times g(y) \\ &= \sum\limits_{y=1}^n g(y)\sum\limits_{x\times y \leq n} f(x)\\ &= \sum\limits_{y=1}^n g(y)\times F(\lfloor \frac{n}{y}\rfloor)\\ \end{aligned}

也就是说，我们有：

g(1)F(n) = \sum\limits_{i = 1}^n f*g(i) - \sum\limits_{y=2}^n g(y)\times F(\lfloor \frac{n}{y}\rfloor)

这里如果我们假设 f *g 的前缀和和 g 的区间和都可以 O(1) 计算（这也是我们应用杜教筛的一个条件），那么计算 F，我们需要的就是计算后面 \sum\limits_{y=2}^n g(y)\times F(\lfloor \frac{n}{y}\rfloor) ，那么这个显然可以数论分块，如果直接数论分块，我们递归计算 F(\lfloor \frac{n}{y}\rfloor)就是 \mathcal{O}(\sqrt n) 次的了。

好，我们解决了单次计算的循环次数问题，那么我们发现每次递归的 F(x) 的取值实质上是 1,2\dots \sqrt n ，这个如果直接记忆化暴力计算，通过调和级数我们可以知道复杂度是 \mathcal{O}(n^\frac{3}{4}) 的。

但是我们发现，对于较小的 n 通过 \mathcal{O}\sqrt n 来计算是血亏的，我们运用一个类似于根号分治的思想，对于 1\dots n^c 的答案用 \mathcal{O}(n^c) 的时间预处理出来，通过一些积分魔术我们可以得到当 c = \frac{2}{3} 的时候复杂度最优秀，为 \mathcal{O}(n^{ \frac{2}{3}}) 。

记忆化的部分，我们可以直接记录一个 ans[x] 代表 \frac{n}{x} 时候的答案。对于计算 f *g 以及 g 的前缀和，我们发现只需要在 \mathcal{O}(\sqrt n) 的时间内算完就不影响复杂度，实际上你也可以套一层杜教筛（

例题

P4213 【模板】杜教筛（Sum)

给定一个正整数，求

ans_1=\sum_{i=1}^n\varphi(i)

ans_2=\sum_{i=1}^n \mu(i)​

1 \leq n \lt 2^{31}

对于 \varphi ，我们有 \varphi(n) * 1 = id，这里 1 就是我们的 g，id 就是我们的 f*g，这两个的求和是非常简单的。

对于 \mu，我们有 \mu(n) * 1 = \epsilon(n) = [n = 1]，这同样非常容易求和。

参考代码：

const int N = 5e6 + 10;
const int M = 5e3 + 10;

ll phi[N],maxn,maxm,n,mu[N],pri[N],tot,prephi[M],premu[M];
bool vis[N];

void init() {
    phi[1] = mu[1] = 1;
    for (int i = 2;i <= N-10;++i) {
        if (!vis[i]) pri[++tot] = i,mu[i] = -1,phi[i] = i - 1;
        for (int j = 1;j <= tot && pri[j] * i <= N - 10;++j) {
            vis[pri[j] * i] = 1;
            if (i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
            phi[i * pri[j]] = phi[i] * phi[pri[j]];
            mu[i * pri[j]] = -mu[i];
        }
    }
    for (int i = 1;i <= N-10;++i) mu[i] += mu[i-1],phi[i] += phi[i-1];
}//在O(n^c)的时间内进行预处理。

ll getphi(ll x) {
    if (x < N-10) return phi[x];
    ll now = n / x;
    if (prephi[now]) return prephi[now];
    //注意这边，我们记忆化的实际上是n / x的值，因为我们递归调用到的总是一个形如n/x的形式，所以这里的数组可以开的相对小一点
    ll ans = x * (x + 1) / 2;//\sum f*g
    for (ll l = 2,r;l <= x;l = r + 1) {
        r = x / (x / l);
        ans -= (r - l + 1) * getphi(x / l);//-\sum g
    }
    return prephi[now] = ans;//记得记忆化
}

ll getmu(ll x) {
    if (x < N-10) return mu[x];
    ll now = n / x;
    if (premu[now]) return premu[now];
    ll ans = 1;
    for (ll l = 2,r;l <= x;l = r + 1) {
        r = x / (x / l);
        ans -= (r - l + 1) * getmu(x / l);
    }
    return premu[now] = ans;
}

P3768 简单的数学题

求

\left(\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\right) \bmod p

n \leq 10^{10}

\begin{aligned} &\sum_{i=1}^n\sum_{j=1}^n ij \gcd(i,j)\\ &= \sum_{i=1}^n\sum_{j=1}^n\sum_{d=1}^nd ij [\gcd(i,j) = d]\\ &= \sum_{i=1}^{n/d}\sum_{j=1}^{n/d}\sum_{d=1}^nd^3 ij [\gcd(i,j) = 1]\\ &= \sum_{d = 1} ^ n d^3 \sum\limits_{i=1}^{n/d}\sum\limits_{j=1}^{n/d}ij \sum\limits_{k | i,k | j} \mu(k)\\ &= \sum_{d = 1} ^ n d^3 \sum\limits_{k} k^2\mu(k) \sum\limits_{i=1}^{n/dk}\sum\limits_{j=1}^{n/dk}ij \\ &= \sum_{d = 1} ^ n d^3 \sum\limits_{k} k^2\mu(k) (\sum_{i=1}^{n/dk} i) ^2 \\ &\text{Let sum(x)} = \sum\limits_{i=1}^x i,T = dk\\ &= \sum_{T=1}^n sum(\frac{n}{T})^2 T^2\sum_{d|T} d\mu(\frac{T}{d}) \end{aligned}

那么我们要求的就是 T^2\sum_{d|T}d\mu(\frac{T}{d})，然后我们发现后面这个东西再 Dirchlet 卷积一下有 \sum_{d|T}d\mu(\frac{T}{d}) = \varphi(T)，所以就变成了 T^2\varphi(T) ，太神奇啦家人们！

那么我们总的式子就变成了：

\sum_{T=1}^n sum(\frac{n}{T})^2 T^2\varphi(T)

别急，按杜教筛套路继续推下去，我们令 f(i) = i^2 \varphi(i) ， F(i) = \sum\limits_{i=1}^n f(i)

我们来考虑这个 f 怎么卷比较可做！因为有 \varphi * 1 = id，所以我们直接设 g(x) = x^2 ，这样我们就易得 f * g (x) = x ^3。这就好起来了！所以我们的 F 就好做了：

\begin{aligned} g(1)F(n) &= \sum\limits_{i = 1}^n f*g(i) -\sum\limits_{i=2}^n g(i)\times F(\lfloor \frac{n}{i}\rfloor) \\ F(n) &= \sum\limits_{i=1}^n i^3 - \sum\limits_{i=2}^n i^2 \times F(\lfloor \frac{n}{i}\rfloor) \end{aligned}

这个是典型的杜教筛式，所以我们就数论分块里面再套一个杜教筛就OK了。

代码：

const int N = 8000000 + 10;
const int M = 1e5 + 10;

ll phi[N],maxn,maxm,n,mu[N],pri[N],tot,premu[M];
ll mod;
bool vis[N];
unordered_map <ll,ll> prephi;

void init() {
    phi[1] = mu[1] = 1;
    for (int i = 2;i <= N-10;++i) {
        if (!vis[i]) pri[++tot] = i,mu[i] = -1,phi[i] = i - 1;
        for (int j = 1;j <= tot && pri[j] * i <= N - 10;++j) {
            vis[pri[j] * i] = 1;
            if (i % pri[j] == 0) {
                phi[i * pri[j]] = phi[i] * pri[j];
                break;
            }
            phi[i * pri[j]] = phi[i] * phi[pri[j]];
            mu[i * pri[j]] = -mu[i];
        }
    }
    for (ll i = 1;i <= N-10;++i) phi[i] = (phi[i-1] + phi[i] * i % mod * i  % mod) % mod;
}

ll inv2,inv6;

ll sum2(ll x) {
    x %= mod;
    return x * (x + 1) % mod * inv2 % mod; 
}

ll sum3(ll x) {
    x %= mod;
    return x * (x + 1) % mod * (2 * x % mod + 1) %mod *inv6 %mod;
}

ll dusieve(ll x) {
    //printf("%lld ",x);
    if (x < N-10) return phi[x];
    ll now = n / x;
    if (prephi[x]) return prephi[x];
    ll ans = sum2(x) * sum2(x) % mod;
    for (ll l = 2,r;l <= x;l = r + 1) {
        r = x / (x / l);
        ans = ((ans - (((sum3(r) - sum3(l-1)) % mod + mod) % mod * dusieve(x / l))) % mod + mod) % mod;
    }
    return prephi[x] = (ans + mod) % mod;
}

ll ffpow(ll a,ll b) {
    a %= mod;
    ll ans = 1;
    for (;b;b >>= 1) {
        if (b & 1) ans = ans * a % mod;
        a = a * a % mod;
    }
    return ans;
}

signed main() {
    read(mod,n);
    inv2 = ffpow(2,mod - 2);inv6 = ffpow(6,mod - 2);
    init();
    ll ans = 0;
    for (ll l = 1,r;l <= n;l = r + 1) {
        r = n / (n / l);
        ll a1 = sum2(n / l) * sum2(n / l) % mod;
        ll a2 = ((dusieve(r) - dusieve(l - 1)) % mod + mod) % mod;
        ans = (ans + (a1 * a2) % mod) % mod;
    }
    if (ans < 0) ans += mod;
    printf("%lld\n",ans % mod);
    return 0;
}