-- |
-- Module:      Math.NumberTheory.Primes.Factorisation.Montgomery
-- Copyright:   (c) 2011 Daniel Fischer
-- Licence:     MIT
-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
--
-- Factorisation of 'Integer's by the elliptic curve algorithm after Montgomery.
-- The algorithm is explained at
-- <http://programmingpraxis.com/2010/04/23/modern-elliptic-curve-factorization-part-1/>
-- and
-- <http://programmingpraxis.com/2010/04/27/modern-elliptic-curve-factorization-part-2/>
--

{-# 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
  ( -- *  Complete factorisation functions
    -- ** Functions with input checking
    factorise
  --   -- * Partial factorisation
  , smallFactors
  --   -- ** Single curve worker
  , 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' n@ produces the prime factorisation of @n@. @'factorise' 0@ is
--   an error and the factorisation of @1@ is empty. Uses a 'StdGen' produced in
--   an arbitrary manner from the bit-pattern of @n@.
--
-- __Warning:__ there are no guarantees of any particular
-- order of prime factors, do not expect them to be ascending.
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)

----------------------------------------------------------------------------------------------------
--                                    Factorisation wrappers                                      --
----------------------------------------------------------------------------------------------------

-- | A wrapper around 'curveFactorisation' providing a few default arguments.
--   The primality test is 'bailliePSW', the @prng@ function - naturally -
--   'randomR'. This function also requires small prime factors to have been
--   stripped before.
stdGenFactorisation :: Maybe Integer     -- ^ Lower bound for composite divisors
                    -> StdGen            -- ^ Standard PRNG
                    -> Maybe Int         -- ^ Estimated number of digits of smallest prime factor
                    -> Integer           -- ^ The number to factorise
                    -> [(Integer, Word)] -- ^ List of prime factors and exponents
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' is the driver for the factorisation. Its performance (and success)
--   can be influenced by passing appropriate arguments. If you know that @n@ has no prime divisors
--   below @b@, any divisor found less than @b*b@ must be prime, thus giving @Just (b*b)@ as the
--   first argument allows skipping the comparatively expensive primality test for those.
--   If @n@ is such that all prime divisors must have a specific easy to test for structure, a
--   custom primality test can improve the performance (normally, it will make very little
--   difference, since @n@ has not many divisors, and many curves have to be tried to find one).
--   More influence has the pseudo random generator (a function @prng@ with @6 <= fst (prng k s) <= k-2@
--   and an initial state for the PRNG) used to generate the curves to try. A lucky choice here can
--   make a huge difference. So, if the default takes too long, try another one; or you can improve your
--   chances for a quick result by running several instances in parallel.
--
--   'curveFactorisation' @n@ requires that small (< 65536) prime factors of @n@
--   have been stripped before. Otherwise it is likely to cycle forever.
--
--   'curveFactorisation' is unlikely to succeed if @n@ has more than one (really) large prime factor.
--
curveFactorisation
  :: forall g.
     Maybe Integer                  -- ^ Lower bound for composite divisors
  -> (Integer -> Bool)              -- ^ A primality test
  -> (Integer -> g -> (Integer, g)) -- ^ A PRNG
  -> g                              -- ^ Initial PRNG state
  -> Maybe Int                      -- ^ Estimated number of digits of the smallest prime factor
  -> Integer                        -- ^ The number to factorise
  -> [(Integer, Word)]              -- ^ List of prime factors and exponents
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
          -- All factors (both @pfs@ and @cfs@), are pairwise coprime. This is
          -- because 'repFact' returns either a single factor, or output of 'workFact'.
          -- In its turn, 'workFact' returns either a single factor,
          -- or concats 'repFact's over coprime integers. Induction completes the proof.
          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)]
                -- Since all @cs@ are coprime, we can factor each of
                -- them and just concat results, without summing up
                -- powers of the same primes in different elements.
                ([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

-- | @'largePFPower' bd n@ produces the pair @(b,k)@ with the largest
--   exponent @k@ such that @n == b^k@, where @bd > 1@ (it is expected
--   that @bd@ is much larger, at least @1000@ or so), @n > bd^2@ and @n@
--   has no prime factors @p <= bd@, skipping the trial division phase
--   of @'highestPower'@ when that is a priori known to be superfluous.
--   It is only present to avoid duplication of work in factorisation
--   and primality testing, it is not expected to be generally useful.
--   The assumptions are not checked, if they are not satisfied, wrong
--   results and wasted work may be the consequence.
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)

----------------------------------------------------------------------------------------------------
--                                         The workhorse                                          --
----------------------------------------------------------------------------------------------------

-- | @'montgomeryFactorisation' n b1 b2 s@ tries to find a factor of @n@ using the
--   curve and point determined by the seed @s@ (@6 <= s < n-1@), multiplying the
--   point by the least common multiple of all numbers @<= b1@ and all primes
--   between @b1@ and @b2@. The idea is that there's a good chance that the order
--   of the point in the curve over one prime factor divides the multiplier, but the
--   order over another factor doesn't, if @b1@ and @b2@ are appropriately chosen.
--   If they are too small, none of the orders will probably divide the multiplier,
--   if they are too large, all probably will, so they should be chosen to fit
--   the expected size of the smallest factor.
--
--   It is assumed that @n@ has no small prime factors.
--
--   The result is maybe a nontrivial divisor of @n@.
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
    -- Small step: for each prime p <= b1
    -- multiply point 'p0' by the highest power p^k <= b1.
    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
      -- If small step did not succeed, perform a big step.
      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

-- | The implementation follows the algorithm at p. 6-7
-- of <http://www.hyperelliptic.org/tanja/SHARCS/talks06/Gaj.pdf Implementing the Elliptic Curve Method of Factoring in Reconfigurable Hardware>
-- by K. Gaj, S. Kwon et al.
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 ]

-- | Same as map (id *** flip multiply p) [from, thn .. to],
-- but calculated in more efficient way.
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)

-- primes, compactly stored as a bit sieve
primeStore :: [PrimeSieve]
primeStore :: [PrimeSieve]
primeStore = Integer -> [PrimeSieve]
psieveFrom Integer
7

-- generate list of primes from arrays
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' n@ finds all prime divisors of @n > 1@ up to 2^16 by trial division and returns the
--   list of these together with their multiplicities, and a possible remaining factor which may be composite.
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##) -- 65536 * 65536 - 1
        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)

-- | For a given estimated decimal length of the smallest prime factor
-- ("tier") return parameters B1, B2 and the number of curves to try
-- before next "tier".
-- Roughly based on http://www.mersennewiki.org/index.php/Elliptic_Curve_Method#Choosing_the_best_parameters_for_ECM
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)