# 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:

A_{M} = pi ^{.} ( 1 - sigma_{(n=0 to infinity)}( n ^{.} b_{n}^{2} ) )

Where b_{n} are the coefficients of the Laurent
series. The calculation of the coefficients is discussed below.

Because the b_{n} 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 = e^{e 10N}

Which gives 6.39 × 10^{11} terms to get the first digit, and an
inconceivable 1.13 × 10^{118} 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)}(b_{n} / w^{n})

This expands to:

psi(w) = w - 1 / 2 + 1 / (8 w) - 1 / (4 w^{2}) + 15 / (128 w^{3})
+ 0 / w^{4} - 47 / (1024 w^{5}) - 1 / (16 w^{6}) + 987 / (32768 w^{7})
+ 0 / w^{8} - 3673 / (262144 w^{9}) + 1 / (32 w^{10})
- 61029 / (4194304 w^{11}) . . .

The coefficients of the series are
b_{0} = -1/2,
b_{1} = 1/8,
b_{2} = -1/4,
b_{3} = 15/128,
b_{4} = 0,
b_{5} = -47/1024, etc. The terms b_{n} are the same in the
formula for the Mandelbrot set area which was given above:

A_{M} = pi ^{.} ( 1 - sigma_{(n=0 to infinity)}( n ^{.} b_{n}^{2} ) )

The terms b_{n} are defined recursively in terms of three other sets
of terms called w_{(n,m)}, a_{j}, and u_{(n,k)}. The
sets w and u are infinite 2-dimensional arrays (kind of like
Pascal's triangle) and the other two 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:

b_{n} = -1/2, if n = 0

otherwise b_{n} = ( - w_{(n,n+1)} ) / n

w_{(n,m)} = 0, if n = 0

otherwise w_{(n,m)} = a_{m-1} + w_{(n-1,m)} + sigma_{(j = 0 ... m-2)}(a_{j} ^{.} w_{(n-1,m-j-1)})

a_{j} = u_{(0,j+1)}

u_{(n,k)} = 1, if 2^{n}-1 = k

otherwise u_{(n,k)} = sigma_{(j = 0 .. k)}(u_{(n-1,j)} ^{.} u_{(n-1,k-j)}), if 2^{n}-1 > k

otherwise u_{(n,k)} = 0, if 2^{n+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
b_{n} 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 b _{n} forn 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 10^{61} 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)}, a_{j}, 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 b _{n} forn 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
n^{2} 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 n^{2} 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 10^{118} terms, and that would mean
storing 10^{236} 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 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-2016. Mu-ency index

This page was written in the "embarrassingly readable" markup language RHTF, and was last updated on 2016 Jan 02. s.11