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

发表于 | 分类于 | 阅读次数:
字数统计: 3,126 字

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

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;
	}
}