diff --git a/app/FactorizationUI.hs b/app/FactorizationUI.hs index f09d15d..32db30e 100644 --- a/app/FactorizationUI.hs +++ b/app/FactorizationUI.hs @@ -1,7 +1,11 @@ module FactorizationUI (run) where +import System.Exit (exitSuccess) import Utils (askChoice, askNumber) -import Factorization (fermatFactorization, pollardPminus1, dixon) +import Factorization + ( fermatFactorization, pollardPminus1, dixon + , gnfsWith, GNFSParams(..), defaultParams + ) run :: IO () run = do @@ -9,13 +13,15 @@ run = do putStrLn "1) Fermat factorization" putStrLn "2) Pollard p-1" putStrLn "3) Dixon" + putStrLn "4) General Number Field Sieve (GNFS)" - choice <- askChoice 3 + choice <- askChoice 4 case choice of 1 -> fermatUI 2 -> pollardUI 3 -> dixonUI + 4 -> gnfsUI _ -> error "Impossible" fermatUI :: IO () @@ -50,3 +56,45 @@ dixonUI = do putStrLn ("n = " ++ show n) putStrLn ("p = " ++ show p) putStrLn ("q = " ++ show q) + +gnfsUI :: IO () +gnfsUI = do + putStrLn "General Number Field Sieve" + putStrLn "Factors a composite n (not a prime or a prime power) into two non-trivial factors." + putStrLn "Practical for n up to roughly 30 digits; the running time grows quickly beyond that." + n <- askNumber "Enter an integer n > 1:" + + putStrLn "Parameters:" + putStrLn "1) Automatic (recommended)" + putStrLn "2) Custom" + mode <- askChoice 2 + + let auto = defaultParams n + params <- case mode of + 1 -> return auto + _ -> do + putStrLn ("Automatic values for this n: degree = " ++ show (paramDegree auto) + ++ ", bound B = " ++ show (paramRatBound auto) + ++ ", sieve half-width = " ++ show (paramHalfWidth auto)) + d <- askNumber "Enter the polynomial degree d (odd, 3 or 5):" + if d /= 3 && d /= 5 + then do + putStrLn "The degree must be 3 or 5." + exitSuccess + else return () + b <- askNumber "Enter the factor base bound B (used for both bases):" + w <- askNumber "Enter the sieve half-width (a ranges over [-w, w]):" + return auto { paramDegree = fromInteger d + , paramRatBound = b + , paramAlgBound = b + , paramHalfWidth = fromInteger w + } + + putStrLn "Running GNFS (this may take a while)..." + case gnfsWith params n of + Left err -> + putStrLn ("Error: " ++ err) + Right (p, q) -> do + putStrLn ("n = " ++ show n) + putStrLn ("p = " ++ show p) + putStrLn ("q = " ++ show q) diff --git a/haskell-math.cabal b/haskell-math.cabal index a01e005..c68c9bb 100644 --- a/haskell-math.cabal +++ b/haskell-math.cabal @@ -6,7 +6,7 @@ version: 0.1.0.0 license: BSD-3-Clause license-file: LICENSE author: Sam HADOW -maintainer: sam.hadow@inbox.lv +maintainer: sam.hadow@proton.me category: Math build-type: Simple extra-doc-files: README.md @@ -32,10 +32,20 @@ library Factorization.FermatFactorization Factorization.PollardPminus1 Factorization.Dixon + Factorization.GNFS ModularSquareRoot.TonelliShanks + GNFS.FactorBase + GNFS.FiniteField + GNFS.PolySelect + GNFS.Polynomial + GNFS.Sieve + GNFS.SquareRoot + GNFS.Util build-depends: base ^>=4.18.2.1, - random ^>=1.2 + random ^>=1.2, + containers ^>=0.6.7, + array ^>=0.5.8 hs-source-dirs: src default-language: GHC2021 diff --git a/src/Factorization.hs b/src/Factorization.hs index bd7b5ef..7798101 100644 --- a/src/Factorization.hs +++ b/src/Factorization.hs @@ -1,5 +1,6 @@ -module Factorization ( fermatFactorization, pollardPminus1, dixon) where +module Factorization ( fermatFactorization, pollardPminus1, dixon, gnfs, gnfsWith, GNFSParams(..), defaultParams) where import Factorization.FermatFactorization import Factorization.PollardPminus1 import Factorization.Dixon +import Factorization.GNFS diff --git a/src/Factorization/Dixon.hs b/src/Factorization/Dixon.hs index a150d9e..0555287 100644 --- a/src/Factorization/Dixon.hs +++ b/src/Factorization/Dixon.hs @@ -1,7 +1,9 @@ module Factorization.Dixon (dixon) where +import Data.Bits (testBit) +import Data.List (transpose) import Primes.ErathosteneSieve (erathosteneSieve) -import LinearAlgebra.GF2 (fromBools, gaussianEliminationMask, maskToIndicesInt) +import LinearAlgebra.GF2 (fromBools, gaussianEliminationAllMasks) import ModularArithmeticUtils (modExp, modMul, countDiv) @@ -9,7 +11,7 @@ factorOverBase :: [Integer] -> Integer -> Maybe [Integer] factorOverBase primes a = go primes a [] where go [] 1 acc = Just (reverse acc) - go [] _ acc = Nothing + go [] _ _ = Nothing go (p:ps) num acc = let (e, num') = countDiv p num in go ps num' (e : acc) @@ -34,20 +36,23 @@ dixon n b maxX numCollected = length pairs in if numCollected <= numCols then Left $ "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 -> Left "No linear dependency found" - Just mask -> - let indices = maskToIndicesInt mask - in if null indices - then Left "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 Left "non trivial factors not found, try larger b or maximum x" - else Right (g, n `div` g) + 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 + + -- 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" diff --git a/src/Factorization/GNFS.hs b/src/Factorization/GNFS.hs new file mode 100644 index 0000000..8de6bc9 --- /dev/null +++ b/src/Factorization/GNFS.hs @@ -0,0 +1,91 @@ +module Factorization.GNFS + ( gnfs + , gnfsWith + , GNFSParams(..) + , defaultParams + ) where + +import Data.Bits (testBit) +import Data.Maybe (listToMaybe, mapMaybe) + +import GNFS.FactorBase +import GNFS.PolySelect +import GNFS.Sieve (sieveRelations) +import GNFS.SquareRoot (extractFactor) +import GNFS.Util (bitLength, isPrimeDet, isqrt, perfectPower) +import LinearAlgebra.GF2 (gaussianEliminationAllMasks) +import Primes.ErathosteneSieve (erathosteneSieve) + +-- | Tunable parameters of the number field sieve. +data GNFSParams = GNFSParams + { paramDegree :: Int -- ^ degree d of the number field polynomial (odd: 3 or 5) + , paramRatBound :: Integer -- ^ rational factor base bound + , paramAlgBound :: Integer -- ^ algebraic factor base bound + , paramQCCount :: Int -- ^ number of quadratic-character columns + , paramHalfWidth :: Int -- ^ sieve a over [-w, w] for each b + , paramMaxB :: Integer -- ^ give up after sieving b = 1 .. maxB + , paramSlack :: Double -- ^ log-sieve tolerance (natural-log units) + , paramExtra :: Int -- ^ relations collected beyond the number of columns + } deriving Show + +-- | Heuristic parameters from the size of n. The factor base bound was fitted empirically (a bit more than the asymptotic L-notation optimum, which undershoots badly at these small sizes): ln B grows by about 0.058 per bit of n. +defaultParams :: Integer -> GNFSParams +defaultParams n = GNFSParams + { paramDegree = if bits <= 130 then 3 else 5 + , paramRatBound = bound + , paramAlgBound = bound + , paramQCCount = min 64 (max 16 (bits `div` 3)) + , paramHalfWidth = fromInteger (bound * 3 `div` 2) + , paramMaxB = 200000 + , paramSlack = 1.5 + , paramExtra = 40 + } + where + bits = bitLength n + bound = max 300 (round (exp (6.9 + 0.058 * fromIntegral (bits - 40) :: Double))) + +-- | Factor n with the General Number Field Sieve using heuristic parameters. +-- Returns a non-trivial factorisation n = p * q (p, q not necessarily prime). +gnfs :: Integer -> Either String (Integer, Integer) +gnfs n = gnfsWith (defaultParams n) n + +gnfsWith :: GNFSParams -> Integer -> Either String (Integer, Integer) +gnfsWith params n + | n < 4 = Left "n must be at least 4" + | even n = Right (2, n `div` 2) + | Just p <- smallFactor = Right (p, n `div` p) + | isPrimeDet n = Left "n is prime" + | Just (r, _) <- perfectPower n = Right (r, n `div` r) + | n < 10 ^ (12 :: Int) = trialDivision + | otherwise = + case selectPolynomial (paramDegree params) n of + Nothing -> Left "no suitable polynomial found" + Just np -> run np + where + -- tiny inputs: GNFS parameters degenerate, plain trial division is instant + trialDivision = + case [ d | d <- [3, 5 .. isqrt n], n `mod` d == 0 ] of + (d : _) -> Right (d, n `div` d) + [] -> Left "n is prime" + + smallFactor = listToMaybe [ p | p <- erathosteneSieve 1000, n `mod` p == 0, p < n ] + + run np = + let fb = buildFactorBase np (paramRatBound params) (paramAlgBound params) (paramQCCount params) + ncols = fbNumCols fb + stream = sieveRelations fb (paramHalfWidth params) (paramSlack params) (paramMaxB params) + attempt extra + | length rels < need = Left "sieve space exhausted before enough relations were found" + | otherwise = + case listToMaybe (mapMaybe tryMask masks) of + Just g -> Right (g, n `div` g) + Nothing + | extra >= 8 * paramExtra params -> Left "no dependency produced a factor" + | otherwise -> attempt (extra * 2) + where + need = ncols + extra + rels = take need stream + rows = map (relationRow fb) rels + masks = gaussianEliminationAllMasks rows + tryMask mk = extractFactor n np [ r | (i, r) <- zip [0 :: Int ..] rels, testBit mk i ] + in attempt (paramExtra params) diff --git a/src/GNFS/FactorBase.hs b/src/GNFS/FactorBase.hs new file mode 100644 index 0000000..0b0a981 --- /dev/null +++ b/src/GNFS/FactorBase.hs @@ -0,0 +1,114 @@ +module GNFS.FactorBase + ( Relation(..) + , FactorBase(..) + , buildFactorBase + , relationRow + ) where + +import Data.Bits (setBit) +import Data.List (foldl') +import qualified Data.Map.Strict as M + +import GNFS.FiniteField (rootsModP) +import GNFS.PolySelect (NFSPoly(..)) +import GNFS.Polynomial +import GNFS.Util (primesFrom, modInv) +import LinearAlgebra.GF2 (BitVec) +import ModularArithmeticUtils (legendre) +import Primes.ErathosteneSieve (erathosteneSieve) + +-- | A relation: the pair (a, b) with a - b*m smooth over the rational base and +-- a - b*alpha smooth over the algebraic base, with the full exponent data. +data Relation = Relation + { relA :: !Integer + , relB :: !Integer + , relRat :: [(Integer, Int)] -- prime -> exponent of |a - b m| + , relAlg :: [((Integer, Integer), Int)] -- (p, r) -> valuation of a - b alpha + } deriving Show + +data FactorBase = FactorBase + { fbPoly :: NFSPoly + , fbRat :: [Integer] -- rational base: primes up to the bound + , fbAlg :: [(Integer, Integer)] -- algebraic base: (p, r), f(r) = 0 mod p, simple root + , fbAlgPrimes :: [Integer] -- distinct primes occurring in fbAlg + , fbAlgLifts :: [(Integer, [(Integer, Integer)])] + -- per prime p: for each root r, the Hensel lifts + -- (p^k, r_k) for k = 1, 2, ... (lazy, infinite) + , fbReal :: [Double] -- real roots of f (sign columns) + , fbQC :: [(Integer, Integer)] -- quadratic-character primes (q, s) + , fbColRat :: M.Map Integer Int + , fbColAlg :: M.Map (Integer, Integer) Int + , fbNumCols :: Int + } + +-- | Build the three factor bases (rational, algebraic, quadratic characters). +buildFactorBase :: NFSPoly -> Integer -> Integer -> Int -> FactorBase +buildFactorBase np ratBound algBound nQC = fb + where + f = npCoeffs np + f' = derivP f + rat = erathosteneSieve ratBound + -- primes whose roots are all simple; primes with a repeated root are dropped + algPrimeRoots = + [ (p, rs) + | p <- erathosteneSieve algBound + , let rs = rootsModP p f + , not (null rs) + , all (\r -> evalP f' r `mod` p /= 0) rs + ] + alg = [ (p, r) | (p, rs) <- algPrimeRoots, r <- rs ] + qc = take nQC + [ (q, s) + | q <- primesFrom (max algBound ratBound + 1) + , let rs = rootsModP q f + , s <- rs + , evalP f' s `mod` q /= 0 + ] + reals = realRoots f + -- Hensel lifting of a simple root r of f mod p to roots mod p^k + lifts p r0 = go p r0 + where + go q rk = + let q' = q * p + inv = case modInv (evalP f' rk) q' of + Just i -> i + Nothing -> error "Hensel lift: derivative not invertible" + rk' = (rk - evalP f rk * inv) `mod` q' + in (q, rk) : go q' rk' + nRat = length rat + nAlg = length alg + fb = FactorBase + { fbPoly = np + , fbRat = rat + , fbAlg = alg + , fbAlgPrimes = map fst algPrimeRoots + , fbAlgLifts = [ (p, lifts p r) | (p, r) <- alg ] + , fbReal = reals + , fbQC = qc + , fbColRat = M.fromList (zip rat [0 ..]) + , fbColAlg = M.fromList (zip alg [0 ..]) + , fbNumCols = 1 + nRat + nAlg + length reals + length qc + } + +-- | Parity vector of a relation. Column layout: +-- [0] sign of a - b m | rational primes | algebraic ideals | real-place signs | quadratic characters +relationRow :: FactorBase -> Relation -> BitVec +relationRow fb rel = foldl' setBit 0 cols + where + a = relA rel + b = relB rel + m = npM (fbPoly fb) + nRat = length (fbRat fb) + nAlg = length (fbAlg fb) + offAlg = 1 + nRat + offReal = offAlg + nAlg + offQC = offReal + length (fbReal fb) + signCol = [ 0 | a - b * m < 0 ] + ratCols = [ 1 + fbColRat fb M.! p | (p, e) <- relRat rel, odd e ] + algCols = [ offAlg + fbColAlg fb M.! pr | (pr, e) <- relAlg rel, odd e ] + realCols = [ offReal + j | (j, th) <- zip [0 ..] (fbReal fb) + , fromInteger a < fromInteger b * th ] + qcCols = [ offQC + j | (j, (q, s)) <- zip [0 ..] (fbQC fb) + , let v = (a - b * s) `mod` q + , legendre v q == q - 1 ] + cols = signCol ++ ratCols ++ algCols ++ realCols ++ qcCols diff --git a/src/GNFS/FiniteField.hs b/src/GNFS/FiniteField.hs new file mode 100644 index 0000000..a2faa6d --- /dev/null +++ b/src/GNFS/FiniteField.hs @@ -0,0 +1,79 @@ +module GNFS.FiniteField + ( irreducibleModP + , rootsModP + , fieldSqrt + , fieldNorm + ) where + +import GNFS.Polynomial +import GNFS.Util (primeDivisors) +import ModularArithmeticUtils (factorOutTwos) + +-- | Rabin's irreducibility test for f over F_p. +irreducibleModP :: Integer -> Poly -> Bool +irreducibleModP p f0 + | d < 1 = False + | otherwise = frobPow d == x && all coprime (primeDivisors d) + where + f = modCoeffs p f0 + d = degP f + x = remP p f [0, 1] + frobPow k = iterate (\u -> powMod p f u p) x !! k + coprime q = gcdP p (subP (frobPow (d `div` q)) x) f == [1] + +-- | All roots of f in F_p. +rootsModP :: Integer -> Poly -> [Integer] +rootsModP p f0 + | p < 100 = [ r | r <- [0 .. p - 1], evalP f0 r `mod` p == 0 ] + | otherwise = split (gcdP p (subP (powMod p f [0, 1] p) [0, 1]) f) + where + f = modCoeffs p f0 + -- g is a product of distinct linear factors + split g + | degP g <= 0 = [] + | degP g == 1 = [ negate (head g) `mod` p ] + | otherwise = case [ (h, q) | c <- [0 .. p - 1] + , let t = powMod p g [c, 1] ((p - 1) `div` 2) + , let h = gcdP p (subP t [1]) g + , degP h > 0, degP h < degP g + , let (q, _) = divModP p h g ] of + ((h, q) : _) -> split h ++ split q + [] -> error "rootsModP: splitting failed" + +-- | Square root in F_p[x]/(f) (f irreducible of degree d) by Tonelli-Shanks. +-- Returns Nothing if the element is not a square. +fieldSqrt :: Integer -> Poly -> Poly -> Maybe Poly +fieldSqrt p f g0 + | null g = Just [] + | powMod p f g halfQ /= [1] = Nothing + | otherwise = Just (go s c0 t0 r0) + where + g = remP p f (modCoeffs p g0) + d = length f - 1 + q = p ^ d + halfQ = (q - 1) `div` 2 + (qq, s) = factorOutTwos (q - 1) + mul = mulMod p f + sq u = mul u u + pw = powMod p f + -- a quadratic non-residue + z = head [ cand | cand <- map (\c -> [c]) [2 .. p - 1] ++ [ [c0, 1] | c0 <- [0 .. p - 1] ] + , pw cand halfQ == [p - 1] ] + c0 = pw z qq + t0 = pw g qq + r0 = pw g ((qq + 1) `div` 2) + go m c t r + | t == [1] = r + | otherwise = + let i = length (takeWhile (/= [1]) (iterate sq t)) -- least i with t^(2^i) = 1 + b = iterate sq c !! fromInteger (m - toInteger i - 1) + c' = sq b + in go (toInteger i) c' (mul t c') (mul r b) + +-- | Norm of u from F_p[x]/(f) down to F_p, as an element of [0, p). +fieldNorm :: Integer -> Poly -> Poly -> Integer +fieldNorm p f u = + case powMod p f u ((p ^ d - 1) `div` (p - 1)) of + [] -> 0 + (c:_) -> c `mod` p + where d = length f - 1 diff --git a/src/GNFS/PolySelect.hs b/src/GNFS/PolySelect.hs new file mode 100644 index 0000000..de352b7 --- /dev/null +++ b/src/GNFS/PolySelect.hs @@ -0,0 +1,47 @@ +module GNFS.PolySelect + ( NFSPoly(..) + , selectPolynomial + ) where + +import Data.Maybe (listToMaybe) +import GNFS.FiniteField (irreducibleModP) +import GNFS.Polynomial +import GNFS.Util (iroot) + +-- | A monic polynomial f of degree d and an integer m with f(m) = n. +data NFSPoly = NFSPoly + { npDegree :: Int + , npCoeffs :: Poly -- ^ little-endian, monic + , npM :: Integer + } deriving Show + +-- | Base-m polynomial selection. Returns a monic polynomial of the requested degree, irreducible over Q (certified by irreducibility modulo some small prime), with f(m) = n. +-- Balanced digits are tried first (smaller coefficients => smaller norms). +selectPolynomial :: Int -> Integer -> Maybe NFSPoly +selectPolynomial d n = listToMaybe + [ NFSPoly d f m + | m <- takeWhile (\m' -> m' > 1 && 2 * m' ^ d > n) [m0, m0 - 1 ..] + , digits <- [balanced m (plainDigits m n), plainDigits m n] + , let f = trimP digits + , length f == d + 1 + , last f == 1 + , irreducibleSomewhere f + ] + where + m0 = iroot d n + smallPrimes = [3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73] + irreducibleSomewhere f = any (\p -> irreducibleModP p f) smallPrimes + +plainDigits :: Integer -> Integer -> [Integer] +plainDigits _ 0 = [] +plainDigits m x = let (q, r) = x `divMod` m in r : plainDigits m q + +-- | Convert digits in [0, m) to balanced digits in (-m/2, m/2], propagating carries. +balanced :: Integer -> [Integer] -> [Integer] +balanced m = go 0 + where + go 0 [] = [] + go c [] = [c] + go c (x:xs) = + let v = x + c + in if 2 * v > m then (v - m) : go 1 xs else v : go 0 xs diff --git a/src/GNFS/Polynomial.hs b/src/GNFS/Polynomial.hs new file mode 100644 index 0000000..6959e4c --- /dev/null +++ b/src/GNFS/Polynomial.hs @@ -0,0 +1,158 @@ +module GNFS.Polynomial + ( Poly + , trimP + , degP + , addP, subP, negP, scaleP, mulP + , evalP + , evalHomog + , derivP + , reduceMonic + , modCoeffs + , remP, divModP + , gcdP + , mulMod + , powMod + , realRoots + ) where + +import GNFS.Util (modInv) + +-- | Polynomial over Z (or F_p), coefficients little-endian: [a0, a1, a2, ...]. +-- All functions return a trimmed list (no trailing zeros); the zero polynomial is []. +type Poly = [Integer] + +trimP :: Poly -> Poly +trimP = reverse . dropWhile (== 0) . reverse + +-- | Degree (-1 for the zero polynomial). +degP :: Poly -> Int +degP p = length (trimP p) - 1 + +addP :: Poly -> Poly -> Poly +addP a b = trimP (go a b) + where + go (x:xs) (y:ys) = (x + y) : go xs ys + go xs [] = xs + go [] ys = ys + +negP :: Poly -> Poly +negP = map negate + +subP :: Poly -> Poly -> Poly +subP a b = addP a (negP b) + +scaleP :: Integer -> Poly -> Poly +scaleP c = trimP . map (* c) + +mulP :: Poly -> Poly -> Poly +mulP [] _ = [] +mulP _ [] = [] +mulP (x:xs) ys = addP (scaleP x ys) (0 : mulP xs ys) + +evalP :: Poly -> Integer -> Integer +evalP f x = foldr (\c acc -> c + x * acc) 0 f + +-- | Homogeneous evaluation b^d * f(a/b) = sum c_i a^i b^(d-i), with d = deg f. +evalHomog :: Poly -> Integer -> Integer -> Integer +evalHomog f a b = sum [ c * a ^ i * b ^ (d - i) | (i, c) <- zip [0 :: Int ..] f ] + where d = length f - 1 + +derivP :: Poly -> Poly +derivP [] = [] +derivP (_:cs) = trimP (zipWith (*) [1 ..] cs) + +-- | Reduce u modulo a MONIC polynomial f over Z (exact, no division needed). +reduceMonic :: Poly -> Poly -> Poly +reduceMonic f u = go (trimP u) + where + d = length f - 1 + go v + | length v <= d = v + | otherwise = + let c = last v + k = length v - 1 - d + in go (trimP (subP v (scaleP c (replicate k 0 ++ f)))) + +-- --- +-- Arithmetic over F_p[x] + +modCoeffs :: Integer -> Poly -> Poly +modCoeffs p = trimP . map (`mod` p) + +-- | Quotient and remainder of u by g (g /= 0) over F_p. +divModP :: Integer -> Poly -> Poly -> (Poly, Poly) +divModP p g0 u0 = go (modCoeffs p u0) [] + where + g = modCoeffs p g0 + dg = length g - 1 + lcI = case modInv (last g) p of + Just i -> i + Nothing -> error "divModP: non-invertible leading coefficient" + go v q + | length v <= dg = (q, v) + | otherwise = + let c = (last v * lcI) `mod` p + k = length v - 1 - dg + v' = modCoeffs p (subP v (scaleP c (replicate k 0 ++ g))) + q' = addP q (replicate k 0 ++ [c]) + in go v' q' + +remP :: Integer -> Poly -> Poly -> Poly +remP p g u = snd (divModP p g u) + +-- | Monic gcd over F_p. +gcdP :: Integer -> Poly -> Poly -> Poly +gcdP p a0 b0 = monic (go (modCoeffs p a0) (modCoeffs p b0)) + where + go a [] = a + go a b = go b (remP p b a) + monic [] = [] + monic g = case modInv (last g) p of + Just i -> modCoeffs p (scaleP i g) + Nothing -> g + +-- | Multiplication in F_p[x]/(f). +mulMod :: Integer -> Poly -> Poly -> Poly -> Poly +mulMod p f u v = remP p f (mulP u v) + +-- | Exponentiation in F_p[x]/(f) (e >= 0). +powMod :: Integer -> Poly -> Poly -> Integer -> Poly +powMod p f u0 e0 = go (remP p f u0) e0 [1] + where + go _ 0 acc = remP p f acc + go b e acc = + let acc' = if odd e then mulMod p f acc b else acc + in go (mulMod p f b b) (e `div` 2) acc' + +-- --- +-- Real roots (floating point; used only for sign bookkeeping) + +-- | All real roots of a polynomial with simple roots, as Doubles, by recursing on the derivative's roots (Rolle) and bisecting between them. +realRoots :: Poly -> [Double] +realRoots f0 + | n <= 0 = [] + | n == 1 = [ negate (head fd) / (fd !! 1) ] + | otherwise = concat (zipWith bisect pts (tail pts)) + where + f = trimP f0 + n = length f - 1 + fd = map fromInteger f :: [Double] + crit = realRoots (derivP f) + bound = 1 + maximum (map abs (init fd)) / abs (last fd) + pts = [negate bound] ++ crit ++ [bound] + ev x = foldr (\c acc -> c + x * acc) 0 fd + bisect lo hi + | slo == 0 = [] -- (cannot happen for simple roots; avoid duplicates) + | shi == 0 = [hi] + | slo /= shi = [go (200 :: Int) lo hi] + | otherwise = [] + where + slo = signum (ev lo) + shi = signum (ev hi) + go 0 l h = (l + h) / 2 + go k l h = + let mid = (l + h) / 2 + sm = signum (ev mid) + in if sm == 0 then mid + else if sm == slo then go (k - 1) mid h + else go (k - 1) l mid diff --git a/src/GNFS/Sieve.hs b/src/GNFS/Sieve.hs new file mode 100644 index 0000000..27fd8a2 --- /dev/null +++ b/src/GNFS/Sieve.hs @@ -0,0 +1,124 @@ +{-# LANGUAGE BangPatterns #-} +module GNFS.Sieve + ( sieveRelations + , sieveLine + ) where + +import Control.Monad (forM_, unless) +import Control.Monad.ST (ST) +import Data.Array.Base (unsafeAt, unsafeRead, unsafeWrite) +import Data.Array.ST (STUArray, newArray, runSTUArray) +import Data.Array.Unboxed (UArray) +import qualified Data.Map.Strict as M + +import GNFS.FactorBase +import GNFS.PolySelect (NFSPoly(..)) +import GNFS.Polynomial (evalHomog) +import GNFS.Util (modInv) +import ModularArithmeticUtils (countDiv) + +-- | Lazy stream of relations: lines b = 1 .. maxB, with a in [-halfA, halfA]. +sieveRelations :: FactorBase -> Int -> Double -> Integer -> [Relation] +sieveRelations fb halfA slack maxB = concatMap (sieveLine fb halfA slack) [1 .. maxB] + +-- | Line-sieve a fixed b over a in [-halfA, halfA], return the verified relations. +-- Log-sums over the factor-base primes (not prime powers) are used as a filter; every survivor is then fully trial-factored, so 'slack' only trades speed against yield. +sieveLine :: FactorBase -> Int -> Double -> Integer -> [Relation] +sieveLine fb halfA slack b = concatMap check [0 .. 2 * halfA] + where + np = fbPoly fb + f = npCoeffs np + d = npDegree np + m = npM np + hi = 2 * halfA + bD = fromInteger b :: Double + bmD = fromInteger (b * m) :: Double + fD = map fromInteger f :: [Double] + + -- prime powers up to this limit are sieved (larger ones hit at most one position) + powLimit = max 1000 (4 * toInteger halfA) + + -- (modulus as Int, log p, start index) for the rational side: a = b*m (mod p^k) + ratSpecs = + [ (qI, log (fromInteger p :: Double), start) + | p <- fbRat fb + , b `mod` p /= 0 + , q <- takeWhile (<= powLimit) (iterate (* p) p) + , let qI = fromInteger q :: Int + , let r = fromInteger ((b `mod` q) * (m `mod` q) `mod` q) :: Int + , let start = (r + halfA) `mod` qI + ] + -- same for the algebraic side: a = b*r_k (mod p^k) for the Hensel-lifted roots + algSpecs = + [ (qI, log (fromInteger p :: Double), start) + | (p, lifts) <- fbAlgLifts fb + , b `mod` p /= 0 + , (q, rk) <- takeWhile ((<= powLimit) . fst) lifts + , let qI = fromInteger q :: Int + , let r = fromInteger ((b `mod` q) * rk `mod` q) :: Int + , let start = (r + halfA) `mod` qI + ] + + runSieve :: [(Int, Double, Int)] -> UArray Int Double + runSieve specs = runSTUArray $ do + arr <- newArray (0, hi) 0 :: ST s (STUArray s Int Double) + forM_ specs $ \(p, lp, start) -> do + let go !i = unless (i > hi) $ do + v <- unsafeRead arr i + unsafeWrite arr i (v + lp) + go (i + p) + go start + return arr + + ratArr = runSieve ratSpecs + algArr = runSieve algSpecs + + check :: Int -> [Relation] + check i = + let a = i - halfA + aD = fromIntegral a :: Double + nr = abs (aD - bmD) + in if nr > 0 && unsafeAt ratArr i >= log nr - slack + then + let fv = sum [ c * aD ^ k * bD ^ (d - k) | (k, c) <- zip [0 :: Int ..] fD ] + in if fv /= 0 && unsafeAt algArr i >= log (abs fv) - slack + then verify (toInteger a) + else [] + else [] + + verify :: Integer -> [Relation] + verify a + | gcd a b /= 1 = [] + | nRat == 0 = [] + | otherwise = + case (factorRat, factorAlg) of + (Just re, Just ae) -> [Relation a b re ae] + _ -> [] + where + nRat = a - b * m + fVal = evalHomog f a b + factorRat = trialFactor (fbRat fb) (abs nRat) + factorAlg = algFactor (fbAlgPrimes fb) (abs fVal) [] + + algFactor [] rest acc = if rest == 1 then Just (reverse acc) else Nothing + algFactor (p:ps) rest acc = + let (e, rest') = countDiv p rest + in if e == 0 + then algFactor ps rest' acc + else case modInv b p of + Nothing -> Nothing + Just bi -> + let r = (a * bi) `mod` p + in if M.member (p, r) (fbColAlg fb) + then algFactor ps rest' (((p, r), fromInteger e) : acc) + else Nothing + +-- | Trial division over a list of primes; Nothing unless n factors completely. +trialFactor :: [Integer] -> Integer -> Maybe [(Integer, Int)] +trialFactor primes n0 = go primes n0 [] + where + go _ 1 acc = Just (reverse acc) + go [] _ _ = Nothing + go (p:ps) n acc = + let (e, n') = countDiv p n + in if e == 0 then go ps n acc else go ps n' ((p, fromInteger e) : acc) diff --git a/src/GNFS/SquareRoot.hs b/src/GNFS/SquareRoot.hs new file mode 100644 index 0000000..235ea86 --- /dev/null +++ b/src/GNFS/SquareRoot.hs @@ -0,0 +1,99 @@ +module GNFS.SquareRoot + ( algebraicSqrt + , extractFactor + , normElement + ) where + +import Control.Monad (guard) +import Data.List (foldl') +import Data.Maybe (listToMaybe) +import qualified Data.Map.Strict as M + +import GNFS.FactorBase (Relation(..)) +import GNFS.FiniteField (irreducibleModP, fieldSqrt, fieldNorm) +import GNFS.PolySelect (NFSPoly(..)) +import GNFS.Polynomial +import GNFS.Util (bitLength, isqrt, modInv, primesFrom) +import ModularArithmeticUtils (modExp) + +-- | Norm N(u(alpha)) of an element of Z[alpha] = Z[x]/(f): the determinant of the multiplication-by-u matrix. +normElement :: Poly -> Poly -> Integer +normElement f u = detInt [ pad (reduceMonic f (mulP u (replicate i 0 ++ [1]))) | i <- [0 .. d - 1] ] + where + d = length f - 1 + pad v = take d (v ++ repeat 0) + +detInt :: [[Integer]] -> Integer +detInt [] = 1 +detInt (r : rs) = sum [ sgn j * c * detInt (map (dropAt j) rs) | (j, c) <- zip [0 :: Int ..] r, c /= 0 ] + where + sgn j = if even j then 1 else -1 + dropAt j row = take j row ++ drop (j + 1) row + +forceP :: Poly -> Poly +forceP p = foldr seq () p `seq` p + +-- | Couveignes' algebraic square root. +-- +-- Given a set of relations whose product prod (a - b alpha) is (hopefully) a square in the number field, compute beta in Z[alpha] with +-- beta^2 = f'(alpha)^2 * prod (a - b alpha). +-- The square root is computed modulo many primes p for which f is irreducible (so Z[alpha]/p is the field F_{p^d}); since d is odd the two roots +-r have norms of opposite sign, and we fix the sign by matching the (known) norm of beta. Candidates are lifted by CRT and verified exactly. +algebraicSqrt :: NFSPoly -> [Relation] -> Maybe Poly +algebraicSqrt np rels + | prodF <= 0 || s * s /= prodF = Nothing + | otherwise = loop 1 (replicate d 0) (primesFrom (2 ^ (50 :: Int))) + where + f = npCoeffs np + d = npDegree np + fp = derivP f + fp2 = reduceMonic f (mulP fp fp) + + prodLin = foldl' (\acc r -> forceP (reduceMonic f (mulP acc [relA r, negate (relB r)]))) [1] rels + gamma2 = reduceMonic f (mulP fp2 prodLin) + + prodF = product [ evalHomog f (relA r) (relB r) | r <- rels ] + s = isqrt prodF + target = normElement f fp * s -- N(beta), up to the global sign of beta + + maxCoef = maximum (1 : map abs gamma2) + cap = bitLength maxCoef + d * (bitLength (maximum (map abs f)) + 8) + 128 + + loop :: Integer -> [Integer] -> [Integer] -> Maybe Poly + loop _ _ [] = Nothing + loop modulus cs (p : ps) + | bitLength modulus > cap = Nothing + | not (irreducibleModP p f) = loop modulus cs ps + | t == 0 = loop modulus cs ps + | otherwise = + case fieldSqrt p f (modCoeffs p gamma2) of + Nothing -> Nothing + Just r0 -> + let nr = fieldNorm p f r0 + in if nr == t then step p r0 + else if nr == (p - t) `mod` p then step p (modCoeffs p (negP r0)) + else Nothing + where + t = target `mod` p + step q r = + let rs = take d (r ++ repeat 0) + invM = case modInv modulus q of { Just i -> i; Nothing -> 1 } + cs' = zipWith (\c x -> c + modulus * (((x - c) * invM) `mod` q)) cs rs + modulus' = modulus * q + lifted = trimP [ if 2 * c > modulus' then c - modulus' else c | c <- cs' ] + in if reduceMonic f (mulP lifted lifted) == gamma2 + then Just lifted + else loop modulus' cs' ps + +-- | From a set of relations that form a dependency, try to produce a factor of n. +extractFactor :: Integer -> NFSPoly -> [Relation] -> Maybe Integer +extractFactor n np rels = do + let f = npCoeffs np + m = npM np + expMap = M.fromListWith (+) [ (p, e) | r <- rels, (p, e) <- relRat r ] + guard (all even (M.elems expMap)) + let y = foldl' (\acc (p, e) -> acc * modExp p (toInteger (e `div` 2)) n `mod` n) 1 (M.toList expMap) + beta <- algebraicSqrt np rels + let x = evalP beta m `mod` n -- phi(beta) + z = evalP (derivP f) m * y `mod` n -- f'(m) * y + guard ((x * x - z * z) `mod` n == 0) + listToMaybe [ g | g <- [gcd (x - z) n, gcd (x + z) n], g > 1, g < n ] diff --git a/src/GNFS/Util.hs b/src/GNFS/Util.hs new file mode 100644 index 0000000..24cebdf --- /dev/null +++ b/src/GNFS/Util.hs @@ -0,0 +1,106 @@ +{-# LANGUAGE BangPatterns #-} +module GNFS.Util + ( bitLength + , isqrt + , iroot + , perfectPower + , isPrimeDet + , primesFrom + , modInv + , primeDivisors + ) where + +import Data.Bits (shiftR, countLeadingZeros) +import Data.Word (Word64) +import ModularArithmeticUtils (modExp, factorOutTwos) + +-- | Number of bits of a positive Integer (0 for n <= 0). +bitLength :: Integer -> Int +bitLength n + | n <= 0 = 0 + | otherwise = go 0 n + where + go :: Int -> Integer -> Int + go !acc x + | x >= 18446744073709551616 = go (acc + 64) (x `shiftR` 64) + | otherwise = acc + (64 - countLeadingZeros (fromInteger x :: Word64)) + +-- | floor(sqrt n) +isqrt :: Integer -> Integer +isqrt n + | n < 0 = error "isqrt: negative" + | n < 2 = n + | otherwise = go (2 ^ ((bitLength n + 1) `div` 2)) + where + go x = let y = (x + n `div` x) `div` 2 + in if y >= x then x else go y + +-- | floor(n^(1/k)) +iroot :: Int -> Integer -> Integer +iroot k n + | n < 2 = n + | k <= 1 = n + | otherwise = go (2 ^ ((bitLength n + k - 1) `div` k)) + where + k' = toInteger k + go x = let y = ((k' - 1) * x + n `div` (x ^ (k - 1))) `div` k' + in if y >= x then x else go y + +-- | If n = r^k for some k >= 2, return (r, k). +perfectPower :: Integer -> Maybe (Integer, Int) +perfectPower n = go 2 + where + maxK = bitLength n + go k + | k > maxK = Nothing + | otherwise = + let r = iroot k n + in if r > 1 && r ^ k == n then Just (r, k) else go (k + 1) + +-- | Strong probable-prime test with the first 13 prime bases. +-- Deterministic for n < 3.3e24, and an overwhelmingly reliable test above. +isPrimeDet :: Integer -> Bool +isPrimeDet n + | n < 2 = False + | n `elem` smallPs = True + | any (\p -> n `mod` p == 0) smallPs = False + | otherwise = all witness smallPs + where + smallPs = [2,3,5,7,11,13,17,19,23,29,31,37,41] + (d, s) = factorOutTwos (n - 1) + witness a = + let x = modExp a d n + in x == 1 || x == n - 1 + || any (== n - 1) (take (fromInteger s - 1) (tail (iterate (\y -> y * y `mod` n) x))) + +-- | Increasing list of primes >= m. +primesFrom :: Integer -> [Integer] +primesFrom m + | m <= 2 = 2 : primesFrom 3 + | even m = filter isPrimeDet [m + 1, m + 3 ..] + | otherwise = filter isPrimeDet [m, m + 2 ..] + +-- | Modular inverse, if it exists. +modInv :: Integer -> Integer -> Maybe Integer +modInv a m = + let (g, x, _) = egcd (a `mod` m) m + in if g == 1 then Just (x `mod` m) else Nothing + +egcd :: Integer -> Integer -> (Integer, Integer, Integer) +egcd a 0 = (a, 1, 0) +egcd a b = + let (q, r) = a `divMod` b + (g, x, y) = egcd b r + in (g, y, x - q * y) + +-- | Distinct prime divisors of a small positive Int. +primeDivisors :: Int -> [Int] +primeDivisors = go 2 + where + go p k + | k < 2 = [] + | p * p > k = [k] + | k `mod` p == 0 = p : go (p + 1) (strip p k) + | otherwise = go (p + 1) k + strip p k | k `mod` p == 0 = strip p (k `div` p) + | otherwise = k diff --git a/src/LinearAlgebra/GF2.hs b/src/LinearAlgebra/GF2.hs index 398da80..611b9aa 100644 --- a/src/LinearAlgebra/GF2.hs +++ b/src/LinearAlgebra/GF2.hs @@ -6,14 +6,14 @@ module LinearAlgebra.GF2 , fromBools , toBools -- Gaussian elimination - , gaussianEliminationMask -- :: nCols -> [BitVec] -> Maybe Mask - , gaussianEliminationIndices -- :: nCols -> [BitVec] -> Maybe [Int] + , gaussianEliminationAllMasks -- :: [BitVec] -> [Mask] , maskToIndicesInt -- :: Mask -> [Int] ) where import Data.Bits import Data.List (findIndex, foldl') import Data.Maybe (listToMaybe, fromMaybe) +import qualified Data.Map.Strict as M -- vector over GF(2) represented as an Integer bitmask. The bit i corresponding to the column i. type BitVec = Integer @@ -33,47 +33,25 @@ toBools :: Int -> BitVec -> [Bool] toBools nCols v = [ testBit v i | i <- [0 .. nCols - 1] ] -- Gaussian elimination over GF(2) using Integers --- Gaussian elimination on a list of BitVec rows (each row is an Integer bitmask of length nCols) and return a mask (as Integer) whose set bits indicate which input rows XOR to produce a vector full of 0 (or Nothing if not enough rows to find a combination). -gaussianEliminationMask :: Int -> [BitVec] -> Maybe Mask -gaussianEliminationMask nCols rows - | null rows = Nothing - | otherwise = findZeroMask (eliminate 0 augmented) +-- Gaussian elimination on a list of BitVec rows (each row is an Integer bitmask of length nCols) and return a mask (as Integer) whose set bits indicate which input rows XOR to produce a vector full of 0. +gaussianEliminationAllMasks :: [BitVec] -> [Mask] +gaussianEliminationAllMasks rows = go M.empty (zip rows (map bit [0 ..])) [] where - m :: Int - m = length rows + go _ [] acc = reverse acc + go basis ((v, mk) : rest) acc = + case reduce basis v mk of + (0, mk') -> go basis rest (mk' : acc) + (v', mk') -> go (M.insert (lowBit v') (v', mk') basis) rest acc - -- initial mask = 1 << i - augmented :: [(BitVec, Mask)] - augmented = zip rows (map (\i -> bit i) [0..m-1]) + reduce basis v mk + | v == 0 = (0, mk) + | otherwise = case M.lookup (lowBit v) basis of + Nothing -> (v, mk) + Just (bv, bm) -> reduce basis (v `xor` bv) (mk `xor` bm) - -- Top level elimination - eliminate :: Int -> [(BitVec, Mask)] -> [(BitVec, Mask)] - eliminate col mat - | col >= nCols = mat - | otherwise = - case findPivot col mat of - Nothing -> eliminate (col + 1) mat - Just pivotIdx -> - let (pivotRow, rest) = removeAt pivotIdx mat - (pvVec, pvMask) = pivotRow - rest' = map (xorIfHasBit pvVec pvMask col) rest - in eliminate (col + 1) (pivotRow : rest') + lowBit :: Integer -> Integer + lowBit v = v .&. negate v - findPivot :: Int -> [(BitVec, Mask)] -> Maybe Int - findPivot col mat = findIndex (\(v,_) -> testBit v col) mat - - xorIfHasBit :: BitVec -> Mask -> Int -> (BitVec, Mask) -> (BitVec, Mask) - xorIfHasBit pvVec pvMask col (v, mask) - | testBit v col = (v `xor` pvVec, mask `xor` pvMask) - | otherwise = (v, mask) - - findZeroMask :: [(BitVec, Mask)] -> Maybe Mask - findZeroMask mat = fmap snd (listToMaybe (filter (\(v,_) -> v == 0) mat)) - - --- wrapper -gaussianEliminationIndices :: Int -> [BitVec] -> Maybe [Int] -gaussianEliminationIndices nCols rows = fmap maskToIndicesInt (gaussianEliminationMask nCols rows) -- helpers diff --git a/test/GF2test.hs b/test/GF2test.hs index 45f7bdf..30d1b82 100644 --- a/test/GF2test.hs +++ b/test/GF2test.hs @@ -1,43 +1,76 @@ module GF2test (run) where +import Data.Bits (xor) +import Data.List (foldl') import Test.HUnit import LinearAlgebra.GF2 ( fromBools - , gaussianEliminationIndices + , gaussianEliminationAllMasks + , maskToIndicesInt ) +-- XOR of the rows selected by a mask +xorSelected :: [Integer] -> [Int] -> Integer +xorSelected rows idxs = foldl' xor 0 [ rows !! i | i <- idxs ] + run :: IO () run = do putStrLn "GF2 Gaussian elimination test" - let rowsBools = - [ [True, False, True] - , [False, True, True] - , [True, True, True] - , [True, True, False] - ] + -- one dependency: rows 0, 1 and 3 XOR to zero + let rowsBits = + map fromBools + [ [True, False, True] + , [False, True, True] + , [True, True, True] + , [True, True, False] + ] + masks = gaussianEliminationAllMasks rowsBits - nCols = length (head rowsBools) - rowsBits = map fromBools rowsBools - maybeIndices = gaussianEliminationIndices nCols rowsBits + assertEqual "expected exactly one dependency" 1 (length masks) + assertEqual "expected rows 0,1,3" [[0,1,3]] (map maskToIndicesInt masks) - case maybeIndices of - Nothing -> assertFailure "expected rows 0,1,3" - Just rows -> assertEqual "expected rows 0,1,3" rows [0,1,3] + -- no dependency case: the three rows are independent + let rowsBits2 = + map fromBools + [ [True, False, True] + , [False, True, True] + , [True, True, True] + ] + assertEqual "Expected no dependency" [] (gaussianEliminationAllMasks rowsBits2) - -- no depency case - let rowsBools2 = - [ [True, False, True] - , [False, True, True] - , [True, True, True] - ] - nCols2 = length (head rowsBools2) - rowsBits2 = map fromBools rowsBools2 - maybeIndices2 = gaussianEliminationIndices nCols2 rowsBits2 + -- several dependencies: 5 rows of rank 2 give a null space of dimension 3 + let rowsBits3 = + map fromBools + [ [True, False, True] -- a + , [False, True, True] -- b + , [True, True, False] -- a xor b + , [True, False, True] -- a again + , [False, False, False] -- zero row + ] + masks3 = gaussianEliminationAllMasks rowsBits3 + idxs3 = map maskToIndicesInt masks3 - assertEqual "Expected no dependency" Nothing maybeIndices2 + assertEqual "expected rows - rank = 3 dependencies" 3 (length masks3) + assertBool "every dependency must XOR to zero" + (all (\ix -> xorSelected rowsBits3 ix == 0) idxs3) + assertBool "no dependency may be empty" (all (not . null) idxs3) + assertBool "the zero row (index 4) must be a dependency on its own" ([4] `elem` idxs3) + + -- the dependencies must be linearly independent (distinct and non-zero masks) + assertEqual "dependencies must be pairwise distinct" (length masks3) (length (dedup masks3)) + + -- a repeated row is a dependency of size two + let rowsBits4 = map fromBools [ [True, False], [False, True], [True, False] ] + assertEqual "duplicate row dependency" [[0,2]] + (map maskToIndicesInt (gaussianEliminationAllMasks rowsBits4)) + + -- empty input + assertEqual "no rows, no dependency" [] (gaussianEliminationAllMasks []) putStrLn "GF2 test done." - + where + dedup [] = [] + dedup (x : xs) = x : dedup (filter (/= x) xs)