基础函数:

// 最大公约数,欧几里得定理
int gcd(int a, int b)
{
    return b?gcd(b, a % b): a;
}
// 拓展欧几里得定理
// 求解ax + by = gcd(a,b)
int ext_gcd(int a, int b, int &x, int &y)
{
    int tmp, ret;
    if(!b)
    {
        x = 1;
        y = 0;
        return a;
    }
    ret = ext_gcd(b, a % b, x, y);
    tmp = x;
    x = y;
    y = tmp - (a / b) * y;
    return ret;
}
//交换数值
void swap(int &a, int &b)
{
    a ^= b ^= a ^= b;
}

/**
 * a的b次方Mod c
 * 参数为整数
 * 使用时注意修改类型
 */
int PowerMod(int a, int b, int c)
{
    int tp = 1;
    while (b)
    {
        if (b & 1)
            tp = (tp * a) % c;
        a = (a * a) % c;
        b >>= 1;
    }
    return tp;
}

1.欧拉函数

Ψ（n） = 小于n且与n互质的数的个数

int eular(int n)
{
    int res = 1, i;
    for(i = 2; i * i < res; i ++)
    {
        if(n % i == 0)
        {
            n /= i;
            res *= (i - 1);
            while(n % i == 0)
                n /= i, res *= i;
        }
    }
    if(n > 1)
        res *= n - 1;
    return res;
}

2.欧拉定理

若a与n互质（即GCD（a，n） = 1），则a^Ψ（n） = 1 （mod n）a^{\varphi(n)} \equiv 1 \pmod n

欧拉函数的一个定理：Ψ（n）= n – sum{Ψ（x）| 其中 n % x == 0}

3.排列组合

/**
 * 排列组合数（素数表示法）
 * 注意传入的数组先初始化为0
 * 复杂度：O(nlog(n)),n为素数个数
 */
// 全排列
// 参数: A(n), p[] 传出数的数组表示指针
// 返回值：结果包含的素数个数
int Arrangement(int n, int p[])
{
    int t, i;
    for(i = 0; i < num_prime && prime[i] <= n; i ++)
    {
        t = n;
        while(t)
            p[i] += t / prime[i], t /= prime[i];
    }
    return i;
}
// 排列
// 参数: A(m,n),m >= n , p[] 传出数的数组表示指针
// 返回值：结果包含的素数个数
int A_Cache[2][MAXP];//缓存项
int Arrangement_A(int m, int n, int p[])
{
    int i;
    ::memset(A_Cache, 0, sizeof(A_Cache));
    int r = Arrangement(m, A_Cache[0]);
    Arrangement(n, A_Cache[1]);
    for(i = 0; i < num_prime; i ++)
        p[i] = A_Cache[0][i] - A_Cache[1][i];
    return r;
}
// 组合
// 参数: C(m,n),m >= n , p[] 传出数的数组表示指针
// 返回值：结果包含的素数个数
int C_Cache[3][MAXP];//缓存项
int Combination(int m, int n, int p[])
{
    int i;
    memset(C_Cache, 0, sizeof(C_Cache));
    int r = Arrangement(m, C_Cache[0]);
    Arrangement(n, C_Cache[1]);
    Arrangement(m - n, C_Cache[2]);
    for(i = 0; i < num_prime; i ++)
        p[i] = C_Cache[0][i] - C_Cache[1][i] - C_Cache[2][i];
    return r;
}

// 取模计算:参数: mod为取模的值，其他参数同上
// 全排列取模
int Arrangement_Mod(int n, int p[], int mod)
{
    int res = 1;
    int i, r = Arrangement(n, p);
    for(i = 0; i < r; i ++)
        res = (res * PowerMod(prime[i], p[i], mod)) % mod;
    return res;
}
// 排列取模
int Arrangement_A_Mod(int m, int n, int p[], int mod)
{
    int res = 1;
    int i, r = Arrangement_A(m, n, p);
    for(i = 0; i < r; i ++)
        res = (res * PowerMod(prime[i], p[i], mod)) % mod;
    return res;
}
// 组合取模
int Combination_Mod(int m, int n, int p[], int mod)
{
    int res = 1;
    int i, r = Combination(m, n, p);
    for(i = 0; i < r; i ++)
        res = (res * PowerMod(prime[i], p[i], mod)) % mod;
    return res;
}

/**
 * 大组合数对素数取模
 * 复杂度: O(n),n素数大小
 */

// 大整数组合取模
const int MAXMOD = 10007;
int C_Map[MAXMOD + 1] = {0};
void BC_Init(int n, int mod)//初始化调用
{
    int i;
    C_Map[1] = 1;
    for(i = 2; i < n; ++ i)
        C_Map[i] = (C_Map[i - 1] * i) % mod;
}
int BC_Combination(int m, int n, int mod)
{
    int x, y;
    if(m < n)
        return 0;
    else
    {
        if(n == 0)
            return 1;
        int mm, nn;
        mm = C_Map[m];
        nn = (C_Map[n] * C_Map[m - n]) % mod;
        int shr = gcd(mm, nn);
        mm /= shr;
        nn /= shr;
        ext_gcd(nn, mod, x, y);
        return (x + mod) * mm % mod;//注意越界
    }
}
//调用的函数
int BN_Combination(int m, int n, int mod)
{
    int m_Num[50] = {0}, n_Num[50] = {0}, i, res = 1;
    //转换为mod进制数
    for(i = 0; m; i ++)
        m_Num[i] = m % mod, m /= mod;
    int t = i;
    for(i = 0; n; i ++)
        n_Num[i] = n % mod, n /= mod;
    for(i = t - 1; i >= 0; i --)
        res = (res * BC_Combination(m_Num[i], n_Num[i], mod)) % mod;//注意越界
    return res;
}

4.分数类+高斯消元

/**
 * 高斯消元和与之配合的分数类
 * 高斯消元复杂度O(n^3),n为未知数个数
 * #include <cmath>
 */

/**
 * 分数类(注意越界)
 */
struct mark
{
    int c, m;
    mark(){}
    mark(int x):c(x), m(1){}
    mark(int _c, int _m)
    {
        int d = gcd(::abs(_c), ::abs(_m));
        c = _c / d;
        m = _m / d;
        if(c < 0 && m < 0)
            c *= -1, m *= -1;
    }

    mark operator+(const mark &r) const
    {
        return mark(c * r.m + m * r.c, m * r.m);
    }

    mark operator-(const mark &r) const
    {
        return mark(c * r.m - m * r.c, m * r.m);
    }

    mark operator*(const mark &r) const
    {
        return mark(c * r.c, m * r.m);
    }

    mark operator/(const mark &r) const
    {
        return mark(c * r.m, m * r.c);
    }

    mark pow(int t)
    {
        mark tp = (*this);
        mark a = mark(1, 1);
        while(t > 0)
        {
            if(t & 1)
                a = a * tp;
            tp = tp * tp;
            t >>= 1;
        }
        return a;
    }

    bool operator==(const mark &r) const
    {
        return ((*this) - r).c == 0;
    }
    //其他判断类似
    bool operator<(const mark &r) const
    {
        mark rk = (*this) - r;
        return rk.c * rk.m < 0;//小心越界
    }
    bool operator>=(const mark &r) const
    {
        return !((*this) < r);//小心越界
    }
};

mark abs(const mark &x)
{
    return mark(::abs(x.c), ::abs(x.m));
}

/**
 * 高斯消元(求解：a[i][j] * x[j] = b[j])
 * 复杂度: O(n^3)
 * 可导入分数类(修改类型时修改zero函数,mark类型和abs函数即可)
 */
struct gauss_mat
{
    static const int maxn = 100;//最大未知数数量
    mark mat[maxn][maxn + 1];//增广矩阵
    mark x[maxn];//解集

    //浮点型和分数类型数据注意修改这里
    bool zero(const mark &x)
    {
        return x.c == 0;
    }

    gauss_mat(){}
    //构造sum{a[i][j] * x[j]} = b[i]
    gauss_mat(mark a[][maxn], mark b[], int _n)
    {
        for(int i = 0; i < _n; i ++)
        {
            for(int j = 0; j < _n; j ++)
                mat[i][j] = a[i][j];
            mat[i][_n] = b[i];
        }
    }

    //获取解
    mark& operator[](int p)
    {
        return x[p];
    }

    //交换行
    void swapR(int r1, int r2, int n)
    {
        for(int i = 0; i <= n; i ++)
        {
            mat[r1][i] = mat[r1][i] + mat[r2][i];
            mat[r2][i] = mat[r1][i] - mat[r2][i];
            mat[r1][i] = mat[r1][i] - mat[r2][i];
        }
    }

    //高斯消元（整数）
    //返回0为有无穷解或无解，返回1有唯一解并计算答案，返回-1无解
    bool gauss(int n)
    {
        int i, j, k, pj;
        for(i = 0; i < n; i ++)
        {
            pj = i;
            //注意类型
            mark p = mat[i][i];
            for(j = i + 1; j < n; j ++)
                if(::abs(p) < ::abs(mat[j][i]))
                    p = mat[j][i], pj = j;
            if(zero(p))
                return false;
            else if(i != pj )
                swapR(i, pj, n);
            for(j = i + 1; j < n; j ++)
            {
                /*
                //-----以下三选一-----
                //参数都是int(易越界)
                int d = gcd(::abs(mat[j][i]), ::abs(p));
                int lp = p / d, rp = mat[j][i] / d;
                d = (mat[j][i] * p < 0)? 1: -1;
                for(k = i; k <= n; k ++)
                    mat[j][k] = lp * mat[j][k] + d * rp * mat[i][k];
                //------int------

                //参数是double
                double d = mat[j][i] / p;
                for(k = i; k <= n; k ++)
                    mat[j][k] = mat[j][k] - d * mat[i][k];
                //------double------
                */
                //参数是mark(为了不易越界，整数建议用分数表示)
                //注意可以删除一些mark标签
                mark d = mat[j][i] / p;
                for(k = i; k <= n; k ++)
                    mat[j][k] = mat[j][k] - d * mat[i][k];
            }
        }

        for(i = n - 1; i >= 0; i --)
        {
            //注意能被整除，否则需要改变类型
            if(zero(mat[i][i]))
                return false;
            mark p = mat[i][i];
            mark sum = mat[i][n];
            for(j = i + 1; j < n; j ++)
                sum = sum -  mat[i][j] * x[j];
            x[i] = sum / p;
        }
        return true;
    }
};