module Math.NumberTheory.Primes.Factorisation.Montgomery
(
factorise
, defaultStdGenFactorisation
, factorise'
, stepFactorisation
, defaultStdGenFactorisation'
, smallFactors
, stdGenFactorisation
, curveFactorisation
, montgomeryFactorisation
, findParms
) where
#include "MachDeps.h"
import System.Random
import Control.Monad.State.Strict
#if __GLASGOW_HASKELL__ < 709
import Control.Applicative
import Data.Word
#endif
import Data.Bits
import Data.IntMap (IntMap)
import qualified Data.IntMap as IM
import Data.List (foldl')
import Data.Maybe
import GHC.TypeNats.Compat
import Math.NumberTheory.Curves.Montgomery
import Math.NumberTheory.Moduli.Class
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+5 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 s `modulo` fromInteger m of
InfMod{} -> error "impossible case"
SomeMod sm -> case montgomeryFactorisation b1 b2 sm 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 :: KnownNat n => Word -> Word -> Mod n -> Maybe Integer
montgomeryFactorisation b1 b2 s = case newPoint (getVal s) n of
Nothing -> Nothing
Just (SomePoint p0) -> do
let q = foldl (flip multiply) p0 smallPowers
z = pointZ q
fromIntegral <$> case gcd n z of
1 -> case gcd n (bigStep q b1 b2) of
1 -> Nothing
g -> Just g
g -> Just g
where
n = getMod s
smallPrimes = takeWhile (<= b1) (2 : 3 : 5 : list primeStore)
smallPowers = map findPower smallPrimes
findPower p = go p
where
go acc
| acc <= b1 `quot` p = go (acc * p)
| otherwise = acc
bigStep :: (KnownNat a24, KnownNat n) => Point a24 n -> Word -> Word -> Integer
bigStep q b1 b2 = rs
where
n = pointN q
b0 = b1 b1 `rem` wheel
qks = zip [0..] $ map (\k -> multiply k q) wheelCoprimes
qs = enumAndMultiplyFromThenTo q b0 (b0 + wheel) b2
rs = foldl' (\ts (_cHi, p) -> foldl' (\us (_cLo, pq) ->
us * (pointZ p * pointX pq pointX p * pointZ pq) `rem` n
) ts qks) 1 qs
wheel :: Word
wheel = 210
wheelCoprimes :: [Word]
wheelCoprimes = [ k | k <- [1 .. wheel `div` 2], k `gcd` wheel == 1 ]
enumAndMultiplyFromThenTo
:: (KnownNat a24, KnownNat n)
=> Point a24 n
-> Word
-> Word
-> Word
-> [(Word, Point a24 n)]
enumAndMultiplyFromThenTo p from thn to = zip [from, thn .. to] progression
where
step = thn from
pFrom = multiply from p
pThen = multiply thn p
pStep = multiply step p
progression = pFrom : pThen : zipWith (\x0 x1 -> add x0 pStep x1) progression (tail progression)
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 :: (Ord a, Num b) => [(a, b)] -> [(a, b)] -> [(a, b)]
merge xs [] = xs
merge [] ys = ys
merge xxs@(x@(p, k) : xs) yys@(y@(q, m) : ys)
= case p `compare` q of
LT -> x : merge xs yys
EQ -> (p, k + m) : merge xs ys
GT -> y : merge xxs ys
mergeAll :: (Ord a, Num b) => [[(a, b)]] -> [(a, b)]
mergeAll = \case
[] -> []
[xs] -> xs
(xs : ys : zss) -> merge (merge xs ys) (mergeAll zss)
testParms :: IntMap (Word, Word, Word)
testParms = IM.fromList
[ (12, ( 400, 40000, 10))
, (15, ( 2000, 200000, 25))
, (20, ( 11000, 1100000, 90))
, (25, ( 50000, 5000000, 300))
, (30, ( 250000, 25000000, 700))
, (35, ( 1000000, 100000000, 1800))
, (40, ( 3000000, 300000000, 5100))
, (45, ( 11000000, 1100000000, 10600))
, (50, ( 43000000, 4300000000, 19300))
, (55, ( 110000000, 11000000000, 49000))
, (60, ( 260000000, 26000000000, 124000))
, (65, ( 850000000, 85000000000, 210000))
, (70, (2900000000, 290000000000, 340000))
]
findParms :: Int -> (Word, Word, Word)
findParms digs = maybe (wheel, 1000, 7) snd (IM.lookupLT digs testParms)