数论

1、分数的四则运算:用结构体表示

2、高精度计算

3、数论

 3.1 最小公倍数 GCD

int gcd(int a,int b){
	//辗转相除
	if(b==0) return a;
	else return gcd(b,a%b); 
} 
int gcd(int a,int b){
	return !b ? a:gcd(b,a%b);
}

最大公约数  LCM

lcm(a,b)=a*b/gcd(a,b)

 3.2  快速幂

快速幂及扩展的矩阵快速幂,快速计算a^n的值

两种方法

//分治 
int fastpow(int a,int n){
	if(n==1) return a;
	int temp=fastpow(a,n/2);
	if(n%2==1) return temp*temp*a;
	else return temp*temp;
}
//位运算,加上取模 
int fastpow(int a,int n,int mod){
	int base=a;
	int res=1;
	while(n){
		if(n&1) res=res*base%mod;
		base=base*base%mod;
		n>>=1;
	}
	return res;
} 

矩阵快速幂,但是难点是把递推关系转化为矩阵的关系,例如:斐波那契数列

struct matrix{
	int m[maxn][maxn];
	matrix(){
		memset(m,0,sizeof(m));
	}
}; 
matrix multi(matrix a,matrix b){
	matrix res;
	for(int i=0;i<maxn;i++){
		for(int j=0;j<maxn;j++){
			for(int k=0;k<maxn;k++){
				res.m[i][j]=(res.m[i][j]+a.m[i][k]*b.m[k][j])%mod;
			}
		}
	}
	return res;
}
matrix fastm(martix a,int n){
	matrix res;
	for(int i=0;i<maxn;i++) 
	res.m[i][i]=1;  //初始化为单位矩阵 就像res=1
	while(n){
		if(n&1) res=multi(res,a);
		a=multi(a,a);
		n>>=1;
	} 
	return res;
}

  

3.3  求素数

普通的做法(试除法),即从2~sqrt(n)之内判断的----n<=10^12 以内可以接受   O(√N)

判断:[2,sqrt(n)]内的所有数-----[2,sqrt(n)]内的所有素数 

bool is_prime(int n){
	if(n==1) return false;
	for(int i=2;i*i<=n;i++){
		if(n%i==0) return false;
	}
	return true;
}

埃氏筛法:O(NlongNlongN)

因为空间用到了is[maxn],所以当maxn=1e7时是可以接受的,不然空间就太大了

优化(1)用来做筛的数到sqrt(n)就可以了(2)for(int j=i+i 可以改为j=i*i,因为前面已经筛过了

int prime[maxn];
bool is[maxn];
int ans=0;
void findd(){
	for(int i=2;i*i<=maxn;i++){  //是< 
		if(is[i]==0){
			prime[ans++]=i;
			for(int j=i*i;j<maxn;j+=i){
				p[j]=1;
			}
		}
	}
}

线性筛模板

void prim(){
    for(int i=2;i<=n;i++){
        if(!vis[i]){
            prime[++ans]=i;
        }
        for(int j=1;j<=ans&&prime[j]*i<=n;j++){
                vis[prime[j]*i]=1;
                if(i%prime[j]==0) break;
        }
    }
}

  

大区间素数

埃氏筛法能解决n<=1e7的问题,但是如果n更大,那么就可以扩展位大区间素数,例如求区间[a,b]内的素数,a<b<=1e12, b-a<=1e6

先用埃氏筛法求出[2,sqrt(b)]范围内的素数,再用这些素数去筛空间[a,b]内的素数

 题目:1619: 【例 1】Prime Distance

3.4    扩展欧几里得算法与二元一次方程的整数解

给出整数a,b,n,问方程ax+by=n什么时候有整数解

!!!有解的充要条件是gcd(a,b)可以整除n

如果确定有解,那么解题方法是先找到一个可行解(x0,y0),然后用通解公式

x=x0+bt

y=y0-at求解,那么问题就是如何求(x0,y0)

扩展欧几里得方法:当方程符合ax+by=gcd(a,b)时,可以用扩展欧几里得算法求(x0,y0)

void extend_gcd(int a,int b,int &x,int &y){
	if(b==0){
		x=1;y=0;
		return;
	}
	extend_gcd(b,a%b,x,y);
	int temp=x;
	x=y;
	y=temp-(a/b)*y;
}

求出来的x,y就是所得到x0,y0,那么通解公式就是:

x=x0+b/gcd*t   gcd=gcd(a,b)

y=y0-a/gcd *t

x对应的最小非负整数解为(x%(b/gcd))+b/gcd)%(b/gcd)

那么问题ax+by=n的整数解

利用上面的结果,步骤如下:

(1)判断ax+by=n是否有整数解,也就是gcd(a,b)能不能整除n

(2)扩展欧几里得求得x0,y0

(3)在ax0+by0=gcd(a,b)两边同时乘以n/gcd(a,b)

(4)对照ax+by=n,得到解为:

x0'=x0*n/gcd(a,b)     y0'=y0*n/gcd(a,b)

那么通解公式为:

x=x0'+b/gcd*k

y=y0'-a/gcd *k

运用场合:(1)求解不定方程(2)求解模的逆元(3)求解同余方程

3.5  同余与逆元

a%m=b%m,那么a和b对m同余,m称为同余的模,记为m|(a-b),即a-b是m的整数倍,同余的符号记为a=b(mod m) (等号三横)

逆元:给出a,m,求方程ax=1(mod m),即ax除以m的余数为1

有解的条件是gcd(a,m)=1,即a,m互素,等价于求ax+my=1

方程ax=1(mod m)的一个解x,称x为a模m的逆元。这样的x很多,都称为逆

int mod_inverse(int a,int m){
	int x,y;
	extend_gcd(a,m,x,y);
	return (m+x%m)%m;  //可能是负数,所以需要处理 
}

  逆元的应用:除法的模,求(a/b)mod m,a,b很大,但是求除法之后再取模会损失精度

把除法的模运算转换成了乘法模运算

(a/b)modm=((a/b)modm)((bk)modm)=((a/b)*(bk))modm=akmodm

一元线性同余方程

ax=b(mod m),即ax-b是m的整数倍,得到ax+my=b,当且仅当gcd(a,m)能整除b是有解

当gcd(a,m)=b时,能够运用扩展欧几里得方法直接求解ax+my=b

当gcd(a,m)!=b 时,就需要运用逆元

如果a'是a模m的逆,则a'a=1mod(m),那么在ax=b(mod m)两边同时乘以a',那么就是x=a'b(mod m)

步骤:求解ax+my=b   同余方程为ax=b(mod m)

(1)有解的条件:gcd(a,m)能够整除b

(2)求ax=1(mod m)得逆元a',等价于扩展欧几里得算法求出ax+my=1

(3)一个特解是x=a'b

(4)代入方程ax+my=b,求解y

中国剩余定理(孙子定理)

https://www.cnblogs.com/wkfvawl/p/9633188.html

 但是这样的中国剩余定理有一个要求就是m1,m2,m3...mn是互素的

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
using namespace std;
const int maxn=1010;
const int INF=0x3fffffff;
typedef long long LL;
typedef unsigned long long ull;
LL n;
//中国剩余定理(孙子定理)
//https://www.cnblogs.com/wkfvawl/p/9633188.html 
LL aa[12],bb[12];
void extend_gcd(LL a,LL b,LL &x,LL &y){
    if(b==0){
        x=1;
        y=0;return;
    }
    extend_gcd(b,a%b,x,y);
    LL tmp=x;
    x=y;
    y=tmp-(a/b)*y;
} 
LL gcd(LL x,LL y){
    if(y==0) return x;
    else return gcd(y,x%y);
}

int main(){
    scanf("%lld",&n);
    
    LL mod=1;
    for(int i=1;i<=n;i++){
        scanf("%lld %lld",&aa[i],&bb[i]);
        mod*=aa[i];
    } 
    LL ans=0;
    LL x,y;
    for(int i=1;i<=n;i++){
        LL p=mod/aa[i];
        extend_gcd(p,aa[i],x,y);
        //cout<<x<<endl;
        ans=(ans+bb[i]*x*p)%mod;
    }
    cout<<(ans+mod)%mod;
    /*   超时到怀疑人生哈哈哈哈 
    for(LL i=mm;;i++){
        bool flag=0;
        for(int j=1;j<=n;j++){
            if(!judge(i-bb[j],aa[j])) {
                flag=1;break;
            }
        }
        if(!flag){
            cout<<i<<endl;break;
        }
    }
    */
return 0;
}
View Code

扩展中国剩余定理

扩展中国剩余定理跟中国剩余定理没半毛钱关系,一个是用扩展欧几里得,一个是用构造

https://www.cnblogs.com/zwfymqz/p/8425731.html

判定有无解:(两个一次处理)

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
using namespace std;
const int maxn=1e5+10;
const int INF=0x3fffffff;
typedef long long LL;
typedef unsigned long long ull;
//变形是需要判断有没有解
//中国剩余定理的要求就是:m1,m2....mn两两互素
//但是如果不是两两互素的话,那么就需要扩展中国剩余定理
//https://www.cnblogs.com/zwfymqz/p/8425731.html 
LL n;
//x=c1(mod m1)
LL c[maxn],m[maxn],x,y;
LL gcd(LL a,LL b){
    if(b==0) return a;
    return gcd(b,a%b);
}
LL extend_gcd(LL a,LL b,LL &x,LL &y){   //不仅算了一元二次方程的解,还算出来了公约数 
    if(b==0){
        x=1;y=0;return a;
    }
    LL r=extend_gcd(b,a%b,x,y);
    LL tmp=x;
    x=y;
    y=tmp-(a/b)*y;
    return r;
}
LL inv(LL a,LL b){  //a在模b下的逆
    LL r=extend_gcd(a,b,x,y);  //最大公约数 
    
    while(x<0) x+=b; //取最小非负的逆 
    return x;
} 
int main(){
    while(~scanf("%lld",&n)){
        for(int i=1;i<=n;i++){
            scanf("%lld %lld",&m[i],&c[i]);
        }
        bool flag=0;
        for(int i=2;i<=n;i++){  //两两合并位1个,并继续扩展 
            LL m1=m[i-1],m2=m[i];
            LL c1=c[i-1],c2=c[i];
            LL t=gcd(m1,m2);
            if((c2-c1)%t!=0){
                flag=1;break;  //必须余数为0 因为式子中要出现相除 
            }
            m[i]=(m1*m2)/t;
            c[i]=(inv(m1/t,m2/t)*(c2-c1)/t)%(m2/t)*m1+c1;
            c[i]=(c[i]%m[i]+m[i])%m[i];   //还有这一步,别忘了 
        } 
        printf("%lld
",flag? -1:c[n]);
    }
return 0;
}
View Code

3.6   质因子分解

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

const int maxn=100010;
bool is_prime(int x){
	if(x==1) return 0;
	for(int i=2;i*i<=x;i++){
		if(x%i==0) return 0;
	}
	return 1;
}
//2147483646
int prime[maxn],ans=0;
struct node{
	int x,num; //用结构体存 
}fac[10];  //10个就够了 
void find_prime(){  //求素数表 
	for(int i=1;i<maxn;i++){
	if(is_prime(i)) prime[ans++]=i; 
}
}
int main(){
	find_prime();
	int n,num=0;
	cin>>n;
	if(n==1) {
		cout<<"1"<<endl;
		return 0; //特判 
	}
	cout<<n<<"=";
	int sqr=(int)(sqrt(1.0*n));
	for(int i=0;prime[i]<=sqr&&i<ans;i++){
		if(n%prime[i]==0){
			fac[num].x=prime[i];
			fac[num].num=0;
			while(n%prime[i]==0){
				fac[num].num++;
				n/=prime[i];
			}
			num++;
		}
		if(n==1) break; //提前退出 
	}
	if(n!=1){ //最后不等于1 ,说明有一个大于sqrt(n)的因子
	fac[num].x=n;
	fac[num++].num=1; 
		
	}
	for(int i=0;i<num;i++){
		if(i>0) cout<<"*";
		cout<<fac[i].x;
		if(fac[i].num!=1) cout<<"^"<<fac[i].num;
	}
	return 0;
}

ps.如果进行因子分解后,得到各因子的个数位e1,e2...ek,那么因子个数就是(e1+1)*(e2+1)*...*(ek+1)

因子之和就是(1+p1+p1^2+...+p1^e1)*(1+p2+p2^2+...p2^e2)*...*(1+pk+pk^2+...+pk^ek)

3.7   关于n!的问题

1、求n!中有多少个质因子p

可以枚举,但是很慢

n!中有(n/p+n/p^2+.....)个质因子p.

eg.10!中有5个2^1,2个2^2,1个2^3,所以质因子2的个数位5+2+1=8个,O(logN)

int cal(int n,int p){
	int ans=0;
	while(n){
		ans+=n/p;
		n/=p;
	} 
	return ans;
} 

应用:计算n!的末尾有几个0,末尾0的个数等于n!中因子10的个数,等于质因子5的个数,计算cal(n,5)就可以了

递归版本:

int cal(int n,int p){
	if(n<p) return 0;
	return n/p+cal(n/p,p); 
}

3.8   组合数的计算

一、计算Cnm

方法一:通过定义计算,即使用Long long也只能承受n<=20;

方法二:递推或者递归,单次计算不会超过O(N^2),记忆化

res[n][m]=js(n-1,m)+js(n-1,m-1)

res[i][j]=res[i][i-j]

long long res[67][67]; //n=67,m=33时溢出
long long js(long long n,long long m){
	if(m==0||m==n) return 1;
	if(res[n][m]!=0) return res[n][m];
	return res[n][m]=js(n-1,m)+js(n-1,m-1);
} 

//计算整张表
void js(){
	int n=60;
	for(int i=0;i<=n;i++){
		res[i][0]=res[i][i]=1;
	}
	for(int i=2;i<=n;i++){
		for(int j=1;j<=i/2;j++){
			res[i][j]=res[i-1][j]+res[i-1][j-1];
			res[i][i-j]=res[i][j];
		}
	}
} 

方法三:定义式变形计算

ans=ans*(n-m+i)/i   i=1~m、

//根据定义式,划分为计算m次
//n=62,m=31时溢出 long long js(long long n,long long m){ long long ans=0; for(long long i=1;i<=m;i++){ ans=ans*(n-m+i)/i; } return ans; }

二、计算Cnmmodp

第一种:递推或者递归

int res[1010][1010];  //m<=n<=1000 ,对p的素性和大小没有限制
int js(int n,int m,int p){
	if(m==0||m==n) return 1;
	if(res[n][m]!=0) return res[n][m];
	return res[n][m]=(js(n-1,m)+js(n-1,m-1))%p;
}

第二种:定义式计算

将组合数进行质因子分解,那么Cn,m=p1^c1*p2^c2...*pn^cn%p,通过快速幂计算每一组pi^ci%p,然后相乘取模就可以了

对Cn,m进行质因子分解:遍历不超过n的所有质数,计算n!,m!,(n-m)!中分别含质因子pi的个数x,y,z,那么Cn,m中含质因子pi的个数为x-y-z

这里就用到了前面的求n!里面的质因子p的个数cal(int n,int p)

m<=n<=10^6,对p的素性没有限制

int prime[maxn]; //素数表(已经得到的,埃氏筛法) 
int js(int n,int m,int p){
	int ans=1;
	for(int i=0;prime[i]<=n;i++){
		int c=cal(n,prime[i])-cal(m,prime[i])-cal(n-m,prime[i]);
		ans=ans*binpow(prime[i],c,p)%p;  //快速幂计算prime[i]^c%p 
	}
}

  

高次同余方程   BSGS

#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<stack>
#include<cstdio>
#include<queue>
#include<map>
#include<vector>
#include<set>
using namespace std;
const int maxn=1010;
const int INF=0x3fffffff;
typedef long long LL;
typedef unsigned long long ull;
//1.y^zmod p
//2. x*y=z(mod p)
//3. y^x=z(mod p)

LL   js1(LL y,LL z,LL p){
    LL tmp=1;
    while(z){
        if(z&1) tmp=tmp*y%p;
        y=y*y%p;
        z>>=1;
    }
    return tmp%p;
}
void extend_gcd(LL a,LL b,LL &x,LL &y){
        if(b==0){
        x=1;
        y=0;return;
    }
    extend_gcd(b,a%b,x,y);
    LL tmp=x;
    x=y;
    y=tmp-(a/b)*y;
}
LL gcd(LL a,LL b){
    if(b==0 ) return a;
    else return gcd(b,a%b);
}
//x*y=z(mod p)  ax=b(mod m)
//ax+my=b
//y*x+p*k=z  已知y,p,z 
//求出x0之后
//最小非负解:((x=x0*z/gcd)+p/gcd)%(p/gcd)
void js2(LL y,LL z,LL p){
    LL A=y,B=z,M=p;
    LL xx,yy;
    extend_gcd(A,M,xx,yy);
    //gcd(a,m)%b==0
    LL d=gcd(y,p);
    if(z%d){
        printf("Orz, I cannot find x!
");return;
    }
    xx=((xx*z/d)%(p/d)+p/d)%(p/d);
    printf("%lld
",xx);
    //return (xx%M+M)%M;
}
//y^x=z(mod p)
//BSGS可以解决高次同余方程:
//https://www.cnblogs.com/kamimxr/p/11555986.html

void js3(LL a,LL ans,LL p){
    map<LL,LL> myhash;
    ans%=p;
    int tmp=sqrt(p)+1;
    for(int i=0;i<tmp;i++){
        //js1其实就是快速幂 
        myhash[(ans)*js1(a,i,p)%p]=i;
    }
    a=js1(a,tmp,p)%p;
    if(a==0&&ans==0){
        printf("1
");return;
    }
    if(a==0&&ans!=0){
        printf("Orz, I cannot find x!
");
        return;
    }
    for(int i=0;i<=tmp;i++){
        if(myhash.find(js1(a,i,p))!=myhash.end()&&(i*tmp-myhash[js1(a,i,p)]>=0)){
            printf("%d
",i*tmp-myhash[js1(a,i,p)]);
            return;
        }
    }
    printf("Orz, I cannot find x!
");
} 
int main(){
    int t,k;
    LL y,z,p;
    while(~scanf("%d %d",&t,&k)){
        if(k==1){
            for(int i=1;i<=t;i++){
                scanf("%lld %lld %lld",&y,&z,&p);
                printf("%lld
",js1(y,z,p));
            }
        }
        else if(k==2){
            for(int i=1;i<=t;i++){
                scanf("%lld %lld %lld",&y,&z,&p);
                js2(y,z,p);
            }
        }
        else if(k==3){
            ////y^x=z(mod p)
            for(int i=1;i<=t;i++){
                scanf("%lld %lld %lld",&y,&z,&p);
                js3(y,z,p);
            }
        }
    }
return 0;
}
View Code
原文地址:https://www.cnblogs.com/shirlybaby/p/12264223.html