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

快速幂 (反复平方法) 求 `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, b 为正整数, 求最大公约数 d, 同时求出系数 x, y 满足 a x + b y = d.

求模 `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)`.

试除法判断素数

遍历所有因数

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

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` 的非平凡因子.

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

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`

求模 `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`.

用 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)`.