【转】数论模板

1.扩展的欧几里德定理

//拓展欧几里得定理，求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;
}