算法模板总结——数学部分

为省赛整理的常用算法模板,数学部分

0-20的阶乘

1
2
3
4
5
const long long fac[21]={1,1,2,6,24,120,720,5040,40320,362880,
3628800,39916800,479001600,6227020800,
87178291200,1307674368000,20922789888000,
355687428096000,6402373705728000,121645100408832000,
2432902008176640000};

错排公式

有n个元素组成的排列,若所有的元素都不在自己原来的位置上,这样的排列记为D(n),有

D(n)=(n-1)[D(n-1)+D(n-2)]
理解博客


最小公倍数lcm(a,b)和最大公因数gcd(a,b)

1
2
3
4
5
6
7
8
int gcd(int a,int b)
{
return b==0?a:gcd(b,a%b);
}
int lcm(int a,int b)
{
return a*b/gcd(a,b);
}

扩展欧几里得原理

求解 ax + by = gcd(a; b) 这里得到的是一组 (x, y) 的可行解,(x + kb, y - ka) 为解集
如果是ax + by = d,这样的话得到的解就是 x0=x*(d/gcd(a,b) ) y0=y*(d/gcd(a,b) )
最后的解集就是:(x0 + k*(b/gcd(a,b) ), y0 - k*(a/gcd(a,b) ) )

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int e_gcd(int a,int b,int &x,int &y)
{
int d=a;
if(!b)
{
x=1;
y=0;
}
else
{
d=e_gcd(b,a%b,y,x);
y-=(a/b)*x;
}
return d;
}


快速幂取模算法

1
2
3
4
5
6
7
8
9
10
11
ll quick_mod(ll a,ll b,ll c)
{
ll ans=1;
while(b>0)
{
if(b&1) ans=(ans*a)%c;
a=(a*a)%c;
b>>=1;
}
return ans;
}

高精度问题

充分利用整数的取值范围,每4位进行处理,减小时间复杂度

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
struct BigInt
{
const static int nlen=4; //控制数组中的每一个数字的长度,为了乘法运算不溢出设定为4
const static int mod=10000; //每个数字大小设定
short n[1000],len; //存放数字的数组以及数组的长度
BigInt()//没有赋值时初始化为0
{
memset(n,0,sizeof(n));
len=1;
}
BigInt(int num)//数字为其赋值时,将数字4位4位存放在数组当中
{
len=0;
while(num>0)
{
n[len++]=num%mod;
num/=mod;
}
}
BigInt(const char *s) //字符串赋值时
{
int l=strlen(s);
len=l%nlen==0?l/nlen:l/nlen+1;//确定数组长度
int index=0;
for(int i=l-1;i>=0;i-=nlen)//每次处理数组中的一位
{
int tmp=0;
int j=i-nlen+1;
if(j<0) j=0;//最后面数字的处理
for(int k=j;k<=i;k++)
tmp=tmp*10+s[k]-'0';
n[index++]=tmp;
}
}
BigInt operator+(const BigInt &b)const //加法操作
{
BigInt res;
res.len=max(len,b.len); //确定位数
for(int i=0;i<res.len;i++)
{
res.n[i]+=(i<len?n[i]:0)+(i<b.len?b.n[i]:0); //对象位置相加
res.n[i+1]+=res.n[i]/mod; //进位处理
res.n[i]=res.n[i]%mod;
}
if(res.n[res.len]>0)res.len++; //最后的结果多出一位时
return res;
}
BigInt operator*(const BigInt &b)const //乘法操作
{
BigInt res;
for(int i=0;i<len;i++) //模拟过程
{
int up=0; //进位存储
for(int j=0;j<b.len;j++)
{
int tmp=n[i]*b.n[i]+up+res.n[i+j];
res.n[i+j]=tmp%mod;
up=tmp/mod;
}
if(up!=0) //处理一遍以后还有进位
res.n[i+b.len]=up;
}
res.len=len+b.len; //先取到位数可能最大的值
while(res.n[res.len-1]==0&&res.len>1)res.len--;
return res;
}
void show() const //输出时的逆序输出
{
printf("%d",n[len-1]);
for(int i=len-2;i>=0;i--)
printf("%04d",n[i]); //注意一定要加04,确保输出四位
printf("\n");
}
};

素数

埃式筛选法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
bool is_prime[MAXSIZE];
int prime[MAXSIZE]; //收集素数
void sieve(int n)
{
int p=0;
memset(is_prime,0,sizeof(is_prime));
is_prime[0]=is_prime[1]=false;
for(int i=2;i<n;i++)
{
if(is_prime[i])
{
prime[p++]=i;
for(int j=2*i;j<n;j+=i)
is_prime[j]=false;
}
}
}

欧拉筛法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
int prime[maxn];
bool is_prime[maxn];
int num; //表示素数个数
void init()
{
num=0;
memset(is_prime,true,sizeof(is_prime));
is_prime[0]=is_prime[1]=false;
for(int i=2;i<MAX;i++)
{
if(is_prime[i])prime[num++]=i;
for(int j=0;j<num && i*prime[j]<MAX;j++)
{
is_prime[i*prime[j]]=false;
if(i%prime[j]==0)break; //保证每一个合数都被他的最小质因数排除
}
}
}

区间筛法

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
bool is_prime[maxp];
int prime[maxp],num;
bool is_prime2[maxn]; //大区间内的素数集合
ll prime2[maxn],num2; //大区间内的素数集合
void prime_init()
{
num=0;
memset(is_prime,true,sizeof(is_prime));
is_prime[0]=is_prime[1]=false;
for(int i=2;i<maxp;i++)
{
if(is_prime[i]) prime[num++]=i;
for(int j=0;j<num && i*prime[j]<maxp;j++)
{
is_prime[i*prime[j]]=false;
if(i%prime[j]==0) break;
}
}
}
void bigger_prime(ll L,ll R)
{
memset(is_prime2,1,sizeof(is_prime2));
ll mul;
for(int i=0;i<num && prime[i]<=R;i++)
{
if(prime[i]<=L)
mul=(L-prime[i])/prime[i]; //获得相差的倍数
else
mul=2;
while(mul*prime[i]<L || mul<=1) mul++; //修正倍数值,不能够等于一
for(ll j=mul*prime[i];j<=R;j+=prime[i])
if(j>=L)
is_prime2[j-L]=0;
}
num2=0;
for(int i=0;i<=R-L;i++)
{
if(is_prime2[i])
{
if(i+L==1) continue; //注意要排除1这个不是素数的数
prime2[num2++]=i+L;
}
}
}

大素数判断和分解(Miller-Rabin素数测试和Pollard Rho 大整数分解)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
//Miller_Rabin 算法进行素数测试
//速度快,而且可以判断 <2^63的数
const int S=20; //随机算法判定次数
//计算 (a*b)%c 加法快速幂
ll mul_mod(ll a,ll b,ll c)
{
a%=c;
b%=c;
ll ret=0;
while(b)
{
if(b&1)
ret+=a,ret%=c;
a<<=1;
if(a>=c)a%=c;
b>>=1;
}
return ret;
}
// 计算x^n %c
ll pow_mod(ll x,ll n,ll mod)
{
if(n==1) return x%mod;
x%=mod;
ll tmp=x;
ll ret=1;
while(n)
{
if(n&1) ret=mul_mod(ret,tmp,mod);
tmp=mul_mod(tmp,tmp,mod);
n>>=1;
}
return ret;
}
//以a为基,n-1=x*2^t a^(n-1)=1(mod n) 验证n是不是合数
//一定是合数返回true,不一定返回false
bool check(ll a,ll n,ll x,ll t)
{
ll ret=pow_mod(a,x,n);
ll last=ret;
for(int i=1;i<=t;i++)
{
ret=mul_mod(ret,ret,n);
if(ret==1 && last!=1 && last!=n-1) return true;
last=ret;
}
if(ret!=1) return true;
return false;
}
// Miller_Rabin()算法素数判定
//是素数返回true.(可能是伪素数,但概率极小)
//合数返回false;
bool Miller_Rabin(ll n)
{
if(n<2) return false;
if(n==2) return true;
if((n&1)==0) return false; //偶数
ll x=n-1,t=0;
while(!(x&1))
{
x>>=1;
t++;
}
for(int i=0;i<S;i++)
{
ll a=rand()%(n-1)+1;
if(check(a,n,x,t))
return false;
}
return true;
}
// pollard_rho 算法进行质因数分解
ll factor[100]; //分解结果
int tol; //分解个数
ll gcd(ll a,ll b)
{
if(a==0) return 1;
if(a<0) return gcd(-a,b);
return b==0?a:gcd(b,a%b);
}
ll Pollard_rho(ll x,ll c)
{
ll i=1,k=2;
ll x0=rand()%x;
ll y=x0;
while(1)
{
i++;
x0=(mul_mod(x0,x0,x)+c)%x;
ll d=gcd(y-x0,x);
if(d!=1&&d!=x) return d;
if(y==x0) return x;
if(i==k){y=x0;k+=k;}
}
}
//对n进行素因子分解
void findfac(ll n)
{
if(Miller_Rabin(n))
{
factor[tol++]=n;
return;
}
ll p=n;
while(p>=n)p=Pollard_rho(p,rand()%(n-1)+1);
findfac(p);
findfac(n/p);
}

素性测试

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
bool is_prime(int n)
{
for(int i=2;i*i<=n;i++)
{
if(n%i==0)
return false;
}
return true;
}
vector<int> divisor(int n) //枚举约数及其个数
{
vector<int> res;
for(int i=2;i*i<=n;i++)
{
if(n%i==0)
{
res.push_back(i);
if(n/i!=i)res.push_back(n/i);
}
}
return res;
}
map<int,int> prime_factor(int n) //整数分解
{
map<int,int> res;
for(int i=2;i*i<=n;i++)
{
while(n%i==0)
{
res[i]++;
n/=i;
}
}
if(n!=0)res[n]=1;
return res;
}

求一个数的因子个数

原理:
若一个整型数能够分解质因数!$ M=p_1^{X_1}*p_2^{X_2}*...*p_n^{X_n} $
则这个数的因子个数为 !$ W=(x_1+1)(x_2+1)...(x_n+1) $
根据这个原理我们就可以快速得到数的因子数表

模板:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
void init_factor_table(int n)
{
fill(factor_table, factor_table + n + 1, 1);
for (int i = 2; i <= n; ++i)
{
if (factor_table[i] == 1) //从最小的素数开始计算
{
for (int j = i; j <= n; j += i)
{
int k = 0;
for (int m = j; m % i == 0; m /= i, ++k); //确定当前i的指数
factor_table[j] *= k + 1;
}
}
}
}


中国剩余定理

给出一个数被三个数取摸后的余数,求这个数。设这三个数为a,b,c;取摸之后的余数为a1,b1,c1。求解步骤:

  1. 对于每一个除数a,找出另两个数b,c的一个最小的满足m%a=1的公倍数m
  2. 用这个每个除数对应的余数乘以刚才得到的m,得到的三个数相加得到M
  3. 再用M除以a,b,c的最小公倍数,得到的结果就是想要的结果
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
int CRT(int n,int a[],int mod[])
{
int tem=1,ret=0,t,y;
for(int i=0;i<n;i++)
tem*=mod[i];
for(int i=0;i<n;i++)
{
int w=tem/mod[i];
t=e_gcd(mod[i],w,t,y); //计算出t w模mi下的的逆元
ret=(ret+y*a[i]*w)%tem;
}
while(ret<=0)
ret+=tem;
return ret;
}

高斯消元

原理解释

模板初级版本

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
void swaprow(double a[][MAX],int m,int maxrowe,int n) //交换m和主元行maxrowe
{
for(int k=m;k<=n+1;k++)
{
swap(a[m][k],a[maxrowe][k]);
}
}
void selectcole(double a[][MAX],int n) //选择列主元并进行消元
{
int i,j,k,maxrowe;
double temp;
for(j=1;j<=n;j++)
{
maxrowe=j;
for(i=j;i<=n;i++)
if(abs(a[i][j])>abs(a[maxrowe][j]))
maxrowe=i;
if(maxrowe!=j)
swaprow(a,j,maxrowe,n); //与最大主元所在的行进行交换
//消元,构成三角形矩阵
for(i=j+1;i<=n;i++)
{
temp=a[i][j]/a[j][j];
for(k=j;k<=n+1;k++)
a[i][k]-=a[j][k]*temp;
}
}
}
void gauss(double a[][MAX] ,int n) //最后的答案为a[][n+1]
{
selectcole(a,n); //构成上三角
for(int i=n;i>=1;i--)
{
for(int j=i+1;j<=n;j++)
a[i][n+1]-=a[i][j]*a[j][n+1];
a[i][n+1]/=a[i][i];
}
}
a数组为构成的增广矩阵,最后的解集为a[][n+1]
```
***
## 欧拉函数
小于n且与n互质的数的个数=n*(1-1/P1)*(1-1/P2)….*(1-1/Pn),其中Pn为不同的质因数
``` cpp
typedef long long ll;
map<ll,ll> prime_factor(ll t)
{
map<ll,ll> ret;
for(int i=2;i*i<=t;i++)
{
while(t%i==0) { ++ret[i];t/=i; }
}
if(t!=1) ret[t]=1;
return ret;
}
ll Euler (ll t)
{
ll ret=t;
map<ll,ll> fac=prime_factor(t);
for(map<ll,ll>::iterator i=fac.begin();i!=fac.end();i++)
ret=ret*(i->first-1)/i->first;
return ret;
}


Graham扫描法求凸包

从最左下方的开始,沿图的边界进行扫描,然后进行取舍

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
struct point
{
double x,y;
}p[maxn],t[maxn];
int n; //点的个数
//得到相应的叉乘
double get_X(point a,point b,point c) //ab到ac的叉乘计算
{
return (b.x-a.x)*(c.y-a.y)-(b.y-a.y)*(c.x-a.x);
}
double len(point a,point b)
{
return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}
bool cmp(point &a,point &b)
{
double pp=get_X(p[0],a,b);
if(pp>0) return true;
if(pp<0) return false;
return len(p[0],a)<len(p[0],b);
}
int Graham() //返回凸包所含有的点的个数
{
if(n<=2) return 0;
else
{
for(int i=1;i<n;i++) //找到起始点
{
if(p[i].x<p[0].x) swap(p[i],p[0]);
else if(p[i].x==p[0].x && p[i].y<p[0].y) swap(p[i],p[0]);
}
sort(p+1,p+n,cmp);
t[0]=p[0];
t[1]=p[1];
int top=1; //栈顶位置
for(int i=2;i<n;i++)
{
while(t>0 && get_X(t[top-1],t[top],p[i])<=0) top--; //考虑的点在右侧的时候将栈顶的点弹出
top++;
t[top]=p[i];
}
return top;
}
}