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 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.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.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