mprime.cpp

   1 // Bignum prime test (returns 1 if prime, 0 if not)
   2
   3 // Uses Algorithm P (probabilistic primality test) from p. 395 of
   4 // "The Art of Computer Programming, Volume 2" by Donald E. Knuth.
   5
   6 #include "stdafx.h"
   7 #include "defs.h"
   8
   9 static int mprimef(unsigned int *, unsigned int *, int);
  10
  11 int
  12 mprime(unsigned int *n)
  13 {
  14         int i, k;
  15         unsigned int *q;
  16
  17         // 1?
  18
  19         if (MLENGTH(n) == 1 && n[0] == 1)
  20                 return 0;
  21
  22         // 2?
  23
  24         if (MLENGTH(n) == 1 && n[0] == 2)
  25                 return 1;
  26
  27         // even?
  28
  29         if ((n[0] & 1) == 0)
  30                 return 0;
  31
  32         // n = 1 + (2 ^ k) q
  33
  34         q = mcopy(n);
  35
  36         k = 0;
  37         do {
  38                 mshiftright(q);
  39                 k++;
  40         } while ((q[0] & 1) == 0);
  41
  42         // try 25 times
  43
  44         for (i = 0; i < 25; i++)
  45                 if (mprimef(n, q, k) == 0)
  46                         break;
  47
  48         mfree(q);
  49
  50         if (i < 25)
  51                 return 0;
  52         else
  53                 return 1;
  54 }
  55
  56 //-----------------------------------------------------------------------------
  57 //
  58 //      This is the actual implementation of Algorithm P.
  59 //
  60 //      Input:          n               The number in question.
  61 //
  62 //                      q               n = 1 + (2 ^ k) q
  63 //
  64 //                      k
  65 //
  66 //      Output:         1               when n is probably prime
  67 //
  68 //                      0               when n is definitely not prime
  69 //
  70 //-----------------------------------------------------------------------------
  71
  72 static int
  73 mprimef(unsigned int *n, unsigned int *q, int k)
  74 {
  75         int i, j;
  76         unsigned int *t, *x, *y;
  77
  78         // generate x
  79
  80         t = mcopy(n);
  81
  82         while (1) {
  83                 for (i = 0; i < MLENGTH(t); i++)
  84                         t[i] = rand();
  85                 x = mmod(t, n);
  86                 if (!MZERO(x) && !MEQUAL(x, 1))
  87                         break;
  88                 mfree(x);
  89         }
  90
  91         mfree(t);
  92
  93         // exponentiate
  94
  95         y = mmodpow(x, q, n);
  96
  97         // done?
  98
  99         if (MEQUAL(y, 1)) {
 100                 mfree(x);
 101                 mfree(y);
 102                 return 1;
 103         }
 104
 105         j = 0;
 106
 107         while (1) {
 108
 109                 // y = n - 1?
 110
 111                 t = msub(n, y);
 112
 113                 if (MEQUAL(t, 1)) {
 114                         mfree(t);
 115                         mfree(x);
 116                         mfree(y);
 117                         return 1;
 118                 }
 119
 120                 mfree(t);
 121
 122                 if (++j == k) {
 123                         mfree(x);
 124                         mfree(y);
 125                         return 0;
 126                 }
 127
 128                 // y = (y ^ 2) mod n
 129
 130                 t = mmul(y, y);
 131                 mfree(y);
 132                 y = mmod(t, n);
 133                 mfree(t);
 134
 135                 // y = 1?
 136
 137                 if (MEQUAL(y, 1)) {
 138                         mfree(x);
 139                         mfree(y);
 140                         return 0;
 141                 }
 142         }
 143 }