{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ScopedTypeVariables #-}
module Math.NumberTheory.Primes.Testing.Probabilistic
( isPrime
, millerRabinV
, bailliePSW
, isStrongFermatPP
, isFermatPP
, lucasTest
) where
import Data.Bits
import Data.Mod
import Data.Proxy
import GHC.Base
import GHC.Integer.GMP.Internals
import GHC.TypeNats (KnownNat, SomeNat(..), someNatVal)
import Math.NumberTheory.Moduli.JacobiSymbol
import Math.NumberTheory.Primes.Small
import Math.NumberTheory.Roots
import Math.NumberTheory.Utils
isPrime :: Integer -> Bool
isPrime :: Integer -> Bool
isPrime Integer
n
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = Integer -> Bool
isPrime (-Integer
n)
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
2 = Bool
False
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
4 = Bool
True
| Bool
otherwise = Int -> Integer -> Bool
millerRabinV Int
0 Integer
n
Bool -> Bool -> Bool
&& Integer -> Bool
bailliePSW Integer
n
millerRabinV :: Int -> Integer -> Bool
millerRabinV :: Int -> Integer -> Bool
millerRabinV Int
k Integer
n
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = Int -> Integer -> Bool
millerRabinV Int
k (-Integer
n)
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
2 = Bool
False
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
4 = Bool
True
| Bool
otherwise = [Integer] -> Bool
go [Integer]
smallPrimes
where
go :: [Integer] -> Bool
go (Integer
p:[Integer]
ps)
| Integer
pInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
p Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
> Integer
n = Bool
True
| Bool
otherwise = (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
p Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
/= Integer
0) Bool -> Bool -> Bool
&& [Integer] -> Bool
go [Integer]
ps
go [] = (Integer -> Bool) -> [Integer] -> Bool
forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer -> Integer -> Bool
isStrongFermatPP Integer
n) (Int -> [Integer] -> [Integer]
forall a. Int -> [a] -> [a]
take Int
k [Integer]
smallPrimes)
smallPrimes :: [Integer]
smallPrimes = (Word16 -> Integer) -> [Word16] -> [Integer]
forall a b. (a -> b) -> [a] -> [b]
map Word16 -> Integer
forall a. Integral a => a -> Integer
toInteger ([Word16] -> [Integer]) -> [Word16] -> [Integer]
forall a b. (a -> b) -> a -> b
$ Word16 -> Word16 -> [Word16]
smallPrimesFromTo Word16
forall a. Bounded a => a
minBound Word16
forall a. Bounded a => a
maxBound
isStrongFermatPP :: Integer -> Integer -> Bool
isStrongFermatPP :: Integer -> Integer -> Bool
isStrongFermatPP Integer
n Integer
b
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = [Char] -> Bool
forall a. HasCallStack => [Char] -> a
error [Char]
"isStrongFermatPP: negative argument"
| Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
1 = Bool
False
| Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
2 = Bool
True
| Bool
otherwise = case Natural -> SomeNat
someNatVal (Integer -> Natural
forall a. Num a => Integer -> a
fromInteger Integer
n) of
SomeNat (Proxy n
_ :: Proxy t) -> Mod n -> Bool
forall (n :: Nat). KnownNat n => Mod n -> Bool
isStrongFermatPPMod (Integer -> Mod n
forall a. Num a => Integer -> a
fromInteger Integer
b :: Mod t)
isStrongFermatPPMod :: KnownNat n => Mod n -> Bool
isStrongFermatPPMod :: Mod n -> Bool
isStrongFermatPPMod Mod n
b = Mod n
b Mod n -> Mod n -> Bool
forall a. Eq a => a -> a -> Bool
== Mod n
0 Bool -> Bool -> Bool
|| Mod n
a Mod n -> Mod n -> Bool
forall a. Eq a => a -> a -> Bool
== Mod n
1 Bool -> Bool -> Bool
|| Word -> Mod n -> Bool
forall t. (Eq t, Num t) => t -> Mod n -> Bool
go Word
t Mod n
a
where
m :: Mod n
m = -Mod n
1
(Word
t, Natural
u) = Natural -> (Word, Natural)
forall a. Integral a => a -> (Word, a)
shiftToOddCount (Natural -> (Word, Natural)) -> Natural -> (Word, Natural)
forall a b. (a -> b) -> a -> b
$ Mod n -> Natural
forall (m :: Nat). Mod m -> Natural
unMod Mod n
m
a :: Mod n
a = Mod n
b Mod n -> Natural -> Mod n
forall (m :: Nat) a.
(KnownNat m, Integral a) =>
Mod m -> a -> Mod m
^% Natural
u
go :: t -> Mod n -> Bool
go t
0 Mod n
_ = Bool
False
go t
k Mod n
x = Mod n
x Mod n -> Mod n -> Bool
forall a. Eq a => a -> a -> Bool
== Mod n
m Bool -> Bool -> Bool
|| t -> Mod n -> Bool
go (t
k t -> t -> t
forall a. Num a => a -> a -> a
- t
1) (Mod n
x Mod n -> Mod n -> Mod n
forall a. Num a => a -> a -> a
* Mod n
x)
isFermatPP :: Integer -> Integer -> Bool
isFermatPP :: Integer -> Integer -> Bool
isFermatPP Integer
n Integer
b = case Natural -> SomeNat
someNatVal (Integer -> Natural
forall a. Num a => Integer -> a
fromInteger Integer
n) of
SomeNat (Proxy n
_ :: Proxy t) -> (Integer -> Mod n
forall a. Num a => Integer -> a
fromInteger Integer
b :: Mod t) Mod n -> Integer -> Mod n
forall (m :: Nat) a.
(KnownNat m, Integral a) =>
Mod m -> a -> Mod m
^% (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Mod n -> Mod n -> Bool
forall a. Eq a => a -> a -> Bool
== Mod n
1
bailliePSW :: Integer -> Bool
bailliePSW :: Integer -> Bool
bailliePSW Integer
n = Integer -> Integer -> Bool
isStrongFermatPP Integer
n Integer
2 Bool -> Bool -> Bool
&& Integer -> Bool
lucasTest Integer
n
lucasTest :: Integer -> Bool
lucasTest :: Integer -> Bool
lucasTest Integer
n
| Integer -> Bool
forall a. Integral a => a -> Bool
isSquare Integer
n Bool -> Bool -> Bool
|| Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 = Bool
False
| Integer
d Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 = Bool
True
| Bool
otherwise = Integer
uo Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
|| Word -> Integer -> Integer -> Bool
forall t. (Eq t, Num t) => t -> Integer -> Integer -> Bool
go Word
t Integer
vo Integer
qo
where
d :: Integer
d = Bool -> Integer -> Integer
find Bool
True Integer
5
find :: Bool -> Integer -> Integer
find !Bool
pos Integer
cd = case Integer -> Integer -> JacobiSymbol
forall a. (Integral a, Bits a) => a -> a -> JacobiSymbol
jacobi (Integer
n Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
cd) Integer
cd of
JacobiSymbol
MinusOne -> if Bool
pos then Integer
cd else (-Integer
cd)
JacobiSymbol
Zero -> if Integer
cd Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
n then Integer
1 else Integer
0
JacobiSymbol
One -> Bool -> Integer -> Integer
find (Bool -> Bool
not Bool
pos) (Integer
cdInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
2)
q :: Integer
q = (Integer
1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
d) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
4
(Word
t,Integer
o) = Integer -> (Word, Integer)
forall a. Integral a => a -> (Word, a)
shiftToOddCount (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1)
(Integer
uo, Integer
vo, Integer
qo) = Integer -> Integer -> Integer -> (Integer, Integer, Integer)
testLucas Integer
n Integer
q Integer
o
go :: t -> Integer -> Integer -> Bool
go t
0 Integer
_ Integer
_ = Bool
False
go t
s Integer
vn Integer
qn = Integer
vn Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
|| t -> Integer -> Integer -> Bool
go (t
st -> t -> t
forall a. Num a => a -> a -> a
-t
1) ((Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n) ((Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n)
testLucas :: Integer -> Integer -> Integer -> (Integer, Integer, Integer)
testLucas :: Integer -> Integer -> Integer -> (Integer, Integer, Integer)
testLucas Integer
n Integer
q (S# Int#
i#) = Int -> (Integer, Integer, Integer)
look (Word -> Int
forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
2)
where
j :: Int
j = Int# -> Int
I# Int#
i#
look :: Int -> (Integer, Integer, Integer)
look Int
k
| Int -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Int
j Int
k = Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer
1 Integer
1 Integer
1 Integer
q
| Bool
otherwise = Int -> (Integer, Integer, Integer)
look (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
go :: Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go Int
k Integer
un Integer
un1 Integer
vn Integer
qn
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = (Integer
un, Integer
vn, Integer
qn)
| Int -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Int
j Int
k = Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer
u2n1 Integer
u2n2 Integer
v2n1 Integer
q2n1
| Bool
otherwise = Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go (Int
kInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer
u2n Integer
u2n1 Integer
v2n Integer
q2n
where
u2n :: Integer
u2n = (Integer
unInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n1 :: Integer
u2n1 = (Integer
un1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n2 :: Integer
u2n2 = ((Integer
un1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
un)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n :: Integer
v2n = (Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n1 :: Integer
v2n1 = ((Integer
un1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- (Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
q)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
un)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n :: Integer
q2n = (Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n1 :: Integer
q2n1 = (Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
q) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
testLucas Integer
n Integer
q (Jp# BigNat
bn#) = Int# -> (Integer, Integer, Integer)
test (Int#
s# Int# -> Int# -> Int#
-# Int#
1#)
where
s# :: Int#
s# = BigNat -> Int#
sizeofBigNat# BigNat
bn#
test :: Int# -> (Integer, Integer, Integer)
test Int#
j# = case BigNat -> Int# -> GmpLimb#
indexBigNat# BigNat
bn# Int#
j# of
GmpLimb#
0## -> Int# -> (Integer, Integer, Integer)
test (Int#
j# Int# -> Int# -> Int#
-# Int#
1#)
GmpLimb#
w# -> Int# -> Word -> Int -> (Integer, Integer, Integer)
look (Int#
j# Int# -> Int# -> Int#
-# Int#
1#) (GmpLimb# -> Word
W# GmpLimb#
w#) (Word -> Int
forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
look :: Int# -> Word -> Int -> (Integer, Integer, Integer)
look Int#
j# Word
w Int
i
| Word -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Word
w Int
i = Int#
-> Word
-> Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go Int#
j# Word
w (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Integer
1 Integer
1 Integer
1 Integer
q
| Bool
otherwise = Int# -> Word -> Int -> (Integer, Integer, Integer)
look Int#
j# Word
w (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
go :: Int#
-> Word
-> Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go Int#
k# Word
w Int
i Integer
un Integer
un1 Integer
vn Integer
qn
| Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = if Int# -> Bool
isTrue# (Int#
k# Int# -> Int# -> Int#
<# Int#
0#)
then (Integer
un,Integer
vn,Integer
qn)
else Int#
-> Word
-> Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go (Int#
k# Int# -> Int# -> Int#
-# Int#
1#) (GmpLimb# -> Word
W# (BigNat -> Int# -> GmpLimb#
indexBigNat# BigNat
bn# Int#
k#)) (Word -> Int
forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) Integer
un Integer
un1 Integer
vn Integer
qn
| Word -> Int -> Bool
forall a. Bits a => a -> Int -> Bool
testBit Word
w Int
i = Int#
-> Word
-> Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go Int#
k# Word
w (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer
u2n1 Integer
u2n2 Integer
v2n1 Integer
q2n1
| Bool
otherwise = Int#
-> Word
-> Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go Int#
k# Word
w (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer
u2n Integer
u2n1 Integer
v2n Integer
q2n
where
u2n :: Integer
u2n = (Integer
unInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n1 :: Integer
u2n1 = (Integer
un1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n2 :: Integer
u2n2 = ((Integer
un1Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
un)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n :: Integer
v2n = (Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n1 :: Integer
v2n1 = ((Integer
un1 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- (Integer
2Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
q)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
un)Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
vnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n :: Integer
q2n = (Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qn) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n1 :: Integer
q2n1 = (Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
qnInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
*Integer
q) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
testLucas Integer
_ Integer
_ Integer
_ = [Char] -> (Integer, Integer, Integer)
forall a. HasCallStack => [Char] -> a
error [Char]
"lucasTest: negative argument"