/*==============================================================================

 File Name:

 ppFactor.c

 Description:

 Routines for integer factorization and primality prime testing.

 Functions:

 factor
 is_probably_prime
 is_almost_surely_prime

 LEGAL

 Primpoly Version 16.1  A Program for Computing Primitive Polynomials.
 Copyright (C) 19992021 by Sean Erik O'Connor. All Rights Reserved.

 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.

 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.

 You should have received a copy of the GNU General Public License
 along with this program. If not, see .

 The author's address is artificer!AT!seanerikoconnor!DOT!freeservers!DOT!com
 with !DOT! replaced by . and the !AT! replaced by @

==============================================================================*/
/*
 Include Files 
*/
#include
#include /* for rand() function */
#include "Primpoly.h"
/*==============================================================================
 factor 
================================================================================
DESCRIPTION
Factor n into primes. Record all the distinct prime factors and how
many times each occurs. Return the number of primes  1.
INPUT
n (int, n >= 1) The integer to be factored.
OUTPUT
primes (bigint *) List of distinct prime factors.
count (int *) List of how many times each factor occurred.
When n = 1, t = 0, and primes[ 0 ] = count[ 0 ] = 1.
RETURNS
t (int) Number of prime factors  1.
Prime factors and their multiplicites are in locations 0 to t.
EXAMPLE
2
For n = 156 = 2 * 3 * 13 we have
k primes[ k ] count[ k ]

0 2 2
1 3 1
2 13 1
METHOD
Method described by D.E. Knuth, ART OF COMPUTER PROGRAMMING, vol. 2, 2nd ed.,

in Algorithm A, pgs. 364365. The running time is O( max \/ p , p )
t1 t
where p is the largest prime divisor of n and p is the next largest.
t t1
(1) First divide out all multiples of 2 and 3 and count them.
(2) Next, divide n by all integers d >= 5 except multiples of 2 and 3.
(3) Halt either when all prime factors have been divided out (leaving n = 1)
or when the current value of n is prime. The stopping test
(d >  n/d  AND r != 0)
 
detects when n is prime.
In the example above, we divided out 2's and 3's first, leaving
n = 13. We continued with trial divisor d = 3. But now the stopping
test was activated because 5 >  13/5  = 2, and remainder = 3 (nonzero).
 
n = 13 itself is the last prime factor of 156. We return t = 2. There
are 2 + 1 = 3 distinct primes.
If we start with d = 5, and add 2 and 4 in succession, we will run through all
the integers except multiples of 2 and 3. To see this, observe the pattern:
Integers: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
Mult. of 2: x x x x x x x x
Mult. of 3: x x x x x
d: x x x x x
Separation: 2 4 2 4
The sequence of divisors d, includes all the primes, so is safe to use for
factorization testing.
Theorem. Suppose no divisor d in the sequence divides n. Suppose at some point,
q < d with r != 0. Then n must be prime.
Proof. We begin with an integer n which has all powers of 2 and 3 removed.
Assume, to the contrary, that n is composite: n = p ... p
t 1 t
n >= p since p is the smallest prime factor.
1 1
2
n >= p since n is composite, so has at least 2 prime factors.
1
2
>= (d + 1) since p > d implies p >= (d + 1). p > d because we assumed
1 1 1
that no d in the sequence divided n, so d couldn't be any of the prime divisors
p
i
2 2 2
= d + 2d + 1 = d + d + d + 1 > d + d
2
> d + n mod d since n mod d < d
2
>  n / d  d + n mod d because our test said q < d, so  n / d  d < d
   
= n. So we get the contradiction n > n. Thus n must be prime.

Note that this is sharper than the bound d > \/ n
Theorem. Conversely, if n is a prime, no d will divide it, so at some point,
d will be large enough that q =  n / d  < d. r != 0 since n is prime.
 
Theorem. The factoring algorithm works.
Proof. If n == 1 we exit immediately. If n is composite, we divide out all powers
of 2 and 3 (if any). Since d runs through all possible prime divisors, we divide
these out also. Composite trial d causes no problem; prime factors of d are divided
out of n before we try to divide by d, so such a d cannot divide n.
We end in one of two ways (1) All divisors of n have been divided out in which case n = 1.
(2) n is a prime, so the stopping test is activiated and n != 1 is recorded as a prime
divisor.
BUGS
Can be slow when n is a prime. We could do a probabilistic test for
the primality of n at the statement, "n_is_prime = (r != 0 && q < d) ;"
which would speed up the test.
Or we could implement the Pollard rho algorithm instead.

 Function Call 
*/
int
factor( bigint n, bigint * primes, int * count )
{
/*
 Local Variables 
*/
int
t = 0, /* Array index for primes and count. */
cnt = 0, /* Counter for factors of 2 and 3. */
new_d = YES, /* YES if current divisor is different from
the previous one. */
n_is_prime = NO, /* YES if current value of n is prime. */
flag = YES ; /* Alternately YES and NO. */
bigint
q, /* The quotient, floor( n / d ) */
r, /* The remainder, n mod d */
d = 5 ; /* Trial divisor and its first value. */
/*
 Function Body 
*/
/* Handle special cases for speed.
These are cases where r = (p^n  1)/(p1) is a large prime or has large prime
factors and p^n is less than the maximum size allowed by this program.
(I don't have all these cases).
*/
if (n == 1)
{
primes[ 0 ] = count[ 0 ] = 1 ;
return( 0 ) ;
}
else if (n == (bigint)4611686018427387903ULL) /* (2^62  1)/(2  1) */
{
primes[ 0 ] = 3 ;
count [ 0 ] = 1 ;
primes[ 1 ] = 715827883ULL ;
count [ 1 ] = 1 ;
primes[ 2 ] = 2147483647ULL ;
count [ 2 ] = 1 ;
return( 2 ) ;
}
else if (n == (bigint)2305843009213693951ULL) /* (2^61  1)/(2  1) */
{
primes[ 0 ] = (bigint)2305843009213693951ULL ;
count [ 0 ] = 1 ;
return( 0 ) ;
}
while( n % 2 == 0 ) /* Remove factors of 2. */
{
n /= 2 ;
++cnt ;
}
if (cnt != 0)
{
primes[ t ] = 2 ;
count[ t++ ] = cnt ;
}
cnt = 0 ;
while( n % 3 == 0 ) /* Remove factors of 3. */
{
n /= 3 ;
++cnt ;
}
if (cnt != 0)
{
primes[ t ] = 3 ;
count[ t++ ] = cnt ;
}
/*
Factor the rest of n. Continue until n = 1 (all factors divided out)
or until n is prime (so n itself is the last prime factor).
*/
do {
/* Integer divide to get quotient and remainder: n = qd + r */
q = n / d ;
r = n % d ;
n_is_prime = (r != 0 && q < d) ; /* n is prime. */
if (r == 0)
{
n = q ; /* Divide d out of n. */
if (new_d) /* New prime factor. */
{
primes[ t ] = d ;
count[ t++ ] = 1 ;
new_d = NO ;
}
else
++count[ t1 ] ; /* Same divisor again. Increment its count. */
}
else {
d += (flag ? 2 : 4) ; /* d did not divide n. Try a new divisor. */
flag = !flag ;
new_d = YES ;
}
} while ( ! n_is_prime && (n != 1)) ;
if (n == 1) /* All factors were divided out. */
return( t  1 ) ;
else { /* Current value of n is prime. It is the last prime factor. */
primes[ t ] = n ;
count[ t ] = 1 ;
return( t ) ;
}
} /* ========================= end of function factor ======================= */
/*
4
Phi[ 5  1 ] / 4
4 4
5  1 = 624 = 2 3 13
Phi = 624 (1  1/2) (1  1/3)(1  1/13) = 192
192 / 4 = 48
Handle special cases.
Warning: can take a long time due to factorization.
*/
bigint
EulerPhi( bigint n )
{
int prime_count = 0 ;
int i = 0 ;
bigint phi ;
/* Factor n into primes. */
bigint primes[ MAXNUMPRIMEFACTORS ] ;
int count [ MAXNUMPRIMEFACTORS ] ;
/* Error return and special cases. */
if (n <= 0)
return 0 ;
else if (n == 1)
return 1 ;
/* Factor n >= 2 into distinct primes. */
prime_count = factor( n, primes, count ) ;
/* Compute Euler phi[ n ] = */
/* */
/* num prime factors */
/*  */
/* n   (1  1/p ) */
/* i = 1 i */
phi = n ;
for (i = 0 ; i <= prime_count ; ++i)
{
phi /= primes[ i ] ;
phi *= (primes[ i ]  1) ;
}
return phi ;
}
/*==============================================================================
 is_probably_prime 
================================================================================
DESCRIPTION
Test whether the number n is probably prime. If n is composite
we are correct always. If n is prime, we will be wrong no more
than about 1/4 of the time on average, probably less for
any x and n.
INPUT
n (int, n >= 0) Number to test for primality.
x (int, 1 < x < n for n > 6) A random integer.
RETURNS
YES If n is probably prime.
NO If n is definitely not prime or n < 0,
or x is out of range.
EXAMPLE
METHOD
MillerRabin probabilistic primality test, described by Knuth.
If n = 1 + 2^k q is prime, and x^q mod n != 1 the sequence,
{ x^q mod n, x^(2q) mod n, ..., x^(2^k q) mod n}
will end with 1 and the element in the sequence just before
1 first appears = n1, since y^2 = 1 (mod p) satisfies
y = +1 or y = 1 only. The probability the algorithm fails is
bounded above by about 1/4 for all n.
BUGS
None.

 Function Call 
*/
int is_probably_prime( int n, int x )
{
/*
 Local Variables 
*/
int
j,
k = 0,
nm1,
q = 0,
y = 0 ;
/*
 Function Body 
*/
/* Not prime if input is out of range. */
if (n < 0)
return NO ;
/* Handle special cases. */
if (n == 0  n == 1  n == 4)
return NO ;
if (n == 2  n == 3  n == 5)
return YES ;
/* Even numbers aren't prime. */
if (n % 2 == 0)
return NO ;
/* Return not prime if x is out of range. */
if (x <= 1  x >= n)
return NO ;
/* Factor out powers of 2 to get n = 1 + 2^k q, q odd. */
nm1 = n  1 ;
k = 0 ;
while( nm1 % 2 == 0)
{
nm1 /= 2 ;
++k ;
}
q = nm1 ;
/* y = x^q (mod n) */
y = power_mod( x, q, n ) ;
for (j = 0 ; j < k ; ++j)
{
if ( (j == 0 && y == 1)  (y == n1))
return YES ;
else if (j > 0 && y == 1)
return NO ;
/* Compute y^2 (mod n) */
y = power_mod( y, 2, n ) ;
}
/* Shouldn't get here, but do it anyway for safety. */
return NO ;
} /* =============== end of function is_probably_prime ===================== */
/*==============================================================================
 is_almost_surely_prime 
================================================================================
DESCRIPTION
Test whether the number n is almost surely prime. If n is
composite, we always return NO. If n is prime, the
probability of returning YES successfully is
NUM_PRIME_TEST_TRIALS
1  (1/4)
INPUT
n (int, n >= 0)
RETURNS
YES If n is probably prime with a very high degree of
probability.
NO If n is definitely not prime.
1 If inputs are out of range.
EXAMPLE
METHOD
Use the MillerRabin probabilitic primality test for lots of
random integers x. If the test shows n is composite for any
given x, is is nonprime for sure.
If it passes the the primality test for several random values
of x, the probability is it prime is approximately
NUM_PRIME_TEST_TRIALS
1  2
assuming the random number generator is really random.
BUGS
The system pseudorandom number generator needs to be a good one
(i.e. passes the spectral test and other statistical tests).

 Function Call 
*/
int is_almost_surely_prime( int n )
{
/*
 Local Variables 
*/
int
trial = 0,
x = 3 ;
/*
 Function Body 
*/
/* Seed the randomnumber generator. */
srand( 314159 );
for (trial = 1 ; trial <= NUM_PRIME_TEST_TRIALS ; ++trial)
{
/* Generate a new random integer such that 1 < x < n. */
x = rand() % n ;
if (x <= 1) x = 3 ;
/* Definitely not prime. */
if (is_probably_prime( n, x ) == NO)
return NO ;
}
/* Almost surely prime because it passed the probably prime test
* above many times.
*/
return YES ;
} /* ============== end of function is_almost_surely_prime ================ */