{-
Copyright 2012 Ken Takusagawa

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 <http://www.gnu.org/licenses/>.

 -}
module Primality (isPrime
                 ,primes
                 ,miller_rabin_1
                 ,modular_exponentiation
                 , miller_rabin_isPrime_fast
                 , miller_rabin_isPrime
                 , miller_rabin_many
                 ) where
{
import Control.Exception;
-- import Debug.Trace(trace);
type Prime = Integer;
isPrime :: Prime -> Bool;
isPrime x = if (x<0) then (isPrime $ negate x) else
                if (x<2) then False else
                    testDivisibility (primes ()) x;
testDivisibility :: [Prime] -> Prime -> Bool;
testDivisibility primeList x = and (map (notDivisible x) (takeWhile (lessThanSquare x) primeList));
isPrime1 :: Prime -> Bool;
isPrime1 = testDivisibility smallPrimes;
lessThanSquare :: Prime -> Prime -> Bool;
lessThanSquare s x = x*x<=s;
notDivisible :: Prime -> Prime -> Bool;
notDivisible big small = (mod big small) /= 0;
smallPrimes :: [Prime];
smallPrimes = 3:(filter isPrime1 [5,7..]);
primes :: () -> [Prime];
primes _ =2:(filter isPrime1 [3,5..]);

modular_exponentiation :: Integer -> Integer -> Integer -> Integer;
modular_exponentiation a n modulo_p = 
    case (compare n 1) of 
      { (EQ) -> a;
        (GT) -> (case (mod n 2) of 
                 { (0) -> (let {
                           result :: Integer;
                           result = (modular_exponentiation a (div n 2) modulo_p)
                         }
                         in (modular_square result modulo_p));
                   (_) -> (mod ((*) a (modular_exponentiation a ((-) n 1) modulo_p)) modulo_p)
                 })
      };
modular_square :: Integer -> Integer -> Integer;
modular_square result modulo_p = (mod ((*) result result) modulo_p);

data Miller_rabin_result = Probable_prime | Composite | N(Integer) deriving (Show);

modular_square_miller_rabin :: Miller_rabin_result -> Integer -> Miller_rabin_result;
modular_square_miller_rabin x modulo_p = let {
  answer = modular_square_miller_rabin' x modulo_p;
  } in {- trace ("modular_square_miller_rabin " ++ (show x) ++ " = " ++ (show answer)) $ -} answer;
modular_square_miller_rabin' :: Miller_rabin_result -> Integer -> Miller_rabin_result;
modular_square_miller_rabin' x modulo_p = (case x of {
  (N(n)) -> (N (modular_square n modulo_p));
  (_) -> x
});

data Modular_classify = Negative_unity | Unity | Number;

pow_2_special::Integer->Integer->Integer->Miller_rabin_result;
pow_2_special a n modulo_p = 
case (mod n 2) of 
  {
    (0)->(let {
            result::Miller_rabin_result;
            result = (modular_square_miller_rabin (pow_2_special a (div n 2) modulo_p) modulo_p)
          }
          in (case result of 
                {
                  (N(b))->(case (classify b modulo_p) of 
                             {
                               (Negative_unity)->Probable_prime;
                               (Unity)-> Composite ; -- nontrivial sqrt of unity (div n 2)
                               (Number)->result
                             });
                  (_)->result
                }));
    (_)->(let {
            result::Integer;
            result = (modular_exponentiation a n modulo_p)
          }
          in (case (classify result modulo_p) of 
                {
                  (Negative_unity)->Probable_prime;
                  (Unity)->Probable_prime;
                  (Number)->(N result)
                }))
  };

classify::Integer->Integer->Modular_classify;
classify n modulo_p = 
    (case ((==) n ((-) modulo_p 1)) of 
       {
         (True)->Negative_unity;
         (_)->(case ((==) n 1) of 
                 {
                   (True)->Unity;
                   (_)->Number
                 })
       });

miller_rabin_1::Integer->Integer->Bool;
miller_rabin_1 n a = ((assert ((<) 2 n))((assert ((==) 1 (mod n 2)))((case (compare a ((-) n 1)) of {
  (LT)->(case (pow_2_special a ((-) n 1) n) of {
  (Probable_prime)->True;
  (Composite)->False;
  (N(1))->True;
  (_)->False
});
  (_)->True
}))));

miller_rabin_isPrime_fast :: Integer -> Bool;
miller_rabin_isPrime_fast n =
    if ( n < 118670087467 ) then
        if (n == 3215031751) then
            False
        else
            miller_rabin_many [2,3,5,7] n
    else error "too big miller_rabin_isPrime_fast";
        -- and $ map (miller_rabin_1 n) bases;
        
-- If n < 25,000,000,000 is a 2, 3, 5 and 7-SPRP, then either n  = 3,215,031,751 or n is prime [PSW80]. (This is actually true for n  < 118,670,087,467 [Jaeschke93].)

miller_rabin_isPrime :: Integer -> Bool;
miller_rabin_isPrime = miller_rabin_many miller_rabin_bases;
miller_rabin_many :: [Integer] -> Integer -> Bool;
miller_rabin_many bases n = and $ map (miller_rabin_1 n) bases;
miller_rabin_bases :: [Integer];
miller_rabin_bases = take 20 $ primes ();

}