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.
#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) ; }
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?