{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE MagicHash #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE UnboxedTuples #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
module Math.NumberTheory.Primes.Factorisation.Montgomery
(
factorise
, smallFactors
, montgomeryFactorisation
, findParms
) where
import Control.Arrow
import Control.Monad.Trans.State.Lazy
import Data.Array.Base (bounds, unsafeAt)
import Data.Bits
import Data.IntMap (IntMap)
import qualified Data.IntMap as IM
import Data.List (foldl')
import Data.Maybe
import Data.Mod
import Data.Proxy
#if __GLASGOW_HASKELL__ < 803
import Data.Semigroup
#endif
import Data.Traversable
import GHC.Exts
import GHC.Integer.GMP.Internals hiding (integerToInt, wordToInteger)
import GHC.Natural
import GHC.TypeNats (KnownNat, SomeNat(..), natVal, someNatVal)
import System.Random
import Math.NumberTheory.Curves.Montgomery
import Math.NumberTheory.Euclidean.Coprimes (splitIntoCoprimes, unCoprimes)
import Math.NumberTheory.Logarithms (integerLogBase')
import Math.NumberTheory.Roots
import Math.NumberTheory.Primes.Sieve.Eratosthenes (PrimeSieve(..), psieveFrom)
import Math.NumberTheory.Primes.Sieve.Indexing (toPrim)
import Math.NumberTheory.Primes.Small
import Math.NumberTheory.Primes.Testing.Probabilistic
import Math.NumberTheory.Utils hiding (splitOff)
import Math.NumberTheory.Utils.FromIntegral
factorise :: Integral a => a -> [(a, Word)]
factorise :: a -> [(a, Word)]
factorise a
0 = [Char] -> [(a, Word)]
forall a. HasCallStack => [Char] -> a
error [Char]
"0 has no prime factorisation"
factorise a
n' = ((Natural, Word) -> (a, Word)) -> [(Natural, Word)] -> [(a, Word)]
forall a b. (a -> b) -> [a] -> [b]
map ((Natural -> a) -> (Natural, Word) -> (a, Word)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first Natural -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral) [(Natural, Word)]
sfs [(a, Word)] -> [(a, Word)] -> [(a, Word)]
forall a. Semigroup a => a -> a -> a
<> ((Integer, Word) -> (a, Word)) -> [(Integer, Word)] -> [(a, Word)]
forall a b. (a -> b) -> [a] -> [b]
map ((Integer -> a) -> (Integer, Word) -> (a, Word)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first Integer -> a
forall a. Num a => Integer -> a
fromInteger) [(Integer, Word)]
rest
where
n :: a
n = a -> a
forall a. Num a => a -> a
abs a
n'
([(Natural, Word)]
sfs, Maybe Natural
mb) = Natural -> ([(Natural, Word)], Maybe Natural)
smallFactors (a -> Natural
forall a b. (Integral a, Num b) => a -> b
fromIntegral' a
n)
sg :: StdGen
sg = Int -> StdGen
mkStdGen (a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n Int -> Int -> Int
forall a. Bits a => a -> a -> a
`xor` Int
0xdeadbeef)
rest :: [(Integer, Word)]
rest = case Maybe Natural
mb of
Maybe Natural
Nothing -> []
Just Natural
m -> Maybe Integer
-> StdGen -> Maybe Int -> Integer -> [(Integer, Word)]
stdGenFactorisation (Integer -> Maybe Integer
forall a. a -> Maybe a
Just (Integer -> Maybe Integer) -> Integer -> Maybe Integer
forall a b. (a -> b) -> a -> b
$ Integer
65536 Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Integer
65536) StdGen
sg Maybe Int
forall a. Maybe a
Nothing (Natural -> Integer
forall a. Integral a => a -> Integer
toInteger Natural
m)
stdGenFactorisation :: Maybe Integer
-> StdGen
-> Maybe Int
-> Integer
-> [(Integer, Word)]
stdGenFactorisation :: Maybe Integer
-> StdGen -> Maybe Int -> Integer -> [(Integer, Word)]
stdGenFactorisation Maybe Integer
primeBound =
Maybe Integer
-> (Integer -> Bool)
-> (Integer -> StdGen -> (Integer, StdGen))
-> StdGen
-> Maybe Int
-> Integer
-> [(Integer, Word)]
forall g.
Maybe Integer
-> (Integer -> Bool)
-> (Integer -> g -> (Integer, g))
-> g
-> Maybe Int
-> Integer
-> [(Integer, Word)]
curveFactorisation Maybe Integer
primeBound Integer -> Bool
bailliePSW (\Integer
m -> (Integer, Integer) -> StdGen -> (Integer, StdGen)
forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Integer
6, Integer
m Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Integer
2))
curveFactorisation
:: forall g.
Maybe Integer
-> (Integer -> Bool)
-> (Integer -> g -> (Integer, g))
-> g
-> Maybe Int
-> Integer
-> [(Integer, Word)]
curveFactorisation :: Maybe Integer
-> (Integer -> Bool)
-> (Integer -> g -> (Integer, g))
-> g
-> Maybe Int
-> Integer
-> [(Integer, Word)]
curveFactorisation Maybe Integer
primeBound Integer -> Bool
primeTest Integer -> g -> (Integer, g)
prng g
seed Maybe Int
mbdigs Integer
n
| Integer
n Integer -> Integer -> Bool
forall a. Eq a => a -> a -> Bool
== Integer
1 = []
| Integer -> Bool
ptest Integer
n = [(Integer
n, Word
1)]
| Bool
otherwise = State g [(Integer, Word)] -> g -> [(Integer, Word)]
forall s a. State s a -> s -> a
evalState (Integer -> Int -> State g [(Integer, Word)]
fact Integer
n Int
digits) g
seed
where
digits :: Int
digits :: Int
digits = Int -> Maybe Int -> Int
forall a. a -> Maybe a -> a
fromMaybe Int
8 Maybe Int
mbdigs
ptest :: Integer -> Bool
ptest :: Integer -> Bool
ptest = (Integer -> Bool)
-> (Integer -> Integer -> Bool) -> Maybe Integer -> Integer -> Bool
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Integer -> Bool
primeTest (\Integer
bd Integer
k -> Integer
k Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
<= Integer
bd Bool -> Bool -> Bool
|| Integer -> Bool
primeTest Integer
k) Maybe Integer
primeBound
rndR :: Integer -> State g Integer
rndR :: Integer -> State g Integer
rndR Integer
k = (g -> (Integer, g)) -> State g Integer
forall (m :: * -> *) s a. Monad m => (s -> (a, s)) -> StateT s m a
state (Integer -> g -> (Integer, g)
prng Integer
k)
perfPw :: Integer -> (Integer, Word)
perfPw :: Integer -> (Integer, Word)
perfPw = (Integer -> (Integer, Word))
-> (Integer -> Integer -> (Integer, Word))
-> Maybe Integer
-> Integer
-> (Integer, Word)
forall b a. b -> (a -> b) -> Maybe a -> b
maybe Integer -> (Integer, Word)
forall a. Integral a => a -> (a, Word)
highestPower (Integer -> Integer -> (Integer, Word)
largePFPower (Integer -> Integer -> (Integer, Word))
-> (Integer -> Integer) -> Integer -> Integer -> (Integer, Word)
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Integer -> Integer
forall a. Integral a => a -> a
integerSquareRoot) Maybe Integer
primeBound
fact :: Integer -> Int -> State g [(Integer, Word)]
fact :: Integer -> Int -> State g [(Integer, Word)]
fact Integer
1 Int
_ = [(Integer, Word)] -> State g [(Integer, Word)]
forall (m :: * -> *) a. Monad m => a -> m a
return [(Integer, Word)]
forall a. Monoid a => a
mempty
fact Integer
m Int
digs = do
let (Word
b1, Word
b2, Word
ct) = Int -> (Word, Word, Word)
findParms Int
digs
Factors [(Integer, Word)]
pfs [(Integer, Word)]
cfs <- Integer -> Word -> Word -> Word -> State g Factors
repFact Integer
m Word
b1 Word
b2 Word
ct
case [(Integer, Word)]
cfs of
[] -> [(Integer, Word)] -> State g [(Integer, Word)]
forall (m :: * -> *) a. Monad m => a -> m a
return [(Integer, Word)]
pfs
[(Integer, Word)]
_ -> do
[[(Integer, Word)]]
nfs <- [(Integer, Word)]
-> ((Integer, Word) -> State g [(Integer, Word)])
-> StateT g Identity [[(Integer, Word)]]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [(Integer, Word)]
cfs (((Integer, Word) -> State g [(Integer, Word)])
-> StateT g Identity [[(Integer, Word)]])
-> ((Integer, Word) -> State g [(Integer, Word)])
-> StateT g Identity [[(Integer, Word)]]
forall a b. (a -> b) -> a -> b
$ \(Integer
k, Word
j) ->
((Integer, Word) -> (Integer, Word))
-> [(Integer, Word)] -> [(Integer, Word)]
forall a b. (a -> b) -> [a] -> [b]
map ((Word -> Word) -> (Integer, Word) -> (Integer, Word)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second (Word -> Word -> Word
forall a. Num a => a -> a -> a
* Word
j)) ([(Integer, Word)] -> [(Integer, Word)])
-> State g [(Integer, Word)] -> State g [(Integer, Word)]
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Integer -> Int -> State g [(Integer, Word)]
fact Integer
k (if [(Integer, Word)] -> Bool
forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Integer, Word)]
pfs then Int
digs Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
5 else Int
digs)
[(Integer, Word)] -> State g [(Integer, Word)]
forall (m :: * -> *) a. Monad m => a -> m a
return ([(Integer, Word)] -> State g [(Integer, Word)])
-> [(Integer, Word)] -> State g [(Integer, Word)]
forall a b. (a -> b) -> a -> b
$ [[(Integer, Word)]] -> [(Integer, Word)]
forall a. Monoid a => [a] -> a
mconcat ([(Integer, Word)]
pfs [(Integer, Word)] -> [[(Integer, Word)]] -> [[(Integer, Word)]]
forall a. a -> [a] -> [a]
: [[(Integer, Word)]]
nfs)
repFact :: Integer -> Word -> Word -> Word -> State g Factors
repFact :: Integer -> Word -> Word -> Word -> State g Factors
repFact Integer
1 Word
_ Word
_ Word
_ = Factors -> State g Factors
forall (m :: * -> *) a. Monad m => a -> m a
return Factors
forall a. Monoid a => a
mempty
repFact Integer
m Word
b1 Word
b2 Word
count =
case Integer -> (Integer, Word)
perfPw Integer
m of
(Integer
_, Word
1) -> Integer -> Word -> Word -> Word -> State g Factors
workFact Integer
m Word
b1 Word
b2 Word
count
(Integer
b, Word
e)
| Integer -> Bool
ptest Integer
b -> Factors -> State g Factors
forall (m :: * -> *) a. Monad m => a -> m a
return (Factors -> State g Factors) -> Factors -> State g Factors
forall a b. (a -> b) -> a -> b
$ Integer -> Word -> Factors
singlePrimeFactor Integer
b Word
e
| Bool
otherwise -> (Word -> Word) -> Factors -> Factors
modifyPowers (Word -> Word -> Word
forall a. Num a => a -> a -> a
* Word
e) (Factors -> Factors) -> State g Factors -> State g Factors
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Integer -> Word -> Word -> Word -> State g Factors
workFact Integer
b Word
b1 Word
b2 Word
count
workFact :: Integer -> Word -> Word -> Word -> State g Factors
workFact :: Integer -> Word -> Word -> Word -> State g Factors
workFact Integer
1 Word
_ Word
_ Word
_ = Factors -> State g Factors
forall (m :: * -> *) a. Monad m => a -> m a
return Factors
forall a. Monoid a => a
mempty
workFact Integer
m Word
_ Word
_ Word
0 = Factors -> State g Factors
forall (m :: * -> *) a. Monad m => a -> m a
return (Factors -> State g Factors) -> Factors -> State g Factors
forall a b. (a -> b) -> a -> b
$ Integer -> Word -> Factors
singleCompositeFactor Integer
m Word
1
workFact Integer
m Word
b1 Word
b2 Word
count = do
Integer
s <- Integer -> State g Integer
rndR Integer
m
case Natural -> SomeNat
someNatVal (Integer -> Natural
forall a. Num a => Integer -> a
fromInteger Integer
m) of
SomeNat (Proxy n
_ :: Proxy t) -> case Word -> Word -> Mod n -> Maybe Integer
forall (n :: Nat).
KnownNat n =>
Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation Word
b1 Word
b2 (Integer -> Mod n
forall a. Num a => Integer -> a
fromInteger Integer
s :: Mod t) of
Maybe Integer
Nothing -> Integer -> Word -> Word -> Word -> State g Factors
workFact Integer
m Word
b1 Word
b2 (Word
count Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
1)
Just Integer
d -> do
let cs :: [(Integer, Word)]
cs = Coprimes Integer Word -> [(Integer, Word)]
forall a b. Coprimes a b -> [(a, b)]
unCoprimes (Coprimes Integer Word -> [(Integer, Word)])
-> Coprimes Integer Word -> [(Integer, Word)]
forall a b. (a -> b) -> a -> b
$ [(Integer, Word)] -> Coprimes Integer Word
forall a b.
(Eq a, GcdDomain a, Eq b, Num b) =>
[(a, b)] -> Coprimes a b
splitIntoCoprimes [(Integer
d, Word
1), (Integer
m Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`quot` Integer
d, Word
1)]
([Factors] -> Factors)
-> StateT g Identity [Factors] -> State g Factors
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap [Factors] -> Factors
forall a. Monoid a => [a] -> a
mconcat (StateT g Identity [Factors] -> State g Factors)
-> StateT g Identity [Factors] -> State g Factors
forall a b. (a -> b) -> a -> b
$ [(Integer, Word)]
-> ((Integer, Word) -> State g Factors)
-> StateT g Identity [Factors]
forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [(Integer, Word)]
cs (((Integer, Word) -> State g Factors)
-> StateT g Identity [Factors])
-> ((Integer, Word) -> State g Factors)
-> StateT g Identity [Factors]
forall a b. (a -> b) -> a -> b
$
\(Integer
x, Word
xm) -> if Integer -> Bool
ptest Integer
x
then Factors -> State g Factors
forall (f :: * -> *) a. Applicative f => a -> f a
pure (Factors -> State g Factors) -> Factors -> State g Factors
forall a b. (a -> b) -> a -> b
$ Integer -> Word -> Factors
singlePrimeFactor Integer
x Word
xm
else Integer -> Word -> Word -> Word -> State g Factors
repFact Integer
x Word
b1 Word
b2 (Word
count Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
1)
data Factors = Factors
{ Factors -> [(Integer, Word)]
_primeFactors :: [(Integer, Word)]
, Factors -> [(Integer, Word)]
_compositeFactors :: [(Integer, Word)]
}
singlePrimeFactor :: Integer -> Word -> Factors
singlePrimeFactor :: Integer -> Word -> Factors
singlePrimeFactor Integer
a Word
b = [(Integer, Word)] -> [(Integer, Word)] -> Factors
Factors [(Integer
a, Word
b)] []
singleCompositeFactor :: Integer -> Word -> Factors
singleCompositeFactor :: Integer -> Word -> Factors
singleCompositeFactor Integer
a Word
b = [(Integer, Word)] -> [(Integer, Word)] -> Factors
Factors [] [(Integer
a, Word
b)]
instance Semigroup Factors where
Factors [(Integer, Word)]
pfs1 [(Integer, Word)]
cfs1 <> :: Factors -> Factors -> Factors
<> Factors [(Integer, Word)]
pfs2 [(Integer, Word)]
cfs2
= [(Integer, Word)] -> [(Integer, Word)] -> Factors
Factors ([(Integer, Word)]
pfs1 [(Integer, Word)] -> [(Integer, Word)] -> [(Integer, Word)]
forall a. Semigroup a => a -> a -> a
<> [(Integer, Word)]
pfs2) ([(Integer, Word)]
cfs1 [(Integer, Word)] -> [(Integer, Word)] -> [(Integer, Word)]
forall a. Semigroup a => a -> a -> a
<> [(Integer, Word)]
cfs2)
instance Monoid Factors where
mempty :: Factors
mempty = [(Integer, Word)] -> [(Integer, Word)] -> Factors
Factors [] []
mappend :: Factors -> Factors -> Factors
mappend = Factors -> Factors -> Factors
forall a. Semigroup a => a -> a -> a
(<>)
modifyPowers :: (Word -> Word) -> Factors -> Factors
modifyPowers :: (Word -> Word) -> Factors -> Factors
modifyPowers Word -> Word
f (Factors [(Integer, Word)]
pfs [(Integer, Word)]
cfs)
= [(Integer, Word)] -> [(Integer, Word)] -> Factors
Factors (((Integer, Word) -> (Integer, Word))
-> [(Integer, Word)] -> [(Integer, Word)]
forall a b. (a -> b) -> [a] -> [b]
map ((Word -> Word) -> (Integer, Word) -> (Integer, Word)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second Word -> Word
f) [(Integer, Word)]
pfs) (((Integer, Word) -> (Integer, Word))
-> [(Integer, Word)] -> [(Integer, Word)]
forall a b. (a -> b) -> [a] -> [b]
map ((Word -> Word) -> (Integer, Word) -> (Integer, Word)
forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second Word -> Word
f) [(Integer, Word)]
cfs)
largePFPower :: Integer -> Integer -> (Integer, Word)
largePFPower :: Integer -> Integer -> (Integer, Word)
largePFPower Integer
bd Integer
n = Word -> Integer -> (Integer, Word)
rawPower Word
ln Integer
n
where
ln :: Word
ln = Int -> Word
intToWord (Integer -> Integer -> Int
integerLogBase' (Integer
bdInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
+Integer
1) Integer
n)
rawPower :: Word -> Integer -> (Integer, Word)
rawPower :: Word -> Integer -> (Integer, Word)
rawPower Word
mx Integer
n = case Integer -> Integer -> Maybe Integer
forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot Integer
4 Integer
n of
Just Integer
r -> case Word -> Integer -> (Integer, Word)
rawPower (Word
mx Word -> Word -> Word
forall a. Integral a => a -> a -> a
`quot` Word
4) Integer
r of
(Integer
m,Word
e) -> (Integer
m, Word
4Word -> Word -> Word
forall a. Num a => a -> a -> a
*Word
e)
Maybe Integer
Nothing -> case Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactSquareRoot Integer
n of
Just Integer
r -> case Word -> Integer -> (Integer, Word)
rawOddPower (Word
mx Word -> Word -> Word
forall a. Integral a => a -> a -> a
`quot` Word
2) Integer
r of
(Integer
m,Word
e) -> (Integer
m, Word
2Word -> Word -> Word
forall a. Num a => a -> a -> a
*Word
e)
Maybe Integer
Nothing -> Word -> Integer -> (Integer, Word)
rawOddPower Word
mx Integer
n
rawOddPower :: Word -> Integer -> (Integer, Word)
rawOddPower :: Word -> Integer -> (Integer, Word)
rawOddPower Word
mx Integer
n
| Word
mx Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
3 = (Integer
n,Word
1)
rawOddPower Word
mx Integer
n = case Integer -> Maybe Integer
forall a. Integral a => a -> Maybe a
exactCubeRoot Integer
n of
Just Integer
r -> case Word -> Integer -> (Integer, Word)
rawOddPower (Word
mx Word -> Word -> Word
forall a. Integral a => a -> a -> a
`quot` Word
3) Integer
r of
(Integer
m,Word
e) -> (Integer
m, Word
3Word -> Word -> Word
forall a. Num a => a -> a -> a
*Word
e)
Maybe Integer
Nothing -> Word -> Integer -> (Integer, Word)
badPower Word
mx Integer
n
badPower :: Word -> Integer -> (Integer, Word)
badPower :: Word -> Integer -> (Integer, Word)
badPower Word
mx Integer
n
| Word
mx Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
< Word
5 = (Integer
n,Word
1)
| Bool
otherwise = Word -> Word -> Integer -> [Word] -> (Integer, Word)
forall t a.
(Integral t, Integral a) =>
a -> a -> t -> [a] -> (t, a)
go Word
1 Word
mx Integer
n ((Word -> Bool) -> [Word] -> [Word]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
mx) ([Word] -> [Word]) -> [Word] -> [Word]
forall a b. (a -> b) -> a -> b
$ (Word -> Word -> Word) -> Word -> [Word] -> [Word]
forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl Word -> Word -> Word
forall a. Num a => a -> a -> a
(+) Word
5 ([Word] -> [Word]) -> [Word] -> [Word]
forall a b. (a -> b) -> a -> b
$ [Word] -> [Word]
forall a. [a] -> [a]
cycle [Word
2,Word
4])
where
go :: a -> a -> t -> [a] -> (t, a)
go !a
e a
b t
m (a
k:[a]
ks)
| a
b a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
k = (t
m,a
e)
| Bool
otherwise = case a -> t -> Maybe t
forall a b. (Integral a, Integral b) => b -> a -> Maybe a
exactRoot a
k t
m of
Just t
r -> a -> a -> t -> [a] -> (t, a)
go (a
ea -> a -> a
forall a. Num a => a -> a -> a
*a
k) (a
b a -> a -> a
forall a. Integral a => a -> a -> a
`quot` a
k) t
r (a
ka -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
ks)
Maybe t
Nothing -> a -> a -> t -> [a] -> (t, a)
go a
e a
b t
m [a]
ks
go a
e a
_ t
m [] = (t
m,a
e)
montgomeryFactorisation :: KnownNat n => Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation :: Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation Word
b1 Word
b2 Mod n
s = case Integer -> Integer -> Maybe SomePoint
newPoint (Natural -> Integer
forall a. Integral a => a -> Integer
toInteger (Mod n -> Natural
forall (m :: Nat). Mod m -> Natural
unMod Mod n
s)) Integer
n of
Maybe SomePoint
Nothing -> Maybe Integer
forall a. Maybe a
Nothing
Just (SomePoint Point a24 n
p0) -> do
let q :: Point a24 n
q = (Point a24 n -> Word -> Point a24 n)
-> Point a24 n -> [Word] -> Point a24 n
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl ((Word -> Point a24 n -> Point a24 n)
-> Point a24 n -> Word -> Point a24 n
forall a b c. (a -> b -> c) -> b -> a -> c
flip Word -> Point a24 n -> Point a24 n
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply) Point a24 n
p0 [Word]
smallPowers
z :: Integer
z = Point a24 n -> Integer
forall (a24 :: Nat) (n :: Nat). Point a24 n -> Integer
pointZ Point a24 n
q
case Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n Integer
z of
Integer
1 -> case Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
gcd Integer
n (Point a24 n -> Word -> Word -> Integer
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Point a24 n -> Word -> Word -> Integer
bigStep Point a24 n
q Word
b1 Word
b2) of
Integer
1 -> Maybe Integer
forall a. Maybe a
Nothing
Integer
g -> Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
g
Integer
g -> Integer -> Maybe Integer
forall a. a -> Maybe a
Just Integer
g
where
n :: Integer
n = Natural -> Integer
forall a. Integral a => a -> Integer
toInteger (Mod n -> Natural
forall (n :: Nat) (proxy :: Nat -> *).
KnownNat n =>
proxy n -> Natural
natVal Mod n
s)
smallPowers :: [Word]
smallPowers
= (Word -> Word) -> [Word] -> [Word]
forall a b. (a -> b) -> [a] -> [b]
map Word -> Word
findPower
([Word] -> [Word]) -> [Word] -> [Word]
forall a b. (a -> b) -> a -> b
$ (Word -> Bool) -> [Word] -> [Word]
forall a. (a -> Bool) -> [a] -> [a]
takeWhile (Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
b1) (Word
2 Word -> [Word] -> [Word]
forall a. a -> [a] -> [a]
: Word
3 Word -> [Word] -> [Word]
forall a. a -> [a] -> [a]
: Word
5 Word -> [Word] -> [Word]
forall a. a -> [a] -> [a]
: [PrimeSieve] -> [Word]
list [PrimeSieve]
primeStore)
findPower :: Word -> Word
findPower Word
p = Word -> Word
go Word
p
where
go :: Word -> Word
go Word
acc
| Word
acc Word -> Word -> Bool
forall a. Ord a => a -> a -> Bool
<= Word
b1 Word -> Word -> Word
forall a. Integral a => a -> a -> a
`quot` Word
p = Word -> Word
go (Word
acc Word -> Word -> Word
forall a. Num a => a -> a -> a
* Word
p)
| Bool
otherwise = Word
acc
bigStep :: (KnownNat a24, KnownNat n) => Point a24 n -> Word -> Word -> Integer
bigStep :: Point a24 n -> Word -> Word -> Integer
bigStep Point a24 n
q Word
b1 Word
b2 = Integer
rs
where
n :: Integer
n = Point a24 n -> Integer
forall (a24 :: Nat) (n :: Nat).
KnownNat n =>
Point a24 n -> Integer
pointN Point a24 n
q
b0 :: Word
b0 = Word
b1 Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
b1 Word -> Word -> Word
forall a. Integral a => a -> a -> a
`rem` Word
wheel
qks :: [(Integer, Point a24 n)]
qks = [Integer] -> [Point a24 n] -> [(Integer, Point a24 n)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] ([Point a24 n] -> [(Integer, Point a24 n)])
-> [Point a24 n] -> [(Integer, Point a24 n)]
forall a b. (a -> b) -> a -> b
$ (Word -> Point a24 n) -> [Word] -> [Point a24 n]
forall a b. (a -> b) -> [a] -> [b]
map (Word -> Point a24 n -> Point a24 n
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
`multiply` Point a24 n
q) [Word]
wheelCoprimes
qs :: [(Word, Point a24 n)]
qs = Point a24 n -> Word -> Word -> Word -> [(Word, Point a24 n)]
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Point a24 n -> Word -> Word -> Word -> [(Word, Point a24 n)]
enumAndMultiplyFromThenTo Point a24 n
q Word
b0 (Word
b0 Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Word
wheel) Word
b2
rs :: Integer
rs = (Integer -> (Word, Point a24 n) -> Integer)
-> Integer -> [(Word, Point a24 n)] -> Integer
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (\Integer
ts (Word
_cHi, Point a24 n
p) -> (Integer -> (Integer, Point a24 n) -> Integer)
-> Integer -> [(Integer, Point a24 n)] -> Integer
forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (\Integer
us (Integer
_cLo, Point a24 n
pq) ->
Integer
us Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* (Point a24 n -> Integer
forall (a24 :: Nat) (n :: Nat). Point a24 n -> Integer
pointZ Point a24 n
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Point a24 n -> Integer
forall (a24 :: Nat) (n :: Nat). Point a24 n -> Integer
pointX Point a24 n
pq Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- Point a24 n -> Integer
forall (a24 :: Nat) (n :: Nat). Point a24 n -> Integer
pointX Point a24 n
p Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Point a24 n -> Integer
forall (a24 :: Nat) (n :: Nat). Point a24 n -> Integer
pointZ Point a24 n
pq) Integer -> Integer -> Integer
forall a. Integral a => a -> a -> a
`rem` Integer
n
) Integer
ts [(Integer, Point a24 n)]
qks) Integer
1 [(Word, Point a24 n)]
qs
wheel :: Word
wheel :: Word
wheel = Word
210
wheelCoprimes :: [Word]
wheelCoprimes :: [Word]
wheelCoprimes = [ Word
k | Word
k <- [Word
1 .. Word
wheel Word -> Word -> Word
forall a. Integral a => a -> a -> a
`div` Word
2], Word
k Word -> Word -> Word
forall a. Integral a => a -> a -> a
`gcd` Word
wheel Word -> Word -> Bool
forall a. Eq a => a -> a -> Bool
== Word
1 ]
enumAndMultiplyFromThenTo
:: (KnownNat a24, KnownNat n)
=> Point a24 n
-> Word
-> Word
-> Word
-> [(Word, Point a24 n)]
enumAndMultiplyFromThenTo :: Point a24 n -> Word -> Word -> Word -> [(Word, Point a24 n)]
enumAndMultiplyFromThenTo Point a24 n
p Word
from Word
thn Word
to = [Word] -> [Point a24 n] -> [(Word, Point a24 n)]
forall a b. [a] -> [b] -> [(a, b)]
zip [Word
from, Word
thn .. Word
to] [Point a24 n]
progression
where
step :: Word
step = Word
thn Word -> Word -> Word
forall a. Num a => a -> a -> a
- Word
from
pFrom :: Point a24 n
pFrom = Word -> Point a24 n -> Point a24 n
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply Word
from Point a24 n
p
pThen :: Point a24 n
pThen = Word -> Point a24 n -> Point a24 n
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply Word
thn Point a24 n
p
pStep :: Point a24 n
pStep = Word -> Point a24 n -> Point a24 n
forall (a24 :: Nat) (n :: Nat).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply Word
step Point a24 n
p
progression :: [Point a24 n]
progression = Point a24 n
pFrom Point a24 n -> [Point a24 n] -> [Point a24 n]
forall a. a -> [a] -> [a]
: Point a24 n
pThen Point a24 n -> [Point a24 n] -> [Point a24 n]
forall a. a -> [a] -> [a]
: (Point a24 n -> Point a24 n -> Point a24 n)
-> [Point a24 n] -> [Point a24 n] -> [Point a24 n]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (Point a24 n -> Point a24 n -> Point a24 n -> Point a24 n
forall (n :: Nat) (a24 :: Nat).
KnownNat n =>
Point a24 n -> Point a24 n -> Point a24 n -> Point a24 n
`add` Point a24 n
pStep) [Point a24 n]
progression ([Point a24 n] -> [Point a24 n]
forall a. [a] -> [a]
tail [Point a24 n]
progression)
primeStore :: [PrimeSieve]
primeStore :: [PrimeSieve]
primeStore = Integer -> [PrimeSieve]
psieveFrom Integer
7
list :: [PrimeSieve] -> [Word]
list :: [PrimeSieve] -> [Word]
list [PrimeSieve]
sieves = [[Word]] -> [Word]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Word
off Word -> Word -> Word
forall a. Num a => a -> a -> a
+ Int -> Word
forall a. Num a => Int -> a
toPrim Int
i | Int
i <- [Int
0 .. Int
li], UArray Int Bool -> Int -> Bool
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> Int -> e
unsafeAt UArray Int Bool
bs Int
i]
| PS Integer
vO UArray Int Bool
bs <- [PrimeSieve]
sieves, let { (Int
_,Int
li) = UArray Int Bool -> (Int, Int)
forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
bounds UArray Int Bool
bs; off :: Word
off = Integer -> Word
forall a. Num a => Integer -> a
fromInteger Integer
vO; }]
smallFactors :: Natural -> ([(Natural, Word)], Maybe Natural)
smallFactors :: Natural -> ([(Natural, Word)], Maybe Natural)
smallFactors = \case
NatS# Word#
0## -> [Char] -> ([(Natural, Word)], Maybe Natural)
forall a. HasCallStack => [Char] -> a
error [Char]
"0 has no prime factorisation"
NatS# Word#
n# -> case Word# -> (# Word#, Word# #)
shiftToOddCount# Word#
n# of
(# Word#
0##, Word#
m# #) -> Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord Word#
m# Int
1
(# Word#
k#, Word#
m# #) -> (Natural
2, Word# -> Word
W# Word#
k#) (Natural, Word)
-> ([(Natural, Word)], Maybe Natural)
-> ([(Natural, Word)], Maybe Natural)
forall a b. a -> ([a], b) -> ([a], b)
<: Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord Word#
m# Int
1
NatJ# BigNat
n -> case BigNat -> (Word, BigNat)
shiftToOddCountBigNat BigNat
n of
(Word
0, BigNat
m) -> BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat BigNat
m Int
1
(Word
k, BigNat
m) -> (Natural
2, Word
k) (Natural, Word)
-> ([(Natural, Word)], Maybe Natural)
-> ([(Natural, Word)], Maybe Natural)
forall a b. a -> ([a], b) -> ([a], b)
<: BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat BigNat
m Int
1
where
a
x <: :: a -> ([a], b) -> ([a], b)
<: ~([a]
l,b
b) = (a
xa -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
l,b
b)
!(Ptr Addr#
smallPrimesAddr#) = Ptr Word16
smallPrimesPtr
goBigNat :: BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat :: BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat !BigNat
m i :: Int
i@(I# Int#
i#)
| Int# -> Bool
isTrue# (BigNat -> Int#
sizeofBigNat# BigNat
m Int# -> Int# -> Int#
==# Int#
1#)
= Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord (BigNat -> Word#
bigNatToWord BigNat
m) Int
i
| Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
smallPrimesLength
= ([], Natural -> Maybe Natural
forall a. a -> Maybe a
Just (BigNat -> Natural
NatJ# BigNat
m))
| Bool
otherwise
= let p# :: Word#
p# =
#if MIN_VERSION_base(4,16,0)
word16ToWord#
#endif
(Addr# -> Int# -> Word#
indexWord16OffAddr# Addr#
smallPrimesAddr# Int#
i#) in
case BigNat
m BigNat -> Word# -> (# BigNat, Word# #)
`quotRemBigNatWord` Word#
p# of
(# BigNat
mp, Word#
0## #) ->
let (# Word
k, BigNat
r #) = Word -> BigNat -> (# Word, BigNat #)
forall a. Num a => a -> BigNat -> (# a, BigNat #)
splitOff Word
1 BigNat
mp in
(Word# -> Natural
NatS# Word#
p#, Word
k) (Natural, Word)
-> ([(Natural, Word)], Maybe Natural)
-> ([(Natural, Word)], Maybe Natural)
forall a b. a -> ([a], b) -> ([a], b)
<: BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat BigNat
r (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
where
splitOff :: a -> BigNat -> (# a, BigNat #)
splitOff !a
k BigNat
x = case BigNat
x BigNat -> Word# -> (# BigNat, Word# #)
`quotRemBigNatWord` Word#
p# of
(# BigNat
xp, Word#
0## #) -> a -> BigNat -> (# a, BigNat #)
splitOff (a
k a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) BigNat
xp
(# BigNat, Word# #)
_ -> (# a
k, BigNat
x #)
(# BigNat, Word# #)
_ -> BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat BigNat
m (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
goWord :: Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord :: Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord Word#
1## !Int
_ = ([], Maybe Natural
forall a. Maybe a
Nothing)
goWord Word#
m# !Int
i
| Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
smallPrimesLength
= if Int# -> Bool
isTrue# (Word#
m# Word# -> Word# -> Int#
`leWord#` Word#
4294967295##)
then ([(Word# -> Natural
NatS# Word#
m#, Word
1)], Maybe Natural
forall a. Maybe a
Nothing)
else ([], Natural -> Maybe Natural
forall a. a -> Maybe a
Just (Word# -> Natural
NatS# Word#
m#))
goWord Word#
m# i :: Int
i@(I# Int#
i#)
= let p# :: Word#
p# =
#if MIN_VERSION_base(4,16,0)
word16ToWord#
#endif
(Addr# -> Int# -> Word#
indexWord16OffAddr# Addr#
smallPrimesAddr# Int#
i#) in
if Int# -> Bool
isTrue# (Word#
m# Word# -> Word# -> Int#
`ltWord#` (Word#
p# Word# -> Word# -> Word#
`timesWord#` Word#
p#))
then ([(Word# -> Natural
NatS# Word#
m#, Word
1)], Maybe Natural
forall a. Maybe a
Nothing)
else case Word#
m# Word# -> Word# -> (# Word#, Word# #)
`quotRemWord#` Word#
p# of
(# Word#
mp#, Word#
0## #) ->
let !(# Word#
k#, Word#
r# #) = Word# -> Word# -> (# Word#, Word# #)
splitOff Word#
1## Word#
mp# in
(Word# -> Natural
NatS# Word#
p#, Word# -> Word
W# Word#
k#) (Natural, Word)
-> ([(Natural, Word)], Maybe Natural)
-> ([(Natural, Word)], Maybe Natural)
forall a b. a -> ([a], b) -> ([a], b)
<: Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord Word#
r# (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
where
splitOff :: Word# -> Word# -> (# Word#, Word# #)
splitOff Word#
k# Word#
x# = case Word#
x# Word# -> Word# -> (# Word#, Word# #)
`quotRemWord#` Word#
p# of
(# Word#
xp#, Word#
0## #) -> Word# -> Word# -> (# Word#, Word# #)
splitOff (Word#
k# Word# -> Word# -> Word#
`plusWord#` Word#
1##) Word#
xp#
(# Word#, Word# #)
_ -> (# Word#
k#, Word#
x# #)
(# Word#, Word# #)
_ -> Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord Word#
m# (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1)
testParms :: IntMap (Word, Word, Word)
testParms :: IntMap (Word, Word, Word)
testParms = [(Int, (Word, Word, Word))] -> IntMap (Word, Word, Word)
forall a. [(Int, a)] -> IntMap a
IM.fromList
[ (Int
12, ( Word
400, Word
40000, Word
10))
, (Int
15, ( Word
2000, Word
200000, Word
25))
, (Int
20, ( Word
11000, Word
1100000, Word
90))
, (Int
25, ( Word
50000, Word
5000000, Word
300))
, (Int
30, ( Word
250000, Word
25000000, Word
700))
, (Int
35, ( Word
1000000, Word
100000000, Word
1800))
, (Int
40, ( Word
3000000, Word
300000000, Word
5100))
, (Int
45, ( Word
11000000, Word
1100000000, Word
10600))
, (Int
50, ( Word
43000000, Word
4300000000, Word
19300))
, (Int
55, ( Word
110000000, Word
11000000000, Word
49000))
, (Int
60, ( Word
260000000, Word
26000000000, Word
124000))
, (Int
65, ( Word
850000000, Word
85000000000, Word
210000))
, (Int
70, (Word
2900000000, Word
290000000000, Word
340000))
]
findParms :: Int -> (Word, Word, Word)
findParms :: Int -> (Word, Word, Word)
findParms Int
digs = (Word, Word, Word)
-> ((Int, (Word, Word, Word)) -> (Word, Word, Word))
-> Maybe (Int, (Word, Word, Word))
-> (Word, Word, Word)
forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Word
wheel, Word
1000, Word
7) (Int, (Word, Word, Word)) -> (Word, Word, Word)
forall a b. (a, b) -> b
snd (Int -> IntMap (Word, Word, Word) -> Maybe (Int, (Word, Word, Word))
forall a. Int -> IntMap a -> Maybe (Int, a)
IM.lookupLT Int
digs IntMap (Word, Word, Word)
testParms)