mandelbrot-fourier.c

// Compute the coefficients of the Jungreis function, i.e., the
// Fourier coefficients of the harmonic parametrization of the
// boundary of the Mandelbrot set, using the formulae given in
// following paper: John H. Ewing & Glenn Schober, "The area of the
// Mandelbrot set", Numer. Math. 61 (1992) 59-72 (esp. formulae (7)
// and (9)).  (Note that their numerical values in table 1 give the
// coefficients of the inverse series.)

// The coefficients betatab[m+1][0] are the b_m such that
// z + sum(b_m*z^-m) defines a biholomorphic bijection from the
// complement of the unit disk to the complement of the Mandelbrot
// set.  The first few values are:
//   {-1/2, 1/8, -1/4, 15/128, 0, -47/1024, -1/16, 987/32768, 0, -3673/262144}

// Formula: b_m = beta(0,m+1) where beta(0,0)=1,
// and by downwards induction on n we have:
// beta(n-1,m) = (beta(n,m) - sum(beta(n-1,j)*beta(n-1,m-j), j=2^n-1..m-2^n+1) - beta(0,m-2^n+1))/2 if m>=2^n-1, 0 otherwise
// (beta(n,m) is written betatab[m][n] in this C program).

#include <stdio.h>
#include <math.h>
#include <float.h>
#include <assert.h>

#ifndef NBCOEFS
#define NBCOEFS 8193
#endif

#ifndef SIZEN
#define SIZEN 14
#endif

#if ((NBCOEFS+4)>>SIZEN) != 0
#error Please increase SIZEN
#endif

double betatab[NBCOEFS+1][SIZEN];

int
main (void)
{
  for ( int n=0 ; n<SIZEN ; n++ )
    betatab[0][n] = 1;
  for ( int m=1 ; m<=NBCOEFS ; m++ )
    {
      for ( int n=SIZEN ; n>=1 ; n-- )
	{
	  if ( m < (1<<n)-1 )
	    betatab[m][n-1] = 0;
	  else
	    {
	      double v;
	      assert (n < SIZEN);
	      v = betatab[m][n];
	      for ( int k=((1<<n)-1) ; k<=m-(1<<n)+1 ; k++ )
		v -= betatab[k][n-1] * betatab[m-k][n-1];
	      v -= betatab[m-(1<<n)+1][0];
	      betatab[m][n-1] = v/2;
	    }
	}
      printf ("%d\t%.17e\t%a\n", m-1, betatab[m][0], betatab[m][0]);
    }
}