diff --git a/app/FactorizationUI.hs b/app/FactorizationUI.hs index 32db30e..dc4b6fb 100644 --- a/app/FactorizationUI.hs +++ b/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) diff --git a/src/Factorization.hs b/src/Factorization.hs index 7798101..ac51e7e 100644 --- a/src/Factorization.hs +++ b/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 diff --git a/src/Factorization/Dixon.hs b/src/Factorization/Dixon.hs index 0555287..524b671 100644 --- a/src/Factorization/Dixon.hs +++ b/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