ACM数论模板

w2w原创

1.欧几里得定理

int gcd(int a, int b){

    return b==0?a:gcd(b, a%b);

}

2.扩展欧几里得（求ax+by = gcd(a,b)的特解)

void e_gcd(LL a, LL b, LL &d, LL &x, LL &y){

    if(b==0){

        x = 1; y = 0; d = a;

    }

    else{

        e_gcd(b, a%b, d, y, x);

        y-= x*(a/b);

    }

}

3.中国剩余定理

同余方程组

x ≡ a1(mod m1)

x ≡ a2(mod m2)

... ...

x ≡ ak(mod mk)

方程组所有的解的集合就是：

x1 = N1*a1 + N2*a2 + ... + Nk*ak

 

其中 Ni mod mi = 1，Ni = ai * ti , 可用欧几里得扩展定理求 ti. 其中M = m1*m2*m3····*mn;

   //互质版
　　#include <iostream> 

    using namespace std; 

    //参数可为负数的扩展欧几里德定理 

    void exOJLD(int a, int b, int &x, int &y){ 

        //根据欧几里德定理 

        if(b == 0){//任意数与0的最大公约数为其本身。 

            x = 1; 

            y = 0; 

        }else{ 

            int x1, y1; 

            exOJLD(b, a%b, x1, y1); 

            if(a*b < 0){//异号取反 

                x = - y1; 

                y = a/b*y1 - x1; 

            }else{//同号 

                x = y1; 

                y = x1 - a/b* y1; 

            } 

        } 

    } 

    //剩余定理 

    int calSYDL(int a[], int m[], int k){ 

        int N[k];//这个可以删除 

        int mm = 1;//最小公倍数 

        int result = 0; 

        for(int i = 0; i < k; i++){ 

            mm *= m[i]; 

        } 

        for(int j = 0; j < k; j++){ 

            int L, J; 

            exOJLD(mm/m[j], -m[j], L, J); 

            N[j] = m[j] * J + 1;//1 

            N[j] = mm/m[j] * L;//2 【注】1和2这两个值应该是相等的。 

            result += N[j]*a[j]; 

        } 

        return (result % mm + mm) % mm;//落在(0, mm)之间，这么写是为了防止result初始为负数，本例中不可能为负可以直接 写成：return result%mm;即可。 

    } 

     

     

    int main(){ 

        int a[3] = {2, 3, 2}; 

        int m[3] = {3, 5, 7};     

        cout<<"结果:"<<calSYDL(a, m, 3)<<endl; 

    } 

　　//不互质版
　　    /**
    中国剩余定理（不互质）
    */  
    #include <iostream>  
    #include <cstdio>  
    #include <cstring>  
    using namespace std;  
    typedef long long LL;  
    LL Mod;  
      
    LL gcd(LL a, LL b)  
    {  
        if(b==0)  
            return a;  
        return gcd(b,a%b);  
    }  
      
    LL Extend_Euclid(LL a, LL b, LL&x, LL& y)  
    {  
        if(b==0)  
        {  
            x=1,y=0;  
            return a;  
        }  
        LL d = Extend_Euclid(b,a%b,x,y);  
        LL t = x;  
        x = y;  
        y = t - a/b*y;  
        return d;  
    }  
      
    //a在模n乘法下的逆元，没有则返回-1  
    LL inv(LL a, LL n)  
    {  
        LL x,y;  
        LL t = Extend_Euclid(a,n,x,y);  
        if(t != 1)  
            return -1;  
        return (x%n+n)%n;  
    }  
      
    //将两个方程合并为一个  
    bool merge(LL a1, LL n1, LL a2, LL n2, LL& a3, LL& n3)  
    {  
        LL d = gcd(n1,n2);  
        LL c = a2-a1;  
        if(c%d)  
            return false;  
        c = (c%n2+n2)%n2;  
        c /= d;  
        n1 /= d;  
        n2 /= d;  
        c *= inv(n1,n2);  
        c %= n2;  
        c *= n1*d;  
        c += a1;  
        n3 = n1*n2*d;  
        a3 = (c%n3+n3)%n3;  
        return true;  
    }  
      
    //求模线性方程组x=ai(mod ni),ni可以不互质  
    LL China_Reminder2(int len, LL* a, LL* n)  
    {  
        LL a1=a[0],n1=n[0];  
        LL a2,n2;  
        for(int i = 1; i < len; i++)  
        {  
            LL aa,nn;  
            a2 = a[i],n2=n[i];  
            if(!merge(a1,n1,a2,n2,aa,nn))  
                return -1;  
            a1 = aa;  
            n1 = nn;  
        }  
        Mod = n1;  
        return (a1%n1+n1)%n1;  
    }  
    LL a[1000],b[1000];  
    int main()  
    {  
        int i;  
        int k;  
        while(scanf("%d",&k)!=EOF)  
        {  
            for(i = 0; i < k; i++)  
                scanf("%I64d %I64d",&a[i],&b[i]);  
            printf("%I64d ",China_Reminder2(k,b,a));  
        }  
        return 0;  
    } 

 

 

4.欧拉函数（求一个数前面的所有与这个数互质的数的个数）

  Euler函数表达通式：euler(x)=x(1-1/p1)(1-1/p2)(1-1/p3)(1-1/p4)…(1-1/pn),其中p1,p2……pn为x的所有素因数，x是不为0的整数。euler(1)=1（唯一和1互质的数就是1本身）。

 

Euler函数有几个性质：

　　1.如果q,p互质，则Euler（p*q） = Euler（p）*Euler（q）；

　　2.如果 a = p^k，则Euler（a） = p^k - p^k-1;

 

    //直接求解欧拉函数 

    int euler(int n){ //返回euler(n)  

         int res=n,a=n; 

         for(int i=2;i*i<=a;i++){ 

             if(a%i==0){ 

                 res=res/i*(i-1);//先进行除法是为了防止中间数据的溢出  

                 while(a%i==0) a/=i; 

             } 

         } 

         if(a>1) res=res/a*(a-1); 

         return res; 

    } 

   

    //线性筛选欧拉函数O(n)用到了一下性质：

    //(1) 若(N%a==0 && (N/a)%a==0) 则有:E(N)=E(N/a)*a;

    //(2) 若(N%a==0 && (N/a)%a!=0) 则有:E(N)=E(N/a)*(a-1);  

    //注意：如果范围过大 可能不适宜开数组来做

    int euler[maxN], vis[maxN], prime[maxN/5], e[maxN], cnt = 0;
    void make_euler(){
        memset(vis, 0, sizeof(vis));
        euler[1] = 1;
        for(int i=2; i<maxN ; ++i){
            if(vis[i] == 0){
                prime[cnt++] = i;
                euler[i] = i-1;
            }
            for(int j=0 ; j<cnt && i*prime[j] < maxN; ++j){
                vis[i*prime[j]] = 1;
                if( i%prime[j] == 0){
                    euler[i*prime[j]] = euler[i] *prime[j];
                    break;
                }
                else euler[i*prime[j]] = euler[i] *(prime[j]-1);
            }
        }
    }
   

 

 

 5.求N以前N的约数个数

   约数个数的性质，对于一个数N，N=p1^a1 + p2^a2 + ... + pn^an。其中p1 ,p2, p3... pn是N的质因数，a1 ,a2, a2,...an为相应的指数，则
                                                           div_num[N]=(p1+1)*(p2+1)*(p3+1)* ... *(pn+1);
结合这个算法的特点，在程序中如下运用：
  对于div_num：

(1)如果i|prime[j] 那么 div_num[i*prime[j]]=div_sum[i]/(e[i]+1)*(e[i]+2)                  //最小素因子次数加1
(2)否则 div_num[i*prime[j]]=div_num[i]*div_num[prime[j]]                                     //满足积性函数条件

  对于e：

(1)如果i|pr[j]  e[i*pr[j]]=e[i]+1; //最小素因子次数加1
(2)否则 e[i*pr[j]]=1;              //pr[j]为1次

 

 

    #include<string.h> 

    #include<iostream>

    #define M 100000 

    using namespace std;

    int prime[M/3],e[M],div_num[M];           // e[i]表示第i个素数因子的个数 

    bool flag[M]; 

    void get_prime() 

    { 

        int i,j,k; 

        memset(flag,false,sizeof(flag)); 

        k=0; 

        for(i=2;i<M;i++){ 

            if(!flag[i]){                             

                prime[k++]=i; 

                e[i]=1; 

                div_num[i]=2;                       //素数的约数个数为2 

            } 

            for(j=0;j<k&&i*prime[j]<M;j++){ 

                flag[i*prime[j]]=true;             

                    if(i%prime[j]==0){ 

                        div_num[i*prime[j]]=div_num[i]/(e[i]+1)*(e[i]+2); 

                        e[i*prime[j]]=e[i]+1; 

                        break; 

                    } 

                    else{ 

                        div_num[i*prime[j]]=div_num[i]*div_num[prime[j]]; 

                        e[i*prime[j]]=1; 

                    } 

            }

        } 

    } 

 

 

6.莫比乌斯函数

　　一个讲得比较清楚的PPT：http://wenku.baidu.com/link?url=UARIPTGHjN78vIzedWT2iwICudBIbsuZ5WMrYwJJjp2P5x7hUvtvSoVKiW7a92GiiF7aCJu1FYid2eB5iM9Wh-hW2Bfd1UfJgrstX7nZnrm

　　线性筛打表莫比乌斯函数：

 

int mob[maxN], vis[maxN], prime[maxN], cnt=0;

void make_mobius(){

    mob[1] = 1;

    memset(vis, 0, sizeof(vis));

    for(int i = 2; i<maxN ; ++i){

        if(!vis[i]){

            mob[i] = -1;

            prime[cnt++] = i;

        }

        for(int j= 0; j<cnt && i*prime[j] < maxN ; ++j){

            vis[i*prime[j]] = 1;

            if(i%prime[j] == 0){

                mob[i*prime[j]] = 0;

                break;

            }

            else mob[i*prime[j]] = -mob[i];

        }

    }

}

 

7.卢卡斯定理

 //卢卡斯定理  C(m, n)%p  

 LL Lucas(LL m, LL n, LL p){

     LL res = 1;

     while(n && m){

         LL n1 = n%p, m1 = m%p;

         //费马小定理求逆元 

         res = res* fac[n1]* qsm(fac[m1]*fac[n1-m1]%p, p-2, p)%p;

         //扩展欧几里得求逆元 

         //res = res* fac[n1]* reverse(fac[m1],p)* reverse(fac[n1-m1],p)%p;

         n /= p;

         m /= p;

     }

     return (res%p + p)%p;

 }

 

 

 

8.卡特兰数

 

 

递推公式

 

9.伯努利数

 

   伯努利数满足条件，且有

    

     那么继续得到

    

    void init_ber(){

        ber[0] = 1;

        for(int i = 1 ; i<maxn; ++i){

            LL ans = 0;

            for(int j = 0 ; j<i ; ++j)

                ans = (ans + comb[i+1][j]*ber[j])%mod;

            ans = -ans*inv(i+1, mod)%mod;

            ber[i] = (ans%mod + mod)%mod;

        }

    }

 

10.乘法逆元

 

因为可能会很大，超过int范围，所以在快速幂时要二分乘法。

逆元打表 ( O（n）)：

 

     inv[1] = 1;  

     for(int i=2;i<N;i++)  

     {  

         if(i >= MOD) break;  

         inv[i] = (MOD - MOD / i) * inv[MOD % i]% MOD;  

     }  

 

 

11. miller rabin算法 pollard rho算法（概率高效判断素数，求因子）

#include <stdlib.h>

#include <string.h>

#include <algorithm>

#include <iostream>

#include <math.h>

 

const int Times=10;

const int N=5500;

 

using namespace std;

typedef long long LL;

 

LL ct,cnt,c,x,y;

LL fac[N],num[N],arr[N];

 

LL gcd(LL a,LL b)

{

    return b? gcd(b,a%b):a;

}

 

LL multi(LL a,LL b,LL m)

{

    LL ans=0;

    while(b)

    {

        if(b&1)

        {

            ans=(ans+a)%m;

            b--;

        }

        b>>=1;

        a=(a+a)%m;

    }

    return ans;

}

 

LL quick_mod(LL a,LL b,LL m)

{

    LL ans=1;

    a%=m;

    while(b)

    {

        if(b&1)

        {

            ans=multi(ans,a,m);

            b--;

        }

        b>>=1;

        a=multi(a,a,m);

    }

    return ans;

}

 

bool Miller_Rabin(LL n)

{

    if(n==2) return true;

    if(n<2||!(n&1)) return false;

    LL a,m=n-1,x,y;

    int k=0;

    while((m&1)==0)

    {

        k++;

        m>>=1;

    }

    for(int i=0; i<Times; i++)

    {

        a=rand()%(n-1)+1;

        x=quick_mod(a,m,n);

        for(int j=0; j<k; j++)

        {

            y=multi(x,x,n);

            if(y==1&&x!=1&&x!=n-1) return false;

            x=y;

        }

        if(y!=1) return false;

    }

    return true;

}

 

LL Pollard_rho(LL n,LL c)

{

    LL x,y,d,i=1,k=2;

    y=x=rand()%(n-1)+1;

    while(true)

    {

        i++;

        x=(multi(x,x,n)+c)%n;

        d=gcd((y-x+n)%n,n);

        if(1<d&&d<n) return d;

        if(y==x) return n;

        if(i==k)

        {

            y=x;

            k<<=1;

        }

    }

}

 

void find(LL n,int c)

{

    if(n==1) return;

    if(Miller_Rabin(n))

    {

        fac[ct++]=n;

        return ;

    }

    LL p=n;

    LL k=c;

    while(p>=n) p=Pollard_rho(p,c--);

    find(p,k);

    find(n/p,k);

}

 

void dfs(LL dept, LL product=1)

{

    if(dept==cnt)

    {

        arr[c++]=product;

        return;

    }

    for(int i=0; i<=num[dept]; i++)

    {

        dfs(dept+1,product);

        product*=fac[dept];

    }

}

 

void Solve(LL n)

{

    ct=0;

    find(n,120);

    sort(fac,fac+ct);

    num[0]=1;

    int k=1;

    for(int i=1; i<ct; i++)

    {

        if(fac[i]==fac[i-1])

            ++num[k-1];

        else

        {

            num[k]=1;

            fac[k++]=fac[i];

        }

    }

    cnt=k;

}

 

const int M=1000005;

bool prime[M];

int p[M];

int k1;

 

void isprime()

{

    k1=0;

    int i,j;

    memset(prime,true,sizeof(prime));

    for(i=2;i<M;i++)

    {

        if(prime[i])

        {

            p[k1++]=i;

            for(j=i+i;j<M;j+=i)

            {

                prime[j]=false;

            }

        }

    }

}

 

int main()

{

    LL n,t,record,tmp;

    isprime();

    while(cin>>n)

    {

        LL ans=-1;

        record=0;

        while(true)

        {

            tmp=1;

            if(Miller_Rabin(n))

            {

                Solve(n-1);

                c=0;

                dfs(0,1);

                sort(arr,arr+c);

                bool flag=0;

                for(int i=0; i<c; i++)

                {

                    if(quick_mod(10,arr[i],n)==1)

                    {

                        tmp=arr[i];

                        break;

                    }

                    if(i==c-2) flag=1;

                }

                if(flag)

                {

                    if(ans<tmp) record=n;

                    cout<<record<<endl;

                    break;

                }

            }

            else

            {

                bool flag=false;

                LL tmp1=(LL)sqrt(n*1.0);

                if(tmp1*tmp1==n&&Miller_Rabin(tmp1))

                {

                    x=tmp1;

                    y=2;

                    flag=true;

                }

                else

                {

                    LL cnt1=0,rea=n;

                    for(int i=0;i<k1;i++)

                    {

                        if(rea%p[i]==0)

                        {

                            x=p[i];

                            while(rea%p[i]==0)

                            {

                                rea/=p[i];

                                cnt1++;

                            }

                            break;

                        }

                    }

                    if(rea==1) flag=true;

                    y=cnt1;

                }

                if(flag)

                {

                    Solve(x-1);

                    c=0;

                    dfs(0,1);

                    sort(arr,arr+c);

                    bool flag=0;

                    for(int i=0; i<c; i++)

                    {

                        if(quick_mod(10,arr[i],x)==1)

                        {

                            tmp=1;

                            for(int j=0; j<y-1; j++)

                                tmp*=x;

                            tmp*=(x-1);

                            break;

                        }

                        if(i==c-2) flag=1;

                    }

                    if(flag)

                    {

                        if(ans<tmp)

                        {

                            ans=tmp;

                            record=n;

                        }

                    }

                }

            }

            n--;

        }

    }

    return 0;

}

 

12.快速乘快速幂

 

LL qsm(LL a, LL n, LL m){

    LL ans = 1;

    while(n>0){

        if(n&1){

            ans = (a*ans)%m;

        }

        a = (a*a)%m;

        n>>=1;

    }

    return ans;

}

long long q_mul( long long a, long long b, long long mod ) //快速计算 (a*b) % mod

{

    long long ans = 0;  // 初始化

    while(b)                //根据b的每一位看加不加当前a

    {

        if(b & 1)           //如果当前位为1

        {

            b--;               //也可不要，方便理解而已

            ans =(ans+ a)%mod;   //ans+=a

        }

        b /= 2;                         //b向前移位

        a = (a + a) % mod;          //更新a

    }

    return ans;

}

 

 

13.博弈

奇异局势： ak =[k（1+√5）/2]，bk= ak + k  （k=0，1，2，…,n 方括号表示取整函 数）。