数论模板（个人模板）
基础函数:
// 最大公约数,欧几里得定理
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;
}
};Last updated