brainwagon

It’s 3/14 again, and that means that it’s π day! Huzzah. This year, I thought I’d try implementing a way of computing π which was entirely new to me: finding π hiding inside the Mandelbrot set.

David Boll made a posting back in 1991 to sci.math:

I posted this to alt.fractals about a month ago, and it occurred to me 
  that readers of this newsgroup would probably be interested also. 
  I was trying to show that the 'neck' of the mandelbrot set at (-.75,0) is 
  actually of zero thickness. Accordingly, I wrote a quickie program and 
  started checking the number of iterations that points of the form 
  (-.75,dy) went thru before blowing up (with dy being a small number). 
  Here's a table of results: 
            dy    # of iterations 
            .1             33 
            .01           315 
            .001         3143 
            .0001       31417 
            .00001     314160 
  Notice anything peculiar about the # of iterations? How about the product 
  of the # of iterations with dy? 

Pretty neat. Surf over to his page on the subject for more information and for Gerald Edgar’s proof. I wrote my own implementation which you can find below.
[sourcecode lang=”C”]
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <inttypes.h>
#include <complex.h>

/*
* _______
* ( _ ) _
* | | | | _ __ ___ __ _ _ __ __| |
* | | | | | ‘_ ` _ \ / _` | ‘_ \ / _` |
* | | | | | | | | | | (_| | | | | (_| |
* |_| |_| |_| |_| |_|\__,_|_| |_|\__,_|
*
* A program that computes approximations to pi by iterating over
* the Mandelbrot set. Not terrifically useful really, there are
* much faster ways, but this was a way which I had not seen
* before, so I thought it was worth archiving in programmatic form.
*
* Written for pi-day, 3/14/2011 by Mark VandeWettering.
*
* References:
* https://home.comcast.net/~davejanelle/mandel.html
*/

uint64_t
mand(double complex c, uint64_t limit)
{
uint64_t i ;
double complex z ;

z = 0 ;
for (i=1; i<limit; i++) {
z = z*z + c ;
if (cabs(z) > 2.0) break ;
}
return i ;
}

main(int argc, char *argv[])
{
int i ;
uint64_t limit = 10LL ;
uint64_t iter ;
double complex c ;
double y = 1. ;

if (argc == 1) {
printf("computing some approximations to pi using the mandelbrot set.\n") ;
printf(" (close approximations might take a few seconds, warning…)\n") ;
for (i=0; i<9; i++) {
c = -0.75 + y * I ;
iter = mand(c, limit) ;
printf("%9" PRIu64 " iterations, pi is %.8lf, +/- %.8lf\n", iter, 10.*iter/limit, y) ;
y /= 10 ;
limit *= 10LL ;
}
}
exit(0) ;
}
[/sourcecode]

Bonus points: This program when run produces a result which isn’t correct during the last iteration (the Mandelbrot iteration runs to limit, instead of stopping where it should). I suspect that it’s got a bug in it somewhere, and isn’t using enough precision. I’ve just got back from a 24 hour trip to London, a visit from a locksmith (don’t ask) and now 12 hours of sleep, but my brain isn’t functioning right now. Can anyone else spot the issue?