Math Algorithm

常见的数学算法包括：快速幂、矩阵运算、欧拉函数与欧拉定理、欧几里得算法及其扩展、中国剩余定理（CRT）、高斯消元、乘法逆元、组合数求解、质数筛选、博弈论模型及其他杂项技巧。

快速幂

快速幂的目的就是做到快速求幂，假设我们要求 a^b,按照朴素算法就是把 a 连乘 b 次，这样一来时间复杂度是O(b)也即是O(n)级别，快速幂能做到O(log n).

例如当 b=11 时, 将 b 写成二进制格式 1011，即1 * 2³ + 0 * 2² + 1 * 2¹ + 1 * 2⁰：

a¹¹ = a^{(2⁰ + 2¹ + 2³)}

int qpow(int a, int b) {
int ans = 1, base = a;
  while (b != 0) {
      if ((b & 1) != 0) ans *= base;
      base *= base;
      b >>= 1;
  }
  return ans;
}

快速幂模块

另附快速乘法 (mod p)意义下

long qmul(long a,long b,long mod) {
  long res=0;
   while(b>0) {
       if((b&1)==1) res=(res+a)%mod;
       a=(a+a)%mod;
       b>>=1;
   }
   return res;
}

矩阵

矩阵快速幂矩阵快速幂就是把快速幂里的乘法换成矩阵乘法 求解 A = mat^K mat为原矩阵，A是答案

static long n, N, k, P = (long)1e9+7;
static long[][] mat, A;
private static void qmpow() {
    A = new long[N][N];
    for (int i = 0;i < n;i++) A[i][i] = 1; //A首先构造成一个单位矩阵（相当于 1）
    while (k != 0){
        if ((k & 1) != 0) A = mul(A, mat); //把快速幂中的乘法换成矩阵乘法就行
        mat = mul(mat, mat);
        k >>= 1;
    }
}
private static long[][] mul(long[][] x, long[][] y) {
    long[][] res = new long[x.length][y[0].length];
    for (int i = 0;i < x.length;i++) for (int j = 0;j < y[0].length;j++){
    for (int k = 0;k < y.length;k++)
        res[i][j] = (res[i][j]+x[i][k]*y[k][j]%P)%P;
    }
    return res;
}

矩阵加速（数列）板子

那么我们可以维护两个我们首先要确定目标矩阵。下面这个矩阵就是我想要的矩阵.

并且有递推式 F[i] = F[i-1]+F[i-3] 那么这个矩阵每个元素就可以表示成如下

通过每一项的系数可以得出M_i与M_i − 1递推关系

$记$

则写个两三项可发现 F[i] = A^i − 3[0][0] + A^i − 3[0][1] + A^i − 3[0][2]

用矩阵快速幂求出A^i − 3即可

欧拉函数&欧拉定理

欧拉函数 Euler’s totient function，即 φ(n) ，表示的是小于等于 n 和 n 互质的数的个数。其中 p₁, p₂, p₃, ⋯, p_k 是 n 的所有不重复的素因子。对应下方代码中 ans = ans/i * (i − 1);

int euler_phi(int n) {
     // 分解质因数过程中计算
   int ans = n;
   for (int i = 2; i * i <= n; i++)
        if (n % i == 0) {
            ans = ans / i * (i - 1);
            while (n % i == 0) n /= i;
        }
   if (n > 1) ans = ans / n * (n - 1);
   return ans;
}

欧拉定理若 gcd (a, m) = 1 ，则 a^φ(m) ≡ 1( mod m) 。
拓展欧拉定理

筛法求欧拉函数注意到在线性笑中，每一个合数都是被最小的质因子笑掉。比如设 p₁ 是 n 的最小质因子，，那么线性篮的过程中 n 通过 n^′ × p₁ 篮掉。观察线性篮的过程，我们还需要处理两个部分，下面对 n^′ mod p₁ 分情况讨论。如果 n^′ mod p₁ = 0 ，那么 n^′ 包含了 n 的所有质因子。

那如果 n^′ mod p₁ ≠ 0 呢，这时 n^′ 和 p₁ 是互质的，根据欧拉函数性质，我们有:

static void euler(int n){
    phi[1] = 1;
    for (int i = 2;i <= n;i++){
        if (is[i] == 0){
            prime[++cnt] = i;
            phi[i] = i-1;
        }
        for (int j = 1;i*prime[j] <= n;j++){
            is[i*prime[j]] = 1;
            if (i % prime[j] != 0){
                phi[i*prime[j]] = phi[i]*(prime[j]-1);
            }else{
                phi[i*prime[j]] = phi[i]*prime[j];
                break;
            }
        }
    }
}

欧几里得定理

辗转相除法

最大公因数 GCD(a, b)
最小公倍数 LCM(a, b) = a/GCD(a, b)*b

1
2
3

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

1
2
3

int lcm(int a, int b){
	return a/gcd(a, b)*b;
}

Binary GCD

a, b 均为偶数。对于这种情况，显然 2 是公约数之一，直接将两个数都除以 2 ，再递归下去求即可，
a, b 中有且仅有一个偶数。不妨设 a 是偶数，那么显然 2 不是公约数之一，直接将 a 除以 2 求即可,
a, b 均为奇数。不妨设 a > b ，那么有。设 a = k₁gcd (a, b), b = k₂gcd (a, b) 且 gcd (k₁, k₂) = 1 ，先考虑证明 gcd (a − b, b) = gcd (a, b) a − b = (k₁ − k₂)gcd (a, b), b = k₂gcd (a, b) ，那么有gcd (a − b, b) = gcd (k₁ − k₂, k₂) ⋅ gcd (a, b) ，即证 gcd (k₁ − k₂, k₂) = 1 考虑反证法，假设 gcd (k₁ − k₂, k₂) = m 其中 m ∈ ℕ^*, m > 1 ，再设 k₁ − k₂ = t₁m, k₂ = t₂m 就有 k₁ = (t₁ + t₂)m ，那么 gcd (k₁, k₂) = m ⋅ gcd (t₁ + t₂, t₂) > 1 与 k₁, k₂ 互质矛盾，得证 gcd (k₁ − k₂, k₂) = 1

int bgcd(int a, int b){
  if (a == 0) return b;
  if (b == 0) return a;
  int p = 0;
  while (((a|b)&1) == 0){
      a >>= 1; b >>= 1; p++;
  }
  if (a > b) return bgcd(a-b, b)<<p;
  else return bgcd(a, b-a)<<p;
}

拓展欧几里得

贝祖定理：a、b 是整数，那么一定存在整数 x、y 使得 ax+by=gcd(a,b)。

如果 ax+by=1 有解，那么gcd(a,b)=1；

解方程：ax+by=m (a, b, m 均为正整数)

首先解 ax+by=gcd(a, b). 递归边界：b = 0, a = gcd(a, b); 显然 x = 1, y = 0. 递归推导：当前层为 a, b; 求解 ax+by=gcd(a, b). 由递归下一层 b, a%b, 可得 bx₁+(a%b)y₁=gcd(a, b).

由于 a%b = a - (a/b)b. —(3) 联立(2)(3) 解得 ay₁ + b(x₁ - (a/b)y₁) = gcd(a, b). 对比 (1) 式解得 x = y₁, y = x₁-a/b*y₁.

static int x, y;
static int exgcd(int a, int b){
    if (b == 0){
        x = 1; y = 0;
        return a;
    }
    int gcd = exgcd(b, a%b);
    int tp = x;
    x = y;
    y = tp-a/b*y;
    return gcd;
}

解 ax+by=m 若方程 ax+by=m 的一组整数解为(x，y). 则它的任意整数解都可以写成(x + k _ b’，y - k _ a’)。其中 a’ = a / gcd(a，b), b’ = b / gcd(a，b)，k 取任意整数

ax+by=m 有解的充分条件 ==gcd(a，b) | m== //m 能整除gcd(a, b)

设（x₁， y₁）是 ax+by=gcd(a, b) 解

则 ax+by=m 的解为 x = x₁c/gcd(a，b), y = y₁c/gcd(a，b),

public static void main(String[] args){
    //a*x+b*y=m
    int d = exgcd(a, b, x, y);
    if (m % d != 0) out.println("无解");
    else{
        x = m/d*x; // x为一个特解
        y = m/d*y; // y为一个特解
        int kx = b/gcd, ky = a/gcd; //即为上面的 b' 和 a'
        int minx = (x%kx+kx)%kx; // 求此时x的最小非负整数解
        int miny = (y%ky+ky)%ky; // 求此时y的最小非负整数解
        int maxx = (c-miny*b)/a; // 求y为最小非负整数，此时x的最大整数解
        int maxy = (c-minx*a)/b; // 求x为最小非负整数，此时y的最大整数解
    }
}

青蛙约会

中国剩余定理（CRT）

求解如下形式的一元线性同余方程组==最小正整数解==。(注意b₁,b₂,⋯,b_k两两不一定互质，拓展版)

设两个方程分别是 x ≡ a₁( mod m₁)、x ≡ a₂( mod m₂) ；将它们转化为不定方程: x = m₁p + a₁ = m₂q + a₂ ，其中 p, q 是整数，则有 m₁p − m₂q = a₂ − a₁ 由裴蜀定理，当 a₂ − a₁ 不能被 gcd (m₁, m₂) 整除时，无解；其他情况下，可以通过扩展欧几里得算法解出来一组可行解 (p, q) ；则原来的两方程组成的模方程组的解为 x ≡ b( mod M) ，其中 b = m₁p + a₁，M = lcm (m₁, m₂)

static void solve() throws Exception {
    //  x = py+r => x==r(mod p)
    n = ni();
    long p1 = ni(), r1 = ni();
    for (int i = 2;i <= n;i++){
        long p2 = ni(), r2 = ni();
        long gcd = exgcd(p1, p2);
        if (abs(r2-r1)%gcd != 0){
            out.println(-1);
            return;
        }
        x = x*(r2-r1)/gcd;
        long mod = Math.abs(p2)/gcd;
        x = (x%mod+mod)%mod;
        r1 = x*p1+r1;
        p1 = p1/gcd*p2;
    }
    out.println(r1);
}

高斯消元（解一次方程组）

1.选择一个尚未被选过的未知数作为主元，选择一个包含这个主元的方程。

2.通过加减消元，消掉其它方程中的这个未知数。

3.重复以上步骤，直到把每一行都变成只有一项有系数。（转换为对角矩阵）

4.输出时记得除以方程主元的系数。

static double[][] a = new double[N][N]; //double 存
static void solve() throws Exception {
    n = ni();
    for (int i = 1;i <= n;i++) // n行
        for (int j = 1;j <= n+1;j++) // n+1列
            a[i][j] = ni();
    for (int i = 1;i <= n;i++){
        int max = i;
        for (int j = i+1;j <= n;j++) //找到每一列中 最大的那个数
            if (abs(a[j][i]) > abs(a[max][i])) max = j;
        for (int j = 1;j <= n+1;j++){//交换
            double tp = a[i][j];
            a[i][j] = a[max][j];
            a[max][j] = tp;
        }
        if (a[i][i] == 0){
            out.println("不存在唯一解");
            return;
        }
        for (int j = 1;j <= n;j++){// 消除系数
            if (i == j) continue;
            double tp = a[j][i]/a[i][i];
            for (int k = i+1;k <= n+1;k++)
                a[j][k] -= a[i][k]*tp;
        }
    }
    for (int i = 1;i <= n;i++)
        //最终每一个未知数解要除以它的系数（a[i][i]）
        out.printf("%.2f\n", a[i][n+1]/a[i][i]);
}

乘法逆元

定义：如果一个线性同余方程 ax ≡ 1 ( p) ，则 x 称为 a p 的逆元，记作。所以我们可以知道：a 除以一个数模 p，等于 a 乘这个数的乘法逆元模 p.

快速幂 费马小定理（==此时 p 一定为素数==）

int qpow(long a, long b) {
    int ans = 1;
    a = (a % p + p) % p; //让a变成正数（模p的意义下）
    while (b != 0) {
        if (b & 1) ans = (a * ans) % p;
        a = (a * a) % p;
        b >>= 1;
    }
    return ans;
}

线性算法求一连串数字的逆元（==此时 p 一定是素数==）

用于求一连串数字对于一个 mod p的逆元。洛谷 P3811

只能用这种方法，别的算法都比这些要求一串要慢。

首先我们有一个, 1⁻¹ ≡ 1 (mod p)

然后设 p = k * i + r, (1 < r < i < p)也就是 k是 p/i的商，r 是余数。

再将这个式子放到 (mod p)意义下就会得到：

k * i + r ≡ 0 (mod p)

然后乘上i⁻¹, r⁻¹就可以得到:

于是，我们就可以从前面推出当前的逆元了。

1
2
3

inv[1] = 1;
for(int i = 2; i <= n; ++ i)
    inv[i] = (p - p / i) * inv[p % i] % p;

组合数(模P意义下)

打表（杨辉三角）

for(int i=0;i<=n;i++){
    c[i][0]=c[i][i]=1;
    for(int j=1;j<i;j++){
        c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
    }
}

预处理 O(nm) 查询 O(1)

乘法逆元+快速幂+阶乘

int pow(int x,int y){...} //快速幂
//求逆元 这里用欧拉定理，只能适用p为素数;可以用拓欧求逆元。
int inv(int x,int p){
    return pow(x,p-2);
}
fac[0]=1;
for(int i=1;i<=n;i++)
    fac[i]=fac[i-1]*i; //预处理出阶乘
return ((fac[n]*inv(fac[m],p))%p*inv(fac[n-m],p))%p;

预处理 O(n) 查询 O(log p)

分解质因数

static void solve() throws Exception{
    long t = nl();
    // cnt[i][0]:第i个质数 cnt[i][1]:第i个质因数出现的次数
    int cnt[][] = new int[(int)1e6][2], len = 0;
    for (int i = 2;(long)i*i <= t;i++){
        if (t % i == 0){
            cnt[len][0] = i;
            while (t % i == 0){
                cnt[len][1]++;
                t /= i;
            }len++;
        }
    }
}

质数

埃氏筛 O(Nlog₂log₂N)

bool numlist[100000001]; //true:合数 false:素数
int prime[20000001]; //存储从小到大的素数（2，3，5，···）
int cnt; //存储素数数量
void work(int n){
    for(int i=2; i<=n; i++){
        if(numlist[i]==false){
            prime[++cnt] = i ;
            for(int j=i; i*j<=n; j++)
                numlist[i*j] = true;
        }
    }
    return;
}

欧拉筛 O(N)

bool numlist[100000001];
int prime[20000001], cnt;
void work(int n){
    for(int i=2; i<=n; i++){
        if(numlist[i]==false) prime[++cnt] = i ;
        for(int j=1; j<=cnt&&i*prime[j]<=n; j++){
            numlist[i*prime[j]] = true ;
            if(i%prime[j]==0)
                break;
        }
    }
    return;
}

博弈论

nim 游戏

nim 游戏的规则是这样的：地上有 n 堆石子，每人每次可从任意一堆石子里取出任意多枚石子扔掉，可以取完，不能不取。每次只能从一堆里取。最后没石子可取的人就输了。假如甲是先手，且告诉你这 n 堆石子的数量，他想知道是否存在先手必胜的策略。

1
2
3

int t = 0;
while (n-- > 0) t ^= ni();
out.println((t==0?"No":"Yes"));

杂项

三角形面积

海伦公式：假设在平面内，有一个三角形，边长分别为 a、b、c ，三角形的面积 S 可由以下公式求得:

而公式里的 p 为半周长 (周长的一半) :

向量积：

任意多边形面积

对于任意一个多边形，如果已知其各个顶点的坐标 A₁(x₁, y₁), A₂(x₂, y₂), …, A_n(x_n, y_n) ，那么这个多边形的面积为:

其中 x_n + 1 = x₁, y_n + 1 = y₁。