module Math.NumberTheory.Primes.Factorisation.Montgomery
(
factorise
, defaultStdGenFactorisation
, factorise'
, stepFactorisation
, defaultStdGenFactorisation'
, smallFactors
, stdGenFactorisation
, curveFactorisation
, montgomeryFactorisation
, findParms
) where
#include "MachDeps.h"
import GHC.Base
#if __GLASGOW_HASKELL__ < 705
import GHC.Word
#endif
import System.Random
import Control.Monad.State.Strict
#if __GLASGOW_HASKELL__ < 709
import Control.Applicative
#endif
import Data.Bits
import Data.Maybe
import Math.NumberTheory.Logarithms
import Math.NumberTheory.Logarithms.Internal
import Math.NumberTheory.Powers.General (highestPower, largePFPower)
import Math.NumberTheory.Powers.Squares (integerSquareRoot')
import Math.NumberTheory.Primes.Sieve.Eratosthenes
import Math.NumberTheory.Primes.Sieve.Indexing
import Math.NumberTheory.Primes.Testing.Probabilistic
import Math.NumberTheory.Unsafe
import Math.NumberTheory.Utils
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
factorise' :: Integer -> [(Integer,Int)]
factorise' n = defaultStdGenFactorisation' (mkStdGen $ fromInteger n `xor` 0xdeadbeef) n
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 < (m1) then k else error "Curves exhausted",k+1)) 6 Nothing r
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
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
stdGenFactorisation :: Maybe Integer
-> StdGen
-> Maybe Int
-> Integer
-> [(Integer,Int)]
stdGenFactorisation primeBound sg digits n
= curveFactorisation primeBound bailliePSW (\m -> randomR (6,m2)) sg digits n
curveFactorisation :: Maybe Integer
-> (Integer -> Bool)
-> (Integer -> g -> (Integer,g))
-> g
-> Maybe Int
-> Integer
-> [(Integer,Int)]
curveFactorisation primeBound primeTest prng seed mbdigs n
| ptest n = [(n,1)]
| otherwise = evalState (fact n digits) seed
where
digits = fromMaybe 8 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)
perfPw = case primeBound of
Nothing -> highestPower
Just bd -> largePFPower (integerSquareRoot' bd)
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 = case perfPw m of
(_,1) -> workFact m b1 b2 count
(b,e)
| ptest b -> return ([(b,e)],[])
| otherwise -> do
(as,bs) <- workFact b b1 b2 count
return $ (mult e as, mult e bs)
workFact m b1 b2 count
| count < 0 = return ([],[(m,1)])
| otherwise = do
s <- rndR m
case montgomeryFactorisation m b1 b2 s of
Nothing -> workFact m b1 b2 (count1)
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 (count1)
(cp,cc) <- if ptest cof
then return ([(cof,1)],[])
else repFact cof b1 b2 (count1)
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 (count1)
(cp,cc) <- if ptest c'
then return ([(c',1)],[])
else repFact c' b1 b2 (count1)
(gp,gc) <- if ptest g
then return ([(g,2)],[])
else dbl <$> repFact g b1 b2 (count1)
return (mergeAll [dp,cp,gp], dc ++ cc ++ gc)
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 (k1) (dbl pt)
p2 = dbln l2 p0
#if WORD_SIZE_IN_BITS == 64
mul a b c = (a*b) `quot` c
#else
mul a b c = fromInteger ((toInteger a * b) `quot` c)
#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
data Curve = C !Integer !Integer
data Point = P !Integer !Integer
montgomeryData :: Integer -> Integer -> (Curve, Point)
montgomeryData n s = (C an ad4, P x z)
where
u = (s*s5) `mod` n
v = (4*s) `mod` n
d = (vu)
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
add :: Integer -> Point -> Point -> Point -> Point
add n (P x0 z0) (P x1 z1) (P x2 z2) = P x3 z3
where
a = (x1z1)*(x2+z2)
b = (x1+z1)*(x2z2)
c = a+b
d = ab
x3 = (c*c*z0) `rem` n
z3 = (d*d*x0) `rem` n
double :: Integer -> Curve -> Point -> Point
double n (C an ad4) (P x z) = P x' z'
where
r = x+z
s = xz
u = r*r
v = s*s
t = uv
x' = (ad4*u*v) `rem` n
z' = ((ad4*v+t*an)*t) `rem` n
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)
primeStore :: [PrimeSieve]
primeStore = psieveFrom 7
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 :: 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)
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)
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