跳转至

# 杜教筛

杜教筛被用于处理一类数论函数的前缀和问题．对于数论函数 𝑓 f，杜教筛可以在低于线性时间的复杂度内计算 𝑆(𝑛) =∑𝑛 𝑖=1 𝑓(𝑖)S(n)=∑i=1 n f(i)．

## 算法思想

我们想办法构造一个 𝑆(𝑛)S(n) 关于 𝑆(⌊𝑛 𝑖⌋)S(⌊n i⌋) 的递推式．

对于任意一个数论函数 𝑔 g，必满足：

𝑛∑𝑖=1(𝑓∗𝑔)(𝑖)=𝑛∑𝑖=1∑𝑑∣𝑖 𝑔(𝑑)𝑓(𝑖 𝑑)=𝑛∑𝑖=1 𝑔(𝑖)𝑆(⌊𝑛 𝑖⌋)∑i=1 n(f∗g)(i)=∑i=1 n∑d∣i g(d)f(i d)=∑i=1 n g(i)S(⌊n i⌋)

其中 𝑓 ∗𝑔 f∗g 为数论函数 𝑓 f 和 𝑔 g 的 狄利克雷卷积．

略证

𝑔(𝑑)𝑓(𝑖 𝑑)g(d)f(i d) 就是对所有 𝑖 ≤𝑛 i≤n 的做贡献，因此变换枚举顺序，枚举 𝑑 d,𝑖 𝑑 i d（分别对应新的 𝑖,𝑗 i,j）

𝑛∑𝑖=1∑𝑑∣𝑖 𝑔(𝑑)𝑓(𝑖 𝑑)=𝑛∑𝑖=1⌊𝑛/𝑖⌋∑𝑗=1 𝑔(𝑖)𝑓(𝑗)=𝑛∑𝑖=1 𝑔(𝑖)⌊𝑛/𝑖⌋∑𝑗=1 𝑓(𝑗)=𝑛∑𝑖=1 𝑔(𝑖)𝑆(⌊𝑛 𝑖⌋)∑i=1 n∑d∣i g(d)f(i d)=∑i=1 n∑j=1⌊n/i⌋g(i)f(j)=∑i=1 n g(i)∑j=1⌊n/i⌋f(j)=∑i=1 n g(i)S(⌊n i⌋)

那么可以得到递推式：

𝑔(1)𝑆(𝑛)=𝑛∑𝑖=1 𝑔(𝑖)𝑆(⌊𝑛 𝑖⌋)−𝑛∑𝑖=2 𝑔(𝑖)𝑆(⌊𝑛 𝑖⌋)=𝑛∑𝑖=1(𝑓∗𝑔)(𝑖)−𝑛∑𝑖=2 𝑔(𝑖)𝑆(⌊𝑛 𝑖⌋)g(1)S(n)=∑i=1 n g(i)S(⌊n i⌋)−∑i=2 n g(i)S(⌊n i⌋)=∑i=1 n(f∗g)(i)−∑i=2 n g(i)S(⌊n i⌋)

假如我们可以构造恰当的数论函数 𝑔 g 使得：

1. 可以快速计算 ∑𝑛 𝑖=1(𝑓 ∗𝑔)(𝑖)∑i=1 n(f∗g)(i)；
2. 可以快速计算 𝑔 g 的前缀和，以用数论分块求解 ∑𝑛 𝑖=2 𝑔(𝑖)𝑆(⌊𝑛 𝑖⌋)∑i=2 n g(i)S(⌊n i⌋)．

则我们可以在较短时间内求得 𝑔(1)𝑆(𝑛)g(1)S(n)．

注意

无论数论函数 𝑓 f 是否为积性函数，只要可以构造出恰当的数论函数 𝑔 g, 便都可以考虑用杜教筛求 𝑓 f 的前缀和．

如考虑 𝑓(𝑛) =i 𝜑(𝑛)f(n)=i φ(n), 显然 𝑓 f 不是积性函数，但可取 𝑔(𝑛) =1 g(n)=1, 从而：

𝑛∑𝑘=1(𝑓∗𝑔)(𝑘)=i 𝑛(𝑛+1)2∑k=1 n(f∗g)(k)=i n(n+1)2

计算 ∑𝑘≤𝑚(𝑓 ∗𝑔)(𝑘)∑k≤m(f∗g)(k) 和 ∑𝑘≤𝑚 𝑔(𝑘)∑k≤m g(k) 的时间复杂度均为 𝑂(1)O(1), 故可以考虑使用杜教筛．

## 时间复杂度

令 𝑅(𝑛) ={⌊𝑛 𝑘⌋:𝑘=2,3,…,𝑛}R(n)={⌊n k⌋:k=2,3,…,n}．利用数论分块的 性质 可知，对任意的 𝑚 ∈𝑅(𝑛)m∈R(n)，都有 𝑅(𝑚) ⊆𝑅(𝑛)R(m)⊆R(n)．也就是说，使用记忆化之后，只需要对所有 𝑘 ∈𝑅(𝑛)k∈R(n) 计算一次 𝑆(𝑘)S(k) 就可以得到 𝑅(𝑛)R(n) 的值．而这些点的数目 |𝑅(𝑛)| =𝑂(√𝑛)|R(n)|=O(n)．

设计算 ∑𝑛 𝑖=1(𝑓 ∗𝑔)(𝑖)∑i=1 n(f∗g)(i) 和 ∑𝑛 𝑖=1 𝑔(𝑖)∑i=1 n g(i) 的时间复杂度均为 𝑂(1)O(1). 设计算 𝑆(𝑛)S(n) 的时间复杂度为 𝑇(𝑛)T(n), 则：

𝑇(𝑛)=∑𝑘∈𝑅(𝑛)𝑇(𝑘)=Θ(√𝑛)+⌊√𝑛⌋∑𝑘=1 𝑂(√𝑘)+⌊√𝑛⌋∑𝑘=2 𝑂(√𝑛 𝑘)=𝑂(∫√𝑛 0(√𝑥+√𝑛 𝑥)d 𝑥)=𝑂(𝑛 3/4).T(n)=∑k∈R(n)T(k)=Θ(n)+∑k=1⌊n⌋O(k)+∑k=2⌊n⌋O(n k)=O(∫0 n(x+n x)d x)=O(n 3/4).

若我们可以预处理出一部分 𝑆(𝑘)S(k), 其中 𝑘 =1,2,…,𝑚 k=1,2,…,m，𝑚 ≥⌊√𝑛⌋m≥⌊n⌋．设预处理的时间复杂度为 𝑇 0(𝑚)T 0(m)，则此时的 𝑇(𝑛)T(n) 为：

𝑇(𝑛)=𝑇 0(𝑚)+∑𝑘∈𝑅(𝑛);𝑘>𝑚 𝑇(𝑘)=𝑇 0(𝑚)+⌊𝑛/𝑚⌋∑𝑘=1 𝑂(√𝑛 𝑘)=𝑂(𝑇 0(𝑚)+∫𝑛/𝑚 0√𝑛 𝑥 d 𝑥)=𝑂(𝑇 0(𝑚)+𝑛√𝑚).T(n)=T 0(m)+∑k∈R(n);k>m T(k)=T 0(m)+∑k=1⌊n/m⌋O(n k)=O(T 0(m)+∫0 n/m n x d x)=O(T 0(m)+n m).

若 𝑇 0(𝑚) =𝑂(𝑚)T 0(m)=O(m)（如线性筛），由均值不等式可知：当 𝑚 =Θ(𝑛 2/3)m=Θ(n 2/3) 时，𝑇(𝑛)T(n) 取得最小值 𝑂(𝑛 2/3)O(n 2/3).

伪证一例

设计算 𝑆(𝑛)S(n) 的复杂度为 𝑇(𝑛)T(n), 则有：

𝑇(𝑛)=Θ(√𝑛)+𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2 𝑇(⌊𝑛 𝑖⌋)⎞⎟ ⎟⎠T(n)=Θ(n)+O(∑i=2⌊n⌋T(⌊n i⌋)) 𝑇(⌊𝑛 𝑖⌋)=Θ(√𝑛 𝑖)+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛/𝑖⌋∑𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋)⎞⎟ ⎟ ⎟⎠=𝑂(√𝑛 𝑖)T(⌊n i⌋)=Θ(n i)+O(∑j=2⌊n/i⌋T(⌊n i j⌋))=O(n i)

其中，𝑂(∑⌊√𝑛/𝑖⌋𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋))O(∑j=2⌊n/i⌋T(⌊n i j⌋)) 视作高阶无穷小，从而可以舍去．故：

𝑇(𝑛)=Θ(√𝑛)+𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2√𝑛 𝑖⎞⎟ ⎟⎠=𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=1√𝑛 𝑖⎞⎟ ⎟⎠=𝑂(∫√𝑛 0√𝑛 𝑥 d 𝑥)=𝑂(𝑛 3/4)T(n)=Θ(n)+O(∑i=2⌊n⌋n i)=O(∑i=1⌊n⌋n i)=O(∫0 n n x d x)=O(n 3/4)
Bug

问题在于「视作高阶无穷小，从而可以舍去」这一处．我们将 𝑇(⌊𝑛 𝑖⌋)T(⌊n i⌋) 代入 𝑇(𝑛)T(n) 的式子里，有：

𝑇(𝑛)=Θ(√𝑛)+𝑂⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2√𝑛 𝑖⎞⎟ ⎟⎠+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋)⎞⎟ ⎟ ⎟⎠=𝑂(√𝑛+∫√𝑛 0√𝑛 𝑥 d 𝑥)+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋)⎞⎟ ⎟ ⎟⎠=𝑂(𝑛 3/4)+𝑂⎛⎜ ⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋)⎞⎟ ⎟ ⎟⎠T(n)=Θ(n)+O(∑i=2⌊n⌋n i)+O(∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊n i j⌋))=O(n+∫0 n n x d x)+O(∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊n i j⌋))=O(n 3/4)+O(∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊n i j⌋))

我们考虑 ⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋)∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊n i j⌋) 这部分，不难发现：

⌊√𝑛⌋∑𝑖=2⌊√𝑛/𝑖⌋∑𝑗=2 𝑇(⌊𝑛 𝑖 𝑗⌋)=Ω⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2 𝑇⎛⎜ ⎜⎝⎢ ⎢ ⎢⎣𝑛 𝑖⋅⌊√𝑛 𝑖⌋−1⎥ ⎥ ⎥⎦⎞⎟ ⎟⎠⎞⎟ ⎟⎠=Ω⎛⎜ ⎜⎝⌊√𝑛⌋∑𝑖=2 𝑇(⌊√𝑛 𝑖⌋)⎞⎟ ⎟⎠∑i=2⌊n⌋∑j=2⌊n/i⌋T(⌊n i j⌋)=Ω(∑i=2⌊n⌋T(⌊n i⋅⌊n i⌋−1⌋))=Ω(∑i=2⌊n⌋T(⌊n i⌋))

由于没有引入记忆化，因此上式中的 𝑇(⌊√𝑛 𝑖⌋)T(⌊n i⌋) 仍然是 Ω((𝑛 𝑖)1/4)Ω((n i)1/4) 的，进而所谓的「高阶无穷小」部分是不可以舍去的．

实际上杜教筛的亚线性时间复杂度是由记忆化保证的．只有使用了记忆化之后才能保证不会出现那个多重求和的项．

## 例题

### 问题一

P4213【模板】杜教筛（Sum）

求 𝑆 1(𝑛) =∑𝑛 𝑖=1 𝜇(𝑖)S 1(n)=∑i=1 n μ(i) 和 𝑆 2(𝑛) =∑𝑛 𝑖=1 𝜑(𝑖)S 2(n)=∑i=1 n φ(i) 的值，1 ≤𝑛 <2 3 1 1≤n<2 31.

莫比乌斯函数前缀和欧拉函数前缀和

我们知道：

𝜖=[𝑛=1]=𝜇∗1=∑𝑑∣𝑛 𝜇(𝑑)ϵ=[n=1]=μ∗1=∑d∣n μ(d)𝑆 1(𝑛)=𝑛∑𝑖=1 𝜖(𝑖)−𝑛∑𝑖=2 𝑆 1(⌊𝑛 𝑖⌋)=1−𝑛∑𝑖=2 𝑆 1(⌊𝑛 𝑖⌋)S 1(n)=∑i=1 n ϵ(i)−∑i=2 n S 1(⌊n i⌋)=1−∑i=2 n S 1(⌊n i⌋)

时间复杂度的推导见 时间复杂度 一节．

对于较大的值，需要用 `map`/`unordered_map` 存下其对应的值，方便以后使用时直接使用之前计算的结果．

当然也可以用杜教筛求出 𝜑(𝑥)φ(x) 的前缀和，但是更好的方法是应用莫比乌斯反演．

莫比乌斯反演杜教筛
𝑛∑𝑖=1 𝑛∑𝑗=1[g c d(𝑖,𝑗)=1]=𝑛∑𝑖=1 𝑛∑𝑗=1∑𝑑∣𝑖,𝑑∣𝑗 𝜇(𝑑)=𝑛∑𝑑=1 𝜇(𝑑)⌊𝑛 𝑑⌋2∑i=1 n∑j=1 n[gcd(i,j)=1]=∑i=1 n∑j=1 n∑d∣i,d∣j μ(d)=∑d=1 n μ(d)⌊n d⌋2

由于题目所求的是 ∑𝑛 𝑖=1∑𝑖 𝑗=1[g c d(𝑖,𝑗) =1]∑i=1 n∑j=1 i[gcd(i,j)=1], 所以我们排除掉 𝑖 =1,𝑗 =1 i=1,j=1 的情况，并将结果除以 2 2 即可．

观察到，只需求出莫比乌斯函数的前缀和，就可以快速计算出欧拉函数的前缀和了．时间复杂度 𝑂(𝑛 2 3)O(n 2 3).

求 𝑆(𝑛) =∑𝑛 𝑖=1 𝜑(𝑖)S(n)=∑i=1 n φ(i).

同样的，𝜑 ∗1 =i d φ∗1=id, 从而：

𝑆(𝑛)=𝑛∑𝑖=1 𝑖−𝑛∑𝑖=2 𝑆(⌊𝑛 𝑖⌋)=1 2 𝑛(𝑛+1)−𝑛∑𝑖=2 𝑆(⌊𝑛 𝑖⌋)S(n)=∑i=1 n i−∑i=2 n S(⌊n i⌋)=1 2 n(n+1)−∑i=2 n S(⌊n i⌋)
代码实现

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    #include <cstring>
    #include <iostream>
    #include <map>
    using namespace std;
    constexpr int MAXN = 2000010;
    long long T, n, pri[MAXN], cur, mu[MAXN], sum_mu[MAXN];
    bool vis[MAXN];
    map<long long, long long> mp_mu;

    long long S_mu(long long x) {  // 求mu的前缀和
      if (x < MAXN) return sum_mu[x];
      if (mp_mu[x]) return mp_mu[x];  // 如果map中已有该大小的mu值，则可直接返回
      long long ret = (long long)1;
      for (long long i = 2, j; i <= x; i = j + 1) {
        j = x / (x / i);
        ret -= S_mu(x / i) * (j - i + 1);
      }
      return mp_mu[x] = ret;  // 路径压缩，方便下次计算
    }

    long long S_phi(long long x) {  // 求phi的前缀和
      long long ret = (long long)0;
      long long j;
      for (long long i = 1; i <= x; i = j + 1) {
        j = x / (x / i);
        ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
      }
      return (ret - 1) / 2 + 1;
    }

    int main() {
      cin.tie(nullptr)->sync_with_stdio(false);
      cin >> T;
      mu[1] = 1;
      for (int i = 2; i < MAXN; i++) {  // 线性筛预处理mu数组
        if (!vis[i]) {
          pri[++cur] = i;
          mu[i] = -1;
        }
        for (int j = 1; j <= cur && i * pri[j] < MAXN; j++) {
          vis[i * pri[j]] = true;
          if (i % pri[j])
            mu[i * pri[j]] = -mu[i];
          else {
            mu[i * pri[j]] = 0;
            break;
          }
        }
      }
      for (int i = 1; i < MAXN; i++)
        sum_mu[i] = sum_mu[i - 1] + mu[i];  // 求mu数组前缀和
      while (T--) {
        cin >> n;
        cout << S_phi(n) << ' ' << S_mu(n) << '\n';
      }
      return 0;
    }

### 问题二

「LuoguP3768」简单的数学题

大意：求

𝑛∑𝑖=1 𝑛∑𝑗=1 𝑖⋅𝑗⋅g c d(𝑖,𝑗)(m o d 𝑝)∑i=1 n∑j=1 n i⋅j⋅gcd(i,j)(mod p)

其中 𝑛 ≤1 0 1 0,5 ×1 0 8 ≤𝑝 ≤1.1 ×1 0 9 n≤10 10,5×10 8≤p≤1.1×10 9,𝑝 p 是质数．

利用 𝜑 ∗1 =i d φ∗1=id 做莫比乌斯反演化为：

𝑛∑𝑑=1 𝐹 2(⌊𝑛 𝑑⌋)⋅𝑑 2 𝜑(𝑑)∑d=1 n F 2(⌊n d⌋)⋅d 2 φ(d)

其中 𝐹(𝑛) =1 2 𝑛(𝑛 +1)F(n)=1 2 n(n+1)

对 ∑𝑛 𝑑=1 𝐹(⌊𝑛 𝑑⌋)2∑d=1 n F(⌊n d⌋)2 做数论分块，𝑑 2 𝜑(𝑑)d 2 φ(d) 的前缀和用杜教筛处理：

𝑓(𝑛)=𝑛 2 𝜑(𝑛)=(i d 2⁡𝜑)(𝑛)f(n)=n 2 φ(n)=(id 2⁡φ)(n)𝑆(𝑛)=𝑛∑𝑖=1 𝑓(𝑖)=𝑛∑𝑖=1(i d 2⁡𝜑)(𝑖)S(n)=∑i=1 n f(i)=∑i=1 n(id 2⁡φ)(i)

需要构造积性函数 𝑔 g，使得 𝑓 ×𝑔 f×g 和 𝑔 g 能快速求和．

单纯的 𝜑 φ 的前缀和可以用 𝜑 ∗1 φ∗1 的杜教筛处理，但是这里的 𝑓 f 多了一个 i d 2 id 2，那么我们就卷一个 i d 2 id 2 上去，让它变成常数：

𝑆(𝑛)=𝑛∑𝑖=1((i d 2⁡𝜑)∗i d 2)(𝑖)−𝑛∑𝑖=2 i d 2⁡(𝑖)𝑆(⌊𝑛 𝑖⌋)S(n)=∑i=1 n((id 2⁡φ)∗id 2)(i)−∑i=2 n id 2⁡(i)S(⌊n i⌋)

化一下卷积：

((i d 2⁡𝜑)∗i d 2)(𝑖)=∑𝑑∣𝑖(i d 2⁡𝜑)(𝑑)i d 2⁡(𝑖 𝑑)=∑𝑑∣𝑖 𝑑 2 𝜑(𝑑)(𝑖 𝑑)2=∑𝑑∣𝑖 𝑖 2 𝜑(𝑑)=𝑖 2∑𝑑∣𝑖 𝜑(𝑑)=𝑖 2(𝜑∗1)(𝑖)=𝑖 3((id 2⁡φ)∗id 2)(i)=∑d∣i(id 2⁡φ)(d)id 2⁡(i d)=∑d∣i d 2 φ(d)(i d)2=∑d∣i i 2 φ(d)=i 2∑d∣i φ(d)=i 2(φ∗1)(i)=i 3

再化一下 𝑆(𝑛)S(n):

𝑆(𝑛)=𝑛∑𝑖=1((i d 2⁡𝜑)∗i d 2)(𝑖)−𝑛∑𝑖=2 i d 2⁡(𝑖)𝑆(⌊𝑛 𝑖⌋)=𝑛∑𝑖=1 𝑖 3−𝑛∑𝑖=2 𝑖 2 𝑆(⌊𝑛 𝑖⌋)=(1 2 𝑛(𝑛+1))2−𝑛∑𝑖=2 𝑖 2 𝑆(⌊𝑛 𝑖⌋)S(n)=∑i=1 n((id 2⁡φ)∗id 2)(i)−∑i=2 n id 2⁡(i)S(⌊n i⌋)=∑i=1 n i 3−∑i=2 n i 2 S(⌊n i⌋)=(1 2 n(n+1))2−∑i=2 n i 2 S(⌊n i⌋)

分块求解即可．

代码实现

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    35
    36
    37
    38
    39
    40
    41
    42
    43
    44
    45
    46
    47
    48
    49
    50
    51
    52
    53
    54
    55
    56
    57
    58
    59
    60
    61
    62
    63
    64
    65
    66
    67
    68
    69
    70
    71
    72
    73
    74
    75
    // 不要为了省什么内存把数组开小,会卡80
    #include <cmath>
    #include <iostream>
    #include <map>
    using namespace std;
    constexpr int N = 5e6, NP = 5e6, SZ = N;
    long long n, P, inv2, inv6, s[N];
    int phi[N], p[NP], cnt, pn;
    bool bp[N];
    map<long long, long long> s_map;

    long long ksm(long long a, long long m) {  // 求逆元用
      long long res = 1;
      while (m) {
        if (m & 1) res = res * a % P;
        a = a * a % P, m >>= 1;
      }
      return res;
    }

    void prime_work(int k) {  // 线性筛phi，s
      bp[0] = bp[1] = true, phi[1] = 1;
      for (int i = 2; i <= k; i++) {
        if (!bp[i]) p[++cnt] = i, phi[i] = i - 1;
        for (int j = 1; j <= cnt && i * p[j] <= k; j++) {
          bp[i * p[j]] = true;
          if (i % p[j] == 0) {
            phi[i * p[j]] = phi[i] * p[j];
            break;
          } else
            phi[i * p[j]] = phi[i] * phi[p[j]];
        }
      }
      for (int i = 1; i <= k; i++)
        s[i] = (1ll * i * i % P * phi[i] % P + s[i - 1]) % P;
    }

    long long s3(long long k) {  // 立方和
      return k %= P, (k * (k + 1) / 2) % P * ((k * (k + 1) / 2) % P) % P;
    }

    long long s2(long long k) {  // 平方和
      return k %= P, k * (k + 1) % P * (k * 2 + 1) % P * inv6 % P;
    }

    long long calc(long long k) {  // 计算S(k)
      if (k <= pn) return s[k];
      if (s_map[k]) return s_map[k];  // 对于超过pn的用map离散存储
      long long res = s3(k), pre = 1, cur;
      for (long long i = 2, j; i <= k; i = j + 1)
        j = k / (k / i), cur = s2(j),
        res = (res - calc(k / i) * (cur - pre) % P) % P, pre = cur;
      return s_map[k] = (res + P) % P;
    }

    long long solve() {
      long long res = 0, pre = 0, cur;
      for (long long i = 1, j; i <= n; i = j + 1) {
        j = n / (n / i);
        cur = calc(j);
        res = (res + (s3(n / i) * (cur - pre)) % P) % P;
        pre = cur;
      }
      return (res + P) % P;
    }

    int main() {
      cin.tie(nullptr)->sync_with_stdio(false);
      cin >> P >> n;
      inv2 = ksm(2, P - 2), inv6 = ksm(6, P - 2);
      pn = (long long)pow(n, 0.666667);  // n^(2/3)
      prime_work(pn);
      cout << solve();
      return 0;
    }

### 参考资料

1. 任之洲，2016，《积性函数求和的几种方法》，2016 年信息学奥林匹克中国国家队候选队员论文
2. 杜教筛的时空复杂度分析 - riteme.site

> 本页面最近更新：2026/1/7 08:56:54，更新历史
> 发现错误？想一起完善？ 在 GitHub 上编辑此页！
> 本页面贡献者：StudyingFather, Tiphereth-A, hsfzLZH1, Ir1d, Enter-tainer, sshwy, c-forrest, Marcythm, MegaOwIer, Henry-ZHR, Xeonacid, Backl1ght, Great-designer, huayucaiji, kenlig, ksyx, Menci, Nanarikom, nanmenyangde, ouuan, purple-vine, shawlleyw, Sshwy
> 本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用