mdbtxt1
mdbtxt2
Proceed to Safety

# Laurent Series

Robert P. Munafo, 2000 Apr 19.

In general, a Laurent series is an infinite series, like a Taylor series, whose terms can be used to approximate some function. One important difference from Taylor series is that the Laurent series method can be used for functions that are less well-behaved, including non-differentiable functions.

In Mu-Ency "Laurent series" refers to a specific Laurent Series described by Ewing and Schober. The function it models is the conformal map psi mapping the unit disk onto the Mandelbrot Set.

Ewing and Schober showed that the area of the Mandelbrot Set could be computed according to the formula:

AM = pi . ( 1 - sigma(n=0 to infinity)( n . bn2 ) )

Where bn are the coefficients of the Laurent series. The calculation of the coefficients is discussed below.

Because the bn values are squared, the sum always increases, and therefore the value of the whole expression always decreases, as successive finite sum approximations are computed. Thus, for a finite n, the formula gives an upper bound to the area.

Attempts to estimate the area using the pixel counting method yield estimates between 1.50 and 1.55. When Ewing and Schober evaluated this series to 240,000 terms, it looked as though the series was going to converge somewhere between 1.60 and 1.65. The series converges so slowly that Ewing and Schober did not try to determine the value of the limit any more precisely than this.

The following table shows the convergence behavior of the series. Scott Huddleston compiled most of the values in this table; the last (for 240,000 terms) is from the Ewing and Schober paper. Each entry in the table represents a doubling of the number of terms used to generate the estimate (except the last two, which are approximately double):

 number of terms area upper bound delta 1.50659 + e / ln(n) (see below) 128 2.02781 2.06718 o 0.05029 256 1.97752 1.99710 o 0.05001 512 1.92751 1.94260 o 0.03217 1024 1.89534 1.89900 o 0.04073 2048 1.85461 1.86333 o 0.02009 4096 1.83452 1.83360 u 0.02836 8192 1.80616 1.80844 o 0.01980 16384 1.78636 1.78688 o 0.01953 32768 1.76683 1.76820 o 0.01346 65536 1.75337 1.75184 u 0.01490 115232 1.73847 1.73997 o 0.01097 240000 1.72750 1.72615 u

Notice how the first differences go up and down: the sequence converges erratically. Also notice how, even though we keep doubling the number of terms, each doubling brings the area estimate down almost as much as the previous doubling did. That means if we want to use this method to get even two good digits of the Mandelbrot set's area, we'll have to double many more times.

When this was discussed in 1992, many believed that the limit of this series is greater than the limit of the pixel counting method. Ewing and Schober suggested either that the formula based on the Laurent series might be wrong (due most probably to poor understanding of the boundary of the Mandelbrot set and how the conformal map approaches it), or that perhaps the boundary has positive area. (It is known that the boundary's Hausdorff dimension is 2.0, but that doesn't guarantee that it has positive area.)

However, it seems most likely given the first-difference data, that the series does converge to the same value as given by pixel-counting, but that the rate of convergence is extremely slow. A good formula for the difference between a Ewing/Schober estimate and the pixel-counting area is:

difference = e / ln(number-of-terms)

The value of 1.50659 + e/ln(n) is shown in the 4th column of the table above. Each value is marked with an o if it is an overestimate of column 2, (that is, higher than the actual Laurent series sum), and a u if it is an underestimate. As you can see, it starts out consistently higher, but after the 5th row of the table it continually goes back and forth between being an overestimate and being an underestimate. Thus, it appears to be a good statistical model of the behavior of the Laurent series sum.

Based on that, the number of terms needed to get a higher bound which agrees with the actual Mandelbrot set area to within N digits past the decimal point is

number-of-terms = ee 10N

Which gives 6.39 × 1011 terms to get the first digit, and an inconceivable 1.13 × 10118 terms to get the second!

The Coefficients and Details of Computation

This is the definition of the conformal mapping psi:

psi(w) = w + sigma(n = 0 to infinity)(bn / wn)

This expands to:

psi(w) = w - 1 / 2 + 1 / (8 w) - 1 / (4 w2) + 15 / (128 w3) + 0 / w4 - 47 / (1024 w5) - 1 / (16 w6) + 987 / (32768 w7) + 0 / w8 - 3673 / (262144 w9) + 1 / (32 w10) - 61029 / (4194304 w11) . . .

The coefficients of the series are b0 = -1/2, b1 = 1/8, b2 = -1/4, b3 = 15/128, b4 = 0, b5 = -47/1024, etc. The terms bn are the same in the formula for the Mandelbrot set area which was given above:

AM = pi . ( 1 - sigma(n=0 to infinity)( n . bn2 ) )

The terms bn are defined recursively in terms of three other sets of terms called w(n,m), aj, and u(n,k). The sets w and u are infinite 2-dimensional arrays (kind of like Pascal's triangle) and the other two (b and a) are ordinary infinite series. All of their values are rational numbers, with powers of 2 in the denominator. Here are the definitions Ewing and Schober gave for the four sets of terms:

bn = -1/2, if n = 0
otherwise bn = ( - w(n,n+1) ) / n

w(n,m) = 0, if n = 0
otherwise w(n,m) = am-1 + w(n-1,m) + sigma(j = 0 ... m-2)(aj . w(n-1,m-j-1))

aj = u(0,j+1)

u(n,k) = 1, if 2n-1 = k
otherwise u(n,k) = sigma(j = 0 .. k)(u(n-1,j) . u(n-1,k-j)), if 2n-1 > k
otherwise u(n,k) = 0, if 2n+1-1 > k
otherwise u(n,k) = (u(n+1,k) - sigma(j = 1 ... k-1)(u(n,j) . u(n,k-j))) / 2

(These definitions were adapted from the Maple code given by Gerald Edgar and listed below)

Recursion and the Importance of Storing Values

As you can see, the values of w(n,m) and u(n,k) are usually computed in terms of several other w(n,m) and/or u(n,k) values. When this happens, smaller n, m, n and k are used and the recursive process repeats until they all manage to get down to one of the trivial cases (like w(n,m) when n=0).

This means that if you try to write a computer program that computes bn terms, and defines these four sets of terms as functions that call themselves and each other, the program will take exponentially longer as you go to higher values of n. Here is a table showing how many times one of the four functions will be invoked using this "simple recursive" approach:

 to compute bn for n from 0 to: function calls ratio series sum 0 1 3.1415 15.00 1 15 3.0925 7.13 2 107 2.6998 6.06 3 648 2.5703 5.57 4 3608 2.5703 5.28 5 19063 2.5372 5.12 6 97570 2.4636 5.01 7 488768 2.4437 4.98 8 2411908 2.4437

Clearly that won't work — each time you add a term it's taking about 4 times as long as all the previous terms put together. The exponential rate is decreasing, but not very quickly. The first entry in the table above is for 128 terms. Even if the average exponential rate for the first 128 terms turned out to be as low as 3.00, the calculation would require something like 1061 steps!

This is a common problem in numerical mathematics, and the solution is to tell the computer to remember each value of w(n,m), aj, and u(n,k) the first time it is computed. The recursive function will not have to be evaluated the second and subsequent times, saving many steps. Here is the same table for the first 8 terms, using this improvement:

 to compute bn for n from 0 to: function calls ratio series sum 0 1 3.1415 11.00 1 11 3.0925 1.82 2 20 2.6998 1.55 3 31 2.5703 1.42 4 44 2.5703 1.34 5 59 2.5372 1.37 6 81 2.4636 1.25 7 101 2.4437 1.22 8 123 2.4437 16 381 2.2730 32 1295 2.1830 72 5857 2.0928

Now, to get n terms, you only have to perform a little more than n2 steps. That is certainly a significant improvement over the exponential pattern. However, there is a substantial drawback in the fact that you have to store the terms. The w(n,m) and u(n,k) terms create the greatest impact, because they are square arrays: to compute n terms, you end up computing n2 values of w(n,m) and u(n,k). That means that when they computed 240,000 terms of the series, Ewing and Schober's computer had to store about 57,600,000,000 terms of w(n,m) and an equal number of u(n,k).

Given this, and the fact that the series converges so slowly, it is clear that the Ewing-Schober sequence has little practical value as a numerical algorithm for computing the area of the Mandelbrot set. In the example above where we considered computing 2 digits past the decimal point, it took 10118 terms, and that would mean storing 10236 values in memory, give or take a few.

Gerald Edgar has provided Maple code to compute the coefficients, using two different methods. The first method is described as being equivalent to the method of Ewing and Schober:

u := proc(n,k) local j; option remember; if 2^n-1=k then 1; elif 2^n-1>k then sum('u(n-1,j)*u(n-1,k-j)',j=0..k); elif 2^(n+1)-1>k then 0; else (u(n+1,k)-sum('u(n,j)*u(n,k-j)',j=1..k-1))/2; fi; end; a := proc(j) option remember; u(0,j+1) end; w := proc(n,m) local j; option remember; if n=0 then 0; else a(m-1)+w(n-1,m)+sum('a(j)*w(n-1,m-j-1)',j=0..m-2); fi; end; b := proc(m) option remember; if m=0 then -1/2; else -w(m,m+1)/m; fi; end;

Adam Majewski provides this translation into maxima:

u[n,k] := block( [j], if 2^n-1=k then 1 else if 2^n-1>k then sum(u[n-1,j]*u[n-1,k-j],j,0,k) else if 2^(n+1)-1>k then 0 else (u[n+1,k]-sum(u[n,j]*u[n,k-j],j,1,k-1))/2 );   a[j] := u[0,j+1];   w[n,m] := block( [j], if n=0 then 0 else a[m-1]+w[n-1,m]+sum(a[j]*w[n-1,m-j-1],j,0,m-2) );   b(m) := if m=0 then -1/2 else -w[m,m+1]/m ;

Gerald Edgar's second method (also in Maple) is described as being more efficient:

u := proc(n,k) local j; option remember; if 2^n-1=k then 1; elif 2^n-1>k then sum('u(n-1,j)*u(n-1,k-j)',j=0..k); elif 2^(n+1)-1>k then 0; else (u(n+1,k)-sum('u(n,j)*u(n,k-j)',j=1..k-1))/2; fi; end; a := proc(j) option remember; u(0,j+1) end; v := proc(k,j) local l; option remember; if k=1 then -a(j-2)-sum('v(1,l)*a(j-l-1)',l=2..j-1); else v(1,j-k+1) + v(k-1,j-1)+sum('v(1,l)*v(k-1,j-l)',l=2..j-k); fi; end; b := proc(j) local k; option remember; if j=0 then -a(0); else -a(j)-sum('v(k,j)*b(k)',k=1..j-1); fi; end;

Adam Majewski also provides this numerical version in maxima:

betaF[n,m]:=block ( [nnn:2^(n+1)-1], if m=0 then 1.0 else if ((n>0) and (m < nnn)) then 0.0 else (betaF[n+1,m]- sum(betaF[n,k]*betaF[n,m-k], k, nnn, m-nnn) -betaF[0,m-nnn])/2.0 );   b(m):=betaF[0,m+1];

References

Acknowledgments

Area estimates table (Laurent Series method):

The 240000 estimate is from the Ewing & Schober article.

The rest are from: Scott Huddleston (scott at crl labs tek com)

Opinions on the discrepancy: Stan Isaacs (isaacs at hpcc01 hp com)

Ewing & Schober reference: Jeff Shallit (shallit at graceland waterloo edu)

Laurent series formula and Maple code: Gerald Edgar (edgar at math ohio-state edu)

Here is Scott Huddleston's complete table of area estimates: # terms: area estimate 72: area < 2.09288 128: area < 2.02781 180: area < 2.01237 256: area < 1.97752 360: area < 1.94961 512: area < 1.92751 720: area < 1.91255 1024: area < 1.89534 1440: area < 1.87172 2048: area < 1.85461 2880: area < 1.84576 4096: area < 1.83452 5760: area < 1.81649 8192: area < 1.80616 11520: area < 1.79642 16384: area < 1.78636 23040: area < 1.7747 32768: area < 1.76683 46080: area < 1.75936 65536: area < 1.75337 92160: area < 1.74321 115232: area < 1.73847

Here is the beginning of the Laurent series, in both reduced and unreduced forms:

 n reduced unreduced 0 -1 / 2 -1 / 2 1 1 / 8 1 / 8 2 -1 / 4 -8 / 32 3 15 / 128 15 / 128 4 0 / 1 0 / 512 5 -47 / 1024 -94 / 2048 6 -1 / 16 -512 / 8192 7 987 / 32768 987 / 32768 8 0 / 1 0 / 131072 9 -3673 / 262144 -7346 / 524288 10 1 / 32 65536 / 2097152 11 -61029 / 4194304 -122058 / 8388608 12 0 / 1 0 / 33554432 13 -689455 / 33554432 -2757820 / 134217728 14 -21 / 512 -22020096 / 536870912 15 59250963 / 2147483648 16 0 17 -164712949 / 17179869184 18 39 / 2048 19 -2402805839 / 274877906944 20 -1 / 64 21 -4850812329 / 2199023255552 22 29 / 2048 23 -18151141041 / 70368744177664 24 0 25 3534139462275 / 562949953421312 26 -1039 / 131072 27 -22045971176589 / 9007199254740992 28 -1 / 256 29 -750527255965871 / 72057594037927936 30 -4579 / 524288 31 54146872254247683 / 9223372036854775808 32 0 33 -155379776183158669 / 73786976294838206464 34 2851 / 1048576 35 -6051993294029466699 / 1180591620717411303424 36 -1 / 1024 37 7704579806709870203 / 9444732965739290427392 38 92051 / 16777216 39 -403307733528668035403 / 302231454903657293676544 40 0 41 1650116480759617184697 / 2417851639229258349412352 42 -229813 / 67108864 43 36124726440442241978477 / 38685626227668133590597632 44 -41 / 4096 45 -225851495844149964787753 / 309485009821345068724781056 46 564373 / 67108864 47 -35761228458796476847725737 / 19807040628566084398385987584 48 0 49 362376876750551361794705823 / 158456325028528675187087900672 50 -29407003 / 8589934592 51 -6510398483578238274151194427 / 2535301200456458802993406410752 52 33 / 8192 53 74815618913797220433481657203 / 20282409603651670423947251286016 54 -30057875 / 34359738368 55 -698617278028915809388280344009 / 649037107316853453566312041152512 56 0 57 -8675905413734991085610532783493 / 5192296858534827628530496329220096 58 -27868893 / 68719476736 59 -375687870961637050293461860951517 / 83076749736557242056487941267521536 60 1 / 4096 61 -1418434432207399687114226995905967 / 664613997892457936451903530140172288 62 -11847286243 / 1099511627776 63 1084116104452462070609082665064238307 / 170141183460469231731687303715884105728 64 0

From the Mandelbrot Set Glossary and Encyclopedia, by Robert Munafo, (c) 1987-2022.     Mu-ency index

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2019 Feb 01. s.27