GNFS implementation + GF2 gaussian elimination fix

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)

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
 

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

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
+         else
-                  allExps = map snd pairs
+           let bitvecs = map (fromBools . map (\e -> e `mod` 2 == 1) . snd) pairs
-                  bitvecs = map (fromBools . map (\e -> e `mod` 2 == 1)) allExps
+               -- a basis of ALL dependencies between the rows (left null space over GF(2))
-              in case gaussianEliminationMask numCols bitvecs of
+               masks = gaussianEliminationAllMasks bitvecs
-                   Nothing -> Left "No linear dependency found"
+
-                   Just mask ->
+               -- Try one dependency: x = product of the chosen x_i, z = sqrt of the product of the chosen (x_i^2 mod n).
-                     let indices = maskToIndicesInt mask
+               -- Each dependency gives a non-trivial factor with probability about 1/2, so we move on to the next one when it fails.
-                     in if null indices
+               tryMask mask =
-                        then Left "Empty dependency"
+                 let chosen   = [ pr | (i, pr) <- zip [0 :: Int ..] pairs, testBit mask i ]
-                        else let productX = foldl (\acc idx -> modMul acc (xs !! idx) n) 1 indices
+                     productX = foldl (\acc (x, _) -> modMul acc x n) 1 chosen
-                                 fs = [ sum [ (allExps !! idx !! j) | idx <- indices ] `div` 2 | j <- [0 .. numCols - 1] ]
+                     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)
+                     z        = foldl (\acc (p, f) -> if f == 0 then acc else modMul (modExp p f n) acc n)
-                                 diff = productX - z
+                                      1 (zip base fs)
-                                 g = gcd (abs diff) n
+                     g        = gcd (abs (productX - z)) n
-                             in if g == 1 || g == n
+                 in if g == 1 || g == n then Nothing else Just g
-                                then Left "non trivial factors not found, try larger b or maximum x"
+           in case masks of
-                                else Right (g, n `div` g)
+                [] -> 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"

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)

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

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

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

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

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)

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 ]

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

src/LinearAlgebra/GF2.hs

@@ -6,14 +6,14 @@ module LinearAlgebra.GF2
   , fromBools
   , toBools
     -- Gaussian elimination
-  , gaussianEliminationMask     -- :: nCols -> [BitVec] -> Maybe Mask
+  , gaussianEliminationAllMasks -- :: [BitVec] -> [Mask]
-  , gaussianEliminationIndices  -- :: nCols -> [BitVec] -> Maybe [Int]
   , 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).
+-- 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.
-gaussianEliminationMask :: Int -> [BitVec] -> Maybe Mask
+gaussianEliminationAllMasks :: [BitVec] -> [Mask]
-gaussianEliminationMask nCols rows
+gaussianEliminationAllMasks rows = go M.empty (zip rows (map bit [0 ..])) []
-  | null rows = Nothing
-  | otherwise = findZeroMask (eliminate 0 augmented)
   where
-    m :: Int
+    go _ [] acc = reverse acc
-    m = length rows
+    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
+    reduce basis v mk
-    augmented :: [(BitVec, Mask)]
+      | v == 0    = (0, mk)
-    augmented = zip rows (map (\i -> bit i) [0..m-1])
+      | 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
+    lowBit :: Integer -> Integer
-    eliminate :: Int -> [(BitVec, Mask)] -> [(BitVec, Mask)]
+    lowBit v = v .&. negate v
-    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')
 
-    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
 

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 =
+    -- 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)
+    assertEqual "expected exactly one dependency" 1 (length masks)
-        rowsBits = map fromBools rowsBools
+    assertEqual "expected rows 0,1,3" [[0,1,3]] (map maskToIndicesInt masks)
-        maybeIndices = gaussianEliminationIndices nCols rowsBits
 
-    case maybeIndices of
+    -- no dependency case: the three rows are independent
-        Nothing   -> assertFailure "expected rows 0,1,3"
+    let rowsBits2 =
-        Just rows -> assertEqual "expected rows 0,1,3" rows [0,1,3]
+          map fromBools
-
-
-    -- 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
 
-    assertEqual "Expected no dependency" Nothing maybeIndices2
+    assertEqual "Expected no dependency" [] (gaussianEliminationAllMasks 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 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)