《算法导论》——数论

1 . 欧几里得算法（递归法球两个数的最大公约数）

　　算法比较简单就直接贴代码了：

　　　　int gcd(int a , int b){ 　　　　　　

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

　　　　}

　　在这个算法的基础上可以该进，gcd(a*n , b*n)=n*gcd(a , b);这个定理给我们提供了思路：可不可以在算法中加如判断减少递归的次数呢？当然可以。

int gac(int m,int n){

　　 if(m==0 || n==0) return (m==0) ? n : m;

　　if(m<n)

　　　　return gac(n,m); //m<n时返回gac（n，m）的最大公约数

　　else{

　　　　　　if(m%2 !=1){

　　　　　　　　if(n%2 !=1) //两个数都是k的倍数时函数返回gac（m/k,n/k）的最大公约数

　　　　　　　　　　return 2*gac(m/2,n/2);

　　　　　　　　else

　　　　　　　　　　return gac(m/2,n); //两个数中若只有一个数是k的倍数则返回gac（m/k,n)的最大公约数

　　　　 　　}

　　　　 　　else{

　　　　　　　　if(n%2!=1) return gac(m,n/2); //两个数中若只有一个数是k的倍数则返回gac（m,n/k)的最大公约数

　　　　　　　　else return gac(n,m-n); //两个数都不是k的倍数则返回两个数中小的数和两个数的差的最大公约数

　　　　 　　　　}

　　　　}

}

2 . 扩展形式的欧几里得算法，再求两个数的最大公约数的时候，还有可以利用的信息，假设有这样的一个方程：d=a*x+b*y ; 怎么求？　观察一下，不难发现可以令 d=gcd(a , b) ， 由于gcd(a , b)=gcd(b , a%b )，令d’=gcd（b ，a%b），那 d’=gcd（b ，a%b）=b*x’+a%b*y’，由于a%b=a -（a/b）*b ，于是有等式 b*x’+(a -（a/b）*b)*y’=d=a*x+b*y成立，在合并同类项：d=a*y'+b*（x’-（a/b）*y'），选择x=y'，y=x’-（a/b）*y'；就是方程的一组解。

　　struct equation{
int d=0;
int x=0;
int y=0;
};

equation extend_euclid(int a , int b){ //利用欧几里得算法求解二元方程的一组解
static equation result; //gcd(a , b)=gcd(b , a % b)
if(b == 0){ //
result.d=a;
result.x=1;
result.y=0;
return result;
}else {
equation temp=extend_euclid(b , a%b);
result.d=temp.d;
result.x=temp.y;
result.y=temp.x - (a / b) * temp.y;
return result;
}
}

3 . 模线性方程：a*x=b（mod n）

　　a.根据《算法导论》中数论的讨论中对这类方程的解法，第一步就是确定该方程有没有解，另d=gcd（a ， n），当d整除b的时候该方程才会有解，并且其中一个特解 x0=x’（b/d） % n ；其中x’为gcd（a ， n）=a*x+n*y的解；

　　b.假设方程有解，并且x'是该方程的一个解，那么该方程有d个解，并且都可以用x'表示出来，x[i]=x'+i*(n/d);其中i从0增长到d-1；

void modular_linear_equation_solver(int a ,int b ,int n){

　　equation result=extend_euclid(a , n) ;

if(b % result.d == 0){

　　　　 int x=(result.x*(b/result.d)) % n;

　　　　 for(int i=0 ; i<result.d ; ++i)

　　　　　　printf("%d ",(x+i*(n/result.d))%n);

　　 }

　　else

　　　　printf("no solutions\n");

}

4. 反复平方法求元素的幂：　

　int repede_pow(int a , int b){
　　 int result=1;
　　for(int i = 0 ; i < 64 ; ++i ){
　　 result *= result ; //反复平方
　　if(b & 1 << i) //考虑到b为奇数的情况
　　 result *= a ;
　　 }
　　 return result ;
　}

5. 反复平方法幂的模运算：

long long int modular_exponentiation(long long a , long long b , long long n){ //反复平方法求元素幂

　　 long long result=1;

　　long long c=0;

　　 long long temp = n ;

　　int bit_size = 0 ;

while(temp){

　　　　++bit_size ;

　　　　 temp >>= 1 ;

　 }

for(int i = bit_size ; i >= 0 ; --i){

　　　　 c <<= 1 ; //通过c的成倍增长，达到反复平方

　　　　 result *= result ;

　　　　 result %= n ;

　　　　 if(n << i & 1){ //考虑到奇数的情况

　　　　　　　　c += 1 ;

　　　　　　　 result *= a ;

　　　　　　　 result %= n ;

　　　　 }

　　}

　　return result ;

}

6 . 好了下面说一下素数判定的问题

　　a.暴力测试；

　　b.miller_rabin测试法；

　　暴力测试方法不用多说，说一下miller_rabin测试方法吧，基于前面的快速幂模运算，和乘积模运算，我们可以多次测试多次测试一个数是不是和数判定这个数是不是素数。

　　可以编程实现：

#include<stdio.h>

long long multi(long long a , long long b , long long n){

　　　　long long result=0;

　　　　 a = a % n;

　　　　while(b){

　　　　　　if(b & 1){

　　　　　　　　result = result + a ;

　　　　　　　　if(result >= n) result -= n;

　　　　　　}

　　　　　 a <<= 1;

　　　　　　 if(a > n) a -= n;

　　　　　　 b >>= 1;

　　　　}

　　　　return result;

}

long long pow_mod(long long a , long long b , long long n){

　　　　if(b == 0)

　　　　　　return 1;

　　　　 if(b == 1)

　　　　　　 return a%n;

　　 long long answer=1;

　　　　while(b){

　　　　　　 if(b & 1) answer = multi(answer , a , n); //answer = answer * a % n;

　　　　　　 b >>= 1; //a = a * a % n;

　　　　　　a = multi(a , a , n);

　　　　}

　　　　return answer;

}

bool witness(long long a , long long n){

　　　　 if(n == 2 || n == a) return 0;

　　　　 if((n & 1) == 0) return 1;

　　　long long u=n-1;

　　　　 while(!(u & 1)) u >>= 1;

　　 long long y,x=pow_mod(a , u , n);

　　　　 while(u != n-1){

　　　　　　　　u <<= 1;

　　　　　　　 y = multi(x , x , n);

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

　　　　　　　　x=y;

　　　　}

　　　　if(y != 1) return 1;

　　　　return 0;

}

bool miller_rabin(long long n){

　　　　 if(n < 2) return 0;

　　　　long long test[]={2 , 3 , 7};

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

　　　　　　　　if(witness(test[i] , n)) return 0;

　　　　}

　　　　 return 1;

}

int main(){

　　　　printf("%lld",pow_mod(7,560,561));

　　　　printf("\n\n\n");

　　　　long long date;

　　　　while(scanf("%lld",&date)){

　　　　　　　　if(miller_rabin(date))

　　　　　　　　　　printf("%lldyes\n\n",date);

　　　　}

}

c.还有一种就是打表的方法：

//写得太好了，每一个函数都是一个模版呀~~~~

#include<cstdio>

#include<cstdlib>

#define gcc 10007

#define MAX ((INT)1<<63)-1

typedef unsigned long long INT;

INT p[10]={2,3,5,7,11,13,17,19,23,29};

inline INT gcd(INT a,INT b) {

　　INT m=1;

　　 if(!b)

　　　　 return a;

　　while(m) {

　　　　m=a%b;

　　　　a=b;

　　　　b=m;

　　}

　　return a;

}

//计算a*b%n

inline INT multi_mod(INT a,INT b,INT mod) {

　　INT sum=0;

　　while(b) {

　　　　if(b&1) sum=(sum+a)%mod;

　　　　a=(a+a)%mod;

　　　　b>>=1;

　　}

　　return sum;

}

//计算a^b%n;

inline INT quickmod(INT a,INT b,INT mod) {

　　　　INT sum=1;

　　　　while(b) {

　　　　　　 if(b&1) sum=multi_mod(sum,a,mod);

　　　　　　a=multi_mod(a,a,mod);

　　　　　　b>>=1;

　　　　 }

　　　　return sum;

}

bool miller_rabin(INT n) {

　　 int i,j,k=0;

　　INT u,m,buf; //将n分解为m*2^k

　　if(n==2)

　　　　return true;

　　if(n<2||!(n&1))

　　　　 return false;

　　 m=n-1;

　　while(!(m&1)) k++,m>>=1;

　　 for(i=0;i<9;i++) {

　　　　if(p[i]>=n)

　　　　　　return true;

　　　　 u=quickmod(p[i],m,n);

　　　　if(u==1)

　　　　　　continue;

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

　　　　　　buf=multi_mod(u,u,n);

　　　　　　 if(buf==1&&u!=1&&u!=n-1)

　　　　　　　　 return false;

　　　　　　 u=buf;

　　　　} //如果p[i]^(n-1)%n!=1那么n为合数

　　　　if(u-1)

　　　　 return false;

　　 }

　　return true;

}

//寻找n的一个因子，该因子并不一定是最小的，所以下面要二分查找最小的那个因子

INT pollard(INT n) {

　　　　 INT i=1;

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

　　　　 INT y=x;

　　　　 INT k=2;

　　　　INT d;

　　　　do {

　　　　　　 i++;

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

　　　　　　if(d>1&&d<n)

　　　　　　return d;

　　　　　　 if(i==k)

　　　　　　　　 y=x,k*=2;

　　　　　　x=(multi_mod(x,x,n)+n-gcc)%n;

　　　　}while(y!=x);

　　　　return n;

}

INT pollard_min(INT n) {

　　　　 INT p,a,b=MAX;

　　　　 if(n==1) return MAX;

　　　　if(miller_rabin(n)) return n;

　　　　p=pollard(n);

　　　　a=pollard_min(p);//二分查找

　　　　INT y=n/p;

　　　　b=pollard_min(y);

　　　　return a<b?a:b;

}

7.中国剩余定理求模线性方程组：