New release of HaskellForMaths

I've just uploaded a new version v0.4.1 of HaskellForMaths, containing three new modules and a couple of other improvements. The additions are as follows:

Math.Algebras.Quaternions

This module was already present: it defines the quaternion algebra on the basis {1,i,j,k}, where multiplication is defined by:
i^2 = j^2 = k^2 = ijk = -1

This is enough information to figure out the full multiplication table. For example:
ijk = -1
=> (ijk)k = -k
=> ij(kk) = -k (associativity of multiplication)
=> ij = k
It turns out that the basis elements i,j,k anti-commute in pairs, eg ij = -ji, etc.

In this release I've added a couple of new things.

First, the quaternions are a division algebra, so I've added a Fractional instance.

Specifically, we can define a conjugation operation on the quaternions (similar to complex conjugation) via
conj (w+xi+yj+zk) = w-xi-yj-zk
Then we can define a quadratic norm via
sqnorm q = q * conj q = w^2+x^2+y^2+z^2
Since the norm is always a scalar, we can define a multiplicative inverse by
q^-1 = conj q / sqnorm q

For example:

$ cabal install HaskellForMaths
$ ghci
> :m Math.Algebras.Quaternions
> (2*i+3*j)^-1 :: Quaternion Q
-2/13i-3/13j

(If you leave out the type annotation, you'll be working in Quaternion Double.)

Second, the quaternions have an interesting role in 3- and 4-dimensional geometry.

Given any non-zero quaternion q, the map x -> q^-1 x q turns out to be a rotation of the 3-dimensional space spanned by {i,j,k}. To multiply rotations together (ie do one then another), just multiply the quaternions. This turns out to be a better way to represent rotations than 3*3 matrices:
- It's more compact - four scalars rather than nine
- They're faster to multiply - 16 scalar multiplications versus 27
- It's more robust against rounding error - whatever quaternion you end up with will still represent a rotation, whereas a sequence of matrix multiplications of rotations might not be quite a rotation any more, due to rounding error.

If you're curious, the function reprSO3 converts a quaternion to the corresponding 3*3 matrix:

> reprSO3 (1+2*i) :: [[Q]]
[[1,0,0],[0,-3/5,-4/5],[0,4/5,-3/5]]

(Exercise: Figure out why we got this matrix.)

Quaternions can also be used to represent rotations of 4-dimensional space - see the documentation.


Math.Algebras.Octonions

This is a new module, providing an implementation of the 8-dimensional non-associative division algebra of octonions. I follow Conway's notation [1], so the octonions have basis {1,e0,e1,e2,e3,e4,e5,e6}, with multiplication defined by:
e_i * e_i = -1, for i in [0..6]
e_i+1 * e_i+2 = e_i+4, where the indices are taken modulo 7.

The octonions are not associative, but they are an inverse loop, so they satisfy x^-1(xy) = y = (yx)x^-1. This is enough to enable us to deduce the full multiplication table from the relations above.

Like the quaternions, the octonions have conjugation and a norm, and multiplicative inverses:

> :l Math.Algebras.Octonions
> (2+i0+2*i3)^-1
2/9-1/9i0-2/9i3

The octonions are an exceptional object in mathematics: there's nothing else quite like them. They can be used to construct various other exceptional objects, such as the root lattice E8, or the Lie group G2. Hopefully I'll be able to cover some of that stuff in a future installment.

[1] Conway and Smith, On Quaternions and Octonions


Math.NumberTheory.Prime

The main function in this module is isPrime :: Integer -> Bool, which tells you whether a number is prime or not. It's implemented using the Miller-Rabin test.

The basic idea of the test is:
- If p is prime, then Zp is a field
- In a field, the equation x^2 = 1 has only two solutions, 1 and -1
- Given an arbitrary b coprime to p, we know from Fermat's little theorem that b^(p-1) = 1 (mod p)
- So if p-1 = q * 2^s, with q odd, then either b^q = 1 (mod p), or there is some r, 0 <= r < s with b^(q*2^r) = -1 (mod p)

The idea of the algorithm is to try to show that p isn't prime by trying to find a b where the above is not true. We take several different values of b at random, and repeatedly square b^q, to see whether we get -1 or not.

The advantage of the Miller-Rabin test, as compared to trial division say, is that it has a fast running time even for very large numbers. For example:

> :m Math.NumberTheory.Prime
> :set +s
> nextPrime $ 10^100
10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000267
(0.09 secs, 14904632 bytes)

The potential disadvantage of the Miller-Rabin test is that it is probabilistic: There is a very small chance (1 in 10^15 in this implementation) that it could just fail to hit on a b which disproves n's primeness, so that it would say n is prime when it isn't. In practice, at those odds it's not worth worrying about.


Math.NumberTheory.Factor

The main function in this module is pfactors :: Integer -> [Integer], which returns the prime factors of a number (with multiplicity). It uses trial division to try to find prime factors less than 10000. After that, it uses the elliptic curve method to try to split what remains. The elliptic curve method relies on some quite advanced maths, but the basic idea is this:
- If p is a prime, then Zp is a field
- Given a field, we can do "arithmetic on elliptic curves" over the field.
- So to factor n, pretend that n is prime, try doing arithmetic on elliptic curves, and wait till something goes wrong.
- It turns out that if we look at what went wrong, we can figure out a non-trivial factor of n

Here it is in action:

> :m Math.NumberTheory.Factor
> pfactors $ 10^30+5
[3,5,4723,1399606163,10085210079364883]
(0.55 secs, 210033624 bytes)
> pfactors $ 10^30+6
[2,7,3919,758405810021,24032284101871]
(2.31 secs, 883504748 bytes)

I love the way it can crunch through 12-digit prime factors with relative ease.