girepository/cmph/miller_rabin.c

   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