Commits

dafis committed 390926f

Elliptic curve factorisation after Montgomery

  • Participants
  • Parent commits 75b2c9f

Comments (0)

Files changed (1)

File Math/NumberTheory/Primes/Factorisation/Montgomery.hs

+-- |
+-- Module:      Math.NumberTheory.Primes.Factorisation.Montgomery
+-- Copyright:   (c) 2011 Daniel Fischer
+-- Licence:     MIT
+-- Maintainer:  Daniel Fischer <daniel.is.fischer@googlemail.com>
+-- Stability:   Provisional
+-- Portability: Non-portable (GHC extensions)
+--
+-- 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/>
+--
+-- The implementation is not very optimised, so it is not suitable for factorising numbers
+-- with only huge prime divisors. However, factors of 20-25 digits are normally found in
+-- acceptable time. The time taken depends, however, strongly on how lucky the curve-picking
+-- is. With luck, even large factors can be found in seconds; on the other hand, finding small
+-- factors (about 10 digits) can take minutes when the curve-picking is bad.
+--
+-- Given enough time, the algorithm should be able to factor numbers of 100-120 digits, but it
+-- is best suited for numbers of up to 50-60 digits.
+{-# LANGUAGE CPP, BangPatterns, MagicHash #-}
+{-# OPTIONS_HADDOCK hide #-}
+module Math.NumberTheory.Primes.Factorisation.Montgomery
+  ( -- *  Complete factorisation functions
+    -- ** Functions with input checking
+    factorise
+  , defaultStdGenFactorisation
+    -- ** Functions without input checking
+  , factorise'
+  , stepFactorisation
+  , defaultStdGenFactorisation'
+    -- * Partial factorisation
+  , smallFactors
+  , stdGenFactorisation
+  , curveFactorisation
+    -- ** Single curve worker
+  , montgomeryFactorisation
+  ) where
+
+#include "MachDeps.h"
+
+import GHC.Base
+import GHC.Word
+import Data.Array.Base
+
+import System.Random
+import Control.Monad.State.Strict
+import Control.Applicative
+import Data.Bits
+import Data.Maybe
+
+import Math.NumberTheory.Logarithms
+import Math.NumberTheory.Logarithms.Internal
+import Math.NumberTheory.Primes.Sieve.Eratosthenes
+import Math.NumberTheory.Primes.Sieve.Indexing
+import Math.NumberTheory.Primes.Testing.Probabilistic
+import Math.NumberTheory.Utils
+
+-- | @'factorise' n@ produces the prime factorisation of @n@, including
+--   a factor of @(-1)@ if @n < 0@. @'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@.
+factorise :: Integer -> [(Integer,Int)]
+factorise n
+    | n < 0     = (-1,1):factorise (-n)
+    | n == 0    = error "0 has no prime factorisation"
+    | n == 1    = []
+    | otherwise = factorise' n
+
+-- | Like 'factorise', but without input checking, hence @n > 1@ is required.
+factorise' :: Integer -> [(Integer,Int)]
+factorise' n = defaultStdGenFactorisation' (mkStdGen $ fromInteger n `xor` 0xdeadbeef) n
+
+-- | @'stepFactorisation'@ is like 'factorise\'', except that it doesn't use a
+--   pseudo random generator but steps through the curves in order.
+--   This strategy turns out to be surprisingly fast, on average it doesn't
+--   seem to be slower than the 'StdGen' based variant.
+stepFactorisation :: Integer -> [(Integer,Int)]
+stepFactorisation n
+    = let (sfs,mb) = smallFactors 100000 n
+      in sfs ++ case mb of
+                  Nothing -> []
+                  Just r  -> curveFactorisation (Just 10000000000) bailliePSW
+                                                (\m k -> (if k < (m-1) then k else error "Curves exhausted",k+1)) 6 Nothing r
+
+-- | @'defaultStdGenFactorisation'@ first strips off all small prime factors and then,
+--   if the factorisation is not complete, proceeds to curve factorisation.
+--   For negative numbers, a factor of @-1@ is included, the factorisation of @1@
+--   is empty. Since @0@ has no prime factorisation, a zero argument causes
+--   an error.
+defaultStdGenFactorisation :: StdGen -> Integer -> [(Integer,Int)]
+defaultStdGenFactorisation sg n
+    | n == 0    = error "0 has no prime factorisation"
+    | n < 0     = (-1,1) : defaultStdGenFactorisation sg (-n)
+    | n == 1    = []
+    | otherwise = defaultStdGenFactorisation' sg n
+
+-- | Like 'defaultStdGenFactorisation', but without input checking, so
+--   @n@ must be larger than @1@.
+defaultStdGenFactorisation' :: StdGen -> Integer -> [(Integer,Int)]
+defaultStdGenFactorisation' sg n
+    = let (sfs,mb) = smallFactors 100000 n
+      in sfs ++ case mb of
+                  Nothing -> []
+                  Just m  -> stdGenFactorisation (Just 10000000000) sg Nothing 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,Int)]  -- ^ List of prime factors and exponents
+stdGenFactorisation primeBound sg digits n
+    = curveFactorisation primeBound bailliePSW (\m -> randomR (6,m-2)) sg digits n
+
+-- | @'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'@ requires that small prime factors have been stripped before. Also, it is
+--   unlikely to succeed if @n@ has more than one (really) large prime factor.
+curveFactorisation :: 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,Int)]               -- ^ List of prime factors and exponents
+curveFactorisation primeBound primeTest prng seed mbdigs n
+    | ptest n   = [(n,1)]
+    | otherwise = evalState (fact n digits) seed
+      where
+        digits = fromMaybe 6 mbdigs
+        mult 1 xs = xs
+        mult j xs = [(p,j*k) | (p,k) <- xs]
+        dbl (u,v) = (mult 2 u, mult 2 v)
+        ptest = case primeBound of
+                  Just bd -> \k -> k <= bd || primeTest k
+                  Nothing -> primeTest
+        rndR k = state (\gen -> prng k gen)
+        fact m digs = do let (b1,b2,ct) = findParms digs
+                         (pfs,cfs) <- repFact m b1 b2 ct
+                         if null cfs
+                           then return pfs
+                           else do
+                               nfs <- forM cfs $ \(k,j) ->
+                                   mult j <$> fact k (if null pfs then digs+4 else digs)
+                               return (mergeAll $ pfs:nfs)
+        repFact m b1 b2 count
+            | count < 0 = return ([],[(m,1)])
+            | otherwise = do
+                s <- rndR m
+                case montgomeryFactorisation m b1 b2 s of
+                  Nothing -> repFact m b1 b2 (count-1)
+                  Just d  -> do
+                      let !cof = m `quot` d
+                      case gcd cof d of
+                        1 -> do
+                            (dp,dc) <- if ptest d
+                                         then return ([(d,1)],[])
+                                         else repFact d b1 b2 (count-1)
+                            (cp,cc) <- if ptest cof
+                                         then return ([(cof,1)],[])
+                                         else repFact cof b1 b2 (count-1)
+                            return (merge dp cp, dc ++ cc)
+                        g -> do
+                            let d' = d `quot` g
+                                c' = cof `quot` g
+                            (dp,dc) <- if ptest d'
+                                         then return ([(d',1)],[])
+                                         else repFact d' b1 b2 (count-1)
+                            (cp,cc) <- if ptest c'
+                                         then return ([(c',1)],[])
+                                         else repFact c' b1 b2 (count-1)
+                            (gp,gc) <- if ptest g
+                                         then return ([(g,2)],[])
+                                         else dbl <$> repFact g b1 b2 (count-1)
+                            return  (mergeAll [dp,cp,gp], dc ++ cc ++ gc)
+
+----------------------------------------------------------------------------------------------------
+--                                         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 :: Integer -> Word -> Word -> Integer -> Maybe Integer
+montgomeryFactorisation n b1 b2 s = go p5 (list primeStore)
+  where
+    l2 = wordLog2' b1
+    b1i = toInteger b1
+    (^~) :: Word -> Int -> Word
+    w ^~ i = w ^ i
+    (e, p0) = montgomeryData n s
+    dbl pt = double n e pt
+    dbln 0 !pt = pt
+    dbln k pt = dbln (k-1) (dbl pt)
+    p2 = dbln l2 p0
+#if WORD_SIZE_IN_BITS == 64
+    mul a b c = (a*b) `quot` c       -- can't overflow, work on Int
+#else
+    mul a b c = fromInteger ((toInteger a * b) `quot` c) -- might overflow if Int is used
+#endif
+    adjust bd ml w
+      | w <= bd     = adjust bd ml (w*ml)
+      | otherwise   = w
+    l3 = mul l2 190537 301994
+    w3 = 3 ^~ l3
+    pw3 = adjust (b1 `quot` 3) 3 w3
+    p3 = multiply n e pw3 p2
+    l5 = mul l2 1936274 4495889
+    w5 = 5 ^~ l5
+    pw5 = adjust (b1 `quot` 5) 5 w5
+    p5 = multiply n e pw5 p3
+    go (P _ 0) _ = Nothing
+    go !pt@(P _ z) (pr:prs)
+      | pr <= b1    = let !lp = integerLogBase' (fromIntegral pr) b1i
+                      in go (multiply n e (pr ^~ lp) pt) prs
+      | otherwise   = case gcd n z of
+                        1 -> lgo (multiply n e pr pt) prs
+                        g -> Just g
+    go (P _ z) _    = case gcd n z of
+                        1 -> Nothing
+                        g -> Just g
+    lgo (P _ 0) _ = Nothing
+    lgo !pt@(P _ z) (pr:prs)
+      | pr <= b2    = lgo (multiply n e pr pt) prs
+      | otherwise   = case gcd n z of
+                        1 -> Nothing
+                        g -> Just g
+    lgo (P _ z) _   = case gcd n z of
+                        1 -> Nothing
+                        g -> Just g
+
+----------------------------------------------------------------------------------------------------
+--                            Helpers, Curves and elliptic arithmetics                            --
+----------------------------------------------------------------------------------------------------
+
+-- A Montgomery curve is given by y^2 = x^3 + (A_n / A_d - 2)*x^2 + x (mod n).
+-- We store A_n and 4*A_d, since A_n occurs with the factor 4 in all formulae.
+data Curve = C !Integer !Integer
+
+-- Point in the projective plane, will be on the curve
+-- A coordinate transformation eliminates the y-coordinate, hence
+-- we store only x and z
+data Point = P !Integer !Integer
+
+-- Get curve and point to start
+-- Input should satisfy 6 <= s < n-1
+montgomeryData :: Integer -> Integer -> (Curve, Point)
+montgomeryData n s = (C an ad4, P x z)
+  where
+    u = (s*s-5) `mod` n
+    v = (4*s) `mod` n
+    d = (v-u)
+    x = (u*u*u) `mod` n
+    z = (v*v*v) `mod` n
+    an = ((d*d)*(d*(3*u+v))) `mod` n
+    ad4 = (16*x*v) `mod` n
+
+-- Addition on the curve, given the modulus n and three points,
+-- p0, p1 and p2, with p0 = p2 - p1, calculate the point p1 + p2.
+-- Note that the addition does not depend on the curve.
+add :: Integer -> Point -> Point -> Point -> Point
+add n (P x0 z0) (P x1 z1) (P x2 z2) = P x3 z3
+  where
+    a = (x1-z1)*(x2+z2)
+    b = (x1+z1)*(x2-z2)
+    c = a+b
+    d = a-b
+    x3 = (c*c*z0) `rem` n
+    z3 = (d*d*x0) `rem` n
+
+-- Double a point on the curve.
+double :: Integer -> Curve -> Point -> Point
+double n (C an ad4) (P x z) = P x' z'
+  where
+    r = x+z
+    s = x-z
+    u = r*r
+    v = s*s
+    t = u-v
+    x' = (ad4*u*v) `rem` n
+    z' = ((ad4*v+t*an)*t) `rem` n
+
+-- Multiply a point on the curve by a Word.
+-- Within Word range, we can use the faster variant going
+-- from high-order bits to low-order.
+multiply :: Integer -> Curve -> Word -> Point -> Point
+multiply n cve (W# w##) p =
+    case wordLog2# w## of
+      l# -> go (l# -# 1#) p (double n cve p)
+  where
+    go 0# !p0 !p1 = case w## `and#` 1## of
+                      0## -> double n cve p0
+                      _   -> add n p p0 p1
+    go i# p0 p1 = case (uncheckedShiftRL# w## i#) `and#` 1## of
+                    0## -> go (i# -# 1#) (double n cve p0) (add n p p0 p1)
+                    _   -> go (i# -# 1#) (add n p p0 p1) (double n cve p1)
+
+{-  Not (yet) needed
+-- Multiply a point on the curve by an Integer.
+multIgr :: Integer -> Curve -> Integer -> Point -> Point
+multIgr n cve k p = go k
+  where
+    go 1 = (p, double n cve p)
+    go m = case m `quotRem` 2 of
+             (q,r) -> let !(!s, l) = go q
+                      in case r of
+                           0 -> (double n cve s, add n p s l)
+                           _ -> (add n p s l, double n cve l)
+-}
+
+-- primes, compactly stored as a bit sieve
+primeStore :: [PrimeSieve]
+primeStore = psieveFrom 7
+
+-- generate list of primes from arrays
+list :: [PrimeSieve] -> [Word]
+list sieves = concat [[off + toPrim i | i <- [0 .. li], unsafeAt bs i]
+                                | PS vO bs <- sieves, let { (_,li) = bounds bs; off = fromInteger vO; }]
+
+-- | @'smallFactors' bound n@ finds all prime divisors of @n > 1@ up to @bound@ by trial division and returns the
+--   list of these together with their multiplicities, and a possible remaining factor which may be composite.
+smallFactors :: Integer -> Integer -> ([(Integer,Int)], Maybe Integer)
+smallFactors bd n = case shiftToOddCount n of
+                      (0,m) -> go m prms
+                      (k,m) -> (2,k) <: if m == 1 then ([],Nothing) else go m prms
+  where
+    prms = tail (primeStore >>= primeList)
+    x <: ~(l,b) = (x:l,b)
+    go m (p:ps)
+        | m < p*p   = ([(m,1)], Nothing)
+        | bd < p    = ([], Just m)
+        | otherwise = case splitOff p m of
+                        (0,_) -> go m ps
+                        (k,r) | r == 1 -> ([(p,k)], Nothing)
+                              | otherwise -> (p,k) <: go r ps
+    go m [] = ([(m,1)], Nothing)
+
+-- helpers: merge sorted lists
+merge :: [(Integer,Int)] -> [(Integer,Int)] -> [(Integer,Int)]
+merge xxs@(x@(p,k):xs) yys@(y@(q,m):ys) = case compare p q of
+                                            LT -> x : merge xs yys
+                                            EQ -> (p,k+m) : merge xs ys
+                                            GT -> y : merge xxs ys
+merge xs [] = xs
+merge _ ys = ys
+
+mergeAll :: [[(Integer,Int)]] -> [(Integer,Int)]
+mergeAll [] = []
+mergeAll [xs] = xs
+mergeAll (xs:ys:zss) = merge (merge xs ys) (mergeAll zss)
+
+-- Parameters for the factorisation, the two b-parameters for montgomery and the number of tries
+-- to use these, depending on the size of the factor we are looking for.
+-- The numbers are roughly based on the parameters listed on Dario Alpern's ECM site.
+testParms :: [(Int,Word,Word,Int)]
+testParms = [ (12, 400, 10000, 10), (15, 2000, 50000, 25), (20, 11000, 150000, 90)
+            , (25, 50000, 500000, 300), (30, 250000, 1500000, 700)
+            , (35, 1000000, 4000000, 1800), (40, 3000000, 12000000, 5100)
+            , (45, 11000000, 45000000, 10600), (50, 43000000, 200000000, 19300)
+            , (55, 80000000, 400000000,30000), (60, 120000000, 700000000, 50000)
+            ]
+
+findParms :: Int -> (Word, Word, Int)
+findParms digs = go (100, 1000, 7) testParms
+  where
+    go p ((d,b1,b2,ct):rest)
+      | digs < d    = p
+      | otherwise   = go (b1,b2,ct) rest
+    go p [] = p