module Factor.Prime
where
import qualified Data.List as List
import qualified Data.Set as Set
import System.Random (RandomGen)
import qualified System.Random as Random
import Factor.Util
type Prime = Integer
primes :: [Prime]
primes :: [Prime]
primes = Prime
2 Prime -> [Prime] -> [Prime]
forall a. a -> [a] -> [a]
: Prime
3 Prime -> [Prime] -> [Prime]
forall a. a -> [a] -> [a]
: Prime
5 Prime -> [Prime] -> [Prime]
forall a. a -> [a] -> [a]
: Prime
-> ((Prime, Prime), Set (Prime, Prime))
-> (Prime, Prime, [Prime])
-> [Prime]
forall a.
(Num a, Ord a) =>
a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
sieveP Prime
7 ((Prime
9,Prime
6),Set (Prime, Prime)
forall a. Set a
Set.empty) ([Prime] -> (Prime, Prime, [Prime])
forall b. Num b => [b] -> (b, b, [b])
nextQ (Int -> [Prime] -> [Prime]
forall a. Int -> [a] -> [a]
drop Int
2 [Prime]
primes))
where
sieveP :: a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
sieveP a
n ((a
pk,a
p),Set (a, a)
ps) (a, a, [a])
qs =
case a -> a -> Ordering
forall a. Ord a => a -> a -> Ordering
compare a
n a
pk of
Ordering
LT -> a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
checkQ a
n ((a
pk,a
p),Set (a, a)
ps) (a, a, [a])
qs
Ordering
EQ -> a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
sieveP (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
2) (a -> a -> Set (a, a) -> ((a, a), Set (a, a))
forall b.
(Ord b, Num b) =>
b -> b -> Set (b, b) -> ((b, b), Set (b, b))
incrP a
pk a
p Set (a, a)
ps) (a, a, [a])
qs
Ordering
GT -> a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
sieveP a
n (a -> a -> Set (a, a) -> ((a, a), Set (a, a))
forall b.
(Ord b, Num b) =>
b -> b -> Set (b, b) -> ((b, b), Set (b, b))
incrP a
pk a
p Set (a, a)
ps) (a, a, [a])
qs
checkQ :: a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
checkQ a
n ((a, a)
pk,Set (a, a)
ps) (a
q2,a
q,[a]
qs) | a
n a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
q2 = a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
sieveP (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
2) ((a, a), Set (a, a))
pks ([a] -> (a, a, [a])
forall b. Num b => [b] -> (b, b, [b])
nextQ [a]
qs)
where pks :: ((a, a), Set (a, a))
pks = a -> a -> Set (a, a) -> ((a, a), Set (a, a))
forall b.
(Ord b, Num b) =>
b -> b -> Set (b, b) -> ((b, b), Set (b, b))
incrP a
q2 (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
q) ((a, a) -> Set (a, a) -> Set (a, a)
forall a. Ord a => a -> Set a -> Set a
Set.insert (a, a)
pk Set (a, a)
ps)
checkQ a
n ((a, a), Set (a, a))
ps (a, a, [a])
qs = a
n a -> [a] -> [a]
forall a. a -> [a] -> [a]
: a -> ((a, a), Set (a, a)) -> (a, a, [a]) -> [a]
sieveP (a
na -> a -> a
forall a. Num a => a -> a -> a
+a
2) ((a, a), Set (a, a))
ps (a, a, [a])
qs
incrP :: b -> b -> Set (b, b) -> ((b, b), Set (b, b))
incrP b
pk b
p Set (b, b)
s = Set (b, b) -> ((b, b), Set (b, b))
forall a. Set a -> (a, Set a)
Set.deleteFindMin ((b, b) -> Set (b, b) -> Set (b, b)
forall a. Ord a => a -> Set a -> Set a
Set.insert (b
pkb -> b -> b
forall a. Num a => a -> a -> a
+b
p, b
p) Set (b, b)
s)
nextQ :: [b] -> (b, b, [b])
nextQ [] = [Char] -> (b, b, [b])
forall a. HasCallStack => [Char] -> a
error [Char]
"ran out of primes"
nextQ (b
q : [b]
qs) = (b
qb -> b -> b
forall a. Num a => a -> a -> a
*b
q, b
q, [b]
qs)
type PrimePower = (Prime,Integer)
multiplyPrimePowers :: [PrimePower] -> [PrimePower] -> [PrimePower]
multiplyPrimePowers :: [(Prime, Prime)] -> [(Prime, Prime)] -> [(Prime, Prime)]
multiplyPrimePowers [(Prime, Prime)]
pks1 [] = [(Prime, Prime)]
pks1
multiplyPrimePowers [] [(Prime, Prime)]
pks2 = [(Prime, Prime)]
pks2
multiplyPrimePowers ((Prime
p1,Prime
k1) : [(Prime, Prime)]
pks1) ((Prime
p2,Prime
k2) : [(Prime, Prime)]
pks2) =
case Prime -> Prime -> Ordering
forall a. Ord a => a -> a -> Ordering
compare Prime
p1 Prime
p2 of
Ordering
LT -> (Prime
p1,Prime
k1) (Prime, Prime) -> [(Prime, Prime)] -> [(Prime, Prime)]
forall a. a -> [a] -> [a]
: [(Prime, Prime)] -> [(Prime, Prime)] -> [(Prime, Prime)]
multiplyPrimePowers [(Prime, Prime)]
pks1 ((Prime
p2,Prime
k2) (Prime, Prime) -> [(Prime, Prime)] -> [(Prime, Prime)]
forall a. a -> [a] -> [a]
: [(Prime, Prime)]
pks2)
Ordering
EQ -> (Prime
p1, Prime
k1 Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
k2) (Prime, Prime) -> [(Prime, Prime)] -> [(Prime, Prime)]
forall a. a -> [a] -> [a]
: [(Prime, Prime)] -> [(Prime, Prime)] -> [(Prime, Prime)]
multiplyPrimePowers [(Prime, Prime)]
pks1 [(Prime, Prime)]
pks2
Ordering
GT -> (Prime
p2,Prime
k2) (Prime, Prime) -> [(Prime, Prime)] -> [(Prime, Prime)]
forall a. a -> [a] -> [a]
: [(Prime, Prime)] -> [(Prime, Prime)] -> [(Prime, Prime)]
multiplyPrimePowers ((Prime
p1,Prime
k1) (Prime, Prime) -> [(Prime, Prime)] -> [(Prime, Prime)]
forall a. a -> [a] -> [a]
: [(Prime, Prime)]
pks1) [(Prime, Prime)]
pks2
productPrimePowers :: [[PrimePower]] -> [PrimePower]
productPrimePowers :: [[(Prime, Prime)]] -> [(Prime, Prime)]
productPrimePowers [] = []
productPrimePowers [[(Prime, Prime)]
pks] = [(Prime, Prime)]
pks
productPrimePowers [[(Prime, Prime)]]
pksl = [(Prime, Prime)] -> [(Prime, Prime)] -> [(Prime, Prime)]
multiplyPrimePowers [(Prime, Prime)]
pks1 [(Prime, Prime)]
pks2
where
pks1 :: [(Prime, Prime)]
pks1 = [[(Prime, Prime)]] -> [(Prime, Prime)]
productPrimePowers [[(Prime, Prime)]]
pksl1
pks2 :: [(Prime, Prime)]
pks2 = [[(Prime, Prime)]] -> [(Prime, Prime)]
productPrimePowers [[(Prime, Prime)]]
pksl2
([[(Prime, Prime)]]
pksl1,[[(Prime, Prime)]]
pksl2) = Int
-> [[(Prime, Prime)]] -> ([[(Prime, Prime)]], [[(Prime, Prime)]])
forall a. Int -> [a] -> ([a], [a])
List.splitAt ([[(Prime, Prime)]] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length [[(Prime, Prime)]]
pksl Int -> Int -> Int
forall a. Integral a => a -> a -> a
`div` Int
2) [[(Prime, Prime)]]
pksl
factor :: [Prime] -> Integer -> ([PrimePower],Integer)
factor :: [Prime] -> Prime -> ([(Prime, Prime)], Prime)
factor [Prime]
_ Prime
0 = ([],Prime
0)
factor [Prime]
ps Prime
n | Prime
n Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
< Prime
0 = ([(Prime, Prime)]
pks, Prime -> Prime
forall a. Num a => a -> a
Prelude.negate Prime
s)
where ([(Prime, Prime)]
pks,Prime
s) = [Prime] -> Prime -> ([(Prime, Prime)], Prime)
factor [Prime]
ps (Prime -> Prime
forall a. Num a => a -> a
Prelude.negate Prime
n)
factor [Prime]
ps0 Prime
n0 | Bool
otherwise = [Prime] -> Prime -> ([(Prime, Prime)], Prime)
go [Prime]
ps0 Prime
n0
where
go :: [Prime] -> Prime -> ([(Prime, Prime)], Prime)
go [Prime]
_ Prime
1 = ([],Prime
1)
go [] Prime
n = ([],Prime
n)
go (Prime
p : [Prime]
ps) Prime
n = (if Prime
k Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
0 then [(Prime, Prime)]
pks else (Prime
p,Prime
k) (Prime, Prime) -> [(Prime, Prime)] -> [(Prime, Prime)]
forall a. a -> [a] -> [a]
: [(Prime, Prime)]
pks, Prime
s)
where
([(Prime, Prime)]
pks,Prime
s) = [Prime] -> Prime -> ([(Prime, Prime)], Prime)
go [Prime]
ps Prime
m
(Prime
k,Prime
m) = Prime -> Prime -> (Prime, Prime)
divPower Prime
p Prime
n
factorProduct :: [Prime] -> [Integer] -> ([PrimePower],[Integer])
factorProduct :: [Prime] -> [Prime] -> ([(Prime, Prime)], [Prime])
factorProduct [Prime]
ps [Prime]
nl = ([[(Prime, Prime)]] -> [(Prime, Prime)]
productPrimePowers [[(Prime, Prime)]]
pksl, [Prime]
sl)
where ([[(Prime, Prime)]]
pksl,[Prime]
sl) = [([(Prime, Prime)], Prime)] -> ([[(Prime, Prime)]], [Prime])
forall a b. [(a, b)] -> ([a], [b])
unzip ([([(Prime, Prime)], Prime)] -> ([[(Prime, Prime)]], [Prime]))
-> [([(Prime, Prime)], Prime)] -> ([[(Prime, Prime)]], [Prime])
forall a b. (a -> b) -> a -> b
$ (Prime -> ([(Prime, Prime)], Prime))
-> [Prime] -> [([(Prime, Prime)], Prime)]
forall a b. (a -> b) -> [a] -> [b]
map ([Prime] -> Prime -> ([(Prime, Prime)], Prime)
factor [Prime]
ps) [Prime]
nl
trialDivision :: Integer -> ([PrimePower],Integer)
trialDivision :: Prime -> ([(Prime, Prime)], Prime)
trialDivision = [Prime] -> Prime -> ([(Prime, Prime)], Prime)
factor [Prime]
primes
primesUnder :: Integer -> Integer
primesUnder :: Prime -> Prime
primesUnder Prime
n = Int -> Prime
forall a. Integral a => a -> Prime
toInteger ([Prime] -> Int
forall (t :: * -> *) a. Foldable t => t a -> Int
length ((Prime -> Bool) -> [Prime] -> [Prime]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (\Prime
p -> Prime
p Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
<= Prime
n) [Prime]
primes))
primesUnderEstimate :: Log2Integer -> Log2Integer
primesUnderEstimate :: Log2Integer -> Log2Integer
primesUnderEstimate = Log2Integer -> Log2Integer
go
where
go :: Log2Integer -> Log2Integer
go Log2Integer
log2n | Log2Integer
log2n Log2Integer -> Log2Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Log2Integer
t = Log2Integer -> Log2Integer
f Log2Integer
log2n
go Log2Integer
log2n = Log2Integer -> Log2Integer -> Log2Integer
forall a. Ord a => a -> a -> a
max Log2Integer
f_t (Log2Integer
log2e Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Log2Integer
log2n Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer -> Log2Integer
log2 Log2Integer
log2n)
f :: Log2Integer -> Log2Integer
f = Prime -> Log2Integer
log2Integer (Prime -> Log2Integer)
-> (Log2Integer -> Prime) -> Log2Integer -> Log2Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Prime -> Prime
primesUnder (Prime -> Prime) -> (Log2Integer -> Prime) -> Log2Integer -> Prime
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Log2Integer -> Prime
exp2Integer
f_t :: Log2Integer
f_t = Log2Integer -> Log2Integer
f Log2Integer
t
t :: Log2Integer
t = Log2Integer
5.0
nthPrime :: Integer -> Integer
nthPrime :: Prime -> Prime
nthPrime Prime
n = [Prime]
primes [Prime] -> Int -> Prime
forall a. [a] -> Int -> a
!! (Prime -> Int
forall a. Num a => Prime -> a
Prelude.fromInteger Prime
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1)
nthPrimeEstimate :: Log2Integer -> Log2Integer
nthPrimeEstimate :: Log2Integer -> Log2Integer
nthPrimeEstimate = Log2Integer -> Log2Integer
go
where
go :: Log2Integer -> Log2Integer
go Log2Integer
log2n | Log2Integer
log2n Log2Integer -> Log2Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Log2Integer
t = Log2Integer -> Log2Integer
f Log2Integer
log2n
go Log2Integer
log2n = Log2Integer -> Log2Integer -> Log2Integer
forall a. Ord a => a -> a -> a
max Log2Integer
f_t (Log2Integer
log2n Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Log2Integer -> Log2Integer
log2 Log2Integer
log2n Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
log2e)
f :: Log2Integer -> Log2Integer
f = Prime -> Log2Integer
log2Integer (Prime -> Log2Integer)
-> (Log2Integer -> Prime) -> Log2Integer -> Log2Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Prime -> Prime
nthPrime (Prime -> Prime) -> (Log2Integer -> Prime) -> Log2Integer -> Prime
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Log2Integer -> Prime
exp2Integer
f_t :: Log2Integer
f_t = Log2Integer -> Log2Integer
f Log2Integer
t
t :: Log2Integer
t = Log2Integer
5.0
smoothUnderLowerBound :: Integer -> Log2Integer -> Log2Integer
smoothUnderLowerBound :: Prime -> Log2Integer -> Log2Integer
smoothUnderLowerBound = (Log2Integer -> Log2Integer)
-> [(Prime, Log2Integer -> Log2Integer)]
-> Prime
-> Log2Integer
-> Log2Integer
forall t t. Ord t => t -> [(t, t)] -> t -> t
go (Log2Integer -> Log2Integer -> Log2Integer
forall a b. a -> b -> a
const Log2Integer
0.0) [(Prime, Log2Integer -> Log2Integer)]
bnds
where
go :: t -> [(t, t)] -> t -> t
go t
_ [] t
_ = [Char] -> t
forall a. HasCallStack => [Char] -> a
error [Char]
"ran out of bounds"
go t
f' ((t
p,t
f) : [(t, t)]
pfs) t
b = if t
b t -> t -> Bool
forall a. Ord a => a -> a -> Bool
< t
p then t
f' else t -> [(t, t)] -> t -> t
go t
f [(t, t)]
pfs t
b
bnds :: [(Prime, Log2Integer -> Log2Integer)]
bnds = Log2Integer
-> Log2Integer -> [Prime] -> [(Prime, Log2Integer -> Log2Integer)]
bndp Log2Integer
0.0 Log2Integer
0.0 [Prime]
primes
bndp :: Log2Integer
-> Log2Integer -> [Prime] -> [(Prime, Log2Integer -> Log2Integer)]
bndp Log2Integer
_ Log2Integer
_ [] = [Char] -> [(Prime, Log2Integer -> Log2Integer)]
forall a. HasCallStack => [Char] -> a
error [Char]
"ran out of primes"
bndp Log2Integer
i Log2Integer
w (Prime
p : [Prime]
ps) = (Prime
p,Log2Integer -> Log2Integer
f) (Prime, Log2Integer -> Log2Integer)
-> [(Prime, Log2Integer -> Log2Integer)]
-> [(Prime, Log2Integer -> Log2Integer)]
forall a. a -> [a] -> [a]
: Log2Integer
-> Log2Integer -> [Prime] -> [(Prime, Log2Integer -> Log2Integer)]
bndp Log2Integer
i' Log2Integer
w' [Prime]
ps
where
f :: Log2Integer -> Log2Integer
f Log2Integer
log2n = Log2Integer
i' Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
* Log2Integer -> Log2Integer
log2 Log2Integer
log2n Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
w'
i' :: Log2Integer
i' = Log2Integer
i Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Log2Integer
1.0
w' :: Log2Integer
w' = Log2Integer
w Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Log2Integer -> Log2Integer
log2 Log2Integer
i' Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Log2Integer -> Log2Integer
log2 (Prime -> Log2Integer
log2Integer Prime
p)
smoothUnder :: Integer -> Log2Integer -> Log2Integer
smoothUnder :: Prime -> Log2Integer -> Log2Integer
smoothUnder Prime
b Log2Integer
log2n =
if Log2Integer
log2n Log2Integer -> Log2Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Log2Integer
log2b then Log2Integer
log2n else Log2Integer -> Log2Integer -> Log2Integer
forall a. Ord a => a -> a -> a
max Log2Integer
lwr (Log2Integer
log2n Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
u Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
* Log2Integer -> Log2Integer
log2 Log2Integer
u)
where
lwr :: Log2Integer
lwr = (if Log2Integer
log2n Log2Integer -> Log2Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Prime -> Log2Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral Prime
b then Log2Integer -> Log2Integer
forall a. a -> a
id
else Log2Integer -> Log2Integer -> Log2Integer
forall a. Ord a => a -> a -> a
max (Prime -> Log2Integer -> Log2Integer
smoothUnderLowerBound Prime
b Log2Integer
log2n)) Log2Integer
log2b
u :: Log2Integer
u = Log2Integer
log2n Log2Integer -> Log2Integer -> Log2Integer
forall a. Fractional a => a -> a -> a
/ Log2Integer
log2b
log2b :: Log2Integer
log2b = Prime -> Log2Integer
log2Integer Prime
b
smoothProb :: Integer -> Width -> Log2Probability
smoothProb :: Prime -> Int -> Log2Integer
smoothProb = (Log2Integer -> Log2Integer) -> Int -> Log2Integer
forall a.
Integral a =>
(Log2Integer -> Log2Integer) -> a -> Log2Integer
prob ((Log2Integer -> Log2Integer) -> Int -> Log2Integer)
-> (Prime -> Log2Integer -> Log2Integer)
-> Prime
-> Int
-> Log2Integer
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Prime -> Log2Integer -> Log2Integer
smoothUnder
where
prob :: (Log2Integer -> Log2Integer) -> a -> Log2Integer
prob Log2Integer -> Log2Integer
f a
w = Log2Integer
f_t' Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Log2Integer -> Log2Integer
log2 (Log2Integer
2.0 Log2Integer -> Log2Integer -> Log2Integer
forall a. Floating a => a -> a -> a
** (Log2Integer -> Log2Integer
f Log2Integer
t Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
f_t') Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
1.0) Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
t'
where
f_t' :: Log2Integer
f_t' = Log2Integer -> Log2Integer
f Log2Integer
t'
t' :: Log2Integer
t' = Log2Integer
t Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer
1.0
t :: Log2Integer
t = a -> Log2Integer
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
w
smoothProbTrials :: Integer -> Width -> Integer -> Log2Probability
smoothProbTrials :: Prime -> Int -> Prime -> Log2Integer
smoothProbTrials Prime
b Int
t Prime
c = Log2Integer -> Log2Integer
log2 (Log2Integer
1.0 Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
- Log2Integer -> Log2Integer
forall a. Floating a => a -> a
Prelude.exp (-Log2Integer
x))
where x :: Log2Integer
x = Log2Integer
2.0 Log2Integer -> Log2Integer -> Log2Integer
forall a. Floating a => a -> a -> a
** (Prime -> Int -> Log2Integer
smoothProb Prime
b Int
t Log2Integer -> Log2Integer -> Log2Integer
forall a. Num a => a -> a -> a
+ Prime -> Log2Integer
log2Integer Prime
c)
type Gfp = Integer
valid :: Prime -> Gfp -> Bool
valid :: Prime -> Prime -> Bool
valid Prime
p Prime
x = Prime
0 Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
<= Prime
x Bool -> Bool -> Bool
&& Prime
x Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
< Prime
p
fromInt :: Prime -> Int -> Gfp
fromInt :: Prime -> Int -> Prime
fromInt Prime
p = Prime -> Prime -> Prime
Factor.Prime.fromInteger Prime
p (Prime -> Prime) -> (Int -> Prime) -> Int -> Prime
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Int -> Prime
forall a. Integral a => a -> Prime
toInteger
fromInteger :: Prime -> Integer -> Gfp
fromInteger :: Prime -> Prime -> Prime
fromInteger Prime
p Prime
n = Prime
n Prime -> Prime -> Prime
forall a. Integral a => a -> a -> a
`mod` Prime
p
toSmallestInteger :: Prime -> Gfp -> Integer
toSmallestInteger :: Prime -> Prime -> Prime
toSmallestInteger Prime
p Prime
x = if Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
x Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
< Prime
x then Prime
x Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
p else Prime
x
negate :: Prime -> Gfp -> Gfp
negate :: Prime -> Prime -> Prime
negate Prime
_ Prime
0 = Prime
0
negate Prime
p Prime
x = Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
x
add :: Prime -> Gfp -> Gfp -> Gfp
add :: Prime -> Prime -> Prime -> Prime
add Prime
p Prime
x Prime
y = Prime -> Prime -> Prime
Factor.Prime.fromInteger Prime
p (Prime
x Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
y)
sum :: Prime -> [Gfp] -> Gfp
sum :: Prime -> [Prime] -> Prime
sum Prime
p = (Prime -> Prime -> Prime) -> Prime -> [Prime] -> Prime
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Prime -> Prime -> Prime -> Prime
add Prime
p) Prime
0
subtract :: Prime -> Gfp -> Gfp -> Gfp
subtract :: Prime -> Prime -> Prime -> Prime
subtract Prime
p Prime
x Prime
y = Prime -> Prime -> Prime -> Prime
add Prime
p Prime
x (Prime -> Prime -> Prime
Factor.Prime.negate Prime
p Prime
y)
multiply :: Prime -> Gfp -> Gfp -> Gfp
multiply :: Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
x Prime
y = Prime -> Prime -> Prime
Factor.Prime.fromInteger Prime
p (Prime
x Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
* Prime
y)
square :: Prime -> Gfp -> Gfp
square :: Prime -> Prime -> Prime
square Prime
p Prime
x = Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
x Prime
x
cube :: Prime -> Gfp -> Gfp
cube :: Prime -> Prime -> Prime
cube Prime
p Prime
x = Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
x (Prime -> Prime -> Prime
square Prime
p Prime
x)
product :: Prime -> [Gfp] -> Gfp
product :: Prime -> [Prime] -> Prime
product Prime
p = (Prime -> Prime -> Prime) -> Prime -> [Prime] -> Prime
forall (t :: * -> *) a b.
Foldable t =>
(a -> b -> b) -> b -> t a -> b
foldr (Prime -> Prime -> Prime -> Prime
multiply Prime
p) Prime
1
multiplyExp :: Prime -> Gfp -> Gfp -> Integer -> Gfp
multiplyExp :: Prime -> Prime -> Prime -> Prime -> Prime
multiplyExp Prime
_ Prime
0 Prime
_ Prime
_ = Prime
0
multiplyExp Prime
_ Prime
z Prime
_ Prime
0 = Prime
z
multiplyExp Prime
_ Prime
_ Prime
0 Prime
_ = Prime
0
multiplyExp Prime
p Prime
z0 Prime
x0 Prime
k0 = Prime -> Prime -> Prime -> Prime
forall t. Integral t => Prime -> Prime -> t -> Prime
go Prime
z0 Prime
x0 Prime
k0
where
go :: Prime -> Prime -> t -> Prime
go Prime
z Prime
_ t
0 = Prime
z
go Prime
z Prime
1 t
_ = Prime
z
go Prime
z Prime
x t
k = Prime -> Prime -> t -> Prime
go Prime
z' Prime
x' t
k'
where
z' :: Prime
z' = if t -> Bool
forall a. Integral a => a -> Bool
even t
k then Prime
z else Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
z Prime
x
x' :: Prime
x' = Prime -> Prime -> Prime
square Prime
p Prime
x
k' :: t
k' = t
k t -> t -> t
forall a. Integral a => a -> a -> a
`div` t
2
exp :: Prime -> Gfp -> Integer -> Gfp
exp :: Prime -> Prime -> Prime -> Prime
exp Prime
p = Prime -> Prime -> Prime -> Prime -> Prime
multiplyExp Prime
p Prime
1
exp2 :: Prime -> Gfp -> Integer -> Gfp
exp2 :: Prime -> Prime -> Prime -> Prime
exp2 Prime
_ Prime
x Prime
0 = Prime
x
exp2 Prime
p Prime
x Prime
k = Prime -> Prime -> Prime -> Prime
exp2 Prime
p (Prime -> Prime -> Prime
square Prime
p Prime
x) (Prime
k Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1)
invertF :: Prime -> Gfp -> Factor Integer Gfp
invertF :: Prime -> Prime -> Factor Prime Prime
invertF Prime
_ Prime
0 = [Char] -> Factor Prime Prime
forall a. HasCallStack => [Char] -> a
error [Char]
"cannot invert zero"
invertF Prime
_ Prime
1 = Prime -> Factor Prime Prime
forall a b. b -> Either a b
Right Prime
1
invertF Prime
p Prime
x = if Prime
g Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
1 then Prime -> Factor Prime Prime
forall a b. b -> Either a b
Right Prime
i else Prime -> Factor Prime Prime
forall a b. a -> Either a b
Left Prime
g
where
(Prime
g,(Prime
_,Prime
t)) = Prime -> Prime -> (Prime, (Prime, Prime))
egcd Prime
p Prime
x
i :: Prime
i = if Prime
t Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
< Prime
0 then Prime
t Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
p else Prime
t
invert :: Prime -> Gfp -> Gfp
invert :: Prime -> Prime -> Prime
invert Prime
p Prime
x = Factor Prime Prime -> Prime
forall f a. Show f => Factor f a -> a
runFactor (Factor Prime Prime -> Prime) -> Factor Prime Prime -> Prime
forall a b. (a -> b) -> a -> b
$ Prime -> Prime -> Factor Prime Prime
invertF Prime
p Prime
x
divideF :: Prime -> Gfp -> Gfp -> Factor Integer Gfp
divideF :: Prime -> Prime -> Prime -> Factor Prime Prime
divideF Prime
p Prime
x Prime
y = do
Prime
z <- Prime -> Prime -> Factor Prime Prime
invertF Prime
p Prime
y
Prime -> Factor Prime Prime
forall (m :: * -> *) a. Monad m => a -> m a
return (Prime -> Factor Prime Prime) -> Prime -> Factor Prime Prime
forall a b. (a -> b) -> a -> b
$ Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
x Prime
z
divide :: Prime -> Gfp -> Gfp -> Gfp
divide :: Prime -> Prime -> Prime -> Prime
divide Prime
p Prime
x Prime
y = Factor Prime Prime -> Prime
forall f a. Show f => Factor f a -> a
runFactor (Factor Prime Prime -> Prime) -> Factor Prime Prime -> Prime
forall a b. (a -> b) -> a -> b
$ Prime -> Prime -> Prime -> Factor Prime Prime
divideF Prime
p Prime
x Prime
y
elements :: Prime -> [Gfp]
elements :: Prime -> [Prime]
elements Prime
p = [Prime
0 .. (Prime
pPrime -> Prime -> Prime
forall a. Num a => a -> a -> a
-Prime
1)]
uniform :: RandomGen r => Prime -> r -> (Gfp,r)
uniform :: Prime -> r -> (Prime, r)
uniform Prime
p = (Prime, Prime) -> r -> (Prime, r)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
Random.randomR (Prime
0, Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1)
uniformNonZero :: RandomGen r => Prime -> r -> (Gfp,r)
uniformNonZero :: Prime -> r -> (Prime, r)
uniformNonZero Prime
p = (Prime, Prime) -> r -> (Prime, r)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
Random.randomR (Prime
1, Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1)
millerRabinTrials :: Int
millerRabinTrials :: Int
millerRabinTrials = Int
100
isPrime :: RandomGen r => Integer -> r -> (Bool,r)
isPrime :: Prime -> r -> (Bool, r)
isPrime Prime
n | Prime
n Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
<= Prime
3 = (,) (Prime
2 Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
<= Prime
n)
isPrime Prime
n | Prime -> Bool
forall a. Integral a => a -> Bool
even Prime
n = (,) Bool
False
isPrime Prime
n = Int -> r -> (Bool, r)
forall t t. (Eq t, Num t, RandomGen t) => t -> t -> (Bool, t)
trials Int
millerRabinTrials
where
trials :: t -> t -> (Bool, t)
trials t
0 t
r = (Bool
True,t
r)
trials t
t t
r = if Prime -> Bool
trial Prime
a then t -> t -> (Bool, t)
trials (t
t t -> t -> t
forall a. Num a => a -> a -> a
- t
1) t
r' else (Bool
False,t
r')
where (Prime
a,t
r') = (Prime, Prime) -> t -> (Prime, t)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
Random.randomR (Prime
2, Prime
n Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1) t
r
trial :: Prime -> Bool
trial Prime
a = Bool -> Bool
not (Bool -> Bool) -> Bool -> Bool
forall a b. (a -> b) -> a -> b
$ Prime -> Prime -> Bool
forall t. (Eq t, Num t) => Prime -> t -> Bool
witness (Prime -> Prime -> Prime -> Prime
Factor.Prime.exp Prime
n Prime
a Prime
m) Prime
k
witness :: Prime -> t -> Bool
witness Prime
x t
0 = Prime
x Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
/= Prime
1
witness Prime
x t
i = if Prime
x2 Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
1 then Prime
x Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
/= Prime
1 Bool -> Bool -> Bool
&& Prime
x Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
/= Prime
n1 else Prime -> t -> Bool
witness Prime
x2 (t
i t -> t -> t
forall a. Num a => a -> a -> a
- t
1)
where x2 :: Prime
x2 = Prime -> Prime -> Prime
square Prime
n Prime
x
(Prime
k,Prime
m) = Prime -> Prime -> (Prime, Prime)
divPower Prime
2 Prime
n1
n1 :: Prime
n1 = Prime
n Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1
nextPrime :: RandomGen r => Integer -> r -> (Prime,r)
nextPrime :: Prime -> r -> (Prime, r)
nextPrime Prime
n r
r = if Bool
b then (Prime
n,r
r') else Prime -> r -> (Prime, r)
forall r. RandomGen r => Prime -> r -> (Prime, r)
nextPrime (Prime
n Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
1) r
r'
where (Bool
b,r
r') = Prime -> r -> (Bool, r)
forall r. RandomGen r => Prime -> r -> (Bool, r)
isPrime Prime
n r
r
previousPrime :: RandomGen r => Integer -> r -> (Prime,r)
previousPrime :: Prime -> r -> (Prime, r)
previousPrime Prime
n r
_ | Prime
n Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
< Prime
2 = [Char] -> (Prime, r)
forall a. HasCallStack => [Char] -> a
error [Char]
"no previous prime"
previousPrime Prime
n r
r = if Bool
b then (Prime
n,r
r') else Prime -> r -> (Prime, r)
forall r. RandomGen r => Prime -> r -> (Prime, r)
previousPrime (Prime
n Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1) r
r'
where (Bool
b,r
r') = Prime -> r -> (Bool, r)
forall r. RandomGen r => Prime -> r -> (Bool, r)
isPrime Prime
n r
r
uniformPrime :: RandomGen r => Int -> r -> (Prime,r)
uniformPrime :: Int -> r -> (Prime, r)
uniformPrime Int
w | Int
w Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
2 = [Char] -> r -> (Prime, r)
forall a. HasCallStack => [Char] -> a
error ([Char] -> r -> (Prime, r)) -> [Char] -> r -> (Prime, r)
forall a b. (a -> b) -> a -> b
$ [Char]
"no primes with width " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
w
uniformPrime Int
2 = Int -> r -> (Prime, r)
forall r. RandomGen r => Int -> r -> (Prime, r)
uniformInteger Int
2
uniformPrime Int
w = r -> (Prime, r)
forall t. RandomGen t => t -> (Prime, t)
pick
where
pick :: t -> (Prime, t)
pick t
r = if Bool
b then (Prime
n,t
r'') else t -> (Prime, t)
pick t
r''
where
(Prime
n,t
r') = Int -> t -> (Prime, t)
forall r. RandomGen r => Int -> r -> (Prime, r)
uniformOddInteger Int
w t
r
(Bool
b,t
r'') = Prime -> t -> (Bool, t)
forall r. RandomGen r => Prime -> r -> (Bool, r)
isPrime Prime
n t
r'
uniformComposite :: RandomGen r => Int -> r -> (Integer,r)
uniformComposite :: Int -> r -> (Prime, r)
uniformComposite Int
w | Int
w Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
3 = [Char] -> r -> (Prime, r)
forall a. HasCallStack => [Char] -> a
error ([Char] -> r -> (Prime, r)) -> [Char] -> r -> (Prime, r)
forall a b. (a -> b) -> a -> b
$ [Char]
"no composites with width " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
w
uniformComposite Int
w = r -> (Prime, r)
forall t. RandomGen t => t -> (Prime, t)
pick
where
pick :: b -> (Prime, b)
pick b
r = if Bool
b then b -> (Prime, b)
pick b
r'' else (Prime
n,b
r'')
where
(Prime
n,b
r') = Int -> b -> (Prime, b)
forall r. RandomGen r => Int -> r -> (Prime, r)
uniformInteger Int
w b
r
(Bool
b,b
r'') = Prime -> b -> (Bool, b)
forall r. RandomGen r => Prime -> r -> (Bool, r)
isPrime Prime
n b
r'
isResidue :: Prime -> Integer -> Bool
isResidue :: Prime -> Prime -> Bool
isResidue Prime
p Prime
n = Prime -> Prime -> Residue
jacobiSymbol Prime
n Prime
p Residue -> Residue -> Bool
forall a. Eq a => a -> a -> Bool
== Residue
Residue
nonResidue :: Prime -> Integer -> Bool
nonResidue :: Prime -> Prime -> Bool
nonResidue Prime
p Prime
n = Prime -> Prime -> Residue
jacobiSymbol Prime
n Prime
p Residue -> Residue -> Bool
forall a. Eq a => a -> a -> Bool
== Residue
NonResidue
nextResidue :: Prime -> Integer -> Integer
nextResidue :: Prime -> Prime -> Prime
nextResidue Prime
p Prime
n = if Prime -> Prime -> Bool
isResidue Prime
p Prime
n then Prime
n else Prime -> Prime -> Prime
nextResidue Prime
p (Prime
n Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
1)
nextNonResidue :: Prime -> Integer -> Integer
nextNonResidue :: Prime -> Prime -> Prime
nextNonResidue Prime
2 Prime
_ = [Char] -> Prime
forall a. HasCallStack => [Char] -> a
error [Char]
"no non-residues modulo 2"
nextNonResidue Prime
p Prime
n = if Prime -> Prime -> Bool
nonResidue Prime
p Prime
n then Prime
n else Prime -> Prime -> Prime
nextNonResidue Prime
p (Prime
n Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
1)
sqrt :: Prime -> Gfp -> Gfp
sqrt :: Prime -> Prime -> Prime
sqrt Prime
2 = Prime -> Prime
forall a. a -> a
id
sqrt Prime
p =
Prime -> Prime
norm (Prime -> Prime) -> (Prime -> Prime) -> Prime -> Prime
forall b c a. (b -> c) -> (a -> b) -> a -> c
.
(if Prime
r Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
1 then Prime -> Prime
sqrt3Mod4
else if Prime
r Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
2 then Prime -> Prime
sqrt5Mod8
else Prime -> Prime
tonelliShanks)
where
(Prime
r,Prime
s) = Prime -> Prime -> (Prime, Prime)
divPower Prime
2 (Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1)
norm :: Prime -> Prime
norm Prime
x = if Prime
2 Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
* Prime
x Prime -> Prime -> Bool
forall a. Ord a => a -> a -> Bool
< Prime
p then Prime
x else Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
x
sqrt3Mod4 :: Prime -> Prime
sqrt3Mod4 Prime
x = Prime -> Prime -> Prime -> Prime
Factor.Prime.exp Prime
p Prime
x Prime
sqrt3Mod4Exp
sqrt3Mod4Exp :: Prime
sqrt3Mod4Exp = (Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
1) Prime -> Prime -> Prime
forall a. Integral a => a -> a -> a
`div` Prime
4
sqrt5Mod8 :: Prime -> Prime
sqrt5Mod8 Prime
x = Prime -> [Prime] -> Prime
Factor.Prime.product Prime
p [Prime
x, Prime
y, Prime
z Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
1]
where
x2 :: Prime
x2 = Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
2 Prime
x
y :: Prime
y = Prime -> Prime -> Prime -> Prime
Factor.Prime.exp Prime
p Prime
x2 Prime
sqrt5Mod8Exp
z :: Prime
z = Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
x2 (Prime -> Prime -> Prime
square Prime
p Prime
y)
sqrt5Mod8Exp :: Prime
sqrt5Mod8Exp = (Prime
p Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- Prime
5) Prime -> Prime -> Prime
forall a. Integral a => a -> a -> a
`div` Prime
8
tonelliShanks :: Prime -> Prime
tonelliShanks Prime
x =
if Prime -> Prime -> Bool
isResidue Prime
p Prime
x then Prime -> Prime -> Prime -> Prime -> Prime
tonelliShanksLoop Prime
tonelliShanksInit Prime
d Prime
t Prime
r else Prime
0
where
d :: Prime
d = Prime -> Prime -> Prime -> Prime
Factor.Prime.exp Prime
p Prime
x ((Prime
s Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
1) Prime -> Prime -> Prime
forall a. Integral a => a -> a -> a
`div` Prime
2)
t :: Prime
t = Prime -> Prime -> Prime -> Prime
Factor.Prime.exp Prime
p Prime
x Prime
s
tonelliShanksInit :: Prime
tonelliShanksInit = Prime -> Prime -> Prime -> Prime
Factor.Prime.exp Prime
p (Prime -> Prime -> Prime
nextNonResidue Prime
p Prime
2) Prime
s
tonelliShanksLoop :: Prime -> Prime -> Prime -> Prime -> Prime
tonelliShanksLoop Prime
c Prime
d Prime
t Prime
m =
if Prime
t Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
1 then Prime
d else Prime -> Prime -> Prime -> Prime -> Prime
tonelliShanksLoop Prime
b2 Prime
db Prime
tb2 Prime
i
where
i :: Prime
i = Prime -> Prime -> Prime
forall t. Num t => Prime -> t -> t
tonelliShanksMin Prime
t Prime
1
b :: Prime
b = Prime -> Prime -> Prime -> Prime
exp2 Prime
p Prime
c (Prime
m Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
- (Prime
i Prime -> Prime -> Prime
forall a. Num a => a -> a -> a
+ Prime
1))
b2 :: Prime
b2 = Prime -> Prime -> Prime
square Prime
p Prime
b
db :: Prime
db = Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
d Prime
b
tb2 :: Prime
tb2 = Prime -> Prime -> Prime -> Prime
multiply Prime
p Prime
t Prime
b2
tonelliShanksMin :: Prime -> t -> t
tonelliShanksMin Prime
t t
i = if Prime
t2 Prime -> Prime -> Bool
forall a. Eq a => a -> a -> Bool
== Prime
1 then t
i else Prime -> t -> t
tonelliShanksMin Prime
t2 (t
i t -> t -> t
forall a. Num a => a -> a -> a
+ t
1)
where t2 :: Prime
t2 = Prime -> Prime -> Prime
square Prime
p Prime
t