/* Author: Pate Williams (c) 1997 Special number field sieve factorization of r ^ 3 + 2. See "Prime Numbers and Computer Methods for Factorization" by Hans Riesel second edition pages 214-216. */ #include <math.h> #include <stdio.h> #include "lip.h" /* default: -200, +200, +100, 10, 10, 20, 21 */ /* crashes: -500, +500, +250, 15, 15, 30, 31 */ #define A_MIN -200 #define A_MAX +200 #define B_MAX +100 #define FB1_SIZE 10 /* number of rational primes */ #define FB2_SIZE 10 /* number of algebraic primes */ #define FB_SIZE 20 /* FB1_SIZE + FB2_SIZE */ #define FB_SIZE1 21 /* FB_SIZE + 1 */ struct element {long a, b, c;}; struct zelement {verylong a, b, c;}; struct factor {long expon; verylong prime;}; void norm(verylong zN, struct zelement e, verylong *nm) { verylong a = 0, b = 0, c = 0; verylong za = 0, zb = 0, zc = 0; zcopy(e.a, &a); zcopy(e.b, &b); zcopy(e.c, &c); zsexpmod(a, 3, zN, &za); zsexpmod(b, 3, zN, &zb); zsmulmod(zb, 2, zN, &zc); zsubmod(za, zc, zN, &zb); zsexpmod(c, 3, zN, &za); zsmulmod(za, 4, zN, &zc); zaddmod(zb, zc, zN, &za); zmulmod(a, b, zN, &zb); zmulmod(zb, c, zN, &zc); zsmulmod(zc, 6, zN, &zb); zaddmod(za, zb, zN, nm); zfree(&a); zfree(&b); zfree(&c); zfree(&za); zfree(&zb); zfree(&zc); } void multiply(verylong zN, struct zelement e1, struct zelement e2, struct zelement *product) { verylong za = 0, zb = 0, zc = 0; zmulmod(e1.a, e2.a, zN, &za); zmulmod(e1.b, e2.c, zN, &zb); zsmulmod(zb, 2, zN, &zc); zsubmod(za, zc, zN, &zb); zmulmod(e1.c, e2.b, zN, &za); zsmulmod(za, 2, zN, &zc); zsubmod(zb, zc, zN, &product->a); zmulmod(e1.a, e2.b, zN, &za); zmulmod(e1.b, e2.a, zN, &zb); zaddmod(za, zb, zN, &zc); zmulmod(e1.c, e2.c, zN, &za); zsmulmod(za, 2, zN, &zb); zsubmod(zc, zb, zN, &product->b); zmulmod(e1.a, e2.c, zN, &za); zmulmod(e1.b, e2.b, zN, &zb); zaddmod(za, zb, zN, &zc); zmulmod(e1.c, e2.a, zN, &za); zaddmod(zc, za, zN, &product->c); zfree(&za); zfree(&zb); zfree(&zc); } void power(verylong zN, struct zelement x, long b, struct zelement *a) /* returns a = x ^ b */ { struct zelement s = {0, 0, 0}, t = {0, 0, 0}; zone(&a->a); zzero(&a->b); zzero(&a->c); zcopy(x.a, &s.a); zcopy(x.b, &s.b); zcopy(x.c, &s.c); while (b != 0) { if (b & 1) { multiply(zN, *a, s, &t); zcopy(t.a, &a->a); zcopy(t.b, &a->b); zcopy(t.c, &a->c); } b >>= 1; if (b != 0) { multiply(zN, s, s, &t); zcopy(t.a, &s.a); zcopy(t.b, &s.b); zcopy(t.c, &s.c); } } zfree(&s.a); zfree(&s.b); zfree(&s.c); zfree(&t.a); zfree(&t.b); zfree(&t.c); } int rational_factor(long a, long b, verylong zr, long *FB1, long *expon) { long c, i, p; verylong za = 0, zn = 0; zsmul(zr, b, &za); zsadd(za, a, &zn); for (i = 0; i < FB1_SIZE; i++) expon[i] = 0; if (zscompare(zn, 0) < 0) { expon[0] = 1; znegate(&zn); } if (zscompare(zn, 1) <= 0) return 0; for (i = 1; i < FB1_SIZE; i++) { p = FB1[i]; if (zsdiv(zn, p, &za) == 0) { c = 1; zcopy(za, &zn); while (zsdiv(zn, p, &za) == 0) { c++; zcopy(za, &zn); } expon[i] = c; if (zscompare(zn, 1) == 0) { zfree(&za); zfree(&zn); return 1; } } } zfree(&za); zfree(&zn); return 0; } int divide(verylong zN, struct zelement alpha, struct zelement beta, struct zelement *gamma) { int flag = 0; verylong a = 0, b = 0, c = 0, d = 0, e = 0; verylong f = 0, g = 0, h = 0, i = 0, x = 0; verylong aa = 0, ab = 0, ac = 0; verylong ba = 0, bb = 0, bc = 0; verylong dd = 0, ee = 0, ff = 0; verylong za = 0, zb = 0, zc = 0, zd = 0, ze = 0; verylong zf = 0; verylong zN1 = 0; zcopy(alpha.a, &aa); zcopy(alpha.b, &ab); zcopy(alpha.c, &ac); zcopy(beta.a, &ba); zcopy(beta.b, &bb); zcopy(beta.c, &bc); zcopy(aa, &g); zcopy(ab, &h); zcopy(ac, &i); zcopy(ba, &d); zcopy(bb, &e); zcopy(bc, &f); if (zcompare(aa, ba) == 0 && zcompare(ab, bb) == 0 && zcompare(ac, bc) == 0) { zfree(&d); zfree(&e); zfree(&f); zfree(&g); zfree(&h); zfree(&i); zfree(&aa); zfree(&ab); zfree(&ac); zfree(&ba); zfree(&bb); zfree(&bc); zone(&gamma->a); zzero(&gamma->b); zzero(&gamma->c); return 1; } znegate(&aa); znegate(&ab); znegate(&ac); if (zcompare(aa, ba) == 0 && zcompare(ab, bb) == 0 && zcompare(ac, bc) == 0) { zfree(&d); zfree(&e); zfree(&f); zfree(&g); zfree(&h); zfree(&i); zfree(&aa); zfree(&ab); zfree(&ac); zfree(&ba); zfree(&bb); zfree(&bc); zone(&gamma->a); znegate(&gamma->a); zzero(&gamma->b); zzero(&gamma->c); return 1; } zsexpmod(d, 3, zN, &za); zsexpmod(e, 3, zN, &zb); zsexpmod(f, 3, zN, &zc); zsmulmod(zb, 2, zN, &zd); zsmulmod(zc, 4, zN, &ze); zsubmod(za, zd, zN, &zf); zaddmod(zf, ze, zN, &za); zmulmod(d, e, zN, &zb); zmulmod(zb, f, zN, &zc); zsmulmod(zc, 8, zN, &zd); zaddmod(za, zd, zN, &x); zmulmod(d, d, zN, &dd); zmulmod(e, e, zN, &ee); zmulmod(f, f, zN, &ff); zmulmod(dd, g, zN, &za); zmulmod(ee, h, zN, &zb); zsmulmod(zb, 2, zN, &zc); zsubmod(za, zc, zN, &zb); zmulmod(ff, i, zN, &za); zsmulmod(za, 4, zN, &zc); zaddmod(zb, zc, zN, &za); zmulmod(d, f, zN, &zb); zmulmod(zb, h, zN, &zc); zsmulmod(zc, 2, zN, &zd); zaddmod(za, zd, zN, &zb); zmulmod(e, f, zN, &za); zmulmod(za, g, zN, &zc); zsmulmod(zc, 2, zN, &zd); zaddmod(zb, zd, zN, &za); zmulmod(d, e, zN, &zb); zmulmod(zb, i, zN, &zc); zsmulmod(zc, 2, zN, &zd); zaddmod(za, zd, zN, &a); zmulmod(dd, h, zN, &za); zmulmod(ee, i, zN, &zb); zsmulmod(zb, 2, zN, &zc); zsubmod(za, zc, zN, &zb); zmulmod(ff, g, zN, &za); zsmulmod(za, 2, zN, &zc); zsubmod(zb, zc, zN, &za); zmulmod(d, e, zN, &zb); zmulmod(zb, g, zN, &zc); zsubmod(za, zc, zN, &zb); zmulmod(d, f, zN, &za); zmulmod(za, i, zN, &zc); zsmulmod(zc, 2, zN, &zd); zaddmod(zb, zd, zN, &za); zmulmod(e, f, zN, &zb); zmulmod(zb, h, zN, &zc); zsmulmod(zc, 2, zN, &zd); zaddmod(za, zd, zN, &b); zmulmod(dd, i, zN, &za); zmulmod(ee, g, zN, &zb); zaddmod(za, zb, zN, &zc); zmulmod(ff, h, zN, &za); zsmulmod(za, 2, zN, &zb); zsubmod(zc, zb, zN, &za); zmulmod(e, f, zN, &zb); zmulmod(zb, i, zN, &zc); zsmulmod(zc, 2, zN, &zb); zaddmod(za, zb, zN, &zc); zmulmod(d, e, zN, &za); zmulmod(za, h, zN, &zb); zsubmod(zc, zb, zN, &za); zmulmod(d, f, zN, &zb); zmulmod(zb, g, zN, &zc); zsubmod(za, zc, zN, &c); zsadd(zN, - 1, &zN1); if (zscompare(x, 0) != 0 && zscompare(x, 1) != 0 && zcompare(x, zN1) != 0) { zrshift(zN, 1, &zf); if (zcompare(x, zf) > 0) { zsub(x, zN, &za); zcopy(za, &x); } if (zcompare(a, zf) > 0) { zsub(a, zN, &za); zcopy(za, &zb); zabs(&zb); if (zcompare(zb, zN) == 0) zzero(&za); zcopy(za, &a); } if (zcompare(b, zf) > 0) { zsub(b, zN, &za); zcopy(za, &zb); zabs(&zb); if (zcompare(zb, zN) == 0) zzero(&za); zcopy(za, &b); } if (zcompare(c, zf) > 0) { zsub(c, zN, &za); zcopy(za, &zb); zabs(&zb); if (zcompare(zb, zN) == 0) zzero(&za); zcopy(za, &c); } zdiv(a, x, &za, &zb); if (zscompare(zb, 0) == 0) { zcopy(za, &a); zdiv(b, x, &za, &zb); if (zscompare(zb, 0) == 0) { zcopy(za, &b); zdiv(c, x, &za, &zb); if (zscompare(zb, 0) == 0) { zcopy(za, &c); if (zscompare(a, 0) < 0) { zadd(a, zN, &za); zmod(za, zN, &zb); zcopy(zb, &a); } if (zscompare(b, 0) < 0) { zadd(b, zN, &za); zmod(za, zN, &zb); zcopy(zb, &b); } if (zscompare(c, 0) < 0) { zadd(c, zN, &za); zmod(za, zN, &zb); zcopy(zb, &c); } flag = 1; } } } } else { zzero(&a); zzero(&b); zzero(&c); } zcopy(a, &gamma->a); zcopy(b, &gamma->b); zcopy(c, &gamma->c); zfree(&a); zfree(&b); zfree(&c); zfree(&d); zfree(&e); zfree(&f); zfree(&g); zfree(&h); zfree(&i); zfree(&x); zfree(&aa); zfree(&ab); zfree(&ac); zfree(&ba); zfree(&bb); zfree(&bc); zfree(&dd); zfree(&ee); zfree(&ff); zfree(&za); zfree(&zb); zfree(&zc); zfree(&zd); zfree(&ze); zfree(&zf); zfree(&zN1); return flag; } int algebraic_factor(verylong zN, struct zelement *FB2, long *expon, struct zelement e) { int flag = 0; long count, i; static struct element u[8] = {{- 1, 1, - 1}, {5, - 4, 3}, {- 19, 15, - 12}, {73, - 58, 46}, {- 281, 223, - 177}, {1081, - 858, 681}, {- 4159, 3301, - 2620}, {16001, - 12700, 10080}}; struct zelement n = {0, 0, 0}, p = {0, 0, 0}; struct zelement q = {0, 0, 0}, unit = {0, 0, 0}; verylong nm = 0, za = 0, zb = 0, zc = 0; verylong zN1 = 0; for (i = 0; i < FB2_SIZE; i++) expon[i] = 0; zsadd(zN, - 1, &zN1); norm(zN, e, &nm); if (zscompare(nm, 0) == 0) { zfree(&nm); return 0; } zcopy(e.a, &n.a); zcopy(e.b, &n.b); zcopy(e.c, &n.c); for (i = 2; i < FB2_SIZE; i++) { zcopy(FB2[i].a, &p.a); zcopy(FB2[i].b, &p.b); zcopy(FB2[i].c, &p.c); if (divide(zN, n, p, &q)) { count = 1; zcopy(q.a, &n.a); zcopy(q.b, &n.b); zcopy(q.c, &n.c); while (divide(zN, n, p, &q)) { count++; zcopy(q.a, &n.a); zcopy(q.b, &n.b); zcopy(q.c, &n.c); } expon[i] = count; if (zscompare(n.a, 1) == 0 && zscompare(n.b, 0) == 0 && zscompare(n.c, 0) == 0) { expon[0] = 0; zfree(&nm); zfree(&n.a); zfree(&n.b); zfree(&n.c); zfree(&p.a); zfree(&p.b); zfree(&p.c); zfree(&q.a); zfree(&q.c); return 1; } if (zcompare(n.a, zN1) == 0 && zscompare(n.b, 0) == 0 && zscompare(n.c, 0) == 0) { expon[0] = 1; zfree(&nm); zfree(&n.a); zfree(&n.b); zfree(&n.c); zfree(&p.a); zfree(&p.b); zfree(&p.c); zfree(&q.a); zfree(&q.c); return 1; } } } norm(zN, n, &nm); if (zscompare(nm, 1) == 0 || zcompare(nm, zN1) == 0) { for (i = 0; i < 8; i++) { zintoz(u[i].a, &za); zintoz(u[i].b, &zb); zintoz(u[i].c, &zc); if (zcompare(n.a, unit.a) == 0 && zcompare(n.b, unit.b) == 0 && zcompare(n.c, unit.c) == 0) { zfree(&nm); zfree(&za); zfree(&zb); zfree(&zc); zfree(&n.a); zfree(&n.b); zfree(&n.c); zfree(&p.a); zfree(&p.b); zfree(&p.c); zfree(&q.a); zfree(&q.c); zfree(&unit.a); zfree(&unit.b); zfree(&unit.c); return 0; } znegate(&n.a); znegate(&n.b); znegate(&n.c); if (zcompare(n.a, unit.a) == 0 && zcompare(n.b, unit.b) == 0 && zcompare(n.c, unit.c) == 0) { zfree(&nm); zfree(&za); zfree(&zb); zfree(&zc); zfree(&n.a); zfree(&n.b); zfree(&n.c); zfree(&p.a); zfree(&p.b); zfree(&p.c); zfree(&q.a); zfree(&q.c); zfree(&unit.a); zfree(&unit.b); zfree(&unit.c); return 0; } znegate(&n.a); znegate(&n.b); znegate(&n.c); } zcopy(FB2[1].a, &p.a); zcopy(FB2[1].b, &p.b); zcopy(FB2[1].c, &p.c); if (divide(zN, n, p, &q)) { count = 1; zcopy(q.a, &n.a); zcopy(q.b, &n.b); zcopy(q.c, &n.c); if ((zscompare(n.a, 1) != 0 || zcompare(n.a, zN1) != 0) || zscompare(n.b, 0) != 0 || zscompare(n.c, 0) != 0) { while (divide(zN, n, p, &q)) { count++; zcopy(q.a, &n.a); zcopy(q.b, &n.b); zcopy(q.c, &n.c); if ((zscompare(n.a, 1) == 0 || zcompare(n.a, zN1) == 0) && zscompare(n.b, 0) == 0 && zscompare(n.c, 0) == 0) break; } } expon[1] = count; expon[0] = zscompare(n.a, 1) == 0 ? 0 : 1; flag = 1; } } zfree(&nm); zfree(&za); zfree(&zb); zfree(&zc); zfree(&n.a); zfree(&n.b); zfree(&n.c); zfree(&p.a); zfree(&p.b); zfree(&p.c); zfree(&q.a); zfree(&q.c); zfree(&unit.a); zfree(&unit.b); return flag; } void special_number_field_sieve(verylong zN, verylong zr, long *FB1) { int done, factored = 0, found, prime = 0; long count = 0, f_count = 0; long a, b, exp; long D, c[FB_SIZE1], i, j, k, s, X[FB_SIZE1]; long a_expon[FB2_SIZE], r_expon[FB1_SIZE], expon[FB_SIZE]; long e_matrix[FB_SIZE1][FB_SIZE], matrix[FB_SIZE1][FB_SIZE]; struct element fb2[28] = {{0, 0, 0}, {1, 1, 0}, {0, 1, 0}, {- 1, 1, 0}, {1, 0, 1}, {1, 1, - 1}, {1, - 2, 0}, {3, 0, - 1}, {3, - 1, 0}, {1, - 3, 1}, {3, 1, 1}, {1, 3, 0}, {3, 0, 2}, {1, 2, - 2}, {3, - 1, 3}, {1, - 2, 1}, {3, 4, 0}, {1, 0, - 3}, {- 1, 4, - 2}, {3, - 3, - 1}, {- 1, 0, 2}, {1, 1, 2}, {1, 0, 3}, {1, 4, 0}, {5, 0, 2}, {1, - 5, 2}, {- 1, 1, 4}, {3, 4, - 2}}; struct factor f[32] = {{0, 0}}, t = {0, 0}; struct zelement e = {0, 0, 0}, z = {0, 0, 0}; struct zelement product = {0, 0, 0}; struct zelement FB2[28] = {{0, 0, 0}}; verylong za = 0, zb = 0, zc = 0, zd = 0, zp = 0; verylong zx = 0, zy = 0; zone(&z.a); zzero(&z.b); zzero(&z.c); for (i = 0; i < 28; i++) { zintoz(fb2[i].a, &FB2[i].a); zintoz(fb2[i].b, &FB2[i].b); zintoz(fb2[i].c, &FB2[i].c); } for (a = A_MIN; !factored && !prime && a <= A_MAX; a++) { for (b = 0; !factored && !prime && b <= B_MAX; b++) { if (a == 0 && b == 0) continue; if (rational_factor(a, b, zr, FB1, r_expon)) { if (a >= 0) zintoz(a, &e.a); else { zintoz(a, &za); zadd(za, zN, &e.a); } zintoz(b, &e.b); zzero(&e.c); if (algebraic_factor(zN, FB2, a_expon, e)) { for (i = 0; i < FB1_SIZE; i++) e_matrix[count][i] = r_expon[i]; for (i = 0; i < FB2_SIZE; i++) e_matrix[count][i + FB1_SIZE] = a_expon[i]; for (i = 0; i < FB_SIZE; i++) matrix[count][i] = e_matrix[count][i] & 1; count++; printf("\b\b%2ld", count); if (count == FB_SIZE1) { count = 0; for (i = 0; i < FB_SIZE1; i++) c[i] = - 1; for (k = 0; k < FB_SIZE1; k++) { found = 0; j = 0; while (!found && j < FB_SIZE) { found = matrix[k][j] != 0 && c[j] < 0; if (!found) j++; } if (found) { matrix[k][j] = 1; for (i = 0; i < FB_SIZE; i++) { if (i != j ) { D = matrix[k][i]; matrix[k][i] = 0; for (s = k + 1; s < FB_SIZE1; s++) matrix[s][i] = (matrix[s][i] + D * matrix[s][j]) & 1; } } c[j] = k; } else { for (j = 0; j < FB_SIZE1; j++) X[j] = 0; for (j = 0; j < FB_SIZE1; j++) { if (j == k) X[j] = 1; else for (s = 0; s < FB_SIZE; s++) if (c[s] == j) X[j] = matrix[k][s]; } for (i = 0; i < FB_SIZE; i++) expon[i] = 0; for (i = 0; i < FB_SIZE1; i++) { if (X[i] != 0) { for (j = 0; j < FB_SIZE; j++) expon[j] += e_matrix[i][j]; } } zone(&zx); for (i = 1; i < FB1_SIZE; i++) { zintoz(FB1[i], &zp); zsexp(zp, expon[i] / 2, &za); zmulmod(za, zx, zN, &zb); zcopy(zb, &zx); } if (expon[0] / 2 & 1) znegate(&zx); if (expon[FB1_SIZE] / 2 & 1) { znegate(&z.a); znegate(&z.b); znegate(&z.c); } for (i = 1; i < FB2_SIZE; i++) { power(zN, FB2[i], expon[FB1_SIZE + i] / 2, &e); multiply(zN, z, e, &product); zcopy(product.a, &z.a); zcopy(product.b, &z.b); zcopy(product.c, &z.c); } zmulmod(zr, z.b, zN, &za); zmulmod(zr, zr, zN, &zb); zmulmod(zb, z.c, zN, &zc); zaddmod(za, zc, zN, &zb); zaddmod(zb, z.a, zN, &zy); zsub(zx, zy, &za); zgcd(za, zN, &zd); if (zscompare(zd, 1) != 0 && zcompare(zd, zN) != 0) { exp = 0; zmod(zN, zd, &zb); done = zscompare(zb, 0) != 0; while (!done) { zdiv(zN, zd, &za, &zb); done = zscompare(zb, 0) != 0; if (!done ) { exp++; zcopy(za, &zN); } } f[f_count].expon = exp; zcopy(zd, &f[f_count++].prime); factored = zscompare(zN, 1) == 0; if (!factored) prime = zprobprime(zN, 5); } } } } } } } } if (!factored) { f[f_count].expon = 1; zcopy(zN, &f[f_count++].prime); } for (i = 0; i < f_count - 1; i++) { for (j = i + 1; j < f_count; j++) { if (zcompare(f[i].prime, f[j].prime) > 0) { t.expon = f[i].expon; zcopy(f[i].prime, &t.prime); f[i].expon = f[j].expon; zcopy(f[j].prime, &f[i].prime); f[j].expon = t.expon; zcopy(t.prime, &f[j].prime); } } } printf("\nfactors:\n"); for (i = 0; i < f_count; i++) { printf("\t"); zwrite(f[i].prime); if (f[i].expon != 1) printf(" ^ %ld\n", f[i].expon); else printf("\n"); } if (!prime) printf("\tlast factor is composite\n"); zfree(&za); zfree(&zb); zfree(&zc); zfree(&zd); zfree(&zp); zfree(&zx); zfree(&zy); zfree(&e.a); zfree(&e.b); zfree(&e.c); zfree(&z.a); zfree(&z.b); zfree(&z.c); zfree(&product.a); zfree(&product.b); zfree(&product.c); for (i = 0; i < 28; i++) { zfree(&FB2[i].a); zfree(&FB2[i].b); zfree(&FB2[i].c); } for (i = 0; i < f_count; i++) zfree(&f[i].prime); zfree(&t.prime); } int main(void) { long FB1[FB1_SIZE], i; verylong zN = 0, za = 0, zr = 0; zpstart2(); FB1[0] = - 1; for (i = 1; i < FB1_SIZE; i++) FB1[i] = zpnext(); printf("factors the number r ^ 3 + 2\n"); for (;;) { printf("enter r or 0 to quit:\n"); zread(&zr); if (zscompare(zr, 0) == 0) break; zsexp(zr, 3, &za); zsadd(za, 2, &zN); zwrite(zN); if (!zprobprime(zN, 5)) { printf(" is composite\n"); special_number_field_sieve(zN, zr, FB1); } else printf(" is prime\n"); } zfree(&zN); zfree(&za); zfree(&zr); return 0; }