prime factors decomposition

app/FactorizationUI.hs

@@ -2,6 +2,7 @@ module FactorizationUI (run) where
 
 import System.Exit (exitSuccess)
 import Utils (askChoice, askNumber)
+import PrimeFactors (primeFactors, formatFactors)
 import Factorization
     ( fermatFactorization, pollardPminus1, dixon, dixonDefaultParams
     , gnfsWith, GNFSParams(..), defaultParams
@@ -14,14 +15,16 @@ run = do
     putStrLn "2) Pollard p-1"
     putStrLn "3) Dixon"
     putStrLn "4) General Number Field Sieve (GNFS)"
+    putStrLn "5) Prime factors (GNFS)"
 
-    choice <- askChoice 4
+    choice <- askChoice 5
 
     case choice of
         1 -> fermatUI
         2 -> pollardUI
         3 -> dixonUI
         4 -> gnfsUI
+        5 -> primeFactorsUI
         _ -> error "Impossible"
 
 fermatUI :: IO ()
@@ -47,7 +50,7 @@ dixonUI :: IO ()
 dixonUI = do
     putStrLn "Dixon"
     putStrLn "Factors a composite n (not a prime, a perfect square or a prime power)."
-    putStrLn "Practical for n up to roughly 22 digits; 24 digits already takes about a minute."
+    putStrLn "Practical for n up to roughly 22 digits; 24 digits already takes a minute or two."
     n <- askNumber "Enter an integer n > 1:"
 
     putStrLn "Parameters:"
@@ -118,3 +121,17 @@ gnfsUI = do
             putStrLn ("n = " ++ show n)
             putStrLn ("p = " ++ show p)
             putStrLn ("q = " ++ show q)
+
+primeFactorsUI :: IO ()
+primeFactorsUI = do
+    putStrLn "Prime factors (GNFS)"
+    putStrLn "Decomposes n into its prime factors."
+    putStrLn "Large composites are split with the General Number Field Sieve, so numbers with more than about 30 digits may take very long."
+    n <- askNumber "Enter an integer n > 1:"
+
+    putStrLn "Running (this may take a while)..."
+    case primeFactors n of
+        Left err ->
+            putStrLn ("Error: " ++ err)
+        Right factors ->
+            putStrLn (show n ++ " = " ++ formatFactors factors)

haskell-math.cabal

@@ -24,6 +24,7 @@ library
         ModularSquareRoot
         ModularArithmeticUtils
         LinearAlgebra.GF2
+        PrimeFactors
     other-modules:
         Primes.FermatPrimeTest
         Primes.MillerRabin

src/PrimeFactors.hs

@@ -0,0 +1,77 @@
+module PrimeFactors
+  ( primeFactors
+  , formatFactors
+  , expandFactors
+  ) where
+
+import Data.List (intercalate)
+import qualified Data.Map.Strict as M
+
+import Factorization.GNFS (gnfs)
+import GNFS.Util (isPrimeDet, perfectPower)
+import ModularArithmeticUtils (countDiv)
+import Primes.ErathosteneSieve (erathosteneSieve)
+
+-- | A factorisation as a list of (prime, exponent), sorted by increasing prime.
+type Factors = [(Integer, Int)]
+
+-- | Decompose n >= 1 into its prime factors, e.g. 360 -> [(2,3),(3,2),(5,1)].
+-- 1 gives the empty factorisation.
+--
+-- Small primes are removed by trial division; what is left is split recursively with GNFS ('gnfs'): primes are recognised with a strong probable-prime test (deterministic  below about 3.3e24), perfect powers are reduced to their root, and every other composite is cut in two by GNFS and both halves are factored further.
+--
+-- Returns Left if GNFS fails on one of the composite cofactors.
+primeFactors :: Integer -> Either String Factors
+primeFactors n
+  | n < 1     = Left "n must be a positive integer"
+  | otherwise =
+      let (small, m) = stripSmallPrimes n
+      in if m == 1
+           then Right small
+           else merge small <$> factorLarge m
+
+-- | Primes below this bound are removed by trial division before GNFS is used.
+trialBound :: Integer
+trialBound = 10000
+
+stripSmallPrimes :: Integer -> (Factors, Integer)
+stripSmallPrimes n = go (erathosteneSieve trialBound) n []
+  where
+    go _ 1 acc = (reverse acc, 1)
+    go [] m acc = (reverse acc, m)
+    go (p : ps) m acc
+      | p * p > m = if m > 1 then (reverse ((m, 1) : acc), 1) else (reverse acc, 1)
+      | otherwise =
+          case countDiv p m of
+            (0, _)  -> go ps m acc
+            (e, m') -> go ps m' ((p, fromInteger e) : acc)
+
+-- | Factor m > 1 that has no prime factor below 'trialBound'.
+factorLarge :: Integer -> Either String Factors
+factorLarge m
+  | isPrimeDet m = Right [(m, 1)]
+  | Just (r, k) <- perfectPower m =
+      map (\(p, e) -> (p, e * k)) <$> factorLarge r
+  | otherwise =
+      case gnfs m of
+        Left err -> Left ("could not factor " ++ show m ++ ": " ++ err)
+        Right (p, q)
+          | p > 1 && q > 1 && p * q == m -> merge <$> factorLarge p <*> factorLarge q
+          | otherwise -> Left ("GNFS returned an invalid split of " ++ show m)
+
+-- | Combine two factorisations, adding the exponents of equal primes.
+merge :: Factors -> Factors -> Factors
+merge a b = M.toAscList (M.fromListWith (+) (a ++ b))
+
+-- | Render a factorisation as "2^2*3*5": exponents of 1 are omitted, and 1 (the empty
+-- factorisation) is written "1".
+formatFactors :: Factors -> String
+formatFactors [] = "1"
+formatFactors fs = intercalate "*" (map one fs)
+  where
+    one (p, 1) = show p
+    one (p, e) = show p ++ "^" ++ show e
+
+-- | Multiply a factorisation back out.
+expandFactors :: Factors -> Integer
+expandFactors = product . map (\(p, e) -> p ^ e)