{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Mandelbrot Polynomials and Matrices\n",
"\n",
"```{image} ../Figures/Mandelbrot/explosion32comp.gif\n",
":height: 400px\n",
":alt: Exploding Mandelbrot 32\n",
":align: center\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A note to the student/reader\n",
"\n",
"This unit uses some ideas of the late [Benoit Mandelbrot](https://en.wikipedia.org/wiki/Benoit_Mandelbrot) and some known facts about the so-called _Mandelbrot Set_, together with some of our own ideas, to help you to learn the following:\n",
"\n",
"1. A bit more Python programming, including more practice with the use of _numerical libraries_ from NumPy, _iteration_ and _recursion_; some new things, such as _loop invariants_, _automatic differentiation_, _symbolic computation_, and a bit more about graphics.\n",
"2. A bit more about polynomials: the _cubic formula_ (we won't need the quartic formula, and even our use of the cubic formula is a bit contrived; but we think it's fun), the _cost of numerical solution of polynomials_ (we will point you to the current numerical champion polynomial solver, [MPSolve](https://en.wikipedia.org/wiki/MPSolve)), and the rather necessary-to-know notions of _numerical stability_ and _conditioning_. Most of these topics are not taught as thoroughly as they should be, in one's first numerical analysis course; so this material here should strengthen the results when you do encounter them.\n",
"3. The surprising use of _eigenvalues_ to find roots of polynomials. The first treatment of eigenvalues (typically in a Linear Algebra course, even though the problem is not, strictly speaking, linear) usually runs the other way around, and defines the eigenvalues of a matrix in terms of the _characteristic polynomial_ of the matrix. Indeed, we will teach you a concept not in the textbooks: namely, the concept of a _minimal height companion matrix_ and we will show you one such for the Mandelbrot polynomials.\n",
"4. A bit more about _dynamical systems_, especially about _iteration_ and _composition_.\n",
"5. An excellent approximate formula for the _largest magnitude root_ of the Mandelbrot polynomial.\n",
"6. An honest-to-goodness analytic solution to the Mandelbrot iteration (this is a very new result, published only in 2021), valid for all $c$ _outside_ the Mandelbrot set.\n",
"7. That we don't know everything about Mandelbrot polynomials and matrices, and that you might be able to answer some open questions."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A note to the Instructor\n",
"\n",
"The students tend to really like this unit. Mandelbrot wrote so well that it is quite possible to ask students to read his papers (especially his math education papers). Read them yourself, if you haven't already. \n",
"\n",
"We have put some of our own research into this unit. We believe that the results we have included are accessible (and of course we believe that they are interesting). There are some open problems here that we believe some students (or you!) may be able to close for us. At the time of writing, the most itching one is whether or not the Mandelbrot polynomials are unimodal. This seems like it ought to be easy, but we couldn't do it.\n",
"\n",
"The notion of iteration (discrete dynamical system) is not, apparently, of interest to every mathematician. But it is an extremely powerful idea for computation, and there are many open mathematical problems associated with it. And it is natural to introduce in a computational context.\n",
"\n",
"Rounding errors don't play a terribly large role in this unit, although they do push us to eigenvalue methods and away from polynomials that have large coefficients. This may be problematic for some students to grasp."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Mandelbrot Set and Mandelbrot Polynomials\n",
"\n",
"The [Wikipedia article on the Mandelbrot Set](https://en.wikipedia.org/wiki/Mandelbrot_set) is pretty comprehensive, and readable. We'll sum up the definition here, and try to flesh it out a bit. The definition is based on the following, apparently very simple, function:\n",
"\\begin{equation*}\n",
"f(z) = z^2 + c \n",
"\\end{equation*}\n",
"This is apparently only a simple quadratic in $z$, and depends linearly on the constant $c$. We now consider _what happens when we repeatedly apply this function_, starting at 0 (an apparently simple point). That is, put \n",
"\\begin{equation*}\n",
"z_0 = 0\n",
"\\end{equation*}\n",
"and then define $z_k$ for $k=1$, $k=2$, $k=3$, and so on in turn by the iterative formula (or recursive formula, or _dynamical_ formula)\n",
"\\begin{equation*}\n",
"z_{n+1} = f(z_n) = z_n^2 + c .\n",
"\\end{equation*}\n",
"This gives us $z_1 = z_0^2 + c = 0^2 + c = c$, and then subsequently $z_2 = z_1^2 + c = c^2 + c$, and then $z_3 = (c^2+c)^2 + c = c^4 + 2c^3 + c^2 + c$, and so on."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Pass the parcel, again\n",
"\n",
"In this unit one could play \"Pass the Parcel\" again in class, with the initiator choosing the the value of $c$ and starting with $z_0 = 0$ so $z_1 = c$. Passing each successive $z_k$ on to the next person to compute $z_{k+1}$ makes the connection with the activity from Continued Fractions. Indeed, that game could be played with Newton iteration or Halley iteration for rootfinding, and for the iterations used in the Fractals unit. It's a bit harder to find a \"pass the parcel\" game for Bohemian matrices but it could be done, if you wanted. In the next unit, on the Chaos Game Representation, there is a better game to play, but it is obviously related to this one.\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Several important things are going on here: this looks simple, and in one sense it is, but we will see some surprising things from this. This is an example of a _nonlinear dynamical system_ (perhaps the \"most famous\" example of a nonlinear dynamical system). We are repeatedly composing the function $f(z) = z^2+c$ with itself: $z_2 = f(c)$ and $z_3 = f(f(c))$ and $z_4 = f(f(f(c))$. Sometimes repeated composition is written with a power in brackets: $f^{(n)}(c)$ means $f(f(f(\\cdots(f(c))\\cdots)))$ where there are $n$ copies of $f$ used; this is where the pernicious notation $f^{(-1)}(c)$ comes from, by the way. We will try to avoid this confusing notation. Several things can be deduced about repeated composition. For instance, the result of composing two polynomials together is another polynomial. If $f(c)$ is of degree $m$ and $g(c)$ is of degree $n$ then $f(g(c))$ is of degree $mn$, and so is $g(f(c))$. If we were working with non-polynomials then we would also have to worry about matching up the domains and ranges, but polynomials are great because they always give a finite result output for any finite value input, and so composition always works.\n",
"\n",
"We have written $z_3$ out above as an explicit degree $4$ polynomial (in the variable $c$), but if we actually want to compute a numerical value of $z_3$ given a numerical value for $c$ (say, $c=-1.2$) then it turns out to be better in several ways to just use the iterative formula itself and not the explicit polynomial. This is true even though we could write that explicit polynomial somewhat more efficiently in what is known as [Horner form](https://en.wikipedia.org/wiki/Horner%27s_method): $z_3 = c\\cdot(1 + c\\cdot(1 + c\\cdot(2 + c)))$, which can be evaluated using only three floating-point multiplications and three floating-point additions. Using the iterative formula directly, though, there is no work to compute $z_0$, there is no floating-point work to compute $z_1 = c$, there is one multiplication and one addition to compute $z_2$, and one more multiplication and one more addition to compute $z_3$, showing that $z_3$ can be computed in only two \"flops\" instead of three. Note: one \"flop\" is _defined to be_ one floating-point multiplication (or division) together with one floating-point addition (or subtraction). It's not a fine-tuned measure; it was meant for use in older analyses of the cost of computation. On modern architectures, the notion is not really all that helpful; but you can see here that no matter what, the cost is less if we use the iteration directly. This advantage only increases as the iteration proceeds: to compute $z_4$ you only need one more \"flop,\" i.e. $3$ flops, whereas from its polynomial form you need $2^3 = 8$ flops (even ignoring the work that has to be done simply to write the polynomial out). For $z_{30}(c)$ the iteration takes only 29 flops; but the explicit polynomial is degree $2^{29}$ and would take $2^{29}$ flops to evaluate (actually, it's worse: we would have to use multiple precision to deal with the big integer coefficients). Thus, somehow, the iteration has _compressed_ a very high degree polynomial into a very efficient box."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(MandelbrotActivity-1)="
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Mandelbrot Activity 1\n",
":class: tip\n",
"Read the Wikipedia entry on the Mandelbrot set, and the entry on Mandelbrot. Mandelbrot's _educational_ books and papers are amazing, and amazingly clearly written. If you don't come back from this activity, well, at least you left going in a good direction.\n",
"{ref}`[Our thoughts] `\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Back to Mandelbrot iteration with a symbolic c\n",
"\n",
"Another important thing to notice in the above is that by leaving $c$ symbolic, we have turned Mandelbrot's iteration formula (or recursion formula) into some kind of generator for certain polynomials. These are called (naturally enough) _Mandelbrot polynomials_. We list the first few here, or at least their coefficients."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n",
"0.0 + 1.0·x¹\n",
"0.0 + 1.0·x¹ + 1.0·x²\n",
"0.0 + 1.0·x¹ + 1.0·x² + 2.0·x³ + 1.0·x⁴\n",
"0.0 + 1.0·x¹ + 1.0·x² + 2.0·x³ + 5.0·x⁴ + 6.0·x⁵ + 6.0·x⁶ + 4.0·x⁷ +\n",
"1.0·x⁸\n",
"0.0 + 1.0·x¹ + 1.0·x² + 2.0·x³ + 5.0·x⁴ + 14.0·x⁵ + 26.0·x⁶ + 44.0·x⁷ +\n",
"69.0·x⁸ + 94.0·x⁹ + 114.0·x¹⁰ + 116.0·x¹¹ + 94.0·x¹² + 60.0·x¹³ +\n",
"28.0·x¹⁴ + 8.0·x¹⁵ + 1.0·x¹⁶\n"
]
}
],
"source": [
"import numpy as np\n",
"from numpy.polynomial import Polynomial as Poly\n",
"N = 8\n",
"\n",
"z0 = Poly([0])\n",
"print(z0)\n",
"c = Poly([0,1]) # We can do \"symbolic polynomial arithmetic\" just with these coefficient vectors in Python\n",
"z1 = z0*z0 + c\n",
"print(z1)\n",
"z2 = z1*z1 + c\n",
"print(z2)\n",
"z3 = z2*z2 + c\n",
"print(z3)\n",
"z4 = z3*z3 + c\n",
"print(z4)\n",
"z5 = z4*z4 + c\n",
"print(z5)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Translating that Python notation for polynomials back into something standard, we have\n",
"\\begin{align}\n",
"z_0 &= 0 \\\\\n",
"z_1 &= c \\\\\n",
"z_2 &= c + c^2 \\\\\n",
"z_3 &= c + c^2 + 2c^3 + c^4\\\\\n",
"z_4 &= c + c^2 + 2c^3 + 5c^4 + 6c^5 + 6c^6 + 4c^7 + c^8 \\\\\n",
"z_5 &= c + c^2 + 2c^3 + 5c^4 + 14c^5 + 26c^6 + 44c^7 + 69c^8 + 94c^9 + 114c^{10} + 116c^{11} + 94c^{12} + 60c^{13} + 28c^{14} + 8c^{15} + c^{16} \n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Mandelbrot set is _defined to be_ the set of $c$ for which this iteration remains bounded for all $n$: that is, everything _except_ those values of $c$ for which $|z_n| \\to \\infty$. That's a weird sort of definition, and you have to first reassure yourself that there really are values of $c$ for which the $z_n$ do _not_ remain bounded. Indeed, almost the first thing you would try was $c=1$: then $z_0 = 0$, $z_1 = 1$, $z_2 = 2$, $z_3 = 5$, $z_4 = 26$, and so on; each succeeding number is one more than the square of the previous, and this goes to infinity rapidly, in some sense."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Some Activities with Mandelbrot polynomials with a more mathematical flavour"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(MandelbrotActivity-2)="
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Mandelbrot Activity 2\n",
":class: tip\n",
"Prove by induction that the degree of the $n$th Mandelbrot polynomial is $2^{n-1}$, and that the leading coefficient is $1$.\n",
"{ref}`[Our proof] `\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(MandelbrotActivity-3)="
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Mandelbrot Activity 3\n",
":class: tip\n",
" Prove that $c=-2$ and $c=0$ are both in the Mandelbrot set. Show that all zeros of Mandelbrot polynomials give _periodic orbits_ under the Mandelbrot iteration and are thus in the Mandelbrot set.\n",
"{ref}`[Our proofs] `\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(MandelbrotActivity-4)="
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Mandelbrot Activity 4\n",
":class: tip\n",
"Is $i = \\sqrt{-1}$ in the Mandelbrot set?\n",
"{ref}`[Yes] `\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(MandelbrotActivity-5)="
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Mandelbrot Activity 5\n",
":class: tip\n",
" Prove that the Mandelbrot polynomials are _unimodal_. At the time we write this, this question is open: solve it, and you could publish a paper with your proof ([Maple Transactions](https://www.mapletransactions.org) would be a good place). The word \"unimodal\" just means that the size of the coefficients increases to a peak, and then decays again. As in $z_4$, the peak might be attained by two coefficients, not just one. Seriously, we don't know how to prove this. We think it's true, though.\n",
"{ref}`[We've got nothing.] `\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"(MandelbrotActivity-6)="
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"```{admonition} Mandelbrot Activity 6\n",
":class: tip\n",
"Solve $z_3 = 0$ by hand, as follows. First, divide out the visible root, $c=0$, to get the cubic equation $0 = 1 + c + 2c^2 + c^3$. Then put $c = \\xi - 2/3$ so $c^3 = (\\xi -2/3)^3 = \\xi^3 - 2\\xi^2 + 4\\xi/9 - 8/27$ which will transform the equation into a cubic in $\\xi$ of the form $0 = q + p\\xi + \\xi^3$. Then put $\\xi = u + v$ (introducing two new variables seems like the opposite of progress, but trust us, it will help). This transforms the equation to $(u+v)^3 + p(u+v) + q = 0$, or $u^3 + 3uv(u+v) + p(u+v) + q = 0$. Show that if you can find $u$ and $v$ that simultaneously solve $3uv + p = 0$ and $u^3 + v^3 + q = 0$ then you can solve the original equation. Do so. Your formula at the end should give you three, and only three, roots. Write out your answer explicitly, and compare it with (say) the answer from Wolfram Alpha. \n",
"{ref}`[Our derivation] `\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Programming Mandelbrot polynomial iterations\n",
"Let us write a short program to evaluate Mandelbrot polynomials of arbitrary order. This code shows several features: how to define a procedure (function); how to set a type for a parameter (more of a hint than a requirement, sadly); the use of keyword arguments; and the use of comments."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# This short program will evaluate Mandelbrot polynomials at a given point,\n",
"# or else it will compute the coefficient vector of a Mandelbrot polynomial.\n",
"# Which it will do depends on the type of the argument \"c.\"\n",
"def Mandelbrot(n: int, c):\n",
" z = 0*c # try to inherit the type of c\n",
" for i in range(n):\n",
" z = z*z + c\n",
" return(z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"That code tries to force the variable `z` to have the same data type as `c` does, by multiplying 0 by `c` and storing that type of 0 in the variable `z`. It's probably not necessary, because the very first time the loop is executed `z` gets replaced by something involving `c`. But, in the case `n==0` that loop won't execute at all, and the value of `z` being returned will be 0. Let's see if that works."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n"
]
}
],
"source": [
"z = Mandelbrot(0, 1.2)\n",
"print(z)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0\n"
]
}
],
"source": [
"z = Mandelbrot(0, Poly([0,1]))\n",
"print(z)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ok, so it doesn't. This is a bit disappointing. However, it's pretty harmless in this case, so we just leave the code as it is. To be fair, deciding automatically what _type_ of zero one should return is fairly difficult. Is $0$ different from `0.0 + 0.0j`? How about a zero polynomial of grade 3, like $0 + 0x + 0x^2 + 0x^3$? Or a five hundred by five hundred matrix of zeros? The data type of computational objects really matters, in many cases. But for us, here, let's just move on.\n",
"\n",
"\n",
"Now let's experiment with the code. First, let's use it to compute a polynomial, by calling it with `c` being a polynomial of degree 1."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0 + 1.0·x¹ + 1.0·x² + 2.0·x³ + 5.0·x⁴ + 6.0·x⁵ + 6.0·x⁶ + 4.0·x⁷ +\n",
"1.0·x⁸\n"
]
}
],
"source": [
"z = Mandelbrot(4, Poly([0, 1])) # The above code works with Polys (defined in previous cell)\n",
"print(z)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.10507775999999969\n"
]
}
],
"source": [
"z = Mandelbrot(4, -1.2) # The code also works with numeric input\n",
"print(z)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"-1j\n"
]
}
],
"source": [
"z = Mandelbrot(21, 1.0j) # Will work with complex numbers also\n",
"print(z)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.17381269635430296\n"
]
}
],
"source": [
"z = Mandelbrot(30, -1.2) # Since it is iterative, we can ask for high iteration numbers on numerical input\n",
"print(z)\n",
"\n",
"# Don't try that with c = Poly([0,1]) though; the result would be a vector of length half a billion or so"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0 + 1.0·x¹ + 1.0·x² + 2.0·x³ + 5.0·x⁴ + 14.0·x⁵ + 42.0·x⁶ + 132.0·x⁷ +\n",
"365.0·x⁸ + 950.0·x⁹ + 2398.0·x¹⁰ + 5916.0·x¹¹ + 14290.0·x¹² +\n",
"33708.0·x¹³ + 77684.0·x¹⁴ + 175048.0·x¹⁵ + 385741.0·x¹⁶ + 831014.0·x¹⁷ +\n",
"1749654.0·x¹⁸ + 3598964.0·x¹⁹ + 7228014.0·x²⁰ + 14162220.0·x²¹ +\n",
"27049196.0·x²² + 50323496.0·x²³ + 91143114.0·x²⁴ + 160617860.0·x²⁵ +\n",
"275276716.0·x²⁶ + 458591432.0·x²⁷ + 742179284.0·x²⁸ + 1166067016.0·x²⁹ +\n",
"1777171560.0·x³⁰ + 2625062128.0·x³¹ + 3754272037.0·x³² +\n",
"5193067630.0·x³³ + 6939692682.0·x³⁴ + 8948546308.0·x³⁵ +\n",
"11120136162.0·x³⁶ + 13299362332.0·x³⁷ + 15286065700.0·x³⁸ +\n",
"16859410792.0·x³⁹ + 17813777994.0·x⁴⁰ + 17999433372.0·x⁴¹ +\n",
"17357937708.0·x⁴² + 15941684776.0·x⁴³ + 13910043524.0·x⁴⁴ +\n",
"11500901864.0·x⁴⁵ + 8984070856.0·x⁴⁶ + 6609143792.0·x⁴⁷ +\n",
"4562339774.0·x⁴⁸ + 2943492972.0·x⁴⁹ + 1766948340.0·x⁵⁰ + 981900168.0·x⁵¹ +\n",
"502196500.0·x⁵² + 234813592.0·x⁵³ + 99582920.0·x⁵⁴ + 37945904.0·x⁵⁵ +\n",
"12843980.0·x⁵⁶ + 3807704.0·x⁵⁷ + 971272.0·x⁵⁸ + 208336.0·x⁵⁹ +\n",
"36440.0·x⁶⁰ + 4976.0·x⁶¹ + 496.0·x⁶² + 32.0·x⁶³ + 1.0·x⁶⁴\n"
]
}
],
"source": [
"zbig = Mandelbrot(7, Poly([0, 1])) # n=8 is ghastly already; n=9 is twice as ghastly\n",
"print(zbig)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"0.0 + 1.0·x¹ + 1.0·x² + 2.0·x³ + 5.0·x⁴ + 14.0·x⁵ + 26.0·x⁶ + 44.0·x⁷ +\n",
"69.0·x⁸ + 94.0·x⁹ + 114.0·x¹⁰ + 116.0·x¹¹ + 94.0·x¹² + 60.0·x¹³ +\n",
"28.0·x¹⁴ + 8.0·x¹⁵ + 1.0·x¹⁶\n"
]
}
],
"source": [
"# If you forget the ordering, you can use keyword arguments (\"kwargs\")\n",
"zforgot = Mandelbrot(c=Poly([0,1]), n=5)\n",
"print(zforgot)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"One advantage of using Polynomials from NumPy is that they \"know how to find their own roots\". Let's test that out; the Mandelbrot polynomials turn out to be tough customers, though (which is why they are used as test problems for [MPSolve](https://en.wikipedia.org/wiki/MPSolve), not at all coincidentally—we learned about them first from the MPSolve test suite)."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfQAAAHWCAYAAACBsnu3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAU0klEQVR4nO3dX6zkZ33f8c/XXvCFFYlSG2xszFJpZWKiBtGRizc3RoHUtlAdoiSyb/CmFyuicF8jLFL1Ziv1DkFxfIEWpAraG4NVzsb8kSoTrVB8FtnE7uKysVbyai1sSOUUERlt/fTijOGwPmfP2TOzZ8bfeb2kozN/fvN7Hj/6rd6eOTO/qTFGAIA3t6sWPQEAYHaCDgANCDoANCDoANCAoANAA4IOAA3MJehV9aWqeqmqntnm/jur6pWqemr689l5jAsAbDgwp/0cT/L5JF+5xDbfG2N8bE7jAQCbzOUZ+hjjiST/MI99AQCXbz//hn5HVT1dVSeq6v37OC4AtDevl9x38oMk7xlj/Lyq7kny9SSHttqwqo4mOZok11577b963/vet09TBIDFOnXq1E/HGNfv5bE1r3O5V9XBJP9jjPE7u9j2bJLJGOOnl9puMpmM9fX1ucwPAJZdVZ0aY0z28th9ecm9qm6oqppevn067s/2Y2wAWAVzecm9qr6a5M4k11XVuSR/meQtSTLGeDjJHyf586q6kOSfktw3fM0bAMzNXII+xrh/h/s/n42PtQEAV4AzxQFAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA6yIY2unc+tDJ3Js7fSip8IVIOgAK+L4ybN59cJrOX7y7KKnwhUg6AAr4sjhg7nmwFU5cvjgoqfCFVBjjEXPYVuTyWSsr68vehoAsC+q6tQYY7KXx3qGDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADcwl6VX2pql6qqme2ub+q6nNVdaaqflhVH5zHuADAhnk9Qz+e5K5L3H93kkPTn6NJvjincQGAzCnoY4wnkvzDJTa5N8lXxobvJ3lbVd04j7EBgP37G/pNSV7YdP3c9LY3qKqjVbVeVesvv/zyvkwOAN7s9ivotcVtW34rzBjjkTHGZIwxuf7666/wtACYB9+1vnj7FfRzSd696frNSc7v09gAXGG+a33x9ivojyX5xPTd7h9K8soY48V9GhuAK8x3rS/eXL4Pvaq+muTOJNcl+UmSv0zyliQZYzxcVZXk89l4J/wvkvzZGGPHLzr3fegArJJZvg/9wDwmMMa4f4f7R5K/mMdYAMAbOVMcADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgAcJmOrZ3OrQ+dyLG104ueyq8IOgBcpuMnz+bVC6/l+Mmzi57Krwg6AFymI4cP5poDV+XI4YOLnsqv1Bhj0XPY1mQyGevr64ueBgDsi6o6NcaY7OWxnqEDQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA+zBMp7Lm9Um6AB7sIzn8ma1CTrAHizjubxZbc7lDgBLwrncAWDFCToANDCXoFfVXVX1XFWdqaoHt7j/zqp6paqemv58dh7jAgAbDsy6g6q6OskXknw0ybkkT1bVY2OM/3XRpt8bY3xs1vEAgDeaxzP025OcGWM8P8b4ZZKvJbl3DvsFAHZpHkG/KckLm66fm952sTuq6umqOlFV75/DuADA1MwvuSepLW67+LNwP0jynjHGz6vqniRfT3Joy51VHU1yNEluueWWOUwPAPqbxzP0c0neven6zUnOb95gjPGPY4yfTy+vJXlLVV231c7GGI+MMSZjjMn1118/h+kBQH/zCPqTSQ5V1Xur6q1J7kvy2OYNquqGqqrp5dun4/5sDmMDAJnDS+5jjAtV9akkjye5OsmXxhjPVtUnp/c/nOSPk/x5VV1I8k9J7hvLfIo6AHiTcepXAFgSTv0KACtO0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0W5Nja6dz60IkcWzu96KkADQg6LMjxk2fz6oXXcvzk2UVPBWhA0GFBjhw+mGsOXJUjhw8ueipAA74+FQCWhK9PBYAVJ+gA0ICgA0ADgg4AM1iWj6AKOgDMYFk+giroADCDZfkIqo+tAcCS8LE1AFhxgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANzCXoVXVXVT1XVWeq6sEt7q+q+tz0/h9W1QfnMS4AsGHmoFfV1Um+kOTuJLclub+qbrtos7uTHJr+HE3yxVnHBQB+bR7P0G9PcmaM8fwY45dJvpbk3ou2uTfJV8aG7yd5W1XdOIexmTq2djq3PnQix9ZOL9W+Os0F2Bv/jvfHPIJ+U5IXNl0/N73tcrdJklTV0apar6r1l19+eQ7TWw3HT57Nqxdey/GTZ5dqX53mAuyNf8f7Yx5Bry1uG3vYZuPGMR4ZY0zGGJPrr79+5smtiiOHD+aaA1flyOGDS7WvTnMB9sa/4/1RY2zZ1d3voOqOJP9hjPFvptc/nSRjjGObtvmrJP9zjPHV6fXnktw5xnjxUvueTCZjfX19pvkBwJtFVZ0aY0z28th5PEN/MsmhqnpvVb01yX1JHrtom8eSfGL6bvcPJXllp5gDALt3YNYdjDEuVNWnkjye5OokXxpjPFtVn5ze/3CStST3JDmT5BdJ/mzWcQGAX5s56EkyxljLRrQ33/bwpssjyV/MYywA4I2cKQ4AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAWAGx9ZO59aHTuTY2umFzkPQAWAGx0+ezasXXsvxk2cXOg9BB4AZHDl8MNccuCpHDh9c6DxqjLHQCVzKZDIZ6+vri54GAOyLqjo1xpjs5bGeoQNAA4IOAA0IOizIsrwzFuhB0GFBluWdsUAPgg4LsizvjAV68C53AFgS3uUOACtO0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGDszy4Kp6e5L/luRgkrNJ/nSM8X+22O5skv+b5P8lubDX89QCAFub9Rn6g0m+O8Y4lOS70+vb+fAY4wNiDgDzN2vQ703y5enlLyf5wxn3BwDswaxBf+cY48Ukmf5+xzbbjSTfqqpTVXV0xjEBgIvs+Df0qvpOkhu2uOszlzHO740xzlfVO5J8u6p+NMZ4YpvxjiY5miS33HLLZQwBAKtrx6CPMT6y3X1V9ZOqunGM8WJV3ZjkpW32cX76+6WqejTJ7Um2DPoY45EkjyTJZDIZO/8nAACzvuT+WJIHppcfSPKNizeoqmur6rdev5zkD5I8M+O4AMAmswb9PyX5aFX9OMlHp9dTVe+qqrXpNu9M8jdV9XSSv03yzTHGX884LgCwyUyfQx9j/CzJ729x+/kk90wvP5/kd2cZBwC4NGeKA4AGBB0AGhB0AGhA0AH24Nja6dz60IkcWzu96KlAEkEH2JPjJ8/m1Quv5fjJs4ueCiQRdIA9OXL4YK45cFWOHD646KlAkqTGWN6TsU0mk7G+vr7oaQDAvqiqU3v9VlLP0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdAC4TMt46l9BB4DLtIyn/hV0ALhMy3jqX6d+BYAl4dSvALDiBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGhB0AGhA0AGgAUEHgAYEHQAaEHQAaEDQAaABQQeABgQdABqYKehV9SdV9WxVvVZVk0tsd1dVPVdVZ6rqwVnGBGD5HFs7nVsfOpFja6cXPZWVNesz9GeS/FGSJ7bboKquTvKFJHcnuS3J/VV124zjArBEjp88m1cvvJbjJ88ueiora6agjzFOjzGe22Gz25OcGWM8P8b4ZZKvJbl3lnEBWC5HDh/MNQeuypHDBxc9lZV1YB/GuCnJC5uun0vyr7fbuKqOJjmaJLfccsuVnRkAc/Hpe347n77ntxc9jZW2Y9Cr6jtJbtjirs+MMb6xizFqi9vGdhuPMR5J8kiSTCaTbbcDAH5tx6CPMT4y4xjnkrx70/Wbk5yfcZ8AwCb78bG1J5Mcqqr3VtVbk9yX5LF9GBcAVsasH1v7eFWdS3JHkm9W1ePT299VVWtJMsa4kORTSR5PcjrJfx9jPDvbtAGAzWZ6U9wY49Ekj25x+/kk92y6vpZkbZaxAIDtOVMcADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOsCKOLZ2Orc+dCLH1k4veipcAYIOsCKOnzybVy+8luMnzy56KlwBgg6wIo4cPphrDlyVI4cPLnoqXAE1xvJ+5fhkMhnr6+uLngYA7IuqOjXGmOzlsZ6hA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQg6ADQgKADQAOCDgANCDoANCDoANCAoANAA4IOAA0IOgA0IOgA0ICgA0ADgg4ADQg6ADQwU9Cr6k+q6tmqeq2qJpfY7mxV/V1VPVVV67OMCQC80YEZH/9Mkj9K8le72PbDY4yfzjgeALCFmYI+xjidJFU1n9kAAHuyX39DH0m+VVWnquroPo0JACtjx2foVfWdJDdscddnxhjf2OU4vzfGOF9V70jy7ar60RjjiW3GO5rk9ei/WlXP7HKMVXZdEn/O2Jl12j1rtTvWafes1e7cutcH7hj0McZH9rrzTfs4P/39UlU9muT2JFsGfYzxSJJHkqSq1scY277Zjg3WaXes0+5Zq92xTrtnrXZnljeOX/GX3Kvq2qr6rdcvJ/mDbLyZDgCYk1k/tvbxqjqX5I4k36yqx6e3v6uq1qabvTPJ31TV00n+Nsk3xxh/Pcu4AMBvmvVd7o8meXSL288nuWd6+fkkv7vHIR7Z++xWinXaHeu0e9Zqd6zT7lmr3dnzOtUYY54TAQAWwKlfAaCBpQl6Vf3nqvpRVf2wqh6tqrdts93Kn0b2Mtbqrqp6rqrOVNWD+zzNhXNq4t27jLVa9WPq7VX17ar68fT3P9tmu5U8pnY6PmrD56b3/7CqPriIeS6DXazVnVX1yvQYeqqqPrvTPpcm6Em+neR3xhj/Msn/TvLpS2z74THGB1b4IxA7rlVVXZ3kC0nuTnJbkvur6rZ9neXivX5q4i0/InmRVT+mdlwrx1SS5MEk3x1jHEry3en17azUMbXL4+PuJIemP0eTfHFfJ7kkLuPf0vemx9AHxhj/caf9Lk3QxxjfGmNcmF79fpKbFzmfZbbLtbo9yZkxxvNjjF8m+VqSe/drjstgjHF6jPHcoufxZrDLtVr5Yyob/71fnl7+cpI/XNxUls5ujo97k3xlbPh+krdV1Y37PdElcEX+LS1N0C/y75Kc2OY+p5H9Tdut1U1JXth0/dz0Nt7IMbU7jqnknWOMF5Nk+vsd22y3isfUbo4Px9CG3a7DHVX1dFWdqKr377TTWb9t7bLs5jSyVfWZJBeS/NdtdrPr08i+mc1hrbb6xpx2H2nY71MTv5nNYa1W/pi6jN2sxDF1kd0cHytxDO3CbtbhB0neM8b4eVXdk+Tr2fhTxbb2Neg7nUa2qh5I8rEkvz+2+Tzd5ZxG9s1sDmt1Lsm7N12/Ocn5+c1wOez3qYnfzOawVit/TFXVT6rqxjHGi9OXil/aZh8rcUxdZDfHx0ocQ7uw4zqMMf5x0+W1qvovVXXdpb6GfGlecq+qu5L8+yT/dozxi222cRrZ7G6tkjyZ5FBVvbeq3prkviSP7dcc3ywcU5fFMbXx3/vA9PIDSd7wysYKH1O7OT4eS/KJ6bvdP5Tkldf/hLFidlyrqrqhauO7yavq9mz0+meX3OsYYyl+kpzJxt8Unpr+PDy9/V1J1qaX/0WSp6c/z2bjpcKFz30Z12p6/Z5svAv+71dxrZJ8PBv/J/xqkp8kedwxtfe1ckyNJPnn2Xh3+4+nv9/umPqN9XnD8ZHkk0k+Ob1c2Xh3998n+bskk0XPeYnX6lPT4+fpbLz5+fBO+3SmOABoYGlecgcA9k7QAaABQQeABgQdABoQdABoQNABoAFBB4AGBB0AGvj/I9uPFQaI1b0AAAAASUVORK5CYII=\n",
"text/plain": [
"