{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE DataKinds #-}
{-# LANGUAGE KindSignatures #-}
{-# LANGUAGE LambdaCase #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# OPTIONS_GHC -fno-warn-type-defaults #-}
{-# OPTIONS_HADDOCK hide #-}
module Math.NumberTheory.Primes.Factorisation.Montgomery
(
factorise
, defaultStdGenFactorisation
, factorise'
, stepFactorisation
, defaultStdGenFactorisation'
, smallFactors
, stdGenFactorisation
, curveFactorisation
, montgomeryFactorisation
, findParms
) where
import Control.Arrow
import Control.Monad.Trans.State.Lazy
import System.Random
import Data.Bits
import Data.IntMap (IntMap)
import qualified Data.IntMap as IM
import Data.List (foldl')
import Data.Maybe
#if __GLASGOW_HASKELL__ < 803
import Data.Semigroup
#endif
import Data.Traversable
import GHC.TypeNats.Compat
import Math.NumberTheory.Curves.Montgomery
import Math.NumberTheory.Euclidean.Coprimes (splitIntoCoprimes, unCoprimes)
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.Primes.Types (unPrime)
import Math.NumberTheory.Unsafe
import Math.NumberTheory.Utils
factorise :: Integer -> [(Integer, Word)]
factorise n
| abs n == 1 = []
| n < 0 = factorise (-n)
| n == 0 = error "0 has no prime factorisation"
| otherwise = factorise' n
factorise' :: Integer -> [(Integer, Word)]
factorise' n = defaultStdGenFactorisation' (mkStdGen $ fromInteger n `xor` 0xdeadbeef) n
stepFactorisation :: Integer -> [(Integer, Word)]
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 :: StdGen -> Integer -> [(Integer, Word)]
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, Word)]
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, Word)]
stdGenFactorisation primeBound sg digits n
= curveFactorisation primeBound bailliePSW (\m -> randomR (6,m-2)) sg digits n
curveFactorisation
:: forall g.
Maybe Integer
-> (Integer -> Bool)
-> (Integer -> g -> (Integer, g))
-> g
-> Maybe Int
-> Integer
-> [(Integer, Word)]
curveFactorisation primeBound primeTest prng seed mbdigs n
| n == 1 = []
| ptest n = [(n, 1)]
| otherwise = evalState (fact n digits) seed
where
digits :: Int
digits = fromMaybe 8 mbdigs
ptest :: Integer -> Bool
ptest = maybe primeTest (\bd k -> k <= bd || primeTest k) primeBound
rndR :: Integer -> State g Integer
rndR k = state (prng k)
perfPw :: Integer -> (Integer, Word)
perfPw = maybe highestPower (largePFPower . integerSquareRoot') primeBound
fact :: Integer -> Int -> State g [(Integer, Word)]
fact 1 _ = return mempty
fact m digs = do
let (b1, b2, ct) = findParms digs
Factors pfs cfs <- repFact m b1 b2 ct
case cfs of
[] -> return pfs
_ -> do
nfs <- forM cfs $ \(k, j) ->
map (second (* j)) <$> fact k (if null pfs then digs + 5 else digs)
return $ mconcat (pfs : nfs)
repFact :: Integer -> Word -> Word -> Word -> State g Factors
repFact 1 _ _ _ = return mempty
repFact m b1 b2 count =
case perfPw m of
(_, 1) -> workFact m b1 b2 count
(b, e)
| ptest b -> return $ singlePrimeFactor b e
| otherwise -> modifyPowers (* e) <$> workFact b b1 b2 count
workFact :: Integer -> Word -> Word -> Word -> State g Factors
workFact 1 _ _ _ = return mempty
workFact m _ _ 0 = return $ singleCompositeFactor m 1
workFact m b1 b2 count = 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 (count - 1)
Just d -> do
let cs = unCoprimes $ splitIntoCoprimes [(d, 1), (m `quot` d, 1)]
fmap mconcat $ flip mapM cs $
\(x, xm) -> if ptest x
then pure $ singlePrimeFactor x xm
else repFact x b1 b2 (count - 1)
data Factors = Factors
{ _primeFactors :: [(Integer, Word)]
, _compositeFactors :: [(Integer, Word)]
}
singlePrimeFactor :: Integer -> Word -> Factors
singlePrimeFactor a b = Factors [(a, b)] []
singleCompositeFactor :: Integer -> Word -> Factors
singleCompositeFactor a b = Factors [] [(a, b)]
instance Semigroup Factors where
Factors pfs1 cfs1 <> Factors pfs2 cfs2
= Factors (pfs1 <> pfs2) (cfs1 <> cfs2)
instance Monoid Factors where
mempty = Factors [] []
mappend = (<>)
modifyPowers :: (Word -> Word) -> Factors -> Factors
modifyPowers f (Factors pfs cfs)
= Factors (map (second f) pfs) (map (second f) cfs)
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, Word)], 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 = map unPrime $ 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)
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)