Dixon factorization

app/FactorizationUI.hs

@@ -1,19 +1,21 @@
 module FactorizationUI (run) where
 
 import Utils (askChoice, askNumber)
-import Factorization (fermatFactorization, pollardPminus1)
+import Factorization (fermatFactorization, pollardPminus1, dixon)
 
 run :: IO ()
 run = do
     putStrLn "Factorization"
     putStrLn "1) Fermat factorization"
     putStrLn "2) Pollard p-1"
+    putStrLn "3) Dixon"
 
-    choice <- askChoice 2
+    choice <- askChoice 3
 
     case choice of
         1 -> fermatUI
         2 -> pollardUI
+        3 -> dixonUI
         _ -> error "Impossible"
 
 fermatUI :: IO ()
@@ -34,3 +36,13 @@ pollardUI = do
             putStrLn ("Other factor: " ++ show (n `div` factor))
         Nothing ->
             putStrLn "Failed to find a factor (try a different B or method)."
+
+dixonUI :: IO ()
+dixonUI = do
+    n <- askNumber "Enter an integer n > 1:"
+    b <- askNumber "Enter the bound B:"
+    x <- askNumber "Enter the max random integer x:"
+    let (p, q) = dixon n b x
+    putStrLn ("n = " ++ show n)
+    putStrLn ("p = " ++ show p)
+    putStrLn ("q = " ++ show q)

haskell-math.cabal

@@ -31,6 +31,7 @@ library
         Primes.ErathosteneSieve
         Factorization.FermatFactorization
         Factorization.PollardPminus1
+        Factorization.Dixon
         ModularSquareRoot.TonelliShanks
     build-depends:
         base ^>=4.18.2.1,

src/Factorization.hs

@@ -1,4 +1,5 @@
-module Factorization ( fermatFactorization, pollardPminus1) where
+module Factorization ( fermatFactorization, pollardPminus1, dixon) where
 
 import Factorization.FermatFactorization
 import Factorization.PollardPminus1
+import Factorization.Dixon

src/Factorization/Dixon.hs

@@ -0,0 +1,53 @@
+module Factorization.Dixon (dixon) where
+
+import Primes.ErathosteneSieve (erathosteneSieve)
+import LinearAlgebra.GF2 (fromBools, gaussianEliminationMask, maskToIndicesInt)
+import ModularArithmeticUtils (modExp, modMul, countDiv)
+
+
+factorOverBase :: [Integer] -> Integer -> Maybe [Integer]
+factorOverBase primes a = go primes a []
+  where
+    go [] 1 acc = Just (reverse acc)
+    go [] _ acc = Nothing
+    go (p:ps) num acc = let (e, num') = countDiv p num
+                        in go ps num' (e : acc)
+
+
+dixon :: Integer -> Integer -> Integer -> (Integer, Integer)
+dixon n b maxX
+  | n <= 3    = error "n must be greater than 3"
+  | even n    = (2, n `div` 2)
+  | otherwise =
+      let base = erathosteneSieve b
+          numCols = length base
+          floorSqrt m
+            | m < 2 = m
+            | otherwise = go m
+            where go k = let q = (k + m `div` k) `div` 2
+                         in if q < k then go q else q
+          startX = floorSqrt n + 1
+          pairs = [(x, exps) | x <- [startX .. maxX],
+                               let a = (x * x) `mod` n,
+                               a > 1,
+                               Just exps <- [factorOverBase base a]]
+          numCollected = length pairs
+      in if numCollected <= numCols
+         then error $ "Not enough smooth squares found up to " ++ show maxX
+         else let xs = map fst pairs
+                  allExps = map snd pairs
+                  bitvecs = map (fromBools . map (\e -> e `mod` 2 == 1)) allExps
+              in case gaussianEliminationMask numCols bitvecs of
+                   Nothing -> error "No linear dependency found"
+                   Just mask ->
+                     let indices = maskToIndicesInt mask
+                     in if null indices
+                        then error "Empty dependency"
+                        else let productX = foldl (\acc idx -> modMul acc (xs !! idx) n) 1 indices
+                                 fs = [ sum [ (allExps !! idx !! j) | idx <- indices ] `div` 2 | j <- [0 .. numCols - 1] ]
+                                 z = foldl (\acc (p, f) -> if f == 0 then acc else modMul (modExp p f n) acc n) 1 (zip base fs)
+                                 diff = productX - z
+                                 g = gcd (abs diff) n
+                             in if g == 1 || g == n
+                                then error "non trivial factors not found, try larger b or maximum x"
+                                else (g, n `div` g)

src/ModularArithmeticUtils.hs

@@ -1,4 +1,4 @@
-module ModularArithmeticUtils (modExp, modMul, legendre, factorOutTwos, jacobi) where
+module ModularArithmeticUtils (modExp, modMul, legendre, factorOutTwos, jacobi, countDiv) where
 
 import Data.Bits (testBit, shiftR)
 
@@ -58,3 +58,11 @@ jacobi2 m =
     3 -> -1
     5 -> -1
     _ -> error "jacobi2: unexpected (m mod 8) for odd m"
+
+
+-- divide q by m as many times
+-- returns the number of division made and the remainder
+countDiv :: Integer -> Integer -> (Integer, Integer)
+countDiv q m
+  | m < 2 || m `mod` q /= 0 = (0, m)
+  | otherwise = let (ee, mm) = countDiv q (m `div` q) in (ee + 1, mm)