cube/zenzicubi.co
cube/zenzicubi.co
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

haskellify finite-field.1

Browse Source
This commit is contained in:
queue-miscreant 2025-08-01 04:30:27 -05:00
parent 86778da52a
commit 9385797c63
1 changed files with 164 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

228
posts/finite-field/1/index.qmd
View File

@ -19,15 +19,25 @@ categories:
.cayley-table th:nth-child(1) {
color: green;
}
.red {
color: red;
}
.green {
color: green;
}
</style>
```{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.
<!-- TODO: python for typesetting, haskell for implementation? -->
| 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 | <span class="green"> 5 </span> |
| $x^3 + x + 1$ | 11 | 7 |
| $x^3 + x^2 + 1$ | 13 | 11 |
| $x^4 + x + 1$ | 19 | 13 |
| $x^4 + x^3 + 1$ | <span class="red"> 25 </span> | <span class="green"> 17 </span> |
| $x^4 + x^3 + x^2 + x + 1$ | 31 | 19 |
| $x^5 + x^2 + 1$ | 37 | <span class="green"> 23 </span> |
| $x^5 + x^3 + 1$ | 41 | <span class="green"> 29 </span> |
| $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*)
</summary>
```{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
```
</details>
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 = "<span style=\"color: green\">" ++ x ++ "</span>"
redSpan x = "<span style=\"color: red\">" ++ x ++ "</span>"
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 uses |+| operation, which is not defined in this article -->
```{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]]
```
</details>
@ -588,8 +669,6 @@ Haskell implementation of the companion matrix
</summary>
```{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]
```
</details>