pollard.cpp

   1 // Factor using the Pollard rho method
   2
   3 #include "stdafx.h"
   4 #include "defs.h"
   5
   6 static unsigned int *n;
   7
   8 void
   9 factor_number(void)
  10 {
  11         int h;
  12
  13         save();
  14
  15         p1 = pop();
  16
  17         // 0 or 1?
  18
  19         if (equaln(p1, 0) || equaln(p1, 1) || equaln(p1, -1)) {
  20                 push(p1);
  21                 restore();
  22                 return;
  23         }
  24
  25         n = mcopy(p1->u.q.a);
  26
  27         h = tos;
  28
  29         factor_a();
  30
  31         if (tos - h > 1) {
  32                 list(tos - h);
  33                 push_symbol(MULTIPLY);
  34                 swap();
  35                 cons();
  36         }
  37
  38         restore();
  39 }
  40
  41 // factor using table look-up, then switch to rho method if necessary
  42
  43 // From TAOCP Vol. 2 by Knuth, p. 380 (Algorithm A)
  44
  45 void
  46 factor_a(void)
  47 {
  48         int k;
  49
  50         if (MSIGN(n) == -1) {
  51                 MSIGN(n) = 1;
  52                 push_integer(-1);
  53         }
  54
  55         for (k = 0; k < 10000; k++) {
  56
  57                 try_kth_prime(k);
  58
  59                 // if n is 1 then we're done
  60
  61                 if (MLENGTH(n) == 1 && n[0] == 1) {
  62                         mfree(n);
  63                         return;
  64                 }
  65         }
  66
  67         factor_b();
  68 }
  69
  70 void
  71 try_kth_prime(int k)
  72 {
  73         int count;
  74         unsigned int *d, *q, *r;
  75
  76         d = mint(primetab[k]);
  77
  78         count = 0;
  79
  80         while (1) {
  81
  82                 // if n is 1 then we're done
  83
  84                 if (MLENGTH(n) == 1 && n[0] == 1) {
  85                         if (count)
  86                                 push_factor(d, count);
  87                         else
  88                                 mfree(d);
  89                         return;
  90                 }
  91
  92                 mdivrem(&q, &r, n, d);
  93
  94                 // continue looping while remainder is zero
  95
  96                 if (MLENGTH(r) == 1 && r[0] == 0) {
  97                         count++;
  98                         mfree(r);
  99                         mfree(n);
 100                         n = q;
 101                 } else {
 102                         mfree(r);
 103                         break;
 104                 }
 105         }
 106
 107         if (count)
 108                 push_factor(d, count);
 109
 110         // q = n/d, hence if q < d then n < d^2 so n is prime
 111
 112         if (mcmp(q, d) == -1) {
 113                 push_factor(n, 1);
 114                 n = mint(1);
 115         }
 116
 117         if (count == 0)
 118                 mfree(d);
 119
 120         mfree(q);
 121 }
 122
 123 // From TAOCP Vol. 2 by Knuth, p. 385 (Algorithm B)
 124
 125 int
 126 factor_b(void)
 127 {
 128         int k, l;
 129         unsigned int *g, *one, *t, *x, *xprime;
 130
 131         one = mint(1);
 132         x = mint(5);
 133         xprime = mint(2);
 134
 135         k = 1;
 136         l = 1;
 137
 138         while (1) {
 139
 140                 if (mprime(n)) {
 141                         push_factor(n, 1);
 142                         mfree(one);
 143                         mfree(x);
 144                         mfree(xprime);
 145                         return 0;
 146                 }
 147
 148                 while (1) {
 149
 150                         if (esc_flag) {
 151                                 mfree(one);
 152                                 mfree(n);
 153                                 mfree(x);
 154                                 mfree(xprime);
 155                                 stop("esc");
 156                         }
 157
 158                         // g = gcd(x' - x, n)
 159
 160                         t = msub(xprime, x);
 161                         MSIGN(t) = 1;
 162                         g = mgcd(t, n);
 163                         mfree(t);
 164
 165                         if (MEQUAL(g, 1)) {
 166                                 mfree(g);
 167                                 if (--k == 0) {
 168                                         mfree(xprime);
 169                                         xprime = mcopy(x);
 170                                         l *= 2;
 171                                         k = l;
 172                                 }
 173
 174                                 // x = (x ^ 2 + 1) mod n
 175
 176                                 t = mmul(x, x);
 177                                 mfree(x);
 178                                 x = madd(t, one);
 179                                 mfree(t);
 180                                 t = mmod(x, n);
 181                                 mfree(x);
 182                                 x = t;
 183
 184                                 continue;
 185                         }
 186
 187                         push_factor(g, 1);
 188
 189                         if (mcmp(g, n) == 0) {
 190                                 mfree(one);
 191                                 mfree(n);
 192                                 mfree(x);
 193                                 mfree(xprime);
 194                                 return -1;
 195                         }
 196
 197                         // n = n / g
 198
 199                         t = mdiv(n, g);
 200                         mfree(n);
 201                         n = t;
 202
 203                         // x = x mod n
 204
 205                         t = mmod(x, n);
 206                         mfree(x);
 207                         x = t;
 208
 209                         // xprime = xprime mod n
 210
 211                         t = mmod(xprime, n);
 212                         mfree(xprime);
 213                         xprime = t;
 214
 215                         break;
 216                 }
 217         }
 218 }
 219
 220 void
 221 push_factor(unsigned int *d, int count)
 222 {
 223         p1 = alloc();
 224         p1->k = NUM;
 225         p1->u.q.a = d;
 226         p1->u.q.b = mint(1);
 227         push(p1);
 228         if (count > 1) {
 229                 push_symbol(POWER);
 230                 swap();
 231                 p1 = alloc();
 232                 p1->k = NUM;
 233                 p1->u.q.a = mint(count);
 234                 p1->u.q.b = mint(1);
 235                 push(p1);
 236                 list(3);
 237         }
 238 }