数论基础
本文内容：
- 二分法实现 sqrt
- 快速幂
- GCD 和 LCM
- 分数四则运算
- 素数表
- 质因子分解
- 大整数运算
- n! 的质因子个数问题
- 组合数求解
- 组合数模 p 求解
- Lucas 定理
1. 求sqrt(2)的值
基于零点定理。
考虑函数 \(f(x)=x^2\) ，给定 \(l=1, r=2\)，则：
l = mid if f(mid)<2
r = mid if f(mid)>2
重复二分直到 l 和 r 足够靠近（即小于设置的精度 EPS）
代码：
#include <iostream>
#include <math.h>
using namespace std;
const double eps = 1e-8;
double square(double x)
{
return x * x;
}
int main()
{
double l = 1, r = 1.5;
double mid;
while ((r - l) >= eps)
{
mid = (l + r) / 2;
if (square(mid) > 2)
{
r = mid;
}
else
{
l = mid;
}
}
cout << mid << endl;
cout << sqrt(2) << endl;
}
2. 快速幂
给定 \(a, b, m\) ，求出 \(a^b\%m\) 的值（模 \(m\) 是为了防止溢出）。
比较暴力的做法就是循环 b 次乘以 a 。
伪代码：
for i=1 to b
{
a = (a*a)%m
}
快速幂解法：
-
递归快速幂：当 \(b\) 是偶数，那么先计算 \(k = a^{b/2}\)，显然对于 \(k\) 可以采取递归处理，那么最后只需要计算 \(k^2\)。如果 \(b\) 是奇数，那么 \(b-1\) 是偶数，即计算 \(a*f(a,b-1,m)\)。算法复杂度为 \(O(logb)\) 。
-
迭代快速幂：将 \(b\) 看作是二进制数，例如 \(b=13=(1101)_2\) ，显然有 \(a^{13} = a^8×a^4×a^1\) ，由于二进制数中每个 bit 代表的十进制都可以由前一个 bit 的十进制计算出来（例如，\(a^8 = a^4 × a^4\)）。迭代快速幂主要思想是遍历 \(b\) 的每一个 bit 位，同时记住当前 bit 的十进制，累积之。算法复杂度 \(O(32)\) 。
需要特别注意的地方：
- 输入
b=1 - 输入
a>m
代码：
/*
Input:a,b,m
Output: a^(b)%m
*/
#include <iostream>
using namespace std;
uint64_t binaryPowerRec(uint64_t a, uint64_t b, uint64_t m)
{
if (b == 0)
return 1;
if (b == 1)
return a;
else if (b & 1)
{
//return (binaryPowerRec(a, b / 2, m) * binaryPowerRec(a, b / 2, m) * a) % m;
return (binaryPowerRec(a, b - 1, m) * a) % m;
}
else
{
int k = binaryPowerRec(a, b / 2, m);
return (k * k) % m;
}
}
uint64_t binaryPowerItor(uint64_t a, uint64_t b, uint64_t m)
{
if (b == 0)
return 1;
uint64_t k = 1;
while (b)
{
if (b & 1)
{
k = (k * a) % m;
}
a = (a * a) % m;
b >>= 1;
}
return k;
}
int main()
{
uint64_t a, b, m;
a = 31, b = 31, m = 32;
a = a % m; //if a>m
cout << binaryPowerRec(a, b, m) << endl;
cout << binaryPowerItor(a, b, m) << endl;
}
3. GCD 和 LCM
GCD 是最大公约数， LCM 是最小公倍数。
GCD 基于欧几里得算法：
证明也比较简单:
| GCD证明 |
|---|
| 设 \(a=kb+r\)，其中 \(k\) 是 \(a\) 除以 \(b\) 的商，\(r\) 是 \(a\) 除以 \(b\) 的余数。 显然有 \(r = a - kb\) 成立。 设 \(d\) 是 \(a\) 和 \(b\) 的一个公约数， 则 \(d\) 也是 \(r = a - kb = a \% b\) 的一个约数。 因此：\(d\) 是 \(a,b\) 的一个公约数，也是 \(b,a\%b\) 的一个公约数。 由于 \(d\) 的任意性，所以 { \(a,b\) 的公约数} \(\subseteq\) { \(b, a\%b\) 的公约数}。 同理可证：{ \(b, a\%b\) 的公约数} \(\subseteq\) { \(a,b\) 的公约数} 因此：{ \(a,b\) 的公约数} \(=\) { \(b, a\%b\) 的公约数} 显然：\(max(\{ a,b 的公约数\}) = max(\{b,a\%b的公约数\})\) 即：\(gcd(a,b) = gcd(b, a\%b)\) |
代码实现：
#include <iostream>
using namespace std;
//最大公约数， require a>b
int gcd(int a, int b)
{
return (b == 0) ? a : gcd(b, a % b);
}
//最小公倍数
int lcm(int a, int b)
{
return a / gcd(a, b) * b;
}
int main()
{
cout << gcd(10, 3) << endl;
cout << lcm(10, 5) << endl;
cout << lcm(lcm(3, 9), 12) << endl;
cout << lcm(-9, -3) << endl;
}
注意 LCM ，如果是 a*b/gcd(a,b) ，a*b 可能会溢出。
4. 分数四则运算
实现分数的四则运算，《算法笔记》一书比较基础（确实不该看这部分内容，太简单了）。就当复习一下C++基础知识吧。
/*
1. 当分数为负，约定分子为负
2. 当分数为0， up=1,down=0
3. 分子分母不可约
*/
#include <iostream>
#include <algorithm>
using namespace std;
//最大公约数
int gcd(int a, int b)
{
return (b == 0) ? a : gcd(b, a % b);
}
//最小公倍数
int lcm(int a, int b)
{
return a / gcd(a, b) * b;
}
class Fraction
{
public:
int up, down;
Fraction(int u = 0, int d = 1)
{
up = u, down = d;
reduction();
}
void reduction()
{
if (down < 0)
{
up = -up, down = -down;
}
if (up == 0)
{
down = 1;
return;
}
int a = abs(up), b = abs(down);
if (a < b)
swap(a, b);
int k = gcd(a, b);
up /= k;
down /= k;
}
void print()
{
cout << up << '/' << down << endl;
}
void operator=(const Fraction &f)
{
up = f.up, down = f.down;
}
Fraction operator+(const Fraction &f) const
{
Fraction c(up * f.down + down * f.up, down * f.down);
return c;
}
Fraction operator-(const Fraction &f) const
{
Fraction c(up * f.down - down * f.up, down * f.down);
return c;
}
Fraction operator*(const Fraction &f) const
{
Fraction c(up * f.up, down * f.down);
return c;
}
Fraction operator/(const Fraction &f) const
{
Fraction c(up * f.down, down * f.up);
return c;
}
};
int main()
{
Fraction a(1, 2), b(1, 4);
Fraction c;
c = a + b, c.print();
c = a - b, c.print();
c = b - a, c.print();
c = a * b, c.print();
c = a / b, c.print();
}
5. 素数表
求解素数的经典算法就不说了。
在这里的要求是输出 \(0 - n\) 范围内的素数，普通算法复杂度为 \(O(n\sqrt{n})\) 。
还有另一类经典算法叫「筛选法」，\(Eratosthenes\) 筛选法和欧拉筛选法。前者复杂度 \(O(nln(n))\) ，后者复杂度为 \(O(n)\) 。
\(Eratosthenes\) 筛选法举例：
2 3 4 5 6 7 8 9 10 (2 是素数，那么去除 2 的倍数)
2 3 5 7 9 (3 是素数，那么去除 3 的倍数)
2 3 5 7 (5 是素数，以此类推)
若果一个数没有被前面的数「筛选掉」（例如上面的 5 ），那么这个数一定是素数。
算法复杂度：
代码：
#include <iostream>
#include <vector>
#define N 100
using namespace std;
vector<int> table;
bool nonprime[N] = {false};
void findPrime()
{
for (int i = 2; i < N; i++)
{
if (!nonprime[i])
{
table.push_back(i);
for (int j = i + i; j < N; j += i)
nonprime[j] = true;
}
}
}
int main()
{
findPrime();
for (int x : table)
cout << x << " ";
}
6. 质因子分解
将一个正整数 \(n\) 写成若干个质数的乘积的形式。
对于任意正整数 \(n\) ，如果它存在 1 和本身之外的因子，那么一定是在 \(\sqrt{n}\) 的左右成对出现。
在这里，就可推出二级结论：一个整数 \(n\) ，如果存在 \([2,n]\) 内的因子，那么最有只有一个质因子大于 \(\sqrt{n}\)。
规定输入与输输出如下：
Input: n = 12
Output: 12 = 2^2 * 3
在这里，使用结构体 struct Factor 存储一个「质因子项」，其包括质因子和它出现的次数。
struct Factor
{
int x;
int count;
Factor(int k = 0, int c = 0)
{
x = k, count = c;
}
};
完整代码：
vector<Factor> decomposeFactor(int n)
{
vector<Factor> result;
int k = sqrt(n);
for (int i = 2; i <= k; i++)
{
if (n % i == 0)
{
Factor f(i, 0);
while (n % i == 0)
{
f.count++;
n /= i;
}
result.push_back(f);
}
}
if (n != 1)
{
result.push_back(Factor(n, 1));
}
return result;
}
int main()
{
int n = 97532468;
cin >> n;
if (n == 1 || n == 0)
{
cout << n << "=" << n << endl;
return 0;
}
auto v = decomposeFactor(n);
int len = v.size();
cout << n << "=";
cout << v.front().x;
if (v.front().count != 1)
{
cout << "^" << v.front().count;
}
for (int i = 1; i < len; i++)
{
cout << "*" << v[i].x;
if (v[i].count != 1)
{
cout << '^' << v[i].count;
}
}
cout << endl;
return 0;
}
7. 大整数运算
实现大整数 \(A \quad op \quad B\) 的运算，比如「9999999999999999999 * 9」这样的等式，比较经典的模拟型的题目。
数据结构：
class BigInteger
{
vector<int> bits;//存放各个位的数字
//实现API
BigInteger operator+(BigInteger &b);
BigInteger operator-(BigInteger &b);
BigInteger operator*(int k);
BigInteger operator/(int k);
}
对于加法和减法，有 2 种处理办法：
- 对于A和B，先「对准长度」，再从低位往高位处理（本次采取此法）
- 对于A和B，先对较短的一部分相加，再处理剩余部分
加减法要特别注意的地方就是：处理进位和借位。
加法代码：
BigInteger operator+(BigInteger &b)
{
// 方法0：开始对准长度
// 方法1：先算短的，再加剩余的
BigInteger c;
int maxlen = max(bits.size(), b.bits.size());
b.bits.resize(maxlen), bits.resize(maxlen);
int t = 0, flag = 0;
for (int i = 0; i < maxlen; i++)
{
t = bits[i] + b.bits[i] + flag;
c.bits.push_back(t % 10);
flag = t / 10;
}
c.bits.push_back(flag);
format(), b.format(), c.format();
return c;
}
减法代码：
BigInteger operator-(BigInteger &b)
{
//require: a>=b
BigInteger c;
int maxlen = max(bits.size(), b.bits.size());
bits.resize(maxlen), b.bits.resize(maxlen);
int t = 0, flag = 0;
for (int i = 0; i < maxlen; i++)
{
t = bits[i] - b.bits[i] - flag;
if (t >= 0)
c.bits.push_back(t), flag = 0;
else
c.bits.push_back(t + 10), flag = 1;
}
format(), b.format(), c.format();
return c;
}
对于乘除法，我们不实现 BigInteger op BigInteger ，而是先实现 BigInteger op int。
对于乘法：
1 4 7
× 3 5
---------
2 4 5
1 4 0
3 5
---------
5 1 4 5
过程：
- \(7 × 35 = 245\)，当前位为 \(5\)，进位 \(24\)
- \(4×35+24=164\)，当前位为 \(4\)，进位为 \(16\)
- \(1×35+16=51\) ，当前位为 \(1\)，进位为 \(5\)
- 当前位为进位 \(5\)
乘法代码：
BigInteger operator*(int k)
{
BigInteger c;
auto &v = bits;
int len = v.size();
int t = 0, flag = 0;
for (int i = 0; i < len; i++)
{
t = v[i] * k + flag;
c.bits.push_back(t % 10);
flag = (t / 10);
}
while (flag)
{
c.bits.push_back(flag % 10);
flag = flag / 10;
}
return c;
}
现在考虑如何实现 BigInteger1 * BigInteger2 :
BigInteger c("0");
for k in BigInteger2.bits: # high bits to low bits
c = (c*10) + (BigInteger1 * k)
其中，BigInteger = BigInteger * int 和 BigInteger = BigInteger + BigInteger 均已实现。
完整代码:
BigInteger operator*(BigInteger b)
{
BigInteger c("0");
BigInteger t;
int len = b.bits.size();
for (int i = len - 1; i >= 0; i--)
{
t = ((*this) * b.bits[i]);
c = (c * 10) + t;
}
return c;
}
下面考虑 BigInteger = BigInteger / int 的实现。
除法的过程：
_0_1_7_6_
7/ 1 2 3 4
7
-----------
5 3
4 9
-----------
4 4
4 2
-----------
2
算法主要思想：
- 1 和 7 比较，不够除，该位商为0，余数 \(r=1\)
- 当前被除数为 \(r=10r+2=12\) ，商为1，余数 \(r=5\)
- 当前被除数为 \(r=10r+3=53\) ，商为7，余数 \(r=4\)
- 当前被除数为 \(r=10r+4=44\) ，商为6，余数 \(r=2\)
除法代码：
BigInteger operator/(int k)
{
BigInteger c;
int len = bits.size();
int r = 0;
auto &v = bits;
c.bits.resize(len);
for (int i = len - 1; i >= 0; i--)
{
r = r * 10 + v[i];
if (r < k)
c.bits[i] = 0;
else
{
c.bits[i] = r / k;
r = r % k;
}
}
c.format();
return c;
}
至于如何实现 BigInteger = BigInteger1 / BigInteger2 ？暂时没想到好的办法，初步想法是把 BigInteger2 分解为 2 个不超过 int的因子，即：BigInteger2 = int1 * int2 ，然后令：
BigInteger c = BigInteger1
c = c / int1
c = c / int2
return c
显然，当除数为一个超过 int 范围的质数的时候，此法无能为力，因此不是一个好的算法。如果有好的想法，欢迎留言。
8. n! 的质因子个数问题
现考虑一个关于 \(n!\) 的问题：求 \(n!\) 中有多少个值为 \(p\) 的质因子。
比如：\(6!\) 中有 4 个质因子 2 。
对于这个问题，比较直观的做法是计算从 [1,n] 的每个数字各有多少个质因子 \(p\) ，然后累计结果。时间复杂度为 \(O(nlogn)\) 。
代码如下：
int calFactor(int n, int p)
{
int ans = 0;
for (int i = 2; i <= n; i++)
{
int t = i;
while (t % p == 0)
{
ans++;
t /= p;
}
}
return ans;
}
现在给出一种 \(O(logn)\) 的解法。
考虑 \(10!\) 中质因子 2 的个数。显然：
- \(10!\) 中有因子 \(2^1\) 的个数为 \(10/2=5\) .
- \(10!\) 中有因子 \(2^2\) 的个数为 \(10/4=2\) .
- \(10!\) 中有因子 \(2^3\) 的个数为 \(10/8=1\) .
所以 \(10!\) 中质因子 2 的个数为 5+2+1=8。
实际上，该过程可以推广为：\(n!\) 中有 $(\frac{n}{p} + \frac{n}{p^2} + \frac{n}{p^3} + ...) $ 个质因子 \(p\) 。其中除法均为向下取整。
(该推论可以从“分数约分”的角度来证明)。
代码实现：
int calFactor2(int n, int p)
{
int ans = 0;
while (n)
{
ans += (n / p);
n /= p;
}
return ans;
}
利用该算法，可以快速计算出 \(n!\) 末尾有多少个零：末尾 0 的个数等于 \(n!\) 中因子 10 的个数，而这又等于质因子 5 的个数（为什么不是质因子 2 的个数呢？因此 4, 6, 8 中也含有 2 ）。亦即计算：calFactor2(n,5)。
9. 组合数求解
众所周知，组合数具有下列性质：
当然我们可以按照式子 1 去暴力求解，时间复杂度为 \(O(n)\) ，但是显然这种方式需要存储三个阶乘数值，容易溢出。
但其实，组合数还有下面这一条性质：
顺便一提该公式的含义：\(C_n^{m}\) 表示从 \(n\) 个里取 \(m\) 个，现在如果第 \(n\) 个必取出，那么有 \(C_{n-1}^{m-1}\) 的取法；如果第 \(n\) 个必然不取出，那么就有 \(C_{n-1}^{m}\) 种取法。
在计算机编程中，我们是十分喜欢递推公式的，因为有递推公式就意味着可以通过循环迭代进行优化。
对于上式的第一项，递归边界是 n==m；对于第二项，递归的边界是 m==0。
递归形式实现：
//方法1：递推公式-递归形式
uint64_t CombineRec(int n, int m)
{
return (m == 0 || m == n) ? 1 : CombineRec(n - 1, m) + CombineRec(n - 1, m - 1);
}
显然，\(C_{n-1}^{m-1}\) 和 \(C_{n-1}^{m}\) 必然包含了重复的计算。（类似于斐波那契数列的过程）
通过二维数组记录计算过程，递归的第二种形式：
//方法1：递推公式-消除重复计算
const int N = 100;
uint64_t table[N][N] = {{0}};
uint64_t combineRec2(int n, int m)
{
if (n == m || m == 0)
return 1;
if (table[n][m])
return table[n][m];
return (table[n][m] = combineRec2(n - 1, m) + combineRec2(n - 1, m - 1));
}
计算所有结果并存放到二维数组：
//方法1：递推公式-迭代形式
uint64_t caltable(int n, int m)
{
for (int i = 0; i < N; i++)
table[i][0] = table[i][i] = 1;
for (int i = 2; i < N; i++)
{
for (int j = 0; j <= (i / 2); j++)
{
table[i][j] = table[i][i - j] = table[i - 1][j] + table[i - 1][j - 1];
}
}
return table[n][m];
}
此外，我们可以组合数定义式有以下规律：
显然可以发现 k*(n-m+i)/i 这个规律。
//方法2：定义式变形
uint64_t combineNumber(int n, int m)
{
uint64_t k = 1;
for (int i = 1; i <= m; i++)
{
k = k * (n - m + i) / i;
}
return k;
}
10. 组合数模 p 求解
基于上一小节的组合数求解，很容易改写出以下算法：
int table[MAX][MAX] = {{0}};
const int p = 13;
int C1(int n, int m)
{
if (m == 0 || m >= n)
return 1;
if (table[n][m])
return table[n][m];
return (table[n][m] = (C1(n - 1, m) + C1(n - 1, m - 1)) % p);
}
int C2(int n, int m)
{
for (int i = 0; i < n; i++)
table[i][0] = table[i][i] = 1;
for (int i = 2; i <= n; i++)
{
for (int j = 0; j <= (i / 2); j++)
{
table[i][j] = table[i][i - j] = (table[i - 1][j] + table[i - 1][j - 1]) % p;
}
}
return table[n][m];
}
上述算法对 \(p\) 的素数性没有要求，且允许 \(p\) 在 int 的范围内即可。
前面写了这么多，总不能到这里就结束了，对吧。下面接着看基于质因数分解的方法。
现在将组合数 \(C_n^{m}\) 进行质因子分解，设分解结果为：
那么：
\(p_i^{c_i} \% p\) 可通过快速幂算法来实现，时间复杂度是 \(O(1)\) 。
那么怎么将 \(C_n^{m}\) 进行质因子分解呢？考虑到 \(C_n^{m} = \frac{n!}{m!(n-m)!}\) ，公式中存在三个阶乘，这就是上面提到的 n!的质因子个数问题 。
遍历不超过 \(n\) 的所有质数 \(p_i\) ，然后计算 \(n!, m!, (n-m)!\) 中分别含质因子 \(p_i\) 的个数 \(x_i, y_i, z_i\)，就可以得到 \(C_n^{m}\) 中质因子 \(p_i\) 的个数为 \(x_i-y_i-z_i\) 。
先来回忆一下需要用到的辅助函数。
- \(n!\) 的质因子 \(p\) 的个数，复杂度 \(O(log_p{n})\)
int countFactorial(int n, int p)
{
int k = 0;
while (n)
{
k += n / p;
n /= p;
}
return k;
}
- 快速幂，复杂度 \(O(1)\)
uint64_t power(int a, int b, int mod)
{
if (b == 0)
return 1;
uint64_t n = 1;
while (b)
{
if (b & 1)
n = (n * a) % mod;
b >>= 1;
a = (a * a) % mod;
}
return n;
}
- 计算 \(n\) 以内的素数表，复杂度 \(O(nln(n))\)
vector<int> getPrimeTable(int n)
{
bool nonprime[MAX] = {0};
vector<int> v;
for (int i = 2; i <= n; i++)
{
if (!nonprime[i])
{
v.push_back(i);
for (int j = i + i; j <= n; j += i)
nonprime[j] = true;
}
}
return v;
}
- 组合数模 \(p\)
uint64_t C3(int n, int m)
{
vector<int> primes = getPrimeTable(n);
int x, y, z;
uint64_t ans = 1;
for (int p : primes)
{
x = countFactorial(n, p);
y = countFactorial(m, p);
z = countFactorial(n - m, p);
ans = ans * power(p, x - y - z, MODP) % MODP;
}
return ans;
}
假设素数表里面有 \(k(k≤n)\) 个素数。
现在来计算一下总的时间复杂度：\(O(nln(n)) + O(klog_p{n}) = O(nlogn)\)
11. Lucas 定理
Lucas 定理其实还是用于求解组合数模 \(p\) 的问题。
如果 \(p\) 是素数，将 \(m,n\) 表示为 \(p\) 进制：
Lucas 定理：
例如，对于 \(C_8^3 \% 5\) 来说，\(n=8, m=3, p=5\)，那么：
所以：
代码实现：
uint64_t lucas(int n, int m)
{
if (m == 0)
return 1;
return (C3(n % MODP, m % MODP) * lucas(n / MODP, m / MODP)) % MODP;
}
对于 C3(n % MODP, m % MODP) 这一项，显然是对应于公式的 \(C_{n_0}^{m_0}\)，其次是递归项 lucas(n / MODP, m / MODP) 表示的是 \(p\) 进制表示中，从低位到高位，可利用 2 进制这一特殊情况去代入理解。
