maximus:book.git
3 years agoDecoupling rendering and colouring.
Claude Heiland-Allen [Sat, 21 Dec 2013 03:10:43 +0000 (03:10 +0000)]
Decoupling rendering and colouring.

Unfortunately the colour scheme we chose looked quite bad, and it took
a few minutes to render all the images.  To reduce waiting time, we can
save the raw data from rendering to a file, and colour it later.

We rename our mandelbrot.c program to render.c, and strip out the colour
type handling.  We delete the render_grey() from mandelbrot_imp.c, and
rename render_colour() back to render(), using a colour mapping that
preserves information.  With 3 bytes, each with 8 bits, we can represent
24bit integers, enough range to represent over 16 million iterations.
The >> operator is bitwise shift right, here we use it to extract the
individual bytes from the escape time.  We rely on putchar() truncating
to the least significant byte in the resulting integer.

We write a new program to remap the colours in the PPM output from the
render program.  Our colour.c reads the image header from its standard
input using scanf(), which works a bit like printf() in reverse.  We
check that we read the header successfully by confirming that we parsed
two items, the width and height, and that the next character after the
maxval of 255 is a newline.  Then we loop over all the pixels, repacking
the RGB values into a single integer.

We now use the HSV colour space, allowing us to specify colours by hue,
saturation and value.  Hue is circular, from red at 0, with increasing
values running through the rainbow yellow green blue violet, and then
back through purple and pink to red again at 1, where the cycle repeats.
Saturation blends colour intensity from greyscale at 0 to full blast at
1, and value controls brightness: black at 0 and full strength at 1.
We use single precision float for the calculations: 24bits is more than
enough, given each colour channel has only 8bits.  We need to scale the
output values from hsv2rgb() from [0..1] to [0..255], and we clamp them
to the output range to make sure no overflow glitches occur.

We add the new rules to the Makefile, and the new binaries to the
.gitignore file, and rebuild to check everything works:

    $ make
    $ time for exp in $(seq -w 0 15)
      do
        ./render \
          -1.759937295271429505091916757 \
           0.012591790526496335594898047 \
           2e-${exp} 2 > ${exp}-render.ppm
      done
    real    3m51.110s
    user    3m50.706s
    sys     0m0.288s
    $ time for i in ??-render.ppm
      do
        ./colour < ${i} > ${i}-colour.ppm
      done
    real    0m2.190s
    user    0m2.080s
    sys     0m0.072s
    $

As we suspected, rendering takes a lot longer than colouring, and having
separated the two processes we no longer have to rerender when we adjust
the colour scheme.

3 years agoColour images.
Claude Heiland-Allen [Sat, 21 Dec 2013 01:32:05 +0000 (01:32 +0000)]
Colour images.

So far our images have progressed from black and white to greyscale,
with the grey level based on the escape time.  We're writing PGM files,
part of the PNM family of formats, and this family include PPM which
supports colour.  The text header for PPM is very similar to PGM, with
the magic number P6 in place of PGM's P5.  With the maxval of 255 that
we are using, each pixel has 3 bytes, for the red, green and blue colour
channels (in that order, tightly packed).

We rename our existing render to render_grey, and write a new
render_colour to output the modified header and pixel data.  Between
calculating the escape time and outputting the pixel data, we need
to turn the linear value into a point in 3D colour space.  About the
simplest way to do this is to vary the speed at which the output colour
channels change with respect to the input value.  We define a colour()
function in the mandelbrot.h header because it doesn't depend on the
floating point type used for calculations, and make the red channel
change more quickly and the blue channel chnage more slowly.

In the main program we handle another command line argument, to switch
between greyscale rendering and colour rendering.  With 3 number types
and 2 colour types, we have ended up with 6 render functions to choose
from in our now-nested switch statements.  We increase the maximum
iteration limit to show more detail near the main body of the Mandelbrot
set, and switch the default number type to double because it is both
faster and more precise than float (though we still need to choose the
float type argument if we want to specify a colour type).

Now let's render a sequence of images zooming in towards a miniature
copy within the Mandelbrot set, to compare colour with greyscale:

    $ time for exp in $(seq -w 0 15)
      do
        ./mandelbrot \
          -1.759937295271429505091916757 0.012591790526496335594898047 \
          2e-${exp} 2 0 > ${exp}-grey.pgm
        ./mandelbrot \
          -1.759937295271429505091916757 0.012591790526496335594898047 \
          2e-${exp} 2 1 > ${exp}-colour.ppm
      done
    real    7m41.621s
    user    7m40.925s
    sys     0m0.492s

This is a bash shell script, which loops over the sequence of numbers from
0 to 15, and each pass through the loop it runs our program twice, using the
number as radius exponent, once for greyscale and once for colour.  We end
up with 32 images to browse through.

3 years agoGeneric programming with preprocessor macros.
Claude Heiland-Allen [Sat, 21 Dec 2013 00:41:12 +0000 (00:41 +0000)]
Generic programming with preprocessor macros.

In supporting double precision to allow deeper zooming without pixelation,
we duplicated a lot of code, just to change the number type used.  Now we'll
get rid of the code duplication, leaving us with one implementation of each
algorithm.  We move one copy of the implementation into mandelbrot_imp.c,
using the name FTYPE wherever ``float'' or ``double'' is needed.  We also
wrap each function that has variants for different number types in FNAME().

Now, FTYPE and FNAME() are not defined by the C standard, so the code won't
compile.  We'll define them ourselves, in a header file mandelbrot.h.  This
begins and ends with an idiomatic header guard, to avoid errors if the file
would be included more than once.  Next it includes the library headers that
mandelbrot_imp.c needs.  We to define FTYPE to be a floating point type,
like float, whose variant function suffix is f.  Then we include the
implementation code.  Now we can undefine FTYPE and FNAME(), and repeat this
stanza for the other floating types available in C: double and long double.

This technique uses the C preprocessor, which the compiler runs early on in
the compilation pipeline (first the preprocessor performs file inclusion and
macro expansion, then the C compiler generates assembly source for the
target architecture, which the assembler translates to machine code, which
is then linked into a program).  Macro expansion is similar to search and
replace within a text editor, with macros able to take parameters.
The ## operator in a macro joins two tokens together to form a new token.

Now our main program file mandelbrot.c includes our header file, and its
main() parses the command line arguments and calls the appropriate function
from our mini-library.  We use strtold() now instead of atof() to avoid
losing precision before we even begin.  The switch() statement is an
alternative to nested if statements.  Finally we add the extra files to the
prerequisites of our Makefile target, so that ``make'' will rebuild the
program when we update any of them.

Now we can benchmark the implementation with different number types:

    $ time ./mandelbrot 0 0 2 0 > float.pgm
    real    0m0.296s
    user    0m0.288s
    sys     0m0.004s
    $ time ./mandelbrot 0 0 2 1 > double.pgm
    real    0m0.130s
    user    0m0.128s
    sys     0m0.000s
    $ time ./mandelbrot 0 0 2 2 > long-double.pgm
    real    0m0.525s
    user    0m0.528s
    sys     0m0.000s
    $

These results are surprising: float is actually slower than double,
despite having less precision.

3 years agoDouble precision.
Claude Heiland-Allen [Fri, 20 Dec 2013 23:34:25 +0000 (23:34 +0000)]
Double precision.

We've been using C's floating point number type ``float''.  Floating
point numbers have a mantissa and an exponent, similar to scientific
notation for decimal numbers like 1.234e-5.  The mantissa has a fixed
width, which means there are only so many significant figures that
can be represented exactly.  The ``float'' type is precise to about
7 decimal digits, having 24 bits (to convert from bits to digits we
multiply by \(log 2 / log 10\), and the exponent says where the
significant bits begin (leading 0s are not significant).

Recall our last pixellated image:

    $ ./mandelbrot -0.8031  0.1780  0.000001 > mandelbrot.pgm

Using the device coordinate transformation, the distance between
adjacent pixels in the image is 0.000001/(720/2) = 2.78e-9.  Comparing
this with the image center shows us why the image is pixelated:
-0.8031 = -0.8031000 to 7 significant figures (working in decimal).
The neighbouring representable numbers with 7 significant figures are
-0.8031001 and -0.8030999.  Each of these differs from the center
coordinate by 0.0000001 = 1e-7, which is 36 times larger the pixel
spacing we computed as 2.78e-9.  This means that for every 36 pixels
in the horizontal direction, there's only one ``float'' value that
can be used: not surprisingly, if 36 pixels use the same input value,
they'll give the same output value.

The situation is better in the vertical direction: 0.1780 is around
4.5 times smaller than -0.8031.  The nature of binary floating point
means representable values are more closely spaced the smaller the
absolute value: if you halve the absolute value, the spacing halves too.
For example with 7 significant figures in decimal, the unit of least
precision for 1.000000 is 0.000001 = 1e-6, and the unit of least
precision for 1000000 = 1e6 = 1.

Now we know why the pixelation occurs (we don't have enough precision,
our mantissa isn't wide enough), and why it's so much worse in one
direction than the other (floating point number precision is relative
to the absolute size of the value).

C has more than one floating point number type, and the ``float''
we have been using is a ``single precision'' type.  If we switch to
the ``double precision'' type, imaginatively called ``double'', then
our pixellation problem will be solved, at least for the moment.  And
sure enough, after replicating our functions to make them use double
precision, we get a much clearer image:

    $ ./mandelbrot -0.8031  0.1780  0.000001 0 > float.pgm
    $ ./mandelbrot -0.8031  0.1780  0.000001 1 > double.pgm
    $ display float.pgm double.pgm

We kept the single precision floating point version and renamed it
with an ``f'' suffix, as is the idiom in C, and use atoi() to read
the 4th command line argument as a flag to tell our program which
implementation to use.

3 years agoRefactoring.
Claude Heiland-Allen [Fri, 20 Dec 2013 23:15:02 +0000 (23:15 +0000)]
Refactoring.

The last two images looked somehow off -- grainy and pixelated.  We'll
investigate this soon, but first let's reorganize our program, breaking
it down into smaller pieces.  So far the work is all performed in the
main() function, which loops over all the pixels, transforms the device
coordinates, and calculates the escape time for the point.  These three
tasks are unrelated, so we split each into its own function, respectively
render(), coordinate(), and calculate().  We define render() after the
other two, because C needs to have declarations before use, and render()
uses both coordinate() and calculate().  This refactoring should make it
easier to adapt and extend our program in the future.

3 years agoEscape time colouring.
Claude Heiland-Allen [Fri, 20 Dec 2013 21:45:58 +0000 (21:45 +0000)]
Escape time colouring.

So far we've been colouring in a binary fashion: points possibly maybe
inside the Mandelbrot set are black, and points definitely outside are
white.  We're also breaking out of the loop early as soon as we know
a point has escaped.  We can add the information about ``how many
iterations it took before we knew the point escaped'' to our images
instead of throwing it away.  This technique is called ``escape time
colouring''.

First we need a variable to hold the final iteration count.  We initialize
it to zero, because this is really an impossible value: z_0 is always 0,
and 0 is always less than 2.  We correct our loop counter initialization:
within our loop we perform the iteration step before we check whether the
iterate escaped, so the n we need at that point starts from 1.  If the
point escaped we store the final n before breaking out of the loop.
After the loop finishes, we can output the final n directly.

Now the images look even more interesting: the isolated speckles with
copies of the whole seem to be connected to the main body by intricate
branches and spirals:

    $ ./mandelbrot -0.75    0       1.25     > mandelbrot.pgm
    $ ./mandelbrot -1.76    0       0.05     > mandelbrot.pgm
    $ ./mandelbrot -1.765   0.023   0.01     > mandelbrot.pgm
    $ ./mandelbrot -1.76648 0.04173 0.0005   > mandelbrot.pgm
    $ ./mandelbrot -0.8031  0.1780  0.1      > mandelbrot.pgm
    $ ./mandelbrot -0.8031  0.1780  0.01     > mandelbrot.pgm
    $ ./mandelbrot -0.8031  0.1780  0.001    > mandelbrot.pgm
    $ ./mandelbrot -0.8031  0.1780  0.0001   > mandelbrot.pgm
    $ ./mandelbrot -0.8031  0.1780  0.00001  > mandelbrot.pgm
    $ ./mandelbrot -0.8031  0.1780  0.000001 > mandelbrot.pgm

3 years agoAvoiding square roots.
Claude Heiland-Allen [Fri, 20 Dec 2013 21:29:39 +0000 (21:29 +0000)]
Avoiding square roots.

Our program is still on the slow side.  The cabs() function we call
to find the absolute value of each iterate is performing a square root:

    \[ |x + i y| = \sqrt{x^2 + y^2} \]

This square root is unnecessary, because if we know both sides
of a comparison are non-negative, we can square both sides and
get the same result.  Moreover, because the escape radius doesn't
change, we can lift its squaring outside all the loops.

Because we'll use this square root avoidance technique more than once,
we define a function once, then call it by name.  Our cabs2() needs to
extract the real and imaginary parts, which we can do using the functions
creal() and cimag().

Now our program runs twice as fast as the previous version:

    $ time ./mandelbrot -0.75 0 1.25 > mandelbrot.pgm
    real    0m0.864s
    user    0m0.864s
    sys     0m0.000s

3 years agoEscaping a loop early: ``break''.
Claude Heiland-Allen [Fri, 20 Dec 2013 21:15:57 +0000 (21:15 +0000)]
Escaping a loop early: ``break''.

When we increased the resolution, our program took noticeably longer to
render each image.  We can speed it up.  Our program currently keeps
iterating to the maximum iterations limit for every pixel, only then
checking if it didn't exceed the escape radius.  But if any iterate along
the way exceeded the escape radius, we know that it's already well on the
way to infinity.  If we check each time through the loop, then we can
stop looping, with a ``break'' statement -- it makes the control flow
pass immediately to the statement following the loop body.

This simple optimization makes our program four times faster for the
initial view:

    $ time ./mandelbrot -0.75 0 1.25 > mandelbrot.pgm
    real    0m6.276s
    user    0m6.260s
    sys     0m0.016s
    $ time ./mandelbrot -0.75 0 1.25 > mandelbrot.pgm
    real    0m1.593s
    user    0m1.584s
    sys     0m0.004s
    $

3 years agoImage output.
Claude Heiland-Allen [Fri, 20 Dec 2013 20:37:50 +0000 (20:37 +0000)]
Image output.

The last image had curious speckles near the bottom, but our text-based
image resolution is too small to make them out properly.  It's also upside
down: mathematical convention dictates that complex numbers are visualized
with the imaginary part increasing in the upwards direction.  We can fix
the flip by negating the imaginary part when transforming from device
coordinates.

We want to create images with pixels instead of characters.  There are lots
of image formats out there, but one of the simplest is a family of related
formats called PNM, the ``portable anymap file format''.  We'll use the
portable graymap format PGM.  PGM files start with a text header, containing
metadata about the image.  First are the two characters P5 meaning ``this is
a PGM image'', followed by a blank space (typically a newline character).
Then come the width and height of the image, as decimal text, separated by
spaces.  The header ends with the ``maxval'' as decimal text, which is the
value that will represent white (zero represents black), followed by a single
newline character.  The header is then followed by the image data, in rows
from top to bottom with each row from left to right.  With a maxval of 255,
each pixel takes one byte.

We can increase the size to something a lot bigger than we had before, no
longer limited by our command line terminal width.  We no longer need to
adjust the aspect ratio because PNM formats have square pixels.  We can
use printf() for the image header text.  Because C treats bytes and
characters interchangably we can still use putchar(), but now we use
0 and 255 for a black and white image.  PNM formats have no separator between
rows (there is no need because the width is fixed by the header), so we
remove the one we output for text images.

Now our program writes image data, but our terminal expects text.  We
should ``redirect'' the output from our program to a file, so that we can
open it with an image viewer (such as ``display'' from the ImageMagick suite).
Command line shell redirection uses the ``>'' character to send the output
to a file (overwriting any pre-existing file, so be careful).  PGM files
usually have a .pgm extension.

    $ ./mandelbrot -0.75 0 1.25 > mandelbrot.pgm
    $ display mandelbrot.pgm

What's that blob on the left?

    $ ./mandelbrot -1.75 0 0.125 > mandelbrot.pgm
    $ display mandelbrot.pgm
    $ ./mandelbrot -1.76 0 0.03 > mandelbrot.pgm
    $ display mandelbrot.pgm

It looks like our original image of the whole Mandelbrot set.  What about
those speckles from before?

    $ ./mandelbrot -0.1 0.85 0.25 > mandelbrot.pgm
    $ display mandelbrot.pgm

They also look like (shrunken, rotated) copies of our original image of the
whole Mandelbrot set.

3 years agoConverting strings to numbers.
Claude Heiland-Allen [Fri, 20 Dec 2013 20:02:26 +0000 (20:02 +0000)]
Converting strings to numbers.

The region of interest in the image we have so far looks a bit small
and to the left of the frame.  We could edit the program to change the
center and radius and recompile it, but that would tire quickly.  Instead
we'll use the command line arguments.  The arguments are passed to main()
as strings, but we need numbers.  The C library header ``stdlib.h''
declares a function ``atof()'' to read a representation of a number in
decimal format from a string and convert it to a floating point number.

However, we should first check that we have been given arguments before
trying to use them.  If there is no string where we try to read a string,
we get ``undefined behaviour'': anything from a program crash (if we're
lucky) to carrying on with incorrect data and generating meaningless
output.  We want three floating point numbers from three arguments, for
the real and imaginary parts of the center, and the radius of the view;
we get the program name as the first argument, so we need to be sure that
we have more than three arguments.  In case fewer arguments were provided,
we leave the center and radius unchanged at the default values we set
previously.

Now we can zoom in to see more of the detail of the Mandelbrot set:

    $ ./mandelbrot -0.75 0.00 1.25
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ----------------------------------------------xxx-----------------------
    ----------------------------------------------xxx-----------------------
    ----------------------------------------xxxxxxxxxxxxx-xxx---------------
    ----------------------------------------xxxxxxxxxxxxxxx-----------------
    --------------------------------------xxxxxxxxxxxxxxxxxxx---------------
    -----------------------------xxxxx---xxxxxxxxxxxxxxxxxxxx---------------
    ----------------------------xxxxxxxx-xxxxxxxxxxxxxxxxxxx----------------
    --------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx------------------
    ----------------------------xxxxxxxx-xxxxxxxxxxxxxxxxxxx----------------
    -----------------------------xxxxx---xxxxxxxxxxxxxxxxxxxx---------------
    --------------------------------------xxxxxxxxxxxxxxxxxxx---------------
    ----------------------------------------xxxxxxxxxxxxxxx-----------------
    ----------------------------------------xxxxxxxxxxxxx-xxx---------------
    ----------------------------------------------xxx-----------------------
    ----------------------------------------------xxx-----------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    $ ./mandelbrot -1.50 0.00 0.50
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ----------------------------------------------------------x------------x
    ---------------------------------------------------xxxxxxxxxxxxx-------x
    --------------------------------------------------xxxxxxxxxxxxxxxxx---xx
    -------------------------------------------------xxxxxxxxxxxxxxxxxxx--xx
    ------------------------------------------------xxxxxxxxxxxxxxxxxxxxx-xx
    -------------------------------------------xxxx-xxxxxxxxxxxxxxxxxxxxx-xx
    --------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx
    -------------------------------------------xxxx-xxxxxxxxxxxxxxxxxxxxx-xx
    ------------------------------------------------xxxxxxxxxxxxxxxxxxxxx-xx
    -------------------------------------------------xxxxxxxxxxxxxxxxxxx--xx
    --------------------------------------------------xxxxxxxxxxxxxxxxx---xx
    ---------------------------------------------------xxxxxxxxxxxxx-------x
    ----------------------------------------------------------x------------x
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    $ ./mandelbrot -0.10 0.85 0.25
    xxx---------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx---------------
    ------------xxx---xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----xx---------------
    -------------------x--x---x-xxxxxxxxxxxxxxx-x-x---x---------------------
    -----------------------------xxxxxxxxxxx--------------------------------
    ---------------------------xxxxxxxxxxxxxx-------------------------------
    --------------------------xxxxxxxxxxxxxxxxx-----------------------------
    --------------------------xxxxxxxxxxxxxxxxx-----------------------------
    -------------------------xxxxxxxxxxxxxxxxxx-----------------------------
    --------------------------xxxxxxxxxxxxxxxxx-----------------------------
    --------------------------xxxxxxxxxxxxxxx-------------------------------
    ------------------------------xxxxxxxxx---------------------------------
    ---------------------------------xxxx-----------------------------------
    ---------------------------------xxxx-----------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    -----------------------------------------x------------------------------
    ------------------------------------------------------------------------
    -------------------------------x----------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    $

3 years agoThe Mandelbrot set.
Claude Heiland-Allen [Fri, 20 Dec 2013 19:22:59 +0000 (19:22 +0000)]
The Mandelbrot set.

Pick a complex number \(c\).  Starting from \(z_0 = 0\), form a sequence
of ``iterates'' by repeating \(z_{n+1} = z_n^2 + c\).  For some values of
\(c\), the iterates \(z_n\) get bigger and bigger, escaping towards
infinity.  These iterates are unbounded: for any size you pick, there'll
be some later iterate that is bigger still.  For other values of \(c\),
the iterates remain bounded.  The Mandelbrot set is a map of the parameter
plane, telling us which \(c\) values will give unbounded iterates, and
which remain bounded.  Mathematically, the Mandelbrot set consists of the
\(c\) values that remain bounded.  What does it look like?

It's possible to show that if an iterate gets bigger than 2, then it'll
escape to infinity.  So the whole Mandelbrot set must be inside a circle
\(|c| <= 2\).  One way to check that a point is not in the Mandelbrot set
is to keep iterating and see if it escapes.  But for any particular \(c\)
value, we don't yet know if it is in the Mandelbrot set or not -- if it
turns out that it is in the set, we would keep iterating forever.  So
that we don't have to wait beyond the end of time to see our picture of
the Mandelbrot set, we choose to give up after a certain maximum iteration
count.  If the point still hasn't escaped, it might or might not be in
the Mandelbrot set.  But if the point escaped, we can be sure the point
is not in the Mandelbrot set.  We can draw these two cases with two
different characters, giving us our first approximation to the Mandelbrot
set.

    $ ./mandelbrot
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ----------------------------------xx------------------------------------
    -------------------------------xxxxxxxxx--------------------------------
    -----------------------------xxxxxxxxxxx--------------------------------
    -----------------------xxxxxxxxxxxxxxxxxx-------------------------------
    --------------xxxxxxxxxxxxxxxxxxxxxxxxx---------------------------------
    -----------------------xxxxxxxxxxxxxxxxxx-------------------------------
    -----------------------------xxxxxxxxxxx--------------------------------
    -------------------------------xxxxxxxxx--------------------------------
    ----------------------------------xx------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    $

Mathematical variables are immutable.  In maths, \(z = z^2 + c\) is an
equation, and you solve it by finding a value for \(z\) that makes the
value on both sides the same.  For example, for \(c = 0\) there are two
solutions: \(z = 0\) and \(z = 1\).  In the C programming language, variables
can change as the program runs.  We already saw this with the loop counters,
in our use of the increment operator {\tt ++i}.  In C, {\tt =} is an
``assignment'' operator, not mathematical equality (for which C uses
{\tt ==}).  Variables in C are like boxes that contain values.  On the right
hand side of {\tt =}, values are taken out of the boxes and used to calculate
the result of the expression.  Then {\tt =} takes that result and stores it
into the box on the left hand side, overwriting the previously stored value.
The use of z as an ``l-value'' on the left hand side of {\tt z = z * z + c}
has a different meaning to the usees of z as an ``r-value'' on the right
hand side.

3 years agoComplex numbers.
Claude Heiland-Allen [Fri, 20 Dec 2013 16:39:44 +0000 (16:39 +0000)]
Complex numbers.

Different kinds of numbers can be thought of as a tower.
Counting gives us the whole numbers, then the desire for
an identity element for addition gives us zero:

    \[ 0 + x = x + 0 = x \]

Inverting addition gives subtraction, which needs negative
numbers:

    \[ x + (-x) = (-x) + x = 0 \]

Repeated addition gives multiplication, with identity 1,
and inverting multiplication gives division which needs
fractions:

    \[ n \times (1 / n) = (1 / n) \times n = 1 \]

Repeated multiplication gives powers, like \(n^2 = n \times n\),
and fractional powers give irrational numbers like
\(2^\frac{1}{2} = \sqrt{2} = 1.41421\ldots \) which can't
be represented by a fraction.

Taking the square root of a negative number can't give
any of the kinds of numbers we've seen so far, because
it's easy to work out that a negative number multiplied
by a negative number gives a positive number, and so
does a positive number multiplied by a positive number.
So we make up a new imaginary number \(i\) so that

    \[ i \times i = -1 \]

and then we have complex numbers with two parts: the sum
of a real part and an imaginary part.  The rules for
complex number arithmetic follow from the definition of \(i\)
and the rules for regular arithmetic (which is associative,
commutative, and distributive):

    \[ (a + b) + c = a + (b + c)
       (a \times b) \times c = a \times (b \times c) \]
    \[ a + b = b + a
       a \times b = b \times a\]
    \[ a \times (b + c) = (a \times b) + (a \times c) \]

Happily the C99 language has library support for complex
numbers in ``complex.h'', so we can include it and then
use complex floating point types, using {\tt I} for \(i\),
and regular arithmetic operators.  We now need to link
with the C maths library (``libm'') which contains the
implementations of the functions declared in the header,
so we add -lm to our Makefile rule.  Now we can draw our
circle \(|c|<r\) using the ``cabs()'' function.

The whole numbers without zero are written \(\mathbb{N}^+\), whole
numbers with zero \(\mathbb{N}\) (the N is for natural). The whole
numbers with negative numbers (the integers) are written \(\mathbb{Z}\)
(the Z is for Zahlen, the German word for numbers).  Fractions (the
rational numbers) are written \(\mathbb{Q}\) (the Q is for quotient, the
result of division).

Adding roots of polynomials gives the algebraic numbers, and
(a subset of) the complex numbers.  It's possible to show that
there are still more numbers, for example \(\pi\) and \(e\) (the
base of natural logarithms) are transcendental: eventually you
reach the continuum of the real number line written \(\mathbb{R}\),
and the complex plane written \(\mathbb{C}\).

3 years agoAspect ratio.
Claude Heiland-Allen [Fri, 20 Dec 2013 12:05:05 +0000 (12:05 +0000)]
Aspect ratio.

The previous program's output failed to look like a circle
for two reasons.  First, the width and height are unequal,
so scaling each device coordinate separately to [-1,1]
will distort the shape.  We correct this distortion by
choosing one of the axes to scale to [-1,1], we'll choose
the vertical direction becasue our image is wider than it
is tall and we want to have the whole circle in the image.
Dividing both coordinates by the same height/2.0 will make
the y axis range [-1,1] and the x axis range [-k,k] with
k > 1.

The second reason our circle doesn't look circular is that
the characters displayed are not square.  This does depend
on the particular font used, but common terminal fonts have
characters roughly half as wide as they are tall.  This makes
the spacing between characters is half as much in the horizontal
direction than the vertical direction.  We correct the spacing
by multiplying by the aspect ratio.

Now our circle looks much more circular:

    $ ./mandelbrot
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    ----------------------------xxxxxxxxxxxxxxxxx---------------------------
    -------------------------xxxxxxxxxxxxxxxxxxxxxxx------------------------
    ----------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxx---------------------
    ---------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx--------------------
    -------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx------------------
    ------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----------------
    ------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    ------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----------------
    ------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----------------
    -------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx------------------
    ---------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx--------------------
    ----------------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxx---------------------
    -------------------------xxxxxxxxxxxxxxxxxxxxxxx------------------------
    ----------------------------xxxxxxxxxxxxxxxxx---------------------------
    ------------------------------------------------------------------------
    $

3 years agoUser coordinates.
Claude Heiland-Allen [Fri, 20 Dec 2013 11:28:31 +0000 (11:28 +0000)]
User coordinates.

Device coordinates are integer pairs in \([0,width)\times[0,height)\).
We now transform the device coordinates to ``user coordinates'' giving
us a frame of reference more appropriate to our needs.  Suppose we want
to draw a filled circle defined by \(x^2 + y^2 < r^2\).  This circle is
centered on (0,0) in user space, but the center of device space is
(width/2, height/2).  Let's pick \(r < 1\), then the whole circle will
be inside a box \([-1,1]\times[-1,1]\).  We transform the device
coordinates (i,j) to user coordinates (x,y), to make drawing the circle
easier.

We need fractional numbers now, for which C has the ``float'' type.
Floating point arithmetic stores a number of significant digits, and
a scaling exponent (like scientific notation).  We use floating point
constants ``2.0'', otherwise ``/'' performs integer division, and we
also exploit the fact that C ``promotes'' the other operands to floating
point as necessary.

With this coordinate transformation, our circle looks like this:

    $ ./mandelbrot
    ------------------------------------------------------------------------
    ------------------------------------------------------------------------
    -----------------------xxxxxxxxxxxxxxxxxxxxxxxxxxx----------------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    --------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-------------
    -----------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------
    ---------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx--------
    -------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx------
    ------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----
    -----xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----
    ----xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx---
    ----xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx---
    ----xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx---
    -----xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----
    ------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-----
    -------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx------
    ---------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx--------
    -----------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------
    --------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx-------------
    -----------------xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx----------------
    -----------------------xxxxxxxxxxxxxxxxxxxxxxxxxxx----------------------
    ------------------------------------------------------------------------
    $

3 years agoDevice coordinates.
Claude Heiland-Allen [Fri, 20 Dec 2013 10:56:29 +0000 (10:56 +0000)]
Device coordinates.

Western text is written left to right in rows from top
to bottom.  For text-based images, this gives an origin
(0,0) at the top left, with the next character to the
right at (1,0) and the first character in the line below
at (0,1).  We can generate a table of coordinates with
two nested for loops: the outer one loops over rows, and
the inner one loops over columns within each row.

We would print the coordinates in a table, but so far we
only know how to print strings.  The C standard library
includes the function ``printf()'' for formatted output.
Its first parameter is a format string, which contains
conversion specifiers prefixed by ``\%''.  You can read
the full list in its reference manual by running
``man 3 printf''.

We want to print our int coordinates in decimal, so we put
``\%d'' in the format string, once for each coordinate.
printf() is variadic: the number of parameters and the type
of each parameter depends on the format string.  We want to
have each row on a separate line, so we need to emit a newline
character `` '\char`\\n' '' at the end of each row.  The
``putchar()'' function prints a single character.

Our program now outputs a table of device coordinates:

    $ ./mandelbrot
    (0,0)(1,0)(2,0)(3,0)(4,0)(5,0)(6,0)(7,0)(8,0)(9,0)
    (0,1)(1,1)(2,1)(3,1)(4,1)(5,1)(6,1)(7,1)(8,1)(9,1)
    (0,2)(1,2)(2,2)(3,2)(4,2)(5,2)(6,2)(7,2)(8,2)(9,2)
    (0,3)(1,3)(2,3)(3,3)(4,3)(5,3)(6,3)(7,3)(8,3)(9,3)
    (0,4)(1,4)(2,4)(3,4)(4,4)(5,4)(6,4)(7,4)(8,4)(9,4)
    (0,5)(1,5)(2,5)(3,5)(4,5)(5,5)(6,5)(7,5)(8,5)(9,5)
    (0,6)(1,6)(2,6)(3,6)(4,6)(5,6)(6,6)(7,6)(8,6)(9,6)
    (0,7)(1,7)(2,7)(3,7)(4,7)(5,7)(6,7)(7,7)(8,7)(9,7)
    (0,8)(1,8)(2,8)(3,8)(4,8)(5,8)(6,8)(7,8)(8,8)(9,8)
    (0,9)(1,9)(2,9)(3,9)(4,9)(5,9)(6,9)(7,9)(8,9)(9,9)
    $

3 years agoPrinting all arguments: a ``for'' loop.
Claude Heiland-Allen [Thu, 19 Dec 2013 19:48:36 +0000 (19:48 +0000)]
Printing all arguments: a ``for'' loop.

So far the program just prints its name, and ignores
any other arguments.  Suppose we want to print all of
them.  The idiomatic way to do this is to write a loop.

The ``for'' loop syntax has three clauses.  In the
first we declare and initialize a counter variable.
This happens once, before any looping happens.

In the second clause we check that it's within the
valid range of argument indices.  If the condition
is false, then the control flow skips beyond the
end of the loop body without executing it.

If the condition clause is true, then the loop body is
executed.  In our program, we print the i'th command
line argument.  When the control flow reaches the end
of the loop body, the third clause is executed before
going back to the start of the loop: in our program we
increment the counter.  This process repeats so long as
the condition clause remains true.

An example invocation of our program might look like

    $ ./mandelbrot one two three
    ./mandelbrot
    one
    two
    three
    $

where ``\$'' is the command prompt.

3 years agoInput and output.
Claude Heiland-Allen [Thu, 19 Dec 2013 19:27:09 +0000 (19:27 +0000)]
Input and output.

A program that returns a status code isn't
very interesting.  To interact with the
outside world needs input from the world and
output to the world.  The C standard library
supports this within the ``stdio.h'' header
file.  One function declared in the header is
``puts()'', which outputs a string (a sequence
of characters in memory, terminated by the null
character `` '\char`\\0' ''.

A program can also take input from its
command line arguments, these are passed as
parameters to ``main'': the first parameter
is the integer number of command line arguments,
and the second contains the arguments themselves.
The collection of arguments and the count
includes the program name as the first argument.

The count has the type ``int'', a small integer:
the mathematical integers extend to infinity, but
computers are based on small fixed-width numbers
and C keeps at a fairly low level of abstraction
above the hardware.

In C, ``*'' denotes a pointer type.  Pointers are
references to areas of memory containing data.
The arguments are made up of characters stored
in consecutive memory locations, and they are
passed to the program as a pointer to consecutive
memory locations, each containing another pointer
to the first character in each argument.  This
gives the type ``char **'', the names ``argc''
and ``argv'' are just a common convention..

Getting the data from the memory location referred
to by a pointer is called ``dereferencing'', and in
C you can use ``*'' (as a unary operator this time)
to get the location pointed to, or ``[n]'' to get
the n'th item, with each item stored consecutively
in memory.  In C, counting starts from 0 for the
first item.

The program now prints the name used to invoke it
to its standard output stream.

3 years agoIgnoring generated files.
Claude Heiland-Allen [Thu, 19 Dec 2013 19:03:16 +0000 (19:03 +0000)]
Ignoring generated files.

A version control system allows you to keep track
of the change history of a collection of files
over time.  Git is one of the more popular tools.
The ``git status'' command gives a brief summary
of what has changed since the last ``git commit'',
for example, listing files that haven't been added
to the repository.  Some of these files might be
temporary or system-specific -- for example, compiled
programs.  Adding them to the ``.gitignore'' file
makes the status output less noisy, allowing us to
focus on important changes.

3 years agoBuilding a C program.
Claude Heiland-Allen [Thu, 19 Dec 2013 18:52:34 +0000 (18:52 +0000)]
Building a C program.

A computer can't run a C program without translating
it to machine code for the CPU, a task performed by a
compiler.  The GNU Compiler Collection includes the
gcc compiler for C.  We'll be using the C99 standard
dialect of C.  We'll enable every warning, so that we
get hints about possible bugs.  We'll set optimization
to a high level, so that we don't have to wait so long
for the program to run (it might take longer to compile).

A Makefile contains rules for building targets from
prerequisites.  Running ``make'' will check if the
default target is up to date compared to its
prerequisites, and if not it will run the commands listed
in the rules.  This happens recursively, so that if you
change one file in a large project, only the affected
path through the tree of dependencies will be rebuilt.

3 years agoA minimal C program.
Claude Heiland-Allen [Thu, 19 Dec 2013 18:46:43 +0000 (18:46 +0000)]
A minimal C program.

Every program needs an entry point -- in the programming
language C this is a function called ``main''.  When the
program has finished running, it exits, returning a status
code to the operating system to signal success (with 0) or
failure (non-zero values).