算法模板 数学

矩阵

矩阵快速幂

测试：洛谷P3390, P1939
struct mat{
    int a[sz][sz];//注意开long long
    mat(){memset(a, 0, sizeof(a));}
    mat operator - (const mat &b)const{
        mat res;
        for(int i = 0; i < sz; i++)
            for(int j = 0; j < sz; j++)
                res.a[i][j] = (a[i][j] - b.a[i][j] + mod) % mod;
        return res;
    }
    mat operator + (const mat &b)const{
        mat res;
        for(int i = 0; i < sz; i++)
            for(int j = 0; j < sz; j++)
                res.a[i][j] = (a[i][j] + b.a[i][j]) % mod;
        return res;
    }
    mat operator * (const mat &b)const{
        mat res;
        for(int i = 0; i < sz ; i++)
            for(int j = 0; j < sz; j++)
                for(int k = 0; k < sz; k++)
                    res.a[i][j] = (res.a[i][j] + a[i][k] * b.a[k][j] % mod) % mod;
        return res;
    }
    mat operator ^ (int x)const{
        mat res, base;
        for(int i = 0; i < sz; i++) res.a[i][i] = 1;
        for(int i = 0; i < sz; i++)
            for(int j = 0; j < sz; j++)
                base.a[i][j] = a[i][j] % mod;
        while(x){
            if(x & 1) res = res * base;
            base = base * base;
            x >>= 1;
        }
        return res;
    }
};

矩阵树定理

/*
洛谷P6178
n个点，m条边，生成树T的权值为T所有边权的乘积
生成树为根为1的外向树
0为无向图，1为有向图
*/
#include<bits/stdc++.h>

using namespace std;
typedef long long ll;
const int N = 310;
const int mod = 1e9 +7;
int a[N][N];
ll ans;

int qpow(int a, int b){
    int res = 1;
    while(b){
        if(b & 1) res = (ll)res * a % mod;
        a = (ll)a * a % mod;
        b >>= 1;
    }
    return res;
}

int gauss(int n, int m){
    int r, c;
    for(r = 2, c = 2; r <= n && c <= m; r++, c++){
        int t = r;
        for(int i = r; i <= n; i++){
            if(a[i][c]){
                if(i != t) ans = mod - ans;
                t = i;
                break;
            }
        }
        if(t != r){
            for(int i = c; i <= m; i++)
                swap(a[r][i], a[t][i]);
        }
        int inv = qpow(a[r][c], mod - 2);
        for(int i = r + 1; i <= n; i++){
            if(a[i][c] == 0) continue;
            int x = (ll)a[i][c] * inv % mod;
            for(int j = c; j <= m; j++)
                a[i][j] = (a[i][j] - (ll)a[r][j] * x % mod + mod) % mod;
        }
    }
}

int main(){
    int n, m, t;
    cin >> n >> m >> t;
    while(m--){
        int u, v, w;
        cin >> u >> v >> w;
        if(t == 0){//无向图
            a[u][u] = (a[u][u] + w) % mod;
            a[v][v] = (a[v][v] + w) % mod;
            a[u][v] = (a[u][v] - w + mod) % mod;
            a[v][u] = (a[v][u] - w + mod) % mod;
        }
        else{//有向图
            a[v][v] = (a[v][v] + w) % mod;
            a[u][v] = (a[u][v] - w + mod) % mod;
        }
    }
    ans = 1;
    gauss(n, n);
    for(int i = 2; i <= n; i++)
        ans = ans * a[i][i] % mod;
    cout << (ans + mod) % mod << endl;
    return 0;
}

随机素数测试和大数分解

//测试洛谷P4718
struct Pollard_Rho{
    const static int N = 1e6 + 10;
    const static int S = 8;
    ll p[N], vis[N], fac[N], sz, tot; //多组输入注意初始化tot = 0
    ll mul(ll a, ll b, ll mod) {
        return (__int128) a * b % mod;
    }
    //如果代码超时且需要分解大数的质因数可以用，否则不要用
    void getprime(int n){
        int cnt = 0;
        sz = n;
        vis[1] = 1;
        for(int i = 2; i <= n; i++){
            if(!vis[i]){
                p[++cnt] = i;
                vis[i] = i;
            }
            for(int j = 1; j <= cnt && p[j] <= n / i; j++){
                int t = i * p[j];
                vis[t] = p[j];
                if(i % p[j] == 0) break;
            }
        }
    }
    ll qpow(ll a, ll b, ll mod) {
        ll res = 1;
        while(b){
            if(b & 1) res = mul(res, a, mod);
            a = mul(a, a, mod);
            b >>= 1;
        }
        return res;
    }
    bool check(ll a, ll n) {
        ll t = 0, u = n - 1;
        while (!(u & 1))
            t++, u >>= 1;
        ll x = qpow(a, u, n), xx = 0;
        while (t--) {
            xx = mul(x, x, n);
            if (xx == 1 && x != 1 && x != n - 1)
                return false;
            x = xx;
        }
        return xx == 1;
    }
    bool miller(ll n) {
        if(n == 2) return true;
        if(n < 2 || !(n & 1)) return false;
        if(n <= sz) return vis[n] == n;
        for(int i = 0; i <= S; ++i) {
            if(!check(rand() % (n - 1) + 1, n))
                return false;
        }
        return true;
    }
    ll gcd(ll a, ll b){
        return b == 0 ? a : gcd(b, a % b);
    }
    ll Abs(ll x) {
        return x < 0 ? -x : x;
    }
    ll Pollard_rho(ll n){
        ll s = 0, t = 0, c = rand() % (n - 1) + 1, v = 1, ed = 1;
        while(1) {
            for(int i = 1; i <= ed; ++i) {
                t = (mul(t, t, n) + c) % n;
                v = mul(v, Abs(t - s), n);
                if(i % 127 == 0) {
                    ll d = gcd(v, n);
                    if(d > 1) return d;
                }
            }
            ll d = gcd(v, n);
            if(d > 1) return d;
            s = t;
            v = 1;
            ed <<= 1;
        }
    }
    void getfactor(ll n){//得到所有的质因子(可能有重复的
        if(n <= sz) {
            while(n != 1)
                fac[tot++] = vis[n], n /= vis[n];
            return;
        }
        if(miller(n)) {
            fac[tot++] = n;
        }
        else{
            ll d = n;
            while(d >= n)
                d = Pollard_rho(n);
            getfactor(d);
            getfactor(n / d);
        }
    }
};

基于值域预处理的快速gcd

//基于值域预处理的快速gcd
//预处理O(值域) 查询O(1)
//测试：洛谷P5435
const int N = 1e6, T = 1e3;
int f[N + 10][3], pre[T + 10][T + 10], p[N + 10], cnt;
bool vis[N + 10];

void init(){
	f[1][0] = f[1][1] = f[1][2] = 1;
	for(int i = 2; i <= N; i++){
		if(!vis[i]){
			f[i][0] = f[i][1] = 1;
			f[i][2] = i;
			p[++cnt] = i;
		}
		for(int j = 1; j <= cnt && p[j] <= N / i; j++){
			int t = i * p[j];
			vis[t] = true;
			f[t][0] = f[i][0] * p[j];
			f[t][1] = f[i][1];
			f[t][2] = f[i][2];
			if(f[t][0] > f[t][1]) swap(f[t][0], f[t][1]);
			if(f[t][1] > f[t][2]) swap(f[t][1], f[t][2]);
			if(i % p[j] == 0) break;
		}
	}
	for(int i = 0; i <= T; i++)
		pre[0][i] = pre[i][0] = i;
	for(int i = 1; i <= T; i++)
		for(int j = 1; j <= i; j++)
			pre[i][j] = pre[j][i] = pre[j][i % j];
}

int gcd(int a, int b){
	int res = 1;
	for(int i = 0; i < 3; i++){
		int t;
		if(f[a][i] > T){
			if(b % f[a][i] == 0) t = f[a][i];
			else t = 1;
		}
		else t = pre[f[a][i]][b % f[a][i]];
		b /= t;
		res *= t;
	}
	return res;
}

二次剩余

struct Num{
    ll x, y;
    Num(){}
    Num(int x, int y):x(x), y(y){}
};
ll w;
Num mul(Num a, Num b, ll p){//复数乘法
    Num res(0, 0);
    res.x = ((a.x * b.x % p + a.y * b.y % p * w % p) % p + p) % p;
    res.y = ((a.x * b.y % p + a.y * b.x % p) % p + p) % p;
    return res;
}

ll qpow(ll a, ll b, ll p){
    ll res = 1;
    while(b){
        if(b & 1) res = res * a % p;
        a = a * a % p;
        b >>= 1;
    }
    return res;
}

ll qpow(Num a, ll b, ll p){
    Num res(1, 0);
    while(b){
        if(b & 1) res = mul(res, a, p);
        a = mul(a, a, p);
        b >>= 1;
    }
    return res.x % p;
}
//p必须是素数
//返回一个解x，另一个解为p - x
//测试：洛谷P5491
ll Cipolla(ll n, ll p){
    n %= p;
    if(p == 2) return n;
    if(qpow(n, (p - 1) / 2, p) == p - 1) return -1;//无解
    ll a;
    while(true){
        a = rand() % p;
        w = ((a * a % p - n) % p + p) % p;
        if(qpow(w, (p - 1) / 2, p) == p - 1) break;
    }
    Num x(a, 1);
    return qpow(x, (p + 1) / 2, p);
}