n 次方根取整 求 `|__ root n a __|`, 其中整数 `a ge 0`, `n ge 2`. 使用牛顿迭代法 `x = ((n-1)x + a//x^(n-1))//n`, 其中除法均向下取整.

说明: `235^2 = 55225`, `236^2 = 55696`, 因此 `|__ sqrt 55555 __| = 235`

快速幂 (反复平方法) 求 `a^n`, 先将指数 `n` 写为二进制: `n = (n_k cdots n_1 n_0)_2` `= sum_(i=0)^k n_i 2^i`, 于是 `a^n = prod_("where "n_i = 1) a^(2^i)`.

在算法每次迭代的开头, 变量 `a = 5, 5^2, 5^4, 5^8`, 对应到 `14` 的二进制表示 `0*1 + 1*2 + 1*4 + 1*8` 中系数为 `1` 的项, 结果为 `5^14 = 5^2 * 5^4 * 5^8`.

辗转相除法 设 a, b 为正整数, 求最大公约数 d, 同时求出系数 x, y 满足 a x + b y = d.

75 = (1, 0)
32 = (0, 1)
11 = (1, -2)
10 = (-2, 5)
1 = (3, -7)

求模 `n` 乘法逆元 若 `a x -= 1 (mod n)`, 则称 `x` 为 `a` 模 `n` 的乘法逆元.
若 `(a, n) = 1`, 则存在 `x, y` 使得 `a x + n y = 1`, 即 `a x -= 1 (mod n)`. 因此应用 gcdExtended 就可获得乘法逆.

线性求逆 求 `1` 到 `p-1` 的各整数模素数 `p` 的乘法逆元.
设 `0 lt a lt p`, `x = p mod a`, `y = |__p // a__|`. 则 `0 lt x lt p`, 说明 `x` 模 `p` 有逆元. 计算知 `x + y a = p -= 0 (mod p)`, `x -= -y a (mod p)`, `x a^-1 -= -y (mod p)`, `a^-1 = x^-1 (p-y) (mod p)`.

前缀积求逆 求整数 `a_1, cdots, a_n` 模素数 `p` 的乘法逆元, 其中每个整数都与 `p` 互素. 这是求逆元的一种离线算法. 记第 `i` 个前缀积为 `"pre"_i -= prod_(j=1)^i a_j (mod p)`. 注意若 `(a, p) = (b, p) = 1`, 则 `a^-1 b^-1 -= (a b)^-1 (mod p)`. 于是前缀积的逆等于逆的前缀积, 即 `"pre"_i^-1 -= prod_(j=1)^i a_j^-1 (mod p)`. 利用这一性质设计算法如下: 先使用 gcdExtended 计算所有数的乘积模 `p` 的逆元 `"pre"_n^-1`, 然后递推 `a_i^-1 -= "pre"_i^-1 * "pre"_(i-1) (mod p)`,
`"pre"_(i-1)^-1 -= "pre"_i^-1 * a_i (mod p)`. 时间复杂度为 `O(n + log p)`.

试除法判断素数

143 = 11 * 13 不是素数

遍历所有因数

30 的所有因数为 1, 2, 3, 5, 6, 10, 15, 30

因数分解 (暴力算法) 设 n 是正整数, 返回 n 的所有素因子及其次数.

252 = 2^2 * 3^2 * 7

linux 系统可以用 factor 命令分解整数.

Pollard Rho因数分解, O(n^(1/4)), 1975 首先生成数列 `x_n`: 选取种子 `x_0`, `x_(n+1) = f(x_n) (mod n)`, 其中 `f(x) = x^2 + 1`. 设 `p | n`, 寻找整数 `x_i, x_j` 满足 `x_i = x_j (mod p)` 且 `x_i != x_j (mod n)`, 于是 `gcd(x_i - x_j, n)` 是 `n` 的非平凡因子.

65537 -> prime
9379 -> 113
32 -> boom!

Eratosthenes 筛法 求 n 以内的所有素数. 注意到素数的倍数一定是合数, 算法先将 bool 数组 flag 初始化为 false, 一旦发现 j 是合数, 就令 flag[j] = true. 时间复杂度为 `O(n log log n)`.

100 以内的素数为 2,3,5,7,11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97.

n `pi(10^n) = 10^n` 以内的素数个数
7 664,579
8 5,761,455
9 50,847,534

Euler 线性筛 求 N 以内素数的一种线性时间的算法. 在 Eratosthenes 算法中, 一个合数可能被反复判断: 如 12 被 2*6 排除, 也被 3*4 排除. 现在规定每个合数都被自己的最小素因子排除, 从而 12 只能由 2 来排除, 等等. 由于每个合数只被排除一次, 所以时间复杂度为 `O(N)`.
算法使用的数组 flag[N+1] 同上. 另有数组 primes[M] 用于保存已经得到的素数. 可以参考上表确定合适的 M 值, 如 N = 1e9 时可取 M = 50847534.

算法的关键一步: 若素数 `p | i`, 则对任意素数 `q gt p`, `q` 都不可能是 `k = i q` 的最小素因子, 所以算法内循环要在 i % p == 0 时 break.
现在证明这个算法确实使得每个合数都被自己的最小素因子排除. 由于每个合数都可以写为自己的最小素因子和一个正整数之积, 故只需证明对 `i = 2, 3, 4, cdots`, 任何形如 `k = i p` 的合数都被 `p` 排除, 其中 `p` 是 `k` 的最小素因子.
`i = 2` 时, 要使素数 `p` 是 `k = i p` 的最小素因子, 只有 `p = 2`. 显然算法将 `4 = 2 xx 2` 排除, 因此 `i = 2` 时结论成立.
`i gt 2` 时, 假设结论对所有小于 `i`, 大于等于 `2` 的正整数均成立, 由于算法的外层和内层循环分别枚举整数 `i` 和素数 `p`, 所以合数 `k = i p` 必然在某一时刻被排除. 现在设 `p` 是 `k = i p` 的最小素因子, 但 `k = j q` 却被 `q` 所排除, 其中 `q gt p` 是一素数. 设 `k_1 = j p`, 显然 `p` 也是 `k_1` 的最小素因子. 注意到 `j lt i`, 由归纳假设, `k_1` 应该由 `p` 排除. 然而由 `p | k = j q` 和 `p, q` 是不同素数知 `p | j`, 按算法, 内循环应当执行 break 语句. 换言之, 不可能使内层循环继续迭代, 从而排除合数 `k = j q`. 矛盾. 所以 `k` 必定由它的最小素因子排除.
由数学归纳法, 结论对所有 `i = 2, 3, 4, cdots` 成立.

已知正整数 `n`, 求最大的正整数 `m | n`, 且 `m` 无平方因子, 即不存在素数 `p` 使得 `p^2 | m`.

为避免对 `n` 作因子分解, 可以求出 `2` 到 `sqrt n` 间所有素数的乘积 `M`, 然后 `m = gcd(M, n)`. 这里 `M` 是一个很大的数; 如果 `n` 的范围已知, 可以通过查表得到 `M`. 如果 `M` 过大, 可以考虑每三个素数乘积求一次 gcd, 再将结果相乘, 即 `m = gcd(2*3*5, n) * gcd(7*11*13, n) * cdots`

CRC 算法 (Cyclic Redundancy Check, 循环冗余校验) 是用于通信的校验算法. 该算法的输入是一个 `n+1` 位二进制串 `b` 和 一个 `bbb F_2` 上的 `n` 次既约多项式 `f`.
所谓 `bbb F_2` 上的 `n` 次多项式是指, 系数为 0 或 1, 系数上的加法和乘法是模 2 的 `n` 次多项式. `bbb F_2` 上所有次数不超过 `n` 的多项式共有 `2^n` 个, 它们构成一个域 `bbb F_(2, n)(x)`, 可以进行四则运算. 如乘法为 `(x^3 + x^2 + 1)(x^3 + x + 1)` `= x^6 + x^5 + x^4 + 3x^3 + x^2 + x + 1` `= x^6 + x^5 + x^4 + x^3 + x^2 + x + 1`. 以上过程简记为 1101 × 1011 = 111 1111. 又, `bbb F_(2, n)(x)` 上的加减法恰好是异或运算: 1101 + 1011 = 1101 - 1011 = 0110. 所谓 `f` 是既约的, 是指 `f` 在 `bbb F_2` 上没有根, 这里就是说 0 和 1 都不是 `f` 的根. 这容易办到, 只要保证 `f` 有奇数个系数为 1 的项, 且常数项等于 1 即可.
当我们计算乘法所得的结果超过 `n` 次, 这个结果需要模 `f`, 使它落在 `bbb F_(2, n)(x)` 中. 比如, 当 `n = 4`, 取既约多项式 `f = x^4 + x + 1`, 则上面的结果化为 `x^6 + x^5 + x^4 + x^3 + x^2 + x + 1` `= (x^4+x+1)(x^2+x+1) + x^2 + x` 即, 111 1111 模 `f =` 1 0011 的最终结果是 `x^2 + x =` 110. 这个求余数的过程很重要, 我们用竖式再算一遍 (注意减法就是异或运算, 不需要借位): 
              111 --- 商 x^2 + x + 1 不重要, 可以不写
        ---------
1 0011 / 111 1111
         100 11
         -------
          11 001
          10 011
          -------
           1 0101
           1 0011
           ------
              110 --- 余数 x^2 + x

现在可以介绍 CRC 的计算过程了: `CRC(b, f) = x^n b(x) mod f(x) + x^n b(x)`. 将字节串 `b` 视作多项式, 首先将它乘以 `x^n`, 这相当于在串尾补 `n` 个 0. 然后计算余数 `r(x) = x^n b(x) mod f(x)`, 这是一个次数不超过 `n` 的多项式. 最后将 `x^n b(x)` 与 `r(x)` 相加, 即用 `r(x)` 替换串尾的 `n` 个 0 即可.
将数据 `b(x)` 和校验帧 `r(x)` 相连, 一起发送到接收方, 接收方只需验证 `CRC(b, f) mod f(x) -= 0`, 就通过了 CRC 校验. 若余数不为零, 则说明通信传输出现了问题.

数字 `n` 称为 CRC 算法的位宽 (width), 比如 CRC-8 的位宽等于 8, CRC-32 的位宽等于 32 等. 多项式 `f` 称为 CRC 算法的生成多项式. 如 CRC-8 的生成多项式是 `x^8 + x^2 + x + 1`, CRC-32 的生成多项式是 `x^32 + x^26 + x^23 + x^22 + x^16 + x^12 + x^11 + x^10 + x^8 + x^7 + x^5 + x^4 + x^2 + x + 1`. `f` 用二进制表示时, 由于首项系数必为 1, 只记录 `n` 个低次项即可. 如 CRC-8 和 CRC-32 的生成多项式分别记为 0x07 和 0x04C1 1DB7

下面这个计算 CRC-16 的例子来自 CSDN: 
#include <stdio.h>
#include <stdint.h>

const uint16_t POLY = 0x1021;

/* @param data 待校验的数据
 * @param len 数据包含的字节数
 * @param crc 初始的 CRC 余数
 */
uint16_t crc16(unsigned char *data, int len, uint16_t crc) {
  while (len--) {
    crc = crc ^ (*data++ << 8); // fetch byte from memory, XOR into CRC top byte
    // loop for 8 bits
    for (int i = 0; i < 8; i++) {
      if (crc & 0x8000) // if b15 is set
        crc = (crc << 1) ^ POLY;
      else // if b15 is off
        crc <<= 1;
    }
    crc &= 0xffff; // ensure CRC remains 16-bit value
  }
  return crc;
}

int main() {
  uint16_t c1, c2;

  // 计算 data 的 CRC-16
  unsigned char data[] = {'1', '2', '3', '4', '5', '6', '7', '8', '9'};
  c1 = crc16(data, 9, 0xffff);
  printf("%04x\n", c1); // 29b1

  // 分两段计算 CRC-16
  c2 = crc16(data, 4, 0xffff);
  c2 = crc16(data + 4, 5, c2);
  printf("%04x\n", c2); // 29b1

  return 0;
}

这里是 CRC 在线计算.

求模 `p` 的平方根, `p -= 3(mod 4)` 设 `p` 是模 4 余 3 的素数, `n` 为模 `p` 的二次剩余, 记 `s = (p-1)/2`, 由 Legendre 符号知道 `n^s -= 1 (mod p)`. 这里 `s` 是奇数, 容易验证 `+-n^((s+1)//2)` 就是 `n` 模 `p` 的平方根.

对于模 4 余 1 的素数, 我们有如下算法:

    Cipolla 算法: 求模 `p` 的平方根
    [来自 olderciyuan] 设 `p` 为奇素数, `n` 为模 `p` 的二次剩余, 则 `n` 模 `p` 恰有两个平方根. Cipolla 算法能在 `O(log p)` 时间内求出这两个平方根.
  1. 取合适的 `a in bbb F_p`, 使得 `a^2 - n` 是二次非剩余. 由于模 `p` 的二次非剩余恰有一半, 我们只要随机选取, 然后用 Euler 判别法, 若 `(a^2-n)^((p-1)//2) = -1`, 即可判断它是二次非剩余.
  2. 将 `a^2-n` 的一个平方根 `i = sqrt(a^2 - n)` 加到 `bbb F_p` 中, 得到扩域 `bbb F_p(i)`. 这个域的特征为 `p`, 其中的元素形如 `A + Bi`, `A, B in bbb F_p`.
  3. 可以证明 `(a+i)^(p+1) = n`, 因此 `n` 在扩域 `bbb F_p(i)` 中有一对平方根 `+-(a+i)^((p+1)//2)`. 事实上, 这对平方根属于 `bbb F_p`, 因此它们就是 `n` 模 `p` 的平方根.
  4. 上述过程中应用快速幂算法, 总的时间复杂度为 `O(log p)`.
  1. 先证 `i^p = -i`. 事实上 `i * i^(p-1)` `= i * (a^2-n)^((p-1)//2)` `= -i`. 最后一个等式成立是因为 `a^2-n` 为二次非剩余.
  2. Fermat 小定理: `AA a in bbb F_p`, `a^p = a`.
  3. 由于 `bbb F_p(i)` 的特征为 `p`, 故 `AA a, b in bbb F_p(i)` 有 `(a+b)^p = a^p + b^p`.
  4. 现在可以算出: `(a+i)^(p+1)` `= (a+i)(a+i)^p` `= (a+i)(a^p+i^p)` `= (a+i)(a-i)` `= a^2-i^2` `= n`.
  5. 下证 `n` 的平方根必定属于 `bbb F_p`. 设 `(A+B i)^2 = n`, 即 `A^2 + B i^2 + 2A B i = n`. 比较两边虚部得 `A B = 0`, 因此 `A = 0` 或 `B = 0`.
    若 `A = 0`, 有 `B^2 i^2 = n`, 但 `B^2, n` 均为二次剩余, `i^2` 为二次非剩余, 矛盾. 故 `B = 0`, 即 `A + B i in bbb F_p`.

说明: `5^2 -= 8 (mod 17)`, 因此 8 模 17 的平方根为 `+-5`. 我们只显示其中一个根 (最小正根).

用 sympy 求模 `p` 的平方根: sympy.ntheory.residue_ntheory.sqrt_mod(a, p, all_roots=False)

中国剩余定理 给定两两互素的正整数 `n_1, cdots, n_k` 和任意 `k` 个整数 `a_1, cdots, a_k`, 则线性同余方程组 `{ x -= a_1 (mod n_1); cdots; x -= a_k (mod n_k); :}` 在模 `N` 的意义下存在唯一解 `x`. 这里 `N = n_1 * cdots * n_k`.
事实上, 记 `M_i = prod_(j != i) n_j`, 又记 `M_i^-1` 是 `M_i` 模 `n_i` 的逆, 则 `x -= sum_(i=1)^k a_i * M_i^-1 * M_i (mod N)`. 验证: 上式两边模 `n_i`, 由于除了 `M_i` 这一项外, 各项都是 `n_i` 的倍数, 我们得到 `x -= a_i * M_i^-1 * M_i -= a_i (mod n_i)`.

说明: `x -= 2 (mod 5), x -= 3 (mod 11), x -= 5 (mod 17)` 的解为 `x -= 872 (mod 935).