#include #include #include #include /* This program has nothing to do with the solution to the actual problem. * Instead it tries to show how primality testing and factoring of numbers * up to 2^63 can be done in a reasonably efficient way. The program makes * use of a probabilistic algorithms, i.e. there is a very small possibility * that the result is in fact wrong. However, the probability can easily be * reduced to a value that is much smaller than the probability of a * hardware failure. */ /* Multiplication of two unsigned long long modulo mod. We do it this way to avoid overflows whenever possible. This should work reliably as long as x, y and mod are less than 2^63, i.e. suitable for a signed long long. */ unsigned long long mulmod (unsigned long long x, unsigned long long y, unsigned long long mod) { unsigned long long ret = 0ULL; x %= mod; y %= mod; while (y) { if (y & 1ULL) { ret += x; ret %= mod; } y >>= 1; x <<= 1; x %= mod; } return ret; } /* Use mulmod to calculate x*x modulo mod */ unsigned long long sq (unsigned long long x, unsigned long long mod) { return mulmod (x,x,mod); } /* Fast exponentiation modulo mod: Return value is equal to base to the * power of exp modulo mod. We use mulmod so this should work without * overflows for values of exp, mod and base up to 2^63-1 * References: * http://de.wikipedia.org/wiki/Bin%C3%A4re_modulo-Exponentiation */ unsigned long long expmod (unsigned long long base, unsigned long long exp, unsigned long long mod) { unsigned long long res = 1ULL; unsigned long long mp = base; while (exp) { if (exp % 2ULL) { res = mulmod (res, mp, mod); } exp >>= 1; mp = sq (mp, mod); } return res; } /* Greatest common divisor of a and b. */ unsigned long long ggt (unsigned long long a, unsigned long long b) { while (a && b) { if (a > b) { a %= b; } else { b %= a; } } return a+b; } /* Run a single instance of the miller/rabin primality test with base a * on the prime candidate p. If this function returns false the number * is definitely not a prime. If the function returns true p is prime * with a probability of about 75%. See references for miller_isprime below. */ int miller_isprime_1 (unsigned long long p, unsigned long long a) { unsigned long long ex = p-1ULL; unsigned long long r; unsigned long long t; unsigned long long oldt; if (p < 2) return 0; if (p < 4) return 1; r = 0; while (ex % 2ULL == 0) { ex /= 2; r++; } t = expmod (a, ex, p); if (t == 1ULL || t == p-1ULL) return 1; while (r > 1) { oldt = t; t = sq (t, p); if (t == p-1ULL) return 1; r--; } return 0; } /* Prime numbers up to 101. */ unsigned long long ps[] = { 2,3,5,7,11,13,17,19,23,29,31,37,41,43, 47,53,59,61,67,71,73,79,83,89,97,101,0 }; /* miller/rabin primality test on p. If the function returns false * p is not a prime. If the function returns true the value is very * probably prime. * References: * http://de.wikipedia.org/wiki/Miller-Rabin-Test * http://mathworld.wolfram.com/Rabin-MillerStrongPseudoprimeTest.html */ int miller_isprime (unsigned long long p) { int i; for (i=0; ps[i]; ++i) { if (!miller_isprime_1 (p, ps[i])) return 0; } return 1; } /* Pollard-Rho Method for factorization. The function an integral divisor * of p. However, there's a small probability that the function returns * p even if p is composite. This function shouldn't be called with a * prime number argument. It's worst case run time for a composite number * in the same order as the naive factoring algorithm. Note that the return * value is not guaranteed to be prime. * References: * http://de.wikipedia.org/wiki/Pollard-Rho-Methode * http://mathworld.wolfram.com/PollardRhoFactorizationMethod.html */ unsigned long long pollard_rho (unsigned long long p) { unsigned long long x = 2; unsigned long long y = 2; unsigned long long d = 1; unsigned long long c = 2; while (d == 1LL) { unsigned long long N = p; x = sq(x,p) + c; x %= p; y = sq(y,p) + c; y %= p; y = sq(y,p) + c; y %= p; if (x > y) { d = x-y; } else { d = y-x; } while (d && N) { if (d > N) { d %= N; } else { N %= d; } } d += N; if (d == p) break; } return d; } /* Return the smallest prime factor of p, returns p if p is itself prime. */ unsigned long long factor_naive (unsigned long long p) { unsigned long long d; if (p % 2ULL == 0) return 2ULL; for (d=3; d*d <= p; ++d) { if (p % d == 0) return d; } return p; } /* We do the following for each number in the input file: - Run the Miller/Rabin test and believe it if it says that the number is actually prime. This is the part where probability comes into play. If we really want to be sure, we could run the naive factoring at this point to verify that the number is really prime. - Try to factor the number using the Pollard-rho method. - If Pollard-rho fails to find a divisor use the naive method to find a factor. Should happen very rarely! */ int main (int argc, char ** argv) { char * fn = "INPUT"; if (argc >= 2) fn = argv[1]; int fd = open ("INPUT", O_RDONLY); if (fd < 0) { perror ("open"); exit (1); } long long Nr; while (read (fd, &Nr, sizeof (Nr))) { unsigned long long N; unsigned long long d; if (Nr < 0LL) Nr = -Nr; N = Nr; if (miller_isprime (N)) { printf ("%lld is prime\n", Nr); continue; } d = pollard_rho (N); /* This is almost never called. */ if (d == N) d = factor_naive (N); if (d == N) { printf ("%lld is prime\n", Nr); } else { long long di = d; printf ("%lld=%lld*%lld\n", Nr, di, Nr/di); } } return 0; }