前言

本文总结了一些与 OI 有关的简单数论内容。 由于作者水平有限，部分结论的证明过程省略。

一. 公约数/扩展欧几里得算法

这里将正整数 a,b 的最大公约数记作 gcd(a,b) . 根据欧几里得算法，有 gcd(a,b)=gcd(b,a\mod b) . 令 d=gcd(a,b) ，由裴蜀定理可知，不定方程 ax+by=d 有正整数解。设 bx_0+(a\mod b)y_0=d 拆开取模，得 bx_0+(a-\lfloor \frac{a}{b}\rfloor b)y_0=d 整理得 ay_0+b(x_0-\lfloor \frac{a}{b}\rfloor y_0)=d 即 \left\{\begin{matrix}x_1=y_0\\y_1=x_0-\lfloor \frac{a}{b}\rfloor y_0\end{matrix}\right. 迭代求解即可，复杂度 O(\log n) .

二. 中国剩余定理

有物不知其数，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

上面的问题出自《孙子算经》。用现代数学的角度看这个问题，可以归纳为求解以下形式的同余方程组： \left\{\begin{matrix}x\equiv a_1(\mod m_1)\\x\equiv a_2(\mod m_2)\\x\equiv a_3(\mod m_3)\\\cdots\\x\equiv a_k(\mod m_k)\end{matrix}\right. 其中 m_1,m_2,m_3,\cdots,m_k 为两两互质的正整数。求 x 的最小正整数解和方程的通解。

对于这个方程组，我们尝试构造这样一组解： x=\sum_{i=1}^na_it_i 而 t_i 的选取，只需满足以下关系：

• t_i\equiv 1(\mod m_i)
• t_i\equiv 0(\mod m_j),j\neq i 令 M=\prod_{i=1}^k m_i ，即 m_1\cdots m_k 的最小公倍数。那么 t_i 取方程 \frac M{M_i} t_i\equiv 1(\mod m_i) 的最小非负整数解即可满足上述关系。t_i 可以用扩展欧几里得算法求出。将其代入前面的式子，可以得到方程的一组解 x ，而通解就是 x+iM(i\in\mathbb{Z}). 特别的，方程的最小整数解为 (x\mod M+M))\mod M.

扩展中国剩余定理

此定理适用于 m_1,m_2,m_3,\cdots,m_k 不一定两两互质的情况。其实它与普通 CRT 使用的方法并不一样，但似乎更容易理解一些。 假设我们已经求出了前 p 个同余方程组成的方程组的通解为 x+iM ，现在要求加入第 p+1 个同余方程组的通解 x+iM'. 显然前 p+1 个方程的一个解也是前 p 个方程的一个解，因此 M' 一定是 M 的倍数，记为 M'=tM. 那么 t 应为以下同余方程的最小正整数解： x+tM\equiv a_{p+1}(\mod m_{p+1}) 若此方程无解，则原方程组无解；否则，我们可以用扩展欧几里得算法求出 t ，从而求出 M'. 如此反复进行 k 次 ex-GCD ，最终可以将 k 个同余方程合并为一个。这样我们就完成了求解。

三. Lucas 定理

Lucas 定理主要用于求组合数对给定质数取模的结果，即 C_n^m\mod p ，其中 n,m,p 同阶，p 为质数。这时不能预处理阶乘再求逆了，因为当 n\geq p 时组合数的分母上会含有 p 这个因子，所以它模 p 余 0. Lucas 定理的内容为以下公式： C_n^m=C_{\lfloor \frac n p\rfloor}^{\lfloor \frac m p\rfloor}\cdot C_{n\mod p}^{m\mod p} 其证明略为复杂，再此省略。同样迭代求解直到 n<p 即可。

扩展 Lucas

此定理适用于 p 不为质数的情况。 先假设 p 可以表示成一个质数的幂次的形式，即 p=p_0^{\alpha} . 则答案为 \frac{\frac{n!}{p_0^x}}{\frac{m!}{p_0^y}\frac{(n-m)!}{p_0^z}}p_0^{x-y-z}\mod p_0^\alpha 其中 x,y,z 分别表示 n!,m!,(n-m)! 中质因子 p_0 的次数。 求一个数中某个质因子的次数可以这么做：设递推函数 g(n,p)=g(\lfloor \frac n {p}\rfloor,p)+\lfloor \frac n {p}\rfloor ，则 g(n,p) 即为 n 中质因子 p 的次数。 若 p 不能表示成质数的幂次，则将 p 质因数分解，得 p=p_1^{\alpha_1}p_2^{\alpha_2}\cdots p_k^{\alpha_k} 对于每个质因子，求出 \left\{\begin{matrix}C_n^m\mod p_1^{\alpha 1}\\C_n^m\mod p_2^{\alpha 2}\\\cdots\\ C_n^m\mod p_k^{\alpha k}\end{matrix}\right. 并用中国剩余定理将其合并即可。

四. 原根

定义

设 m 是正整数，a 是整数，若 a 模 m 的阶等于 \varphi(m)，则 a 为模 m 的一个原根。 ——百度百科

数论阶

若 a,m 为整数且互质，满足 a^r\equiv 1(\mod m) 的最小整数 r，称为 a 模 m 的阶，通常记作 \delta_m(a)=r . 由欧拉定理，有 a^{\varphi(m)}\equiv 1(\mod m) ，因此对于任意与 m 互质的整数 a ，都有 \delta_m(a)\leq \varphi(m).

由以上定义可知，不等式等号成立时，a 为模 m 的一个原根，a^1,a^2,a^3,...,a^{\varphi(m)} 在模 m 意义下都不为 1 且构成了 m 的一个缩系。当 m 为质数时，它们取遍了 [1,m-1] 内每个整数恰好一次。这使原根具有了一些优秀的性质和应用，比如将模意义下的乘法转化为加法。此外，要判定一个数是否有原根，可以使用以下方法：

• 模 m 有原根的充要条件是 m=1,2,4,p,2p,p^n ，其中 p 是奇素数，n 是任意正整数
• 若模一个数 m 有原根，则他有 \varphi(\varphi(m)) 个原根

求 m 的任意一个原根

一般来说，一个数最小的原根会比较小，大约是 O(m^{\frac 1 4}) 的。考虑从 2 开始枚举，对于一个数 a ，依次判断 a^1...a^{\varphi(m)} 在模 m 意义下是否为 1 . 这样做的效率比较慢，因为我们做了很多无用的判断。若 a^n\not\equiv 1(\mod m) ，则对于 n 的任意一个因数 d 同样有 a^d\not\equiv 1(\mod m) . 因此，我们只需枚举 \varphi(m) 的每个质因子 p_i ，判断 a^{\frac{\varphi(m)}{p_i}} 模 m 是否为 1 就足够了。这样做的复杂度是 O(\sqrt m+m^{\frac 1 4}\log^2 m) 的。

代码

int get_pr(int n)
{
    static int f[MAXN];
    int t=n-1, k=0;
    for (int i=2; i*i<n; i++)
        if (t%i==0)
        {
            f[++k]=i;
            while (t%i==0) t/=i;
        }
    for (int i=2; i<n; i++)
    {
        bool flag=1;
        for (int j=1; j<=k; j++)
            if (qpow(i, (n-1)/f[j])==1)
            {
                flag=0;
                break;
            }
        if (flag) return i;
    }
}

五. BSGS 算法

求以下同余方程的正整数解： a^x\equiv b(\mod p) 其中 a 与 p 互质。 令 x=qm-r . 有 a^{qm}=a^rb 在 [0,m-1] 范围内枚举 r 并将其值用哈希或 map 存下，在 [0,\lfloor \frac p m\rfloor] 范围内枚举 q 并查看等号右边是否有与其相等的值。容易证明当 m=\sqrt p 时复杂度最优，为 O(\sqrt p).

扩展 BSGS

适用于 a 不与 p 互质的情况。 考虑将这个问题转化为 a 与 p 互质的情况。令 d=gcd(a,p) 将方程写为以下形式： a^x+kp=b 由裴蜀定理可知方程有解的条件为 d|b. 因此将等式两边同时除 d ： \frac a d a^{x-1}+k\frac p d=\frac z d 一直这样处理，直到 gcd(a,p')=1 为止。设递归了 c 此，则方程可化为： \frac a {d^c}\cdot y^{x-c}+k\cdot \frac p {d^c}=\frac b {d^c} 即 \frac a {d^c}\cdot y^{x-c}=\frac b {d^c}(\mod \frac p {d^c}) 使用 BSGS 算法求解后回代即可。复杂度仍为 O(\sqrt p).

六. 整除分块

与其说整除分块是一种算法，不如说它是一个技巧。有时我们需要计算这样的式子： \sum_{i=1}^n f(i)\lfloor\frac{n}{i}\rfloor 根据经验，我们发现 \lfloor\frac{n}{i}\rfloor 只有 O(\sqrt n) 种取值。对于每种取值，i 都会有一个连续的范围。假设函数 f 的前缀和可以预处理后快速求出，那么我们可以枚举 \lfloor\frac{n}{i}\rfloor 的所有取值，并和 f 的连续一段的和相乘。具体可以见代码：

// s[i] 为函数 f(i) 的前缀和 
int ans=0;
for (int i=1, j; i<=n; i=j+1)
{
    j=n/(n/i);
    ans+=(s[j]-s[i-1])*(n/i);
}
printf("%d\n", ans);

七. 莫比乌斯反演

莫比乌斯函数

定义莫比乌斯函数 \mu(n) 如下： 首先将 n 质因数分解，表示成 n=\prod_{i=1}^{k}{p_i}^{\alpha_i} 则 \mu(n)=\left\{\begin{matrix}1&n=1\\(-1)^k&\forall i,\alpha_i=1\\0&\exists i,\alpha_i>1\end{matrix}\right. 用语言描述，就是当 n=1 时， \mu(n)=1 ；当 n 没有幂次大于平方的质因子（即 n 的所有质因子都不同）时，若 n 有奇数个质因子，则 mu(n)=-1 ，若 n 有偶数个质因子，则 mu(n)=1 ；否则 n 有幂次大于 2 的质因子，此时函数值 \mu(n)=0 。综上所述，莫比乌斯函数只有 1, -1, 0 这三种取值。

性质

• 性质一：对于任意一对互质的正整数 p,q ， \mu(pq)=\mu(p)\mu(q) 。所以 莫比乌斯函数是积性函数 。这一性质说明莫比乌斯函数可以被线性筛。

证明： 根据定义，可知若 p,q 中一个包含幂次大于平方的质因子，则它的莫比乌斯函数值为 0 ，那么这两个函数值的乘积也为 0 。否则，因为 p,q 互质，所以它们没有相同的质因子。设 p,q 的质因子数分别为 f(p),f(q) ，则 f(pq)=f(p)+f(q) 因为莫比乌斯函数的取值由质因子数的奇偶性决定，容易推出 \mu(pq)=\mu(p)\mu(q) 注意，对于不互质的正整数对，由于存在相同的质因子，此命题不成立。所以莫比乌斯函数是积性函数，但不是完全积性函数。

• 性质二：对于任意正整数 n ， \sum_{d|n}\mu(d)=[n=1] 。（这里的 [n=1] 表示当且仅当括号内条件成立时值为 1 ，否则为 0 。下文同。）

证明： 根据定义，当 n=1 时命题显然成立。 当 n>1 时， 将 n 按上文形式分解质因数。当 a_i>1 时函数值为 0 ，对整个式子的值没有贡献，可以不用考虑。我们只统计有多少个质因子的指数为 1 。那么上述式子可以化为： \sum_{i=0}^k C_k^i*(-1)^i 由二项式定理得： \sum_{i=0}^k\begin{pmatrix}k\\i\\\end{pmatrix}(-1)^i1^{k-i}=[(-1)+1]^k=0

• 性质三：对于任意正整数 n ， \sum_{d|n}\frac{\mu(d)}d=\frac{\varphi(n)}n

这条性质用到的地方不多，因此证明暂时省略。

狄利克雷卷积

狄利克雷卷积是一种数论函数卷积。若数论函数 f,g,h 满足 h(n)=\sum_{d|n}f(d)g(\frac n d) 则称函数 f 为函数 g 和函数 h 的狄利克雷卷积，记作 h(n)=(f*g)(n) ，下文简单地记作 h=f*g 。可以证明狄利克雷卷积满足交换律和结合律。此外，两个积性函数的狄利克雷卷积仍是一个积性函数。

一些基本的完全积性函数

为了下文展开方便，下面给出一些完全积性函数的记号，这些函数都满足对于任意正整数 a,b ，f(ab)=f(a)f(b) .

• 单位元函数 \varepsilon(n)=[n=1]
• 常函数 I(n)=1
• 恒等函数 id(n)=n

容易看出，上文中莫比乌斯函数的性质二可记作 \mu*I=\varepsilon ，性质三可记作 \mu*id=\varphi.

莫比乌斯反演定理

f(n) 和 g(n) 为数论函数，若它们满足 f(n)=\sum_{d|n}g(d) 则有 g(n)=\sum_{d|n}\mu(d)f(\frac n d)

证明： 原式可以表示为狄利克雷卷积的形式，即 f=g*I. 在等式两边同时卷上 \mu ，得 f*\mu=g*I*\mu 由 \mu*I=\varepsilon ，得 f*\mu=g 展开即结果式。

例题

Luogu3455&BZOJ1101 【POI2007 ZAP-Queries】

题意

T 组数据，每组数据给定 n,m ，求 \sum_{i=1}^n\sum_{j=1}^m[gcd(i,j)=k] T\leq 50000,n,m\leq 50000.

题解

不妨设 n\leq m. 原始可以化为 =\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}[gcd(i,j)=1] 利用莫比乌斯函数性质二，将 [gcd(i,j)=1] 换成约数 \mu 和的形式，得 =\sum_{i=1}^{\lfloor\frac{n}{k}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{k}\rfloor}\sum_{d|gcd(i,j)}\mu(d) 枚举 d ，那么根据上式，i,j 应分别在 \lfloor\frac{n}{k}\rfloor,\lfloor\frac{m}{k}\rfloor 的范围内且都被 d 整除，这样的 (i,j) 对数是可以快速计算的。于是式子化为： =\sum_{d=1}^n\mu(d)\lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor 观察到后面的 \lfloor\frac{n}{kd}\rfloor\lfloor\frac{m}{kd}\rfloor 是可以整除分块的，因此我们可以线性筛出 \mu ，这样就可以做到单组数据 O(\sqrt n) 的时间复杂度。

八. 杜教筛

杜教筛是一种求积性函数前缀和的低于线性的筛法。我们有时要求 S(n)=\sum_{i=1}^n f(i) ，其中 f 为一个积性函数。但 n 可能有 10^10 之大，线性筛的复杂度不能承受。此时我们就可以使用这种算法。 首先构造两个积性函数 g,h ，使 h=f*g ，且 g,h 的前缀和可以快速计算. 展开得 h(i)=\sum_{d|i}f(\frac i d)g(d) 对式子两边求和，得 \sum_{i=1}^n h(i)=\sum_{i=1}^n\sum_{d|i}f(\frac i d)g(d) 将 d 提到前面，式子化为 \sum_{i=1}^n h(i)=\sum_{d=1}^n g(d)\sum_{i=1}^{\lfloor\frac n d\rfloor}f(i) =\sum_{d=1}^n g(d)S(\lfloor\frac n d\rfloor) 将 g(1) 那一项单独提出来，得 \sum_{i=1}^n h(i)=g(1)S(n)+\sum_{d=2}^n g(d)S(\lfloor\frac n d\rfloor) 上式就是杜教筛的套路式。移项，得 S(n)=\frac{\sum_{i=1}^n h(i)-\sum_{i=2}^n g(i)S(\lfloor\frac n i\rfloor)}{g(1)} 由于 g 是积性函数，g(1)=1 ，所以 S(n)=\sum_{i=1}^n h(i)-\sum_{i=2}^n g(i)S(\lfloor\frac n i\rfloor) 发现 S(\lfloor\frac n i\rfloor) 那部分可以整除分块，于是递归去求 S(i) ，这样做的复杂度是 O(n^{\frac 3 4}) 的。但是可以在 i\leq n^{\frac 2 3} 的时候直接线性筛预处理，并开个哈希/映射以下存下 i>n^{\frac 2 3} 的答案，这样复杂度可以优化到 O(n^{\frac 2 3}).

例题

Luogu3768 【简单的数学题】

题意

给定 n,p ，求 (\sum_{i=1}^n\sum_{j=1}^n ijgcd(i,j))\mod p . n\leq 10^{10},5\times 10^8\leq p \leq 1.1\times 10^9

题解

先根据莫比乌斯反演的老套路，设 gcd(i,j)=k ，把 k 提出来： Ans=\sum_{k=1}^n k\sum_{i=1}^n\sum_{j=1}^n ij[gcd(i,j)=k] =\sum_{k=1}^n k^3\sum_{i=1}^{\lfloor\frac n k\rfloor}\sum_{j=1}^{\lfloor\frac n k\rfloor} ij[gcd(i,j)=1] 对 [gcd(i,j)=k] 莫比乌斯反演： =\sum_{k=1}^n k^3\sum_{i=1}^{\lfloor\frac n k\rfloor}\sum_{j=1}^{\lfloor\frac n k\rfloor} ij\sum_{d\mid i,j}\mu(d) 把 d 提出来，i,j 都是 d 的倍数： =\sum_{k=1}^n k^3\sum_{d=1}^{\lfloor\frac n k\rfloor}\mu(d)d^2\sum_{i=1}^{\lfloor\frac n{kd}\rfloor}\sum_{j=1}^{\lfloor\frac n{kd}\rfloor} ij 设 S(n)=\sum_{i=1}^n i ，则有： =\sum_{k=1}^n k^3\sum_{d=1}^{\lfloor\frac n k\rfloor}\mu(d)d^2 S^2(\lfloor\frac n{kd}\rfloor) 令 T=kd ，先枚举 T ，再枚举 k ： =\sum_{T=1}^n\sum_{k\mid T}k^3\mu\left(\frac T k\right)\left(\frac T k\right)^2 S(\lfloor\frac n{T}\rfloor) 约分并把只与 T 有关的项提到前面： =\sum_{T=1}^nS(\lfloor\frac n{T}\rfloor)T^2\sum_{k\mid T}k\mu\left(\frac T k\right) 考虑莫比乌斯函数的第三条性质即 \mu*id=\varphi ，式子最后的 \sum_{k\mid T}k\mu\left(\frac T k\right) 可以化为 \varphi(T) ： =\sum_{T=1}^nS(\lfloor\frac n{T}\rfloor)T^2\varphi(T) 前面的 S(\lfloor\frac n{T}\rfloor) 可以整除分块，而后面的 T^2\varphi(T) 显然是一个积性函数。考虑怎么用杜教筛筛出它的前缀和。 设 f(n)=n^2\varphi(n) ，h=f*g 。展开得： h(n)=\sum_{d\mid n}d^2\varphi(d)g(\frac n d) 我们令 g(i)=i^2 ，就可以把式子右边的 d^2 消掉： h(n)=n^2\sum_{d\mid n}\varphi(d) 而欧拉函数 \varphi 有一个众所周知的性质：\sum_{d\mid n}\varphi(d)=n ，即 \varphi*I=id. 于是： h(n)=n^3 而 h(n) 的前缀和可以用以下公式计算： \sum_{i=1}^n i^3=\left[\frac{n(n+1)}2\right]^2 这样，我们就可以在低于线性的时间内计算出 f(n) 的前缀和了。

代码

#include <cstdio>
#include <cmath>
#include <cstring>
#include <tr1/unordered_map>
using namespace std;
using namespace std;
typedef long long ll;
const int MAXN=5E6;
int p, lim, inv2, inv6; ll n;
int prim[MAXN], phi[MAXN];
int sum[MAXN];
bool book[MAXN];
tr1::unordered_map<ll, int> hash;
void exgcd(int a, int b, int& x, int& y)
{
    if (b==0) x=1, y=0;
    else exgcd(b, a%b, y, x), y-=a/b*x;
}
inline int pre(ll n)
{
    n=n%p;
    return n*(n+1)%p*inv2%p;
}
inline int square(ll n)
{
    n=n%p;
    return n*(n+1)%p*(2*n+1)%p*inv6%p;
}
inline int cube(ll n)
{
    return 1ll*pre(n)*pre(n)%p;
}
void sieve()
{
    memset(book, 0, sizeof book);
    phi[1]=1;
    for (int i=2, cp=0; i<=lim; i++)
    {
        if (!book[i]) prim[++cp]=i, phi[i]=i-1;
        for (int j=1; j<=cp&&i*prim[j]<=lim; j++)
        {
            book[i*prim[j]]=1;
            if (i%prim[j]==0)
            {
                phi[i*prim[j]]=phi[i]*prim[j];
                break;
            }
            phi[i*prim[j]]=phi[i]*(prim[j]-1);
        }
    }
    for (int i=1; i<=lim; i++) sum[i]=(sum[i-1]+1ll*i*i%p*phi[i]%p)%p;
}
int djs(ll n)
{
    if (n<=lim) return sum[n];
    if (hash.count(n)) return hash[n];
    int res=cube(n);
    for (ll d=2, k; d<=n; d=k+1)
    {
        k=n/(n/d);
        res=(res-1ll*(square(k)-square(d-1))*djs(n/d)%p+p)%p;
    }
    return hash[n]=res;
}
int main()
{
//  freopen("P3768.in", "r", stdin);
//  freopen("P3768.out", "w", stdout);
    scanf("%d%lld", &p, &n);
    int t;
    exgcd(2, p, inv2, t);
    exgcd(6, p, inv6, t);
    inv2=(inv2+p)%p; inv6=(inv6+p)%p;
    lim=pow(n, 2.0/3);
    sieve();
    int ans=0;
    for (ll i=1, j; i<=n; i=j+1)
    {
        j=n/(n/i);
        int s=pre(n/i);
        ans=(ans+1ll*s*s%p*(djs(j)-djs(i-1))%p+p)%p;
    }
    printf("%d\n", ans);
    return 0;
}

九. Powerful Number 优化积性函数前缀和

这种求积性函数前缀和的方法和杜教筛十分相似。先以一道题为例：

已知 f(x) 为积性函数，满足 f(p^q)=p^k ，其中 p 为质数，k 为定值。求 f(x) 的前 n 项和。

先定义 Powerful Number 为每个质因子的幂次都 \geq 2 的正整数。则有结论：n 以内 Powerful Number 的数量是 O(\sqrt n) 的。

证明： 由于每个质因子间是独立的，易证所有 Powerful Number 都可以表示成 a^2b^3 的形式。考虑枚举 a 则 n 以内这种数的数量为 O\left(\sum_{a=1}^{\sqrt n}\left(\frac n{a^2}\right)^{\frac 1 3}\right) =O\left(\int_1^{\sqrt n}\left(\frac n{x^2}\right)^{\frac 1 3}dx\right) =O(\sqrt n) 考虑构造数论函数 g,h ，使 f=g*h ，且 g(p)=f(p). 可以令 g(x)=x^k ，这样当 p 为质数时，有 f(p)=g(1)h(p)+g(p)h(1) 由于 g(1)=h(1)=1 ，可以得出 h(p)=0 . 又因为 h 是积性函数，所以仅当 x 为 Powerful Number 时，h(x) 才不为 0. 考虑求解答案。答案为 \sum_{i=1}^n f(i)=\sum_{ij\leq n}g(i)h(j) =\sum_{i=1}^nh(i)\sum_{j=1}^{\lfloor\frac n i\rfloor}g(j) 后面的 g 的前缀和可以用插值等求自然数幂和的方法快速求出。现在唯一的问题就是求所有 Powerful Number x 的 h(x). 由于 f,g 都是积性函数，所以 h 也是积性函数。因此我们考虑求 h(p^q). 由 f=g*h ，有以下等式： f(p^q)=\sum_{i=0}^q g(p^i)h(p^{q-i}) p^k=\sum_{i=0}^q p^{ik}h(p^{q-i}) 令 h(p^q)=p^k-p^{2k} ，等式成立。这样我们就可以在搜索 Powerful Number 时顺便把它们的 h(x) 的值也求出来了。

十. min-25 筛

以下内容部分引用自租酥雨的博客，并加入了一些自己的理解。

min-25 筛与杜教筛相似，也是一种求积性函数前缀和的方法。一般来说，它的适用范围更加广泛，设 p 为质数，则它对要求的积性函数只有两条要求：

• f(p) 可以用简单多项式表示
• f(p^k) 可以快速计算

预处理

设 P 为质数集合，我们需要预处理出对于每个 x=\lfloor\frac n i\rfloor ，\sum_{i=1}^x [i\in P] f(i) 的值。 考虑模拟埃氏筛法的过程，先将所有 i 视为质数，代入 f(p) 的多项式，求出 \sum_{i=1}^n f(i) . 然后，枚举每个质数，筛掉所有以它为最小质因子的数，最终得到答案。但是这样的复杂度是 O(n\ln n) ，我们无法承受。考虑到 f 是积性函数的性质，我们可以用 DP 的方法解决这个问题。 设 g(n,j)=\sum_{i=1}^n[i\in P\;or\;lpf(i)>P_j]f(i) ，其中 lpf(i) 表示 i 的最小质因子。也就是说，g(n,j) 就是对于在 [1,n] 范围内对前 j 个质数做了埃氏筛法后没有筛到的数，将这些数视为质数代入 f(p) 的多项式求出的值之和。我们发现，要求的 \sum_{i=1}^x [i\in P] f(i) 就是 g(x,|P|) ，其中 |P| 是 \sqrt n 以内质数集的大小。

考虑 g(n,j) 的转移，分两种情况：

1. P_j^2>n。此时运行的第j次已经不会再筛掉任何数了（因为第 j 次运行中筛掉的最小的数是 P_j^2，所以此时 g(n,j)=g(n,j−1).
2. P_j^2\leq n。这时候我们就要考虑哪些数被筛掉了。被筛掉的数一定含有质因子 P_j ，且除掉 P_j 后最小的质因子会大于等于 P_j 。考虑减去 f(P_j)\cdot g(\lfloor\frac n{P_j}\rfloor,j−1)，但在 g(\lfloor\frac n{P_j}\rfloor,j−1) 中多减去了 \sum_{i=1}^{j-1}f(P_i) 这些最小质因子小于 P_j 的函数值，所以再把它们加上就好了。 所以总结起来就是： g(n,j)=\left\{\begin{matrix}g(n,j-1) & P_j^2>n\\g(n,j-1)-f(P_j)\left[g(\lfloor\frac n{P_j}\rfloor,j-1)-\sum_{i=1}^{j-1}f(P_i)\right]&P_j^2\leq n\end{matrix}\right.

观察 DP 方程，发现转移中需要用到的第一维状态就是 \lfloor\frac n i\rfloor 的所有取值，共 O(\sqrt n) 种。而最终我们需要的只有第二维为 |P| 的状态，其他状态都是过程量，并无实际意义。因此可以用滚动数组将空间优化到 O(\sqrt n).

复杂度证明： 考虑对于每个 x=\lfloor\frac n i\rfloor ，它只会在枚举 P_j\leq\sqrt n 时产生贡献。由素数定理，\sqrt n 内的素数有 O\left(\frac{\sqrt n}{\log\sqrt n}\right) 个。因此，时间复杂度可估计为 \sum_{i=1}^{\sqrt n}O\left(\frac{\sqrt i}{\log\sqrt i}\right)+\sum_{i=1}^{\sqrt n}O\left(\frac{\sqrt \frac n i}{\log\sqrt\frac n i}\right) =O\left(\int_1^{\sqrt n}\frac{\sqrt \frac n x \mathrm{d}x}{\log\sqrt\frac n x}\right) =O\left(\frac{n^{\frac 3 4}}{\log n}\right)

求解

我们设 S(n,j)=\sum_{i=1}^n[lpf(i)\geq P_j]f(i) ，也就是所有满足最小质因子大于等于 P_j 的 f 值之和。 那么最终的答案就是 S(n,1)+f(1) . 鉴于质数的答案我们已经算出来了，是 g(n,|P|)-\sum_{i=1}^{j-1}f(P_i)。（因为要保证最小质因子大于等于Pj所以要把小于它的质数减掉） 考虑合数。我们枚举这个合数的最小质因子及其出现次数，然后直接乘即可。 S(n,j)=g(n,|P|)-\sum_{i=1}^{j-1}f(P_i)+\sum_{k=j}^{p_k^2\leq n}\sum_{e=1}^{p_k^{e+1}\leq n}f(p_k^e)\cdot S(\left\lfloor\frac{n}{p_k^e}\right\rfloor,k+1)+f(p_k^{e+1})

解释一下这个式子，设有一合数为 x ，则 p_k 枚举的是 x 的最小质因子，而 e 枚举的是这个质因子在 x 中的幂次（即 p_k^e\mid x,p_k^{e+1}\nmid x ），则除去这个质因子后它的最小质因子应不小于 k+1. 但是这样会有表示为一个质数的幂次形式的合数没有被考虑到，因此应额外加上 f(p_k^{e+1}). 这一部分看似很暴力，实际上复杂度也被证明是 O\left(\frac{n^{\frac 3 4}}{\log n}\right) 的。具体证明过程请参考朱震霆的 2018 年国家集训队论文《一些特殊的数论函数求和问题》。

例题

LOJ6053 【简单的函数】

题意

定义积性函数 f(x) 满足 f(1)=1,f(p^c)=p\oplus c ，其中 \oplus 表示按位异或，求 f(x) 的前 n 项和模 10^9+7 的值。 n\leq 10^{10} .

题解

要应用 min-25 筛，首先应把 f(p) 表示为简单多项式的形式。观察到对于所有奇素数 p ，都有 f(p)=p\oplus 1=p-1 . 这个函数虽然不满足积性，但可以将它拆成两个积性函数 g(p)=p 和 h(p)=1 相减。 于是我们对于每个 x=\lfloor\frac n i\rfloor ，预处理出 \sum_{i=1}^x [i\in P] i 和 \sum_{i=1}^x [i\in P] 1 再将两者作差就可以了。 值得一提的是，本题也可以用上面介绍的 Powerful Number 的方法求解，效率更高，在此不再展开说明。

代码

#include <cstdio>
#include <cmath>
#include <cstring>
typedef long long ll;
const int MAXN=110000;
const int MAXM=2*MAXN;
const int MOD=1E9+7;
ll n; int m, sq;
int cp, pri[MAXN], sp[MAXM];
bool book[MAXN];
ll w[MAXM];
int id1[MAXM], id2[MAXM];
int g[MAXM], h[MAXM];
inline int id(ll i)
{
    return i<=sq?id1[i]:id2[n/i];
}
void get_prime(int lim)
{
    memset(book, 0, sizeof book);
    cp=0;
    for (int i=2; i<=lim; i++)
    {
        if (!book[i]) pri[++cp]=i, sp[cp]=(sp[cp-1]+i)%MOD;
        for (int j=1; j<=cp&&i*pri[j]<=lim; j++)
        {
            book[i*pri[j]]=1;
            if (i%pri[j]==0) break;
        }
    }
}
void dp()
{
    for (ll i=1, j; i<=n; i=j+1)
    {
        j=n/(n/i); w[++m]=n/i;
        if (w[m]<=sq) id1[w[m]]=m;
        else id2[n/w[m]]=m;
        g[m]=((w[m]+2)%MOD)*((w[m]-1)%MOD)%MOD;
        if (g[m]&1) g[m]+=MOD; g[m]/=2;
        h[m]=(w[m]-1)%MOD;
    }
    for (int j=1; j<=cp; j++)
        for (int i=1; i<=m&&1ll*pri[j]*pri[j]<=w[i]; i++)
        {
            g[i]=((g[i]-1ll*pri[j]*(g[id(w[i]/pri[j])]-sp[j-1]))%MOD+MOD)%MOD;
            h[i]=(h[i]-(h[id(w[i]/pri[j])]-(j-1))+MOD)%MOD;
        }
}
int calc(ll i, int j)
{
    if (i<pri[j]) return 0;
    int res=((g[id(i)]-h[id(i)])-(sp[j-1]-(j-1))+MOD)%MOD;
    for (int k=j; k<=cp&&1ll*pri[k]*pri[k]<=i; k++)
    {
        ll p; int e;
        for (p=pri[k], e=1; p*pri[k]<=i; p*=pri[k], e++)
            res=(res+1ll*(pri[k]^e)*calc(i/p, k+1)+(pri[k]^e+1))%MOD;
    }
    return res;
}
int main()
{
//  freopen("loj6053.in", "r", stdin);
//  freopen("loj6053.out", "w", stdout);
    scanf("%lld", &n);
    sq=sqrt(n);
    get_prime(sq);
    dp();
    printf("%d\n", n==1?1:(calc(n, 1)+3)%MOD);
}

十一. min-25 筛的两种优化

1. min-25 本人在博客中提到的优化