/* from http://www.mindspring.com/~pate/course/chap08/squfof.c PLB 20AUG2010 */
/*
Author: Pate Williams (c) 1997
"Algorithm 8.7.2 (Shanks's SQUFOF). Given an odd
integer N, this algorithm tries to find a
non-trivial factor of N." -Henri Cohen-
See "A Course in Computational Algebraic Number
Theory" by Henri Cohen page 437.
*/
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
struct factor {long expon, prime;};
long exp_mod(long x, long b, long n)
/* returns x ^ b mod n */
{
long a = 1l, s = x;
while (b != 0) {
if (b & 1l) a = (a * s) % n;
b >>= 1;
if (b != 0) s = (s * s) % n;
}
return a;
}
long gcd(long a, long b)
{
long r;
a = labs(a);
b = labs(b);
while (b > 0)
r = a % b, a = b, b = r;
return a;
}
int RabinMiller(long N)
/* returns 0 if N is composite 1 otherwise */
{
long N1 = N - 1, a, b, c = 20, e, q = N1, t = 0;
if (N == 2) return 1;
while (!(q & 1)) q >>= 1, t++;
while (c > 0) {
do a = rand() % N; while (a == 0);
e = 0;
b = exp_mod(a, q, N);
if (b != 1) {
while (b != 1 && b != N1 && e <= t - 2) {
b = (b * b) % N;
e++;
}
if (b != N1) return 0;
}
c--;
}
return 1;
}
int square_test(long n, long *s)
{
int square = 1;
long q11[11], q63[63], q64[64], q65[65];
long k, r, t;
for (k = 0; k < 11; k++) q11[k] = 0;
for (k = 0; k < 6; k++) q11[(k * k) % 11] = 1;
for (k = 0; k < 63; k++) q63[k] = 0;
for (k = 0; k < 32; k++) q63[(k * k) % 63] = 1;
for (k = 0; k < 64; k++) q64[k] = 0;
for (k = 0; k < 32; k++) q64[(k * k) % 64] = 1;
for (k = 0; k < 65; k++) q65[k] = 0;
for (k = 0; k < 33; k++) q65[(k * k) % 65] = 1;
t = n % 64;
r = n % 45045l;
if (q64[t] == 0) square = 0;
if (square && q63[r % 63] == 0) square = 0;
if (square && q65[r % 65] == 0) square = 0;
if (square && q11[r % 11] == 0) square = 0;
if (square) *s = sqrt(n); else *s = 0;
return square;
}
void reduce(long A, long B, long C,
long *a, long *b, long *c)
{
long D, l, r, s;
*a = A;
*b = B;
*c = C;
D = B * B - 4 * A * C;
s = sqrt(D);
l = labs(*c);
if (l < s / 2)
r = - *b + 2 * l * ((*b + s) / (2 * l));
else
r = - *b + 2 * l;
*a = *c;
*b = r;
*c = (r * r - D) / (4 * *c);
}
void reduction(long A, long B, long C,
long *a, long *b, long *c)
{
long D, l, r, s;
*a = A;
*b = B;
*c = C;
D = B * B - 4 * A * C;
s = sqrt(D);
L1:
if (*b < 0) *b = - *b, *a = - *a, *c = - *c;
if (*b > labs(s - 2 * labs(*a)) && *b <= s) return;
l = labs(*c);
if (l < s / 2)
r = - *b + 2 * l * ((*b + s) / (2 * l));
else
r = - *b + 2 * l;
*a = *c;
*b = r;
*c = (r * r - D) / (4 * *c);
goto L1;
}
int SQUFOF(long *N, long *f, long *e)
{
int found;
long A, B, C, D, L, Q[1024], Q_size, a, b, b1, c;
long ap, bp, cp, d, i, j, s, z;
if (RabinMiller(*N)) return - 2;
if (square_test(*N, f)) {
*e = 2;
*N = 1;
return - 1;
}
if (*N % 4 == 1) {
D = *N;
d = sqrt(*N);
b = (d - 1) / 2;
b = (b << 1) + 1;
}
else {
D = 4 * *N;
d = sqrt(D);
b = 2 * (d / 2);
}
a = 1;
c = (b * b - D) / 4;
Q_size = 0;
i = 0;
L = sqrt(d);
ap = a;
bp = b;
cp = c;
L4:
reduce(ap, bp, cp, &A, &B, &C);
i++;
if (i & 1) goto L7;
if (square_test(labs(A), &z)) {
a = z;
found = 0;
for (j = 0; !found && j < Q_size; j++)
found = a == Q[j];
if (!found) goto L8;
}
if (A == 1) return 0;
L7:
if (labs(A) <= L) {
Q[Q_size] = A;
Q_size++;
}
ap = A;
bp = B;
cp = C;
goto L4;
L8:
s = gcd(a, gcd(B, D));
if (s > 1) {
*f = s;
*e = 2;
*N /= s * s;
return 1;
}
else
reduction(a, - B, a * C, &a, &b, &c);
L9:
b1 = b;
ap = a;
bp = b;
cp = c;
reduce(a, b, c, &a, &b, &c);
if (a == 1) {
a = ap;
b = bp;
c = cp;
goto L10;
}
if (b1 != b) goto L9;
L10:
if (a == 1) return 0;
if (a & 1) {
*e = 1;
*f = labs(a);
if (*f == 1 || *f == *N) return 0;
*N /= *f;
return 1;
}
*e = 1;
*f = labs(a / 2);
if (*f == 1 || *f == *N) return 0;
*N /= *f;
return 1;
}
int main(void)
{
int flag;
long N, count, i, j, n;
struct factor f[32], t;
for (;;) {
/* printf("N = "); PLB 20AUG2010 */
/* scanf("%ld", &N); PLB 20AUG2010 */
N = 633003781; /* PLB 20AUG2010 */
if (N <= 0) break;
if (N & 1) {
n = N;
flag = SQUFOF(&N, &f[0].prime, &f[0].expon);
if (flag == - 2) printf("%ld is prime\n", N);
else if (flag == - 1)
printf("%ld = %ld * %ld\n", f[0].prime * f[0].prime,
f[0].prime, f[0].prime);
else if (flag == 0) printf("%ld not factored\n", N);
else {
count = 1;
printf("%ld is composite\nfactors:\n", n);
for (;;) {
flag = SQUFOF(&N, &f[count].prime, &f[count].expon);
if (flag == - 1) break;
if (flag == 1) count++;
if (flag == - 2) {
f[count].expon = 1;
f[count].prime = N;
count++;
break;
}
if (flag == 0) break;
if (N == 1) break;
}
if (flag != 0) {
for (i = 0; i < count - 1; i++)
for (j = i + 1; j < count; j++)
if (f[i].prime > f[j].prime)
t = f[i], f[i] = f[j], f[j] = t;
for (i = 0; i < count; i++) {
printf("%ld", f[i].prime);
if (f[i].expon > 1) printf(" ^ %ld\n", f[i].expon);
else printf("\n");
}
}
else
printf("error - can't factor number\n");
}
}
else
printf("%ld must be odd!\n", N);
}
return 0;
}