开发: C++知识库 Java知识库 JavaScript Python PHP知识库人工智能区块链大数据移动开发嵌入式开发工具数据结构与算法开发测试游戏开发网络协议系统运维
教程: HTML教程 CSS教程 JavaScript教程 Go语言教程 JQuery教程 VUE教程 VUE3教程 Bootstrap教程 SQL数据库教程 C语言教程 C++教程 Java教程 Python教程 Python3教程 C#教程
数码: 电脑笔记本显卡显示器固态硬盘硬盘耳机手机 iphone vivo oppo 小米华为单反装机图拉丁

-> 数据结构与算法 -> 算法竞赛备赛（六）：数论 -> 正文阅读

[数据结构与算法]算法竞赛备赛（六）：数论

作者:_14

title: 机试-数论
date: 2022-2-28 15:03:06
tags:

算法
C++
数论
categories:
算法
cover: graph.png
copyright_author: codeslogan
copyright_url: https://codeslogan.github.io

前言

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

质数(素数)

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

质数的判定

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

时间复杂度： $O (s q r t (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 (s q r t (n))$

最好时间复杂度： $O (l o g (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 (n l o g (n))$

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

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

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

$=n*\lim_{x \to +\infty}{\frac{1}{2}+\frac{1}{3}+\frac{1}{4}+...+\frac{1}{n}}$

$= n ? l n (n) < n ? l o g (n) = n l o g (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 (n l o g (l o g (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只会被最小质因子筛掉

i % primes[j] == 0

$p r i m e s [j]$ 一定是 $i$ 的最小质因子， $p r i m e s [j]$ 一定是 $p r i m e s [j] ? i$ 的最小质因子
i % primes[j] != 0

$p r i m e s [j]$ 一定小于 $i$ 的所有质因子， $p r i m e s [j]$ 也一定是 $p r i m e s [j] ? 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 (s q r t (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)*...*(\alpha_n+1)$

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

数据范围

$\le n \le 100$

$\le a_i \le 2*10^9$

输入格式

输出格式

#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 (l o g (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}$

那么欧拉函数公式如下：

$\phi(n)=N*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(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*p_2}+\frac{N}{p_2*p_3}+...$

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

$N-\frac{N}{p_1}-\frac{N}{p_2}-...-\frac{N}{p_k}+\frac{N}{p_1*p_2}+\frac{N}{p_2*p_3}+...-\frac{N}{p_1*p_2*p_3}-\frac{N}{p_2*p_3*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中每个数欧拉函数的和。这里使用到了前面线性筛质数的方法

若i为质数，那么i的欧拉函数 $\phi(i)$ 为i-1
若i不为质数，求 $\phi(primes[j]*i)$ 分两种情况
1. 若i % primes[j] == 0，意味着primes[j]是i的一个质因子
  
  此时 $\phi(i)=i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})$
  
  则 $\phi(primes[j]*i)=primes[j]*i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})$
  
  对比两个式子发现，后者只比前者多乘了一个 $p r i m e s [j]$ ，因此我们得到公式为
  
  $\phi(primes[j]*i)=primes[j]*\phi(i)$
2. 若i % primes[j] != 0，意味着意味着primes[j]不是i的一个质因子
  
  此时 $\phi(i)=i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})$
  
  则 $\phi(primes[j]*i)=primes[j]*i*(1-\frac{1}{p_1})*(1-\frac{1}{p_2})*...*(1-\frac{1}{p_k})*(1-\frac{1}{p_j})$
  
  对比两个式子，因为此时的primes[j]不是i的一个质因子，所以在求primes[j]*i时，需要比原先式子多乘了primes[j]和(1-1/pj)，化简得最后结果如下
  
  $\phi(primes[j]*i)=p[j]*\phi(i)*(1-\frac{1}{p_j})=\phi(i)*(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 $

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

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

证明

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

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

$aa_1,aa_2,aa_3,...,aa_{\phi(n)}$

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

$a_1*a_2*a_3*...*a_{\phi(n)} \equiv aa_1*aa_2*aa_3*...*aa_{\phi(n)}$

$a_1*a_2*a_3*...*a_{\phi(n)} \equiv a^{\phi(n)}(a_1*a_2*a_3*...*a_{\phi(n)})$

约去相同部分，得

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

因此我们得出结论

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

应用

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

快速幂

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

以 $l o g (k)$ 的时间复杂度，快速求出 $a^k \ mod\ p$ 的结果

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

此算法主要的时间花在了预处理 $O (l o g (k))$ 和k的填二进制拆解 $O (l o g (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*x_i+b_i*y_i=gcd(a_i,b_i)$

推导如下：

$by+(a\ mod\ b)x=d$

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

$\Rightarrow ax+b(y-[\frac{a}{b}]*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*x_i \equiv b_i(mod\ m_i)$

推导

$a_i*x_i \equiv b_i(mod\ m_i)$

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

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

令 $y_i\prime = y_i$

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

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

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

所以可以求得 $y\prime =d$

$\Rightarrow a*x*\frac{b}{d} + m*y\prime*\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*M_i^{-1}*M_i)\ mod \ M$

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

$\begin{cases} x\equiv a_1(mod \ m_1)\\ x\equiv a_2(mod \ m_2) \\ \cdots \cdots \\ x\equiv a_k(mod \ m_k)\end{cases}$

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

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

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

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

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

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

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

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

$M = 3 ? 5 ? 7 = 105$

首先求解 $\equiv 2 (mod \ 3)$

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

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

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

同理，可以求出 $\equiv 3 (mod \ 5)$ 和 $\equiv 2 (mod \ 7)$

$\ mod \ M = 233 \ mod \ 105 = 23$

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

扩展CRT

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

$\begin{cases} x\equiv m_1(mod \ a_1)\\ x\equiv m_2(mod \ a_2) \\ \end{cases}$

化简得，

$\begin{cases} x= k_1a_1+m_1\\ x= k_2a_2+m_2 \\ \end{cases}$

所以 $k_1a_1+m_1=k_2a_2+m_2$

$\Rightarrow k_1a_1-k_2a_2=m_2-m_1$ ，这个式子有解则满足 $gcd(a_1,a_2) | m_2-m_1$

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

$\begin{cases} k_1= k_1+k\frac{a_2}{d}\\ k_2= k_2+k\frac{a_1}{d} \\ \end{cases}$

$x=k_1a_1+m_1$

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

$\Rightarrow x=a_1k_1+m_1+k\frac{a_1a_2}{d}$

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

$\Rightarrow x=a_1k_1+m_1+k[a_1a_2]$

通过对比最前面化简后的式子，我们可以把 $a_1k_1+m_1$ 看成m， $k[a_1a_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;
}