haskellify finite-field.2

posts/finite-field/2/MplIHaskell.hs

@@ -0,0 +1 @@
+../../number-number/1/MplIHaskell.hs

posts/finite-field/2/Previous.hs

@@ -0,0 +1,225 @@
+module Previous where
+
+import Control.Applicative (liftA2)
+import Data.Array
+import Data.List (unfoldr, transpose)
+import Data.Tuple (swap)
+
+zipAdd :: Num a => [a] -> [a] -> [a]
+zipAdd []     ds     = ds
+zipAdd cs     []     = cs
+zipAdd (c:cs) (d:ds) = (c + d):zipAdd cs ds
+
+{------------------------------------------------------------------------------
+ -  ___     _                     _      _
+ - | _ \___| |_  _ _ _  ___ _ __ (_)__ _| |___
+ - |  _/ _ \ | || | ' \/ _ \ '  \| / _` | (_-<
+ - |_| \___/_|\_, |_||_\___/_|_|_|_\__,_|_/__/
+ -            |__/
+ ------------------------------------------------------------------------------}
+
+-- A polynomial is its ascending list of coefficients (of type a)
+newtype Polynomial a = Poly { coeffs :: [a] } deriving Functor
+
+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, Num a) => Eq (Polynomial a) where
+  (==) as bs = all (==0) $ coeffs $ as - bs
+
+
+-- Interpret a number's base-b expansion as a polynomial
+asPoly :: Int -> Int -> Polynomial Int
+-- Build a list with f, which returns either Nothing
+-- or Just (next element of list, next argument to f)
+asPoly b = Poly . unfoldr f where
+  -- Divide x by b. Emit the remainder and recurse with the quotient.
+  f x | x /= 0    = Just $ swap $ divMod x b
+  -- If there's nothing left to divide out, terminate
+      | otherwise = Nothing
+
+-- Horner evaluation of a polynomial at the integer b
+evalPoly :: Int -> Polynomial Int -> Int
+-- Start with the highest coefficient
+-- Multiply by b at each step and add the coefficient of the next term
+evalPoly b (Poly p) = foldr (\y acc -> acc*b + y) 0 p
+
+-- Divide the polynomial ps by qs (coefficients in descending degree order)
+synthDiv' :: (Eq a, Num a) => [a] -> [a] -> ([a], [a])
+synthDiv' ps qs
+  | head qs /= 1 = error "Cannot divide by non-monic polynomial"
+  | otherwise   = splitAt deg $ doDiv ps deg
+  where
+    -- Negate the denominator and ignore leading term
+    qNeg = map negate $ tail qs
+    -- The degree of the result, based on the degrees of the numerator and denominator
+    deg   = max 0 (length ps - length qs + 1)
+    -- Pluck off the head of the list and add a shifted and scaled version of
+    -- qs to the tail of the list. Repeat this d times
+    doDiv xs    0  = xs
+    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 (Poly p) (Poly q) = (Poly $ reverse quot, Poly $ reverse rem)
+  where (quot, rem)        = synthDiv' (reverse p) (reverse q)
+
+
+{------------------------------------------------------------------------------
+ -  ___     _                     _      _    ___
+ - | _ \___| |_  _ _ _  ___ _ __ (_)__ _| |  / __| ___ __ _ _  _ ___ _ _  __ ___ ___
+ - |  _/ _ \ | || | ' \/ _ \ '  \| / _` | |  \__ \/ -_) _` | || / -_) ' \/ _/ -_|_-<
+ - |_| \___/_|\_, |_||_\___/_|_|_|_\__,_|_|  |___/\___\__, |\_,_\___|_||_\__\___/__/
+ -            |__/                                       |_|
+ ------------------------------------------------------------------------------}
+
+-- 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]
+
+-- All monic polynomials with coefficients mod n, ordered by degree
+allMonics :: Int -> [Polynomial Int]
+allMonics n = concat [monics n d | d <- [1..]]
+
+-- All irreducible monic polynomials with coefficients mod n
+irreducibles :: Int -> [Polynomial Int]
+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".
+  -- If any give the zero polynomial, then x is a multiple of an irreducible
+  notMultiple x is = and [not $ all (==0) $ coeffs $ remModN x i | i <- is]
+  -- Sieve out by notMultiple
+  go is (x:xs)
+    | notMultiple x is = x:go (x:is) xs
+    | otherwise        = go is xs
+
+
+{------------------------------------------------------------------------------
+ -  __  __      _       _
+ - |  \/  |__ _| |_ _ _(_)__ ___ ___
+ - | |\/| / _` |  _| '_| / _/ -_|_-<
+ - |_|  |_\__,_|\__|_| |_\__\___/__/
+ -
+ ------------------------------------------------------------------------------}
+
+newtype Matrix a = Mat { unMat :: Array (Int, Int) a } deriving (Functor, Eq)
+
+asScalar = Mat . listArray ((0,0),(0,0)) . pure
+
+transposeM :: Matrix a -> Matrix a
+transposeM (Mat m) = Mat $ ixmap (bounds m) swap m
+
+dotMul :: Num a => [a] -> [a] -> a
+dotMul = (sum .) . zipWith (*)
+
+matMul' xs ys = reshape (length ys) $ liftA2 dotMul xs ys
+matMul  xs ys = matMul' xs $ transpose ys
+
+-- 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
+
+-- Maps
+mapRange :: Ix i => (i -> e) -> (i, i) -> [(i, e)]
+mapRange g r = map (\x -> (x, g x)) $ range r
+
+instance Num a => Num (Matrix a) where
+  (+) (Mat x) (Mat y)
+    | (0,0) == (snd $ bounds x) = Mat $ fmap ((x!(0,0))+) y --adding scalars
+    | otherwise                 = Mat $ zipWithArr (+) x y
+  (*) x y
+    | (0,0) == (snd $ bounds $ unMat x) = fmap (((unMat x)!(0,0))*) y --multiplying scalars
+    | otherwise                         = toMatrix $ matMul' (fromMatrix x) (fromMatrix $ transposeM y)
+  abs = fmap abs
+  negate = fmap negate
+  signum = undefined
+  fromInteger = asScalar . fromInteger --"scalars"
+
+reshape :: Int -> [a] -> [[a]]
+reshape n = unfoldr (reshape' n) where
+  reshape' n x = if null x then Nothing else Just $ splitAt n x
+
+-- 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
+
+-- ...and its inverse
+fromMatrix :: Matrix a -> [[a]]
+fromMatrix (Mat m) = let (_,(_,n)) = bounds m in reshape (n+1) $ elems m
+
+-- 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..]]
+
+-- Companion matrix
+companion :: Polynomial Int -> Matrix Int
+companion (Poly ps)
+  | last ps' /= 1 = error "Cannot find companion matrix of non-monic polynomial"
+  | otherwise     = Mat $ array ((0,0), (n-1,n-1)) $ lastRow ++ shiftI where
+    -- The degree of the polynomial, as well as the size of the matrix
+    n       = length ps' - 1
+    -- Remove trailing 0s from 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))
+    -- 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'
+
+mDims = snd . bounds . unMat
+
+-- trace of a matrix: sum of diagonal entries
+mTrace :: (Num a) => Matrix a -> a
+mTrace a = sum [(unMat a)!(i,i) | i <- [0..(uncurry min $ mDims a)]]
+
+-- Faddeev-Leverrier algorithm for computing the determinant, characteristic polynomial, and inverse matrix
+-- this relies on dividing the trace by the current iteration count, so it is only suitable for Integral types
+fadeev :: (Num a, Integral a) => Matrix a -> ([a], a, Matrix a)
+fadeev a = fadeev' [1] (eye $ n+1) 1 where
+  n = snd $ mDims a
+  fadeev' rs m i
+    | all (==0) (unMat nextm) = (replicate (n - (fromIntegral i) + 1) 0 ++ tram:rs , -tram, m)
+    | otherwise               = fadeev' (tram:rs) nextm (i+1) where
+      am    = a * m
+      tram  = -mTrace am `div` i
+      nextm = am + (fmap (tram *) $ eye (n+1))
+
+-- determinant from Faddeev-Leverrier algorithm
+determinant :: (Num a, Integral a) => Matrix a -> a
+determinant = (\(_,x,_) -> x) . fadeev
+-- characteristic polynomial from Faddeev-Leverrier algorithm
+charpoly :: (Num a, Integral a) => Matrix a -> Polynomial a
+charpoly    = Poly . (\(x,_,_) -> x) . fadeev

posts/finite-field/2/extra.qmd

@@ -0,0 +1,281 @@
+Rather than redefining evaluation for each of these cases,
+  we should map our polynomial into a structure compatible with how we want to evaluate it.
+Essentially, this means that from a polynomial in the base structure,
+  we can derive polynomials in these other structures.
+In particular, we can either have a matrix of polynomials or a polynomial in matrices.
+
+<!-- TODO: notes about functoriality of `fmap`ping eval vs -->
+:::: {layout-ncol="2"}
+::: {}
+$$
+\begin{align*}
+  p &: K[x]
+  \\
+  p(x) &= x^n + p_{n-1}x^{n-1} + ...
+  \\
+  \phantom{= p} & + p_1 x + p_0
+\end{align*}
+$$
+:::
+
+::: {}
+$x$ is a scalar indeterminate
+
+```haskell
+p :: Polynomial K
+```
+:::
+::::
+
+:::: {layout-ncol="2"}
+::: {}
+$$
+\begin{align*}
+  P &: (K[x])^{m \times m}
+  \\
+  P(x I) &= (x I)^n + (p_{n-1})(x I)^{n-1} + ...
+  \\
+  & + p_1(x I)+ p_0 I
+\end{align*}
+$$
+:::
+
+::: {}
+$x$ is a scalar indeterminate, $P(x I)= p(x) I$ is a matrix of polynomials in $x$
+
+```haskell
+asPolynomialMatrix
+  :: Polynomial K -> Matrix (Polynomial K)
+
+pMat :: Matrix (Polynomial K)
+pMat = asPolynomialMatrix p
+```
+:::
+::::
+
+:::: {layout-ncol="2"}
+::: {}
+$$
+\begin{align*}
+  \hat P &: K^{m \times m}[X]
+  \\
+  \hat P(X) &= X^n + (p_{n-1}I)X^{n-1} + ...
+  \\
+  & + (p_1 I) X + (p_0 I)
+\end{align*}
+$$
+:::
+
+::: {}
+$X$ is a matrix indeterminate, $\hat P(X)$ is a polynomial over matrices
+
+```haskell
+asMatrixPolynomial
+  :: Polynomial K -> Polynomial (Matrix K)
+
+pHat :: Polynomial (Matrix K)
+pHat = asMatrixPolynomial p
+```
+:::
+::::
+
+
+### Cayley-Hamilton Theorem
+
+When evaluating the characteristic polynomial of a matrix *with* that matrix,
+  something strange happens.
+Continuing from the previous article, using $x^2 + x + 1$ and its companion matrix, we have:
+
+$$
+\begin{gather*}
+  p(x) = x^2 + x + 1 \qquad C_{p} = C
+  = \left( \begin{matrix}
+       0 &  1 \\
+      -1 & -1
+    \end{matrix} \right)
+  \\ \\
+  \hat P(C) = C^2 + C + (1 \cdot I)
+    = \left( \begin{matrix}
+      -1 & -1 \\
+       1 &  0
+    \end{matrix} \right)
+    + \left( \begin{matrix}
+       0 &  1 \\
+      -1 & -1
+    \end{matrix} \right)
+    + \left( \begin{matrix}
+      1 & 0 \\
+      0 & 1
+    \end{matrix} \right)
+  \\ \\
+  = \left( \begin{matrix}
+      0 & 0 \\
+      0 & 0
+    \end{matrix} \right)
+\end{gather*}
+$$
+
+The result is the zero matrix.
+This tells us that, at least in this case, the matrix *C* is a root of its own characteristic polynomial.
+By the [Cayley-Hamilton theorem](https://en.wikipedia.org/wiki/Cayley%E2%80%93Hamilton_theorem),
+  this is true in general, no matter the degree of *p*, no matter its coefficients,
+  and importantly, no matter the choice of field.
+
+This is more powerful than it would otherwise seem.
+For one, factoring a polynomial "inside" a matrix turns out to give the same answer
+  as factoring a polynomial over matrices.
+
+:::: {layout-ncol="2"}
+::: {}
+
+$$
+\begin{gather*}
+  P(xI) = \left( \begin{matrix}
+      x^2 + x + 1 &           0 \\
+                0 & x^2 + x + 1
+    \end{matrix}\right)
+  \\ \\
+  = (xI - C)(xI - C')
+  \\ \\
+  = \left( \begin{matrix}
+      x &    -1 \\
+      1 & x + 1
+    \end{matrix} \right)
+    \left( \begin{matrix}
+      x - a &    -b \\
+         -c & x - d
+    \end{matrix} \right)
+  \\ \\
+  \begin{align*}
+    x(x - a) + c &= x^2 + x + 1
+    \\
+    \textcolor{green}{x(-b) - (x - d)} &\textcolor{green}{= 0}
+    \\
+    \textcolor{blue}{(x - a) + (x + 1)(-c)} &\textcolor{blue}{= 0}
+    \\
+    (-b) + (x + 1)(x - d) &= x^2 + x + 1
+  \end{align*}
+  \\ \\
+    \textcolor{green}{(-b -1)x +d = 0} \implies b = -1, ~ d = 0 \\
+    \textcolor{blue}{(1 - c)x - a - c = 0} \implies c = 1, ~ a = -1
+  \\ \\
+  C' =
+    \left( \begin{matrix}
+      -1 & -1 \\
+       1 &  0
+    \end{matrix} \right)
+\end{gather*}
+$$
+:::
+
+::: {}
+$$
+\begin{gather*}
+  \hat P(X) = X^2 + X + 1I
+  \\[10pt]
+  = (X - C)(X - C')
+  \\[10pt]
+  = X^2 - (C + C')X + CC'
+  \\[10pt]
+  \implies
+  \\[10pt]
+  C + C' = -I, ~ C' = -I - C
+  \\[10pt]
+  CC' = I, ~ C^{-1} = C'
+  \\[10pt]
+  C' = \left( \begin{matrix}
+      -1 & -1 \\
+       1 &  0
+    \end{matrix} \right)
+\end{gather*}
+$$
+:::
+::::
+
+It's important to not that a matrix factorization is not unique.
+*Any* matrix with a given characteristic polynomial can be used as a root of that polynomial.
+Of course, choosing one root affects the other matrix roots.
+
+
+### Moving Roots
+
+All matrices commute with the identity and zero matrices.
+A less obvious fact is that all of the matrix roots *also* commute with one another.
+By the Fundamental Theorem of Algebra,
+  [Vieta's formulas](https://en.wikipedia.org/wiki/Vieta%27s_formulas) state:
+
+$$
+\begin{gather*}
+  \hat P(X)
+    = \prod_{[i]_n} (X - \Xi_i)
+    = (X - \Xi_0) (X - \Xi_1)...(X - \Xi_{n-1})
+  \\
+  = \left\{ \begin{align*}
+      & \phantom{+} X^n
+      \\
+      & - (\Xi_0 + \Xi_1 + ... + \Xi_{n-1}) X^{n-1}
+      \\
+      & + (\Xi_0 \Xi_1+ \Xi_0 \Xi_2 + ... + \Xi_0 \Xi_{n-1} + \Xi_1 \Xi_2 + ... + \Xi_{n-2} \Xi_{n-1})X^{n-2}
+      \\
+      & \qquad \vdots
+      \\
+      & + (-1)^n \Xi_0 \Xi_1 \Xi_2...\Xi_n
+    \end{align*} \right.
+  \\
+  = X^n -\sigma_1([\Xi]_n)X^{n-1} + \sigma_2([\Xi]_n)X^{n-2} + ... + (-1)^n \sigma_n([\Xi]_n)
+\end{gather*}
+$$
+
+The product range \[*i*\]~*n*~ means that the terms are ordered from 0 to *n* - 1 over the index given.
+On the bottom line, *σ* are
+  [elementary symmetric polynomials](https://en.wikipedia.org/wiki/Elementary_symmetric_polynomial)
+  and \[*Ξ*\]~*n*~ is the list of root matrices from *Ξ*~*0*~ to Ξ~*n-1*~.
+
+By factoring the matrix with the roots in a different order, we get another factorization.
+It suffices to only focus on *σ*~2~, which has all pairwise products.
+
+$$
+\begin{gather*}
+  \pi \in S_n
+  \\
+  \qquad
+    \pi \circ \hat P(X) = \prod_{\pi ([i]_n)} (X - \Xi_i)
+  \\ \\
+  = X^n
+    - \sigma_1 \left(\pi ([\Xi]_n) \vphantom{^{1}} \right)X^{n-1} +
+    + \sigma_2 \left(\pi ([\Xi]_n) \vphantom{^{1}} \right)X^{n-2} + ...
+    + (-1)^n \sigma_n \left(\pi ([\Xi]_n) \vphantom{^{1}} \right)
+  \\ \\
+  \\ \\
+  (0 ~ 1) \circ \hat P(X) = (X - \Xi_{1}) (X - \Xi_0)(X - \Xi_2)...(X - \Xi_{n-1})
+  \\
+  = X^n + ... + \sigma_2(\Xi_1, \Xi_0, \Xi_2, ...,\Xi_{n-1})X^{n-2} + ...
+  \\ \\ \\ \\
+  \begin{array}{}
+    e & (0 ~ 1) & (1 ~ 2) & ... & (n-2 ~~ n-1)
+    \\ \hline
+    \textcolor{red}{\Xi_0 \Xi_1} & \textcolor{red}{\Xi_1 \Xi_0} & \Xi_0 \Xi_1  & & \Xi_0 \Xi_1
+    \\
+    \Xi_0 \Xi_2 & \Xi_0 \Xi_2 & \Xi_0 \Xi_2 & & \Xi_0 \Xi_2
+    \\
+    \Xi_0 \Xi_3  & \Xi_0 \Xi_3 & \Xi_0 \Xi_3 & & \Xi_0 \Xi_3
+    \\
+    \vdots & \vdots & \vdots & & \vdots
+    \\
+    \Xi_0 \Xi_{n-1}  & \Xi_0 \Xi_{n-1} & \Xi_{0} \Xi_{n-1} & & \Xi_{n-1} \Xi_0
+    \\
+    \textcolor{green}{\Xi_1 \Xi_2}  & \Xi_1 \Xi_2 & \textcolor{green}{\Xi_2 \Xi_1}  & & \Xi_1 \Xi_2
+    \\
+    \vdots & \vdots & \vdots & & \vdots
+    \\
+    \textcolor{blue}{\Xi_{n-2} \Xi_{n-1}}  & \Xi_{n-2} \Xi_{n-1} & \Xi_{n-2} \Xi_{n-1} & & \textcolor{blue}{\Xi_{n-1} \Xi_{n-2}}
+  \end{array}
+\end{gather*}
+$$
+
+<!-- TODO: permutation -->
+The "[path swaps](/posts/permutations/1/)" shown commute only the adjacent elements.
+By contrast, the permutation (0 2) commutes *Ξ*~0~ past both *Ξ*~1~ and *Ξ*~2~.
+But since we already know *Ξ*~0~ and *Ξ*~1~ commute by the above list,
+  we learn at this step that *Ξ*~0~ and *Ξ*~2~ commute.
+This can be repeated until we reach the permutation (0 *n*-1) to prove commutativity between all pairs.