道长的算法笔记：数论基础汇总

(一) 素数判定与筛选

　给定一个正整数 \(N\)，如果存在一个数 \(T\)，\(T\) 满足 \((2\leq T \leq N -1)\) 则称 \(N\) 是一个合数，如果不存在这样这样的因数 \(T\)，则称\(N\) 质数。简单来说，一个数\(N\) 如何仅能被 \(1\) 与 \(N\) 本身整除，则称这个数字是质数，或称素数(Prime Number)；数论的大多算法均以素数作为基础。

　接下来我们讨论精确的正整数的素数判定算法，以及给定正整数\(N\)，试求\([1,N]\)之内共有多少素数的算法，也称素数筛法。素数判定实际存在更快的、近似的算法，出错概率极小，密码学领域经常使用，但在算法题中用得很少。

(1.1) 素数判定

　根据素数定义，只要循环判对 \(2\leq i \leq N-1\) 范围之内的正整数是否存在某个数能够整除 \(N\) 即可，但当给定数值很大时，比如 \(N>2^{30}\)，这种算法效率较差。由于因数 \(d|N\) 则其必有 \(\frac{N}{d}|N\)，因而我们只需扫描 \(2\leq i \leq \sqrt{N}\) 这个范围即可，另外需要注意，若把判定条件写成 \(i<\text{sqrt(N)}\) 每次循环都要计算一趟开方，效率会很低，如果写成 \(i \times i \leq n\) 又很容易溢出，故应改成 \(i \leq n / i\)，这种根据定义的素数判别法，我们称为试除法。试除法也被用在素数个数的求解上面。

bool is_prime(int n){
    // 由于0与1即非质数也非合数，故应特判
    if(n < 2) return false;
    for(int i = 2; i <= n / i; i++){
        if(n % i == 0){
            return false;
        }
    }
    return true;
}

(1.2) 素数筛法

(1.2.1) 朴素筛法

　我们使用一个数组 \(\text{excl}\) 记录那些已被排除在素数可能之外的正整数，也即那些已被筛掉的正整数，反之，未被排除则为素数。筛法的思想很简单，某个数的倍数必然是合数，以此类推，故当遍历第 \(i\) 项的时候，如果这一项仍未被筛掉则必为宿舍，将其纳入 \(prime\) 数组，算法时间复杂度\(O(N\log N)\)。

　然而，这种朴素的做法尚有优化的余地，因为其重复的筛掉了同一个数，例如遇到合数的时候，合数 \(6\) 筛掉的 \(12\) 实际也是素数 \(2\)、素数 \(3\) 倍数，

int excl[MAXN], prime[MAXN], cnt;
int prime_sifter(int n){
  for(int i = 2; i <= n; i++){
        if(!excl[i]) prime[cnt++] = i;
        for(int j = i + i; j <= n; j += i){
            excl[j] = true;
        }
    }
    return cnt;// 1~n 范围之内的素数个数
}

(1.2.2) 埃氏筛法

　朴素方法的改良方案是仅当遇到素数才做一趟筛除操作，素数的倍数必然是一个合数，外层遍历 \(2 \leq i \leq n\)，当遍历遇到素数则将其其倍数筛掉，此法即为埃氏筛法，时间复杂度\(O(N\log\log N)\)，已经非常接近线性了，但是这种做法仍有优化的余地。因为只对素数做筛除仍然会有重复筛除的操作，例如 $ 6=2 \times 3$，遍历遇到素数 \(2\) 与 \(3\) 都会筛一趟数字 \(6\)，根据质因数分解定理克制，任意一个正整数均可写为 \(N=p_1^{\alpha_1}p_2^{\alpha_2}p_3^{\alpha_3}...p_k^{\alpha_k}\) 这种形式，因而我们只需要使用其最小质因数筛它便可避免重复筛除了。下文我们会详细介绍这个质因数分解定理。

int excl[MAXN], prime[MAXN], cnt;
int prime_sifter(int n){
  for(int i = 2; i <= n; i++){
        if(!excl[i]){
            prime[cnt++] = i;
            for(int j = i + i; j <= n; j += i){
                excl[j] = true;
            }
        }
    }
    return cnt;// 1~n 范围之内的素数个数
}

(1.2.3) 欧拉筛法

　只用最小质因数筛除的做法称为欧拉筛，也称线性筛，顾名思义时间复杂度\(O(N)\)。这种做法的思想是构造一单调递增素数序列 \(\text{prime}\)，然后使用最小质因数去筛合数。算法有两处基础细节：

　第一处是内循环判定条件，\(\text{prime[j]} \leq n/i\) 其实是把 \(\text{prime[j] * i} \leq n\) 移项处理了，此处是为了避免溢出，同时我们无需在判断条件中写上 \(j \leq cnt\)，这是最坏的情况就是素数序列 \(\text{prime[0,cnt-2]}\) 之内都找不到 \(i\) 最小质因子，说明此时 \(i\) 必定是一个质数，那么 \(i\) 肯定会在 \(\text{prime[cnt-1]}\) 这个位置出现，然后顺理成章的跳出循环。

　第二处是内循环跳出条件，因为素数序列是单调递增的，\(i \% \text{prime[j]} \neq 0\) 意味着 \(i\) 最小质因数尚未找到， 此时仅能确定 \(\text{prime[j]}\) 就是 \(i*\text{prime[j]}\)这一项的最小质因数，筛之，但当 \(i \% \text{prime[j]} = 0\)，说明已经找到了 \(i\) 最小质因数 \(\text{prime[j]}\)，同时\(\text{prime[j]}\) 也是 \(i*\text{prime[j]}\) 这一项的最小质因子，故在后面的步骤中， \(i*\text{prime[j]}\) 肯定会被 \(\text{prime[j]}\) 筛掉，为了避免重复筛除，此轮中我们直接跳出。

int excl[MAXN], prime[MAXN], cnt;
int prime_sifter(int n){
  for(int i = 2; i <= n; i++){
        if(!excl[i]) prime[cnt++] = i;
        for(int j = 0; prime[j] <= n/i; j++){
            excl[prime[j] * i] = true;
            if(i % prime[j] == 0) break;
        }
    }
    return cnt;// 1~n 范围之内的素数个数
}

(二) 约数

(2.1) 试除法计算约数个数

　给定若干正整数\(A_i\)，对于每个整数\(A_i\)，按照从小到大的顺序输出它的所有约数。采用试除法，我们扫描一趟即可，但当 \(A_i\) 非常大的时候，这种做法非常低效，参考素数判定的优化思路，假如因数 \(d|N\)，那么必有 \(\frac{N}{d}|N\)，因而我们只需扫描 \(2\leq i \leq \sqrt{N}\) 这个范围即可，单个数据\(N\)计算约数的时间复杂度 \(O(\sqrt{N})\)，另外一种计算约数个数的方法是通过质因数分解。

　分解正整数得到各个素数项的指数\(\alpha_i\)，再算 \(\prod_{i=1}^k(1+\alpha_i)\)即可，对于多个正整数累乘结果计算约数个数、约数之和均可采用质因数分解定理。

vector<int> divide(int n){
    vector<int> res;
    for(int i = 1; i <= n / i; i++){
        if(n % i == 0){
            res.push_back(i);
            if(i != n / i) res.push_back(n / i);
        }
    }
    sort(res.begin(), res.end());
    return res;
}

(2.2) 质因数分解应用

　质因数分解定理也称算数基本定理，主要是说任意一个大于一的正整数 \(N\), 均可写为下列这种若干素数幂乘积的形式：

\[N=p_1^{\alpha_1}p_2^{\alpha_2}p_3^{\alpha_3}...p_k^{\alpha_k} \]

　这条式子会在接下来被反复用到，是一个非常重要的定理，相关证明会在文末给出，

　质因数分解通过不停对输入的整数 \(N\) 除以素数并在内层循环中不停变化 \(N\) ，以此得到所有素数序列及其质数序列。由于正整数 \(N\) 至多存在一个因数大于 \(\sqrt{N}\)，在外层循环中，我采用了与试除法一样的优化，只遍历 \([2,\sqrt{N}]\) 这个范围，并在循环结束之后另外特判一下。

// 每个整数 n 都可以写成质因数乘积的形式
void get_prime_factor(int n){
    for(int i = 2; i <= n / i; i++){
        if(n % i == 0){ // 这条判断并非必须的，需要具体情况具体分析，此处仅是为了每算一个素因数打印一次
            int s = 0;
            while(n % i == 0){
                 n = n / i; // 注意 n 是一直都在变化的
                 s = s + 1;
            }
            printf("%d %d\n", i, s);
        }
    }
    if(n > 1){ // 每个数至多存在一个大于sqrt(n)正因子
        printf("%d %d\n", n, 1);//然而此处打印的一项未必大于最初的sqrt(n)，例如的 96=2^5 * 3*1，此处打印 3，显然 3 < sqrt(96)
    }
    printf("\n");
}

(2.2.1) 质因数分解计算约数个数

　给定若干正整数\(A_i\)，试着输出这些数的乘积的约数个数，并对 \(10^9+7\) 取模。其实这个问题，与对单个整数求解约数个数的思路基本是一样的，根据质因数分解定理可知，任意一个正整数 \(N\) 均能被唯一地质因数分解，若干正整数相乘自然也不例外，不妨假设：

\[N=\prod_{i=1}^nA_i=p_1^{\alpha_1}p_2^{\alpha_2}p_3^{\alpha_3}...p_k^{\alpha_k} \]

　\(N\) 约数实际就是在这些项中，任意选出若干个因子组合得到的，对于第 \(i\) 项共有 \((1+\alpha_i)\)种可能选择，因而所有可能的约数即为 \(\prod_{i=1}^n(1+\alpha_i)\)，程序如下。

typedef long long llong;
unordered_map<int, int> hash_in;

int n, x, MOD = 1e9 + 7;
int main(){
    scanf("%d", &n);
    for(int i = 0; i < n; i++){
        scanf("%d", &x);
        for(int j = 2; j <= x / j; j++){
            while(x % j == 0){
                x =  x / j;
                hash_in[j]++;
            }
        }
        if(x > 1) hash_in[x]++;
    }
    
    llong ans = 1;
    for(auto &kv: hash_in){
        ans = ans * (kv.second + 1) % MOD;
    }
    printf("%lld\n", ans);
    return 0;
}

(2.2.2) 质因数分解计算约数总数

　给定若干正整数\(A_i\)，打印所有这些数的乘积的约数之和，并对 \(10^9+7\) 取模。此题如果分别算出每个数据的约数再合并，复杂度会很高。不妨设 \(N=\prod_{i=1}^nA_i\)，同样的根据质因数分解定理可知，

\[N=p_1^{\alpha_1}p_2^{\alpha_2}p_3^{\alpha_3}...p_k^{\alpha_k} \]

　其中约数是由于这些质因子组合而成的，根据排列组合的原理，所有质因数之和\(S\) 能被写成,

\[\begin{aligned} S &=\prod_{i = 1}^k (p_{i}^{0}+p_{i}^{1}+p_{i}^{2}+p_{i}^{3}+...+p_{i}^{\alpha_1})\\ &=\prod_{i = 1}^k\sum_{j=0}^{\alpha_i}p_i^j \end{aligned} \]

　由此可见，我们并不需要知道每个约数具体是什么，只要算出 \(N=\prod_{i=1}^nA_i\) 质因数及其指数之后，只要利用全一系数多项式算法计算每一个项数 \((p_{i}^{0}+p_{i}^{1}+p_{i}^{2}+p_{i}^{3}+...+p_{i}^{\alpha_1})\)，再对项数边乘边模即可。

#include<bits/stdc++.h>
using namespace std;

typedef long long LLong;
unordered_map<int,int> hash_in;
 
int n, x;
int MOD = 1e9 + 7;

int main(){
    scanf("%d", &n);
    for(int i = 0; i < n; i++){
        scanf("%d", &x);
        for(int j = 2; j <= x / j; j++){
            while(x % j == 0){
                x = x / j;
                hash_in[j]++;
            }
        }
        if(x > 1) hash_in[x]++;
    }
    LLong ans = 1;
    for(auto &kv: hash_in){
        LLong p = kv.first, alpha = kv.second, t = 1;
        while(alpha--){
            t = (t * p + 1) % MOD; //此处全一多项式相乘算法，非常巧妙
        }
        ans = ans * t % MOD;
    }
    printf("%lld\n", ans);
    return 0;
}

(2.2.3) 质因数分解计算数论欧拉函数

　对于正整数 \(N\)，欧拉函数返回小于等于 \(N\) 范围中与其互质的正整数个数，记作\(\varphi(N)\)，同时规定 \(\varphi(1) = 1\)。欧拉函数是一个积性函数，也就是说，对于互质的 \(M,N\) 满足 \(\varphi(MN)=\varphi(M)\varphi(N)\)，而对任意 \(M,N\) 则有:

\[\varphi(MN)=\cfrac{\varphi(M)\varphi(N)gcd(M,N)}{\varphi(gcd(M,N))} \]

　
　不过此处我们只需要考虑互质的情况即可，根据质因数分解定理展开正整 \(N\) 会有，

\[\varphi(N)=\varphi(p_1^{\alpha_1}p_2^{\alpha_2}p_3^{\alpha_3}...p_k^{\alpha_k})=\prod_{i=1}^k\varphi(p_i^{\alpha_i}) \]

　对于其中的每项 \(\varphi(p_i^{\alpha_i})\)，分析 \([1,p_i^{\alpha_i}]\) 区间，会发现，其中与正整数 \(p_i^{\alpha_i}\) 不互质的包括了 \(1\times p_i\)、 \(2\times p_i\)、\(3\times p_i\)，... ，\(p_i^{\alpha_i-1}\times p_i\)，剩下的 \(p_i^{\alpha_i}-p_i^{\alpha_i-1}\) 这么多项的正整数，均与其互质。因而 \(\varphi(p_i^{\alpha_i}) = p_i^{\alpha_i}(1-\frac{1}{p_i})\)，代入上述式子即得，

\[\begin{align} \varphi(N)&=\prod_{i=1}^kp_i^{\alpha_i}(1-\frac{1}{p_i})\\ &=\prod_{i=1}^kp_i^{\alpha_i}\prod_{i=1}^k(\frac{p_i-1}{p_i})\\ &=N\prod_{i=1}^k(\frac{p_i-1}{p_i}) \end{align} \]

　经过分析，会发现欧拉函数与正整数分解之后的素数序列有关，而与素数序列的指数项无关，因而只要通过质因数分解算法计算素数序列即可。

int varphi(int n){
    int res = n;
    for(int j = 2; j <=  n / j; j++){
        int s = 0;
        if(n % j == 0){
            while( n % j == 0){
                n = n / j;
                s = s + 1;
            }
            res = res / j * (j - 1);
        }
    }
    return (n > 1) ? (res / n * (n - 1)) : (res);
}

(三) 数论欧拉函数

　虽然上文已经大致的介绍了一遍欧拉函数，但是通过质因数分解法求解单个正整数的欧拉函数，复杂度 \(O(\sqrt{N})\) 尚可接受，但是我们经常需要批量的计算欧拉函数而不仅仅是计算某个数的欧拉函数，给定数列，长度 \(M\)，数值平均大小约为 \(N\)，此时复杂度大概 \(O(M\sqrt{N})\)，非常容易超时。

　其实线性筛，也即欧拉筛法在筛数的过程中，能够顺带的算出欧拉函数、莫比乌斯函数等等，复杂度\(O(M)\)，只要在筛数的过程中分类讨论即可。

int phi[MAXN],  excl[MAXN], prime[MAXN], cnt;
void euler(int n){
    phi[1] = 1;
    for(int i = 2; i <= n; i++){
        if(!excl[i]){
            prime[cnt++] = i;
            phi[i] = i - 1;
        }
        for(int j  = 0; prime[j] <= n/i; j++){
            excl[i * prime[j]] = 1;
            if(i % prime[j] == 0){
                phi[i * prime[j]] =  phi[i] * prime[j];
                break;
            }
            phi[i * prime[j]] = phi[i] * (prime[j] - 1);
        }
    }
}

　分析上面的代码，首先是不被筛掉的数据， \(i\) 不被筛掉说明 \(i\) 是素数，欧拉函数直接等于 \(\varphi(i)=i-1\)，接着是被筛的数据，为图方便，不妨把 \(\text{prime[j]}\) 记作 \(p_j\)，筛数的过程中，如果 \(i \% p_j=0\)，说明 \(p_j\) 是 \(i\) 最小质因子，也就是说，\(i*p_j\) 与 \(i\) 质因数分解之后，素数序列是一样，仅指数项不同，二者展开式如下，

\[\begin{align} i &= p_1^{\alpha_1}...p_j^{\alpha_j}...p_k^{\alpha_k}\\ i*p_j &= p_1^{\alpha_1}...p_j^{\alpha_j+1}...p_k^{\alpha_k} \end{align} \]

　欧拉函数与指数项无关，二者欧拉函数分别如下，

\[\begin{align} \varphi(i)&=i\prod_{i=1}^k(1-\cfrac{1}{p_i})\\ &=i\prod_{i=1}^k(\cfrac{p_i-1}{p_i})\\ \varphi(i*p_j)&=i*p_j\prod_{i=1}^k(1-\cfrac{1}{p_i})\\ &=i*p_j\prod_{i=1}^k(\cfrac{p_i-1}{p_i})\\ &=p_j * \varphi(i) \end{align} \]

　如果 \(i \% p_j \neq 0\) 说明素数 \(p_j\) 与 \(i\) 互质，\(\varphi(i*p_j)=\varphi(i)\varphi(p_j)\)，简单来说，\(i*p_j\) 分解所得的素数序列会比 \(i\) 多出一项 \(p_j\)，也就是说，其欧拉函数等于 \(\varphi(i)p_j(1-\frac{1}{p_j})\)，化简可得，\(\varphi(i)(p_j-1)\)。上述代码便是根据分类讨论的结果写出的。

(四) 组合数计算

(4.1) 动态规划计算组合数

　组合数表达式 \(C_i^j\) 是指从 \(i\) 个对象中，取出 \(j\) 个对象可能构成的组合数，有时为免符号上下标太小影响阅读，我们偶尔也把 \(C_i^j\) 记作 \(C(i,j)\)。数据量比较小的时候，能用动态规划的方式打表算出每个组合数，问询的时候直接查表给出答案即可，这种做法的复杂度达到了 \(O(N^2)\) 级别，数据量稍大一点则无法通过了。

int MOD = 1e9 + 7;
int n, a, b, comb[MAXN][MAXN];

void solve(){
    for(int i = 0; i <= 2000; i++){
        for(int j = 0; j <= i; j++){
            if(!j){
                comb[i][j] = 1;
            }else{
                comb[i][j] = (comb[i - 1][j] + comb[i - 1][j - 1]) % MOD;
            }
        }
    }
}

int main(){
    scanf("%d", &n);
    solve();
    for(int i = 0; i < n; i++){
        scanf("%d %d", &a, &b);
        printf("%d\n", comb[a][b]);
    }
    return 0;
}

(4.2) 利用乘法逆元计算组合数

　由于组合数增长速度非常快，很容易溢出，故常见的组合数计算都是带模的运算，例如给定 \(1\leq b \leq a \leq 10^5\)，要求 \(C_a^b \% p\)，如果\(p\)是一个素数，并且 \(p>a>b\)，通常给定的\(p=10^9+7\)，在这种情况之下，利用乘法逆元即可避免除法运算。

\[C_a^b=\cfrac{a!}{b!(a-b)!}=a!*(b!)^{-1}*((a-b)!)^{-1} \]

　其中 \((b!)^{-1}=(((b-1)!)b)^{-1}=((b-1)!)^{-1}*b^{-1}\)，

　由于 \(p\) 是一个足够大的质数， \(1\text{~}b\)、\(1\text{~}(a-b)\) 之间每一个数与之都不存在\(1\)以外的公因数，根据费马小定理可知，我们利用快速幂预处阶乘的理乘法逆元，然后再根据组合数公式计算即可。

#define MAXN 100005
typedef long long LLong;

LLong n, a, b, MOD = 1e9 + 7;
LLong fact[MAXN], invfact[MAXN];

LLong quick_pow(LLong a, LLong b, LLong p){
    LLong ans = 1;
    while(b){
        if(b & 1){
            ans = ans * a % p;
        }
        a = a * a % p;
        b >>= 1;
    }
    return ans;
}

// 预处理一下
void solve(){
    fact[0] = invfact[0] = 1;
    for(int i = 1; i <= 100000; i++){
        fact[i] = fact[i - 1] * i % MOD;
        invfact[i] = invfact[i - 1] * quick_pow(i, MOD - 2, MOD) % MOD;
    }
}

int main(){
    solve();
    scanf("%lld", &n);
    for(int i = 0; i < n; i++){
        scanf("%lld %lld", &a, &b);
        printf("%lld\n", fact[a] * invfact[b] % MOD * invfact[a - b] % MOD);
    }
    return 0;
}

(4.3) 鲁卡斯定理计算组合数

　如果组合数计算题给出的模数 \(p\) 不满足\(p>a>b\)，甚至 \(p\) 可能是一个合数，此时乘法逆元有可能是不存在的。上述方法失效。对于这种情况，我们能够改用 Lucas 定理，也即鲁卡斯定理，其常被计算大整数的组合数，本条定理主要是说，在对素数取模的条件之下会有，

\[C_i^j = C_{\lfloor i/p \rfloor}^{\lfloor j/p \rfloor} \times C_{\lfloor i\%p \rfloor}^{\lfloor j\%p \rfloor} \mod~p \]

下面借助多项式展开的方式给出证明，对于任意质数\(p\)，会有，

\[(1+x)^p\equiv(1+x^p)\mod p \]

　关于这条结论只要使用二项式展开即可证明，我们不做详细讨论，接下来，不妨记 \(a=\lfloor i/p \rfloor\), \(b = \lfloor i\%p \rfloor\)，\(c = \lfloor j/p \rfloor\)，\(d = \lfloor j\%p \rfloor\)，对于\((1+x)^i\)，我们将其展开，

\[\begin{align} (1+x)^n &=(1+x)^{ap+b}\\ &=(1+x^p)^a(1+x)^b\\ &=\sum_{m=0}^aC_a^mx^{mp}\sum_{n=0}^bC_b^{n}x^n \end{align} \]

　此时考虑\(x^j\) 系数，式子左侧\(C_i^j\)，右侧\(C_a^mC_b^n\)，其中(\(pm+n=j\))，因而 \(C_a^mC_b^n\) 实际就是\(C_a^cC_b^d\)，所有会有 \(C_i^j = C_{\lfloor i/p \rfloor}^{\lfloor j/p \rfloor} \times C_{\lfloor i\%p \rfloor}^{\lfloor j\%p \rfloor}\)，命题得证。根据式子可以发现，鲁卡斯定理能够递归实现。

LLong quick_pow(LLong a, LLong b, LLong p){
    LLong ans = 1;
    while(b){
        if(b&1) ans = ans * a % p;
        a = a * a % p;
        b >>=1;
    }
    return ans;
}

LLong comb(LLong a, LLong b, LLong p){
    LLong ans = 1;
    for(int i = 1, j = a; i <= b; i++, j--){
        ans = ans * j % p;
        ans = ans * quick_pow(i, p - 2, p) % p;
    }
    return ans;
}

LLong lucas(LLong a, LLong b, LLong p){
    if(a < p && b < p)
        return comb(a,b, p);
    return comb(a%p,b%p,p) * lucas(a/p, b/p, p) % p;
}

(4.4) 高精度计算组合数

　如果题目给出的 \(a,b\) 不算特别大，但又没有给出一个模数，其组合数极有可能会溢出，而且是即便长整型变量也无法承受的溢出。这种情况之下，仅能使用高精度算法来计算，为了避免除法运算，我们不妨先对阶乘做一趟质因数分解，再对分解出来的质数进行高精度乘法即可。

#define MAXN 5005
#define PFCT prime_factor
int excl[MAXN], primes[MAXN], cnt;
void get_primes(int n){
    for(int i = 2; i <= n; i++){
        if(!excl[i]) primes[cnt++] = i;
        for(int j = 0; primes[j] <= n / i; j++){
            excl[i * primes[j]] = 1;
            if(i % primes[j] == 0) break;
        }
    }
}

int prime_factor(int n, int p){
    int ans = 0;
    while(n){
        ans += n / p;
        n = n / p;
    }
    return ans;
}

vector<int> mul(vector<int> &a, int b){
    int t = 0;
    vector<int> c;
    for(int i = 0; i < a.size(); i++){
        t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    while(t){
        c.push_back(t % 10);
        t /= 10;
    }
    while(c.size() > 1 && c.back() == 0) c.pop_back();
    return c;
}

int a, b;
int main(){
    // 先做质因数分解,再做高精度相乘
    scanf("%d %d", &a, &b);
    get_primes(a);
    vector<int> ans;
    ans.push_back(1);
    for(int i = 0; i < cnt; i++){
        int alpha = PFCT(a, primes[i]) - PFCT(b, primes[i]) - PFCT(a - b, primes[i]);
        for(int j = 0; j < alpha; j++){
            ans = mul(ans, primes[i]);
        }
    }
    for(int i = ans.size() - 1;  i >= 0; i--){
        printf("%d", ans[i]);
    }
    printf("\n");
    return 0;
}

(五) 常用性质以及些许不严谨证明

　如果 \(a,b\in Z\)，\(b \neq 0\)，如果存在 \(q\in Z\) 使得 \(a=bq\)，则称 \(b\) 整除 \(a\) 或称\(a\) 能被 \(b\) 整除。记作 \(b|a\)，且称 \(a\) 是 \(b\) 倍数，\(b\) 是 \(a\) 因数，下面给出几条非常显然但又非常重要的性质，此处给出的表达并不严谨，若想构建严谨体系的同学请看姜建国的《数论算法》，个人感觉写得非常好，

• 整除具有传递性，\(c|b\)，\(b|a\) 则有 \(c|a\)
• 整除表达式能够进行类似于加减运算的表达，\(c|b\) 且 \(c|a\) 则有 \(c|(a±b)\)
• 整除表达式能够进行类似于线性组合的表达，如果 \(\forall_i\) 均有 \(b|a_i\) 则有 \(b|\sum s_ia_i\)
• 如果 \(b|a\) 且\(a≠0\) 必有\(|b|\leq|a|\)，如果 \(a,b\) 反过来也成立则有 \(|a|=|b|\)
• 对于一个整系数的多项式，若有 \(d|(b-c)\) 则有 \(f(b) - f(c)\)
• 对于素数p，如果 \(p|ab\)，则有 \(p|a\) 或 \(p|b\)

　如何证明我们上面提到质因数分解定理，也即如何证明算数基本定理？

　我们先证存在性，再证唯一性，对于 \(n=2\) 显然成立，假设对于 \(2\leq n<k\) 均成立，则当 \(n=k\)，\(k\) 若是素数显然也成立，若是合数，则其能被写成 \(n=n_1n_2\)，由于 \(n_1<n,n_2<n\) 均为写为质因数分解形式，因而 \(n\) 亦可写为质因数分解形式。

　接下来证明唯一性，假设 \(n\) 存在两种分解，\(n = p_1...p_s = q_1...q_t\)，我们假设这两种分解的素数都是按照从小到大的顺序排列的，

　由于 \(p_1|n\)，因而会有\(p_1|q_1...q_t\)，每个 \(q\) 均为素数，因而必然存在某个\(q_j\) 满足 \(p_1 | q_j\)，反过来也应存在某个\(p_i\) 满足 \(q_1|p_i\)，由于两条分解式都是从小到大排列的，所以我们知道 \(p_i \geq p_1\)，\(q_j\geq q_1\)，又因为某个 \(p\) 与 \(q\) 都是素数，能够整除它们的，仅有\(1\)或其自身，因而 \(p_1=q_j\) 且 \(q_1=p_i\)，根据这个不等关系我们能够列出，

\[\begin{align} p_1&=q_j \geq q_1\\ q_1&=p_i \geq q_1 \end{align} \]

　因而 \(p_1 = q_1\)，类似的，我们能知道 \(p_2 = q_2\)，以此类推，\(p_1...p_s = q_1...q_t\) 每一项均相等，故质因数分解唯一性得证。

支持作者