{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Symbolic Computation: The Pitfalls\n", "\n", "This collection of notebooks is mostly numerical, with not a lot of exact or symbolic computation, with the exception of polynomial computation done by computing with vectors of monomial basis coefficients, so it's really numerical anyway. Why not do symbolic computation? And, for that matter, why is numerical computing (even with all the unexpected behaviour of floating-point arithmetic) so much more popular than symbolic or exact computing?\n", "\n", "This section explores symbolic computation and its pitfalls. We do so from the point of view of experience and with some authority: we have used symbolic computation (usually in Maple, but also in other symbolic languages) for years (decades!) and know it and its benefits well. Caveat: We do not know SymPy so well, and so if we say that SymPy can't do something, we may well be wrong. One of us knows Sage quite well, but we're not using Sage here (yet).\n", "\n", "One of RMC's earliest co-authors, Honglin Ye, put it well when he suggested that not everyone needs numerical methods but that everyone could use symbolic computation.\n", "\n", "Wait. Isn't that contradictory? If everyone could use it, why aren't they?\n", "\n", "There are, we think, a few main obstacles.\n", "\n", "1. Symbolic computation systems are hard to learn how to use well, because there's a lot to learn (indeed, you kind of have to know the math first, too). Look at the SymPy Tutorial for example. It has ten sections, one labeled \"Gotchas\". The SAGEMATH system, which also works with Python, is both more powerful and more complicated: See the SAGEMATH Tutorial to get started there.\n", "2. Some mathematical problems are inherently too expensive to solve in human lifetimes, even with today's computers, and people unfairly blame symbolic computation systems for this.\n", "3. Even if you can solve a problem exactly, with extra effort, that effort might be wasted because the approximate answers are also the \"exact\" answers to similar problems, and those similar problems might be just as good a model of whatever system you were trying to understand. This is especially true if the data is only known approximately.\n", "4. \"Symbolic Computation\" and \"Computer Algebra\" are related terms—about as close as \"Numerical Analysis\" and \"Computational Science\" if that comparison means anything—but the differences are remarkably important, because what gets implemented is usually a Computer Algebra system, whereas what people actually want to use is a symbolic computation system. We'll show you what that means.\n", "5. Symbolic computation systems are hard to implement well. The major systems (Maple, Mathematica, and Matlab) charge money for their products, and get what they ask for; this is because their systems are better than the free ones in many respects, because they have invested significant programmer time to address the inherent difficulties. Free systems, such as SymPy, will do the easy things for you; and we will see that they can be useful. But in reality there's no comparison (although we admit that the SAGEMATH people may well disagree with our opinion).\n", "\n", "All that said, symbolic computation can be extremely useful (and interesting), and is sometimes worth all the bother. Let's look first at what Python and SymPy can do. Later we'll look at what the difficulties are." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "100 factorial is 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000\n", "The floating point value of p is 9.332621544394415e+157\n" ] } ], "source": [ "n = 100\n", "p = 1\n", "for i in range(n): \n", " p = p*(i+1)\n", "print(n, ' factorial is ', p)\n", "print('The floating point value of p is ', 1.0*p)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first thing we see is that Python has, built-in, arbitrary precision integer arithmetic. Yay?" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "720 factorial is 2601218943565795100204903227081043611191521875016945785727541837850835631156947382240678577958130457082619920575892247259536641565162052015873791984587740832529105244690388811884123764341191951045505346658616243271940197113909845536727278537099345629855586719369774070003700430783758997420676784016967207846280629229032107161669867260548988445514257193985499448939594496064045132362140265986193073249369770477606067680670176491669403034819961881455625195592566918830825514942947596537274845624628824234526597789737740896466553992435928786212515967483220976029505696699927284670563747137533019248313587076125412683415860129447566011455420749589952563543068288634631084965650682771552996256790845235702552186222358130016700834523443236821935793184701956510729781804354173890560727428048583995919729021726612291298420516067579036232337699453964191475175567557695392233803056825308599977441675784352815913461340394604901269542028838347101363733824484506660093348484440711931292537694657354337375724772230181534032647177531984537341478674327048457983786618703257405938924215709695994630557521063203263493209220738320923356309923267504401701760572026010829288042335606643089888710297380797578013056049576342838683057190662205291174822510536697756603029574043387983471518552602805333866357139101046336419769097397432285994219837046979109956303389604675889865795711176566670039156748153115943980043625399399731203066490601325311304719028898491856203766669164468791125249193754425845895000311561682974304641142538074897281723375955380661719801404677935614793635266265683339509760000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000\n" ] } ], "source": [ "n = 720\n", "p = 1\n", "for i in range(n): \n", " p = p*(i+1)\n", "print(n, ' factorial is ', p)\n", "# print('The floating point value of p is ', 1.0*p) # Causes OverflowError" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Large integers cost more to manipulate, and the above number is pretty long. But SymPy will do it if you ask. One thing you might want to do is factor those numbers. Or one might just want to know the prime factors." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The prime factors of 720 : [2, 3, 5]\n", "The prime factors of 2601218943565795100204903227081043611191521875016945785727541837850835631156947382240678577958130457082619920575892247259536641565162052015873791984587740832529105244690388811884123764341191951045505346658616243271940197113909845536727278537099345629855586719369774070003700430783758997420676784016967207846280629229032107161669867260548988445514257193985499448939594496064045132362140265986193073249369770477606067680670176491669403034819961881455625195592566918830825514942947596537274845624628824234526597789737740896466553992435928786212515967483220976029505696699927284670563747137533019248313587076125412683415860129447566011455420749589952563543068288634631084965650682771552996256790845235702552186222358130016700834523443236821935793184701956510729781804354173890560727428048583995919729021726612291298420516067579036232337699453964191475175567557695392233803056825308599977441675784352815913461340394604901269542028838347101363733824484506660093348484440711931292537694657354337375724772230181534032647177531984537341478674327048457983786618703257405938924215709695994630557521063203263493209220738320923356309923267504401701760572026010829288042335606643089888710297380797578013056049576342838683057190662205291174822510536697756603029574043387983471518552602805333866357139101046336419769097397432285994219837046979109956303389604675889865795711176566670039156748153115943980043625399399731203066490601325311304719028898491856203766669164468791125249193754425845895000311561682974304641142538074897281723375955380661719801404677935614793635266265683339509760000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 : [2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251, 257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389, 397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541, 547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677, 683, 691, 701, 709, 719]\n" ] } ], "source": [ "from sympy import primefactors \n", "\n", "primefactors_n = primefactors(n) \n", "print(\"The prime factors of {} : {}\".format(n, primefactors_n)) \n", "\n", "primefactors_p = primefactors(p) \n", "print(\"The prime factors of {} : {}\".format(p, primefactors_p)) \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Factoring seems like such a simple problem, and it's so natural to have it implemented in a symbolic computation system. The number 720! is 1747 digits long. Maybe all 1700--odd digits long integers are so easy to factor?\n", "\n", "Um, no. See the discussion [here](https://en.wikipedia.org/wiki/Integer_factorization) to get started. Let's take a modest problem and time it here." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The prime factors of 3000000000238000000004719 : {1000000000039: 1, 3000000000121: 1}\n", "--- 0.9057760238647461 seconds ---\n" ] } ], "source": [ "funny = 3000000000238000000004719\n", "#notfunny = 45000000000000000057990000000000000024761900000000000003506217\n", "from sympy import factorint\n", "import time\n", "start_time = time.time()\n", "factordict = factorint(funny) \n", "print(\"The prime factors of {} : {}\".format(funny, factordict)) \n", "print(\"--- %s seconds ---\" % (time.time() - start_time))\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That factoring of $3000000000238000000004719$ __with the previous version of SymPy__ took between 8 and 11 seconds on this machine (different times if executed more than once) but has since been improved to be just under a second; on this very same machine, Maple's \"ifactor\" command succeeds so quickly that it registers no time taken at all, possibly because it is using a very specialized method; factoring integers is an important feature of symbolic computation and Maple's procedures for it have been a subject of serious research for a long time. Maple's help pages cite three important papers, and tell you that it uses an algorithm called the quadratic sieve. Maple can factor $45000000000000000057990000000000000024761900000000000003506217$ into its three prime factors in about 7.5 seconds on this machine; in contrast, after fifty minutes running trying to factor that with factorint as above, RMC had to hard-restart to get Python's attention. (Haven't tried with the new version, though.)\n", "\n", "That SymPy takes so long to factor integers, in comparsion, suggests that it still isn't using the best methods (the documentation says that it switches between three methods, trial division, Pollard rho, and Pollard p-1) ; and because factoring is such a basic algorithm (an even more basic one is GCD or Greatest Common Divisor) this will have important knock-on effects.\n", "\n", "But factoring, as old an idea as it is, is complicated enough to be used as a basic idea in modern cryptography. The slowness of SymPy is not completely its fault: the problem is hard.\n", "\n", "Let's move on to computing with functions. As previously stated, most supposedly \"symbolic\" systems are really \"algebra\" systems: this means that they work well with polynomials (even multivariate polynomials). A polynomial considered as an algebraic object is isomorphic to a polynomial considered as a function, but the difference in viewpoint can alter the affordances. What the word \"affordance\" means in this context is that \"something can happen with it\": for instance, you can pick out a lowest-degree term; or you can add it to another polynomial; or you can square it; and so on. As a function, you can evaluate it at a particular value for the symbols (variables)." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left\\{- \\sqrt{3}, \\sqrt{3}\\right\\}$" ], "text/plain": [ "{-sqrt(3), sqrt(3)}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from sympy import *\n", "x = symbols('x')\n", "\n", "solveset(Eq(x**2, 3), x)\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left\\{- \\frac{1}{3 \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}} + \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}, - \\frac{\\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}}{2} + \\frac{1}{6 \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}} + i \\left(\\frac{\\sqrt{3}}{6 \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}} + \\frac{\\sqrt{3} \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}}{2}\\right), - \\frac{\\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}}{2} + \\frac{1}{6 \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}} + i \\left(- \\frac{\\sqrt{3} \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}}{2} - \\frac{\\sqrt{3}}{6 \\sqrt{\\frac{1}{2} + \\frac{\\sqrt{93}}{18}}}\\right)\\right\\}$" ], "text/plain": [ "{-1/(3*(1/2 + sqrt(93)/18)**(1/3)) + (1/2 + sqrt(93)/18)**(1/3), -(1/2 + sqrt(93)/18)**(1/3)/2 + 1/(6*(1/2 + sqrt(93)/18)**(1/3)) + I*(sqrt(3)/(6*(1/2 + sqrt(93)/18)**(1/3)) + sqrt(3)*(1/2 + sqrt(93)/18)**(1/3)/2), -(1/2 + sqrt(93)/18)**(1/3)/2 + 1/(6*(1/2 + sqrt(93)/18)**(1/3)) + I*(-sqrt(3)*(1/2 + sqrt(93)/18)**(1/3)/2 - sqrt(3)/(6*(1/2 + sqrt(93)/18)**(1/3)))}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solveset(Eq(x**3+x-1, 0), x)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--- 0.0 seconds ---\n" ] } ], "source": [ "start_time = time.time()\n", "#solveset(Eq(x**4+x-1, 0), x) # Interrupted after about two hours: the code did not succeed\n", "print(\"--- %s seconds ---\" % (time.time() - start_time))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "__Updated__ With the new version of SymPy, that computation succeeds in under 3 seconds. Got the right answers, too.\n", "\n", "Nevertheless. In those original two hours, RMC went and had his dinner; then downloaded a paper by Dave Auckly from the American Mathematical Monthly 2007 which talks about solving the quartic with a pencil (an algebraic geometer's pencil!), read it, and solved the problem by hand, including solving the resolvent cubic by hand, which he already knew how to do. And got it right, too. So there.\n", "\n", "In contrast, Maple (nearly instantaneously) returns—if you force it to by saying you want the explicit solution—the answer\n", "\n", "$$\n", "\\frac{\\sqrt{6}\\, \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{12}+\\frac{\\mathrm{I} \\sqrt{6}\\, \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}+12 \\sqrt{6}\\, \\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}-48 \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}}}{12}\n", ", \n", "\\frac{\\sqrt{6}\\, \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{12}-\\frac{\\mathrm{I} \\sqrt{6}\\, \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}+12 \\sqrt{6}\\, \\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}-48 \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}}}{12}\n", ", \n", "-\\frac{\\sqrt{6}\\, \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{12}+\\frac{\\sqrt{6}\\, \\sqrt{\\frac{-\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}+12 \\sqrt{6}\\, \\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}+48 \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}}}{12}\n", ", \n", "-\\frac{\\sqrt{6}\\, \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{12}-\\frac{\\sqrt{6}\\, \\sqrt{\\frac{-\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}+12 \\sqrt{6}\\, \\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}+48 \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}} \\sqrt{\\frac{\\left(108+12 \\sqrt{849}\\right)^{\\frac{2}{3}}-48}{\\left(108+12 \\sqrt{849}\\right)^{\\frac{1}{3}}}}}}}{12}\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is an example of what Velvel Kahan calls a \"wallpaper expression.\" He also famously said, \"Have you ever asked a computer algebra system a question, and then, as the screensful of answer whizzed past your eyes, said \"I wish I hadn't asked?\"\"\n", "\n", "The use of that exact answer (quickly obtained or not) is questionable. Then, of course, the Abel-Ruffini theorem says that there is no general formula for solving polynomials of degree $5$ or higher in terms of radicals. For degree $5$ polynomials, there is a solution in terms of elliptic functions; again, it's complicated enough that it's of questionable use. Then there is Galois theory which describes the algebraic structures of polynomials. See the interesting historical essay by Nick Trefethen on What we learned from Galois .\n", "\n", "The lesson here is that even when you can solve something exactly, maybe you shouldn't. \n", "\n", "There are some interesting things you can do with univariate polynomials of high degree, including with the algebraic numbers that are their roots. But computation with them isn't so easy. SymPy actually has some quite advanced features for polynomials, including multivariate polynomials." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's reinforce that lesson (\"Even if you can, maybe you shouldn't\") with a little linear algebra. Some of the material of this OER uses determinants of matrices, and introductory textbooks give the formula\n", "\n", "$$\n", "\\det \\left[\\begin{array}{cc}\n", "a_{1,1} & a_{1,2} \n", "\\\\\n", " a_{2,1} & a_{2,2} \n", "\\end{array}\\right]\n", " = a_{1,1} a_{2,2}-a_{1,2} a_{2,1}\n", "$$\n", "\n", "usually with the simpler notation where the matrix is has entries $a$, $b$, $c$, $d$ and the determinant is $ad-bc$.\n", "\n", "The formula for the three-by-three case is more obnoxious:\n", "\n", "$$\n", "\\det \\left[\\begin{array}{ccc}\n", "a_{1,1} & a_{1,2} & a_{1,3} \n", "\\\\\n", " a_{2,1} & a_{2,2} & a_{2,3} \n", "\\\\\n", " a_{3,1} & a_{3,2} & a_{3,3} \n", "\\end{array}\\right] = a_{1,1} a_{2,2} a_{3,3}-a_{1,1} a_{2,3} a_{3,2}-a_{1,2} a_{2,1} a_{3,3}+a_{1,2} a_{2,3} a_{3,1}+a_{1,3} a_{2,1} a_{3,2}-a_{1,3} a_{2,2} a_{3,1}\n", "$$\n", "\n", "This six-term formula, which can be rearranged to make it less costly to evaluate, cannot be made (much) simpler for humans to understand.\n", "\n", "The four-by-four case has twenty-four terms in its determinant, and requires effort to read. Again it can be rearranged to make it less costly to evaluate, but cannot be made much simpler for humans to understand.\n", "\n", "The _general_ formula looks pretty simple, however:\n", "\n", "$$\n", "\\det \\mathbf{A} = \\sum_{\\sigma \\in S_n} \\mathrm{sgn}(\\sigma) \\prod_{i=1}^n a_{i,\\sigma_i}\n", "$$\n", "\n", "The difficulty comes in unpacking the notation. The outer sum is over all $n!$ permutations of the integers $1$, $2$, $\\ldots$, $n$ (so when $n=4$ there are $4! = 24$ terms, when $n=5$ there are $5!=120$ terms, etc). One has to be able to determine the _sign_ of a permutation as well. The factorial number of terms is pretty inescapable.\n", "\n", "Computer algebra systems will quite happily give you these determinantal formulas, written out. But even for a ten-by-ten matrix (which is pretty small in today's terms) there will be $10! = 3,628,800$ terms in the answer, which will occupy several screens. Evaluation of that formula, given numerical values for the $a_{i,j}$, has to consider each of those (more than three million) terms.\n", "\n", "Now, sometimes things can be done with those symbolic expressions. We mentioned above that they can be rearranged to be less costly to evaluate. Here is a \"straight-line program\" to evaluate the five-by-five determinant, found by Maple's codegen[optimize] command.\n", "\n", "\n", " t2 := a[5, 4]; \n", " t22 := a[1, 4]; \n", " t23 := a[1, 3]; \n", " t24 := a[1, 2]; \n", " t8 := a[4, 3]; \n", " t9 := a[4, 2]; \n", " t28 := t23*t9-t24*t8; \n", " t3 := a[5, 3]; \n", " t4 := a[5, 2]; \n", " t29 := t23*t4-t24*t3; \n", " t38 := t3*t9-t4*t8; \n", " t7 := a[4, 4]; \n", " t47 := t2*t28-t22*t38-t29*t7; \n", " t12 := a[3, 4]; \n", " t13 := a[3, 3]; \n", " t14 := a[3, 2]; \n", " t30 := t22*t9-t24*t7; \n", " t31 := t22*t8-t23*t7; \n", " t46 := t12*t28-t13*t30+t14*t31; \n", " t17 := a[2, 4]; \n", " t18 := a[2, 3]; \n", " t19 := a[2, 2]; \n", " t45 := (t18*t9-t19*t8)*t12-(t17*t9-t19*t7)*t13+(t17*t8-t18*t7)*t14; \n", " t39 := t2*t9-t4*t7; \n", " t40 := t2*t8-t3*t7;\n", " t44 := t12*t38-t13*t39+t14*t40; \n", " t43 := t17*t28-t18*t30+t19*t31; \n", " t33 := t2*t24-t22*t4; \n", " t34 := t2*t23-t22*t3; \n", " t42 := t17*t29+t18*t33-t19*t34; \n", " t41 := t17*t38-t18*t39+t19*t40; \n", " t25 := a[1, 1]; \n", " t5 := a[5, 1]; \n", " t32 := t2*t25-t22*t5; \n", " t27 := t23*t5-t25*t3; \n", " t26 := t24*t5-t25*t4; \n", " t20 := a[2, 1]; \n", " t15 := a[3, 1]; \n", " t10 := a[4, 1]; \n", " t1 := (-t42*t15+(t17*t27+t18*t32-t20*t34)*t14+(-t17*t26-t19*t32+t20*t33)*t13\n", " +(t18*t26-t19*t27+t20*t29)*t12)*a[4, 5]+(-t45*t5+t44*t20-t41*t15+((t17*t3-t2*t18)*t14\n", " +(-t17*t4+t19*t2)*t13+(t18*t4-t19*t3)*t12)*t10)*a[1,5]\n", " +(t46*t5-t44*t25-t47*t15+(-t12*t29-t13*t33+t14*t34)*t10)*a[2, 5]\n", " +(t10*t42+t20*t47+t25*t41-t43*t5)*a[3, 5]\n", " +(t45*t25-t46*t20+t43*t15+((-t17*t23+t22*t18)*t14+(t17*t24-t19*t22)*t13\n", " +(-t18*t24+t19*t23)*t12)*t10)*a[5, 5]\n", "\n", "(that was cut-and-pasted in and hand-edited for clarity: hopefully no typos were introduced). The cost of evaluating the original big ugly expression with $120$ terms is 119 additions + 480 multiplications + 600 subscripts.\n", "The cost of evaluating that _straight line program_ above (generated by computer algebra) is only 25 subscripts + 40 assignments + 106 multiplications + 66 additions, in comparison; this is a significant reduction in cost. Also, there were no _divisions_ introduced, which normal floating-point operations do. So, this compression is possibly of interest.\n", "\n", "But sometimes no compression is possible, and the only thing one can say about a huge symbolic answer is \"It is what it is.\"\n", "\n", "The _general_ case of optimal algorithms for linear algebra even over the integers, never mind with symbols, is a topic of much current study. We do not even know the cheapest way to multiply two matrices! And that algorithm gives a key to almost every other algorithm. These problems are _hard_." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Symbolic computation with functions\n", "\n", "Let's try some calculus-like things." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left\\{i \\left(2 n \\pi + \\arg{\\left(x \\right)}\\right) + \\log{\\left(\\left|{x}\\right| \\right)}\\; \\middle|\\; n \\in \\mathbb{Z}\\right\\}$" ], "text/plain": [ "ImageSet(Lambda(_n, I*(2*_n*pi + arg(x)) + log(Abs(x))), Integers)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "y = symbols('y')\n", "\n", "solveset(Eq(exp(y), x), y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "RMC really doesn't like that \"solution\"! It has separated the real and imaginary parts without needing to. A perfectly good answer would be $\\ln_n(x)$, which looks a lot simpler.\n", "\n", "$\\ln_k(z)$, which might not look familiar to you, means $\\ln(z) + 2\\pi i k$. \n", "\n", "Fine. We will live with it. The solution is not actually wrong.\n", "\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\left\\{y\\; \\middle|\\; y \\in \\mathbb{C} \\wedge - x + y e^{y} = 0 \\right\\}$" ], "text/plain": [ "ConditionSet(y, Eq(-x + y*exp(y), 0), Complexes)" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solveset(Eq(y*exp(y), x), y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oh, that's disappointing. See the Wikipedia article on Lambert W to see what that should have been." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[LambertW(x)]" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "solve(Eq(y*exp(y), x), y)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's better, but—like the logarithm above—there should be multiple branches." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Integrals, and the difference between Computer Algebra and Symbolic Computation\n", "\n", "We'll start with a nasty example. You can find nice examples of SymPy and integration in many places, so we will assume that you have seen instances of computer implementations of the fundamental theorem of calculus: to find areas under curves by using antiderivatives. The nasty integral that we will use is, if $x \\ne 0$,\n", "\n", "$$\n", "f(x) = \\frac{e^{-1/x}}{x^2(1+e^{-1/x})^2}\n", "$$\n", "\n", "with the removal of the discontinuity at $x=0$ by defining $f(0)=0$.\n", "\n", "We will try to use both NumPy and SymPy to integrate this (infinitely smooth! Just check the derivatives at zero. The function is _infinitely flat_ there) function on various intervals in the real axis. The function has $f(x)\\ge 0$, and is zero only at $x=0$. Therefore the integral of $f(x)$ from any $a$ to any $b > a$ will be positive. Essentially positive functions have positive area underneath them, end of story. See the figure below. We do avoid computing the function at exactly zero; but that's just laziness on our part." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "from matplotlib import pyplot as plt" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Area under the curve from -1 to 1 is approximately 0.5371959827682197\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD8CAYAAACfF6SlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAZ5klEQVR4nO3dfYxcV3nH8e8Tv4QkkLCJI2piy3HaNcGYysSbOG0l89KFbFpqUxHAqSBOa8tAcbGKQDVCRCgRooAEcktUsJKUmKo4IRXqBoxdTJJGVYm7u63BL8H2JmC8xiV2srELjp3d5Okfc8dcj2d2Znfunftyfh/J8sydOzPPnDvz7JnnnnPG3B0REQnLBVkHICIinafkLyISICV/EZEAKfmLiARIyV9EJEBK/iIiAWop+ZtZn5ntN7NhM9tQ5/bbzeyYme2K/q1JPlQREUnK9GY7mNk04G7g7cAIMGBm/e6+r2bXB9x9XQoxiohIwlrp+d8ADLv70+7+IrAFWJFuWCIikqamPX/gKuBw7PoIsLTOfu82s2XAAeCv3f1w7Q5mthZYC3DJJZcsufbaaycfsYhIwIaGho67+5XtPk4ryb8VDwPfdPczZvZB4H7gbbU7ufsmYBNAT0+PDw4OJvT0IiJhMLNDSTxOK2WfI8Dc2PU50baz3P1Zdz8TXb0HWJJEcCIiko5Wkv8A0G1m881sJrAS6I/vYGazY1eXA08mF6KIiCStadnH3cfNbB2wHZgG3Ofue83sTmDQ3fuBj5rZcmAceA64PcWYRUSkTZbVks6q+YuITJ6ZDbl7T7uPoxm+IiIBUvIXEQmQkr+ISICU/EVEAqTkLyISICV/EZEAKfmLiARIyV9EJEBK/iIiAVLyFxEJkJK/iEiAlPxFRAKk5C8iEiAlfxGRACn5i4gESMlfRCRASv4iIgFS8hcRCZCSv4hIgJT8RUQCND3rAKSzhg6NsnHHAfoWzebBgZ+DGe/tmcu2PUdZ37uAJfO6sg5RcmTo0Ch3Pbz3nPdJ36LZer+UgJJ/AOIJ/4vbf8LoqTF2HznB6KkxAA49+2tGT41x8oUxLr1ohj7UcvY9c/L0OLtGTgC/eZ9U3zu7j5zgEzddqz8EBWXunskT9/T0+ODgYCbPHZKhQ6OsuX+A0VNjdF084+z/n7jp2vN6/idPj7Pr8PN0XTyDe1Zdrw9zoOLvmcVzLqtsrOn5VzsR8feU3jOdYWZD7t7T9uMo+ZfTOT23KKE366XVfuj1LSAs9d4zjRJ6vW+Tes90hpK/TOi2e3fy+MHjk/5A1iaAZd2z2Lx6aQcilqzpPVMMSSV/1fxLJt4jAybdC1syr4vNq5eePdF38vQ4Q4dG1ZMruaFDo5x8YYzFc1/Np9+5cMrvmep777Z7d+obQM5pqGeJVMs2jx88zrY9R9m8eumUP3xL5nVx6UUz2HX4eTbuOJBwpJI3G3ccYNfICS59xfS23jObVy9l256jPH7wOGvuH2Do0GjCkUpSlPxLZOOOA2dPvq3vXdD2463vXcCy7llne3L6IJfP0KFRbrt3J32LZrOse1Zi75vqiWB1HPJLyb8Eaj/ASY26qO3J6YNcPht3HEjkm2Lcknld3LPqenUcck41/4KLj9ABUjnRVu0NJtErlHxJ69hWOw7Vk8i7j5zQUNCcUc+/4JIu9dRT/SAD6sWVRPXbIpBYj78elYDyS8m/oNIq9UykWiLQh7j4OnUsVQLKL5V9Cqr64YV0Sj31qPxTHp08lrUlIOjce1YaU8+/YNIYndGqJfO6WN+7gI07Dqj3VmDV8fidHoev0WP5op5/wWTR48/T80v7sjqG+gaQLy31/M2sz8z2m9mwmW2YYL93m5mbWdtTj+V88VmYWZVeqr03lX6KK+tjuL53AYvnXHZ29rhko2nyN7NpwN3AzcBC4FYzW1hnv1cB64GdSQcpFUnMwmxXtfemIXvFlfUx1OzxfGil538DMOzuT7v7i8AWYEWd/e4CPg+cTjA+Ids6f7OY1HMrjjwdM9X/s9dK8r8KOBy7PhJtO8vMrgPmuvt3J3ogM1trZoNmNnjs2LFJBxuqNGZhJhWTem7Fkadjptnj2Wv7hK+ZXQB8Cbi92b7uvgnYBJUlndt97hDkoc5fj4Z9Fk8ej9n63gWcfGFMq8dmoJWe/xFgbuz6nGhb1auARcBjZvYz4EagXyd9k5GHOn89WdeNZfLyeMxU/89OK8l/AOg2s/lmNhNYCfRXb3T3E+4+y92vdvergSeA5e6uX2ppQx7r/LXyVEOWieX5WKn+n42mZR93HzezdcB2YBpwn7vvNbM7gUF375/4EWQqijCevggxSkWej5XG/2ejpZq/u28FttZsu6PBvm9pP6xw1fslrrzKYw1Z6ivCsVL9v7P0G745U+396HdQJUR6/zeX1G/4am2fHMnryB6RTlH9v3OU/HMkryN7msnzycTQFe3YaPx/52hhtxwpQl22njyfTAxdUY9NUT8LRaKefw506leV0pL1QmHSWFGPjX49Ln3q+edAUXtnVfEPquRL0Y9N0T8beaaef8Z0klekMS3/nB4l/4wV9SRvraKdWCy7shwPLf+QHiX/jBRh+YbJyNOKkVKu46Hhn+lQzT8jZatlanRGvpTpeGj5h3Qo+WegjHX+op9YLJsyHg8t/5AslX0yUJY6v0gnqf6fLPX8O6hIi7aJ5FH1M1Ot/6/vXaAO1BSp599Befw5xqSVZZRJkZX5GGj5h+So599BZToJ10jZTmQXUQjHIITPUtrU8++Aoi/fMBlFXU6gTEI4Blr+oX3q+XdACD2xqjKOMimakI5BSJ+tpKnnn7IyDusUyQst/zB1Sv4p07BOkfRo+OfUqeyTEg3rFOkMDf+cGvX8UxLCsM6JlHm4YV6F2uYa/jk1Sv4pUJ2/XAuLFUXoba76/+Qo+adAdf4whhvmTehtrvr/5Ji7Z/LEPT09Pjg4mMlzpyVe59+256hqjyIdFsJn0MyG3L2n3cfRCd8EacyxSLa0/HPrlPwTpCnnIvmgz2JzqvknIKTlG0SKQMs/NKfkn4DQR1k0EurQwyyorevTZ7MxJf82aVhnY/rgdY7auj4N/2xMyb9NGtbZWOhDDztJbV2fhn82pqGeUxTCkDKRMijbZ1VDPTOmYZ0ixaDhn/Wp7DMFqvOLFI/q/+dS8p8C1flFikf1/3O1lPzNrM/M9pvZsJltqHP7h8xst5ntMrP/MLOFyYeavepwur5Fs3VyTaSAqifGq8s/h/wNoGnyN7NpwN3AzcBC4NY6yf2f3f2N7r4Y+ALwpaQDzYPQl2meKo1BT4/adnK0/PNvtNLzvwEYdven3f1FYAuwIr6Du5+MXb0EyGYIUYpU5586jUFPj9p2alT/by35XwUcjl0fibadw8w+YmZPUen5f7TeA5nZWjMbNLPBY8eOTSXezKjOP3Uag54ete3UqP7fwjh/M7sF6HP3NdH1DwBL3X1dg/3/DLjJ3VdN9LhFGedftjHCIlJR1M92J8f5HwHmxq7PibY1sgX4h3aCyhON5xcpp9DH/7dS9hkAus1svpnNBFYC/fEdzKw7dvWPgYPJhZgNjewRCUOoI4Ca9vzdfdzM1gHbgWnAfe6+18zuBAbdvR9YZ2a9wBgwCkxY8ikC9fhFwhDqN4CWlndw963A1pptd8Qur084rkxpZI9IeNb3LuDkC2NnRwAVof7fDs3wrUMje5Kn8ejJU5smK7QRQEr+Marzp0fj0ZOnNk1eSPV/reoZozp/evSbqslTmyYvpPq/kn9Edf50xX9TVZKhNk1PCPV/lX0iqvOLSFUI9f/ge/7xWX6gr9AiUlHNBdX6f1FmALcq2ORfTfonT4+z6/DzQHlreyIyebX1/5MvjHHpRTNK80cg2LLP2ZO77hrZIyINVUcAYVaq0VVB9vzjJ3c//c6FpfgrLiLpqH4DGDo0yl0P7y3NSeAge/46uZsdTUxqn9owG2U7CRxU8tckruxpYlL71IbZKdMksKDKPprElT1NTGqf2jA7ZZoEFkzy1ySufNDEpPapDbNXhklgpU/+5wzpHDnBsu5ZhTxQIpIf1fr/4wePc9fDews5BLT0NX8N6RSRNBR9CGhpe/71Zu4W6a+yiORb0YeAlrbnX+3xb9tzlM2rlxbmgIhIsRR1CGgpk79O7opIJxVxCGipkn91HP9d39mnSVw5p4lKk6c2y69qCWjbnqNnTwLn/ViVKvnr5G5xaKLS5KnN8q9IJ4FLc8JX6/UUiyYqTZ7aLP+KdBK4ND1/rddTLNUPiY5V69RmxVGEk8CF7/nrx1hEJI/y/mMwhU7+Q4dGWXP/AKOnxoDirrEhIuVTuw7Q7iMnuGfV9bn5A1Doss/GHQcYPTVG18Uz1OMXkVxa37uArotnMHpqLFcloEIm/9qlmfP011REJG7JvC7uWXV97uYBFKrso9/dFZEiyuPvAReq569x/CJSZHmaB1CY5H/OOP4/eYOGvJWAZqy2Tm1VDtVvAJ9+50IWz7ns7DyALOS+7KP1+MtLv6zWOrVVueTh9wBynfzjQzkXz7lMpZ6S0YzV1qmtyqd6LE+eHs9kKKi5e0eeqFZPT48PDg5OuE/15EjXxTM0okdESineyV3WPavpNzszG3L3nnafN5c1fw3lFJFQZDUUNFdlHw3lFJEQZTEUNDfJX/V9EQldJ88DtFT2MbM+M9tvZsNmtqHO7R8zs31m9mMz+4GZzZtsIPGlGjSUU0RCFB8KmvaSEE2Tv5lNA+4GbgYWArea2cKa3f4H6HH33wUeAr7QagCq74uInKsT5wFa6fnfAAy7+9Pu/iKwBVgR38HdH3X3U9HVJ4A5rTx5tdSjH1oPmyYwNaa2CVftT0OuuX8g0fdBK8n/KuBw7PpItK2R1cD36t1gZmvNbNDMBn/y0yPc9Z19WpVT9POEE1DbSHxV0Lse3sv0y6/qTuJxEz3ha2bvB3qAN9e73d03AZsALpzd7fE1etTjD5cmMDWmtpFqCag6EvKCmRddmsTjNp3kZWa/B3zG3W+Krn8SwN0/V7NfL/D3wJvd/ZlmT3zF1a/3f/v3/1TSFxFp0dChUW687o0nx54duazdx2ql7DMAdJvZfDObCawE+uM7mNmbgK8By1tJ/ADzZ12ixC8iMglL5nUx/tyRg0k8VtPk7+7jwDpgO/Ak8KC77zWzO81sebTbF4FXAt8ys11m1t/g4UREJAdaqvm7+1Zga822O2KXexOOS0REUpTLtX1ERCRdSv4iIgFS8hcRCZCSv+SCZrKeT20iaVLyl1zQTNbzqU0kTblZ0lnCppms51ObSJpy/TOOIiJyrlL/jKOIiKRLyV9EJEBK/iIiAVLyFxEJkJK/iEiAlPxFRAKk5C8iEiAlf8kVLWmgNpDOUPKXXNGSBmoD6Qwt7yC5oiUN1AbSGVreQUSkQLS8g4iITJmSv4hIgJT8RUQCpOQvIhIgJX8RkQAp+UsuhTjRKcTXLNlR8pdcCnGiU4ivWbKjSV6SSyFOdArxNUt2NMlLRKRANMlLRESmTMlfRCRASv4iIgFS8hcRCZCSv4hIgJT8RUQCpOQvuRXSjNeQXqvkQ0vJ38z6zGy/mQ2b2YY6ty8zs/82s3EzuyX5MCVEIc14Dem1Sj40neFrZtOAu4G3AyPAgJn1u/u+2G4/B24HPp5GkBKmkGa8hvRaJR9aWd7hBmDY3Z8GMLMtwArgbPJ3959Ft72cQowSqCXzuti8emnWYXRESK9V8qGVss9VwOHY9ZFo26SZ2VozGzSzwWPHjk3lIUREJAEdPeHr7pvcvcfde6688spOPrWIiMS0kvyPAHNj1+dE20REpKBaSf4DQLeZzTezmcBKoD/dsEREJE1Nk7+7jwPrgO3Ak8CD7r7XzO40s+UAZna9mY0A7wG+ZmZ70wxaRETa01LN3923uvsCd/9td/9stO0Od++PLg+4+xx3v8Tdr3D3N6QZtISlzBOgyvzaJN80w1dyr8wToMr82iTf9DOOkntlngBV5tcm+aafcRQRKRD9jKOIiEyZkr+ISICU/EVEAqTkLyISICV/KYQyjocv42uS4lDyl0Io43j4Mr4mKQ6N85dCKON4+DK+JikOjfMXESkQjfMXEZEpU/IXEQmQkr+ISICU/EVEAqTkL4VShrHxZXgNUnxK/lIoZRgbX4bXIMWncf5SKGUYG1+G1yDFp3H+IiIFonH+IiIyZUr+IiIBUvKXwinqaJmixi3lpOQvhVPU0TJFjVvKSaN9pHCKOlqmqHFLOWm0j4hIgWi0j4iITJmSvxRWkU6gFilWCYOSvxRWkU6gFilWCYNO+EphFekEapFilTDohK+ISIHohK9IJM/19DzHJmFT8pfCy3M9Pc+xSdhU85fCy3M9Pc+xSdjU85fCWzKvi/W9C9i440CuyitDh0bZuOMA63sXsGReV9bhiJxDyV9KIY/llTzGJFLVUvI3sz4z229mw2a2oc7tF5rZA9HtO83s6sQjFZnA+t4FLOueRd+i2ZmfYK2e5O1bNJtl3bNU8pFcalrzN7NpwN3A24ERYMDM+t19X2y31cCou/+Oma0EPg+8L42ARepZMq+LzauXctu9O3n84HF2HznBPauu73i5ZejQKGvuH2D01BgAm1cv7ejzi7SqlRO+NwDD7v40gJltAVYA8eS/AvhMdPkh4CtmZp7VJAIJ1vreBew+coLRU2Pc9fBeLr1oBn2LZvPgwM/BjPf2zD3n8rY9Rye8fTL7bttzlJOnxxk9NUbXxTPU45dcazrJy8xuAfrcfU10/QPAUndfF9tnT7TPSHT9qWif4zWPtRZYG11dBOxJ6oWkaBZwvOle2VOcEZt50SXTXnn5a80umGYzLrzEX3553C64YDpAvcu1214+/X/Tp118WUv7nnf72Jlfu7/80ku/eu4X/uILv07zdaJjnrSixPk6d39Vuw/S0aGe7r4J2ARgZoNJzFJLm+JMVhHiNLPB8RPP5DpGKEZbguJMmpklsjRCKyd8jwBzY9fnRNvq7mNm04HLgGeTCFBERJLXSvIfALrNbL6ZzQRWAv01+/QDq6LLtwCPqN4vIpJfTcs+7j5uZuuA7cA04D5332tmdwKD7t4P3At8w8yGgeeo/IFoZlMbcXeS4kxWEeIsQoygOJMWVJyZreopIiLZ0QxfEZEAKfmLiAQo1eRvZu8xs71m9rKZNRxC1Wj5iOgk885o+wPRCec04rzczL5vZgej/8+bFmpmbzWzXbF/p83sXdFtXzezn8ZuW5xVnNF+L8Vi6Y9tT709W2zLxWb2w+i98WMze1/stlTbsp2lSszsk9H2/WZ2U5JxTSHOj5nZvqj9fmBm82K31T3+GcV5u5kdi8WzJnbbquh9ctDMVtXet8NxfjkW4wEzez52W0fa08zuM7NnrDJvqt7tZmZ/F72GH5vZdbHbJt+W7p7aP+D1wOuAx4CeBvtMA54CrgFmAj8CFka3PQisjC5/FfhwSnF+AdgQXd4AfL7J/pdTObF9cXT968AtabblZOIEftVge+rt2UqMwAKgO7r8WuAo8Oq023Ki91psn78EvhpdXgk8EF1eGO1/ITA/epxpGcb51tj778PVOCc6/hnFeTvwlTr3vRx4Ovq/K7rclVWcNfv/FZWBLZ1uz2XAdcCeBrf/EfA9wIAbgZ3ttGWqPX93f9Ld9zfZ7ezyEe7+IrAFWGFmBryNynIRAPcD70op1BXR47f6PLcA33P3UynF08hk4zyrg+3ZNEZ3P+DuB6PLvwCeAa5MIZZadd9rNfvE438I+MOo7VYAW9z9jLv/FBiOHi+TON390dj77wkq8286rZX2bOQm4Pvu/py7jwLfB/pyEuetwDdTiqUhd3+cSqeykRXAZq94Ani1mc1mim2Zh5r/VcDh2PWRaNsVwPPuPl6zPQ2vcfej0eX/BV7TZP+VnP/m+Gz0VezLZnZh4hFWtBrnK8xs0MyeqJam6Fx7TqotzewGKr2xp2Kb02rLRu+1uvtEbXWCStu1ct9Oxhm3mkqPsKre8U9Dq3G+OzqeD5lZdcJoLtszKp/NBx6Jbe5UezbT6HVMqS3bXt7BzHYAv1Xnpk+5+7+2+/hJmSjO+BV3dzNrOP41+kv7RirzHqo+SSXRzaQyBvdvgDszjHOeux8xs2uAR8xsN5UkloiE2/IbwCp3fznanFhbhsDM3g/0AG+ObT7v+Lv7U/UfIXUPA9909zNm9kEq36rellEsrVgJPOTuL8W25ak9E9N28nf33jYfotHyEc9S+VozPeqB1VtWomUTxWlmvzSz2e5+NEpIz0zwUO8Fvu3uY7HHrvZ0z5jZPwIfzzJOdz8S/f+0mT0GvAn4FxJqzyRiNLNLge9S6SQ8EXvsxNqyjsksVTJi5y5V0sp9OxknZtZL5Q/um939THV7g+OfRrJqGqe7x5d5uYfKOaHqfd9Sc9/HEo/wN8/V6rFbCXwkvqGD7dlMo9cxpbbMQ9mn7vIRXjmT8SiV+jpUlo9I65tEfHmKZs9zXj0wSnLVuvq7SG+10qZxmllXtVRiZrOAPwD2dbA9W4lxJvBtKvXLh2puS7Mt21mqpB9YaZXRQPOBbuC/EoxtUnGa2ZuArwHL3f2Z2Pa6xz/DOGfHri4HnowubwfeEcXbBbyDc79NdzTOKNZrqZww/WFsWyfbs5l+4LZo1M+NwImoszS1tkz57PWfUqk/nQF+CWyPtr8W2FpzFvsAlb+mn4ptv4bKB2wY+BZwYUpxXgH8ADgI7AAuj7b3APfE9ruayl/ZC2ru/wiwm0qi+ifglVnFCfx+FMuPov9Xd7I9W4zx/cAYsCv2b3En2rLee41KWWl5dPkVUdsMR211Tey+n4rutx+4OeXPTrM4d0SfqWr79Tc7/hnF+TlgbxTPo8C1sfv+RdTOw8CfZxlndP0zwN/W3K9j7UmlU3k0+myMUDmX8yHgQ9HtRuWHtZ6KYumJ3XfSbanlHUREApSHso+IiHSYkr+ISICU/EVEAqTkLyISICV/EZEAKfmLiARIyV9EJED/D5eMMT8obebJAAAAAElFTkSuQmCC\n", "text/plain": [ "