Double precision.
authorClaude Heiland-Allen <claude@mathr.co.uk>
Fri, 20 Dec 2013 23:34:25 +0000 (23:34 +0000)
committerClaude Heiland-Allen <claude@mathr.co.uk>
Fri, 20 Dec 2013 23:34:25 +0000 (23:34 +0000)
commitecb18f012205e586c683a15c5829d71d06dcdb32
tree9ace1dc2a4403fd777af9a896412ab57053263d8
parent7b71abc5e0df24258b1bfe47369018995061a09f
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.
mandelbrot.c