数论的模板

书上的几个模板,照着打了一遍。。还有一些不解的地方~,有待深入了解。。

//求一次同余方程ax=b(mod m)的最小正整数解
#include<cstdio>
#include<iostream>
using namespace std;
__int64 exgcd(__int64 m,__int64 &x,__int64 n,__int64 &y)
{
	__int64 x1,y1,x0,y0;
	x0=1;y0=0;
	x1=0;y1=1;
	__int64 r=(m%n+n)%n;
	__int64 q=(m-r)/n;
	x=0;y=1;
	while(r)
	{
		x=x0-q*x1;y=y0-q*y1;x0=x1;y0=y1;
		x1=x;y1=y;
		m=n;n=r;r=m%n;
		q=(m-r)/n;
	}
	return n;
}
__int64 f(__int64 a,__int64 b,__int64 m)//用exgcd求解最小正整数解
{
	__int64 x,y;
	__int64 M=exgcd(a,x,m,y);
	if(b%M!=0)
		return -1;
	else
	{
		__int64 s=m/M;
		x=x*(b/M);
		x=(x%s+s)%s;
		return x;//x即为ax=b(mod m)的最小正整数解
	}

}
int main()
{
	
return 0;
}
#include<stdio.h>
#include<stdlib.h>
__int64 ans[1000];
void exgcd(__int64 a,__int64 b,__int64 &d,__int64 &x,__int64 &y)
{
	if(b==0)
	{
		x=1;y=0;
		d=a;
		return;
	}
	else
	{
		exgcd(b,a%b,d,x,y);
		__int64 temp=x;
		x=y;
		y=temp-(a/b)*y;
	}
}
__int64 f(__int64 a,__int64 b,__int64 m)//解一元线性同余方程(<m的所有解)
{
	int i;
	__int64 d,x,y;
	exgcd(a,m,d,x,y);
	if(b%d)
		return -1;//表示无解
	x=x*(b/d)%m;
	for(i=1;i<=d;i++)
		ans[i]=(x+(i-1)*m/d)%m;
	return d;

}
int main()
{
	__int64 d,a,b,m;
	int i;
	while(scanf("%I64d %I64d %I64d",&a,&b,&m)!=EOF)
	{
		d=f(a,b,m);
		if(d==-1)
			printf("no answer\n");
		else
		{
			for(i=1;i<=d;i++)
			printf("%I64d\n",ans[i]);
		}
	}
system("pause");
return 0;
}
//若要求所有的解:ans[i]+m*t是解(t为整数)
/*
x=r1(mod a1)
x=r2(mod a2)
x=r3(mod a3)
……
x=rn(mod an)
令m=[a1,a2,a3,……,an]
解上述一元线性同余方程组小于m的非负整数解
*/
void exgcd(__int64 a,__int64 b,__int64 &d,__int64 &x,__int64 &y)
{
	if(b==0)
	{
		x=1;y=0;
		d=a;
		return;
	}
	else
	{
		exgcd(b,a%b,d,x,y);
		__int64 temp=x;
		x=y;
		y=temp-(a/b)*y;
	}
}
int solve()
{
	int i;
	__int64 a1,a2,r1,r2;
	scanf("%I64d %I64d",&a1,&r1);
	for(i=2;i<=n;i++)
	{
		__int64 a,b,c,d,x0,y0;
		scanf("%I64d %I64d",&a2,&r2);
		a=a1,b=a2,c=r2-r1;
		exgcd(a,b,d,x0,y0);
		if(c%d!=0)
			have=0;
		int t=b/d;
		x0=(x0*(c/d)%t+t)%t;
		r1=a1*x0+r1;
		a1=a1*(a2/d);
	}
	if(!have)
		r1=-1;
	return r1;//返回值即为方程组小于m的非负整数解,返回-1,说明无解
}
//若要求方程组在[0,N](注意[]是闭区间的意思,而非最小公倍数)范围内的解,求法:
//设方程组在m范围内的解为a,则a+m*x(x是自然数)都是方程组的解,求出所有a+m*x(满足a+m*x<=N)即可

//解A的x次方=B(mod C)(0<=x<C)的解
#include<iostream>
#include<cmath>
#define maxn 65535
using namespace std;
struct hash
{
	int a,b,next;
}Hash[maxn*2];
int flg[maxn+66];
int top,idx;
void ins(int a,int b)
{
	int k=b&maxn;
	if(flg[k]!=idx)
	{
		flg[k]=idx;
		Hash[k].next=-1;
		Hash[k].a=a;
		Hash[k].b=b;
		return;
	}
	while(Hash[k].next!=-1)
	{
		if(Hash[k].b==b) return;
		k=Hash[k].next;
	}
	Hash[k].next=++top;
	Hash[top].next=-1;
	Hash[top].a=a;
	Hash[top].b=b;
}
int find(int b)
{
	int k=b&maxn;
	if(flg[k]!=idx) return -1;
	while(k!=-1)
	{
		if(Hash[k].b==b)
			return Hash[k].a;
		k=Hash[k].next;
	}
	return -1;
}
int gcd(int a,int b)
{
	return b==0?a:gcd(b,a%b);
}
int exgcd(int a,int b,int &x,int &y)
{
	int t,ret;
	if(!b) {x=1,y=0;return a;}
	ret=exgcd(b,a%b,x,y);
	t=x,x=y,y=t-a/b*y;
	return ret;
}
int Inval(int a,int b,int n)
{
	int x,y,e;
	exgcd(a,n,x,y);
	e=(__int64)x*b%n;
	return e<0?(e+n):e;
}
int pow_mod(__int64 a,int b,int c)
{
	__int64 ret=1%c;a%=c;
	while(b)
	{
		if(b&1)
			ret=ret*a%c;
		a=a*a%c;
		b=b>>1;
	}
	return ret;
}
int BabyStep(int A,int B,int C)
{
	top=maxn;++idx;
	__int64 buf=1%C,D=buf,K;
	int i,d=0,tmp;
	for(i=0;i<=100;buf=buf*A%C,++i)
		if(buf==B)
			return i;
	while((tmp=gcd(A,C))!=1)
	{
		if(B%tmp) return -1;
		++d;
		C/=tmp;
		B/=tmp;
		D=D*A/tmp%C;
	}
	int M=(int)ceil(sqrt((double)C));
	for(buf=1%C,i=0;i<=M;buf=buf*A%C,++i) ins(i,buf);
	for(i=0,K=pow_mod((__int64)A,M,C);i<=M;D=D*K%C,++i)
	{
		tmp=Inval((int)D,B,C);int w;
		if(tmp>=0 && (w=find(tmp))!=-1)
			return i*M+w+d;
	}
	return -1;
}
int main()
{
	/*
	int A,B,C;
	while(scanf("%d %d %d",&C,&A,&B)!=EOF)
	{
		int tmp=BabyStep(A,B,C);
		if(tmp<0)
			printf("No Answer\n");
		else
			printf("%d\n",tmp);
	}
system("pause");
*/
return 0;
}
//A是矩阵,二分计算A的k次方
#include<cstdlib>
#include<iostream>
#define maxn 101
using namespace std;
typedef struct
{
	int m[maxn][maxn];
}Matrax;
Matrax a,per;
int n,M;//M是??
void init()
{
	int i,j;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			scanf("%d",&a.m[i][j]);
			a.m[i][j]%=M;
			per.m[i][j]=(i==j);
		}
	return;
}
Matrax multi(Matrax a,Matrax b)
{
	Matrax c;
	int k,i,j;
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
		{
			c.m[i][j]=0;
			for(k=0;k<n;k++)
			{
				c.m[i][j]+=a.m[i][k]*b.m[k][j];
			}
			c.m[i][j]%=M;
		}
	return c;
}
Matrax power(int k)
{
	Matrax c,p,ans=per;
	p=a;
	while(k)
	{
		if(k&1)
		{
			ans=multi(ans,p);
			k--;
		}
		else
		{
			k/=2;
			p=multi(p,p);
		}
	}
	return ans;//ans即为A的k次方矩阵
}

中国剩余定理

/*
若m1,m2,m3,……,mr是两两互素的正整数,则同余方程组
x=a1(mod m1)
x=a2(mod m2)
x=a3(mod m3)
……
x=ar(mod mr)
有模M=m1*m2*m3*……*mr的唯一解,即为中国剩余定理
解这样的方程组:

令Mi=M/mi,因为m1,m2,m3,……,mr两两互素,因此(Mi,mi)=1,即MiPi=1(mod mi),那么容易验证a1M1P1+a2M2P2+……+arMrPr就是方程组的解
*/
//求此类同余方程组最小非负整数解的算法实现
int  China(int r)
{
	int M=1;
	int i,Mi,x0,y0,d,ans=0;
	for(i=1;i<=r;i++)
	{
		M*=m[i];
	}
	for(i=1;i<=r;i++)
	{
		Mi=M/m[i];
		exgcd(Mi,m[i],d,x0,y0);
		ans=(ans+Mi*x0*a[i])%M;
	}
	if(ans<0)
		ans+=M;
	return ans;
}

2 thoughts on “数论的模板

发表评论

电子邮件地址不会被公开。 必填项已用 * 标注

*

您可以使用这些 HTML 标签和属性: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>