{- |
module: Factor.Prime
description: Prime integers
license: MIT

maintainer: Joe Leslie-Hurd <joe@gilith.com>
stability: provisional
portability: portable
-}

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

-------------------------------------------------------------------------------
-- The genuine sieve of Eratosthenes
--
-- https://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf
-------------------------------------------------------------------------------

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)

-------------------------------------------------------------------------------
-- Trial division
-------------------------------------------------------------------------------

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

-------------------------------------------------------------------------------
-- Number of primes
--
-- The number of primes at most n is pi(n), which converges to n / log
-- n by the Prime Number Theorem.
--
--   log2(pi(n))
--   ==
--   log2(n / log(n))
--   ==
--   log2(n) - log2(log(n))
--   ==
--   log2(n) - log2(log2(n) / log2(e))
--   ==
--   log2(log2(e)) + log2(n) - log2(log2(n))
--
-- The nth prime function p_n is the inverse to pi and can be
-- estimated by the asymptotic inverse function n * log(n).
--
--   log2(p_n)
--   ==
--   log2(n * log(n))
--   ==
--   log2(n) + log2(log(n))
--   ==
--   log2(n) + log2(log2(n) / log2(e))
--   ==
--   log2(n) + log2(log2(n)) - log2(log2(e))
--
-- https://math.stackexchange.com/questions/606955/estimate-of-nth-prime
-------------------------------------------------------------------------------

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

-------------------------------------------------------------------------------
-- B-smooth numbers
--
-- The probability that a uniform integer in the range [1,N] is
-- B-smooth is given there as u^(−u), where u = log N / log B, giving
-- this estimate for the base-2 log of B-smooth integers at most N:
--
--   log2(N * u^(-u))
--   ==
--   log2(N) - u * log2(u)
--
-- For B fixed and much smaller than N, a better estimate is this
-- lower bound (the volume of a B-dimensional shape):
--
--   log2((1 / pi(B)!) * Prod { log2(N) / log2(p) | p <= B })
--   ==
--   Sum { -log2(i) | 1 <= i <= pi(B) }
--   + pi(B) * log2(log2(N))
--   + Sum { -log2(log2(p)) | p <= B }
--
-- Turning this into a probability:
--
--   log2(Prob(random t-bit number is B-smooth))
--   ==
--   log2((2^(smooth_B(t)) - 2^(smooth_B(t-1))) / 2^(t-1))
--   ==
--   (smooth_B(t-1) + log2(2^(smooth_B(t) - smooth_B(t-1)) - 1)) - (t-1)
--
-- For C independent trials:
--
--   log2(Prob(at least one of C random t-bit numbers is B-smooth))
--   ==
--   log2(1 - (1 - 2 ^ (log2(Prob(random t-bit number is B-smooth)))) ^ C)
--   ==
--   log2(1 - (1 - x / C) ^ C)                          [for some x, see below]
--   ==                                    [in the limit as C goes to infinity]
--   log2(1 - e^(-x))
--
-- where
--
--   x / C == 2 ^ (log2(Prob(random t-bit number is B-smooth)))
--   iff
--   x == 2 ^ (log2(Prob(random t-bit number is B-smooth)) + log2(C))
--
-- https://en.wikipedia.org/wiki/Smooth_number#Distribution
-- https://math.stackexchange.com/questions/22949/probability-of-being-b-smooth
-- "Prime Numbers: A Computational Perspective" by Crandall and Pomerance
-------------------------------------------------------------------------------

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)

-------------------------------------------------------------------------------
-- The field GF(p) of arithmetic modulo a prime
-------------------------------------------------------------------------------

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)

-------------------------------------------------------------------------------
-- Testing for primality using the Miller-Rabin probabilistic primality test
-------------------------------------------------------------------------------

millerRabinTrials :: Int
millerRabinTrials :: Int
millerRabinTrials = Int
100  -- False positive every 2^200 tests (or so)

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'

-------------------------------------------------------------------------------
-- Square roots in GF(p) using the Tonelli-Shanks algorithm
-------------------------------------------------------------------------------

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

{-
map (smoothUnder 20 . fromIntegral) [1..100]
map (smoothUnder 40 . fromIntegral) [1..100]
map (smoothUnder 100 . fromIntegral) [1..100]
map (smoothProb 100) [1..100]
map (smoothProb 1000) [1..100]
map (\t -> smoothProbTrials 1000 t 1000) [1..100]
previousPrime 200504042599 <$> Random.getStdGen
-}