Jun 4, 2015
The latest version of hmatrix supports basic arithmetic operations for integer (Int32 and Int64) elements.
let m = (3><4) [0..] :: Matrix I
m
(3><4) [ 0, 1, 2, 3 , 4, 5, 6, 7 , 8, 9, 10, 11 ]
m <> tr m
(3><3) [ 14, 38, 62 , 38, 126, 214 , 62, 214, 366 ]
This is not very useful at the moment because of the possibility of overflow; integer elements have been included mainly for improved matrix slicing and data manipulation. But once raw integer elements are admitted modular arithmetic can also be easily supported. We have defined a newtype wrapper with a phantom type for statically checked modular operations and the necessary instances of hmatrix classes.
12*7 :: I ./. 15
9
let m = (3><4) [0..] :: Matrix (I ./. 5)
m
(3><4) [ 0, 1, 2, 3 , 4, 0, 1, 2 , 3, 4, 0, 1 ]
m <> tr m
(3><3) [ 4, 3, 2 , 3, 1, 4 , 2, 4, 1 ]
A fundamental result in number theory is that the integers modulo a prime have multiplicative inverses. This is a closely related to Fermat's Little Theorem.
type F n = I ./. n
5 ^ (13-1) :: F 13
1
5 * 5^(13-2) :: F 13
1
Taking advantage of this result we have defined a Fractional instance for the modular types which works for any prime p:
recip 5 :: F 13
8
5*8 :: F 13
1
7/5 :: F 13
4
4*5 :: F 13
7
In order to solve a system of linear equations we only need that the coefficientes belong to a field, so taking advantage of the above Fractional instance we can solve systems with coefficients in the finite field F p. To check this we have defined a general version of Gaussian elimination (actually, LU decomposition and forward/backward substitution) using the tools described here.
linearSolve' a b = luSolve' (luPacked' a) b inv' x = linearSolve' x (ident (rows x))
a = (3><3) [ 1, 2, 3 , 4, 5, 6 , 7, 8, 0 ] :: Matrix (F 17)
b
(3><4) [ 15, 4, 10, 16 , 0, 15, 13, 11 , 15, 13, 11, 9 ]
linearSolve' a b
(3><4) [ 0, 1, 2, 3 , 4, 5, 6, 7 , 8, 9, 10, 11 ]
Using standard real arithmetic we would get a noninteger solution.
An arbitrary logic circuit composed exclusively of XOR gates can be expressed as a linear transformation of vectors in F 2. It is nice that we can use basic linear algebra to obtain the inverse circuit (if it is nonsingular).
dispbin $ circ
1 1 1 1
1
1 1 1 1 1
1 1 1 1
1 1 1
1 1 1 1 1 1
1 1 1 1 1 1 1
1 1 1 1 1 1
1 1 1 1
1 1 1 1
let icirc = inv' circ
dispbin $ icirc
1 1 1 1 1
1 1 1 1 1 1
1 1 1 1
1 1 1 1 1
1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1 1 1 1 1
1 1 1
dispbin $ icirc <> circ
1
1
1
1
1
1
1
1
1
1
A simple vector of elements with no additional structure has similar properties to the elements themselves, replicated in several independent copies. Things become more interesting when we define a multiplication operation to combine and produce vectors of different dimensions. We must only define the result for the canonical basis. A very simple possibility is:
$ e _i \times e _j = e _{i+j} $
This product is associative and $e_0$ is the identity, so we have a monoid. If the basis is written as powers of a variable ($e _k=x^k$) this is the familiar ring of polynomials.
Univariate polynomials can be implemented using arrays of coefficients, with element-by-element addition and multiplication based on convolution.
let y = Poly (fromList [0,1]) :: Poly Double
let a = y^2-3*y+2; b = 3*y^5-7
a*b
Poly (fromList [-14.0,21.0,-7.0,0.0,0.0,6.0,-9.0,3.0])
let c = a*b + 3*y -2
c
Poly (fromList [-16.0,24.0,-7.0,0.0,0.0,6.0,-9.0,3.0])
In conventional notation: $ c = 3.0y^7-9.0y^6+6.0y^5-7.0y^2+24.0y-16.0 $.
If the coefficients belong to a field we can invert the above multiplication using long division, obtaining a quotient and remainder:
let (q,r) = quotRemP c a
(q,r)
(Poly (fromList [-7.0,0.0,0.0,0.0,0.0,3.0]),Poly (fromList [-2.0,3.0]))
$ q = 3.0y^5-7.0 $ and $ r = 3.0y-2.0 $.
a*q+r - c
Poly (fromList [0.0])
(Multiplication of polynomials is equivalent to convolution of coefficient vectors. In principle we could implement division of polynomials by solving a linear system whose coeffient matrix is the expanded convolution operator, similar to the example in this page about sparse systems. This works, but it is probably not what we really want, because it always obtains zero remainder and an exact quotient, possibly of the same degree as the dividend.)
Anyway, the above operations with polynomials work exactly the same with coefficients in a finite field:
let y = Poly (fromList [0,1]) :: Poly (F 13)
let a = y^2-3*y+2; b = 3*y^5-7; c = a*b + 3*y -2
c
Poly (fromList [10,11,6,0,0,6,4,3])
We get the "same" result (mod 13). In conventional notation, $ c = 3y^7+4y^6+6y^5+6y^2+11y+10 $.
Remarkably, we can perform division of this kind of integer polynomials exactly in the same way as if they had real coefficients:
let (q,r) = quotRemP c a
(q,r)
(Poly (fromList [6,0,0,0,0,3]),Poly (fromList [11,3]))
$ q = 3y^5+6 $ and $ r = 3y+11 $.
a*q+r
Poly (fromList [10,11,6,0,0,6,4,3])
That's good, but as in the above example, and similar to standard integer division, polynomial division may leave a remainder. General polynomials over a finite field can be added and multiplied to produce elements of increasing and unbounded degree/dimension.
However, if we redefine the product modulo a given polynomial, we obtain a finite field of nonprime dimension. The modulo operation keeps the dimension/degree bounded, and a multiplicative inverse can again be easily defined (e.g. using powers).
Let's see an example in the finite field $GF(5^6)$ with reducing polynomial $ x^6+x^5+3x^3+x^2+2x+1 $.
let a = 3*x^5+2*x^2+1; b = x^2+4*x^3+4
$ a = 3x^5+2x^2+1 $
$ b = 4x^3+x^2+4 $
let c = a*b
$ c = 2x^4+2x^3+x^2+x $
c/a
FF (Poly (fromList [4,0,1,4]))
$\frac{c}{a} = 4x^3+x^2+4 $
$\frac{c}{b} = 3x^5+2x^2+1 $
$ \frac{1}{a} = 3x^5+3x^4+3x^3+2x^2+3x+2 $
$ \frac{1}{b} = 3x^5+3x^4+3x^3+2x^2+3x+2 $
Wikipedia: Finite field arithmetic.
A Ruiz / 2015-07-09