【算法竞赛备赛】数论

图片无法加载，请点击此处

前言

在考试中，经常会出现一类问题，它们不涉及很深的算法，却和数学息息相关。这样的问题通常难度不大，也不需要特别的数学知识，只要掌握简单的数理逻辑即可。——《算法笔记》胡凡

质数(素数)

质数又称素数，是指除了1和本身外，不能被其它的任何数整除。

质数的判定

如何判断给定的正整数n是否是质数？

时间复杂度：$O(sqrt(n))$

bool isPrime(int n) {
    if (n < 2) return false;
    
    for (int i = 2; i <= n / i; i++) // i <= n / i
        if (n % i == 0) return false;
    
    return true;
}

分解质因数

n中有且仅有一个大于sqrt(n)的质因子。（如果有两个，两者相乘将大于n，明显不符合）

据此，我们可以先枚举小于sqrt(n)的所有质因子

如果最后n仍大于1，那么这个n就是那个大于sqrt(n)的质因子

最坏时间复杂度：$O(sqrt(n))$

最好时间复杂度：$O(log(n))$，（若n为$2^k$，那么在i=2时，所有因子就被分解出来了，所以时间应该是$log(2^k)=k$）

void divide(int n) {
    
    for (int i = 2; i <= n / i; i++) 
        if (n % i == 0) {
            int s = 0;
            while (n % i == 0) {
                s++;
                n /= i;
            }
            printf("%d %d\n", i, s);
        }
    
    if (n > 1) printf("%d %d", n, 1);
}

筛质数

筛选出从1-n之间质数的个数

数据量小于等于$10^6$时，以下两个算法效率不相上下

埃氏筛法思想

一个数是p，如果2~p-1中的质数都不是它的质因数的话，那么p是质数(以此对算法进行优化)

朴素筛法

$O(nlog(n))$

时间复杂度计算，当n=2时，执行n/2次；当n=3时，执行n/3次

$\frac{n}{2}+\frac{n}{3}+\frac{n}{4}+...+\frac{n}{n}$

$=n\ast (\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n})$

$=n\ast {\lim }_{x\to +\infty }\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n}$

$=n\ast ln(n)<n\ast log(n)=nlog(n)$

int primes[N];
bool st[N];

void get_primes(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) primes[cnt++] = i;  
        
        for (int j = i + i; j <= n; j += i) st[j] = true;
    }
}

cout << cnt << endl; // 输出质数个数

优化

$O(nlog(log(n)))$。相当于线性的时间复杂度

1~n中有$\frac{n}{ln(n)}$个质数

int primes[N];
bool st[N];

void get_primes(int n) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) { 
            primes[cnt++] = i;
            // 是质数的倍数的话，判定为合数
            for (int j = i + i; j <= n; j += i) st[j] = true;
        }
    }
}

cout << cnt << endl; // 输出质数个数

线性筛法(*常用)

n只会被最小质因子筛掉

1. i % primes[j] == 0

   $primes[j]$一定是$i$的最小质因子，$primes[j]$一定是$primes[j]\ast i$的最小质因子

2. i % primes[j] != 0

   $primes[j]$一定小于$i$的所有质因子，$primes[j]$也一定是$primes[j]\ast i$的最小质因子

int primes[N]; // 用于存储质数
bool st[N]; // false时为质数

void get_primes(int n ) {
    for (int i = 2; i <= n; i++) {
        if (!st[i]) primes[cnt++] = i;
        for (int j = 0; primes[j] <= n / i; j++) {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break; // primes[j]一定是i的最小质因子
        }
    }
}

约数

12的约数为1,12,3,4,2,6

试除法求约数

$O(sqrt(n))$

vector<int> get_divisors(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); // 避免出现5*5的情况
        }
        
    sort(res.begin(), res.end());
    return res;
}

约数个数

若要求出N的约数个数，首先将N质因数分解，p代表质数，$\alpha$为p的次数

$N=p_1^{{\alpha }_1}+p_2^{{\alpha }_2}+p_3^{{\alpha }_3}+...+p_n^{{\alpha }_n}$

然后求出每个质因数的约数个数，如

$p_1^{{\alpha }_1}:p_1^0,p_1^1,...,p_1^{{\alpha }_1-1},p_1^{{\alpha }_1}$

所以$p_1^{{\alpha }_1}$共${\alpha }_1+1$个约数，同理$p_2^{{\alpha }_2}$共${\alpha }_2+1$个约数……

根据乘法原则，N的约数个数为$({\alpha }_1+1)({\alpha }_2+1)\ast ...\ast ({\alpha }_n+1)$

---

约定$n$个整数$a_i$，请你输出这些数的乘积的约数个数，答案对$10^9+7$取模

数据范围

$1\le n\le 100$

$1\le a_i\le 2\ast 10^9$

输入格式

3
2
6
8

输出格式

12

#include <iostream>
#include <algorithm>
#include <unordered_map>

using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;
int n;

int main() {
    cin >> n;
    
    unordered_map<int, int> primes;
    while (n -- ) {
        int x;
        cin >> x;
        
        // 分解质因数
        for (int i = 2; i <= x/i; i++) 
            while (x % i == 0) {
                x /= i;
                primes[i]++;
            }
        // 记得收尾
        if (x > 1) primes[x]++;
    }
    
    LL res = 1;
    // 代约数个数公式
    for (auto prime: primes) res = res * (prime.second + 1) % mod;
    
    cout << res << endl;
    
    return 0;
}

约数之和

若要求出N的约数之和，首先将N质因数分解，p代表质数，$\alpha$为p的次数

$N=p_1^{{\alpha }_1}+p_2^{{\alpha }_2}+p_3^{{\alpha }_3}+...+p_n^{{\alpha }_n}$

约数之和公式为

$(p_1^0+p_1^1+...+p_1^{{\alpha }_1})(p_2^0+p_2^1+...+p_2^{{\alpha }_2})...(p_n^0+p_n^1+...+p_n^{{\alpha }_n})$

至于为什么是这样计算呢？可以利用乘法分配律进行验证，两两相乘后可以得到N的所有约数(因为这些质因数本就是由约数拆解而来，这步相当于还原了用质数还原了约数)，进行加和就是答案

---

约定$n$个整数$a_i$，请你输出这些数的乘积的约数之和，答案对$10^9+7$取模

#include <iostream>
#include <algorithm>
#include <unordered_map>

using namespace std;
typedef long long LL;
const int mod = 1e9 + 7;

int n;

int main() {
    cin >> n;
    
    unordered_map<int,int> primes;
    // 质因数分解
    while (n -- ) {
        int x;
        cin >> x;
        
        for (int i = 2; i <= x/i; i++)
            while (x % i == 0) {
                x /= i;
                primes[i]++;
            }
        if (x > 1) primes[x]++;
    }
    
    LL res = 1;
    // 约数之和公式
    for(auto prime: primes) {
        int p = prime.first, a = prime.second;
        LL t = 1;
        while(a--) t = (t * p + 1) % mod;
        res = (res * t) % mod;
    }
    
    cout << res << endl;
    
    return 0;
}

最大公约数

欧几里得算法，即辗转相除法，是将两个数辗转相除直到余数为0

$O(log(n))$

int gcd(int a, int b) {
    return b? gcd(b, a % b):a;
}

欧拉函数

在数论，对正整数n，欧拉函数是小于等于n的正整数中与n互质的数的数目

互质：若N个整数的最大公因数是1，则称这N个整数互质。

若要求出N的欧拉函数，首先将N质因数分解，p代表质数，$\alpha$为p的次数

$N=p_1^{{\alpha }_1}+p_2^{{\alpha }_2}+p_3^{{\alpha }_3}+...+p_n^{{\alpha }_n}$

那么欧拉函数公式如下：

$\varphi (n)=N\ast (1-\frac{1}{p_1})\ast (1-\frac{1}{p_2})\ast ...\ast (1-\frac{1}{p_k})$

证明如下（容斥原理）：

N同时代表个数，首先我们减去N中所有质因数倍数的个数

$N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}$

但这样减会出现一个问题，有的数同时是p1、p2的倍数，就会被减去两次，所以我们需要加回来

$N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1\ast p_2}+\frac{N}{p_2\ast p_3}+...$

一些数同时是p1、p2、p3的倍数，那么减去3次的同时，又加上了3次，相当于没有计算，因此

$N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1\ast p_2}+\frac{N}{p_2\ast p_3}+...-\frac{N}{p_1\ast p_2\ast p_3}-\frac{N}{p_2\ast p_3\ast p_4}$

根据以上式子，可以看出一定的规律，将它们进行整理即可得到欧拉函数

求一个数的欧拉函数

int phi(int x)  // 欧拉函数
{
    int res = x;
    for (int i = 2; i <= x / i; i ++ )
        if (x % i == 0)
        {
            res = res / i * (i - 1);
            while (x % i == 0) x /= i;
        }
    if (x > 1) res = res / x * (x - 1);

    return res;
}

筛选法求欧拉

给定一个数n，求出1-n中每个数欧拉函数的和。这里使用到了前面线性筛质数的方法

1. 若i为质数，那么i的欧拉函数$\varphi (i)$为i-1

2. 若i不为质数，求$\varphi (primes[j]\ast i)$分两种情况

   1. 若i % primes[j] == 0，意味着primes[j]是i的一个质因子

      此时$\varphi (i)=i\ast (1-\frac{1}{p_1})\ast (1-\frac{1}{p_2})\ast ...\ast (1-\frac{1}{p_k})$

      则$\varphi (primes[j]\ast i)=primes[j]\ast i\ast (1-\frac{1}{p_1})\ast (1-\frac{1}{p_2})\ast ...\ast (1-\frac{1}{p_k})$

      对比两个式子发现，后者只比前者多乘了一个$primes[j]$，因此我们得到公式为

      $\varphi (primes[j]\ast i)=primes[j]\ast \varphi (i)$

   2. 若i % primes[j] != 0，意味着意味着primes[j]不是i的一个质因子

      此时$\varphi (i)=i\ast (1-\frac{1}{p_1})\ast (1-\frac{1}{p_2})\ast ...\ast (1-\frac{1}{p_k})$

      则$\varphi (primes[j]\ast i)=primes[j]\ast i\ast (1-\frac{1}{p_1})\ast (1-\frac{1}{p_2})\ast ...\ast (1-\frac{1}{p_k})\ast (1-\frac{1}{p_j})$

      对比两个式子，因为此时的primes[j]不是i的一个质因子，所以在求primes[j]*i时，需要比原先式子多乘了primes[j]和(1-1/pj)，化简得最后结果如下

      $\varphi (primes[j]\ast i)=p[j]\ast \varphi (i)\ast (1-\frac{1}{p_j})=\varphi (i)\ast (p_j-1)$

int get_eulers(int n) {
    
    phi[1] = 1;
    for (int i = 2; i <= n; i++) {
        if (!st[i]) {
            primes[cnt++] = i;
            phi[i] = i - 1;
        }
        for (int j = 0; primes[j] <= n / i; j++) {
            st[primes[j]*i] = true;
            if (i % primes[j] == 0) {
                phi[i*primes[j]] = primes[j] * phi[i];
                break;
            }
            phi[i*primes[j]] = (primes[j]-1) * phi[i];
        }
    }
    LL res = 0;
    for (int i = 1; i <= n; i++) res += phi[i];
    
    return res;
}

欧拉数论定理

内容

若$a$与$n$互质（gcd(a,n)=1），则$a^{\phi(n)} \equiv 1 ; mod ; n $

其中$\varphi (n)$称为对模$n$缩系的元素个数

• 当n为质数时，$a^{n-1} \equiv 1 ; mod ; n $（费马小定理）

证明

1~n中，与n互质的数为$a_1,a_2,a_3,...,a_{\varphi (n)}$（注意这里的a与前面的a不同含义）

将上述的数同时乘以a，得

$a{a}_1,a{a}_2,a{a}_3,...,a{a}_{\varphi (n)}$

因a1,a2,a3与n互质并且a与n也互质，所以两个数的乘积也和n互质。将两组分别内部同时相乘，得

$a_1\ast a_2\ast a_3\ast ...\ast a_{\varphi (n)}\equiv a{a}_1\ast a{a}_2\ast a{a}_3\ast ...\ast a{a}_{\varphi (n)}$

$a_1\ast a_2\ast a_3\ast ...\ast a_{\varphi (n)}\equiv a^{\varphi (n)}(a_1\ast a_2\ast a_3\ast ...\ast a_{\varphi (n)})$

约去相同部分，得

$a^{\varphi (n)}\equiv 1(modn)$

因此我们得出结论

$a^{\phi(n)} \equiv 1 ; mod ; n $

应用

这个定理可以用来简化幂的模运算。比如计算$7^{222}$的个位数，实际是求$7^{222}$被10除的余数。7和10互质，且φ(10)=4。由欧拉定理知$7^4\equiv 1(mod10)$。所以$7^{222}=(7^4{)}^{55}\ast 7^2\equiv 49\equiv 9(mod10)$

快速幂

$1\le a,k,p\le 10^9$

以$log(k)$的时间复杂度，快速求出$a^{k}modp$的结果

• 要快速求出上述的结果，我们需要先预处理出$a^{2^0},a^{2^1}...(modp)$的值
• 然后将$a^k$进行二进制分解，$a^k=a^{2^0}\ast a^{2^1}\ast ...=a^{2^0+2^1+...}(modp)$
• 最后用原先预处理的结果，相乘取模得出答案

此算法主要的时间花在了预处理$O(log(k))$和k的填二进制拆解$O(log(k))$上

虽然上述过程有一些复杂，但在最后的算法实现非常简洁

// a^k%p
int qmi(int a, int k, int p) {
    int res = 1;
    
    // 同时处理了两个过程
    while(k) {
        if (k & 1) res = (LL)res * a % p;
        k >>= 1;
        a = (LL)a * a % p;
    }
    
    return res; 
}

扩展欧几里得算法

求出一组$x_i,y_i$，使其满足$a_i\ast x_i+b_i\ast y_i=gcd(a_i,b_i)$

推导如下：

$by+(amodb)x=d$

$\Rightarrow by+(a-[\frac{a}{b}]\ast b)x=d$

$\Rightarrow ax+b(y-[\frac{a}{b}]\ast x)=d$

int exgcd(int a, int b, int& x, int& y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    } 
    
    int d = exgcd(b, a%b, y, x);
    y -= a / b * x;
    
    return d;
}

线性同余方程

题目描述

给定$n$组数据$a_i,b_i,m_i$，对于每组数求出一个$x_i$，使其满足$a_i\ast x_i\equiv b_i(mod{m}_i)$

推导

$a_i\ast x_i\equiv b_i(mod{m}_i)$

$\Rightarrow a_i\ast x_i/m_i=y_i...b_i$

$\Rightarrow a_i\ast x_i-m_i\ast y_i=b_i$

令$y_i\prime =y_i$

$\Rightarrow a\ast x+m\ast y\prime =b$

如果b不是最大公约数d的倍数，那么一定无解

根据扩展欧几里得算法，对于任意的整数a和m，必定存在一组x和y'，使得ax+my'等于gcd(a,m)

所以可以求得$a\ast x+m\ast y\prime =d$

$\Rightarrow a\ast x\ast \frac{b}{d}+m\ast y\prime \ast \frac{b}{d}=b$

所以要求解的xi=x*b/d

中国剩余定理

中国剩余定理(Chinese remainder theorem, CRT)

小学数学时大家都学过这样的问题：有一群人，3个人一组，剩1人，5个人一组，剩3人，10人一组，剩7人。

那么有多少人呢？CRT的出现就是用来解决上述问题

以上的问题可以看成一个线性同余方程组

CRT公式：$x=({\sum }_{i=1}^n{a}_i\ast M_i^{-1}\ast M_i)modM$

以上问题可以转换成以下的线性同余方程组表达

$\{\begin{array}{l}x\equiv a_1(mod{m}_1)\\ x\equiv a_2(mod{m}_2)\\ \cdots \cdots \\ x\equiv a_k(mod{m}_k)\end{array}$

$M=m_1\ast m_2\ast \cdots \ast m_k$

$M_i=M/m_i=\cdots \ast m_{i-1}\ast m_{i+1}\ast \cdots$

此外，还需要求出$M_i$的逆元$M_i^{-1}$，使其满足条件$M_i\ast M_i^{-1}\equiv 1(mod{m}_i)$

而对于本组中的其它i，满足条件$M_j\ast M_j^{-1}\equiv 0(mod{m}_i)$

根据$x=a_i\ast M_i^{-1}\ast M_i$，求出其中的一个值

对每一线性同余方程执行上述过程，最后将结果对M取模，得出最后答案

---

举个例子，三三数之剩二，五五数之剩三，七七数之剩二。问物几何？

即，一个整数除以3余2，除以5余3，除以7余2，求这个整数。

$M=3\ast 5\ast 7=105$

首先求解$x\equiv 2(mod3)$

$M_1=M/m_1=105/3=35$

$M_1^{-1}\ast 35\equiv 1(mod3)$，此步利用扩展欧几里得算法，求出$M_1^{-1}=2$

$x=a_i\ast M_i^{-1}\ast M_i=2\ast 2\ast 35=2\ast 70=140$

同理，可以求出$x\equiv 3(mod5)$和$x\equiv 2(mod7)$

$x=15\ast 2+21\ast 3+70\ast 2=233modM=233mod105=23$

这个余数23就是合乎条件的最小数.

扩展CRT

根据题目所述，这里的$a_1,a_2,\cdots ,a_n$不满足两两互质的条件，所以需要重新推导CRT

$\{\begin{array}{l}x\equiv m_1(mod{a}_1)\\ x\equiv m_2(mod{a}_2)\end{array}$

化简得，

$\{\begin{array}{l}x=k_1{a}_1+m_1\\ x=k_2{a}_2+m_2\end{array}$

所以$k_1{a}_1+m_1=k_2{a}_2+m_2$

$\Rightarrow k_1{a}_1-k_2{a}_2=m_2-m_1$，这个式子有解则满足$gcd(a_1,a_2)|m_2-m_1$

运用扩展欧几里得算法，能够求出k1,k2的通解,非齐次=齐次通解+非齐次特解

$\{\begin{array}{l}{k}_1=k_1+k\frac{a_2}{d}\\ k_2=k_2+k\frac{a_1}{d}\end{array}$

$x=k_1{a}_1+m_1$

$\Rightarrow x=(k_1+k\frac{a_2}{d})a_1+m_1$

$\Rightarrow x=a_1{k}_1+m_1+k\frac{a_1{a}_2}{d}$

其中a1a2/d就是a1a2的最小公倍数，因此

$\Rightarrow x=a_1{k}_1+m_1+k[a_1{a}_2]$

通过对比最前面化简后的式子，我们可以把$a_1{k}_1+m_1$看成m，$k[a_1{a}_2]$看成ka

$\Rightarrow x=x_0+ka$

以这样的形式重复更新合并后来的线性同余方程，可以得到最终答案

#include <iostream>
#include <cmath>

using namespace std;

typedef long long LL;
int n;

LL exgcd(LL a,LL b, LL& x, LL& y) {
    if (!b) {
        x = 1, y = 0;
        return a;
    }
    
    LL d = exgcd(b, a%b, y, x);
    y -= a / b * x;
    return d;
}


int main() {
    cin >> n;
    
    LL a1, m1;
    cin >> a1 >> m1;
    
    bool has_answer = true;
    for (int i = 0; i < n-1; i ++ ) {
        LL a2, m2;
        cin >> a2 >> m2;
        
        LL k1, k2;
        LL d = exgcd(a1, a2, k1, k2);
        if ((m2-m1) % d) { // 当m2-m1不是d的倍数时，无解
            has_answer = false;
            break;
        }
        k1 *= (m2 - m1)/d; // 欧几里得求出的解是以最小公约数d为结果，转换为m2-m1
        LL t = a2 / d;
        k1 = (k1 % t + t) % t; // 最小正整数解
        
        m1 = a1*k1+m1; // 为下一轮循环计算m1,a1
        a1 = abs(a1 / d * a2); // 不能a1*a2就近相乘，否则会溢出
    }
    
    if (has_answer) cout << (m1 % a1 + a1) % a1 << endl;
    else puts("-1");
    
    return 0;
}

204.表达整数的奇怪方式

高斯消元

• 枚举每一列
  1. 找到该列中绝对值最大的那个数所在的行
  2. 将该行移动到最上面（次上、次次上）
  3. 将该行的第一个数变为1
  4. 用这个数把下面行该列的数消成0

组合计数

容斥定理

简单博弈论

参考资料

AcWing算法基础课，当你还在犹豫的时候，别人已经学完了