{-# 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.Exts (Word(..), Int(..), (-#), (<#), isTrue#)
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 forall a. Ord a => a -> a -> Bool
< Integer
0 = Integer -> Bool
isPrime (-Integer
n)
| Integer
n forall a. Ord a => a -> a -> Bool
< Integer
2 = Bool
False
| Integer
n 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 forall a. Ord a => a -> a -> Bool
< Integer
0 = Int -> Integer -> Bool
millerRabinV Int
k (-Integer
n)
| Integer
n forall a. Ord a => a -> a -> Bool
< Integer
2 = Bool
False
| Integer
n 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
pforall a. Num a => a -> a -> a
*Integer
p forall a. Ord a => a -> a -> Bool
> Integer
n = Bool
True
| Bool
otherwise = (Integer
n forall a. Integral a => a -> a -> a
`rem` Integer
p forall a. Eq a => a -> a -> Bool
/= Integer
0) Bool -> Bool -> Bool
&& [Integer] -> Bool
go [Integer]
ps
go [] = forall (t :: * -> *) a. Foldable t => (a -> Bool) -> t a -> Bool
all (Integer -> Integer -> Bool
isStrongFermatPP Integer
n) (forall a. Int -> [a] -> [a]
take Int
k [Integer]
smallPrimes)
smallPrimes :: [Integer]
smallPrimes = forall a b. (a -> b) -> [a] -> [b]
map forall a. Integral a => a -> Integer
toInteger forall a b. (a -> b) -> a -> b
$ Word16 -> Word16 -> [Word16]
smallPrimesFromTo forall a. Bounded a => a
minBound forall a. Bounded a => a
maxBound
isStrongFermatPP :: Integer -> Integer -> Bool
isStrongFermatPP :: Integer -> Integer -> Bool
isStrongFermatPP Integer
n Integer
b
| Integer
n forall a. Ord a => a -> a -> Bool
< Integer
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"isStrongFermatPP: negative argument"
| Integer
n forall a. Ord a => a -> a -> Bool
<= Integer
1 = Bool
False
| Integer
n forall a. Eq a => a -> a -> Bool
== Integer
2 = Bool
True
| Bool
otherwise = case Natural -> SomeNat
someNatVal (forall a. Num a => Integer -> a
fromInteger Integer
n) of
SomeNat (Proxy n
_ :: Proxy t) -> forall (n :: Natural). KnownNat n => Mod n -> Bool
isStrongFermatPPMod (forall a. Num a => Integer -> a
fromInteger Integer
b :: Mod t)
isStrongFermatPPMod :: KnownNat n => Mod n -> Bool
isStrongFermatPPMod :: forall (n :: Natural). KnownNat n => Mod n -> Bool
isStrongFermatPPMod Mod n
b = Mod n
b forall a. Eq a => a -> a -> Bool
== Mod n
0 Bool -> Bool -> Bool
|| Mod n
a forall a. Eq a => a -> a -> Bool
== Mod n
1 Bool -> Bool -> 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) = forall a. Integral a => a -> (Word, a)
shiftToOddCount forall a b. (a -> b) -> a -> b
$ forall (m :: Natural). Mod m -> Natural
unMod Mod n
m
a :: Mod n
a = Mod n
b forall (m :: Natural) 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 forall a. Eq a => a -> a -> Bool
== Mod n
m Bool -> Bool -> Bool
|| t -> Mod n -> Bool
go (t
k forall a. Num a => a -> a -> a
- t
1) (Mod n
x 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 (forall a. Num a => Integer -> a
fromInteger Integer
n) of
SomeNat (Proxy n
_ :: Proxy t) -> (forall a. Num a => Integer -> a
fromInteger Integer
b :: Mod t) forall (m :: Natural) a.
(KnownNat m, Integral a) =>
Mod m -> a -> Mod m
^% (Integer
nforall a. Num a => a -> a -> a
-Integer
1) 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
| forall a. Integral a => a -> Bool
isSquare Integer
n Bool -> Bool -> Bool
|| Integer
d forall a. Eq a => a -> a -> Bool
== Integer
0 = Bool
False
| Integer
d forall a. Eq a => a -> a -> Bool
== Integer
1 = Bool
True
| Bool
otherwise = Integer
uo forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> 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 forall a. (Integral a, Bits a) => a -> a -> JacobiSymbol
jacobi (Integer
n 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 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
cdforall a. Num a => a -> a -> a
+Integer
2)
q :: Integer
q = (Integer
1forall a. Num a => a -> a -> a
-Integer
d) forall a. Integral a => a -> a -> a
`quot` Integer
4
(Word
t,Integer
o) = forall a. Integral a => a -> (Word, a)
shiftToOddCount (Integer
nforall 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 forall a. Eq a => a -> a -> Bool
== Integer
0 Bool -> Bool -> Bool
|| t -> Integer -> Integer -> Bool
go (t
sforall a. Num a => a -> a -> a
-t
1) ((Integer
vnforall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
2forall a. Num a => a -> a -> a
*Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n) ((Integer
qnforall a. Num a => a -> a -> a
*Integer
qn) 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 (forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) 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
| forall a. Bits a => a -> Int -> Bool
testBit Int
j Int
k = Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go (Int
kforall a. Num a => a -> a -> a
-Int
1) Integer
1 Integer
1 Integer
1 Integer
q
| Bool
otherwise = Int -> (Integer, Integer, Integer)
look (Int
kforall 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 forall a. Ord a => a -> a -> Bool
< Int
0 = (Integer
un, Integer
vn, Integer
qn)
| forall a. Bits a => a -> Int -> Bool
testBit Int
j Int
k = Int
-> Integer
-> Integer
-> Integer
-> Integer
-> (Integer, Integer, Integer)
go (Int
kforall 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
kforall a. Num a => a -> a -> a
-Int
1) Integer
u2n Integer
u2n1 Integer
v2n Integer
q2n
where
u2n :: Integer
u2n = (Integer
unforall a. Num a => a -> a -> a
*Integer
vn) forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n1 :: Integer
u2n1 = (Integer
un1forall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n2 :: Integer
u2n2 = ((Integer
un1forall a. Num a => a -> a -> a
-Integer
qforall a. Num a => a -> a -> a
*Integer
un)forall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n :: Integer
v2n = (Integer
vnforall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
2forall a. Num a => a -> a -> a
*Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n1 :: Integer
v2n1 = ((Integer
un1 forall a. Num a => a -> a -> a
- (Integer
2forall a. Num a => a -> a -> a
*Integer
q)forall a. Num a => a -> a -> a
*Integer
un)forall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n :: Integer
q2n = (Integer
qnforall a. Num a => a -> a -> a
*Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n1 :: Integer
q2n1 = (Integer
qnforall a. Num a => a -> a -> a
*Integer
qnforall a. Num a => a -> a -> a
*Integer
q) 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#) (forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) forall a. Num a => a -> a -> a
- Int
1)
look :: Int# -> Word -> Int -> (Integer, Integer, Integer)
look Int#
j# Word
w Int
i
| 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 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
iforall 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 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#)) (forall b. FiniteBits b => b -> Int
finiteBitSize (Word
0 :: Word) forall a. Num a => a -> a -> a
- Int
1) Integer
un Integer
un1 Integer
vn Integer
qn
| 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
iforall 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
iforall a. Num a => a -> a -> a
-Int
1) Integer
u2n Integer
u2n1 Integer
v2n Integer
q2n
where
u2n :: Integer
u2n = (Integer
unforall a. Num a => a -> a -> a
*Integer
vn) forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n1 :: Integer
u2n1 = (Integer
un1forall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
u2n2 :: Integer
u2n2 = ((Integer
un1forall a. Num a => a -> a -> a
-Integer
qforall a. Num a => a -> a -> a
*Integer
un)forall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n :: Integer
v2n = (Integer
vnforall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
2forall a. Num a => a -> a -> a
*Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
v2n1 :: Integer
v2n1 = ((Integer
un1 forall a. Num a => a -> a -> a
- (Integer
2forall a. Num a => a -> a -> a
*Integer
q)forall a. Num a => a -> a -> a
*Integer
un)forall a. Num a => a -> a -> a
*Integer
vnforall a. Num a => a -> a -> a
-Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n :: Integer
q2n = (Integer
qnforall a. Num a => a -> a -> a
*Integer
qn) forall a. Integral a => a -> a -> a
`rem` Integer
n
q2n1 :: Integer
q2n1 = (Integer
qnforall a. Num a => a -> a -> a
*Integer
qnforall a. Num a => a -> a -> a
*Integer
q) forall a. Integral a => a -> a -> a
`rem` Integer
n
testLucas Integer
_ Integer
_ Integer
_ = forall a. HasCallStack => [Char] -> a
error [Char]
"lucasTest: negative argument"