【转】数论模板 1.扩展的欧几里德定理 2.求模乘法的逆 3.快速幂求解a^b 4.求解模方程a^x=b(mod n)。 5.中国剩余定理 6.素数 7.欧拉phi函数 8.求三角形两边向量积 9.求凸包 10.三分法
//拓展欧几里得定理，求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
void gcd(int a,int b,int &d,int &x,int &y)
{
if(!b){d=a;x=1;y=0;}
else{gcd(b,a%b,d,y,x);y-=x*(a/b);}
}
用法1：
求ax+by=c的整数解。ax+by=gcd(a,b)=g的一组解为(x0,y0),则ax+by=c的一组解为(x0*c/g,y0*c/g)。当c不是g的倍数时无整数解
若(x1,y1)是ax+by=c的一组解，则其任意整数解为(x1+k*bb,y1-k*aa)，其中aa=a/gcd(a,b),bb=b/gcd(bb)，k为任意整数
2.求模乘法的逆
//求得a在模n条件下的逆
int inv(int a,int n)
{
int d,x,y;
gcd(a,n,d,x,y);
return d==1?(x+n)%n:-1;
}
或则：
v=pow_mod(a,n-m-1,n);//n为素数，pow_mod(a,n-1,n)=1,费马小定理。所以a^m*a^(n-m-1)=a^(n-1)=1(mod n).
3.快速幂求解a^b
//快速幂求a^b
int pow_mod(int a,int b)
{
int s=1;
while(b)
{
if(b&1)
s=(s*a)%mod;
a=(a*a)%mod;
b=b>>1;
}
return s;
}
4.求解模方程a^x=b(mod n)。
1.n为素数。无解返回-1
//求解模方程a^x=b(mod n)。n为素数，无解返回-1
int log_mod (int a,int b,int n)
{
int m,v,e=1,i;
m=ceil(sqrt(n+0.5));
v=inv(pow_mod(a,m),n);
map<int,int>x;
x[1]=0;
for(i=1;i<m;i++)
{
e=(e*a)%n;
if(!x.count(e))x[e]=i;
}
for(i=0;i<m;i++)
{
if(x.count(b))return i*m+x[b];
b=(b*v)%n;
}
return -1;
}
2.n不是素数。
//hdu 2815 Mod Tree
#include <cstdio>
#include <cstring>
#include <cmath>
#include <map>
#include <iostream>
#include <algorithm>
using namespace std;
#define LL __int64
LL gcd(LL a,LL b)
{
return b==0?a:gcd(b,a%b);
}
//拓展欧几里得定理，求ax+by=gcd(a,b)的一组解(x,y),d=gcd(a,b)
void gcd_mod(LL a,LL b,LL &d,LL &x,LL &y)
{
if(!b){d=a;x=1;y=0;}
else{gcd_mod(b,a%b,d,y,x);y-=x*(a/b);}
}
//求解模方程d*a^(x-c)=b(mod n)。d,a和n互质，无解返回-1
LL log_mod (LL a,LL b,LL n,LL c,LL d)
{
LL m,v,e=1,i,x,y,dd;
m=ceil(sqrt(n+0.5)); //x=i*m+j
map<LL,LL>f;
f[1]=m;
for(i=1;i<m;i++) //建哈希表，保存a^0,a^1,...,a^m-1
{
e=(e*a)%n;
if(!f[e])f[e]=i;
}
e=(e*a)%n;//e=a^m
for(i=0;i<m;i++)//每次增加m次方，遍历所有1<=f<=n
{
gcd_mod(d,n,dd,x,y);//d*x+n*y=1-->(d*x)%n=1-->d*(x*b)%n==b
x=(x*b%n+n)%n;
if(f[x])
{
LL num=f[x];
f.clear();//需要清空，不然会爆内存
return c+i*m+(num==m?0:num);
}
d=(d*e)%n;
}
return -1;
}
int main()
{
LL a,b,n;
while(scanf("%I64d%I64d%I64d",&a,&n,&b)!=EOF)
{
if(b>=n)
{
printf("Orz,I can’t find D!
");
continue;
}
if(b==0)
{
printf("0
");
continue;
}
LL ans=0,c=0,d=1,t;
while((t=gcd(a,n))!=1)
{
if(b%t){ans=-1;break;}
c++;
n=n/t;
b=b/t;
d=d*a/t%n;
if(d==b){ans=c;break;}//特判下是否成立。
}
if(ans!=0)
{
if(ans==-1){printf("Orz,I can’t find D!
");}
else printf("%I64d
",ans);
}
else
{
ans=log_mod(a,b,n,c,d);
if(ans==-1)printf("Orz,I can’t find D!
");
else printf("%I64d
",ans);
}
}
return 0;
}
/*
求解模方程a^x=b(mod n)，n不为素数。拓展Baby Step Giant Step
模板题。
方法：
初始d=1,c=0,i=0;
1.令g=gcd(a,n),若g==1则执行下一步。否则由于a^x=k*n+b;(k为某一整数),则(a/g)*a^k=k*(n/g)+b/g,(b/g为整除，若不成立则无解)
令n=n/g，d=d*a/g，b=b/g,c++则d*a^(x-c)=b(mod n),接着重复1步骤。
2.通过1步骤后，保证了a和d都与n互质，方程转换为d*a^(x-c)=b(mod n)。由于a和n互质，所以由欧拉定理a^phi(n)==1(mod n),(a,n互质)
可知，phi(n)<=n,a^0==1(mod n),所以构成循环，且循环节不大于n。从而推出如果存在解，则必定1<=x<n。(在此基础上我们就可以用
Baby Step Giant Step方法了)
3.令m=ceil(sqrt(n)),则m*m>=n。用哈希表存储a^0,a^1,..,a^(m-1)，接着判断1~m*m-1中是否存在解。
4.为了减少时间，所以用哈希表缩减复杂度。分成m次遍历，每次遍历a^m长度。由于a和d都与n互质，所以gcd(d,n)=1，
所以用拓展的欧几里德定理求得d*x+n*y=gcd(d,n)=1,的一组整数解(x,y)。则d*x+n*y=1-->d*x%n=(d*x+n*y)%n=1-->d*(x*b)%n=((d*x)%n*b%n)%n=b。
所以若x*b在哈希表中存在，值为k，则a^k*d=b(mod n),答案就是ans=k+c+i*m。如果不存在，则令d=d*a^m,i++后遍历下一个a^m，直到遍历a^0到a^(m-1)还未找到，则说明不解并退出。
*/
5.中国剩余定理
互质情况
//中国剩余定理，求得M%A=a,M%B=b,...中的M，其中A,B,C...互质
int china(int a[])
{
int i,j,k,d,ans=0,x,y,M=1;
for(i=0;i<n;i++)
M=M*x[i];
for(i=0;i<n;i++)
{
m=M/x[i];
gcd(x[i],m,d,x,y);
ans=(ans+y*m*a[i])%M;
}
return (ans+M)%M;
}
不互质的情况
int China(int n)
{
int m1,r1,m2,r2,flag=0,i,d,x,y,c,t;
scanf("%d%d",&m1,&r1);
flag=0;
for(i=1;i<n;i++)
{
scanf("%d%d",&m2,&r2);
if(flag)continue;
gcd(m1,m2,d,x,y);//d=gcd(m1,m2);x*m1+y*m2=d;
c=r2-r1;
if(c%d)//对于方程m1*x+m2*y=c，如果c不是d的倍数就无整数解
{
flag=1;
continue;
}
t=m2/d;//对于方程m1x+m2y=c=r2-r1,若(x0,y0)是一组整数解,那么(x0+k*m2/d,y0-k*m1/d)也是一组整数解(k为任意整数)
//其中x0=x*c/d,y0=x*c/d;
x=(c/d*x%t+t)%t;//保证x0是正数，因为x+k*t是解，(x%t+t)%t也必定是正数解(必定存在某个k使得(x%t+t)%t=x+k*t)
r1=m1*x+r1;//新求的r1就是前i组的解，Mi=m1*x+M(i-1)=r2-m2*y(m1为前i个m的最小公倍数);对m2取余时，余数为r2；
//对以前的m取余时，Mi%m=m1*x%m+M(i-1)%m=M(i-1)%m=r
m1=m1*m2/d;
}
if(flag)return -1;
else return r1;
}
6.素数
普通筛选法
const int maxn=10000001;
const int maxc=700000;
int vis[maxn];//vis[i]=0时，i为素数或者1，否则为合数
int prime[maxc];
//筛素数
void sieve(int n)
{
int m=(int)sqrt(n+0.5);//避免浮点误差
memset(vis,0,sizeof(vis));
for(int i=2;i<=m;i++)if(!vis[i])
for(int j=i*i;j<=n;j+=i)vis[j]=1;
}
//生成素数表,存在prime数组中，返回素数个数
int primes(int n)
{
sieve(n);
int t=0;
for(int i=2;i<=n;i++)if(!vis[i])
prime[t++]=i;
return t;
}
线性筛法
：在O（n）时间复杂度内找出所有的素数
void init()//预处理，找出所有1e7以内的素数，以减少查找1e14范围数的因子的时间
{ //现行筛素数的方法，时间复杂度为O(n)
memset(check,false,sizeof(check));
int i,j;
tot=0;
for(i=2;i<=1e7;i++)
{
if(!check[i])prime[tot++]=i;
for(j=0;j<tot;j++)
{
if(i*prime[j]>1e7)break;
check[i*prime[j]]=true;
if(i%prime[j]==0)break;
}
}
//printf("%d
",tot);
//for(i=0;i<20;i++)
// printf("prime[%d]:%d
",i,prime[i]);
}
7.欧拉phi函数
，求不超过n且与n互质的正整数个数
int euler_phi(int n)//求单个欧拉函数
{
int m=(int)sqrt(n+0.5);
int i,ans=n;
for(i=2;i<=m;i++)
if(n%i==0)
{
ans=ans/i*(i-1);
while(n%i==0)n/=i;
}
if(n>1)ans=ans/n*(n-1);
return ans;
}
int phi[maxn];
void euler_phi()
{
int i,j,k;
//欧拉函数，phi[i]表示不超过i的与i互质的整数个数
for(i=2;i<maxn;i++)phi[i]=0;
phi[1]=1;
for(i=2;i<maxn;i++)
if(!phi[i])
for(j=i;j<=maxn;j+=i){
if(!phi[j])phi[j]=j;
phi[j]=phi[j]/i*(i-1);
}
}
8.求三角形两边向量积
，面积S=1/2*|a×b|
int cross(node a,node b,node c)//向量积
{
return (a.x-c.x)*(b.y-c.y)-(b.x-c.x)*(a.y-c.y);//向量ac×bc
}
9.求凸包
，res数组存凸包上的点。
int convex(int n)
{
sort(e,e+n,cmp);//根据x从小到大排序，x相等时根据y从小到大排序
int m=0,i,j,k;
//求得下凸包，逆时针
for(i=0;i<n;i++)
{
while(m>1&&cross(res[m-1],e[i],res[m-2])<=0)m--;
res[m++]=e[i];
}
k=m;
//求得上凸包
for(i=n-2;i>=0;i--)
{
while(m>k&&cross(res[m-1],e[i],res[m-2])<=0)m--;
res[m++]=e[i];
}
if(n>1)m--;//起始点重复。
return m;//返回凸包边界结点数。
}
10.三分法
for(i=0;i<100;i++)
{
m1=l+(r-l)/3;
m2=r-(r-l)/3;
if(find(m1)<find(m2)) r=m2;
else l=m1;
}