1  #include "miller_rabin.h"
       2  
       3  static inline cmph_uint64 int_pow(cmph_uint64 a, cmph_uint64 d, cmph_uint64 n)
       4  {
       5  	cmph_uint64 a_pow = a;
       6  	cmph_uint64 res = 1;
       7  	while(d > 0)
       8  	{
       9  		if((d & 1) == 1)
      10  			res =(((cmph_uint64)res) * a_pow) % n;
      11  		a_pow = (((cmph_uint64)a_pow) * a_pow) % n;
      12  		d /= 2;
      13  	};
      14  	return res;
      15  };
      16  
      17  static inline cmph_uint8 check_witness(cmph_uint64 a_exp_d, cmph_uint64 n, cmph_uint64 s)
      18  {
      19  	cmph_uint64 i;
      20  	cmph_uint64 a_exp = a_exp_d;
      21  	if(a_exp == 1 || a_exp == (n - 1))
      22  		return 1;
      23  	for(i = 1; i < s; i++)
      24  	{
      25  		a_exp = (((cmph_uint64)a_exp) * a_exp) % n;
      26  		if(a_exp == (n - 1))
      27  			return 1;
      28  	};
      29  	return 0;
      30  };
      31  
      32  cmph_uint8 check_primality(cmph_uint64 n)
      33  {
      34  	cmph_uint64 a, d, s, a_exp_d;
      35  	if((n % 2) == 0)
      36  		return 0;
      37  	if((n % 3) == 0)
      38  		return 0;
      39  	if((n % 5) == 0)
      40  		return 0;
      41  	if((n % 7 ) == 0)
      42  		return 0;
      43  	//we decompoe the number n - 1 into 2^s*d
      44  	s = 0;
      45  	d = n - 1;
      46  	do 
      47  	{
      48  		s++;
      49  		d /= 2;
      50  	}while((d % 2) == 0);
      51  
      52  	a = 2;
      53  	a_exp_d = int_pow(a, d, n);
      54  	if(check_witness(a_exp_d, n, s) == 0)
      55  		return 0;
      56  	a = 7;
      57  	a_exp_d = int_pow(a, d, n);
      58  	if(check_witness(a_exp_d, n, s) == 0)
      59  		return 0;
      60  	a = 61;
      61  	a_exp_d = int_pow(a, d, n);
      62  	if(check_witness(a_exp_d, n, s) == 0)
      63  		return 0;
      64  	return 1;
      65  };
      66  
      67