-- |
-- Module:      Math.NumberTheory.Primes.Testing.Certified
-- Copyright:   (c) 2011 Daniel Fischer
-- Licence:     MIT
-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
--
-- Deterministic primality testing.

{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE ViewPatterns        #-}

module Math.NumberTheory.Primes.Testing.Certified
  ( isCertifiedPrime
  ) where

import Data.List (foldl')
import Data.Bits ((.&.))
import Data.Mod
import Data.Proxy
import GHC.Integer.GMP.Internals (powModInteger)
import GHC.TypeNats (SomeNat(..), someNatVal)

import Math.NumberTheory.Roots (integerSquareRoot)
import Math.NumberTheory.Primes (unPrime)
import Math.NumberTheory.Primes.Factorisation.TrialDivision (trialDivisionPrimeTo, trialDivisionTo, trialDivisionWith)
import Math.NumberTheory.Primes.Factorisation.Montgomery (montgomeryFactorisation, smallFactors, findParms)
import Math.NumberTheory.Primes.Testing.Probabilistic (bailliePSW, isPrime, isStrongFermatPP, lucasTest)
import Math.NumberTheory.Primes.Sieve.Eratosthenes (primeList, primeSieve)
import Math.NumberTheory.Utils (splitOff)

-- | @'isCertifiedPrime' n@ tests primality of @n@, first trial division
--   by small primes is performed, then a Baillie PSW test and finally a
--   prime certificate is constructed and verified, provided no step before
--   found @n@ to be composite. Constructing prime certificates can take
--   a /very/ long time, so use this with care.
isCertifiedPrime :: Integer -> Bool
isCertifiedPrime :: Integer -> Bool
isCertifiedPrime Integer
n
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0     = Integer -> Bool
isCertifiedPrime (-Integer
n)
    | Bool
otherwise = Integer -> Bool
isPrime Integer
n Bool -> Bool -> Bool
&& ((Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
bpbd) Bool -> Bool -> Bool
|| PrimalityProof -> Bool
checkPrimalityProof (Integer -> PrimalityProof
certifyBPSW Integer
n))
      where
        bpbd :: Integer
bpbd = Integer
100000000000000000
-- Although it is known that there are no Baillie PSW pseudoprimes below 2^64,
-- use the verified bound 10^17, I don't know whether Gilchrist's result has been
-- verified yet.

-- | A proof of primality of a positive number. The type is
--   abstract to ensure the validity of proofs.
data PrimalityProof
    = Pocklington { PrimalityProof -> Integer
cprime :: !Integer          -- ^ The number whose primality is proved.
                  , PrimalityProof -> Integer
_factorisedPart, PrimalityProof -> Integer
_cofactor :: !Integer
                  , PrimalityProof -> [(Integer, Word, Integer, PrimalityProof)]
_knownFactors :: ![(Integer, Word, Integer, PrimalityProof)]
                  }
    | TrialDivision { cprime :: !Integer        -- ^ The number whose primality is proved.
                    , PrimalityProof -> Integer
_tdLimit :: !Integer }
    | Trivial { cprime :: !Integer              -- ^ The number whose primality is proved.
              }
      deriving Int -> PrimalityProof -> ShowS
[PrimalityProof] -> ShowS
PrimalityProof -> String
(Int -> PrimalityProof -> ShowS)
-> (PrimalityProof -> String)
-> ([PrimalityProof] -> ShowS)
-> Show PrimalityProof
forall a.
(Int -> a -> ShowS) -> (a -> String) -> ([a] -> ShowS) -> Show a
showList :: [PrimalityProof] -> ShowS
$cshowList :: [PrimalityProof] -> ShowS
show :: PrimalityProof -> String
$cshow :: PrimalityProof -> String
showsPrec :: Int -> PrimalityProof -> ShowS
$cshowsPrec :: Int -> PrimalityProof -> ShowS
Show

-- | Check the validity of a 'PrimalityProof'. Since it should be
--   impossible to create invalid proofs by the public interface, this
--   should never return 'False'.
checkPrimalityProof :: PrimalityProof -> Bool
checkPrimalityProof :: PrimalityProof -> Bool
checkPrimalityProof (Trivial Integer
n) = Integer -> Bool
isTrivialPrime Integer
n
checkPrimalityProof (TrialDivision Integer
p Integer
b) = Integer
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
bInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b Bool -> Bool -> Bool
&& Integer -> Integer -> Bool
trialDivisionPrimeTo Integer
b Integer
p
checkPrimalityProof (Pocklington Integer
p Integer
a Integer
b [(Integer, Word, Integer, PrimalityProof)]
fcts) = Integer
b Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
0 Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
b Bool -> Bool -> Bool
&& Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
b Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
pm1 Bool -> Bool -> Bool
&& Integer
a Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== [(Integer, Word, Integer, PrimalityProof)] -> Integer
forall b a c d. (Integral b, Num a) => [(a, b, c, d)] -> a
ppProd [(Integer, Word, Integer, PrimalityProof)]
fcts Bool -> Bool -> Bool
&& ((Integer, Word, Integer, PrimalityProof) -> Bool)
-> [(Integer, Word, Integer, PrimalityProof)] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer, Word, Integer, PrimalityProof) -> Bool
forall b. (Integer, b, Integer, PrimalityProof) -> Bool
verify [(Integer, Word, Integer, PrimalityProof)]
fcts
  where
    pm1 :: Integer
pm1 = Integer
pInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1
    ppProd :: [(a, b, c, d)] -> a
ppProd [(a, b, c, d)]
pps = [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
product [a
pfa -> b -> a
forall a b. (Num a, Integral b) => a -> b -> a
^b
e | (a
pf,b
e,c
_,d
_) <- [(a, b, c, d)]
pps]
    verify :: (Integer, b, Integer, PrimalityProof) -> Bool
verify (Integer
pf,b
_,Integer
base,PrimalityProof
proof) = Integer
pf Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== PrimalityProof -> Integer
cprime PrimalityProof
proof Bool -> Bool -> Bool
&& Integer -> Integer -> Bool
crit Integer
pf Integer
base Bool -> Bool -> Bool
&& PrimalityProof -> Bool
checkPrimalityProof PrimalityProof
proof
    crit :: Integer -> Integer -> Bool
crit Integer
pf Integer
base = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
p (Integer
xInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 Bool -> Bool -> Bool
&& Integer
y Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1
      where
        x :: Integer
x = Integer -> Integer -> Integer -> Integer
powModInteger Integer
base (Integer
pm1 Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
pf) Integer
p
        y :: Integer
y = Integer -> Integer -> Integer -> Integer
powModInteger Integer
x Integer
pf Integer
p

-- | @'isTrivialPrime'@ checks whether its argument is a trivially
--   known prime.
isTrivialPrime :: Integer -> Bool
isTrivialPrime :: Integer -> Bool
isTrivialPrime Integer
n = Integer
n Integer -> [Integer] -> Bool
forall (t :: * -> *) a. (Foldable t, Eq a) => a -> t a -> Bool
`elem` [Integer]
trivialPrimes

-- | List of trivially known primes.
trivialPrimes :: [Integer]
trivialPrimes :: [Integer]
trivialPrimes = [Integer
2,Integer
3,Integer
5,Integer
7,Integer
11,Integer
13,Integer
17,Integer
19,Integer
23,Integer
29]

-- | Certify a small number. This is not exposed and should only
--   be used where correct. It is always checked after use, though,
--   so it shouldn't be able to lie.
smallCert :: Integer -> PrimalityProof
smallCert :: Integer -> PrimalityProof
smallCert Integer
n
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
30    = Integer -> PrimalityProof
Trivial Integer
n
    | Bool
otherwise = Integer -> Integer -> PrimalityProof
TrialDivision Integer
n (Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
1)

-- | @'certify' n@ constructs, for @n > 1@, a proof of either
--   primality or compositeness of @n@. This may take a very long
--   time if the number has no small(ish) prime divisors
certify :: Integer -> Maybe PrimalityProof
certify :: Integer -> Maybe PrimalityProof
certify Integer
n
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
2     = String -> Maybe PrimalityProof
forall a. HasCallStack => String -> a
error String
"Only numbers larger than 1 can be certified"
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
31    = case [Integer] -> Integer -> [(Integer, Word)]
trialDivisionWith [Integer]
trivialPrimes Integer
n of
                    ((Integer
p,Word
_):[(Integer, Word)]
_) | Integer
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
n     -> Maybe PrimalityProof
forall a. Maybe a
Nothing
                              | Bool
otherwise -> PrimalityProof -> Maybe PrimalityProof
forall a. a -> Maybe a
Just (Integer -> PrimalityProof
Trivial Integer
n)
                    [(Integer, Word)]
_ -> String -> Maybe PrimalityProof
forall a. HasCallStack => String -> a
error String
"Impossible"
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
billi = let r2 :: Integer
r2 = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
n Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ Integer
2 in
                  case Integer -> Integer -> [(Integer, Word)]
trialDivisionTo Integer
r2 Integer
n of
                    ((Integer
p,Word
_):[(Integer, Word)]
_) | Integer
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
n       -> Maybe PrimalityProof
forall a. Maybe a
Nothing
                              | Bool
otherwise   -> PrimalityProof -> Maybe PrimalityProof
forall a. a -> Maybe a
Just (Integer -> Integer -> PrimalityProof
TrialDivision Integer
n Integer
r2)
                    [(Integer, Word)]
_ -> String -> Maybe PrimalityProof
forall a. HasCallStack => String -> a
error String
"Impossible"
    | Bool
otherwise = case Natural -> ([(Natural, Word)], Maybe Natural)
smallFactors (Integer -> Natural
forall a. Num a => Integer -> a
fromInteger (Integer -> Integer
forall a. Num a => a -> a
abs Integer
n)) of
                    ([], Just Natural
_) | Bool -> Bool
not (Integer -> Integer -> Bool
isStrongFermatPP Integer
n Integer
2) -> Maybe PrimalityProof
forall a. Maybe a
Nothing
                                 | Bool -> Bool
not (Integer -> Bool
lucasTest Integer
n) -> Maybe PrimalityProof
forall a. Maybe a
Nothing
                                 | Bool
otherwise -> PrimalityProof -> Maybe PrimalityProof
forall a. a -> Maybe a
Just (Integer -> PrimalityProof
certifyBPSW Integer
n)       -- if it isn't we error and ask for a report.
                    ((Natural -> Integer
forall a. Integral a => a -> Integer
toInteger -> Integer
p,Word
_):[(Natural, Word)]
_, Maybe Natural
_)
                      | Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
n -> PrimalityProof -> Maybe PrimalityProof
forall a. a -> Maybe a
Just (Integer -> Integer -> PrimalityProof
TrialDivision Integer
n (Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min Integer
100000 Integer
n))
                      | Bool
otherwise -> Maybe PrimalityProof
forall a. Maybe a
Nothing
                    ([(Natural, Word)], Maybe Natural)
_ -> String -> Maybe PrimalityProof
forall a. HasCallStack => String -> a
error (String
"***Error factorising " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"! Please report this to maintainer of arithmoi.")
      where
        billi :: Integer
billi = Integer
1000000000000

-- | Certify a number known to be not too small, having no small prime divisors and having
--   passed the Baillie PSW test. So we assume it's prime, erroring if not.
--   Since it's presumably a large number, we don't bother with trial division and
--   construct a Pocklington certificate.
certifyBPSW :: Integer -> PrimalityProof
certifyBPSW :: Integer -> PrimalityProof
certifyBPSW Integer
n = Integer
-> Integer
-> Integer
-> [(Integer, Word, Integer, PrimalityProof)]
-> PrimalityProof
Pocklington Integer
n Integer
a Integer
b [(Integer, Word, Integer, PrimalityProof)]
kfcts
  where
    nm1 :: Integer
nm1 = Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1
    h :: Integer
h = Integer
nm1 Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
2
    m3 :: Bool
m3 = Integer -> Int
forall a. Num a => Integer -> a
fromInteger Integer
n Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. (Int
3 :: Int) Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
3
    (Integer
a,[(Integer, Word, Bool)]
pp,Integer
b) = Integer -> (Integer, [(Integer, Word, Bool)], Integer)
findDecomposition Integer
nm1
    kfcts0 :: [(Integer, Word, Integer, PrimalityProof)]
kfcts0 = ((Integer, Word, Bool) -> (Integer, Word, Integer, PrimalityProof))
-> [(Integer, Word, Bool)]
-> [(Integer, Word, Integer, PrimalityProof)]
forall a b. (a -> b) -> [a] -> [b]
map (Integer, Word, Bool) -> (Integer, Word, Integer, PrimalityProof)
forall b.
(Integer, b, Bool) -> (Integer, b, Integer, PrimalityProof)
check [(Integer, Word, Bool)]
pp
    kfcts :: [(Integer, Word, Integer, PrimalityProof)]
kfcts = ([(Integer, Word, Integer, PrimalityProof)]
 -> (Integer, Word, Integer, PrimalityProof)
 -> [(Integer, Word, Integer, PrimalityProof)])
-> [(Integer, Word, Integer, PrimalityProof)]
-> [(Integer, Word, Integer, PrimalityProof)]
-> [(Integer, Word, Integer, PrimalityProof)]
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' [(Integer, Word, Integer, PrimalityProof)]
-> (Integer, Word, Integer, PrimalityProof)
-> [(Integer, Word, Integer, PrimalityProof)]
forall a b c a. [(a, b, c, a)] -> (a, b, c, a) -> [(a, b, c, a)]
force [] [(Integer, Word, Integer, PrimalityProof)]
kfcts0
    force :: [(a, b, c, a)] -> (a, b, c, a) -> [(a, b, c, a)]
force [(a, b, c, a)]
xs t :: (a, b, c, a)
t@(a
_,b
_,c
_,a
prf) = a
prf a -> [(a, b, c, a)] -> [(a, b, c, a)]
`seq` ((a, b, c, a)
t(a, b, c, a) -> [(a, b, c, a)] -> [(a, b, c, a)]
forall a. a -> [a] -> [a]
:[(a, b, c, a)]
xs)
    check :: (Integer, b, Bool) -> (Integer, b, Integer, PrimalityProof)
check (Integer
p,b
e,Bool
byTD) = Integer -> (Integer, b, Integer, PrimalityProof)
go Integer
2
      where
        go :: Integer -> (Integer, b, Integer, PrimalityProof)
go Integer
bs
            | Integer
bs Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
h    = String -> (Integer, b, Integer, PrimalityProof)
forall a. HasCallStack => String -> a
error (Integer -> String
bpswMessage Integer
n)
            | Integer
x Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1    = if Bool
m3 Bool -> Bool -> Bool
&& (Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
2) then (Integer
p,b
e,Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
bs,Integer -> PrimalityProof
Trivial Integer
2) else Integer -> (Integer, b, Integer, PrimalityProof)
go (Integer
bsInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1)
            | Integer
g Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
1    = String -> (Integer, b, Integer, PrimalityProof)
forall a. HasCallStack => String -> a
error (Integer -> String
bpswMessage Integer
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
found Integer
g)
            | Integer
y Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
1    = String -> (Integer, b, Integer, PrimalityProof)
forall a. HasCallStack => String -> a
error (Integer -> String
bpswMessage Integer
n String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
fermat Integer
bs)
            | Bool
byTD      = (Integer
p,b
e,Integer
bs, Integer -> PrimalityProof
smallCert Integer
p)
            | Bool
otherwise = case Integer -> Maybe PrimalityProof
certify Integer
p of
                            Maybe PrimalityProof
Nothing -> String -> (Integer, b, Integer, PrimalityProof)
forall a. HasCallStack => String -> a
error (String
"***Error in factorisation code: " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
p
                                                        String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" was supposed to be prime but isn't.\n"
                                                        String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
"Please report this to the maintainer.\n\n")
                            Just PrimalityProof
ppr ->(Integer
p,b
e,Integer
bs,PrimalityProof
ppr)
              where
                q :: Integer
q = Integer
nm1 Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
p
                x :: Integer
x = Integer -> Integer -> Integer -> Integer
powModInteger Integer
bs Integer
q Integer
n
                y :: Integer
y = Integer -> Integer -> Integer -> Integer
powModInteger Integer
x Integer
p Integer
n
                g :: Integer
g = Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n (Integer
xInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1)

-- | Find a decomposition of p-1 for the pocklington certificate.
--   Usually bloody slow if p-1 has two (or more) /large/ prime divisors.
findDecomposition :: Integer -> (Integer, [(Integer, Word, Bool)], Integer)
findDecomposition :: Integer -> (Integer, [(Integer, Word, Bool)], Integer)
findDecomposition Integer
n = Integer
-> Integer
-> [(Integer, Word, Bool)]
-> [Integer]
-> (Integer, [(Integer, Word, Bool)], Integer)
go Integer
1 Integer
n [] [Integer]
prms
  where
    sr :: Integer
sr = Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot Integer
n
    pbd :: Integer
pbd = Integer -> Integer -> Integer
forall a. Ord a => a -> a -> a
min Integer
1000000 (Integer
srInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
20)
    prms :: [Integer]
prms = (Prime Integer -> Integer) -> [Prime Integer] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Prime Integer -> Integer
forall a. Prime a -> a
unPrime ([Prime Integer] -> [Integer]) -> [Prime Integer] -> [Integer]
forall a b. (a -> b) -> a -> b
$ PrimeSieve -> [Prime Integer]
forall a. Integral a => PrimeSieve -> [Prime a]
primeList (Integer -> PrimeSieve
primeSieve Integer
pbd)
    go :: Integer
-> Integer
-> [(Integer, Word, Bool)]
-> [Integer]
-> (Integer, [(Integer, Word, Bool)], Integer)
go Integer
a Integer
b [(Integer, Word, Bool)]
afs (Integer
p:[Integer]
ps)
        | Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
b     = (Integer
a,[(Integer, Word, Bool)]
afs,Integer
b)
        | Bool
otherwise = case Integer -> Integer -> (Word, Integer)
forall a. (Eq a, GcdDomain a) => a -> a -> (Word, a)
splitOff Integer
p Integer
b of
                        (Word
0,Integer
_) -> Integer
-> Integer
-> [(Integer, Word, Bool)]
-> [Integer]
-> (Integer, [(Integer, Word, Bool)], Integer)
go Integer
a Integer
b [(Integer, Word, Bool)]
afs [Integer]
ps
                        (Word
e,Integer
q) -> Integer
-> Integer
-> [(Integer, Word, Bool)]
-> [Integer]
-> (Integer, [(Integer, Word, Bool)], Integer)
go (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
pInteger -> Word -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Word
e) Integer
q ((Integer
p,Word
e,Bool
True)(Integer, Word, Bool)
-> [(Integer, Word, Bool)] -> [(Integer, Word, Bool)]
forall a. a -> [a] -> [a]
:[(Integer, Word, Bool)]
afs) [Integer]
ps
    go Integer
a Integer
b [(Integer, Word, Bool)]
afs []
        | Integer
a Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
b     = (Integer
a,[(Integer, Word, Bool)]
afs,Integer
b)
        | Integer -> Bool
bailliePSW Integer
b  = (Integer
b,[(Integer
b,Word
1,Bool
False)],Integer
a)   -- Until a Baillie PSW pseudoprime is found, I'm going with this
        | Word
e Word -> Word -> Bool
forall a. Eq a => a -> a -> Bool
== Word
0    = String -> (Integer, [(Integer, Word, Bool)], Integer)
forall a. HasCallStack => String -> a
error (String
"Error in factorisation, " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
p String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" was found as a factor of " String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
b String -> ShowS
forall a. [a] -> [a] -> [a]
++ String
" but isn't.")
        | Bool
otherwise = Integer
-> Integer
-> [(Integer, Word, Bool)]
-> [Integer]
-> (Integer, [(Integer, Word, Bool)], Integer)
go (Integer
aInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
pInteger -> Word -> Integer
forall a b. (Num a, Integral b) => a -> b -> a
^Word
e) Integer
q ((Integer
p,Word
e,Bool
False)(Integer, Word, Bool)
-> [(Integer, Word, Bool)] -> [(Integer, Word, Bool)]
forall a. a -> [a] -> [a]
:[(Integer, Word, Bool)]
afs) []
          where
            p :: Integer
p = Integer -> Int -> Integer -> Integer
findFactor Integer
b Int
8 Integer
6
            (Word
e,Integer
q) = Integer -> Integer -> (Word, Integer)
forall a. (Eq a, GcdDomain a) => a -> a -> (Word, a)
splitOff Integer
p Integer
b

-- | Find a factor of a known composite with approximately digits digits,
--   starting with curve s. Actually, this may loop infinitely, but the
--   loop should not be entered before the heat death of the universe.
findFactor :: Integer -> Int -> Integer -> Integer
findFactor :: Integer -> Int -> Integer -> Integer
findFactor Integer
n Int
digits Integer
s = case Integer
-> Word -> Word -> Word -> Integer -> Either Integer Integer
findLoop Integer
n Word
lo Word
hi Word
count Integer
s of
                          Left Integer
t  -> Integer -> Int -> Integer -> Integer
findFactor Integer
n (Int
digitsInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
5) Integer
t
                          Right Integer
f -> Integer
f
  where
    (Word
lo,Word
hi,Word
count) = Int -> (Word, Word, Word)
findParms Int
digits

-- | Find a factor or say with which curve to continue.
findLoop :: Integer -> Word -> Word -> Word -> Integer -> Either Integer Integer
findLoop :: Integer
-> Word -> Word -> Word -> Integer -> Either Integer Integer
findLoop Integer
_ Word
_  Word
_  Word
0  Integer
s = Integer -> Either Integer Integer
forall a b. a -> Either a b
Left Integer
s
findLoop Integer
n Word
lo Word
hi Word
ct Integer
s
    | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
sInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
2  = Integer -> Either Integer Integer
forall a b. a -> Either a b
Left Integer
6
    | Bool
otherwise = case Natural -> SomeNat
someNatVal (Integer -> Natural
forall a. Num a => Integer -> a
fromInteger Integer
n) of
                    SomeNat (Proxy n
_ :: Proxy t) -> case Word -> Word -> Mod n -> Maybe Integer
forall (n :: Nat).
KnownNat n =>
Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation Word
lo Word
hi (Integer -> Mod n
forall a. Num a => Integer -> a
fromInteger Integer
s :: Mod t) of
                      Maybe Integer
Nothing -> Integer
-> Word -> Word -> Word -> Integer -> Either Integer Integer
findLoop Integer
n Word
lo Word
hi (Word
ctWord -> Word -> Word
forall a. Num a => a -> a -> a
-Word
1) (Integer
sInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1)
                      Just Integer
fct
                        | Integer -> Bool
bailliePSW Integer
fct -> Integer -> Either Integer Integer
forall a b. b -> Either a b
Right Integer
fct
                        | Bool
otherwise -> Integer -> Either Integer Integer
forall a b. b -> Either a b
Right (Integer -> Int -> Integer -> Integer
findFactor Integer
fct Int
8 (Integer
sInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1))

-- | Message in the unlikely case a Baillie PSW pseudoprime is found.
bpswMessage :: Integer -> String
bpswMessage :: Integer -> String
bpswMessage Integer
n = [String] -> String
unlines
                    [ String
"\n***Congratulations! You found a Baillie PSW pseudoprime!"
                    , String
"Please report this finding to the maintainers:"
                    , String
"<daniel.is.fischer@googlemail.com>,"
                    , String
"<andrew.lelechenko@gmail.com>"
                    , String
"The number in question is:\n"
                    , Integer -> String
forall a. Show a => a -> String
show Integer
n
                    , String
"\nOther parties like wikipedia might also be interested."
                    , String
"\nSorry for aborting your programm, but this is a major discovery."
                    ]

-- | Found a factor
found :: Integer -> String
found :: Integer -> String
found Integer
g = String
"\nA nontrivial divisor is:\n" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
g

-- | Fermat failure
fermat :: Integer -> String
fermat :: Integer -> String
fermat Integer
b = String
"\nThe Fermat test fails for base\n" String -> ShowS
forall a. [a] -> [a] -> [a]
++ Integer -> String
forall a. Show a => a -> String
show Integer
b