diff --git a/posts/finite-field/1/index.qmd b/posts/finite-field/1/index.qmd index 3a1c4ea..d18b53c 100644 --- a/posts/finite-field/1/index.qmd +++ b/posts/finite-field/1/index.qmd @@ -19,15 +19,25 @@ categories: .cayley-table th:nth-child(1) { color: green; } - -.red { - color: red; -} -.green { - color: green; -} +```{haskell} +--|echo: false + +import Data.List (intercalate) +import Data.Profunctor (rmap) + +import Colonnade +import qualified Colonnade.Encode as CE +import IHaskell.Display (markdown) + +markdownTable col rows = unlines $ h:h':r where + toColumns = ("| " ++) . (++ " |") . intercalate " | " . foldr (:) [] + h = toColumns $ CE.header id col + h' = [if x == '|' then x else '-' | x <- h] + r = map (toColumns . CE.row id col) rows +``` + [Fields](https://en.wikipedia.org/wiki/Field_%28mathematics%29) are a basic structure in abstract algebra. Roughly, a field is a collection of elements paired with two operations, @@ -162,10 +172,11 @@ This implementation actually works for any base *b*, which is not necessarily pr The only difference is that the coefficients lose "field-ness" for composite *b*. ```{haskell} --- | eval: false +import Data.List (unfoldr) +import Data.Tuple (swap) -- A polynomial is its ascending list of coefficients (of type a) -data Polynomial a = Poly { coeffs :: [a] } +newtype Polynomial a = Poly { coeffs :: [a] } deriving Functor -- Interpret a number's base-b expansion as a polynomial asPoly :: Int -> Int -> Polynomial Int @@ -200,16 +211,16 @@ For GF(2), this is insignificant since 1 is the only nonzero element, but over G $$ \begin{align*} - x^2 + 2 \text{over $GF(5)$} \quad + x^2 + 2 \text{ over GF$(5)$} \quad &\longleftrightarrow \quad {_5 27} \\ - 2x^2 + 4 \text{over $GF(5)$} \quad + 2x^2 + 4 \text{ over GF$(5)$} \quad &\longleftrightarrow \quad {_5 54} \\ - 3x^2 + 1 \text{over $GF(5)$} \quad + 3x^2 + 1 \text{ over GF$(5)$} \quad &\longleftrightarrow \quad {_5 76} \\ - 4x^2 + 3 \text{over $GF(5)$} \quad + 4x^2 + 3 \text{ over GF$(5)$} \quad &\longleftrightarrow \quad {_5 103} \end{align*} $$ @@ -229,8 +240,6 @@ Haskell implementation of monic polynomials over GF(*p*) Again, nothing about this definition depends on the base being prime. ```{haskell} --- | eval: false - -- All monic polynomials of degree d with coefficients mod n monics :: Int -> Int -> [Polynomial Int] monics n d = map (asPoly n) [n^d..2*(n^d) - 1] @@ -247,7 +256,7 @@ Unfortunately, these base-*p* expansions are more difficult to obtain algorithmi and I'll leave this as an exercise to the reader. -### Irreducibles +### Sieving out Irreducibles Over the integers, we can factor a number into primes. To decide if a number is prime, we just divide it (using an algorithm like long division) @@ -265,38 +274,11 @@ However, the converse does not hold. Over GF(2), $$ -(x + 1)^2 = x^2 + 2x + 1 \equiv x^2 + 1 \text{over $GF(2)$} +(x + 1)^2 = x^2 + 2x + 1 \equiv x^2 + 1 \text{ over GF$(2)$} $$ ...but as just mentioned, the right-hand side is irreducible over the integers. -Since we can denote polynomials by numbers, it may be tempting to freely switch - between primes and irreducibles. -However, irreducibles depend on the chosen field and do not generally correspond - to the base-*p* expansion of a prime. - - -| Irreducible over GF(2), *q(x)* | *q*(2) ([OEIS A014580](https://oeis.org/A014580)) | Prime | -|--------------------------------|---------------------------------------------------|---------------------------------| -| $x$ | 2 | 2 | -| $x + 1$ | 3 | 3 | -| $x^2 + x + 1$ | 7 | 5 | -| $x^3 + x + 1$ | 11 | 7 | -| $x^3 + x^2 + 1$ | 13 | 11 | -| $x^4 + x + 1$ | 19 | 13 | -| $x^4 + x^3 + 1$ | 25 | 17 | -| $x^4 + x^3 + x^2 + x + 1$ | 31 | 19 | -| $x^5 + x^2 + 1$ | 37 | 23 | -| $x^5 + x^3 + 1$ | 41 | 29 | -| $x^5 + x^3 + x^2 + x + 1$ | 47 | 31 | - -The red entry in column 2 is not prime. -Dually, the green entries in column 3 do not have binary expansions which correspond - to irreducible polynomials over GF(2). - - -### Dividing and Sieving - Just like integers, we can use [polynomial long division]( https://en.wikipedia.org/wiki/Polynomial_long_division ) with these objects to decide if a polynomial is irreducible. @@ -314,7 +296,10 @@ The algorithm is similar to table-less algorithms for CRCs, but we don't have th We also have to watch out for negation and coefficients other than 1 for when not working over GF(2). ```{haskell} --- | eval: false +zipAdd :: Num a => [a] -> [a] -> [a] +zipAdd [] ds = ds +zipAdd cs [] = cs +zipAdd (c:cs) (d:ds) = (c + d):zipAdd cs ds -- Divide the polynomial ps by qs (coefficients in descending degree order) synthDiv' :: (Eq a, Num a) => [a] -> [a] -> ([a], [a]) @@ -332,7 +317,8 @@ synthDiv' ps qs doDiv (x:xs) d = x:doDiv (zipAdd xs $ map (*x) qNeg) (d - 1) -- Use Polynomial (coefficients in ascending degree order) instead of lists -synthDiv :: (Eq a, Num a) => Polynomial a -> Polynomial a -> (Polynomial a, Polynomial a) +synthDiv :: (Eq a, Num a) => Polynomial a -> Polynomial a + -> (Polynomial a, Polynomial a) synthDiv (Poly p) (Poly q) = (Poly $ reverse quot, Poly $ reverse rem) where (quot, rem) = synthDiv' (reverse p) (reverse q) ``` @@ -347,11 +333,9 @@ Haskell implementation of an irreducible polynomial sieve over GF(*p*) ```{haskell} --- | eval: false - -- All irreducible monic polynomials with coefficients mod n irreducibles :: Int -> [Polynomial Int] -irreducibles n = go [] $ monics n where +irreducibles n = go [] $ allMonics n where -- Divide the polynomial x by i, then take the remainder mod n remModN x i = fmap (`mod` n) $ snd $ synthDiv x i -- Find remainders of x divided by every irreducible in "is". @@ -361,11 +345,59 @@ irreducibles n = go [] $ monics n where go is (x:xs) | notMultiple x is = x:go (x:is) xs | otherwise = go is xs - -print $ take 10 $ irreducibles 2 ``` +Since we can denote polynomials by numbers, it may be tempting to freely switch + between primes and irreducibles. +However, irreducibles depend on the chosen field and do not generally correspond + to the base-*p* expansion of a prime. + +```{haskell} +--|code-fold: true +--|classes: plain + +texifyPoly :: (Num a, Eq a, Show a) => Polynomial a -> String +texifyPoly (Poly xs) = ("$" ++) $ (++ "$") $ texify' $ zip xs [0..] where + texify' [] = "0" + texify' ((c, n):xs) + | all ((==0) . fst) xs = showPow c n + | c == 0 = texify' xs + | otherwise = showPow c n ++ " + " ++ texify' xs + showPow c 0 = show c + showPow 1 1 = "x" + showPow c 1 = show c ++ showPow 1 1 + showPow 1 n = "x^{" ++ show n ++ "}" + showPow c n = show c ++ showPow 1 n + +-- Simple prime sieve +primes = primes' [] 2 where + primes' ps n + | and [n `mod` p /= 0 | p <- ps] = n:primes' (n:ps) (n+1) + | otherwise = primes' ps (n+1) + +somePrimes = takeWhile (<50) primes +someIrr2 = takeWhile (<50) $ map (evalPoly 2) $ irreducibles 2 + +irr2PrimeTable = columns (\(_, f) r -> f r) (\(c, _) -> Headed c) [ + ("Irreducible over GF(2), *q*(x)", texifyPoly . (irreducibles 2 !!)), + ("*q*(2), [OEIS A014580](https://oeis.org/A014580)", + (\x -> if x `notElem` somePrimes then redSpan $ show x else show x) . + (someIrr2 !!)), + ("Prime", + (\x -> if x `notElem` someIrr2 then greenSpan $ show x else show x) . + (primes !!)) + ] where + greenSpan x = "" ++ x ++ "" + redSpan x = "" ++ x ++ "" + +markdown $ markdownTable irr2PrimeTable [0..10] +``` + +The red entry in column 2 is not prime. +Dually, the green entries in column 3 do not have binary expansions which correspond + to irreducible polynomials over GF(2). + Matrices -------- @@ -395,9 +427,14 @@ Laplace expansion is ludicrously inefficient compared to other algorithms, Numeric computation will not be used to keep the arithmetic exact. ```{haskell} --- | eval: false +import Data.Array +newtype Matrix a = Mat { unMat :: Array (Int, Int) a } deriving Functor -newtype Matrix a = Mat { unMat :: Array (Int, Int) a } +-- Simple function for building a Matrix from lists +toMatrix :: [[a]] -> Matrix a +toMatrix l = Mat $ listArray ((0,0),(n-1,m-1)) $ concat l where + m = length $ head l + n = length l determinant :: (Num a, Eq a) => Matrix a -> a determinant (Mat xs) = determinant' xs where @@ -414,7 +451,7 @@ determinant (Mat xs) = determinant' xs where -- Produce the cofactor of row 0, column i cofactor i | xs!(0,i) == 0 = 0 - | otherwise = (parity i) * xs!(0,i) * (determinant' $ minor i) + | otherwise = (parity i) * xs!(0,i) * determinant' (minor i) -- Furthest extent of the bounds, i.e., the size of the matrix (_,(n,_)) = bounds xs -- Build a new Array by eliminating row 0 and column i @@ -454,17 +491,59 @@ Haskell implementation of the characteristic polynomial Since `determinant` was defined for all `Num` and `Eq`, it can immediately be applied if these instances are defined for polynomials. - - ```{haskell} --- | eval: false --- Num instance for polynomials omitted --- instance (Num a, Eq a) => Num (Polynomial a) where --- ... +instance (Eq a, Num a) => (Num (Polynomial a)) where + (+) x@(Poly as) y@(Poly bs) = Poly $ zipAdd as bs + + --convolution + (*) (Poly []) (Poly bs) = Poly [] + (*) (Poly as) (Poly []) = Poly [] + (*) (Poly as) (Poly bs) = Poly $ zeroize $ finish $ foldl convolve' ([], []) as + where convolve' (xs, ys) a = (a:xs, sum (zipWith (*) (a:xs) bs):ys) + finish (xs, ys) = (reverse ys ++) $ finish' xs $ tail bs + finish' xs [] = [] + finish' xs ys = sum (zipWith (*) xs ys):finish' xs (tail ys) + zeroize xs = if all (==0) xs then [] else xs + + -- this definition works + negate = fmap negate + -- but these two might run into some problems + abs = fmap abs + signum = fmap signum + fromInteger 0 = Poly [] + fromInteger a = Poly [fromInteger a] instance Eq a => Eq (Polynomial a) where (==) (Poly xs) (Poly ys) = xs == ys +``` +Then, along with some helper matrix functions, we can build a function for the characteristic polynomial. + +```{haskell} +-- Zero matrix +zero :: Num a => Int -> Matrix a +zero n = Mat $ listArray ((0,0),(n-1,n-1)) $ repeat 0 + +-- Identity matrix +eye :: Num a => Int -> Matrix a +eye n = Mat $ unMat (zero n) // take n [((i,i), 1) | i <- [0..]] + +-- Maps +mapRange :: Ix i => (i -> e) -> (i, i) -> [(i, e)] +mapRange g r = map (\x -> (x, g x)) $ range r + +-- Pointwise application +zipWithArr :: Ix i => (a -> b -> c) + -> Array i a -> Array i b -> Array i c +zipWithArr f a b + | ab == bb = array ab $ map (\x -> (x, f (a!x) (b!x))) $ indices a + | otherwise = error "Array dimension mismatch" where + ab = bounds a + bb = bounds b + +(|+|) (Mat x) (Mat y) = Mat $ zipWithArr (+) x y + +-- Characteristic Polynomial charpoly :: Matrix Int -> Polynomial Int charpoly xs = determinant $ eyeLambda |+| negPolyXs where -- Furthest extent of the bounds, i.e., the size of the matrix @@ -474,7 +553,9 @@ charpoly xs = determinant $ eyeLambda |+| negPolyXs where negPolyXs = fmap (\x -> Poly [-x]) xs -- Identity matrix times lambda (encoded as Poly [0, 1]) eyeLambda :: Matrix (Polynomial Int) - eyeLambda = fmap (\x -> (Poly [x] * Poly [0, 1])) $ eye (n+1) + eyeLambda = (\x -> Poly [x] * Poly [0, 1]) <$> eye (n+1) + +markdown $ texifyPoly $ charpoly $ toMatrix [[1,0],[0,1]] ``` @@ -588,8 +669,6 @@ Haskell implementation of the companion matrix ```{haskell} --- | eval: false - companion :: Polynomial Int -> Matrix Int companion (Poly ps) | last ps' /= 1 = error "Cannot find companion matrix of non-monic polynomial" @@ -600,13 +679,34 @@ companion (Poly ps) ps' = reverse $ dropWhile (==0) $ reverse ps -- Address/value tuples for a shifted identity matrix: -- 1s on the diagonal just above the main diagonal, 0s elsewhere - shiftI = map (\p@(x,y) -> (p, if y == x + 1 then 1 else 0)) $ range ((0,0),(n-2,n-1)) + shiftI = map (\p@(x,y) -> (p, if y == x + 1 then 1 else 0)) $ + range ((0,0),(n-2,n-1)) -- Address/value tuples for the last row of the companion matrix: -- ascending powers of the polynomial lastRow = zipWith (\x y -> ((n-1, x), y)) [0..n-1] $ map negate ps' -- (charpoly . companion) = id :: Polynomial Int -> Polynomial Int ``` + +Applying this to $1 - 2x + x^2$ gives us: + +```{haskell} +--| code-fold: true + +reshape :: Int -> [a] -> [[a]] +reshape n = unfoldr (reshape' n) where + reshape' n x = if null x then Nothing else Just $ splitAt n x + +fromMatrix :: Matrix a -> [[a]] +fromMatrix (Mat m) = let (_,(_,n)) = bounds m in reshape (n+1) $ elems m + +texifyMatrix mat = surround mat' where + mat' = intercalate " \\\\ " $ map (intercalate " & " . map show) $ + fromMatrix mat + surround = ("\\left( \\begin{matrix}" ++) . (++ "\\end{matrix} \\right)") + +markdown $ ("$$" ++) $ (++ "$$") $ texifyMatrix $ companion $ Poly [1, -2, 1] +```