Happy π day!

March 14, 2011 | Math, My Projects | By: Mark VandeWettering

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?

Comments

Pingback from Playing with Pi | Learn Something
Time 3/22/2011 at 12:03 am

[…] I had a pretty good idea how he was going to do it, as he’d recently posted something really cool about a way to use the Mandelbrot set to find pi.  (Sure enough that’s what he did.) […]

Write a comment