mandelbrot-set-fourier.sage

## Compute 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 ctab[m] 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]
## The value 2^(2*m+1)*b_m is an integer: its first few values are:
##   [-1, 1, -8, 15, 0, -94, -512, 987, 0, -7346]

## Two variants are given below: one using exact (integer) computations, and
## one using floating-point computations (which seems to achieve satisfactory
## precision).

## Rational or fp computation: 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 Sage program).

## Integer computation: b_m = B(0,m+1)/2^(2*m+1) where B(0,0)=1/2 (only non-int),
## and B(n,m) = 2^(2*m+3-2^(n+2))*beta(n,m), so:
## B(n-1,m) = 2^(2^(n+1)-1)*B(n,m) - 2^(2^(n+1)-4)*sum(B(n-1,j)*B(n-1,m-j), j=2^n-1..m-2^n+1) - 2*B(0,m-2^n+1) if m>=2^n-1, 0 otherwise
## (B(n,m) is written bigbtab[m][n] in this Sage program).

## EXACT COMPUTATION ##

target = 4097
sizen = ceil(log(target+4)/log(2))

bigbtab = []
for m in range(target+1):
    if m%100 == 0:
        print m
    if m==0:
        bigbtab.append([1/2 if n==0 else None for n in range(sizen)])
    else:
        subtab = [None for n in range(sizen)]
        subtab[sizen-1] = 0
        for n in range(sizen-1,0,-1):
            if m>=2^n-1:
                subtab[n-1] = Integer(2^(2^(n+1)-1)*subtab[n] - 2^(2^(n+1)-4)*sum([bigbtab[k][n-1]*bigbtab[m-k][n-1] for k in range(2^n-1,m-2^n+2)]) - 2*bigbtab[m-2^n+1][0])
            else:
                subtab[n-1] = 0
        bigbtab.append(subtab)

ctab = [bigbtab[m+1][0]/2^(2*m+1) for m in range(target)]

parametric_plot((cos(x)+sum([N(ctab[i])*cos(i*x) for i in range(len(ctab))]), sin(x)-sum([N(ctab[i])*sin(i*x) for i in range(len(ctab))])), (x,0,2*pi), aspect_ratio=1, plot_points=1024).show(xmin=-2.05, xmax=0.55, ymin=-1.3, ymax=1.3)

# sum([parametric_plot((cos(x)+sum([N(ctab[i])*cos(i*x) for i in range(target>>j)]), sin(x)-sum([N(ctab[i])*sin(i*x) for i in range(target>>j)])), (x,0,2*pi), aspect_ratio=1, plot_points=1024, color=Color(j/floor(log(target)/log(2)),0,1 - j/floor(log(target)/log(2)))) for j in range(floor(log(target)/log(2)),0,-1)])

## FLOATING-POINT COMPUTATION ##

target = 4097
sizen = ceil(log(target+4)/log(2))

betatab = []
for m in range(target+1):
    if m%100 == 0:
        print m
    if m==0:
        betatab.append([N(1) if n==0 else None for n in range(sizen)])
    else:
        subtab = [None for n in range(sizen)]
        subtab[sizen-1] = 0
        for n in range(sizen-1,0,-1):
            if m>=2^n-1:
                subtab[n-1] = (subtab[n] - sum([betatab[k][n-1]*betatab[m-k][n-1] for k in range(2^n-1,m-2^n+2)]) - betatab[m-2^n+1][0])/2
            else:
                subtab[n-1] = 0
        betatab.append(subtab)

ctab = [betatab[m+1][0] for m in range(target)]

parametric_plot((cos(x)+sum([ctab[i]*cos(i*x) for i in range(len(ctab))]), sin(x)-sum([ctab[i]*sin(i*x) for i in range(len(ctab))])), (x,0,2*pi), aspect_ratio=1, plot_points=1024).show(xmin=-2.05, xmax=0.55, ymin=-1.3, ymax=1.3)

# sum([parametric_plot((cos(x)+sum([ctab[i]*cos(i*x) for i in range(target>>j)]), sin(x)-sum([ctab[i]*sin(i*x) for i in range(target>>j)])), (x,0,2*pi), aspect_ratio=1, plot_points=1024, color=Color(j/floor(log(target)/log(2)),0,1 - j/floor(log(target)/log(2)))) for j in range(floor(log(target)/log(2)),0,-1)])

######################################################################

## Draw the parametrization itself:
twopiN = N(2*pi)
plot(lambda x: cos(twopiN*x)+sum([ctab2[i]*cos(twopiN*i*x) for i in range(len(ctab))]), (x,0,1)) + plot(lambda x: sin(twopiN*x)-sum([ctab2[i]*sin(twopiN*i*x) for i in range(len(ctab))]), (x,0,1), color="red")

## Draw a conformal grid:
angles = [N(j*pi/60) for j in range(120)] ; arguts = [N(j*pi/60) for j in range(11)] ; sum([parametric_plot((cos(angle)*exp(x)+sum([ctab2[i]*cos(i*angle)*exp(-i*x) for i in range(len(ctab))]), sin(angle)*exp(x)-sum([ctab2[i]*sin(i*angle)*exp(-i*x) for i in range(len(ctab))])), (x,0.,N(pi/6)), color="red") for angle in angles]) + sum([parametric_plot((cos(x)*exp(argut)+sum([ctab2[i]*cos(i*x)*exp(-i*argut) for i in range(len(ctab))]), sin(x)*exp(argut)-sum([ctab2[i]*sin(i*x)*exp(-i*argut) for i in range(len(ctab))])), (x,0.,N(2*pi))) for argut in arguts])

https://commons.wikimedia.org/wiki/File:Jungreis.svg