diff --git a/app/FactorizationUI.hs b/app/FactorizationUI.hs index 40a1f81..eed7eb3 100644 --- a/app/FactorizationUI.hs +++ b/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) diff --git a/haskell-math.cabal b/haskell-math.cabal index 58e5827..a01e005 100644 --- a/haskell-math.cabal +++ b/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, diff --git a/src/Factorization.hs b/src/Factorization.hs index 8419b0d..bd7b5ef 100644 --- a/src/Factorization.hs +++ b/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 diff --git a/src/Factorization/Dixon.hs b/src/Factorization/Dixon.hs new file mode 100644 index 0000000..886bcab --- /dev/null +++ b/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) diff --git a/src/ModularArithmeticUtils.hs b/src/ModularArithmeticUtils.hs index 5fa0cb8..c2c7af6 100644 --- a/src/ModularArithmeticUtils.hs +++ b/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)