算法模板の数学&数论

1、求逆元

int inv(int a) {
    if(a == 1) return 1;
    return (MOD - MOD / a) * inv(MOD % a);
}

2、线性筛法

bool isPrime[MAXN];
int label[MAXN], prime[MAXN];
int n, total;

void makePrime() {
    n = 100000;
    for(int i = 2; i <= n; ++i) {
        if(!label[i]) {
            prime[total++] = i;
            label[i] = total;
        }
        for(int j = 0; j < label[i]; ++j) {
            if(i * prime[j] > n) break;
            label[i * prime[j]] = j + 1;
        }
    }
    for(int i = 0; i < total; ++i) isPrime[prime[i]] = true;
}

3、欧拉函数

void phi_table(int n) {
    //for(int i = 2; i <= n; ++i) phi[i] = 0;
    phi[1] = 1;
    for(int i = 2; i <= n; ++i) if(!phi[i])
        for(int j = i; j <= n; j += i) {
            if(!phi[j]) phi[j] = j;
            phi[j] = phi[j] / i * (i - 1);
        }
}

4、高精度类模板（减法除法的正确性未验证）

#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
using namespace std;

const int MAXN = 100010;

struct bign {
    int len, s[MAXN];

    bign () {
        memset(s, 0, sizeof(s));
        len = 1;
    }
    bign (int num) { *this = num; }
    bign (const char *num) { *this = num; }

    void clear() {
        memset(s, 0, sizeof(s));
        len = 1;
    }

    bign operator = (const int num) {//数字
        char s[MAXN];
        sprintf(s, "%d", num);
        *this = s;
        return *this;
    }
    bign operator = (const char *num) {//字符串
        for(int i = 0; num[i] == '0'; num++) ;  //去前导0
        if(*num == 0) --num;
        len = strlen(num);
        for(int i = 0; i < len; ++i) s[i] = num[len-i-1] - '0';
        return *this;
    }

    bign operator + (const bign &b) const {
        bign c;
        c.len = 0;
        for(int i = 0, g = 0; g || i < max(len, b.len); ++i) {
            int x = g;
            if(i < len) x += s[i];
            if(i < b.len) x += b.s[i];
            c.s[c.len++] = x % 10;
            g = x / 10;
        }
        return c;
    }

    bign operator += (const bign &b) {
        *this = *this + b;
        return *this;
    }

    void clean() {
        while(len > 1 && !s[len-1]) len--;
    }

    bign operator * (const bign &b) {
        bign c;
        c.len = len + b.len;
        for(int i = 0; i < len; ++i) {
            for(int j = 0; j < b.len; ++j) {
                c.s[i+j] += s[i] * b.s[j];
            }
        }
        for(int i = 0; i < c.len; ++i) {
            c.s[i+1] += c.s[i]/10;
            c.s[i] %= 10;
        }
        c.clean();
        return c;
    }
    bign operator *= (const bign &b) {
        *this = *this * b;
        return *this;
    }

    bign operator *= (const int &b) {//使用前要保证>len的位置都是空的
        for(int i = 0; i < len; ++i) s[i] *= b;
        for(int i = 0; i < len; ++i) {
            s[i + 1] += s[i] / 10;
            s[i] %= 10;
        }
        while(s[len]) {
            s[len + 1] += s[len] / 10;
            s[len] %= 10;
            ++len;
        }
        return *this;
    }

    bign operator - (const bign &b) {
        bign c;
        c.len = 0;
        for(int i = 0, g = 0; i < len; ++i) {
            int x = s[i] - g;
            if(i < b.len) x -= b.s[i];
            if(x >= 0) g = 0;
            else {
                g = 1;
                x += 10;
            }
            c.s[c.len++] = x;
        }
        c.clean();
        return c;
    }
    bign operator -= (const bign &b) {
        *this = *this - b;
        return *this;
    }

    bign operator / (const bign &b) {
        bign c, f = 0;
        for(int i = len - 1; i >= 0; i--) {
            f *= 10;
            f.s[0] = s[i];
            while(f >= b) {
                f -= b;
                c.s[i]++;
            }
        }
        c.len = len;
        c.clean();
        return c;
    }
    bign operator /= (const bign &b) {
        *this  = *this / b;
        return *this;
    }

    bign operator % (const bign &b) {
        bign r = *this / b;
        r = *this - r*b;
        return r;
    }
    bign operator %= (const bign &b) {
        *this = *this % b;
        return *this;
    }

    bool operator < (const bign &b) {
        if(len != b.len) return len < b.len;
        for(int i = len-1; i >= 0; i--) {
            if(s[i] != b.s[i]) return s[i] < b.s[i];
        }
        return false;
    }

    bool operator > (const bign &b) {
        if(len != b.len) return len > b.len;
        for(int i = len-1; i >= 0; i--) {
            if(s[i] != b.s[i]) return s[i] > b.s[i];
        }
        return false;
    }

    bool operator == (const bign &b) {
        return !(*this > b) && !(*this < b);
    }

    bool operator != (const bign &b) {
        return !(*this == b);
    }

    bool operator <= (const bign &b) {
        return *this < b || *this == b;
    }

    bool operator >= (const bign &b) {
        return *this > b || *this == b;
    }

    string str() const {
        string res = "";
        for(int i = 0; i < len; ++i) res = char(s[i]+'0') + res;
        return res;
    }
};

istream& operator >> (istream &in, bign &x) {
    string s;
    in >> s;
    x = s.c_str();
    return in;
}

ostream& operator << (ostream &out, const bign &x) {
    out << x.str();
    return out;
}

5、单纯形（UVA 10498）//只能解标准形式

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;

const double EPS = 1e-10;
const int MAXN = 55;
const int INF = 0x3fff3fff;

inline int sgn(double x) {
    return (x > EPS) - (x < -EPS);
}

double A[MAXN][MAXN];
double b[MAXN], c[MAXN];
int N[MAXN], B[MAXN];
int n, m;
double v;

bool init() {
    N[0] = B[0] = v = 0;
    for(int i = 1; i <= n; ++i) N[++N[0]] = i;
    for(int i = 1; i <= m; ++i) B[++B[0]] = n + i;
    return true;
}

void pivot(int l, int e) {
    b[e] = b[l] / A[l][e];
    A[e][l] = 1.0 / A[l][e];
    for(int i = 1; i <= N[0]; ++i) {
        int &x = N[i];
        if(x != e) A[e][x] = A[l][x] / A[l][e];
    }
    for(int i = 1; i <= B[0]; ++i) {
        int &y = B[i];
        b[y] -= A[y][e] * b[e];
        A[y][l] = -A[y][e] * A[e][l];
        for(int j = 1; j <= N[0]; ++j) {
            int &x = N[j];
            if(x != e) A[y][x] -= A[e][x] * A[y][e];
        }
    }
    v += b[e] * c[e];
    c[l] = -A[e][l] * c[e];
    for(int i = 1; i <= N[0]; ++i) {
        int &x = N[i];
        if(x != e) c[x] -= A[e][x] * c[e];
    }
    for(int i = 1; i <= N[0]; ++i) if(N[i] == e) N[i] = l;
    for(int i = 1; i <= B[0]; ++i) if(B[i] == l) B[i] = e;
}

bool simplex() {
    while(true) {
        int e = MAXN;
        for(int i = 1; i <= N[0]; ++i) {
            int &x = N[i];
            if(sgn(c[x]) > 0 && x < e) e = x;
        }
        if(e == MAXN) break;
        double delta = -1;
        int l = MAXN;
        for(int i = 1; i <= B[0]; ++i) {
            int &y = B[i];
            if(sgn(A[y][e]) > 0) {
                double tmp = b[y] / A[y][e];
                if(delta == -1 || sgn(tmp - delta) < 0 || (sgn(tmp - delta) == 0 && y < l)) {
                    delta = tmp;
                    l = y;
                }
            }
        }
        if(l == MAXN) return false;
        pivot(l, e);
    }
    return true;
}

int main() {
    while(scanf("%d%d", &n, &m) != EOF) {
        for(int i = 1; i <= n; ++i) scanf("%lf", &c[i]);
        for(int i = 1; i <= m; ++i) {
            for(int j = 1; j <= n; ++j) scanf("%lf", &A[n + i][j]);
            scanf("%lf", &b[n + i]);
        }
        init();
        simplex();
        printf("Nasa can spend %d taka.
", (int)ceil(v * m));
    }
}

6、高斯消元（-1无解，0无穷解，1唯一解）

int guess_eliminatioin() {
    int rank = 0;
    for(int i = 0, t = 0; i < m && t < n; ++i, ++t) {
        int r = i;
        for(int j = i + 1; j < m; ++j)
            if(mat[r][t].val < mat[j][t].val) r = j;
        if(mat[r][t].isZero()) { --i; continue;}
        else ++rank;
        if(r != i) for(int j = 0; j <= n; ++j) swap(mat[i][j], mat[r][j]);
        for(int j = n; j >= t; --j)
            for(int k = i + 1; k < m; ++k) mat[k][j] -= mat[i][j] * mat[k][t] / mat[i][t];
    }
    for(int i = rank; i < m; ++i)
        if(!mat[i][n].isZero()) return -1;
    if(rank < n) return 0;
    for(int i = n - 1; i >= 0; --i) {
        for(int j = i + 1; j < n; ++j)
            mat[i][n] -= mat[j][n] * mat[i][j];
        mat[i][n] = mat[i][n] / mat[i][i];
    }
    return 1;
}

7、离散对数（小步大步算法）（POJ 3243）

#include <cstdio>
#include <iostream>
#include <algorithm>
#include <cstring>
#include <cmath>
using namespace std;
typedef long long LL;

const int SIZEH = 65537;

struct hash_map {
    int head[SIZEH], size;
    int next[SIZEH];
    LL state[SIZEH], val[SIZEH];

    void init() {
        memset(head, -1, sizeof(head));
        size = 0;
    }

    void insert(LL st, LL sv) {
        LL h = st % SIZEH;
        for(int p = head[h]; ~p; p = next[p])
            if(state[p] == st) return ;
        state[size] = st; val[size] = sv;
        next[size] = head[h]; head[h] = size++;
    }

    LL find(LL st) {
        LL h = st % SIZEH;
        for(int p = head[h]; ~p; p = next[p])
            if(state[p] == st) return val[p];
        return -1;
    }
} hashmap;

void exgcd(LL a, LL b, LL &x, LL &y) {
    if(!b) x = 1, y = 0;
    else {
        exgcd(b, a % b, y, x);
        y -= x * (a / b);
    }
}

LL inv(LL a, LL n) {
    LL x, y;
    exgcd(a, n, x, y);
    return (x + n) % n;
}

LL pow_mod(LL x, LL p, LL n) {
    LL ret = 1;
    while(p) {
        if(p & 1) ret = (ret * x) % n;
        x = (x * x) % n;
        p >>= 1;
    }
    return ret;
}

LL BabyStep_GiantStep(LL a, LL b, LL n) {
    for(LL i = 0, e = 1; i <= 64; ++i) {
        if(e == b) return i;
        e = (e * a) % n;
    }
    LL k = 1, cnt = 0;
    while(true) {
        LL t = __gcd(a, n);
        if(t == 1) break;
        if(b % t != 0) return -1;
        n /= t; b /= t; k = (k * a / t) % n;
        ++cnt;
    }
    hashmap.init();
    hashmap.insert(1, 0);
    LL e = 1, m = LL(ceil(sqrt(n + 0.5)));
    for(int i = 1; i < m; ++i) {
        e = (e * a) % n;
        hashmap.insert(e, i);
    }
    LL p = inv(pow_mod(a, m, n), n), v = inv(k, n);
    for(int i = 0; i < m; ++i) {
        LL t = hashmap.find((b * v) % n);
        if(t != -1) return i * m + t + cnt;
        v = (v * p) % n;
    }
    return -1;
}

int main() {
    LL x, z, k;
    while(cin>>x>>z>>k) {
        if(x == 0 && z == 0 && k == 0) break;
        LL ans = BabyStep_GiantStep(x % z, k % z, z);
        if(ans == -1) puts("No Solution");
        else cout<<ans<<endl;
    }
}

8、母函数

/*
G(x) = 1 + x + x^2 + x^3 + …… = 1 / (1 - x)
G(x) = 1 + 2x + 3x^2 + 4x^3 + …… = 1 / (1 - x)^2

指数型母函数
G(x) = 1 + x + x^2 / 2! + x^3 / 3! + x^4 / 4! + …… = e^x
<1, -1, 1, -1, 1, -1, ……> = e^(-x)
<0, 1, 0, 1, 0, 1, ……> = (e^x - e^(-x)) / 2
<1, 0, 1, 0, 1, 0, ……> = (e^x + e^(-x)) / 2
*/

9、插头DP（URAL1519 括号表示法）

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>
using namespace std;
typedef long long LL;

const int MAXH = 20010;
const int SIZEH = 13131;

struct hash_map {
    int head[SIZEH];
    int next[MAXH], state[MAXH];
    LL value[MAXH];
    int size;

    void init() {
        memset(head, -1, sizeof(head));
        size = 0;
    }

    void insert(int st, LL tv) {
        int h = st % SIZEH;
        for(int i = head[h]; ~i; i = next[i]) {
            if(state[i] == st) {
                value[i] += tv;
                return ;
            }
        }
        value[size] = tv; state[size] = st;
        next[size] = head[h]; head[h] = size++;
    }
} hashmap[2];

hash_map *cur, *last;
int acc[] = {0, -1, 1, 0};

int n, m, en, em;
char mat[14][14];

int getB(int state, int i) {
    i <<= 1;
    return (state >> i) & 3;
}

int getLB(int state, int i) {
    int ret = i, cnt = 1;
    while(cnt) cnt += acc[getB(state, --ret)];
    return ret;
}

int getRB(int state, int i) {
    int ret = i, cnt = -1;
    while(cnt) cnt += acc[getB(state, ++ret)];
    return ret;
}

void setB(int &state, int i, int tv) {
    i <<= 1;
    state = (state & ~(3 << i)) | (tv << i);
}

void update(int x, int y, int state, LL tv) {
    int left = getB(state, y);
    int up = getB(state, y + 1);
    if(mat[x][y] == '*') {
        if(left == 0 && up == 0) cur->insert(state, tv);
        return ;
    }
    if(left == 0 && up == 0) {
        if(x == n - 1 || y == m - 1) return ;
        int newState = state;
        setB(newState, y, 1);
        setB(newState, y + 1, 2);
        cur->insert(newState, tv);
    } else if(left == 0 || up == 0) {
        if(x < n - 1) {
            int newState = state;
            setB(newState, y, up + left);
            setB(newState, y + 1, 0);
            cur->insert(newState, tv);
        }
        if(y < m - 1) {
            int newState = state;
            setB(newState, y, 0);
            setB(newState, y + 1, up + left);
            cur->insert(newState, tv);
        }
    } else {
        int newState = state;
        setB(newState, y, 0);
        setB(newState, y + 1, 0);
        if(left == 1 && up == 1) setB(newState, getRB(state, y + 1), 1);
        if(left == 1 && up == 2 && !(x == en && y == em)) return ;
        if(left == 2 && up == 2) setB(newState, getLB(state, y), 2);
        cur->insert(newState, tv);
    }
}

void findend() {
    for(en = n - 1; en >= 0; --en)
        for(em = m - 1; em >= 0; --em) if(mat[en][em] == '.') return ;
}

LL solve() {
    findend();
    cur = hashmap, last = hashmap + 1;
    last->init();
    last->insert(0, 1);
    for(int i = 0; i < n; ++i) {
        int sz = last->size;
        for(int k = 0; k < sz; ++k) last->state[k] <<= 2;
        for(int j = 0; j < m; ++j) {
            cur->init();
            sz = last->size;
            for(int k = 0; k < sz; ++k)
                update(i, j, last->state[k], last->value[k]);
            swap(cur, last);
        }
    }
    return last->size ? last->value[0] : 0;
}

int main() {
    scanf("%d%d", &n, &m);
    for(int i = 0; i < n; ++i) scanf("%s", mat[i]);
    cout<<solve()<<endl;
}

10、快速傅里叶变换（HDU1402）

#include <cmath>
#include <algorithm>
#include <cstdio>
#include <iostream>
#include <cstring>
#include <complex>
using namespace std;
typedef complex<double> Complex;
const double PI = acos(-1);

void fft_prepare(int maxn, Complex *&e) {
    e = new Complex[2 * maxn - 1];
    e += maxn - 1;
    e[0] = 1;
    for (int i = 1; i < maxn; i <<= 1)
        e[i] = Complex(cos(2 * PI * i / maxn), sin(2 * PI * i / maxn));
    for (int i = 3; i < maxn; i++)
        if ((i & -i) != i) e[i] = e[i - (i & -i)] * e[i & -i];
    for (int i = 1; i < maxn; i++) e[-i] = e[maxn - i];
}
/* f = 1: dft; f = -1: idft */
void dft(Complex *a, int N, int f, Complex *e, int maxn) {
    int d = maxn / N * f;
    Complex x;
    for (int n = N, m; m = n / 2, m >= 1; n = m, d *= 2)
        for (int i = 0; i < m; i++)
            for (int j = i; j < N; j += n)
                x = a[j] - a[j + m], a[j] += a[j + m], a[j + m] = x * e[d * i];
    for (int i = 0, j = 1; j < N - 1; j++) {
        for (int k = N / 2; k > (i ^= k); k /= 2);
        if (j < i) swap(a[i], a[j]);
    }
}

const int MAXN = 131072;
Complex x1[MAXN], x2[MAXN];
char s1[MAXN / 2], s2[MAXN / 2];
int sum[MAXN];

int main() {
    Complex* e = 0;
    fft_prepare(MAXN, e);
    while(scanf("%s%s",s1,s2) != EOF) {
        int n1 = strlen(s1);
        int n2 = strlen(s2);
        int n = 1;
        while(n < n1 * 2 || n < n2 * 2) n <<= 1;
        for(int i = 0; i < n; ++i) {
            x1[i] = i < n1 ? s1[n1 - 1 - i] - '0' : 0;
            x2[i] = i < n2 ? s2[n2 - 1 - i] - '0' : 0;
        }

        dft(x1, n, 1, e, MAXN);
        dft(x2, n, 1, e, MAXN);
        for(int i = 0; i < n; ++i) x1[i] = x1[i] * x2[i];
        dft(x1, n, -1, e, MAXN);
        for(int i = 0; i < n; ++i) x1[i] /= n;

        for(int i = 0; i < n; ++i) sum[i] = round(x1[i].real());
        for(int i = 0; i < n; ++i) {
            sum[i + 1] += sum[i] / 10;
            sum[i] %= 10;
        }

        n = n1 + n2 - 1;
        while(sum[n] <= 0 && n > 0) --n;
        for(int i = n; i >= 0;i--) printf("%d", sum[i]);
        puts("");
    }
}

11、求解线性同余方程组（POJ 2891 exgcd）

void exgcd(LL a, LL b, LL &d, LL &x, LL &y) {
    if(!b) d = a, x = 1, y = 0;
    else {
        exgcd(b, a % b, d, y, x);
        y -= x * (a / b);
    }
}

int main() {
    LL k, a1, a2, r1, r2;
    while(scanf("%I64d", &k) != EOF) {
        bool flag = true;
        scanf("%I64d%I64d", &a1, &r1);
        for(int i = 1; i < k; ++i) {
            scanf("%I64d%I64d", &a2, &r2);
            if(!flag) continue;
            LL r = r2 - r1, d, k1, k2;
            exgcd(a1, a2, d, k1, k2);
            if(r % d) flag = false;
            LL t = a2 / d;
            k1 = (r / d * k1 % t + t) % t;
            r1 = r1 + a1 * k1;
            a1 = a1 / d * a2;
        }
        printf("%I64d
", flag ? r1 : -1);
    }
}