Dixon factorization, automatic parameters

app/FactorizationUI.hs

@@ -3,7 +3,7 @@ module FactorizationUI (run) where
 import System.Exit (exitSuccess)
 import Utils (askChoice, askNumber)
 import Factorization
-    ( fermatFactorization, pollardPminus1, dixon
+    ( fermatFactorization, pollardPminus1, dixon, dixonDefaultParams
     , gnfsWith, GNFSParams(..), defaultParams
     )
 
@@ -45,13 +45,33 @@ pollardUI = do
 
 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."
     n <- askNumber "Enter an integer n > 1:"
-    b <- askNumber "Enter the bound B:"
-    x <- askNumber "Enter the max random integer x:"
 
+    putStrLn "Parameters:"
+    putStrLn "1) Automatic (recommended)"
+    putStrLn "2) Custom"
+    mode <- askChoice 2
+
+    let (autoB, autoX) = dixonDefaultParams n
+    (b, x) <- case mode of
+        1 -> return (autoB, autoX)
+        _ -> do
+            putStrLn ("Automatic values for this n: bound B = " ++ show autoB
+                      ++ ", max x = " ++ show autoX)
+            b' <- askNumber "Enter the bound B:"
+            x' <- askNumber "Enter the max integer x (the search starts at floor(sqrt n) + 1):"
+            return (b', x')
+
+    putStrLn "Running Dixon (this may take a while)..."
     case dixon n b x of
         Left err -> do
             putStrLn ("Error: " ++ err)
+            if mode == 1
+                then putStrLn "Check that n is composite and not a perfect square or a prime power."
+                else return ()
         Right (p, q) -> do
             putStrLn ("n = " ++ show n)
             putStrLn ("p = " ++ show p)

src/Factorization.hs

@@ -1,4 +1,4 @@
-module Factorization ( fermatFactorization, pollardPminus1, dixon, gnfs, gnfsWith, GNFSParams(..), defaultParams) where
+module Factorization ( fermatFactorization, pollardPminus1, dixon, dixonDefaultParams, gnfs, gnfsWith, GNFSParams(..), defaultParams) where
 
 import Factorization.FermatFactorization
 import Factorization.PollardPminus1

src/Factorization/Dixon.hs

@@ -1,4 +1,4 @@
-module Factorization.Dixon (dixon) where
+module Factorization.Dixon (dixon, dixonDefaultParams) where
 
 import Data.Bits (testBit)
 import Data.List (transpose)
@@ -15,44 +15,73 @@ factorOverBase primes a = go primes a []
     go (p:ps) num acc = let (e, num') = countDiv p num
                         in go ps num' (e : acc)
 
+floorSqrt :: Integer -> Integer
+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
 
+-- | Recommended parameters (b, maxX) for factoring n with 'dixon'.
+--
+--   * b    : smoothness bound (the factor base is all primes <= b).
+--   * maxX : upper end of the search range for x. The squares x^2 mod n are tested for x = floor(sqrt n) + 1 .. maxX. 'dixon' stops scanning as soon as it has enough smooth squares, so maxX is only a ceiling and is chosen generously.
+--
+-- Both values were fitted from measurements on balanced semiprimes of 8 to 24 digits:
+--
+--   * ln b = 0.64 * sqrt(ln n * ln ln n)   (the usual shape of the optimal bound for this family of methods; the running time is flat around the optimum, so the constant need not be exact).
+--   * The range actually needed was between 30 and 930 times the number of primes in the factor base (about b / ln b); the ceiling is 200 * (number of digits of n) times that, i.e. roughly 2 to 4 times the worst case measured.
+--
+-- Dixon is only meant for composite n: on a prime (or a prime power) it will scan the whole range before giving up, which can take a while for large n.
+dixonDefaultParams :: Integer -> (Integer, Integer)
+dixonDefaultParams n = (b, floorSqrt n + width)
+  where
+    ln     = log (fromIntegral (max 16 n)) :: Double
+    b      = max 30 (round (exp (0.64 * sqrt (ln * log ln))))
+    cols   = max 1 (round (fromIntegral b / log (fromIntegral b) :: Double)) :: Integer
+    digits = toInteger (length (show n))
+    width  = 200 * digits * cols
+
+-- | Dixon's factorization method.
+--
+-- Scans x = floor(sqrt n) + 1 .. maxX and keeps the x for which x^2 mod n is smooth over the primes <= b. It stops as soon as it has (number of primes + 20) such squares, which guarantees  at least 20 independent dependencies over GF(2); each one gives a non-trivial factor with probability about 1/2. If all of them fail, the search continues with four times as many squares, until the range up to maxX is exhausted.
 dixon :: Integer -> Integer -> Integer -> Either String (Integer, Integer)
 dixon n b maxX
   | n <= 3    = Left "n must be greater than 3"
   | even n    = Right (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 Left $ "Not enough smooth squares found up to " ++ show maxX
-         else
-           let bitvecs = map (fromBools . map (\e -> e `mod` 2 == 1) . snd) pairs
-               -- a basis of ALL dependencies between the rows (left null space over GF(2))
-               masks = gaussianEliminationAllMasks bitvecs
+  | otherwise = attempt 20
+  where
+    base    = erathosteneSieve b
+    numCols = length base
+    startX  = floorSqrt n + 1
 
-               -- Try one dependency: x = product of the chosen x_i, z = sqrt of the product of the chosen (x_i^2 mod n).
-               -- Each dependency gives a non-trivial factor with probability about 1/2, so we move on to the next one when it fails.
-               tryMask mask =
-                 let chosen   = [ pr | (i, pr) <- zip [0 :: Int ..] pairs, testBit mask i ]
-                     productX = foldl (\acc (x, _) -> modMul acc x n) 1 chosen
-                     fs       = [ sum col `div` 2 | col <- transpose (map snd chosen) ]
-                     z        = foldl (\acc (p, f) -> if f == 0 then acc else modMul (modExp p f n) acc n)
-                                      1 (zip base fs)
-                     g        = gcd (abs (productX - z)) n
-                 in if g == 1 || g == n then Nothing else Just g
-           in case masks of
-                [] -> Left "No linear dependency found"
-                _  -> case [ g | mask <- masks, Just g <- [tryMask mask] ] of
-                        (g : _) -> Right (g, n `div` g)
-                        []      -> Left "non trivial factors not found, try larger b or maximum x"
+    -- lazy list of (x, exponents of x^2 mod n over the base) for the smooth squares
+    smooth = [ (x, exps) | x <- [startX .. maxX]
+                         , let a = (x * x) `mod` n
+                         , a > 1
+                         , Just exps <- [factorOverBase base a] ]
+
+    attempt :: Int -> Either String (Integer, Integer)
+    attempt extra
+      | got <= numCols = Left $ "Not enough smooth squares found up to " ++ show maxX
+      | otherwise =
+          case [ g | mask <- gaussianEliminationAllMasks bitvecs, Just g <- [tryMask mask] ] of
+            (g : _) -> Right (g, n `div` g)
+            []
+              | got < numCols + extra -> Left "non trivial factors not found, try larger b or maximum x"
+              | otherwise             -> attempt (extra * 4)
+      where
+        pairs   = take (numCols + extra) smooth
+        got     = length pairs
+        bitvecs = map (fromBools . map (\e -> e `mod` 2 == 1) . snd) pairs
+
+        -- Try one dependency: x = product of the chosen x_i, z = sqrt of the product of the chosen (x_i^2 mod n).
+        -- Each dependency gives a non-trivial factor with probability about 1/2, so we move on to the next one when it fails.
+        tryMask mask =
+          let chosen   = [ pr | (i, pr) <- zip [0 :: Int ..] pairs, testBit mask i ]
+              productX = foldl (\acc (x, _) -> modMul acc x n) 1 chosen
+              fs       = [ sum col `div` 2 | col <- transpose (map snd chosen) ]
+              z        = foldl (\acc (p, f) -> if f == 0 then acc else modMul (modExp p f n) acc n)
+                               1 (zip base fs)
+              g        = gcd (abs (productX - z)) n
+          in if g == 1 || g == n then Nothing else Just g