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?
Pingback: Playing with Pi | Learn Something