-- |
-- 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 :: forall a. Integral a => a -> [(a, Word)]
factorise a
0 = forall a. HasCallStack => [Char] -> a
error [Char]
"0 has no prime factorisation"
factorise a
n' = forall a b. (a -> b) -> [a] -> [b]
map (forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first forall a b. (Integral a, Num b) => a -> b
fromIntegral) [(Natural, Word)]
sfs forall a. Semigroup a => a -> a -> a
<> forall a b. (a -> b) -> [a] -> [b]
map (forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (b, d) (c, d)
first forall a. Num a => Integer -> a
fromInteger) [(Integer, Word)]
rest
  where
    n :: a
n = forall a. Num a => a -> a
abs a
n'
    ([(Natural, Word)]
sfs, Maybe Natural
mb) = Natural -> ([(Natural, Word)], Maybe Natural)
smallFactors (forall a b. (Integral a, Num b) => a -> b
fromIntegral' a
n)
    sg :: StdGen
sg = Int -> StdGen
mkStdGen (forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n 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 (forall a. a -> Maybe a
Just forall a b. (a -> b) -> a -> b
$ Integer
65536 forall a. Num a => a -> a -> a
* Integer
65536) StdGen
sg forall a. Maybe a
Nothing (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 =
  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 -> forall a g. (Random a, RandomGen g) => (a, a) -> g -> (a, g)
randomR (Integer
6, Integer
m 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 :: forall g.
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 forall a. Eq a => a -> a -> Bool
== Integer
1    = []
    | Integer -> Bool
ptest Integer
n   = [(Integer
n, Word
1)]
    | Bool
otherwise = 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 = forall a. a -> Maybe a -> a
fromMaybe Int
8 Maybe Int
mbdigs

        ptest :: Integer -> Bool
        ptest :: Integer -> Bool
ptest = forall b a. b -> (a -> b) -> Maybe a -> b
maybe Integer -> Bool
primeTest (\Integer
bd Integer
k -> Integer
k 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 = 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 = forall b a. b -> (a -> b) -> Maybe a -> b
maybe forall a. Integral a => a -> (a, Word)
highestPower (Integer -> Integer -> (Integer, Word)
largePFPower forall b c a. (b -> c) -> (a -> b) -> a -> c
. 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
_ = forall (m :: * -> *) a. Monad m => a -> m a
return 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
            [] -> forall (m :: * -> *) a. Monad m => a -> m a
return [(Integer, Word)]
pfs
            [(Integer, Word)]
_  -> do
              [[(Integer, Word)]]
nfs <- forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [(Integer, Word)]
cfs forall a b. (a -> b) -> a -> b
$ \(Integer
k, Word
j) ->
                  forall a b. (a -> b) -> [a] -> [b]
map (forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second (forall a. Num a => a -> a -> a
* Word
j)) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Integer -> Int -> State g [(Integer, Word)]
fact Integer
k (if forall (t :: * -> *) a. Foldable t => t a -> Bool
null [(Integer, Word)]
pfs then Int
digs forall a. Num a => a -> a -> a
+ Int
5 else Int
digs)
              forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ forall a. Monoid a => [a] -> a
mconcat ([(Integer, Word)]
pfs 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
_ = forall (m :: * -> *) a. Monad m => a -> m a
return 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   -> forall (m :: * -> *) a. Monad m => a -> m a
return forall a b. (a -> b) -> a -> b
$ Integer -> Word -> Factors
singlePrimeFactor Integer
b Word
e
              | Bool
otherwise -> (Word -> Word) -> Factors -> Factors
modifyPowers (forall a. Num a => a -> a -> a
* Word
e) 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
_ = forall (m :: * -> *) a. Monad m => a -> m a
return forall a. Monoid a => a
mempty
        workFact Integer
m Word
_ Word
_ Word
0 = forall (m :: * -> *) a. Monad m => a -> m a
return 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 (forall a. Num a => Integer -> a
fromInteger Integer
m) of
            SomeNat (Proxy n
_ :: Proxy t) -> case forall (n :: Natural).
KnownNat n =>
Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation Word
b1 Word
b2 (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 forall a. Num a => a -> a -> a
- Word
1)
              Just Integer
d  -> do
                let cs :: [(Integer, Word)]
cs = forall a b. Coprimes a b -> [(a, b)]
unCoprimes forall a b. (a -> b) -> a -> b
$ forall a b.
(Eq a, GcdDomain a, Eq b, Num b) =>
[(a, b)] -> Coprimes a b
splitIntoCoprimes [(Integer
d, Word
1), (Integer
m 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.
                forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap forall a. Monoid a => [a] -> a
mconcat forall a b. (a -> b) -> a -> b
$ forall (t :: * -> *) (m :: * -> *) a b.
(Traversable t, Monad m) =>
t a -> (a -> m b) -> m (t b)
forM [(Integer, Word)]
cs forall a b. (a -> b) -> a -> b
$
                  \(Integer
x, Word
xm) -> if Integer -> Bool
ptest Integer
x
                              then forall (f :: * -> *) a. Applicative f => a -> f a
pure forall a b. (a -> b) -> a -> b
$ Integer -> Word -> Factors
singlePrimeFactor Integer
x Word
xm
                              else forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
fmap ((Word -> Word) -> Factors -> Factors
modifyPowers (forall a. Num a => a -> a -> a
* Word
xm)) (Integer -> Word -> Word -> Word -> State g Factors
repFact Integer
x Word
b1 Word
b2 (Word
count forall a. Num a => a -> a -> a
- Word
1))

data Factors = Factors
  { Factors -> [(Integer, Word)]
_primeFactors     :: [(Integer, Word)]
  , Factors -> [(Integer, Word)]
_compositeFactors :: [(Integer, Word)]
  } deriving (Int -> Factors -> ShowS
[Factors] -> ShowS
Factors -> [Char]
forall a.
(Int -> a -> ShowS) -> (a -> [Char]) -> ([a] -> ShowS) -> Show a
showList :: [Factors] -> ShowS
$cshowList :: [Factors] -> ShowS
show :: Factors -> [Char]
$cshow :: Factors -> [Char]
showsPrec :: Int -> Factors -> ShowS
$cshowsPrec :: Int -> Factors -> ShowS
Show)

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 forall a. Semigroup a => a -> a -> a
<> [(Integer, Word)]
pfs2) ([(Integer, Word)]
cfs1 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 = 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 (forall a b. (a -> b) -> [a] -> [b]
map (forall (a :: * -> * -> *) b c d.
Arrow a =>
a b c -> a (d, b) (d, c)
second Word -> Word
f) [(Integer, Word)]
pfs) (forall a b. (a -> b) -> [a] -> [b]
map (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
bdforall 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 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 forall a. Integral a => a -> a -> a
`quot` Word
4) Integer
r of
                              (Integer
m,Word
e) -> (Integer
m, Word
4forall a. Num a => a -> a -> a
*Word
e)
                  Maybe Integer
Nothing -> case forall a. Integral a => a -> Maybe a
exactSquareRoot Integer
n of
                               Just Integer
r -> case Word -> Integer -> (Integer, Word)
rawOddPower (Word
mx forall a. Integral a => a -> a -> a
`quot` Word
2) Integer
r of
                                           (Integer
m,Word
e) -> (Integer
m, Word
2forall 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 forall a. Ord a => a -> a -> Bool
< Word
3       = (Integer
n,Word
1)
rawOddPower Word
mx Integer
n = case forall a. Integral a => a -> Maybe a
exactCubeRoot Integer
n of
                     Just Integer
r -> case Word -> Integer -> (Integer, Word)
rawOddPower (Word
mx forall a. Integral a => a -> a -> a
`quot` Word
3) Integer
r of
                                 (Integer
m,Word
e) -> (Integer
m, Word
3forall 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 forall a. Ord a => a -> a -> Bool
< Word
5      = (Integer
n,Word
1)
  | Bool
otherwise   = forall {t} {a}.
(Integral t, Integral a) =>
a -> a -> t -> [a] -> (t, a)
go Word
1 Word
mx Integer
n (forall a. (a -> Bool) -> [a] -> [a]
takeWhile (forall a. Ord a => a -> a -> Bool
<= Word
mx) forall a b. (a -> b) -> a -> b
$ forall b a. (b -> a -> b) -> b -> [a] -> [b]
scanl forall a. Num a => a -> a -> a
(+) Word
5 forall a b. (a -> b) -> a -> b
$ 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 forall a. Ord a => a -> a -> Bool
< a
k     = (t
m,a
e)
        | Bool
otherwise = case 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
eforall a. Num a => a -> a -> a
*a
k) (a
b forall a. Integral a => a -> a -> a
`quot` a
k) t
r (a
kforall 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 :: forall (n :: Natural).
KnownNat n =>
Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation Word
b1 Word
b2 Mod n
s = case Integer -> Integer -> Maybe SomePoint
newPoint (forall a. Integral a => a -> Integer
toInteger (forall (m :: Natural). Mod m -> Natural
unMod Mod n
s)) Integer
n of
  Maybe SomePoint
Nothing             -> 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 = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl (forall a b c. (a -> b -> c) -> b -> a -> c
flip forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply) Point a24 n
p0 [Word]
smallPowers
        z :: Integer
z = forall (a24 :: Natural) (n :: Natural). Point a24 n -> Integer
pointZ Point a24 n
q

    case 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 forall a. Integral a => a -> a -> a
gcd Integer
n (forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Point a24 n -> Word -> Word -> Integer
bigStep Point a24 n
q Word
b1 Word
b2) of
        Integer
1 -> forall a. Maybe a
Nothing
        Integer
g -> forall a. a -> Maybe a
Just Integer
g
      Integer
g -> forall a. a -> Maybe a
Just Integer
g
  where
    n :: Integer
n = forall a. Integral a => a -> Integer
toInteger (forall (n :: Natural) (proxy :: Natural -> *).
KnownNat n =>
proxy n -> Natural
natVal Mod n
s)
    smallPowers :: [Word]
smallPowers
      = forall a b. (a -> b) -> [a] -> [b]
map Word -> Word
findPower
      forall a b. (a -> b) -> a -> b
$ forall a. (a -> Bool) -> [a] -> [a]
takeWhile (forall a. Ord a => a -> a -> Bool
<= Word
b1) (Word
2 forall a. a -> [a] -> [a]
: Word
3 forall a. a -> [a] -> [a]
: Word
5 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 forall a. Ord a => a -> a -> Bool
<= Word
b1 forall a. Integral a => a -> a -> a
`quot` Word
p = Word -> Word
go (Word
acc 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 :: forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Point a24 n -> Word -> Word -> Integer
bigStep Point a24 n
q Word
b1 Word
b2 = Integer
rs
  where
    n :: Integer
n = forall (a24 :: Natural) (n :: Natural).
KnownNat n =>
Point a24 n -> Integer
pointN Point a24 n
q

    b0 :: Word
b0 = Word
b1 forall a. Num a => a -> a -> a
- Word
b1 forall a. Integral a => a -> a -> a
`rem` Word
wheel
    qks :: [(Integer, Point a24 n)]
qks = forall a b. [a] -> [b] -> [(a, b)]
zip [Integer
0..] forall a b. (a -> b) -> a -> b
$ forall a b. (a -> b) -> [a] -> [b]
map (forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
`multiply` Point a24 n
q) [Word]
wheelCoprimes
    qs :: [(Word, Point a24 n)]
qs = forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Point a24 n -> Word -> Word -> Word -> [(Word, Point a24 n)]
enumAndMultiplyFromThenTo Point a24 n
q Word
b0 (Word
b0 forall a. Num a => a -> a -> a
+ Word
wheel) Word
b2

    rs :: Integer
rs = forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (\Integer
ts (Word
_cHi, Point a24 n
p) -> forall (t :: * -> *) b a.
Foldable t =>
(b -> a -> b) -> b -> t a -> b
foldl' (\Integer
us (Integer
_cLo, Point a24 n
pq) ->
        Integer
us forall a. Num a => a -> a -> a
* (forall (a24 :: Natural) (n :: Natural). Point a24 n -> Integer
pointZ Point a24 n
p forall a. Num a => a -> a -> a
* forall (a24 :: Natural) (n :: Natural). Point a24 n -> Integer
pointX Point a24 n
pq forall a. Num a => a -> a -> a
- forall (a24 :: Natural) (n :: Natural). Point a24 n -> Integer
pointX Point a24 n
p forall a. Num a => a -> a -> a
* forall (a24 :: Natural) (n :: Natural). Point a24 n -> Integer
pointZ Point a24 n
pq) 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 forall a. Integral a => a -> a -> a
`div` Word
2], Word
k forall a. Integral a => a -> a -> a
`gcd` Word
wheel 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 :: forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Point a24 n -> Word -> Word -> Word -> [(Word, Point a24 n)]
enumAndMultiplyFromThenTo Point a24 n
p Word
from Word
thn Word
to = forall a b. [a] -> [b] -> [(a, b)]
zip [Word
from, Word
thn .. Word
to] [Point a24 n]
progression
  where
    step :: Word
step = Word
thn forall a. Num a => a -> a -> a
- Word
from

    pFrom :: Point a24 n
pFrom = forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply Word
from Point a24 n
p
    pThen :: Point a24 n
pThen = forall (a24 :: Natural) (n :: Natural).
(KnownNat a24, KnownNat n) =>
Word -> Point a24 n -> Point a24 n
multiply Word
thn  Point a24 n
p
    pStep :: Point a24 n
pStep = forall (a24 :: Natural) (n :: Natural).
(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 forall a. a -> [a] -> [a]
: Point a24 n
pThen forall a. a -> [a] -> [a]
: forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith (forall (n :: Natural) (a24 :: Natural).
KnownNat n =>
Point a24 n -> Point a24 n -> Point a24 n -> Point a24 n
`add` Point a24 n
pStep) [Point a24 n]
progression (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 = forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat [[Word
off forall a. Num a => a -> a -> a
+ forall a. Num a => Int -> a
toPrim Int
i | Int
i <- [Int
0 .. Int
li], 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) = forall (a :: * -> * -> *) e i.
(IArray a e, Ix i) =>
a i e -> (i, i)
bounds UArray Int Bool
bs; off :: Word
off = 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## -> 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#) 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) 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
xforall 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 forall a. Ord a => a -> a -> Bool
>= Int
smallPrimesLength
      = ([], forall a. a -> Maybe a
Just (BigNat -> Natural
NatJ# BigNat
m))
      | Bool
otherwise
      = let p# :: Word#
p# =
#if MIN_VERSION_base(4,16,0)
              Word16# -> Word#
word16ToWord#
#endif
              (Addr# -> Int# -> Word16#
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 #) = forall {a}. Num a => a -> BigNat -> (# a, BigNat #)
splitOff Word
1 BigNat
mp in
            (Word# -> Natural
NatS# Word#
p#, Word
k) forall {a} {b}. a -> ([a], b) -> ([a], b)
<: BigNat -> Int -> ([(Natural, Word)], Maybe Natural)
goBigNat BigNat
r (Int
i 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 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 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
_ = ([], forall a. Maybe a
Nothing)
    goWord Word#
m#  !Int
i
      | Int
i 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)], forall a. Maybe a
Nothing)
        else ([], 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)
              Word16# -> Word#
word16ToWord#
#endif
              (Addr# -> Int# -> Word16#
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)], 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#) forall {a} {b}. a -> ([a], b) -> ([a], b)
<: Word# -> Int -> ([(Natural, Word)], Maybe Natural)
goWord Word#
r# (Int
i 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 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 = 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 = forall b a. b -> (a -> b) -> Maybe a -> b
maybe (Word
wheel, Word
1000, Word
7) forall a b. (a, b) -> b
snd (forall a. Int -> IntMap a -> Maybe (Int, a)
IM.lookupLT Int
digs IntMap (Word, Word, Word)
testParms)