123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239 |
- // Factor using the Pollard rho method
- #include "stdafx.h"
- #include "defs.h"
- static unsigned int *n;
- void
- factor_number(void)
- {
- int h;
- save();
- p1 = pop();
- // 0 or 1?
- if (equaln(p1, 0) || equaln(p1, 1) || equaln(p1, -1)) {
- push(p1);
- restore();
- return;
- }
- n = mcopy(p1->u.q.a);
- h = tos;
- factor_a();
- if (tos - h > 1) {
- list(tos - h);
- push_symbol(MULTIPLY);
- swap();
- cons();
- }
- restore();
- }
- // factor using table look-up, then switch to rho method if necessary
- // From TAOCP Vol. 2 by Knuth, p. 380 (Algorithm A)
- void
- factor_a(void)
- {
- int k;
- if (MSIGN(n) == -1) {
- MSIGN(n) = 1;
- push_integer(-1);
- }
- for (k = 0; k < 10000; k++) {
- try_kth_prime(k);
- // if n is 1 then we're done
- if (MLENGTH(n) == 1 && n[0] == 1) {
- mfree(n);
- return;
- }
- }
- factor_b();
- }
- void
- try_kth_prime(int k)
- {
- int count;
- unsigned int *d, *q, *r;
- d = mint(primetab[k]);
- count = 0;
- while (1) {
- // if n is 1 then we're done
- if (MLENGTH(n) == 1 && n[0] == 1) {
- if (count)
- push_factor(d, count);
- else
- mfree(d);
- return;
- }
- mdivrem(&q, &r, n, d);
- // continue looping while remainder is zero
- if (MLENGTH(r) == 1 && r[0] == 0) {
- count++;
- mfree(r);
- mfree(n);
- n = q;
- } else {
- mfree(r);
- break;
- }
- }
- if (count)
- push_factor(d, count);
- // q = n/d, hence if q < d then n < d^2 so n is prime
- if (mcmp(q, d) == -1) {
- push_factor(n, 1);
- n = mint(1);
- }
- if (count == 0)
- mfree(d);
- mfree(q);
- }
- // From TAOCP Vol. 2 by Knuth, p. 385 (Algorithm B)
- int
- factor_b(void)
- {
- int k, l;
- unsigned int *g, *one, *t, *x, *xprime;
- one = mint(1);
- x = mint(5);
- xprime = mint(2);
- k = 1;
- l = 1;
- while (1) {
- if (mprime(n)) {
- push_factor(n, 1);
- mfree(one);
- mfree(x);
- mfree(xprime);
- return 0;
- }
- while (1) {
- if (esc_flag) {
- mfree(one);
- mfree(n);
- mfree(x);
- mfree(xprime);
- stop("esc");
- }
- // g = gcd(x' - x, n)
- t = msub(xprime, x);
- MSIGN(t) = 1;
- g = mgcd(t, n);
- mfree(t);
- if (MEQUAL(g, 1)) {
- mfree(g);
- if (--k == 0) {
- mfree(xprime);
- xprime = mcopy(x);
- l *= 2;
- k = l;
- }
- // x = (x ^ 2 + 1) mod n
- t = mmul(x, x);
- mfree(x);
- x = madd(t, one);
- mfree(t);
- t = mmod(x, n);
- mfree(x);
- x = t;
- continue;
- }
- push_factor(g, 1);
- if (mcmp(g, n) == 0) {
- mfree(one);
- mfree(n);
- mfree(x);
- mfree(xprime);
- return -1;
- }
- // n = n / g
- t = mdiv(n, g);
- mfree(n);
- n = t;
- // x = x mod n
- t = mmod(x, n);
- mfree(x);
- x = t;
- // xprime = xprime mod n
- t = mmod(xprime, n);
- mfree(xprime);
- xprime = t;
- break;
- }
- }
- }
- void
- push_factor(unsigned int *d, int count)
- {
- p1 = alloc();
- p1->k = NUM;
- p1->u.q.a = d;
- p1->u.q.b = mint(1);
- push(p1);
- if (count > 1) {
- push_symbol(POWER);
- swap();
- p1 = alloc();
- p1->k = NUM;
- p1->u.q.a = mint(count);
- p1->u.q.b = mint(1);
- push(p1);
- list(3);
- }
- }
|