Our Imaginary Friends


Having hopefully now convinced you that floating point is the right choice for general purpose numerical computing, provided you put some thought into how you use it, we're ready to start implementing the library.
We shall start with some basic mathematical types and functions upon which we can build more advanced algorithms.
First up are the complex numbers. A complex number \(z = x + iy\) has both a real part \(x\) and an imaginary part \(y\) which is multiplied by \(i = \sqrt{-1}\). We can define arithmetic operations for complex numbers...

Whoa! Hold up there! I'm taking the square root of what now?!

Since multiplying a real number by itself always results in a non-negative number it might not be entirely obvious what taking the square root of a negative number actually means. I could be glib and simply assure you that you don't need to worry about it but I think I can argue from analogy that it's a perfectly reasonable concept.
Imagine that we've defined a symbol \(k\) to mean a multiple of a thousand. Clearly we don't want to mix up the \(k\) part and the unit part of a number so must keep track of them during arithmetic. For example
\[ (2k + 623) + (1k + 491) = (3k + 1104) \]
Now if whenever we we see a number greater than a thousand we replace it with \(k\) plus that number minus one thousand we have a rule for transforming it into a unique, unambiguous representation. In this case
\[ (3k + 1104) = (3k + k + 104) = (4k + 104) \]
If we similarly think of \(i\) as being nothing more than a symbol with the transformation rule that \(i^2 = -1\) we can unambiguously define rules of arithmetic for complex numbers. For example, addition is defined as
\[ z_0 + z_1 = (x_0 + i y_0) + (x_1 + i y_1) = (x_0 + x_1) + i (y_0 + y_1) \]
and multiplication as
\[ \begin{align*} z_0 \times z_1 &= (x_0 + i y_0) \times (x_1 + i y_1)\\ &= (x_0 \times x_1) + (x_0 \times i y_1) + (x_1 \times i y_0) + (i y_0 \times i y_1)\\ &= (x_0 x_1 - y_0 y_1) + i (x_0 y_1 + x_1 y_0) \end{align*} \]
The notion that \(i\) is just a symbol with a rule for transforming it is very much in keeping with the modern view of mathematics as the study of systems of symbol manipulation. Rather than worry about what it all actually means mathematicians are generally more concerned that the rules have interesting consequences and do not lead to any contradictions.

Now division is a slightly trickier operation that is easiest to define in terms of the product of the numerator and the reciprocal of the denominator. The latter is given by
\[ \frac{1}{x + i y} = \frac{x - i y}{x^2 + y^2} \]
as can be seen by multiplying it by \(x + i y\).
\[ \begin{align*} (x + i y) \times \frac{x - i y}{x^2 + y^2} &= \frac{(x + i y) \times (x - i y)}{x^2 + y^2}\\ &= \frac{x^2 - x \times i y + x \times i y - i^2 y^2}{x^2 + y^2}\\ &= \frac{x^2 + y^2}{x^2 + y^2}\\ &= 1 \end{align*} \]
Hence division is defined as
\[ \begin{align*} (x_0 + i y_0) \div (x_1 + i y_1) &= (x_0 + i y_0) \times \frac{x_1 - i y_1}{x_1^2 + y_1^2}\\ &= \frac{(x_0 + i y_0) \times (x_1 - i y_1)}{x_1^2 + y_1^2}\\ &= \frac{(x_0 x_1 - y_0 y_1) + i (x_0 y_1 + x_1 y_0)}{x_1^2 + y_1^2} \end{align*} \]
To define general mathematical operations we shall rely again upon our favourite mathematical tool, Taylor's theorem. To illustrate how, consider the exponential of a complex number
\[ e^{x + i y} = e^x \times e^{i y} \]
Now the \(e^x\) term can be calculated with Math.exp so we just need to work out the value of \(e^{i y}\). Plugging \(i y\) into the Taylor series expansion of the exponential function at zero, also known as the Maclaurin series, we have
\[ \begin{align*} e^{i y} &= 1 + i y + \tfrac{1}{2} (i y)^2 + \tfrac{1}{6} (i y)^3 + \tfrac{1}{24} (i y)^4 + \tfrac{1}{120} (i y)^5 + \dots + \tfrac{1}{n!} (i y)^n + \dots\\ &= 1 + i y - \tfrac{1}{2} y^2 - \tfrac{1}{6} i y^3 + \tfrac{1}{24} y^4 + \tfrac{1}{120} i y^5 + \dots + \tfrac{1}{n!} (i y)^n + \dots\\ &= \left(1 - \tfrac{1}{2} y^2 + \tfrac{1}{24} y^4 + \dots + \tfrac{(-1)^n}{(2n)!} y^{2n} + \dots\right) + \left(i y - \tfrac{1}{6} i y^3 + \tfrac{1}{120} i y^5 + \dots + \tfrac{(-1)^n}{(2n+1)!} i y^{2n+1} + \dots\right)\\ &= \left(1 - \tfrac{1}{2} y^2 + \tfrac{1}{24} y^4 + \dots + \tfrac{(-1)^n}{(2n)!} y^{2n} + \dots\right) + i \left(y - \tfrac{1}{6} y^3 + \tfrac{1}{120} y^5 + \dots + \tfrac{(-1)^n}{(2n+1)!} y^{2n+1} + \dots\right) \end{align*} \]
By gathering the real and imaginary terms separately we can see by inspection that they are the Maclaurin series' for the cosine and sine respectively.
\[ \begin{align*} \cos y &= 1 - \tfrac{1}{2}y^2 + \tfrac{1}{24}y^4 + \dots + \frac{(-1)^n}{(2n)!}y^{2n} + \dots\\ \sin y &= y - \tfrac{1}{6}y^3 + \tfrac{1}{120}y^5 + \dots + \frac{(-1)^n}{(2n+1)!}y^{2n+1} + \dots \end{align*} \]
All three of these series converge for any value of \(y\) so the exponential of a complex number is therefore given by
\[ e^{x + i y} = e^x \cos y + i e^x \sin y \]
This is known as Euler’s formula and, if we set \(x=0\) and \(y=\pi\), yields the surprising identity
\[ e^{i \pi} = -1 \]
linking negative numbers, transcendental irrationals and complex numbers; three extremely important developments in the concept of number.

The realisation that functions could be extended to complex numbers in this way lead to a flurry of activity analysing the properties of functions of complex numbers from which grew the unimaginatively titled subject of complex analysis. Some of its results provide suprising insights into the properties of functions of real numbers but, given that this is a blog about numerical computing and not an undergraduate course in mathematics, we'll not be delving into them.

Figure 1: Complex Numbers As Coordinates
The form of the exponential of a complex number is suggestive of an alternative representation for complex numbers; if we treat the real and imaginary parts as Cartesian coordinates \((x, y)\) then we can uniquely identify them with their distances to the origin, \((0, 0)\), known as the magnitude, \(r\), of the number, and the anticlockwise angle between the line from the origin and the positive \(x\)-axis, known as the argument, \(\theta\), as illustrated by figure 1.
\[ \begin{align*} x + i y &= r (\cos \theta + i \sin \theta) = r e^{i \theta}\\ r &= \sqrt{x^2 + y^2}\\ \theta &= \arctan \tfrac{y}{x} \end{align*} \]
In some sense this provides an interpretation of \(i\) beyond an arbitrary symbol; if we consider positive real numbers to represent forwards and negative real numbers to represent backwards then \(i\) represents left.

Now expressing complex numbers in this form makes some calculations a great deal easier. For example, to raise a complex number to a real power we have
\[ (x + i y)^n = \left(r e^{i \theta}\right)^n = r^n e^{i n \theta} = r^n \cos n \theta + i r^n \sin n \theta \]
Unfortunately, since the sine and cosine functions are periodic (adding \(2\pi\) to their arguments yields the same results) there are an infinite number of representations of this form for each complex number. The usual convention is to choose a specific \(2\pi\) wide range of angles from which we choose the value of the argument; we shall use \(-\pi<\theta\leqslant\pi\).
This can lead to some surprising consequences. For example, consider
\[ \left((-1)^3\right)^\frac{1}{3} \]
Expressing this in the magnitude-argument form we have
\[ \begin{align*} \left((e^{i\pi})^3\right)^\frac{1}{3} &= \left(e^{3i\pi}\right)^\frac{1}{3}\\ &= \left(e^{i\pi}\right)^\frac{1}{3}\\ &= e^{\frac{i\pi}{3}} \end{align*} \]
Putting this back into real and imaginary parts yields
\[ \cos \tfrac{\pi}{3} + i \sin \tfrac{\pi}{3} = 0.5 + \tfrac{\sqrt 3}{2} i \]
and not minus one as we might have expected. The reason for this is that when using complex numbers there are three cubic roots of minus one rather than the one we'd expect when using real numbers. Our calculation just happens to result in a different one than we started with.
We're quite used to this when dealing with square roots of real numbers; when taking the square root of the square of a negative number our calculators and computers will give a positive number. If the behaviour of roots of complex numbers seems surprising it is simply because we are less familiar with them.

A Complex Class

Rather than try to work out how to represent \(i\) as an object, we shall simply define complex numbers as ordered pairs of reals with rules of arithmetic that map onto those of complex numbers. For example
\[ (x_0, y_0) + (x_1, y_1) = (x_0 + x_1, y_0 + y_1)\\ (x_0, y_0) \times (x_1, y_1) = (x_0 x_1 - y_0 y_1, x_0 y_1 + x_1 y_0) \]
The definition of the ak.complex class is illustrated by listing 1.

Listing 1: ak.complex
ak.COMPLEX_T = 'ak.complex';

function Complex(){}
Complex.prototype = {TYPE: ak.COMPLEX_T, valueOf: function(){return ak.NaN;}};

ak.complex = function() {
 var c     = new Complex();
 var arg0  = arguments[0];
 var state = {re: 0, im: 0};

 constructors[ak.nativeType(arg0)](state, arg0, arguments);

 c.re = function() {return state.re;};
 c.im = function() {return state.im;};

 c.toString = function() {
  return '('+state.re.toString()+','+state.im.toString()+'i)';

 c.toExponential = function(d) {
  return '('+state.re.toExponential(d)+','+state.im.toExponential(d)+'i)';

 c.toFixed = function(d) {
  return '('+state.re.toFixed(d)+','+state.im.toFixed(d)+'i)';

 c.toPrecision = function(d) {
  return '('+state.re.toPrecision(d)+','+state.im.toPrecision(d)+'i)';

 return Object.freeze(c);

var constructors = {};

The ordered pair that we use to represent the complex number is here provided by the state object which is populated by a call to a member of the constructors object. Aside from this we have trivial accessors methods re and im for the real and imaginary parts of the number and methods for conversion to string. Being an arithmetic type implemented with floating point I felt it appropriate to implement equivalents of the string conversions supported by JavaScript's Number object.

ak.complex Constructors

Having a state of just two numbers, the constructors for ak.complex are relatively straightforward, as shown in listing 2.

Listing 2: ak.complex Constructors
constructors[ak.NUMBER_T] = function(state, re, args) {
 var arg1 = args[1];

 state.re = Number(re);
 constructors[ak.NUMBER_T][ak.nativeType(arg1)](state, arg1);

constructors[ak.NUMBER_T][ak.NUMBER_T] = function(state, im) {
 state.im = Number(im);

constructors[ak.NUMBER_T][ak.UNDEFINED_T] = function() {

constructors[ak.OBJECT_T] = function(state, obj) {
 state.re = ak.nativeType(obj.re)===ak.FUNCTION_T ? Number(obj.re()) : Number(obj.re);
 state.im = ak.nativeType(obj.im)===ak.FUNCTION_T ? Number(obj.im()) : Number(obj.im);

So if the object is constructed with a single number then that is assigned to its real part and the imaginary part is left at zero, whilst a pair of numbers set the real and imaginary parts in that order. Finally we have construction from an object which initialises the real and imaginary parts with the value of the re and im members, or to the results of them if they're methods.

ak.complex Overloads

As usual we shall implemented the bulk of the functionality with overloads, the simplest of which are given in listing 3.

Listing 3: ak.complex abs And neg
function abs(c) {
 return Math.sqrt(Math.pow(c.re(), 2) + Math.pow(c.im(), 2));
ak.overload(ak.abs, ak.COMPLEX_T, abs);

function neg(c) {
 return ak.complex(-c.re(), -c.im());
ak.overload(ak.neg, ak.COMPLEX_T, neg);

These functions do pretty much the same thing as their equivalents for numbers, neg negating both the real and imaginary parts and abs returning the distance of the complex number from zero when interpreted as Cartesian coordinates, known as the magnitude. For normal numbers the latter simply removes the sign, making the result positive, but for those that is the distance from zero.

For complex numbers we also have the argument, the angle they make with the positive \(x\)-axis when plotted on a graph, provided by a new overload arg. Another new overload is conj which returns the conjugate of a complex number in which the imaginary part is negated. Both are given in listing 4, together with overloads for real numbers.

Listing 4: ak.complex arg And conj
function arg(c) {
 return Math.atan2(c.im(), c.re());

function argR(r) {
 return r>=0 ? 0 : ak.PI;

ak.overload(ak.arg, ak.COMPLEX_T, arg);
ak.overload(ak.arg, ak.NUMBER_T,  argR);

function conj(c) {
 return ak.complex(c.re(), -c.im());

function conjR(r) {
 return r;

ak.overload(ak.conj, ak.COMPLEX_T, conj);
ak.overload(ak.conj, ak.NUMBER_T,  conjR);

Addition proceeds exactly as we have discussed, with mixed arithmetic treating real numbers as complex numbers with an imaginary part of zero, as illustrated by listing 5.

Listing 5: ak.complex add
function add(c0, c1) {
 return ak.complex(c0.re()+c1.re(), c0.im()+c1.im());

function addCR(c, r) {
 return ak.complex(r+c.re(), c.im());

function addRC(r, c) {
 return ak.complex(r+c.re(), c.im());

ak.overload(ak.add, [ak.COMPLEX_T, ak.COMPLEX_T], add);
ak.overload(ak.add, [ak.COMPLEX_T, ak.NUMBER_T],  addCR);
ak.overload(ak.add, [ak.NUMBER_T,  ak.COMPLEX_T], addRC);

Multiplication, division and exponentiation also follow the rules described above as shown, omitting mixed arithmetic overloads, in listing 6.

Listing 6: ak.complex mul, div And exp
function mul(c0, c1) {
 return ak.complex(c0.re()*c1.re()-c0.im()*c1.im(), c0.re()*c1.im()+c0.im()*c1.re());
ak.overload(ak.mul, [ak.COMPLEX_T, ak.COMPLEX_T], mul);

function div(c0, c1) {
 var abs2   = Math.pow(c1.re(), 2) + Math.pow(c1.im(), 2);
 var inv_r1 =  c1.re()/abs2;
 var inv_i1 = -c1.im()/abs2;

 return ak.complex(c0.re()*inv_r1-c0.im()*inv_i1, c0.re()*inv_i1+c0.im()*inv_r1);
ak.overload(ak.div, [ak.COMPLEX_T, ak.COMPLEX_T], div);

function exp(c) {
 var e = Math.exp(c.re());
 return ak.complex(e*Math.cos(c.im()), e*Math.sin(c.im()));
ak.overload(ak.exp, ak.COMPLEX_T, exp);

Now if we're using complex numbers it's convenient to redefine ak overloads that might have invalid results for negative reals but that are perfectly well defined for complex numbers. For example, the square root of a negative number should return a complex number rather than NaN if we have decided to use complex numbers. The square roots of complex and real numbers are given in listing 7.

Listing 7: ak.complex sqrt
function sqrt(c) {
 var r = Math.sqrt(ak.abs(c));
 var a = arg(c) / 2;
 return ak.complex(r*Math.cos(a), r*Math.sin(a));

function sqrtR(r) {
 return (r>=0) ? Math.sqrt(r) : ak.complex(0, Math.sqrt(-r));

ak.overload(ak.sqrt, ak.COMPLEX_T, sqrt);
ak.overload(ak.sqrt, ak.NUMBER_T,  sqrtR);

The full complement of overloads of existing operations, new overloads and replacements for those that are well defined for complex arithmetic are given in listing 8.

Listing 8: ak.complex Overloads
if(!ak.arg)  ak.arg  = function(x) {return ak.arg[ak.type(x)](x)};
if(!ak.conj) ak.conj = function(x) {return ak.conj[ak.type(x)](x)};

ak.overload(ak.abs,  ak.COMPLEX_T, abs);
ak.overload(ak.arg,  ak.COMPLEX_T, arg);
ak.overload(ak.arg,  ak.NUMBER_T,  argR);
ak.overload(ak.conj, ak.COMPLEX_T, conj);
ak.overload(ak.conj, ak.NUMBER_T,  conjR);
ak.overload(ak.cos,  ak.COMPLEX_T, cos);
ak.overload(ak.cosh, ak.COMPLEX_T, cosh);
ak.overload(ak.exp,  ak.COMPLEX_T, exp);
ak.overload(ak.inv,  ak.COMPLEX_T, inv);
ak.overload(ak.log,  ak.COMPLEX_T, log);
ak.overload(ak.log,  ak.NUMBER_T,  logR);
ak.overload(ak.neg,  ak.COMPLEX_T, neg);
ak.overload(ak.sin,  ak.COMPLEX_T, sin);
ak.overload(ak.sinh, ak.COMPLEX_T, sinh);
ak.overload(ak.sqrt, ak.COMPLEX_T, sqrt);
ak.overload(ak.sqrt, ak.NUMBER_T,  sqrtR);
ak.overload(ak.tan,  ak.COMPLEX_T, tan);
ak.overload(ak.tanh, ak.COMPLEX_T, tanh);

ak.overload(ak.add,  [ak.COMPLEX_T, ak.COMPLEX_T], add);
ak.overload(ak.add,  [ak.COMPLEX_T, ak.NUMBER_T],  addCR);
ak.overload(ak.add,  [ak.NUMBER_T,  ak.COMPLEX_T], addRC);
ak.overload(ak.dist, [ak.COMPLEX_T, ak.COMPLEX_T], dist);
ak.overload(ak.dist, [ak.COMPLEX_T, ak.NUMBER_T],  distCR);
ak.overload(ak.dist, [ak.NUMBER_T,  ak.COMPLEX_T], distRC);
ak.overload(ak.div,  [ak.COMPLEX_T, ak.COMPLEX_T], div);
ak.overload(ak.div,  [ak.COMPLEX_T, ak.NUMBER_T],  divCR);
ak.overload(ak.div,  [ak.NUMBER_T,  ak.COMPLEX_T], divRC);
ak.overload(ak.eq,   [ak.COMPLEX_T, ak.COMPLEX_T], eq);
ak.overload(ak.eq,   [ak.COMPLEX_T, ak.NUMBER_T],  eqCR);
ak.overload(ak.eq,   [ak.NUMBER_T,  ak.COMPLEX_T], eqRC);
ak.overload(ak.mul,  [ak.COMPLEX_T, ak.COMPLEX_T], mul);
ak.overload(ak.mul,  [ak.COMPLEX_T, ak.NUMBER_T],  mulCR);
ak.overload(ak.mul,  [ak.NUMBER_T,  ak.COMPLEX_T], mulRC);
ak.overload(ak.ne,   [ak.COMPLEX_T, ak.COMPLEX_T], ne);
ak.overload(ak.ne,   [ak.COMPLEX_T, ak.NUMBER_T],  neCR);
ak.overload(ak.ne,   [ak.NUMBER_T,  ak.COMPLEX_T], neRC);
ak.overload(ak.pow,  [ak.COMPLEX_T, ak.COMPLEX_T], pow);
ak.overload(ak.pow,  [ak.COMPLEX_T, ak.NUMBER_T],  powCR);
ak.overload(ak.pow,  [ak.NUMBER_T,  ak.COMPLEX_T], powRC);
ak.overload(ak.pow,  [ak.NUMBER_T,  ak.NUMBER_T],  powRR);
ak.overload(ak.sub,  [ak.COMPLEX_T, ak.COMPLEX_T], sub);
ak.overload(ak.sub,  [ak.COMPLEX_T, ak.NUMBER_T],  subCR);
ak.overload(ak.sub,  [ak.NUMBER_T,  ak.COMPLEX_T], subRC);

Note that we must check that the dispatch functions for arg and conj haven't been defined by another type before we add them, unlike those that we know for certain have been defined in the core ak library.
The implementations of these overloads can be found in Complex.js

Using ak.complex

Program 1 illustrates some of the basic arithmetic operations described above.

Program 1: Complex Arithmetic

Plotting the results of a function of complex numbers is a bit of a tricky job since we need to convey two pieces of information; either the real and imaginary parts or the magnitude and argument. The latter are generally easier to deal with, as illustrated by program 2.

Program 2: Complex Square Root

Here we have two new chart objects, both attached to the same canvas. The first is a contour plot which represents the result of a function of two arguments with the colour of the pixel plotted at the point they represent. The second is a vector plot which represents starting points with dots and connects them to end points with lines.
The colours used to represent the values in the contour plot are set by the palette method which expects the number of levels as the first argument and an array of colour levels to interpolate between from the lowest, 0, to the highest, 1. Note that the bounds must be three, rather than two, dimensional arrays so that the chart can choose the correct colour for the value at each point.
The vector plot requires an array of ak.ui.pointer objects rather than the points that are expected by simple graphs. These expect the starting point as an array for their first argument and the offsets to the end point as an array for the second.

In program 2 the contour plot represents the magnitude of the result of the function, with darker hues indicating smaller numbers, and the vector plot just their arguments, by giving the pointers constant length.
As we should expect, those values with large negative real parts and small imaginary parts point in the direction of \(\pm i\) and those with large negative real parts and small imaginary parts point in the direction of +1.

Now this is all well and good, but let's face it, there's only one kind of graph of complex numbers that we're really interested in...

Fractals Baby!

Figure 2: What On Earth Was I Thinking?
Beloved of undergraduate mathematicians the world over, fractals hit mainstream popular culture in the late 1980s. In fact, I still have a T-shirt with a black and white print of a fractal.

Made from hemp.

As is tragically shown in figure 2.

Perhaps the most famous fractal is the Mandelbrot set, although I suspect that the majority of my fellow unwashed hippies whose torsos it adorned had no idea how it was defined.
The truly remarkable thing about it is how simple it actually is. At its heart is the quadratic equation
\[ f_c(z) = z^2 + c \]
Now, if for some complex \(c\) and \(z=0\), we repeatedly set \(z = f_c(z)\) and \(z\) never flies off to infinity then \(c\) is a member of the Mandelbrot set. For example, the real value minus one is a member since
\[ \begin{align*} f_{-1}(0) &= \phantom{-}0^2 + -1 = -1\\ f_{-1}(-1) &= -1^2 + -1 = \phantom{-}0 \end{align*} \]
whereas plus one isn't because
\[ \begin{align*} f_{1}(0) &= 0^2 + 1 = 1\\ f_{1}(1) &= 1^2 + 1 = 2\\ f_{1}(2) &= 2^2 + 1 = 5\\ f_{1}(5) &= 5^2 + 1 = 26 \end{align*} \]
and so on.
In general it's quite tricky to actually prove that a complex number is in the Mandelbrot set, so the usual tactic is to approximate it with the set of numbers that don't lead to \(z\) growing above some given magnitude within some given number of assignments of \(z = f_c(z)\). Rather than simply plot the numbers that are in the Mandelbrot set program 3 plots the number of iterations it takes for \(z\) to grow above that magnitude for that familiar psychedelic image.

Program 3: The Mandelbrot Set


Leave a comment