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?

Share Button
Be Sociable, Share!

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