几个初等数论算法

2009-02-05 15:14
//注：这个版本是2009.4.6更新的，原来的算法存在一些小问题。

拓展的欧几里得算法
#include<iostream>
using namespace std;
int extend_Euclid(int& x,int& y,int a,int b);
int main(void){
    int x,y,a,b,d;
   
    cin>>a>>b;
    d = extend_Euclid(x,y,a,b);
   
    cout<<d<<"="<<x<<"*"<<a<<"+"<<y<<"*"<<b<<endl;
   
    system("pause");
    return 0;
}
int extend_Euclid(int& x,int& y,int a,int b){
    if(b==0){
       x = 1;
       y = 0;
       return a;
    }
   
    int d = extend_Euclid(x,y,b,a%b);
    int tx = x;
    x = y;
    y = tx - (a/b)*y;
   
    return d;
}
解线性同余方程
#include<iostream>
using namespace std;
int Extended_Euclid(int a,int b,int& x,int& y);
void solve(int a,int b,int n);
int main(void){
    int a,b,n;
   
    cout<<"Input a,b,n:";
    cin>>a>>b>>n;
   
    solve(a,b,n);
   
    system("pause");
    return 0;
}
void solve(int a,int b,int n){
     int x,y;
     int d = Extended_Euclid(a,n,x,y);
    
     if(b%d!=0){
        cout<<"No solutions!";
        return;
     }
    
     int x0 = x*(b/d) % n;//注：由于C的奇怪的mod运算导致解可能出现负值

                          //。。。。。。虽然也是正确的，但总是不太好

     for(int i=0;i<d;i++)
         cout<<"x"<<i<<"="<<(x0+i*(n/d) )<<endl;
}
int Extended_Euclid(int a,int b,int& x,int& y){
    if(b==0){
       x = 1;
       y = 0;
       return a;
    }
   
    int d = Extended_Euclid(b,a%b,x,y);
    int tx = x;
    x = y;
    y = tx - (a/b)*y;
   
    return d;
}
用中国剩余定理解同余方程组
#include<iostream>
using namespace std;
int Extended_Euclid(int a,int b,int& x,int& y);
void solve(int a[],int n[],int k);
int main(void){
    int k,a[100],n[100];
   
    cout<<"Input the number of formulas:";
    cin>>k;
    cout<<"Notice:n[i] must be a prime number!"<<endl;
    for(int i=0;i<k;i++){
        cout<<"Input a["<<i<<"],n["<<i<<"]:";
        cin>>a[i]>>n[i];
    }
   
    solve(a,n,k);
   
    system("pause");
    return 0;
}
void solve(int a[],int n[],int k){
     int d,x,y,r;
    
     d = ExtendedEuclid(a,n,x,y);
    
     if(b%d!=0){
        cout<<"No Answer!";
        return;
     }
    
     r = ( x*(b/d) )%n;
     if(r<0) r+= n;
     for(int i=0;i<d;i++)
         cout<<( r + i*(n/d) )%n<<"\t";
}
int Extended_Euclid(int a,int b,int& x,int& y){
    if(b==0){
       x = 1;
       y = 0;
       return a;
    }
    
    int d = Extended_Euclid(b,a%b,x,y);
    
    int tx = x;
    x = y;
    y = tx-(a/b);
   
    return d;
}
Miller_Rabin:
#include<iostream>
#include<cstdlib>
#include<ctime>
using namespace std;
const int s = 50;
typedef unsigned long long ULL;
ULL ModularExponentiation(ULL a,ULL b,ULL n);
bool MillerRabin(ULL n);
int main(void){
    ULL n;
   
    cin>>n;
    if(MillerRabin(n) ) cout<<"n should be a prime!"<<endl;
    else cout<<"n is not a prime!"<<endl;
   
    system("pause");
    return 0;
}
ULL ModularExponentiation(ULL a,ULL b,ULL n){//有溢出的危险

    ULL r = 1,t = a;
   
    while(b>0){
          if(b%2) r = (r*t)%n;
          b /= 2;
          t = (t*t)%n;
    }
   
    return r;
}
   
bool MillerRabin(ULL n){
     ULL t,u,a,x,tx;
     int i,j;
    
     t = 0;
     u = n - 1;
     while(u%2==0){
           t++;
           u /= 2;
     }
    
     srand(time(NULL) );
     for(i=0;i<s;i++){
         a = rand()%(n-1)+1;//这里如果是a = rand()%n;，则有可能不小心取到0

         x = ModularExponentiation(a,u,n);
         for(j=0;j<t;j++){
             tx = x;
             x = (x*x)%n;
             if(x==1&&tx!=1&&tx!=n-1) return false;
         }
        
         if(x!=1) return false;
     }
    
     return true;
}