Skip to content

Instantly share code, notes, and snippets.

@ceptreee
Last active June 15, 2019 20:31
Show Gist options
  • Save ceptreee/d4e21358dbd7bd3d0e8640b4792097dd to your computer and use it in GitHub Desktop.
Save ceptreee/d4e21358dbd7bd3d0e8640b4792097dd to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np \n",
"import matplotlib.pyplot as plt \n",
"import sympy as sb \n",
"from IPython.display import Math\n",
"\n",
"x = sb.Symbol('x')\n",
"sb.init_printing()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"# 関数\n",
"f = sb.exp(-x)\n",
"ff = sb.lambdify(x,f)\n",
"\n",
"# 分子の多項式の次数\n",
"n = 2\n",
"# 分母の多項式の次数\n",
"m = 4\n",
"\n",
"# x=0でのn+m+1次までのテイラー展開\n",
"# .removeO()はランダウの記号を取り除く\n",
"f_Taylor = sb.series(f,x,0,n+m+1).removeO()\n",
"ff_Taylor= sb.lambdify(x,f_Taylor)\n",
"\n",
"# テイラー展開の係数\n",
"a = f_Taylor.as_poly().all_coeffs()[::-1]\n",
"a = np.array(a).astype(np.float64)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# パデ近似の分母と分子の係数を求める\n",
"def Pade(a,n,m):\n",
" I1 = np.eye(n+1)\n",
" O1 = np.zeros([m,n+1])\n",
" A1 = np.zeros([n+m+1,m])\n",
"\n",
" for j in range(n+1):\n",
" for i in range(1,m+1):\n",
" if j - i >= 0:\n",
" A1[j,i-1] = -a[j-i]\n",
"\n",
" for j in range(n+1,n+m+1):\n",
" for i in range(1,m+1):\n",
" if j - i >= 0:\n",
" A1[j,i-1] = -a[j-i]\n",
"\n",
" A = np.hstack([np.vstack([I1,O1]),A1])\n",
" c = np.linalg.solve(A,a)\n",
"\n",
" p = c[:n+1]\n",
" q = c[n+1:]\n",
" q = np.hstack([1,q])\n",
" return p,q\n",
"\n",
"p,q = Pade(a,n,m)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# 有理関数\n",
"def RationalFunciton(x,p,q):\n",
" n = len(p)-1\n",
" m = len(q)-1\n",
" P = np.sum([p[i]*x**i for i in range(n+1)],axis=0)\n",
" Q = np.sum([q[i]*x**i for i in range(m+1)],axis=0)\n",
" R = P / Q\n",
" return R,P,Q"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"R,P,Q = RationalFunciton(x,p,q)\n",
"zeros = np.array(sb.solve(P)).astype(complex)\n",
"poles = np.array(sb.solve(Q)).astype(complex)"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1153880b8>"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQEAAAD8CAYAAAB3lxGOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAH4xJREFUeJztnX9MXNedt59jwJhmwAZSg2B4y7ivTQvYxJNJMQqFJN1surtVN61qtV2ldZpWUaRdKa22iTabv9KmahtLSf7ISlHVvHHURF29UbLZV+2u3CYNEFKCjElIjBNwCk49uMZZwDU2eGzG5/1jPGNjj6/NnYH7bc73kUbDXO6ceeByP5x77vlhrLUoiuIuq4IWUBQlWDQEFMVxNAQUxXE0BBTFcTQEFMVxNAQUxXE0BBTFcTQEFMVxNAQUxXEKg/jQa6+91tbX16/oZ54+fZrVq1ev6GcuBfXzj2Q3CMZv7969/2Ot/fhV7WytXfHH9ddfb1eaV199dcU/cymon38ku1kbjB8wYK/yfHTmcqClpSVoBU/Uzz+S3UC+nzMhMDs7G7SCJ+rnH8luIN/PmRAYGxsLWsET9fOPZDeQ7xdIw6CirBRnzpwhHo9z6tSpwBzWrl3Lu+++uyxlr1mzhnA4TFFRke8ynAmBlb4bsVTUzz9ebvF4nNLSUurr6zHGrJzUBSQSCYqLi/NerrWWqakp4vE4kUjEdznOXA5UVFQEreCJ+vnHy+3UqVNUVlYGFgAAhYXL87/WGENlZWXOtRxnQmBwcDBoBU/Uzz9XcgsyAADm5uaWrex8/GzOhICiKNlxJgTKy8uDVvBE/fwj2Q2goKAgaAVPnGkYlN5hQ/38k0+3l96cYOfuEQ4fm6dmXQn33dbA7VtrcyrzYx/72JLfk0wmVyw8nKkJdHd3B63gifr5J19uL705wQMvvsPEsXksMHFsngdefIeX3pzwXeaTTz7Jli1buO6664hEItx888385je/oa2tjWg0yvbt2zlx4gSQusvxgx/8gPb2dp5//nneeusttm3bxpYtW/jSl77EzMxMXn7Oi3EmBKzwqdXVzz/5ctu5e4T5M8lF2+bPJNm5e8R3mffccw+vv/46e/bsIRwOc9ddd/Hwww/z8ssvMzg4SCwW49FHH83sv2bNGnp7e/na177GN7/5TX7605/y9ttvs3nzZh566CHfHl44czkQdAvxlVA//+TL7fCx+SVtXwr33nsvt9xyC+Xl5ezfv58bb7wRSI0wbGtry+z31a9+FYA///nPHDt2jM7OTgB27NjB9u3bc/bIhjMhkP5lSkX9/JMvt5p1JUxkOeFr1pXkVO4LL7zABx98wBNPPMGvf/1rbr31Vn75y19m3feaa67J6bP84MzlwNDQUNAKnqiff/Lldt9tDZQULW6MKykq4L7bGnyXuXfvXh555BGeffZZVq1axbZt23j99dd5//33gVQfgtHR0Uvet3btWsrLy3nttdcA+MUvfrFsQexMTWC5GlXyhfr5J19u6bsA+bw78MQTTzAzM8PNN98MQCwWY9euXXz9618nkUgA8PDDD7Np06ZL3vvMM89wzz33MDc3x4YNG3j66ad9e3jhTAgoytVw+9banG8JXsjTTz/N7OwspaWli7bv2bPnkn0PHjy46PV1113HG2+8kTeXy+HM5UA0Gg1awRP1849kN/DXT2AlcSYEpqeng1bwRP38I9kNYGFhIWgFT5wJgYurWtJQP/9IdoPUbUDJOBMCiqJkx5kQ2LBhQ9AK0Ps4jPcs3jbeA72Py/DzQLKfZDdA9HTo4FAIXNw6Gwi1UXj+zvNBMN6Tel0bleHngWQ/yW4gfxShMyEgorNLpAO270qd+L/7Uep5+y6IdMjw80Cyn2Q3gPn5pXU7vummmxgYGFgmm0vJWwgYYwqMMW8aY36VrzI/kkQ6IPZt6Hkk9RzpCNpISeNxufZRJp81gXuB5ZlSNQ9UVlYGrZBivAcGnoKO+1PP5/7oxPhdBsl+eXPzuFzLhXg8zqc+9Sl27NjBli1b+MpXvsLc3ByvvPIKW7duZfPmzdx1112ZHoQXcrlhx3nlapcq8noAYeAV4BbgV1faP4hlyJLJ5Ip/5iWMdVv700jq+aLXIvw8kOzn5bZ///6lFZY+Jq88vPhY5cDY2JgFbG9vr7XW2m9961v2hz/8oQ2Hw3ZkZMRaa+03vvEN+9hjj1lrre3s7LR79uyxH374of3sZz9rT5w4Ya219ic/+Yl96KGHLik/289IAMuQPQ7cD5zNU3l5p6en58o7LTcTg5k2AOB8G8HEoAw/DyT75dVtGS7XTp48SV1dXWb48B133MErr7xCJBLJjBnYsWPHJT/HG2+8kRl2fN111/HMM8/wwQcf5OxzMTmPHTDGfAE4aq3da4y5yWO/u4G7AWpqaujq6gJSt3dKS0szjTuVlZU0NTVlfiGFhYW0t7czODjI8ePHgdQgjMnJSQ4dOgTAxo0bKS4uZt++fQCsX7+eTZs20dvbC5CZ831gYCBTnWptbSUejzMxkZo1pqGhgYKCAvbv3w9AdXU1kUiEvr4+AEpKSmhtbaW/vz/T0NPW1sb4+DhHjhwBoLGxkWQyychIahKK2tpawuEw/f39AIRC7cQiMfr6+jJVv/b2dkbnr+XE2BhdXV00NzeTSCQ4cOAAAHV1dVRVVWUaisrKyohGo/T29mZ6onV0dDA8PMzU1BSQmm5rdnY2s/JNfX09FRUVmVl5y8vLaWlpobu7G2stxhg6OzsZGhrKDMaJRqNMT09nOuKcPn2amZmZZT9ObW1tSz5OZ8+ezfw9XXyc1q5dy9mzZ0kkEpnf15o1a7DWZo5BUVERq1ev5uTJkxT88XVK9vwc03E/ds/Pma+6nuT/upFQKMSpU6c8yygqKsrMLLxq1SquueYaZmdnOXs29b/RWsv8/Dxzc3Mkk8nM+0+fPs3c3Bxnz54lmUySTCY5efIkp06d4tZbb+VnP/sZaUpLSzPvh1SX5IWFhUvOpyVxtVWGyz2AHwNx4CBwBJgDnvV6TxCXA6+99tqKf+ZSUD//eLkt6XLA43ItF9555x0L2N///vfWWmu/853v2IcfftjW1dXZAwcOWGut3bFjh3388cettecvB44ePbpon5MnT2YuHy4k8MsBa+0D1tqwtbYe+BrwO2vtHbmWm2/a29uDVvBE/fyTNzePy7VcCIVCfPrTn+aZZ55hy5YtTE9P873vfY+nn36a7du3s3nzZlatWsU999yz6H0f//jHM8OOt2zZwrZt23jvvfdycsnK1abF1TyAmxDaMLh3794V/8yloH7+8XJbcsPgMjA8PGybmpqWrfxcawJ5nU/AWtsFdOWzzHyRvk6Vivr5R7IbkGkTkIozPQYVJSg+8YlPZBpDJeJMCMRisaAVPFE//1zJzQY8XfpyTiqSj5/NmRCYnJwMWsET9fOPl9uaNWuYmpoKNAjOnDmzLOXac0uTr1mzJqdynJlj8NChQ3zyk58MWuOyqJ9/vNzC4TDxeJwPP/xwha3Oc+rUqZxP1MuxZs0awuFwTmU4EwKKmxQVFRGJRAJ16OrqYuvWrYE6eOHM5cDGjRuDVvBE/fwj2Q3k+zkTAumuw1JRP/9IdgP5fs6EgORbNKB+uSDZDeT7ORMCiqJkx5kQWL9+fdAKnqiffyS7gXw/E8T901gsZldyDjVILQBRWCj3Zoj6+UeyGwTjZ4zZa629qh5eztQE0mPWpaJ+/pHsBvL9nAkBRVGy40wISL9No37+kewG8v2caRNQFJfQNoEsSA8d9fOPZDeQ7+dMCCzLfO15RP38I9kN5Ps5EwKKomTHmTaB+fl5SkpKVvQzl4L6+UeyGwTjp20CWYjH40EreKJ+/pHsBvL9nAmB9OIVUlE//0h2A/l+zoSAoijZcSYEGhoaglbwRP38I9kN5Ps5EwIFBQVBK3iifv6R7Aby/ZwJgfQCllJRP/9IdgP5fs6EgKIo2XEmBKqrq4NW8ET9/CPZDeT7yZ2JIc8EPe30lVA//0h2g/z5vfTmBDt3j3D42Dw160q477YGbt9am3O5ztQE+vr6glbwRP38I9kN8uP30psTPPDiO0wcm8cCE8fmeeDFd3jpzdz7IDgTAoryl8zO3SPMn0ku2jZ/JsnO3SM5l+1MCEjuWw7qlwuS3SA/foePzS9p+1LIOQSMMXXGmFeNMe8aY4aNMffmbLUMtLa2Bq3gifr5R7Ib5MevZl32ILnc9qWQj5rAAvDP1tpPA9uAfzTGNOah3LzS398ftIIn6ucfyW6QH7/7bmugpGhxp6OSogLuuy333og53x2w1v4J+NO5r2eNMe8CtYCoHhLz87lXm5YT9fOPZDfIj1/6LsBy3B3I63wCxph6oAdottYev+h7dwN3A9TU1Fz/3HPPAbBhwwZKS0sZGhoCoLKykqamJnp6egAoLCykvb2dwcFBjh9PFRmLxZicnOTQoUNAasHH4uLizHJP69evZ9OmTZmpnouLi0kkEoRCocwsL62trcTj8cwIr4aGBgoKCjK9u6qrq4lEIpmW3ZKSElpbW+nv788c1La2NsbHxzly5AgAjY2NJJNJRkZSjTW1tbWEw+HMf4JQKEQsFqOvr49EIgFAe3s7o6OjjI2NEQqFaG5uJpFIcODAAQDq6uqoqqrKTFFVVlZGNBqlt7eXhYUFADo6OhgeHmZqagqAlpYWZmdnGRsbA6C+vp6KigoGBwcBKC8vp6Wlhe7ubqy1GGPo7OxkaGiImZkZAKLRKNPT0xw8eBCA06dPc8MNNyz7cWpra2NgYGBJx+nw4cOsWrVqRY7T0aNHAZZ0nI4dO0YoFFqR45Q+nyoqKq56PoG8hYAxJgR0Az+y1r7otW8Qk4okEgnRs76qn38ku0Ewfis+qYgxpgh4AXjuSgEQFOPj40EreKJ+/pHsBvL98nF3wABPAe9aax/NXWl5SFcFpaJ+/pHsBvL98lETuBH4BnCLMeatc4+/zUO5iqKsAPm4O9ALmDy4LCuNjeLuWi5C/fwj2Q3k+znTYzCZTF55pwBRP/9IdgP5fs6EQPp2kFTUzz+S3UC+nzMhoChKdpwJgdra3HtWLSfq5x/JbiDfz5kQCIfDqS96H4fxnsXfHO9JbQ+QjJ9QJPtJdgP5fs6EQGYQR20Unr/zfBCM96Re10aDUgPcGASzXEh2A/l+zkwvliHSAdt3pU782Ldh4KnU60hHwGKKEgzO1ARCodD5F5GOVAD0PJJ6FhAAi/wEItlPshvI93NmVeJFpC8BtCagfETRVYmzkJnsMR0A23fBLQ+evzS4uLFwhXFhsszlQrIbyPdzJgTS48KZGFz8nz/dRjAxGJQacIGfUCT7SXYD+X7uNQy2f/fSbZEOvRxQnMWZNoGFhQUKC+Vmnvr5R7IbBOOnbQJZGB0dDVrBE/Xzj2Q3kO/nTAik54aTivr5R7IbyPdzJgQURcmOMyHQ3NwctIIn6ucfyW4g38+ZEJB+m0b9/CPZDeT7ORMC6fnhpaJ+/pHsBvL9nAkBRVGy40wI1NXVBa3gifr5R7IbyPdzJgSqqqqCVvBE/fwj2Q3k+zkTAoGOWrwK1M8/kt1Avp8zIaAoSnacCYGysrKgFTxRP/9IdgP5fs4MIFIUl9ABRFno7e0NWsET9fOPZDeQ7+dMCCwsLASt4In6+UeyG8j3cyYEFEXJjjNtAmfPnmXVKrmZp37+kewGwfiteJuAMebzxpgRY8z7xph/yUeZ+WZ4eDhoBU/Uzz+S3UC+X84hYIwpAP4N+BugEfi6MUbcguxTU1NBK3iifv6R7Aby/fJRE/gM8L61dsxaexr4d+Dv81CuoigrQD5CoBY4dMHr+LltomhpaQlawRP1849kN5Dvl48pUE2WbZe0Nhpj7gbuBqipqaGrqwuADRs2UFpaytDQEACVlZU0NTXR05NaDKSwsJD29nYGBwc5fvw4ALFYjMnJSQ4dSmXPxo0bKS4uZt++fQCsX7+eTZs2Ze7PFhcXU1tbyx/+8AdOnDgBQGtrK/F4nImJCQAaGhooKChg//79AFRXVxOJRDILR5SUlNDa2kp/fz/z8/MAtLW1MT4+zpEjRwBobGwkmUwyMjICpJakDofDmQUpQ6EQsViMvr6+zEQT7e3tjI6OEo/HWb16Nc3NzSQSicwY9Lq6OqqqqjL9z8vKyohGo/T29mZuPXV0dDA8PJypdra0tDA7O8vY2BgA9fX1VFRUMDiYWluhvLyclpYWuru7sdZijKGzs5OhoSFmZmYAiEajTE9Pc/DgQQDWrl0LsOzHqa2tjYGBgSUdp6KioozXch+n9HyBSzlOc3NzrF69ekWOU/p8WhLW2pweQBuw+4LXDwAPeL3n+uuvtyvNq6++uuKfuRTUzz+S3awNxg8YsFd5DufjcmAPsNEYEzHGrAa+Bvy/PJSrKMoKkPPlgLV2wRjzT8BuoAD4P9ZacfdE6uvrg1bwRP38I9kN5PvlZVkUa+1/Af+Vj7KWi4qKiqAVPFE//0h2A/l+crtZ5Zl0Y4tU1M8/kt1Avp8zIaAoSnacCYHy8vKgFTxRP/9IdgP5fs4MIFIUl9BJRbLQ3d0dtIIn6ucfyW4g38+ZEAiixrMU1M8/kt1Avp8zIWBMtt7NclA//0h2g78AP20TUJSPHtomkIX0ABOpqJ9/JLuBfD9nQiA96koq6ucfyW4g38+ZEFAUJTvOhEA0Gg1awRP1849kN5Dv50wITE9PB63gifr5R7IbyPdzJgTSM69IRf38I9kN5Ps5EwKKomTHmRDYsGFD0AqeqJ9/JLuBfD9nQmDJky+uMOrnH8luIN/PmRCQ3mFD/fwj2Q3k+zkTAoqiZMeZEKisrEx90fs4jPcs/uZ4T2p7gGT8hCLZT7IbyPdzJgSamppSX9RG4fk7zwfBeE/qdW2wHToyfkKR7CfZDeT7ORMC6ZVyiHTA9l2pE/93P0o9b9+V2h4gGT+hSPaT7Aby/ZwJgUVEOiD2beh5JPUccAAoSpA4EwKFhRcssTDeAwNPQcf9qeeL2wgCYJGfQCT7SXYD+X7uTSqSbgNIXwJc/FpRPgLopCJZyCwAMTG4+IRPtxFMBLtAhPQFKiT7SXYD+X6y6yl5JL1cNu3fvfSbkY7AawEZP6FI9pPsBvL9nKkJKIqSHWdCIBa7qsujwFA//0h2A/l+zoTA5ORk0AqeqJ9/JLuBfD9nQuDQoUNBK3iifv6R7Aby/XIKAWPMTmPMe8aYt40x/2GMWZcvMUVRVoZcawK/BZqttVuAUeCB3JWWh40bNwat4In6+UeyG8j3yykErLW/sdYunHv5BhDOXWl5KC4uDlrBE/Xzj2Q3kO+XzzaBu4D/zmN5eWXfvn1BK3iifv6R7Aby/a7YbdgY8zJQneVbD1pr//PcPg8CMeDL9jIFGmPuBu4GqKmpuf65554DUvOvlZaWZmZfqayspKmpKTPyqrCwkPb2dgYHBzOdLmKxGJOTk5kGl40bN1JcXJz5Za9fv55NmzbR29sLpJI4kUgQCoU4ceIEAK2trcTjcSYmJgBoaGigoKCA/fv3A1BdXU0kEqGvrw+AkpISWltb6e/vZ35+HoC2tjbGx8c5cuQIAI2NjSSTSUZGRgCora0lHA7T398PQCgUIhaL0dfXRyKRAKC9vZ3R0VHGxsYIhUI0NzeTSCQ4cOAAAHV1dVRVVZHuZl1WVkY0GqW3t5eFhVQlrKOjg+HhYaampgBoaWlhdnaWsbExAOrr66moqMj0XCsvL6elpYXu7m6stRhj6OzsZGhoKLNaTjQaZXp6OjNT7unTp7nhhhuW/Ti1tbUxMDCwpON0+PBhVq1atSLH6ejRowBLOk7Hjh0jFAqtyHFKn08VFRVX3W0457EDxpgdwD3A56y1c1fzniDGDuzfv5/GxsYV/cyloH7+kewG+fN76c0Jdu4e4fCxeWrWlXDfbQ3cvrU2675LGTuQUwgYYz4PPAp0Wms/vNr3BRECCwsLokdzqZ9/JLtBfvxeenOCB158h/kzycy2kqICfvzlzVmDYCUHED0BlAK/Nca8ZYx5Msfylo10lVMq6ucfyW6QH7+du0cWBQDA/JkkO3eP5Fx2TvFkrf3fORsoinJFDh+bX9L2peBMj0Hpt2nUzz+S3SA/fjXrSpa0fSm4N6mIovwFIrlN4C8G6aGjfv6R7Ab58bt9ay0//vJmateVYIDadSWXDYClIrdJNc+k7ztLRf38I9kN8ud3+9bavJz0F+NMTUBRlOw40yYwPz9PSUnujSjLhfr5R7IbBOOnbQJZiMfjQSt4on7+kewG8v2cCYF033OpqJ9/JLuBfD9nQkBRlOw4EwINDQ1BK3iifv6R7Aby/ZwJgYKCgqAVPFE//0h2A/l+zoRAevy5VNTPP5LdQL6fMyGgKEp2nAmB6upskyPJQf38I9kN5Ps5EwKRSCRoBU/Uzz+S3UC+nzMhkJ4rUCrq5x/JbiDfz5kQUBQlO86EgOS+5aB+uSDZDeT7OTOASFFcQgcQZSE9p7xU1M8/kt1Avp8zIZBeiEIq6ucfyW4g38+ZEFAUJTvOtAkkEgnRs9Kqn38ku0EwftomkIXx8fGgFTxRP/9IdgP5fs6EQHoxSqmon38ku4F8P2dCQFGU7DgTApJXrQX1ywXJbiDfz5kQSCaTV94pQNTPP5LdQL6fMyEwMpL76q3Lifr5R7IbyPdzJgQURcmOMyFQW5v/5Zvyifr5R7IbyPfLSwgYY75vjLHGmGvzUd5yEA6Hg1bwRP38I9kN5PvlHALGmDrgVuCPuessH9IHcaiffyS7gXy/fNQEHgPuB1a+/7GiKDmT09LkxpgvAhPW2iFjzJX2vRu4G6Cmpoauri4ANmzYQGlpKUNDQwBUVlbS1NRET09PSrCwkPb2dgYHBzl+/DgAsViMyclJDh06BMDGjRspLi5m3759AKxfv55NmzbR29sLQHFxMaFQiIGBgcwy0a2trcTj8cwSUQ0NDRQUFGSmh66uriYSiWSmhiopKaG1tZX+/v7MqLC2tjbGx8czPcIaGxtJJpOZ1uDa2lrC4XDmP0EoFCIWi9HX10cikQCgvb2d0dFR5ubm6Orqorm5mUQiwYEDBwCoq6ujqqoqs8Z9WVkZ0WiU3t5eFhYWAOjo6GB4eJipqSkAWlpamJ2dZWxsDID6+noqKioYHBwEoLy8nJaWFrq7u7HWYoyhs7OToaEhZmZmAIhGo0xPT3Pw4MH08WNmZmbZj1NbW9uSj1NJSUnm72m5j9PRo0cBlnSc0sd2JY5T+nxaClccQGSMeRnINl3qg8C/An9trf2zMeYgELPW/s+VPtTZSUV6H4faKEQ6zm8b74GJQWj/bnBeykeOvA4gstb+lbW2+eIHMAZEgKFzARAGBo0xIudXFjHZY20Unr8zdeJD6vn5O6E2KsPPA8l+kt1Avp/vywFr7TvA+vTrpdQEgiBdrQuUSAds35U68WPfhoGnUq8jHSQ+6ArW7QqI+P1dBsluIN/PmX4CYoh0pAKg55HU84WXBooSAM5MKrKwsEBhYU7toPkhfQlwUU1AjN9lkOwn2Q2C8dNJRbIwOjoatML5ANi+C2558PylwXiPDD8PJPtJdgP5fs6EQPrWTqBMDGb+8wPn2wgmBmX4eSDZT7IbyPeTW4f6KJLtNmCkI/U4d59bUVYaZ2oCzc3NQSt4on7+kewG8v2cCQHpt2nUzz+S3UC+nzMhkO7eKRX1849kN5Dv50wIKIqSHWdCoK6uLmgFT9TPP5LdQL6fMyFQVVUVtIIn6ucfyW4g38+ZEJA+alH9/CPZDeT7ORMCiqJkx5kQKCsrC1rBE/Xzj2Q3kO/nzAAiRXEJHUCUhfQUVlJRP/9IdgP5fs6EQHo+Pqmon38ku4F8P2dCQFGU7DjTJnD27FlWrZKbeernH8luEIyftglkYXh4OGgFT9TPP5LdQL6fMyGQnutdKurnH8luIN/PmRBQFCU7zoRAS0tL0AqeqJ9/JLuBfD9nQmB2djZoBU/Uzz+S3UC+nzMhkF7vTSrq5x/JbiDfz5kQUBQlO4H0EzDGfAh8sMIfey0gcom0c6iffyS7QTB+n7DWfvxqdgwkBILAGDNwtZ0ngkD9/CPZDeT76eWAojiOhoCiOI5LIfCzoAWugPr5R7IbCPdzpk1AUZTsuFQTUBQlC06GgDHm+8YYa4y5NmiXCzHG7DTGvGeMedsY8x/GmHUCnD5vjBkxxrxvjPmXoH0uxBhTZ4x51RjzrjFm2Bhzb9BOF2OMKTDGvGmM+VXQLpfDuRAwxtQBtwJ/DNolC78Fmq21W4BR4IEgZYwxBcC/AX8DNAJfN8Y0Bul0EQvAP1trPw1sA/5RmB/AvcC7QUt44VwIAI8B9wPiGkOstb+x1qbnonoDCAfpA3wGeN9aO2atPQ38O/D3ATtlsNb+yVo7eO7rWVInW22wVucxxoSBvwN+HrSLF06FgDHmi8CEtXYoaJer4C7gvwN2qAUOXfA6jqCT7EKMMfXAVqA/WJNFPE7qH87ZoEW8KAxaIN8YY14GqrN860HgX4G/XlmjxXj5WWv/89w+D5Kq6j63km5ZMFm2iatBGWNCwAvAd621x4P2ATDGfAE4aq3da4y5KWgfLz5yIWCt/ats240xm4EIMGSMgVRVe9AY8xlr7ZGg/dIYY3YAXwA+Z4O/fxsHLlxNMwwcDsglK8aYIlIB8Jy19sWgfS7gRuCLxpi/BdYAZcaYZ621dwTsdQnO9hMwxhwEYtZaMQNPjDGfBx4FOq21HwrwKSTVQPk5YALYA/yDtVbEpHkmlebPANPW2u8G7XM5ztUEvm+t/ULQLtlwqk3gL4AngFLgt8aYt4wxTwYpc66R8p+A3aQa3f6vlAA4x43AN4Bbzv2+3jr3n1dZAs7WBBRFSaE1AUVxHA0BRXEcDQFFcRwNAUVxHA0BRXEcDQFFcRwNAUVxHA0BRXGc/w8GLugbNQKJeQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.close('all')\n",
"plt.figure()\n",
"plt.plot(zeros.real,zeros.imag,'o',label='zero')\n",
"plt.plot(poles.real,poles.imag,'x',label='pole')\n",
"\n",
"plt.grid()\n",
"mx = np.max(np.abs(np.hstack([zeros,poles])))\n",
"\n",
"plt.xlim([-mx,mx])\n",
"plt.ylim([-mx,mx])\n",
"plt.gca().set_aspect(1/plt.gca().get_data_ratio())\n",
"plt.legend()\n"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"xx = np.linspace(0,5,100)\n",
"yy = ff(xx)\n",
"yy_Taylor = ff_Taylor(xx)\n",
"yy_Pade,_,_ = RationalFunciton(xx,p,q)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQoAAAD8CAYAAACPd+p5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnXd4VFX6gN8zk0kvhDQSQoegBAiESFEEVFQEV11BKVYEFV1d689eYF3X3lDZ1VXBXkBdESu4gOLSuwmhFxNCgBBIbzPn98dkxgDJ5E4yM/ceve/zzMNM5s6ZlzOTL+eee873CSklJiYmJp6w6C1gYmJifMxAYWJi0ixmoDAxMWkWM1CYmJg0ixkoTExMmsUMFCYmJs3ik0AhhHhLCHFQCPFLE88LIcRMIcQOIcQmIUSmL97XxMQkMPhqRDEHGOXh+QuAHvW3G4B/+uh9TUxMAoBPAoWU8kfgiIdDLgbekU5WAG2EEMm+eG8TExP/ExSg92kP/NrgcV79zwoaHiSEuAHniIOoUOuAqPAQwmKTCQ8Px2q1UlFRAUBQUBBhYWGUlpa6XkdkZCQVFRXY7XYAIiIiqK2tpaamBoDQ0FCEEFRWVgJgs9kICQmhrKwMAIvFQkRERIvbkFJitVqJiIigvLwch8MBQGRkJNXV1dTW1gIQFhaGlJKqqioAgoODsdlslJeXA2C1WgkPD/dJGy4vgKioKCorK6mrq0NKSUREBHa7nerqagBCQkI89vGJbQCEh4d71UZLPychhPv51n5ODdsoKStn77E6BJCeEu2Tz6miosL9uLWfU0v62Ns2tm3bdlhKmUBzSCl9cgM6A7808dxXwNAGj38ABnhqb0CyRe54oKucM2eOVIHFixfrraAZlVyl9J/vip2HZad7F8iLXlnmszZV61tgjdTw+x2oqx55QIcGj1OB/Z5eYMdKN9thli360q9iviIjI0NvBc2o5Ar+862uc5CWFEmv5Giftala32olUIFiPnB1/dWPwcAxKWWBpxeU1NkASHXscw8hjUzDoZ/RUckV/Oc7LC2B7+8YzhOX9vFZm6r1rVZ8dXn0Q2A50FMIkSeEmCKEmCaEmFZ/yNfALmAH8G/g5ubaLKkVAFyQWs5///tfX2j6lV27dumtoBmVXEEtX5VcvcEnk5lSyonNPC+Bv3jTpsMSTLUUZIXmc9uCzxkzZkyrHE1MGiKl5GBpNYlRIQgh9NYxPIZdmdmmTSw/laViERJbUbZ7JtmodO7cWW8FzajkCv7xPVRWzaB//MCIZ5f4tF3V+lYrhg0U0dHRvLslmBnberBgxU5WrFiht5JH2rZtq7eCZlRyBf/45hY45xISo0J82q5qfasVwwaKiooKbCkZTP9wLdt/PcRnn32mt5JH1q1bp7eCZlRyBf/45h4oAeCUdr674gHq9a1WDBsoAMaOHeu+/+mnn7oXpZiYtJbcA84RxSnJUTqbqEGgVmZ6TVBQEGeffTbdOqVw3ZnJhNtg/fr1ZGYacz9ZbGys3gqaUckV/OPrOvXw9YhCtb7VimFHFGFhYYSEhDDy7BE80G07N3TYy+efztVbq0lUWmijkiv43rfW7mDHQedy7p7tfDuiUK1vtWLYQOFauHL+n8axviqFcFHDkW0/62zVNEuXLtVbQTMquYLvfXcfLqfG7qBj23AiQ3w7qFatb7Vi2EDh4vzzz2d+nnN4OKJdOTk5OTobNY5K8ycquYLvfTu2DeeTG4cw/aJePm0X1OtbrRg+UISHh5Nv6QTAqLZ5fPLxhzobNY5Ki3ZUcgXf+4barAzs0pazT0nyabugXt9qxbCBIirqt3PHc/90GZurk4kSVRzMXmrIqD18+HC9FTSjkiuo5auSqzcYNlA03Ag2ZswYPtsXA8DIlDI2b96sl1aTbNy4UW8FzajkCr73/b+5G3n2u61U1NT5tF1Qr2+1YthA4Uq6AfWnH7aufHakG++sr+Tjjz/W0axxiouL9VbQjEqu4Fvf4vIa5q7N481luwkJsvqsXXf7ivWtVgwbKE5k9CXjGfvyer5YlstHH31kyNMPE+OTvd+5IrNXSjRWy+9zPsEfGDZQhIeHH/d41KhRREc7r37s2rWLtWvX6qHVJEZdCNYYKrmCb31/2X8MgN4pvl1o5UK1vtWKYQOFKx+ii9DQUC6++GLGDT+VT27J5JMP39PJrHGOHPGUW9hYqOQKvvX9Jd8ZKNLbx/iszYao1rdaMWygcCUCbcjEiRO5c6CFy+J2UPvrqpOCiZ7s2bNHbwXNqOQKvvV1nXqk+2lEoVrfasWwgaIxzj33XD7b67xsekmnMpYsWaKvkIlSlFbVsvtwOcFWCz0Szc1g3mDYQBEScnKegKCgIOztMqmWQZwZvo9PP3pHB7PG6dq1q94KmlHJFXznW1Fj55J+KZzbK4ngIP989VXrW60YNlBYrY1fuppwxTUsOtYRi5AkVOQaJvFuwwViRkclV/Cdb1J0KC9O6M+rV/hvwlG1vtWKYQOFq8DJiZx22ml8nRcJwLiOx/jyS2Ok81dpoY1KrqCWr0qu3uCrLNyjhBBb64sQ39fI8x2FEIuFEOvrixSPbsV7kdj7LI7JMPqEFPD9Fx+0Tt7kD8PynUUcLKnSW0NJWh0ohBBW4FWchYh7AROFECduy3sI+ERK2R+YAMxqrt2goKa3/1551TW8uiuVv2/vwTeLfuTAgQMt9vcVcXFxeitoRiVX8I1vZY2dK95YwZAn/0tVrf+ulqnWt1rxxYhiILBDSrlLSlkDfISzKHFDJOC6HhVDM1XCwJm4pim6devGt3vDefiDtew/WMz777/fMnMfkp6erreCZlRyBd/4bjlQgkNC94RIQm2+X7rtQrW+1YovsnY0VoB40AnHTAe+F0LcCkQAIxtrqGGR4sTERPflz65duxIVFeU+/4uLi+Pqq6/mp59+AmDWrFnceeedrF+/npIS53XyrKwsCgsL+fVXp1qPHj0ICQnhl19+wdV+Wloay5YtA5xXWYYMGcKaNWvcxWwHDRpEXl4e+fn5APTs2ROr1erOidGuXTu6dOnC8uXLKSsrIyEhgUGDBrFy5Ur3JOuQIUPYvXu3e9TTq1cv7HY7W7dudXZe+/akpqaycuVKwFnoNisri+XLl7vXkgwdOpRt27Zx8OBBAHr37k11dTXbt28HoEOHDiQlJbFmzRrAmcE8MzOTZcuWuffMDBs2jOzsbIqKiigrK+OMM86gtLTUXbCmc+fOtG3b1p0cNjY2loyMDJYude7WFUIwfPhwNm7c6N7PkJmZyZEjR9xrBxr7nNLT0/nxxx8B5yhx6NChrFu3zqvPqbCw0L19u6Wf05oDzuCQGFTJypUr/fY5LVq0iNDQUJ98TuDMmOXvz0kTWgqUeroBlwFvNHh8FfDyCcfcCdxVf38IkANYPLWblpbmsbhqSUmJ7JCSIGdMHCBnXjdArl692tv6rD5FpeK0KrlK6Rvfv364Tna6d4F8f8Xe1gt5QLW+JYBFirUUIJ4CfFIfmJYDoUC8p0abSwASFRXFmAtG8WCPHdyQupv3337TW2+f4mlOxWio5Aq+8d3w61EAMjr4Z+m2C9X6Viu+CBSrgR5CiC5CiGCck5XzTzhmH3AOgBDiVJyB4pCnRiMjI5t948uvmMzS8k6EiDrCDq/TdU3F0KFDdXtvb1HJFVrvW1xew96iCkJtFnom+Xedg2p9q5VWBwopZR1wC/AdsAXn1Y1sIcTfhBAX1R92F3C9EGIj8CFwbf2wp0maWkfRkOHDh/PpbmdAubJLMfPmzmv5f6SVqFT4RSVXaL3v1kJnouY+7WMIsvp36ZBqfasVXxUp/hpnxfKGP3ukwf0c4Axv2tSy4ctisZDS/wIO1+yhV3Ahzy14n6uuvsqbt/EZrsk5FVDJFVrvO7hrHJumn8eRshofGTWNan2rFcOuzNTKdVOm8EF+CgBnJRwmNzdXZyMTIxIdaqNzfITeGspi2EAREaHtQ01OTmZTdUcALo3fy1tvvOZPrSbJysrS5X1bgkquoJavSq7eYNhAUVtbq/nYSydN4bk9PbhqSQLvvPt+o7ks/E1hYWHA37OlqOQKrfPdV1TB4H/8wH2fbvKhUdOo1rdaMWygqKnRfj55/vnn89J/D/LZj1soPHiIefMCP6npWjCkAiq5Qut81/9azIGSKg6VBuaPh2p9qxXDBgpvsFqt3Hjjje7Hr856VUcbEyOxdq9zdWL/jm10NlEbwwYK1zJYrUydOpVe3VL4/K/9ubV/NRs2bPCTWeP06NEjoO/XGlRyhdb5rtnjDBRZndv6SscjqvWtVgwbKLwtzZaUlMTQM4fzpza7GBu3mzdfe8VPZo3TWEYuo6KSK7Tct7SqltwDJQRZBBmpgRlRqNa3WjFsoGjJKsurp/6F70s7EyzsJJdv4ujRo34waxzXJiYVUMkVWu67ft9RHBJ6t48hLNh/O0YbolrfasWwgaIlnH766Xy61/mXY3KnA7z5xr91NjLRkzX18xNZnWJ1NlEfw+5gsdlsXr9GCMGg8yexfW8ePWyH2LPiP9TV3RGQjTqJiYl+fw9foZIrtNz3gt7tCLNZGdglcIFCtb7Vimhmy4VuDBgwQLakGlhlZSXTJ5/NU6fksLaqPXtOe4yxY8f6wfB46urqlNk5qJIrqOWrkiuAEGKtlLLZVWKGPfVwJSXxlrCwMEK6n8VRRzj9Qvbz0Zx/+tiscVwJcFRAJVdQy1clV28wbKBoDdNuvoUpP8aR9paNeQt+cGcUMvnj8H32AV7+YTvb63eOmrQOwwYKi6XlaikpKYR36M+uvMMAPPfcc77SahKVLoup5Aot8/1sXT7PLdzmXnAVKFTrW60YNlBo3RTWFHfccQcAAkHO2p/YvXu3L7SaZMiQIX5t35eo5Are+zockhW7nTknh3QLbFZs1fpWK4YNFFoS13giMzOT0aPOI+eB7iyfWM3LLzzjI7PGUen0RiVX8N4390ApRytqSYkJpWPbcD9ZNY5qfasVwwYKX1Qqv/3Ou9lfHUa4qCH+6BoOHz7sA7PGaenkqx6o5Are+y7f5RxNDO4W5/UK39aiWt9qxbCBwheMHDmS93Y51/jf0DmfV15+UWcjk0CwfGf9aUfX32cxHj0wbKBo7RwFOBdgnXPpVDZWpxBvKaNq6yJKS/0zCz5o0ImlTIyLSq7gna/dIVmp0/wEqNe3WjFsoPAmcY0nxk8Yz7+2OUcVt/YoYNar/tmCnpeX55d2/YFKruCdb1lVHWefkkhWp1hSYwM7PwHq9a1WAlKkuP6Yy4UQOUKIbCFEs5WFvUlc44mgoCAyz7+G3Nok2luPcmjdf1o9UdoYripVKqCSK3jnGxNu46UJ/Zl30+l+NGoa1fpWKwEpUiyE6AHcD5whpUwHbm/t+3rDNddeyyu5cRxyRFFWWcUbb7wRyLc3MVGeQBUpvh54VUpZDCClPNhco94mrvFEcHAwPYdPpNMzh3htwUaeeuopqqqqfNY+OOtdqoJKrqDdt7rOzjebCzhW6ZvT1pagWt9qJVBFitMAhBA/A1ZgupTy2xMbalikODk52WORYm+L314wejR/f/xxKqsOsn//fmbMmMFjjz3msyLFdXV17Nu3T4kixXV1dYSGhipTpDg2Ntb9XfD0OS3ctI+nV1fRIyGM96/q3ejnBM79QP76nPLy8tyvMYsUe1+keAHwOWADuuAMJm08tdtckeKWMHPmTBkVGSZnTBwgp10ySJaXl/usbZWK06rkKqV23398nSM73btAPv5Vjn+FPKBa32KwIsV5wBdSylop5W5gKxDw5ILXX389V5+XwSNp27k/vZBZrwY2XZ6Jf/lxm3NB3bAeCTqb/P4IVJHi/wBnAQgh4nGeiuzy1GhLEtc0R2hoKOnnXMW22kQ6Bh3h0NrPfbauol27dj5pJxCo5ArafA+WVLGloIRQm4WszvpltFKtb7USqCLF3wFFQogcYDHwf1LKIk/t+msX3pSpU3kp15mF6I4eeTz/7NM+abdLly4+aScQqOQK2nx/3O4cTQzpGkeoLTD5MRtDtb7Vik/WUUgpv5ZSpkkpu0kpH6//2SNSyvn196WU8k4pZS8pZR8p5UfNtemvNfPBwcFkXTiVTdXJtLOWYNn1vXuSsDW4JspUQCVX0Ob747ZDAAxL0/e0Q7W+1YphV2b6k6uvvpoXc5MA+GvXvTz5+AydjUxaS0mV85Ko3oHi94phA0VrEtc0h9Vq5dIp/8fS8k7EiEpsBSvZuXNnq9oMCwvzkZ3/UckVtPnOmTyQVQ+cQ1edK5ar1rdaMWxy3aysLOnPvf1SSq669Hyqi/OYt3QLY8eO1aVmqYmJniifXLe8vNyv7QshuOXevzFv6RYAPv30U3766acWt+daiKMCKrlC8747DpZhlD94qvWtVgwbKBwOh9/fY/DgwUycOBGA3t2Tmfnkwy1+35ZUNtMLlVzBs++uQ2WMfH4pF73ysyGChWp9qxXDBopA8eSTTzLurD6sn1TJA30Kefvtt/VWMvGCRVsKAeiWEBHwbFZ/JAwbKCIjIwPyPh07dqTn4As47Iikf+h+Nn39OseOHfO6HZWSqqrkCp59F+U4L22P7JUUKB2PqNa3WjFsoHBttAkE9z3wMI9vSQXgwVP38fjfHvW6DX9n+fYlKrlC075HymtYs/cINqtguEEui6rWt1oxbKDwVYYrLURGRnL6pX9lRWUH4i1ldC39mezsbK/acO06VAGVXKFp3++zD+CQMKRbPFGhvl/y3xJU61utGDZQBJoJEycwa0cqdimYmryDpx65yxCTYyZNs2BTAQAX9knW2eT3j2EDRaAXrgghuGfGs7x1oAdBwsHUrvu9mtjs1atX8wcZBJVcoXHfqlo7WwpKCLIIzks3xvwEqNe3WjFs2WU9/pr37t2b99sM54NDDh76Yh8lVXdz4YUXEh8f3+xrfVGHJFCo5AqN+4barKx44Bxy9pfQJjxYB6vGUa1vtWLYEYWvU9Vp5aFHZvDAV0Xszj9MUVERd911l6bXubIaqYBKrtC0r81qIaNDmwDbeEa1vtWKYQOFXkRERDBr1izAWbe0Yu9qvvn6a52tTBpSUVNHWXWd3hp/KAwbKIKD9RtOjh49mvETxvPdnX2ZOyKfBW/+3Z3jsSnat28fILvWo5IrnOw7b20eAx5byKuLd+hk1DSq9a1WDBso/JHhyhtenvkyC/OcmcBnpO/l4fvu9nh8ampqILR8gkqucLLvp2vzqK5z0CHABYi1oFrfasWwgcLfm8KaIyEhgcxLbmNZeUfiLWWcbVvDd9+elDjcjUqbgVRyheN9txeWsjHvGFGhQZxnkNWYDVGtb7Vi2EBhBMZPmMDs/WmUyRAujt3Jl//+G0eOHNFb6w/NvHXOkn0X9k3WNeXdHw3DBgqrVf8vgRCCJ174J49u6QrA4+m7uffOWxo9NlB7U3yBSq7wm6/dIfnPemeNlbGZxhziq9a3WjFsoAgPN8b5Z2JiImdOvI9vS7oQY6mgj3U777333knHZWU1m/vDMKjkCr/5LttxmMKSajrHhTOgk36Ztj2hWt9qJWBFiuuPGyeEkEKIZntT7zmKhlzy50uYX9qPR3J7cMecddx0000npc5TKamqSq7wm+/mvKMAXJqZatgt5ar1rVZavTKzQZHic3EW+lkthJgvpcw54bgo4K+AptmeQCSu8Yannn+ZzMxMHA4HZWVlTJw4kWXLlrkv4wZyt2trUckVfvO95ewejOmbQnSoYRcUK9e3WglUkWKAx4CnAX2WXLaSqKgoPvroI2w2G+0T23DbaQ4evO//9Nb6w9ElPoK4SP/UfDFpmoAUKRZC9Ac6SCkXCCGaXJDQsEhxSkqKT4sUN1b8Ni0tzesixQ8++CBDSz7jnKjtxJY4+PDDj0hOdlaHWrlypRJFigGKi4uVKVKcddpA3vjPD3SLsRAaGtqqYtLg3yLF8fHx7u+tWaTYiyLFOEctS4DO9Y+XAFnNtZuent6a2qt+w263y2snXCyPPpwk5aPR8pErBstt27bJ7OxsvdU0o5KrlFK+8c0q2eneBXLau2v0VmkW1foWAxUpjgJ6A0uEEHuAwcD85iY0A5m4xhssFgvPz5rN3RudpeMe7pbLY3dPYe/evTqbaccXldECyWebneUCjXqloyGq9a1W/F6kWEp5TEoZL6XsLKXsDKwALpJS+q9oh5+JjY1l2iOv8NK+NIKEgxcydvL6qy+YiW78QPb+Y2QXOQgPtjJugDHXTvwRCFSRYq8xesWlAQMGEHXGNBaXdSLOUsaMvr/y+N8f01tLE71799ZbQTNv/OTMQXl5VgdD5Z1oCpX61ht8cp1JSvk18PUJP3ukiWNHaGyz9WJ+5ropU7jrtnV0ritj/r5IHn7/UU7tlc7YsWP1VvOIKpfw9h+t5MuN+7EKmDJUjSrhqvSttxh2ZaZeiWu85clnX+SmlV14+H3nmdRVV12FP0sh+gLXLLzReWvZbuockqwkqyF3ijaGKn3rLYYNFKpgs9n44OO57u3FiW3CeWXGre5LUSYtJzE6hKjQIC7oYowM239kDBso9Exc4y1t27Zl9uzZpHVuz8opYbyeuZXpd1xr2J2mHTp0aP4gA3DDsG6suP8czuzdWW8VzajSt95i2EChd+Iabxk8eDD/fvsD5u5PJFjYebnvFm69fhIVFRV6q51EUpLx8jg0RURIkFK+Krl6g2EDhZE2hWlhzZo1DBs2jIRz/4/PiroQJap4IT2bGyZfQU1Njd56x2H0OZRZS3Yw++fdVNU6M1ob3bchKrl6g2EDhaqMnzCBfV2uZGlpBxItJTyRtoGbrp/8u03j7msOHKvipUXbmfFlDjsOlumtY1KPYQOFERLXeEN0dLT7/u133s3isAtZU5lCB+sR7uuwiutvmGqYHbENXY3G09/lUl3n4ILe7ejdPgYwtu+JqOTqDYYNFEZJXKOVzMzM4x4/+tgTfFg5guzqRN7aFs3st+Zw4403GiJYnOhqFDblHeWzdfkEWy3cd8Ep7p8b1bcxVHL1BsMGCtfOQFVw7UJ1IYTgmRdf5cWDQ3lyrnOn3xtvvMENBggWJ7oaASklf/vSudtz8hmd6RQX4X7OiL5NoZKrNxg2UKiwMrMhrm3CDbFYLLz27ze55pprADilcxJXx67mhhumNHp8oNDzvZti/sb9rNlbTFxEMH85u/txzxnRtylUcvUG46YK+p1gsVh48803kVJyVds1DAvfTYeQMqZedw2v/fstQkLMJCwAX250Via/Z1RPokPVujT+R0AY9S93VlaWVOlSk8PhwGJpeoDmcDi44+YpTI35iT5hhyhwxHLn1kxef+NtoqKiAmjavKse2B2S77IPcEHvdiflwzSib1Oo5AoghFgrpWw2h61h/0euzEOqkJ2d7fF5i8XCi/98i/drzmdFWTuSLcXMOmUFN0650p1VKVA056oHVotgdJ/kRpPmGtG3KVRy9QbDBgrVzvVcqcs8IYTgiedfYUnMRL4rSiZWlPPmqf/jkTunsmXLlgBYOtHiGgiOlNdw64frOXDM8wZAo/hqQSVXbzBsoPi9IoTgvoens6/37bzzazJhogZbeT6nn346ixYt0lsvYEgpuWfeRr7cuJ+H/vOL3jomzWDYQKHaOoqMjAyvjr/+xmnEX/ocY76KZdb8DRw9epRRo0Yx85VX/X7Fx1tXf/D2//awaMtBokODmHFxusdjjeCrFZVcvcGwgUK1Jc+lpaVev2b0mDE89tpnpKSkAJDRLYn0vHe55pY7/JqPoyWuvmT5ziL+/pXzVOvJsX1p38ZzNjO9fb1BJVdvMGygUC1TkCudurdkZmayevVqTjsti1ljQjkndAtPxX/BNdNuYueuPb6VrKelrr7g1yMV3Pz+WuockhuGdWV0n+RmX6Onr7eo5OoNhg0UfyRSUlJYuvRH3jmSyf+K25AsjvBupy+Y8+IDvPnx/OYbUITKGjtT3l5NcUUtI3omcO+oU5p/kYkhMGygUG0hUufOnVv1+rCwMF6Z/Qkbut/JrJwogoWdx9p+RfKmF5h019+prPTdqUhrXVtKWLCVS/q3p0diJDMn9sdq0VY/VC/flqCSqzcEpEixEOJOIUSOEGKTEOIHIUSn5tpUbfdo27ZtW92GEIKbb72N/nfMZdrSNpTYbYy2rePR8A8YPOJcd+Wr1uIL15Zy84juzL9lqFerL/X09RaVXL2h1YGiQZHiC4BewEQhRK8TDluPszpYX2AezhqkHjFiZihPuEq8+YIhQ4bw+Lx13JXTj1WHw3hqXQibVi0jMzOTZ559jorq1hVH8qVrcxytqOH6d9awt+i3RERhwd79EQikb2tRydUbAlKkWEq5WErp+s1fgbOamIkH4uLieH3eQtac8iAf/jcXcE7wbvvxIybe/zSzv19r+I1zvx6p4LJ/LWdhTiH3f7ZZbx2TVhCQIsUnMAX4prEnGhYpTkpKMmSR4qaK31ZWVvqlSHGv3n1YsWIFkyZNwl6Yy8v9dhIknuWNH7cw5NthTD69EyP7d/eq+G1lZaXfixT/criO1zbVUlrjoH2k4LIOFSxbtqxFn1NMTIz7u9Dazwn8W6QY+F0WKW71pjAhxGXA+VLKqfWPrwIGSilvbeTYK3FWFRsupfR4/VO1TWH+pra2lqcfn07Eqhe59bQgrAL2OJJ4tO5aSuIzmXH5YPqmttFbk+o6Oy8t2s4/l+5ESjj7lEReuLwfMeHmjlAjEshNYc0VKXYJjQQexFl3tNlFEqotXFm6dKlf27fZbDw4/XHOfnI51/3cgU0HJZ0thbwd/BRTi5/jplf+w5crc3V1lVIy4fUVzFqyE4DbzunBG1dntTpI+LtvfYlKrt7g9yLFAEKI/sBrOIPE77Lcc6DmC/r27cub325gSY9HeWipg/JaGGNdxcP2WUw4O4sZM2ZQXl7Ot78coKis8Xjsa1dXe0IIxg1IpWt8BPOmDeGOc9OwaLwEqqV9FVDJ1RtaPUchpawTQriKFFuBt1xFioE1Usr5wDNAJDC3fhvxPilliwsYG5HGtkf7i6CgIP56x13kXz6Be+69mdMrFvL3n9ZSUeFg+vTpzP34fSovehaLNYjhaQl++0e0AAAdeklEQVSM6ZvMOacmuS9J+sLV4ZCs3H2E2T/vpme7KO46rycAkwZ2ZGxmKqE2313eDmTfthaVXL3BTFzzO2Dp0qX89a9/ZdOmTQAsnhJLUnwbXhBX8b08DRBYLYIBHWM5s0c8V5/emZgw708HKmvsrNl7hJ+2H2bBxv3sr98enhgVwv/uO5sgq2HX75k0gdY5CsMGivT0dKlSEpCNGzfqunPQbrfz9ttvM/PxB/jPn8rp3Mb5S7u2OJJZ9ktYHDkKBxasFsH7lyYzOKs/AM99v5XSqjraxYQSE2YjPNiKEILaOge928fQs50z+9ana/O477NN1Np/+760bxPG2Mz2XDmkE4lRoX77v+ndt96gkitoDxSGzZmpWuIa16UovbBarVx33XWMHz+eV158jqLvn+Ou0xwMiC3jTd5jTcG7vHWkHzUDb6K0+LfRxIer9nG4rPFKZred08MdKJKiQ6lzSPq0j+H0bnGc2yuJzI6xPpmDaA69+9YbVHL1BsMGCpOWERERwb0PPsKRm25h5nNPUr3kn9yeBVnJFtLabqDTI5P4IjKeu+++m6lTpzLjot4UHKuk4FgVZVV1lNc4A3Sw1ULXhN9S5g/s0pZfpp9PRIj5lfkjYthTj/79+8v169frraGZkpISQ1aJKioq4pUXnqbsp39hr6nihRXO0UNoENw/PIKytLFcffP/0bt3b51Nm8aofdsYKrnC7yC5rmqJa44cOaK3QqPExcXx6N+f4qH5+0i65G8kJzvzP0zqY+ORoVYej/uc7OkDufVP/Xjv3XcNmdTYqH3bGCq5eoNhA4VqiWtcy2ONSkxMDPfeey+7d+/m3nvvpSq6G/NyarEIGN/bxssDdnPaypt44sJ23H/LdaxatcowawKM3rcNUcnVGwwbKEz8Q0hICKNGjeK9xTkk/OUb/rJjBI/9VMv+Ugc946387UwYU/ExgwYNolevXjz22GPs3LlTb20TnTHszJRqiWu6du2qt4Jmunbt6t48NHz4cAoKCpjz1pvs/O41xiQXsWC7c0IzNzeXT1+dTtf1TzCrphPth13NJePGB/z/qlrf/h4xJzN9RHFxMbGxsXpraKIpV4fDwY8//sjbb7/N3LlzKS8v55lzQ7j7dGfQLqmWfLWtjnWVKbQdeBnn/WkcmZmZfl+N+HvoW6Oi/GSmaolrvNmyqzdNuVosFkaMGMHs2bMpLCzkww8/JDdqKA8trmVdgZ3oEMHEPjaeGXiIu+yvkv33oaSmpjJlyhTmzZvntzUEv4e+VR3DnnqY6EtERAQTJkxgwoQJFBcXM3/+fF794h3iDy3nwu6C0ztYqayD/fv389Zbb7Hg49nMGBHKHmsXovpcwOnnjGHIkCGEhvpvxaZJ4DBsoAgKMqxao8TFxemtoBlvXWNjY7nmmmu45pprKC0t5ZtvvuGWr+axbOMiwLnfY2RXK9OybEAeDvk66z/8F689IzkY1p3IXueSOfRcBg8eTExMjN999UQlV28w7ByFapvCVKpi7StXu93OqlWr+Oqrr9iybAG9ZC7ndLEyJNVKSNBv8xY1dknMk6VU2wXp6elcOLQvPQaM4LRBgzn11FOb/aPwR+zbQKH8prCePXtKVxoyFViyZAkjRozQW0MT/nItKipi8eLF/PTf7zn6y0K6WgsY3sm53fyst51zTgI4fE8UwVZYV2Bn4yHBkZCOWFMzad93OH37ZZKenn7cKYvZt/5D+U1hJuoRFxfHuHHjGDduHAAHDhzgxx9/ZNmyn+jXbxmbNm0iIUxyrErSJdbCsE5BDOsEkA/kU7lrPlc+V8kX2yRpaWkM7deDU3umUSIjSExMpHv37gQHB+v5X/zDYthAoVoCEJXmVALl2q5dOy6//HIuv/xyAMrKyli9ejUfr1zJljU/4chfSydbMQOSrWS0s9I11kJeicRut7Nlyxaubb+LO9KWcKxKsuX5Z/joiOSwjKU2qiOWdqcS3W0Q3bt3p3v37qSmphqiFoxK3wNvMOyph2pzFCYto7CwkLVr17J+/Xq2blzF6g3ZbN2xCyklT40MYUp/G3HhJ5/zr8q3M+gNZ60Qi4B3/hzOUaKpDktCtOlIWHIabTqmk9qlBx06dCAlJQWbzUzweyLKz1H06tVL+qoyViBYt24dmZmZemtowuiuZWVl5OTksHnzZnJyssldu4zwinza2A/RM85Cj7YWthY5uHeRcz9QxxjB3tujGm3rcIWDSZ9Wsmi3g6SkJEb3bkv/DuFYopMJjutIRGIXYlK6k5SSSrt27UhMTGxVQDF6356I8nMUqu0eddWpUAGju0ZGRjJw4EAGDhwI/DZBWF5ezrZt29i2bRvV27ZxdcoOtm/fzuG8nUz+opjObSx0jhF0bmOhQ4yF1GhBfLiFkmpn0tsDBw7QJ6OYWzqGANvADhQ4b0U/OVhb4CD1vQratm1LQkICd5/mwBIWDWGxBEUlYotOJDQ2mYj49kQldKBNXAJt27YlNjbWPXdi9L5tKYYNFCYmJxIREUH//v3p37//Sc+VlZWxZ88edu3axcY9e/hy3z725e6hrHAPBTIPIQ4ipWRFnp3Z62toF2mhXaSgXaQgMUIQF24hOsQBOLeKHys+wtQJ0UDhb29SUX/Lh5u+quRfa5ylHf98ShD3nhlGhT2I8joL89+LwG4Nw2ELx2GLZKX1NKKioomOjqZrSDERYaGERMUSEtmW0Oi2hEfHERbTloioNkRERBjyFMkngUIIMQp4CWcW7jeklE+e8HwI8A4wACgCxksp93hqMyIiwtPThiMrq9nRm2FQyRW0+UZGRtK7d+8mE/DU1NRQUFBAfn4++fn5bC8oYOn+/RQUFFD4SwGVRfmUFh9CiEqklARZ4OHFVcSHC+LCBHFhFtqGCWLDIDZUUFTx2yl7l1gLg1IEziGKHThaf4PKWsnl//jJfWz2zRH0Smh80vWlldXc/m01NpuN07uE88YFUOuwUoOVOmnFLoKwiyCkCOKj4r5UBrclNDSUwRF5dAg+hggKRgSFIIKCsQSFYLGFUBOawOH4QYSEhBBis5J8bB1BwaFYg0OxhERq/QhaHygaFCk+F2cxoNVCiPlSyoYTDFOAYilldyHEBOApYLyndmtrW1eIN9AUFhYSGam94/VEJVfwjW9wcDCdOnWiU6dOHo+z2+0cOXKEwsJCDh8+zKFDhzh8+DC7Dx/m8OHDFBUVUVRURHFUMT16HKG4uJiPso+y/NdyYkIF0SEQHSKIDhFEBQtOTCm6vsBBSbUkMlgQFiQIt0FEsPPfqvo0sbW1tVhry+jeJgKoq78dzzXv/ocdR5wjoKFjwxjW2wYOoKb+Vs+SPXVc+vZTAEQFQ8n9v2XfyitxaO4/X4wo3EWKAYQQriLFDQPFxcD0+vvzgFeEEEJ6mEmtqWk84atR+fXXX+nWrZveGppQyRUC62u1WklISCAhIUHza6SUlJeXU1xczKJFi0hLS+PYsWPu29OjSykpKaG0tJQfSkr4oqSM0tJSysrK3Lfy8jIqygVWax12u50VeXZ6vVpGSBCEBQlCgiA0SBAaBCFWKCj97Zf8/c21bCi0Y7NASJDAZoFgKwRbBduP/HacQ8KnObXYrBBkgSOV2i9kBKpIsfuY+oJBx4A44HDDgxoWKU5MTFSqSHFZWZlfihSDs/jttm3bOHjQWWSttcVvy8rK/F6k2Jefk5RSmSLFycnJ1NbWEh4ezsCBA73+nM4880w2bNjA/v37qaqqonPnzhw+fJhdu3ZRXV1NTEwMNpuNrrm51NTUYLVaiYuLI7f+cW1tLQkJCezfv5+ysjJqQmoYPTqC8vJySkpKmJFbh8ViwW6315ft3IsWAlKkWAiRXX9MXv3jnfXHFDXVbkZGhlRpy25+fj7t27fXW0MTKrmCWr4quYLxihS7jxFCBAExgMcspKqtzFQpI5dKrqCWr0qu3hCQIsX1j6+pvz8O+K+n+QnAkNmgPeEaKquASq6glq9Krt4QqCLFbwLvCiF24BxJTGjt+5qYmAQOn6yjkFJ+DXx9ws8eaXC/CrjMmzaNuOjEE4mJiXoraEYlV1DLVyVXbzBshg3VzvXS0tL0VtCMSq6glq9Krt5g2EDhuvSlCq7LrCqgkiuo5auSqzcYNlCYmJgYB8MGCpXyDoJap0oquYJaviq5eoNh81GYiWtMTPyPWQAowKgU1FRyBbV8VXL1BsMGCtUS16g0+aqSK6jlq5KrNxg2UJiYmBgHw85RZGZmStdORhWorKwkLCxMbw1NqOQKavmq5Aq/gzkK1RLX5OXl6a2gGZVcQS1flVy9wbA5MxtLXFNbW0teXh5VVVU6GHmmqqqKLVu26K2hCZdraGgoqamphl8un5+fT48ePfTW0IRKrt5g2EDRGHl5eURFRdG5c2fDbUMvLS0lKqrxlPFGo7S0lMjISIqKisjLy6NLly56K5kYHMOeejSsPemiqqqKuLg4wwUJUGuhTUhICEII4uLiDDk6O5GePXvqraAZlVy9wbCBoqlgYMQgAcb1agyXqyrORigVqBWVXL3BsIFCtcQ1KvxldqGSK4BKFeNUcvUGwwYKExMT42DYQGH0mfgTUamKtUqu4MyirQoquXqDYQOF0ScH77rrLnr16sX111/P8OHDPf7y1dTUMGzYMHdK9kBw6NAhJk+eTF5eHtddd91x61KM3rcnotJVGZVcvcGwgcLIa+Z37drFzz//TE5ODv369ePSSy/1eN4fHBzMOeecw8cff+xzl82bN3PhhRcedzt48CAJCQl07NiRu+66i5kzZx43QisvL/e5hz9x1eNQAZVcvUGtMagB2Lp1KyNHjqSurs5dLPfzzz93P3/WWWfxwAMPcO655/LQQw9RUlLCzJkzueSSS7j//vu54oorWvS+u3fv5vbbbyc/Px+LxcK7775Lz5496dOnDwsWLDjp+LKyMnbt2kVQUJBS5QNNDIqUssU3oC2wENhe/29sI8f0A5YD2cAmnAWKm237lFNOkSeSk5Pjvg/47dYcDz74oPz3v/8tq6urZVJSkpRSytLSUimllEuXLpXDhw+X7733nhw9erSsq6uTUkpZV1cn4+PjT2pr6NChMiMj46TbwoUL3cfU1NTIs88+W+7YsUNKKeVXX30lr7322ib9amtr5eTJk+WePXvk008/LRcvXnzc8y7XE/vUqKxYsUJvBc2o5CqllDgz5Tf/u67loCZfDE8D99Xfvw94qpFj0oAe9fdTgAKgTXNtDxgw4KT/lFECxUUXXSRXrVol8/PzZc+ePU96ftiwYTIzM1OWlJQc9/OUlJSTfqaFTz75RCYlJbmDSHp6upw2bZrX7TSGCoHCxH9oDRStPfW4GBhRf/9tYAlwb8MDpJTbGtzfL4Q4CCTgqgvfBEY+j87OziY9PZ3q6mr33ERZWRmRkZFs3ryZgoIC4uPjT1rSXV1dfdKK0zPPPLO+BuTxPPvss4wcORKAjRs38vjjjzNlyhSf+LtcVcFV01UFVHL1htYGiiQpZQGAlLJACOGxqIEQYiAQDOxs4nmPRYqrqqooLS3FarXicDiOm/CMioqivLwch8NZvTk8PJza2lr3bL9r2bLrFzsoKIjQ0FB3G0IIIiMjKS8vd//iRkREuAu/utooKyvDarVit9sJDQ3Fbrdz6NAhbDYbO3bs4IorruCDDz7gnnvu4fPPP+fiiy+murqawsJC4uLiAOdVEFdh24ULFxIcHOwOjBaLhYiICMrKytwe7dq14+uvv2bcuHFYLBa2b99Oenq6e+OczWbDZrO5s4K52mgYgCIjI6msrMRut2O326mrc1bNrqqqYsmSJYYuUlxRUaFMkeLi4mK3a2uLSQNkZGT4vZi0JpobcgCLgF8auV0MHD3h2GIP7SQDW4HBWoY6aWlpJw2TjDBM/t///ifHjRvnfnzdddfJhQsXygMHDsjBgwfL77//XkrpnKsYPHiw+7i5c+fKO++8s0XvWVFRIceOHSvT0tJkRkaGvOKKK1r1f2h4+mOEPm2OE+dYjIxKrlIGbo5iK5AsGwSCJo6LBtYBl2ltOzMz86T/lBG/1OvWrZNXXnmltNvtHo/785//LHNzcwNk5ZmGrkbs0xOpqqrSW0EzKrlKqT1QtHYdRcPiw9cAX5x4QH3h4s+Bd6SUc7U27BrOGZ3+/ftz1llneUwGXFNTwyWXXGKYnYWq9K2L3bt3662gGZVcvaG1geJJ4FwhxHbg3PrHCCGyhBBv1B9zOTAMuFYIsaH+1q+5hlXKcHXddde5Rk6NEhwczNVXXx1AI88EcoWoL3DNG6iASq7e0KrJTCllEXBOIz9fA0ytv/8e8F5r3sfExERfDLuEW6UEpdB4oh2jopIrOK9AqIJKrt5g2EDhaShvRFTyVckV1KrxopKrNxg2UKiWXEWlCUKVXAH3WgYVUMnVGwwbKExMTIyDYQNFcHCw3gpeoVKiHZVcwbkqUhVUcvUGw24zN+KXuaioiHPOcV7kOXDgAFarlYSEBABWrFihuZ26ujri4+M5etTjdhe/oVoQTk1N1VtBMyq5eoNhRxRG3BQWFxfHhg0b2LBhA9OmTeOOO+5wP/bnug9fr3swYt96wrW/QgVUcvUGw44otND5vq+afO4ff+7DpEEdAfhg5T4e+Hxzk8fueXJMq10uv/xyDh48SFVVFXfccQdTp07ltddeY8eOHTzzzDMA/POf/2T37t384x//cL/O4XBw99138/333yOE4NFHH2XcuHEsWrSIJ598kvj4eLKzs9m8uWl/ExN/Y9hAoVp9hNdff50OHTpQUVFBVlYWY8eOZdKkSfTr148nnniCoKAgZs+ezZw5c4573dy5c8nJyWHjxo0cOnSI0047jWHDhgHO05mcnBw6duzoU1eLxbADyUZRaUu8Sq7eYNhAER4e3uwxWkcCkwZ1dI8u/MXrr7/O/PnzAWfpw507d5KVlcWwYcP45ptv6Nq1K1arlV69eh13KrFs2TImTZqE1WqlXbt2DB06lDVr1hAcHMyQIUN8HiTAuX1eJbKymi22bRhUcvUGw/5pUek8etGiRSxevJgVK1awceNG+vbt614HMnXqVObMmcNbb73F5MmTT3qtp8VP/vqFNnLi4sZQKWGtSq7eYNhA4UpAowLHjh0jNjaWsLAwsrOzWb16tfu5M844g507dzJ37lzGjx9/0muHDRvGRx99hN1up7CwkJ9//tnvf5VUW5mp0gIxlVy9wbCnHioxZswYZs2aRUZGBqeccspJqdDGjRtHbm4uMTExJ7123LhxrFixgoyMDIQQPP/88yQmekwUZmIScIRR/7oMGDBArl279rifbdmyhVNPPVUnI8/I+jRkjTFq1Cjuv/9+hg8fHmCrxmnoauQ+dVFXV6dMdTOVXAGEEGullM0OYQ176qHaEK6xvSlFRUWkpaURGxtrmCAB6u2j2bZtW/MHGQSVXL3BsKFPpcQ10PiiqLi4OEN+cVRLXHPw4EFltm+r5OoNhh1RmJiYGAfDBgozcY3/UMkVnGnvVUElV28wbKAw6iRrU6jkq5IrqDVfpZKrNxg2UKg24abSF0QlV8BdQEcFVHL1hlYFCiFEWyHEQiHE9vp/Yz0cGy2EyBdCvNKa99Qbq9VKv3796N27N5dddpnHNP0nMmfOHG655RY/2pmY+IfWjijuA36QUvYAfqh/3BSPAUu1NmzUnAlhYWFs2LCBX375heDgYP71r38Bxsyf0RQquYKzFJ8qqOTqDX4vUgwghBgAJAHfAprWJ2v6Mk8/eaWjmwtfhKz6vRVrZsOC2z20c0yL0kmceeaZbNq0CYDx48eTl5dHVVUVt912GzfccAMAs2fP5oknniA5OZm0tDRCQkIAOHToENOmTWPfvn0AvPjii5xxxhkt8vAW1QJFUlKS3gqaUcnVG/xepFgIYQGeA66ikRogJxyruUhxWFgYja+DdOKQDmqqqqitrcVWXYWnef7S0tLjihS79pk0VqTYdXxdXR0LFizgggsuoLS0lBdffJGEhASsVisDBgzgvPPOo7a2lkcffZRly5YRERHBmDFjyMzMpKamhptvvpkbb7yRYcOGceDAAc4//3zWrFlzXJFi16RjZGQkVVVV7vUPoaGhSCndcw0tKVIcGRmpTJHiwsJC90pSoxcpXrZsmfuqklmk2IsixcAtwD31968FXtFS69CoRYotFovMyMiQGRkZ8pZbbpHV1dVSSinvu+8+2bdvX9m3b18ZHR0tly9fLj///HN51VVXuV/70ksvyb/85S9SSikTEhLc7WRkZMiUlJTjigf7E7NIsf9QyVVK7bVHmx1RSClHNvWcEKJQCJEsnaOJZOBgI4cNAc4UQtwMRALBQogyKaWn+QzDJq5xzVE0ZMmSJSxdupTly5cTHh7OiBEj3Fdtmtr/4XA4WL58uS7rRVRLXBMdHa23gmZUcvUGvxcpllJeIaXsKKXsDNyNs1ixxyAB2hLXGIVjx44RFxdHeHg4ubm57kS7gwYNYsmSJRQVFVFbW8vcub/VaD7vvPN45ZXfLgCdGHz8iWqJazIzM/VW0IxKrt4QiCLFLUKl5CqjRo2iqqqKvn378vDDDzN48GAAkpOTmT59OkOGDGHkyJHHfYlmzpzJmjVr6Nu3L7169XJfPQkEDecuVGDZsmV6K2hGJVdv8HuR4hN+PgeYo7Ht1qj5jcYCWEhICJ999hlRUVEnPTd58uRGM1vFx8fz8ccf+8Xx94ZKm9hUcvUGtU5WTUxMdMGwgaKxv85GRqXsyyq5Au6s5Cqgkqs3GDZQuK5vn4hRT0ma8jUiLlej9uWJZGdn662gGZVcvcGwgaKxc73Q0FCKiooM+QVXqdy93W5HSklRUZESW85di49UQCVXbzBshqvGSE1NJS8vj0OHDumtchJVVVVK/NLBb66hoaG/21qZJr7FsIGisXUUNpuNLl266GDTPMXFxcTGNrl51lCo5ArOZcyqoJKrNxj21EOloTyotTZBJVdQy1clV28wbKBQLbmKa9OOCqjkCmr5quTqDYYNFCYmJsbBsAWAhBClwFa9PbwgHjist4RGVHIFtXxVcgXoKaVsdtGSYSczga1SQwUjoyCEWKOKr0quoJavSq7g9NVynHnqYWJi0ixmoDAxMWkWIweK1/UW8BKVfFVyBbV8VXIFjb6Gncw0MTExDkYeUZiYmBgEM1CYmJg0iyEDhRBilBBiqxBihxCi2fyaeiKEeEsIcVAI8YveLs0hhOgghFgshNgihMgWQtymt5MnhBChQohVQoiN9b4z9HZqDiGEVQixXgixQG+X5hBC7BFCbBZCbGjuMqnh5iiEEFZgG84cnHnAamCilDJHV7EmEEIMA8pwJg02dCnr+kzpyVLKdUKIKGAtcImB+1YAEVLKMiGEDVgG3CalXKGzWpMIIe7EWeQqWkp5od4+nhBC7AGypJTNLhAz4ohiILBDSrlLSlkDfISzhoghkVL+CBzR20MLUsoCKeW6+vulwBagvb5WTVNfesKVpNRWfzPWX7YGCCFSgTFAqxJLGxEjBor2wK8NHudh4C+zqgghOgP9gZX6mnimfii/AWfNmIVSSiP7vgjcAzj0FtGIBL4XQqytr9LXJEYMFI1VzDHsXxEVEUJEAp8Ct0spS/T28YSU0i6l7AekAgOFEIY8vRNCXAgclFKu1dvFC86QUmYCFwB/qT+NbhQjBoo8oGFJ6FRgv04uvzvqz/U/Bd6XUn6mt49WpJRHcRbBHqWzSlOcAVxUf97/EXC2EOI9fZU8I6XcX//vQeBznKf9jWLEQLEa6CGE6CKECAYm4KxIZtJK6icH3wS2SCmf19unOYQQCUKINvX3w4CRQK6+Vo0jpbxfSplaXxFvAvBfKeWVOms1iRAion5CGyFEBHAezprCjWK4QCGlrMNZ2Pg7nJNtn0gpDZvaWAjxIbAc6CmEyBNCTNHbyQNn4Kwqf3b9JbENQojRekt5IBlYLITYhPMPyEIppeEvOypCErBMCLERWAV8JaX8tqmDDXd51MTExHgYbkRhYmJiPMxAYWJi0ixmoDAxMWkWM1CYmJg0ixkoTExMmsUMFCYmJs1iBgoTE5Nm+X+8oUM8o+NEpgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.plot(xx,yy,'k',lw=3,label='$f(x)=' + sb.latex(f)+'$')\n",
"plt.plot(xx,yy_Taylor,'--',label='Taylor')\n",
"plt.plot(xx,yy_Pade,'--',label='Pade')\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.ylim([-0.5,1])\n",
"plt.gca().set_aspect(1/plt.gca().get_data_ratio())\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# Sympyで数式を表示するときに浮動小数点数を丸める\n",
"# https://stackoverflow.com/questions/43804701/round-floats-within-an-expression\n",
"from sympy import preorder_traversal,Float\n",
"def RoundFloat(f,order=2):\n",
" for a in preorder_traversal(f):\n",
" if isinstance(a, Float):\n",
" f = f.subs(a, round(a, order))\n",
" return f"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$f_{\\rm{Taylor}}(x) =\\frac{x^{6}}{720} - \\frac{x^{5}}{120} + \\frac{x^{4}}{24} - \\frac{x^{3}}{6} + \\frac{x^{2}}{2} - x + 1$$"
],
"text/plain": [
"<IPython.core.display.Math object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$$f_{\\rm{Pade}}(x) =\\frac{0.03 x^{2} - 0.33 x + 1.0}{0.03 x^{3} + 0.2 x^{2} + 0.67 x + 1.0}$$"
],
"text/plain": [
"<IPython.core.display.Math object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"text_Taylor = r\"$f_{\\rm{Taylor}}(x) =\"+ sb.latex(f_Taylor)+'$'\n",
"text_Pade = r\"$f_{\\rm{Pade}}(x) =\"+ sb.latex(RoundFloat(R,2))+'$'\n",
"display(Math(text_Taylor))\n",
"display(Math(text_Pade))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x1152691d0>"
]
},
"execution_count": 11,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAisAAAD8CAYAAAC2CuSWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XlUU9fCNvBnJ4GEQEDmUQRkUEFRRrEU1Dpbbau2jm1pFYdr7WBb9baV61Bvp+tQx69WW9tq7XW6FX2tWqtIqYgMigoyKCAiIgoIBEKA5Hx/RFJEhkBOSI7u31quRZKTvZ/sINk5Zw+EYRhQFEVRFEUZKp6+A1AURVEURbWHdlYoiqIoijJotLNCURRFUZRBo50ViqIoiqIMGu2sUBRFURRl0GhnhaIoiqIog8ZKZ4UQ8h0hpJQQcrWNxwkhZCMh5Doh5DIhJICNeimKoiiKevKxdWZlF4Ax7Tw+FoDXw39zAWxjqV6KoiiKop5wrHRWGIaJB1DeziEvAPiRUTkPoAchxJGNuimKoiiKerIJuqkeZwC3mt0uenjfneYHEULmQnXmBRIRP1AiFsLE0hFisRh8Ph+1tbUAAIFAABMTE1RXVzc9D2ZmZqitrYVCoQAAmJqaoqGhAfX19QAAkUgEQghkMhkAwMjICEKhEFKpFADA4/Fgamra5TIYhgGfz4epqSlqamqgVCoBAGZmZpDL5WhoaAAAmJiYgGEY1NXVAQCMjY1hZGSEmpoaAACfz4dYLGaljKZcACCRSCCTydDY2AiGYWBqagqFQgG5XA4AEAqF7bZxyzIAQCwWd6qMrr5PhBD149q+T83LqJLW4GZlIwgAXydzVt6n2tpa9W1t36eutHFny8jJybnPMIwtKIqiDFh3dVZIK/c9ts4/wzDbAWwHgCAnPvPfNxyQ4L0cr7/+uq7zaS0uLg5Dhw7VdwyNcCkroLu8SXllmLr9PPx79sDhhc+wUibX2pYQclPfGSiKojrSXbOBigD0bHbbBUBxe09QgI/eRveRcOqIToOxxd/fX98RNMalrIDu8soblfC2N0M/R3PWyuRa21IURXFBd3VWYgG89nBW0GAAlQzD3GnvCVWNRgAAF2Wh+nS+IWt+Gt7QcSkroLu8Ed62OPleJD6b1J+1MrnWthRFUVzA1tTlvQASAfgQQooIIbMJIfMJIfMfHnIMQB6A6wC+BfCPjsqsalBdORrrUoPTp0+zEVOn8vLy9B1BY1zKCnArL5eyUhRFcQUrY1YYhpneweMMgIWdKVPJM4acIQgS3cY7R/+H8ePHa5WRoppjGAal1XLYSYQgpLUhVRRFUZShMNgVbHv0sMSfUhfwCAOjsgz1DAtD5ebmpu8IGuNSVkA3ee9J5Qj99x8Y+p84VsvlWttSFEVxgcF2VszNzfHTNWOszPHC0fM3cP78eX1HapeVlZW+I2iMS1kB3eTNuqMaW2InEbJaLtfalqIoigsMtrNSW1sLIyd/rNibitxb93Do0CF9R2pXWlqaviNojEtZAd3kzSqpAgD0cWBvJhDAvbalKIriAoPtrADA5MmT1T8fPHhQvXAWRWkrq0R1ZqWPo0TPSSiKoqiOdNeicJ0mEAgwfPhw9O7lhDefdYTYCLh48SICAgxzD0RLS0t9R9AYl7ICusnbdBmI7TMrXGtbiqIoLjDYMysmJiYQCoUYMXwoPuqdi7k9b+J/B/frO1abuLQYGJeyAuznbVAocb1UtfS+jwO7Z1a41rYURVFcYLCdlabFtUZPmIKLdU4Qk3qU5/yl51RtO3v2rL4jaIxLWQH28+bfr0G9QglXKzHMhOyeXORa21IURXGBwXZWmowePRqxRapT9UMdapCZmannRK3j0ngaLmUF2M/raiXGvnlhWDGxH6vlAtxrW4qiKC4w+M6KWCzGbV4vAMAYqyLs++9ePSdqHZcWFuNSVoD9vCIjPkLcrTC8jz2r5QLca1uKoiguMNjOikTy91iCkRNexhW5IySkDqUZZw3y22tkZKS+I2iMS1kBbuXlUlaKoiiuMNjOSvPNC8ePH49DhRYAgBFOUly5ckVfsdqUnp6u7wga41JWgP28H+5Px39OZKO2vpHVcgHutS1FURQXGGxnpbHx7w8SsViM20YeOFTeGz9elOG///2vHpO1rqKiQt8RNMalrAC7eStq6rE/tQg7E/IhFPBZK1ddPsfalqIoigsMtrPS0rgXp2Lypos4nJCFX375xSAvBVGGL6NYtXJtPydz8Hl0fAlFURQXGGxnRSwWP3J7zJgxMDdXzQrKy8tDamqqPmK1yVAXq2sNl7IC7Oa9WlwJAPBzYncxuCZca1uKoiguMNjOikKheOS2SCTCCy+8gCmRfbHvrQDs27tbT8laV15eru8IGuNSVoDdvFdvqzorvs4WrJXZHNfalqIoigsMtrMil8sfu2/69OlYHMLDy9bX0XDrwmMdGn0qKCjQdwSNcSkrwG7epstAvjo6s8K1tqUoiuICg+2stGbkyJE4dFM1pfnFXlLExcXpNxDFKdV1Dci/XwNjPg9ednQDQ4qiKK4w2M6KUCh87D6BQACFQwDkjADPigtx8Jcf9ZCsdR4eHvqOoDEuZQXYy1tbr8CLA50wsp89jAW6+dXnWttSFEVxgcF2Vvj81qeVTpv5Ok5VuoJHGNjWZj2yHos+NV/EztBxKSvAXl57cxE2TBuELTN1NwiWa21LURTFBQbbWamtrW31/uDgYBwrMgMATHGtxJEjR7ozVpu4tBgYl7IC3MrLpawURVFcwUpnhRAyhhCSTQi5TghZ1srjroSQM4SQi4SQy4SQcVrUBTu/YahkTNBfeAcnD/+sXXjqqZF4owylVXX6jkFRFEV1ktadFUIIH8AWAGMB9AMwnRDScjvbTwDsYxhmEIBpALZ2VK5AIGjzsVmvvo4teS74NNcLv52KR0lJSZfzs8Xa2lrfETTGpawAO3ll9QrM3HEeYZ+fRl2D7maRca1tKYqiuICNMyshAK4zDJPHMEw9gF8AvNDiGAZA01xRCwDFHRVqYmLS5mO9e/fG8ZtiLP85FcWlFdizZ0/XkrPI19dX3xE0xqWsADt5r5VUQckAnrZmEBmxv8x+E661LUVRFBe0ffpCc84AbjW7XQQgtMUxKwCcJIQsAmAKYERrBRFC5gKYCwB2dnbqqckeHh6QSCTq8QDW1tZ47bXX8OeffwIAtm7disWLF+PixYuoqlKtoxEUFIS7d+/i1i1VNC8vLwiFQly9ehVN5Xt7eyMhIQGAavZRWFgYUlJSIJVKAQChoaEoKirC7du3AQA+Pj7g8/nIzMwEADg4OMDd3R2JiYmQSqWwtbVFaGgokpKS1AN/w8LCkJ+frz77069fPygUCmRnZ6saz9kZLi4uSEpKAgCYmZkhKCgIiYmJ6rVmwsPDkZOTg9LSUgCAn58f5HI5cnNzAQA9e/aEvb09UlJSAADm5uYICAhAQkKCeo+liIgIZGRkoKysDFKpFM888wyqq6uRl5cHAHBzc4OVlRXS0tIAAJaWlvD398fZs6pdrgkhiIyMRHp6unr/m4CAAJSXl6vXFmntffL19UV8fDwA1dmy8PBwpKWldep9unv3LgghWr1PKSWqDoqdQIakpCSdvU+nTp2CSCRi5X0CAH9/f52/TxRFUYaOaLvHDiHkZQCjGYaZ8/D2qwBCGIZZ1OyYxQ/rWksICQOwE4AfwzDKtsr18fFhmj4oWlNdXQ3fPr0xJ9IVliZA2IL/h6CgIK1eizbi4uIwdOhQvdXfGVzKCrCT951fLuLwpWL8+6X+mBHqyk6wVnCtbQkhqQzD6O8/DkVRlAbYuAxUBKBns9suePwyz2wA+wCAYZhEACIANu0V2vRNui0SiQTjx47Bx17XMdclH3t+2NnZ3Kxqb4yNoeFSVoCdvJduPQAA+PfUzTL7TbjWthRFUVzARmclGYAXIcSdEGIM1QDa2BbHFAJ4DgAIIX2h6qzca69QMzOzDit+ZeYbOFvTC0LSCJP7aXpdcyU8PFxvdXcWl7IC2uetqKnHzbJaiIx48LHX7TooXGtbiqIoLtC6s8IwTCOAtwCcAHANqlk/GYSQVYSQiQ8Pex9ANCEkHcBeAFFMB9ef2lpnpbnIyEgczFd1ama5V+DA/gNdfyFaahpDwAVcygponzf7bjUAoL+zBQR83S4txLW2pSiK4gJWzlkzDHMMwLEW98U0+zkTwDOdKVOTTQp5PB6cBo3F/foC9DO+i7VH9+DV117tTDWsaRowygVcygpon3ewhzUurxiFcmk9S4naxrW2pSiK4gKDXcFWU2/Ono2fbzsBAIbZ3kdWVpaeE1GGyFxkBDcbU33HoCiKorrAYDsrpqaafbA4Ojrislw1u2OSzU18t+MbXcZqkz5nInUWl7IC3MrLpawURVFcYbCdlYaGBo2PnTRjNtYWeOHVOFv8+NMe9boX3enu3bvdXmdXcSkroF3ewrJaDP73H1h28DKLidrGtbalKIriAoPtrNTXaz6+YPTo0fj6dCkOxV/D3dJ7OHCg+wfaNi1qxgVcygpol/firQqUVNXhXnX3dGC51rYURVFcYLCdlc7g8/mYN2+e+vaWrVv0mIYyJKk3Vau4DnLtoeckFEVRVFcZbGelaclyTc2ZMwf9ejvhf28PwqJBcly6dElHyVrn5eXVrfVpg0tZAe3yphSoOitBblZsxWkX19qWoiiKCwy2s9LRCrYt2dvbI/zZSEzokYfJ1vnY+c1mHSVrnVAo7Nb6tMGlrEDX81bXNSCrpAoCHoG/S/ecWeFa21IURXGBwXZWurIa7WtzFuJktRuMiQKONZfx4MEDHSRrXdPGe1zApaxA1/NeLHwAJQP4OVvAxFh3Oy03x7W2pSiK4gKD7ax0xZAhQ3Dwpuob9Bu9SrBzx7d6TkTpU8rD8SpBvSz1nISiKIrShsHuumZkZNTp5xBCEDp6BnJvFsHL6B4Kzv+Kxsb3umVzOTs7O53XwRYuZQW6nnesnwNMjPgIce++zgrX2paiKIoLSAdb9OhNYGAgk5qa2unnyWQyrHhjOL7ok4nUOmcUBK/G5MmTdZDwUY2NjZzZcZdLWQFu5eVSVgAghKQyDENXsqMoyqAZ7GUgqVTapeeZmJhA6DkMD5RiDBQW45dd21hO1rqEhIRuqYcNXMoKcCsvl7JSFEVxhcF2VrQx/x9vYXa8Nby/M8KBo38gJSVF35GobnYyowSb/shF7sMdlymKoijuMtjOCo/X9WhOTk4Q9xyEvKL7AIC1a9eyFatNXJqyyqWsQNfyHkq7jbW/56gXhesuXGtbiqIoLjDYzoqmGxm25b333gMAEBBkpv6J/Px8NmK1KSwsTKfls4lLWYHO51UqGZzPL1M9t7e1LiK1iWttS1EUxQUGOxKwtrZWq+cHBARg3JhRWBuQD1dBBT5Z/xXWbdzKUrrHpaSkcGbHXS5lBTqfN6ukGg9qG+BkIYKrlViHyR7HtbbVRmpqqp1AINgBwA8G/MWHoiiDxxBCKhUKxfdKpXJbYGDgY5sDGmxnRaFQaF3Gu4s/QPGJD9DH6C5sHqTg/v37sLGxYSHd47o6IFgfuJQV6HzexDzVWZXBva07vRKytrjWttoQCAQ7HBwc+tra2lbweDzDnFZIUZTBYxgG9fX1RsXFxYuqqqoCALze8pgn+tvQiBEjsDtPtSfMXLfb2Lxpg54TUd0h8cbDS0Ae3XsJ6CnkZ2trW0U7KhRFaYMQAqFQ2NCrV69KAOGtHWOwnRVtx6wAqgZ4btIcpMudYMOToi77FKqrdTM7JDQ0VCfl6gKXsgKdy6tQMkjS03gVgHttqyUe7ahQFMWWh39PWt0bxWA7Kw0NDayUM3XaVPy/HNXZlUVed7B1yxZWym2pqKhIJ+XqApeyAp3LK61rxPA+dgjqZQkXy+4drwJwr20piqK4gJXOCiFkDCEkmxBynRCyrI1jXiGEZBJCMgghP3dUZn39Y+NrukQgECBg9OvIarCHM/8B7qX9qvXg3dbcvn2b9TJ1hUtZgc7ltRAb4etpg3BgwRAdJmob19qWoiiKC7TurBBC+AC2ABgLoB+A6YSQfi2O8QLwTwDPMAzjC+BdbevtjNejorA5yxr3lBJIZXXYsWNHd1ZPURRFUZQW2DizEgLgOsMweQzD1AP4BcALLY6JBrCFYZgKAGAYprSjQkUiEQvRVIyNjeETOR29vrqHb46m44svvkBdXR1r5QOAj48Pq+XpEpeyAprnlTcq8NuVO6iUsXMJsSu41rZPuuLiYsGQIUO8zM3NB77yyiu9Fi5c6Lxq1SqNdpvs379/35SUFPb+EHVBXV0dGTBgQB+pVEr27NljMWPGDFd95qEofWFj6rIzgFvNbhcBaDnK0BsACCF/QTV4ZgXDMMdbFkQImQtgLgA4OjoiLi4OAODh4QGJRIL09HQAgLW1NXx9fREfH696EQIBwsPDkZaWhqqqKgBAUFAQ7t69i1u3VNHGjhuHT9esgayuFMXFxVi5ciVWr16t3stFKBQiLCwMKSkp6umnoaGhKCoqUp/a9/HxAZ/PR2ZmJgDAwcEB7u7uSExMRGNjIwoLCxEaGoqkpCTIZDIAqkXC8vPzUVJSAgDo168fFAoFsrOzVY3n7AwXFxckJSUBAMzMzBAUFITExETI5XIAQHh4OHJyclBaqurj+fn5QS6XIzc3FwDQs2dP2Nvbq7cVMDc3R0BAABISEtDY2AgAiIiIQEZGBsrKytDY2AiRSITq6mrk5eUBANzc3GBlZYW0tDQAgKWlJfz9/XH27FkwDANCCCIjI5Geno6KCtWqsAEBASgvL0dBQQFr75OXlxeEQiGuXr0KQLWLsaWlpfp3ob336ffLhfgyuQ5etibY86pfq+8ToNo/SlfvU1FRkfo52r5PAODv76/z9+lJFhMT4+Dh4SE/d+5cbnFxsWDgwIH98vLyrmjy3Hfeeafk448/dj5x4sQNXedsi0gkYqKiou5HR0e73r9/X3Ds2DG9ZaEofdJ612VCyMsARjMMM+fh7VcBhDAMs6jZMUcBNAB4BYALgD8B+DEM86Ctcn18fJimP/ps2bRpEz7+aCk+mNAPd2QCrN1zGmIxO4Mw4+LiMHToUFbK0jUuZQU0z/vZb9fwzdk8zI3wwEfj+uo+WCu41rba7Lqcnp5e4O/vf5/tTGzy9vbut3bt2lsTJkyoXr58uX1ubq7ol19+uanJc2tra4mjo6P/5cuXM3r16qWz03WFhYWCl156qXfL+//3v//dcHV1bUxNTRWFhIT4Xr9+/bK7u7v+ThtSVDdIT0+38ff3d2t5PxuXgYoA9Gx22wVAcSvHHGYYpoFhmHwA2QC8WKi7U6Kjo/HaKH/EeOfin753sXXL5u6OQOlQfI7qczPCy1bPSSh9q6urIxKJZGBubq7J1KlTPb29vfv9/vvvFhEREeq1C+bPn+8ycuRIdSdh3rx5LmFhYd5yuZwAgFgsZnx9fWtiY2PN2ci0du1am969e/tKJJKBERERXrdv3xYAwMMOSXbLf66uro0lJSX8RYsWuS5duvT2zp07rdjIQVFcxMZloGQAXoQQdwC3AUwDMKPFMb8CmA5gFyHEBqrLQnntFWpkZMRCtEeJRCL4Pvcqcm6vh7dRKe6l/g/V1QsgkUi0LtvBwYGFhN2DS1kBzfKWVtXh2p0qiIx4CHKz7IZUreNa27KFEBLYXXUxDJPa0TEikYiJi4vLGjVqlE9ZWVk6AFhaWvr7+vqqB6utWLHijpeXV/9z586ZJCQkmJ05c8Y8MTExWygUqk83e3t716Wnp5u0LH/YsGGeKSkpZq3VHRQUJD1z5sz15vctW7bM4ciRI5aHDx++7unpWf/666+7fvDBB8579+5t8yyPTCYj06ZNc9+wYcOtQYMGyQIDA/t8+OGH9ywsLJQdvX6KetJo3VlhGKaREPIWgBNQjUf5jmGYDELIKgApDMPEPnxsFCEkE4ACwIcMw5S1V66udq+dPWcO3pu+B1v6l+I9ryKs+8+X+NfK1VqX6+7uzkK67sGlrIBmeeNzVWdVwjysITJqdU2hbsG1tn2SJScni/v06aNep6C6uppvbm6u/qB3cHBQREdHl0ZFRblLpVJ+fHx8lrW19SP7fEgkEmVJSclj35xadkbac/v2bcHXX3/tmJSUlOnn5ycHgOjo6PuLFi3q1d7zTExMmLi4OHU9ly9fztK0Top60rCyzgrDMMcYhvFmGKY3wzBrHt4X87CjAkZlMcMw/RiG6c8wzC8dlamrPVaMjY0R9PwcXJY7woFfBV7eSfXAVW00Dd7kAi5lBTTLG59zDwAQ4a3fS0Bca9sn2aVLl8R+fn6yptvm5uaKqqqqR/7mBQQE1Obm5prExMTc9vT0fGw8SHV1Nc/CwqJRmxxHjx41b2ho4D3zzDN9JRLJQIlEMvCll17ykkgk2m+ARlFPCYPdyFCXXnvtNURP2YXv/O/gbY+bWLlmJdZ9rZuVbanuUVWn+pzRd2flaaXJpZnudvXqVZMZM2aoz+D26dNHlpmZKYqMjKwFgAsXLpgsXrzYddKkSWU//vijzfz588tblpGTkyOaPn36Y/dHRER4tXcZKD4+Prfpdnl5OX/kyJEVv/32W7uXvimKapvBLrfP4+kuGp/Px6TZH+JsTS9YEBmM7iThxg3tZgSamDx2WdtgcSkroFneXW+E4MJHz8HDRvs9pbTBtbZ9kmVlZYmDgoLUZ1ZGjRpVGR8fLwGA/Px8o5deeslz/fr1N3ft2lWYnZ1tcvTo0UcGr8lkMpKRkWE6YcKEqpZlx8fH59bW1l5s7V/zjgoABAcH154/f16SkJAgBoDy8nLe7t27eyiVdOgJRWnKYDsrbGxk2J7x48fj2+sOeDnOGV/uT8XSpUu1Ko9LG9hxKSugeV47cxEIITpO0z6ute2TqrCwUFBVVcX39/dXD6idN29e2enTpy3Kysr4Y8aM8VqwYMHdmTNnVkokEuXChQtLYmJinJuXsXfv3h6hoaHVbm5uWk0XHjFiRM37779/Z+rUqb3FYvGgvn37+v3222/muvxCRlFPGq3XWdGVvn37MteuXdNpHefPn0dYWJj6dnx8PJ599tkulZWUlMSZDyouZQU6znu9VIretqZ676gA3GvbJ32dlZbeeustZzs7u4aYmJgOB6oNGDCgz86dOwuCg4PZXe6aoqg2tbXOisGOWemOU6SDBw/G9OnTsXfvXvh5OmLj58vxzDOnu3QJqmklVC7gUlag/bx596QYse4s+jtbIPatZ/TeYeFa2z5tNm/erPFOk3T2DUUZjqf+POTnn3+OKcP64+IMGT7qfxc//PCDviNRnXDq2l0AMJgzKxRFURT7DLazYmbW6kB71rm6usJn8FjcV5phkKgYl49tR2VlZafLaX45ydBxKSvQft5Tmaqz+SP62XdXnHZxrW0piqK4wGA7K02bw3WHZR8tx5prLgCAj/sWYs2qf3W6jPz8fLZj6QyXsgJt5y2vqUfKzXIY8QkiDWTKMtfalqIoigsMtrPS0NB9+3WZmZlhyKS3cV7WEzY8KTyq/0JGRkanymjarZcLuJQVaDvvyYwSKBkgrLcNJCL2t2foCq61LUVRFBcYbGelu02bPg1br7tAwRDMcbyOL2Leh6HOlKJUjl6+AwB4vr+jnpNQFEVRumSwnZXuXlyLEIIlK/+D70q8ICBKzPEo7tRg2379+ukwHbu4lBVoPW9dgwLX7lRBwCMY5WsY41UA7rUtRVEUFxjs1GV9nNXw8/PDnh6R+PmeEp8cLkRV3Qd4/vnnYWNj0+FzFQrubPPBpaxA63lFRnyc/+g5ZBZXoYfYWA+pWse1tqUoiuICgz2zUlenn3WYPolZiY/+rwz5t++jrKwM77//vkbPy87O1nEy9nApK9B2XiM+D/49e3RzmvZxrW0piqK4wGA7K/piamqKrVu3AgAICGpvJuO3Y8f0nIpqrra+EVK5VhvhUhRFURxisJ0VY2P9ndofN24cpk6bihOLB2D/0Ns4uvNTVFU9tpfZI5ydndt93JBwKSvweN4DqUUIXP07tpy5rqdEbeNa21JPl4yMDOGSJUscV61aZSeXy7t9FcWNGzdav/XWW84bNmyw7u66Dxw4YD58+HBPXdZx5coV4Zo1a+yioqJ6lpSU8HVVz8GDB82XL19uP23atF76eB/1wWA7K0ZG+p2KumnjJvxeJAIArPS9ieXLPmj3eBcXl+6IxQouZQUez3swtQjyRiV6Won1lKhtXGtb6umyZcsWG7FYrJTL5Twej9ftAwPNzc0VAoGAkUqlOvsgb01qaqrowYMHfDc3N50u4NW/f3+5nZ1dQ1FRkbEuP8MmT55ctXr16rumpqZK2lnRs5qaGr3Wb2tri4AX30FCjStseFIMN0rBiePH2zw+KSmpG9Nph0tZgUfz5t6tRnpRJSQiAUYZyKq1zXGtbaknU11dHRkwYEAfqVRK9uzZYzFjxgxXAJDJZLzJkyc/cHNzkx87dkzS3RmioqIebNiwobixsZFkZmbq5PR5a/UePXrUvLi42OjKlSvixMREVqaattXG0dHRFdHR0fdv3Lih9etrqw4A+OKLL2zHjh1baW5urvuN9AyAwc4GMgRTp03DnKMHMNDzLl6wvIG3vl2F4JAQWFlZ6TvaU+tAWhEA4PkBjhAZdeuXM4qDiouLBVOmTHG/evWq6ZgxYyr27dt3szvqXbhwobO9vb1Guzv379+/7/fff58fFBTE2qwCkUjEREVF3Y+Ojna9f/++4NixYzcAYPbs2WUbN260ra2t5f3nP/8pZqs+TTMcOHDAPDU1VVxUVGTs4eGhk5U/W6u36SxHQUGBMCwsjJXdRlurJzY2VnL+/HnT3Nxc0dq1a4t0UQcAfPnll7ZnzpyRyOVycu/evRpbW9snfhoiMdSFz/r168dkZmbqOwZKS0vxxVsTsbbfNVQqxfigcBi+3fXzY8elpKQgKChIDwk7j0tZgb/zKpQMhnz+B+5WyXFgfhiC3Ayv08i1tiWEpDIM06XA6enpBf7+/vfZzsSmOXPmuNTV1fF2795d2F11FhcXCwYOHNgvLy/vipmZWYd/YHfs2GG5f/9+qxMnTtxgM0dZhDfUAAAgAElEQVRqaqooJCTE9/r165fd3d27b0lwA8jQXfV2Rz2G8D5q49///rftzz//bJOTk2MyYcKE8oMHDxa0d3x6erqNv7+/W8v7DfYykFhsGOMR7Ozs8Oz0ZThe5Q4LXi3683Oxe/fux47j0gcUl7ICf+dNuH4fd6vkcLMWI7CXpZ5TtY5rbfuki4+PN586dWpFd9a5bds26+HDh1dq0lEBgBkzZjw4f/685ObNm6wNcigpKeEvWrTIdenSpbd37typl169LjIsXrzYafHixU7dXa+u6uno9RjC+9geTd4PZ2fnhqVLl955+eWXtfpiw0pnhRAyhhCSTQi5TghZ1s5xUwghDCGkw7/o+h6z0tyLL72I2OqBiMnywnu70rBgwQLcuPHol6DExEQ9pes8LmUF/s57pegBAGBSgAsIMcwxZVxr2ydVXV0dkUgkA3Nzc02mTp3q6e3t3W1LC//+++8WERER1c3vmz9/vsvIkSN7N92eN2+eS1hYmLdcLidisZjx9fWtiY2NNWejfplMRqZNm+a+YcOGWytXriyJjY21rKys7NYvpvrK0F31dkc93d2G7f2OalPu66+//uDVV199YG1trdV6E1qPWSGE8AFsATASQBGAZEJILMMwmS2OkwB4G4BGIxCVSsMaM/TFuk0ICAiAUqmEVCrF9OnTkZCQoJ5i3Z27RGuLS1mBv/O+NdwL4wc4wVxkuEOtuNa2TyqRSMTExcVljRo1yqesrCy9O+vOzs428fX1fWT8yYoVK+54eXn1P3funElCQoLZmTNnzBMTE7OFQiEDAN7e3nXp6emPDfwcNmyYZ0pKillr9QQFBUnPnHl8/r6JiQkTFxenvv/y5ctZ2r+qztFXhu6qtzvq6e427Oh3VN/Y+KsfAuA6wzB5AEAI+QXACwBaDjhZDeBLAO3PATZQEokEv/zyC8LCwmBnaYp3gpX4eNmH+Grd1/qO9lRxtzHVdwSqDW7L/i+wrcc+Gtfn5tyI3vcBYHv8DZt/H8vq1daxBZ+PT236eeS6s31zS6XilvdrIjk5WdynT59aqVRKhg8f7g0AGRkZYl9f31oAOHLkyHV7e/sOByYePXpUEhsba7F9+3aNBkxWV1fzW87QcHBwUERHR5dGRUW5S6VSfnx8fJa1tbW6bolEoiwpKXnsMlBrnRGK0oWOfkf1jY3OijOAW81uFwEIbX4AIWQQgJ4MwxwlhLTZWSGEzAUwFwCcnJwQFxcHAPDw8IBEIkF6uuoLkrW1NXx9fREfH696EQIBwsPDkZaWpl68LSgoCHfv3sWtW6poXl5eEAqFuHr1KgDVWBRvb28kJCQAAIRCIcLCwpCSkgKpVAoACA0NRVFREW7fvg0A8PHxwccff4zwqkN4TpILyyol9u79BY6ODgBU01ZDQ0ORlJQEmUw14DwsLAz5+fkoKSkBoNroTqFQqJdld3Z2houLi3rKq5mZGYKCgpCYmKj+lh4eHo6cnByUlqomFvj5+UEulyM3NxcA0LNnT9jb2yMlJQUAYG5ujoCAACQkJKCxUXXmLSIiAhkZGSgrKwMAVFRUoLq6Gnl5eQAANzc3WFlZIS0tDQBgaWkJf39/nD17FgzDgBCCyMhIpKeno6JCNQQgICAA5eXlKCgo0On7FBQcgh2//oHeFjyIRCKN3ic+n4+mAdoODg5wd3dXX6IxMTHR2ftkY2Oj/r1l433y9/fX+fv0pLp06ZLYz89PZmZmxly4cCEbAPz8/Po2/cw2hUIBPp8Pc3NzRVVV1WOn6wMCAmrXr1/vuG3btnxPT89HBkpWV1fzLCwstDpNTghps7PYGQzDqDuFISEhPsnJya2e2QkICJCmpqY+0pa6yND87FJ9fT0PAL799ls74PGzS12pv6kuXb7W7no92rx3D+9v83e0K/nZpPVsIELIywBGMwwz5+HtVwGEMAyz6OFtHoDTAKIYhikghMQB+IBhmJT2yvXz82OaPrAMiVKpxOyZk7DB6zwseDL863o/zPrXj2hoaODMjruZmZmcyQoAO48nY3VcKcb6OWDbLFb+FuoM19r2SZ4NFB4e7jVjxoyyf/zjH+VN9/n5+fW9evXqNQC4deuWYMqUKR4KhYLY2dk1xMbG5kVHR/ecOXNm+YgRI2oOHz4sOXHihPm4ceOqYmNjLbZu3Vo0efJk9+LiYmNTU1PF/v3785OSksTr1q2z5/F4zPPPP//g3XffLRsyZIj3q6++en/BggXqei9cuGDy/PPPez377LNVd+7cMT537lxO86xDhgzxmj59evmiRYvKmt8fERHh1d5loPj4+Nym2ydOnDBbvny5MwCUlpYaPffcc5U7d+68derUKdMPPvigp5GREePg4FB/4MCBAqFQyCxatMj53LlzZra2tg379+8vkEgkWl97Lysr4w8dOtT7xo0borNnz14LDg6uY7P+psGc69ata3XadVt1AcA333xjtWzZsp4VFRWs99KPHj0qWbNmjaNSqcTChQtLX3vttQeaPK+j19OkvdfFlo5+R1ujaX4AePvtt51u375trM/ZQEUAeja77QKgeXAJAD8AcYSQAgCDAcR2NMi2ocEwZ2jxeDys2/o9Pkh3BwAs752F1R/Mxs2b3bJ8AyuaztBwxaErqs9DQ50B1BzX2vZJlpWVJQ4KCmpzTQ1bW1tFQkJCTkpKSraDg0PD0aNHJbNnzy7btWuXNQDs3r3b+s0331R3Hn766SdLZ2fn+uTk5OyXX365/IsvvrADgMrKSv6JEyduvPvuu2UAMGrUqMr4+Hj1gmv5+flGL730kuf69etv7tq1qzA7O9vk6NGj6sdlMhnJyMgwnTBhwmN7esTHx+fW1tZebO1f844KAIwePVp64cKF7AsXLmQHBgZKJ02aVAEAnp6e9QkJCdnJycnZbm5u9T///HOP5ORkUX5+vjA1NTV7+PDh1Zs2bWJl+XszMzPliRMncseOHauegdWd9bdWF6A663Xo0CFLR0fHejbqaa62tpasW7fO/syZM7lJSUk5mnZUOqOt18WWjn5HtdHQ0IDa2lqiUCiIQqEgtbW1pCuf72x0VpIBeBFC3AkhxgCmAYhtepBhmEqGYWwYhnFjGMYNwHkAEzs6s2LILC0tMT9mM74u9IaAKLHe/wa2b1kPQ12zhssyiiuRUaaE2JiPKYF0KXtKM4WFhYKqqiq+v79/mwutlZaW8seNG9c7ODjY548//rC4deuWcXh4eG1ubq6orKyMf/PmTWHzhdpyc3OFISEhNQAwZMiQ2uvXr4sAYNCgQTV8/t8LFM6bN6/s9OnTFlKplJSXl/PGjBnjtWDBgrszZ86slEgkyoULF5bExMSoN5Hau3dvj9DQ0Go3NzdWvqHJ5XJy6dIl0zFjxkgBwM3NraFpGrVAIGB4PB5z+vRpyejRoysB4MUXX6w8d+4cKx9MQqGQcXJyeuRyVnfW31pdgOqsyqRJkyp4PPYn05w6dcpMJBIpR4wY4Tly5MjehYWFrM8AaOt1sUGT31FtLF261MnU1DRg69atDocPH7YyNTUNWLp0abvTnVuj9TvHMEwjgLcAnABwDcA+hmEyCCGrCCETu1quiQkrKyLrTGBgICTPzMcZaS9Y86RYOeAW1ny6Wt+xNOLn56fvCBrb8Wc+AOCVoJ7oIdbf5paa4lLbPslcXV0b6+vr09o7Vb5z507rESNGVCUnJ2cPHTq0sunLxujRoyujoqJcx48f/8j6LF5eXvKkpCRTADh37pzY09OzDlCdbW3O0dGx8eWXXy5bt26drZWVlTI7Ozvzk08+UZ9yW7Vq1d20tDT1zI4NGzbYr1mz5jYbrxsAfv31V/Pw8PCq5h0oAMjKyjI+deqUxdSpUysrKir4PXr0UACAlZWVoqKiQufLQbNR/7p164o1ueTQvK7GxkYcOHDAas6cOeUdPa8r7ty5Y1RQUCA8derU9Tlz5txftmyZxh/Emr6eJs1fV9fSPk6T39G2aJJ/3bp1xQzDpDb/15nX3ISVHiDDMMcAHGtxX0wbxw7VsEztg+nYm7Nn4/130uDWKEVsoRmW7/kX+vbzxeTJk/UdrV1cmV5b/ECGI+nF4BNgdri7vuNohCttSwGjR4+uioqKcj958qS5iYmJ+g/OnDlzyry8vJy3bt3afOIAZs2aVXHo0KEeQUFBPmKxWLl///685OTkVlev3Lx5s8adD7anpO7fv99y9uzZj4wlKi8v582aNct99+7deSKRiLG0tFQ8ePCA//AxvqWlpU5nfXRn/S3r2rRpk/WUKVPKW3be2GJpadkYHBwsFYlEzMSJE6u++uorB13U0/J16aIOQ2awK9jW1bG2TYZOff6fDViQ5I7le1RXtV599VX1bA9D1TQ7xdB9l5CPRiWDIHu+Qe6w3BqutO3TqmlwLQAMGTJElpOTk/nHH3/cOHr0aN7bb7+tHp/y7LPPVvbs2bMRAJ5//vnq7du3FxkZGeHIkSP5KSkp2fHx8bn29vaKpsf08VpaI5fLSXp6uumoUaOkTfc1NDRg0qRJHjExMcX+/v5yABg2bJj05MmT5gAQGxtrHhYWJm2rTG11Z/2t1ZWZmSnas2eP9bPPPut18+ZNYVRUVM+OyumMiIiImtzcXBOlUonExESxLnZ2bu11PW0MtrPCFUZGRvj5v/vh4qIaT2HXQ4zNKxepp4lSXWdnLoREJMBYd91ttU5Rzf3xxx+mEydO9FyyZEmJvrN0xeHDhyXPPPPMI5eAtm/fbpWenm766aefOoWEhPh8++23liEhITJXV9f6wMBAn99//91i0aJFrM3qioyM9Pzzzz/No6Oj3TZu3GjdnfW3Vte2bdtu//XXX7l//vlnbq9eveS7du261XFJmnNwcFBMmDChIiQkxGfZsmUuq1atYn2DyNZeF9t1GDqD3ciwf//+zJUrV/QdQ2OnTp3CwugoxM8isOTVYO7lAVi385BB7tB848YN9O7du+MDDUCNvBElRTc5k5dLbQs82VOXKYriHs5tZNi0pTdXDB48GN/+8DP2F9vBmCiwacA1LIqegdraWn1He4y9vb2+I2jMVCjgVF4uZaUoiuIKg+2sGNJGhppISUlBREQEbEd+iENl7pCQOqz3zcDcN2aivp71qf1aMfQxNVvjruP7v/JR16Aac2foeZvjUlaKoiiuMNjOCldNnTYNhe6zcLa6J+x4VfjM+xIWRL8BhcJgtlgwaCWVdfj6VC5WHsnE9VKdjfmjKIqiOMRgOyu6mmamK+bmf+/u/u7iD3DG5HmkyJzQk1+OZT0vIHruHIPZSbp5VkPz5YksyBuVGOvnAD9nCwCGnbclLmWlKIriCoPtrIjF3Jiq2iQgIOCR2/9a/Rn2yoYiQ26H73LM8f13uzBv3jyD6LC0zGooLhc9wKG02zDm87BsbB/1/YaatzVcykpRFMUVBttZadpRlyuadm9uQgjBVxu2YENpOD7fr9ohd8eOHZhrAB2WllkNAcMwWHVEtUvyG8+4oZe1qfoxQ8zbFi5lpSiK4gqD7awY6pTqtjQ2Pr67O4/Hwzff7sTrr78OAOjjZo/XLJMxd+7sVo/vLvqsuy2x6cVIuVkBa1NjLBzu+chjhpi3LVzKSlEUxRWsb7hEPYrH42Hnzp1gGAavWqUgQpyPnkIp5rz5Or759jsIhUJ9RzQIR9LvAACWjPGBuYhb09YpiqIo3TLYMysSCSubcHabiIiINh/j8/n4/vvvcUQWgisyW7jz7+Ezt78Q9dp0VFdXd2NKlfay6ss3rwZi68wAvBL0+ErYhpi3LVzKSlEUxRUG21mRyWT6jtApGRkZ7T7O4/GwYdt32FM/GuelDnDkVWBrn/OYN3sWSkq6d2XvjrLqA59HMK6/Iwghjz1miHnbwqWsFEVRXGGwnRWuXfsvKyvr8BhCCD5btxlxFtNxoswRlqQGO/ueQ8ziObh27VqHz2eLJlm7Q3lNPRbtvYiSyvY3rTSUvJrgUlaKWrRokXNgYKDPmDFjPKqrqx/7PGj5+K1btwSDBg3qExwc7DN48GDvmzdvtnnN9uDBg+bLly+3nzZtWi+5XP74t5Au5morQ2eyaaKsrIzfv3//vmKxeFBycrKoszkpdtHG7WaEECxbvgKFfu/ix1uOMCH1MKq5jSFDhuDUqVP6jtdtGIbBkgPpOJJejE9+varvOBT11ElOThbl5+cLU1NTs4cPH169adMm644ed3R0bExJSclKTk7OnjFjRtmWLVts2ip/8uTJVatXr75ramqq7ExnpaNcbWXoTLaWQkJCfFreZ2Zmpjxx4kTu2LFjK7qSk2KXwXZWuLbOir+/f6eOj543HzaT1mL8/1lia+wlPHjwAGPGjMHGzVt0PhOqs1l14YdzBTh1rRTmIgFWvuDb7rGGkFdTXMpKPd1Onz4tGT16dCUAvPjii5Xnzp2TdPS4QCBQL9hZXV3N9/PzkwFAXV0dGTBgQB+pVEr27NljMWPGDFcA+OKLL2zHjh1baW5urvF6DR3laitDa/e3lUsTQqGQcXJyavMUf0c5KXYZbGeFa8vTd2Wg7Ljx47H6m0NwcnICAPj3todv0U94/a33UFfX/qURbehjUG9ziTfK8On/qS57fT55AJx7mLR7vL7zdgaXsj4NiouLBUOGDPEyNzcf+Morr/TSRR0LFy50XrVqlZ0mx/bv379vSkpKq5cUultFRQW/R48eCgCwsrJSVFRU8DV5/Ny5cyYDBgzo8+2339oNHjy4FgBEIhETFRV1Pzo62nXXrl02P/zwQ+GXX35pe+bMGUlmZqbo3r17Gi9J3lGutjK0dn9rubrSVl3NSbHHYDsrcrlc3xE6JS8vr0vPCwgIQHJyMoKDg7B1vAjPia7hC5vDeH3+AtzIK2A35ENdzcqGW+W1+MeeVDQqGcyN8MC4/o4dPkefeTuLS1mfBjExMQ4eHh7yqqqqS/v27bvJdvnFxcWC/fv3Wy9evPieJse/8847JR9//LEz2zm6wtLSUvHgwQM+AJSXl/MtLS0Vmjw+ZMgQ2eXLl7M++eST2ytWrHBoOj4yMlK6b98+m+3btxcaGRlhyZIl944fP54XExNTamtrqy67sLBQEBgY6NPyX2FhoUCTXO1laO3+lrma5ObmGoeEhPiEhIT4XLt2zaTp59zcXGM22o9il8F2Vp4mTk5OOHs2Hj+WB+BcRQ84knL81Oswdm34CDv/G6vveKyR1Ssw+4dkVNQ2YKiPLZaO6dPxkyhKC/Hx8eZTp05tdcwBG7Zt22Y9fPjwSjMzM42u3c6YMePB+fPnJdoO/mTDsGHDpCdPnjQHgNjYWPOwsDBpR4/LZDL12BNLS0uFWCxWAkBJSQl/0aJFrkuXLr29c+dOq/bqdXV1bUxNTc1u+c/V1bVRk1xtZWjt/vZyeXl51V+4cCH7woUL2X379pU1/ezl5VXPRvtR7DLYzgrXFktzc3PT6vkmJibY/P0+XPJcjK2ZEhgTBVZb/R8cL6/HjPc/hUzG3mUhbbN2lYkxHy8OcoaXnRk2Th8EPk+zMXf6ytsVXMr6JKurqyMSiWRgbm6uydSpUz29vb376aKe33//3SIiIkJ97W/+/PkuI0eO7N10e968eS5hYWHeTQNMxWIx4+vrWxMbG6v3HS9DQkJkrq6u9YGBgT6///67xaJFi+4XFhYK3nvvPae2Hv/rr7/EQUFBPqGhod5ff/21/fLly+/KZDIybdo09w0bNtxauXJlSWxsrGVlZWWXP1taqxdQnZF57733nFrLAAAt71+yZEmptrkiIyM9//zzT/Po6Gi3jRs3WjfP0VZOSjcIG4M5CSFjAHwNgA9gB8Mwn7d4fDGAOQAaAdwD8CbDMO2ekh00aBBz8eJFrbN1l6qqKtZ23E1MTMQPH03Dl+E1MOc3IFvhjFdOWmPv99+gXz/t/+aymbUrZPUKmBhrfnlX33k7g0tZAYAQksowTFBXnpuenl7g7+9vsH+gU1NTRaNGjfIpKytL11UdlpaW/r/++mtuZGRkLaA6w+Dl5dX/xIkT2QkJCWY7duywTUxMzLa2tlZfIoiKiuopEAiYHTt2FDUva9iwYZ4pKSlmrdUTFBQkPXPmzHVdvQ6KMhTp6ek2/v7+bi3v13q5fUIIH8AWACMBFAFIJoTEMgyT2eywiwCCGIapJYQsAPAlgKntlVtbW9vewwYnLS0NQ4cOZaWssLAweB9Iw/tzpyLa8Sr+X74Qly8kICAgAKs/XYOFi96GWNj1s8hsZu3Ig9p6fHjgMj4Z31e9OWFnOipA9+bVFpeysm6FRWCbj41cfRPPvK3q2Py10Qa/L297sOuKylT1z1tC+uJetvix+zWQnJws7tOnTy0AZGdnGw8ePLivl5dXnUwm423evPlmUwejneeLPv/8c4eDBw8WtHVMdXU1v/lMFwcHB0V0dHRpVFSUu1Qq5cfHx2c176gAgEQiUZaUlDz2H5h2RiiqbWzsDRQC4DrDMHkAQAj5BcALANSdFYZhzjQ7/jyAWSzU+0SztrbG9gO/Y9uWzdi7cwkA1aDjnPhfML24Di+OGYOokQGtrvhqKG6V1+LNXcnILZWiRt6In6MH6zsS9RS5dOmSuGlaKwCEhoZWHz9+PO/06dOmH330kfNff/2Vq20d5ubmiqqqqkcuLQQEBNSuX7/ecdu2bfmenp4NLZ9TXV3Ns7Cw0GrVS0JI2x1D6onAMEynOudPOjY6K84AbjW7XQQgtJ3jZwP4rbUHCCFzAcwFAHt7e8TFxQEAPDw8IJFIkJ6uOptrbW0NX19fxMfHA1DNrw8PD0daWhqqqqoAAEFBQbh79y5u3VJF8/LyglAoxNWrqgXI7Ozs4O3tjYSEBACqMTJhYWFISUmBVKoaJxUaGoqioiLcvn0bAODj4wM+n4/MTFU/zMHBAe7u7khMTIRMJkNSUhJCQ0ORlJSk3i4gLCwM+fn56iX1+/XrB4VCgezsbFXjOTvDxcUFSUlJAAAzMzMEBQUhMTERcrkc/fz64/z585gxYwYUd7OwaeANCMh/sCP+GsKOR+CNIb0wYpAn7O3tkZKSAgAwNzdHQEAAEhIS1CsBR0REICMjA2VlZZDJZKioqEB1dbV69oqbmxusrKyQlpYGALC0tIS/vz/Onj0LhmFACEFkZCTS09NRUaEarxgQEIDy8nIUFBQ89j5dvd+Iby43oLpeCWczgpd71iIhIaFL75OFhYX6d0Hb9wlQjQ9i+30CgPDwcABQZ/Xz84NcLkduruozsWfPnp16nwDVui26fJ9YpemZj2fevq8+y9KRhRe6vLTz1atXTWbMmPHYksJhYWG1d+7cMb5165ZgypQpHgqFgtjZ2TXExsbmMQyDF154wePBgwcCDw+POgBQKpV44403emZlZZnw+Xz89NNP+b17924AgD59+sgyMzNFTWdpLly4YLJ48WLXSZMmlf3444828+fPL29Zf05Ojmj69OmP3R8REeHV3mWg+Ph4deeKfpBRTxutx6wQQl4GMJphmDkPb78KIIRhmEWtHDsLwFsAIhmGaXduclBQENP0R50CGhoa8OWaFTC9sAGLggXgE6BAaY9/NUahyiYAK18ZjAEuPfQdE/JGBb4+lYttZ2+AYYDhfeyw/pWBsBDrffID1YonecyKjY2N/7Fjx3JCQkJk2dnZxu+8847L8ePH8w4dOmS+bds22yNHjuTx+XzGyMgIb775Zs+JEyc+qKioEKSmpoo3b958+8svv7RNSkoynTJlSkViYqLpxo0bi0+fPm26a9cu6x9//LEQAFasWGGfnZ0t2rt37838/Hyj8PDwPl9++WXhxIkTq93c3Pr/8MMPec8//7x6AK5MJiMODg7+6enpGW5ubo+ddaGop11bY1bYmA1UBKD5VrkuAIpbHkQIGQHgYwATO+qoANxbXOvs2bM6Ld/IyAgfr1iD4Z8n4s2/euJyKQM33l38YPwF5lSsxYLNv+JIUpZeszIMg2nbz2Nr3A0AwDvPeWHHa0Fad1R03bZs4lLWJ1lhYaGgqqqK7+/vr55Gl5SUJAkJCfHZtGmT3fr164tKS0v548aN6x0cHOzzxx9/WNy6dcs4NzdXGBQUVAMAYWFhNQCQkZEhOn78eI+QkBCfpUuXulRWVqoHXc2bN6/s9OnTFmVlZfwxY8Z4LViw4O7MmTMrJRKJcuHChSUxMTGPrKmyd+/eHqGhodVsdFR0ua/PlStXhGvWrLGLiorqWVJS0qlBZu3lOnXqlOnAgQP7BAcH+0yYMMG9vWX4dbW3EAAcPXpUEhYW5h0aGur9448/9gCAEydOmDWtteLm5uY3e/bsx7eA7wS6txC72GicZABehBB3QogxgGkAHlkchBAyCMA3UHVUSlmo0+Doeon8JgMGDMDO45cQ5/UvfHJWiZoGYDz/ApYrtmLa8CCsXLkSNTU1OH61BGXS1vuEbGdtKo8QgimBLvCwMcWB+WF4b6Q3eBpOT9akfC7gUtYnmaura2N9fX2aUChUvyGhoaHVFy5cyD5z5sx1Pz8/+c6dO61HjBhRlZycnD106NBKhmHg5eUlT0tLEwPA+fPnxQDQt2/fuhdeeKHiwoUL2cnJydn79u0raCrT0dGx8eWXXy7bsmWLdXZ2duYnn3yi/vu2atWqu2lpaY98g9iwYYP9mjVrbmv7+nS9r0///v3ldnZ2DUVFRcbNF1JrrrX9dDrK5enpWZ+QkJCdnJyc7ebmVv/zzz+3eTpYk72FupKhtraWrFu3zv7MmTO5SUlJOa+99toDABg9erS0aa2VwMBA6aRJkzRen4fuLaR7WndWGIZphOrSzgkA1wDsYxgmgxCyihAy8eFhXwEwA7CfEHKJEPLkrHT2UHcOdBUIBHj7vfexYHcWltwcij2XG/DerlTU1tZgxYoVCA0ehIW7L2DwZ39gzg/J+N/FIlTV/f1Fjo2sSiWDxBtlmPtjCtb9nqO+f0aIK4698ywCezl9FM4AAB8bSURBVLW7LlSnGPIg4pa4lPVpN3r06KqdO3faPvfcc73LysqMAGDWrFkVOTk5orCwMO9Lly6JAWD69OmVZWVlgtDQUO/Bgwd7b9269ZEPlc2bN9+OiYnR6EvY5cuXs4KDg7VeNKk79vWJjo6uiI6Ovn/jxg2NVnTVJJebm1tD0wJ6AoGA4fF4THsZdLG30KlTp8xEIpFyxIgRniNHjuzdtHJuE7lcTi5dumQ6ZswYKd1byHCwMcAWDMMcA3CsxX0xzX4e0dkyJRJuvW+RkZHdXqezszO27D6Ms2fPwvjG28D9ywCAzUPuwx7vYr3yVZy8FoxT10rB5xEEulriWS8bvDZkSJfqk9UrkHKzHH/m3sfR9GIUV6r+5l669QDvPOcFAZ8HQghERuxukaGPtu0qLmV9mvj4+NQfP378kb0QhgwZIsvJyclseWzL4wDgu+++u9XyPn2qqKjgOzk5NQBt7+vT2uPnzp0zmT9/fq+qqirByZMnc4BH9/W5f/++4NixYzdiY2Ml58+fN83NzRWtXbu2qGX9Xc3VJCsry/jUqVMWn3322Z22MjTtLSSXy8m9e/dqmi/Zr02GO3fuGBUUFAjT0tKyDh8+bL5s2TKnn3/+Wb1n0K+//moeHh5exefzwefzH8ulaVtom5N6FCudFV1omqXBFenp6XrbcTcyMhJpaWn44YcfsHHNR3Azq4GbqAzbsQGpFWbYqngRZ8zG4EJBOVILKxBsIcXgoEEAgLUns1Fd1wgHCxEsTIwgNuaDEIKGRiX8nC3g46DqNB5MLcKyQ5fRoPj7ModzDxNMDnDGrLBeEPB1d7lVn23bWVzKSnGXtvv67Nixw3LFihUOTR/SkZGR0g8//LDX9evXLxsZGWHixInVEydOfGzgYG5urvHMmTPdAaBpPx0A2LNnT76Xl1e9JvvllJeX82bNmuW+e/fuPJFIpP6D0jLDkiVL7i1ZsuSxPZe0zWBpadkYHBwsFYlEzMSJE6u++uorh+aP79+/33L27NnqgeMtc2mao2XulujeQp1jsAN6mqZyckXTNFF94fP5ePPNN/HX5Rv4r+2HWHIGuCtVItBSip02u7GvYibG3fgCwywrUV3x9wSOvRcKsetcAT7/LQv/PHQF7/xyCW/vvYj396fj2JU76uPszUVoVDLo72yBeREeODA/DH8uGYbFo3xgJ9HtJrL6btvO4FJWirv0ta9PR/vpdJSroaEBkyZN8oiJiSn29/dXD6rrzgwRERE1ubm5JkqlEomJiWI3Nzd1DrlcTtLT001HjRol7SgX3VuoexnsmRWqa0xNTbH04xiUL3gLG9d+DnncNrwbBAQ58uBtdQm9YmbgsJkNPvjgA8yZMwcrJ/rhTqUMdyrrIK1rRE29qpNozOfBw9ZUXW6IuxWurhgNUyH9laEofWu+L42NjU3j/v378wsLCwXr16+3W79+fXFrj//111/iZcuWufD5fEYoFDI//fRTQfN9fQYNGiQLDAzs8+GHH96zsLDQeIxIR7kA1eys9evX23l6etalp6ebfvrpp06ffvopoqOjS2fNmvWgOzOsX7++eMKECRUhISE+PB4Pu3btym967uHDhyXPPPNMFZ/PBxttExkZ6Xnt2jXxjRs3RG+++ea9t99+u6x5jtZyUq1jZW8gXXia9wZiU1lZGTav/xLSP/8fFPV1WH9e1ekXCYB/RppC6j0Zr/3jQ/j5+ek5adsMtW1bw6WswJO9zgpFUdyjy3VWdEKh4Nblu/LyxxakNAjW1tb416df4JPYQti/uAqOjo4AgBn9jRATzsca6/8hY0UIFk0YiN0//WSQY4UMtW1bw6WsFEVRXGGwnZWmJcy5omkpc0NlYWGBpUuXIj8/H0uXLkWdeW8cyGwAjwBT/YywKTAfwUkL8NnzDvjnW2/iwoULBrNmiKG3bXNcykpRFMUVBttZoXRDKBRizJgx2H0mE7YLf8PC60Ox+s8GFFcr4WPDx6pngfG1/0VoaCj69euH1atX48YN1mbrURRFUVSnGexoSaFQqO8IneLh4aHvCBrz8PBQb3gXGRmJO3fuYNd3O3HjxDcY71iGo7mqQbZZWVk4uGUFPC5+hq31veAc8RpenDK1218r19r2KcI0bZ5IURSlLaVSSQC0OoDZYM+sNK20yBVcWsSuZVZHR0f88+NPsD3uJiwXHIPSfyZMTVUzgWYNMMLMAUZYG1SMORWfIelDP3w4zhOfrfgnUlNTu+VSEZfb9klGCKmsr6+nO1RSFMUKmUwmIoSUtPaYwXZWamtr9R2hU9LT0/UdQWNtZeXxeBg6dCi+//573L17F3v37kWWJByfnGlA2h0FzIUE0/sb4auQe3hfsQUZn4bDxcUFs2fPxoEDB3S2xsiT0LZPIoVC8X1xcbHpw29DFEVRXaJUKklNTY1JQUGBcWNj48rWjjHYy0CUfpmammLatGmYNm0aKioqEBsb+//bu//YuO/7vuPPN4/k8bdE0vwhk7RERSRnig0t6iKaMS1pS9wmbjCnm+t4QzNlSGBsnYdmw5CkKDYsG7A5CVq0QIN1hevWrou6W7POXtw0sRNTMgGGNkWLtihZpEOxESmKMn9Y5FEUfxw/+4NHWlEk6pgT9f1+7NcDOOhO99XpKZ4gvvX9fu775dvPP8Md73bxmT3Gx2sizC/D+fPneeqpp/juX/0pXz+cw3CklsJf+jQf/8Sv0tbWRk7O1p4wToKzsrLyP2ZmZlreeuutdsCvXaEiEiYrZnZheXn56y0tLd+/3gahHVYyM0Obdl2lpf5cMHOzrcXFxRw5coQjR44wOzvL9773PR5/8a/p7HsZWL0+0Cd3R/hXsSxghBX3x7zxl3/E//yW42LuHgoaH6Cl/QHuvfdetm3btuW9QfKpNV379+9fBI4E3SEiH3yhPSlcLBZzPT09QWekbGVlhYyM0B5V+xm3qjWRSPDaa6/x4osvcrrzuzS6t/lEbYS26gjRzPePDCwmHNuemGUhYezdu5fPtH+Uuv2H+Vjrvdx99903HUw/jF/b2yWdk8KJiNwuoR1WGhoa3JkzZ4LOSFlHRweHDx8OOiMlW9U6OTnJK6+8wqs/+gHvnXyJ3ZExDu1cPTrwD59eXYNkwMRXCsmOQO9Ygr53janoXUSqW6j66CE+ek8Le/fu/ZnDR/rabh0NKyLiA7+OtUiolZaW8vDDD/Pwww8DcOHCBY4dO0Zn56vcc08nb775JmW5jktXHLXFGRzcmcnBnQCjwCjzQy/wG787z/MDjvr6etrvqePuhnpmXD7l5eXs2bOH7OzsIP+IIiISgNAOK76du8GnNTa3q7WyspJHHnmERx55BIB4PM7rr7/OX3V3c7rnVVZGj7Mza5r9OyI0V0bYXZzByIwjkUhw+vRpvlA1xL+r7+DSFcfp3/sWz005JlwxS4V3kVF5N0UfaWXPnj3s2bOH6urqUHzc3ae/ByIivgjtYSDf1qzIL2Z8fJzjx4/zxhtvcKbvNV4/0c+Zd4ZwzvGNT0b54r4sSvN+fg3Ia6MJWp+cAyDD4Jlfy+M9iljIrcC230Xujnq237WX6to6ampquPPOO8nK0ilBrqXDQCLig9AOK42Nje7UqVNBZ6Sst7eXlpaWoDNSEvbWeDzOqVOneOuttzh1qp+3j3eSd3mU7Yl3aSjNoK4kgzOTK3z15dXrR921zfj7L1//ZGwTl1f459+Z5+WzK1RUVPBgUwn7avLIKNpBduld5JfXsu3OPVTcWU1lZSXl5eVpDTVh/9peS8OKiPggtPusfbvq8szMTNAJKQt7a0FBAQcOHODAgQPA+4tW5+bmGBgYYGBggIWBAf7Fne8wODjIxMhP+JfPT7Nrewa7thm7tmdQsy2D6iLjjrwMZhbAOceFCxf4peZpHr8rCgxAAhhbvU2+usLxsRWqn71MSUkJZWVl/IePrZCRWwS5xWQWlpNVVE5O8Q7y76iisKyG7aVllJSUUFxcvL6WJuxfWxERH4V2WBG5Vn5+Pvv27WPfvn0/91w8Hmd4eJihoSH6hof5fz/9KT99e5j4+DBjbgSzizjn+PFIgj99Y5HKggwqC4zKAqM83yjNy6AounpJiqmpKS5NT/GlR4uA8fd/k8vJ2yj86xfn+aOeJQB+7R9k8tX7c7mcyGRuOYMXns0nEcllJSuPlawCuiMfo7CwiKKiInZHp8nPzSFaWEy0oIScohLyikrJ3VZCfuF28vPzdbhKROQat2RYMbNPAX/A6lksn3TOPXHN81HgGWA/MAl8zjk3vNFrrl2bxhexmD970n1qhdR6CwoKaGpqoqmp6brPLy4uMjY2xujoKKOjowyOjXH0/HnGxsYYPznG/OQos9PvYjaPc47MDPiPr1zhjjyjNNcozc2gJNcozoXiHGPy8vuHT2uLM2i901jdVZMA3kveYH7J8ch/e3V92/7fzKex7PoLgf+ge4Ev/90CWVlZfLw2jyc/DUsrERaJsOwiJCyThGXiLJPnpj/KfHYJOTk53Js/Qk32JSwzG8uMYpnZZGRGyciKsphTxsQdrUSjUaJZEXZc6iUzO4dIdg4Z0YJU3wIRkUClPayYWQT4NvAAMAK8bmYvOOeuXnDyRWDaObfHzB4FvgF8bqPXXVpaSjftthofH6egwI9//H1qhVvTm52dzc6dO9m5c+eG2yUSCaamphgfH2diYoJ3332XiYkJzk5MMDExweTkJJOTk0wXTlNXN8X09DTP9b9H17k5tuUYRVEoihpFUaMw28i45kNtb4ytMLPgKMg2cjONvCzIz1798crqxa5ZWloishRnz/Z8YDl5+1lH/vz/8s7U6p6g9n+ay8GmrNVrlS4mb0kdw8v8k6e/AUBhNsz8dtH6cyMz1724qYhI6NyKPSsHgHecc0MAZvYc8BBw9bDyEPCfk/f/GvhDMzO3werexcXFGz0VSufOneMjH/lI0Bkp8akVbm9vJBKhrKyMsrKylH+Nc465uTmmp6d5+eWXqa+v59KlS+u3bz44y8zMDLOzs/xwZobnZ+LMzs4Sj8fXb3NzcS7PGZHIMolEgh+PJGj8dpxoJuRmGtFMyMk0cjIhGoGx2fcHjb94a4kT4wmyMiCaaWRlQHYEsiPG4NT72604+M6pJbIikJkBU/PhXFwvInKtWzGsVAHnrno8ArTeaBvn3LKZXQJKgYmrNzKzx4DHAMrLy+no6ABg9+7dFBYWrl/RtrS0lL1793Ls2LHVP0RmJu3t7fT29q4vcIzFYoyPj3Pu3GpaXV0d0WiUkydPsvb69fX1dHZ2AhCNRmlra6Onp4d4PA5Aa2srIyMjjI6OAtDQ0EAkEmHtU0qVlZXU1tbS1dVFPB6nu7ub1tZWuru7mZ+fB6CtrY2zZ89y4cLqVa8bGxtJJBKsnZ23qqqK6upquru7gdXDGbFYjK6uLhYWVj/t0t7ezsDAABcvXgSgqamJhYUFBgcHAaipqaGiooK1j3oXFRXR0tJCZ2cny8ur/ys/ePAg/f39TE5OEo/HmZ6eZnZ2lqGhIQB27dpFSUkJvb29wOr1gJqbmzl69CjOOcyMQ4cO0dfXt3515ZaWFqamphgeHt7S98k5t/53Id33CSA3N3fL3qcdO3awtLREXl4eBw4c2PT7dP/993PixAnOnz/PlStX2LVrFxMTEwwNDbGwsMC2bdvIyspi99tvs7i4SCQSobS0lLeTj5eWligrK+P8+fPE43EWo4s8+GA+c3NzzMzM8PW3l8nIyCCRSDA7Owv8PSIiYZf2R5fN7NeBX3HOfSn5+PPAAefcv71qm/7kNiPJxz9JbjN5o9dtbm52a9/0fDA6OkpVVVXQGSnxqRX86vWpFfTRZRHxw6244toIUHPV42rg/I22MbNMYBswtdGL+nYG22g0GnRCynxqBb96fWoVEfHFrRhWXgfqzKzWzLKBR4EXrtnmBd6/lPzDwI82Wq8CrO+e98XaYQsf+NQKfvX61Coi4ou016wk16A8Dnyf1Y8uP+Wc6zez/wL0OOdeAP4E+HMze4fVPSqPpvv7ioiIyIfDLTnPinPub4G/vebn/tNV968Av76Z1/TtxFjl5eVBJ6TMp1bwq9enVhERX9yKw0Bbwrdj//X19UEnpMynVvCr16dWERFfhHZYWftYqi/WPgLtA59awa9en1pFRHwR2mFFREREBEI8rGRkhDbtunw6bOVTK/jV61OriIgv0j4p3FaJxWJu7UyfIrI1dFI4EfFBaHdfXL58OeiETfFpsPKpFfzq9alVRMQXoR1WEolE0Amb4tOCYJ9awa9en1pFRHwR2mFFREREBEK8ZqWlpcWtXQHYB/Pz8+Tm5gadkRKfWsGvXp9aQWtWRMQPod2zsrS0FHTCpoyMjASdkDKfWsGvXp9aRUR8EdphZXFxMeiETRkdHQ06IWU+tYJfvT61ioj4IrTDioiIiAiEeFjJyckJOmFTGhoagk5ImU+t4FevT60iIr4I7bBiZkEnbEokEgk6IWU+tYJfvT61ioj4IrTDyvz8fNAJm3Lq1KmgE1LmUyv41etTq4iIL0I7rIiIiIhAiIeVrKysoBM2pbKyMuiElPnUCn71+tQqIuKL0A4rvl29tra2NuiElPnUCn71+tQqIuKL0A4rvl1jpaurK+iElPnUCn71+tQqIuKL0A4rIiIiIpDmsGJmJWb2kpkNJn8svs4295hZl5n1m9mbZva5lMIy/JqjfLoejE+t4FevT60iIr5I60KGZvZNYMo594SZfQ0ods599Zpt6gHnnBs0szuB48Ddzrn3NnrtWCzmenp6fuE2Ebk5XchQRHyQ7u6Lh4Cnk/efBj577QbOuQHn3GDy/nngIlB2sxeem5tLM+326u7uDjohZT61gl+9PrWKiPgiM81fX+GcGwNwzo2ZWflGG5vZASAb+MkNnn8MeAygvLycjo4OAHbv3k1hYSF9fX0AlJaWsnfvXo4dO7b6h8jMpL29nd7eXmZmZgCIxWKMj49z7tw5AOrq6ohGo5w8eZK116+vr6ezsxNY/fRRW1sbPT0964t7W1tbGRkZWb84XUNDA5FIZP3EX5WVldTW1tLV1UU8Hqe7u5vW1la6u7vXT2rX1tbG2bNnuXDhAgCNjY0kEgnOnDkDQFVVFdXV1evf5AoKCojFYnR1dbGwsABAe3s7AwMDXLx4EYCmpiYWFhYYHBwEoKamhoqKCtb2RBUVFdHS0kJnZyfLy8sAHDx4kP7+fiYnJ4nH40xPTzM7O8vQ0BAAu3btoqSkhN7eXgCKi4tpbm7m6NGjOOcwMw4dOkRfXx/T09MAtLS0MDU1xfDw8Ja+T5cvX17/u5Du+wSrh2q26n2anp5eb033fQJobm7e8vdJRCTsbnoYyMxeBq538ojfAZ52zm2/attp59zPrVtJPrcD6ACOOOd+fLOwhoYGt/aNwgcdHR0cPnw46IyU+NQKfvX61Ao6DCQifkh3zcoZ4HByr8oOoMM593NXcjOzIlYHlf/unPvfqbz2/v373fHjx3/httttYWHBm3PD+NQKfvX61AoaVkTED+muWXkBOJK8fwR4/toNzCwb+BvgmVQHFWB917ovzp49G3RCynxqBb96fWoVEfFFusPKE8ADZjYIPJB8jJnFzOzJ5DaPAAeBL5jZieTtnpu98NLSUpppt9faWgcf+NQKfvX61Coi4ou0Ftg65yaBT1zn53uALyXvPws8m87vIyIiIh9eoT3zmm8n12psbAw6IWU+tYJfvT61ioj4IrTDSjoLf4OQSCSCTkiZT63gV69PrSIivgjtsHLlypWgEzbFp49Z+9QKfvX61Coi4ovQDisiIiIiEOJhJTs7O+iETamqqgo6IWU+tYJfvT61ioj4IrTDSlZWVtAJm1JdXR10Qsp8agW/en1qFRHxRWiHFV3IcOv41Ap+9frUKiLii9AOKyIiIiIQ4mElEokEnbApBQUFQSekzKdW8KvXp1YREV+kdSHDrRSLxVxPT0/QGSIfaLqQoYj4ILR7Vnxbs9LV1RV0Qsp8agW/en1qFRHxRWiHlZWVlaATNsWnq0T71Ap+9frUKiLii9AOKyIiIiIQ4jUr+/fvd8ePHw86I2XLy8tkZqZ1EevbxqdW8KvXp1bQmhUR8UNo96z4tjt9YGAg6ISU+dQKfvX61Coi4ovQDitLS0tBJ2zKxYsXg05ImU+t4FevT60iIr4I7bAiIiIiAiEeVnJzc4NO2JSmpqagE1LmUyv41etTq4iIL0I7rIR14e+N+LTGxqdW8KvXp1YREV+Edli5cuVK0AmbMjg4GHRCynxqBb96fWoVEfFFWsOKmZWY2UtmNpj8sXiDbYvMbNTM/jCd31NEREQ+XNLds/I14IfOuTrgh8nHN/JfgaOpvnB2dnaaabdXTU1N0Akp86kV/Or1qVVExBfpDisPAU8n7z8NfPZ6G5nZfqAC+EGqL5yVlZVm2u1VUVERdELKfGoFv3p9ahUR8UW6p9qscM6NATjnxsys/NoNzCwD+F3g88AnNnoxM3sMeAygvLycjo4OAHbv3k1hYSF9fX0AlJaWsnfvXo4dO7b6h8jMpL29nd7eXmZmZgCIxWKMj49z7tw5AOrq6ohGo5w8eZK116+vr6ezsxOAaDRKW1sbPT09xONxAFpbWxkZGWF0dBSAhoYGIpEIp06dAqCyspLa2lq6urqIx+OUlZXR2tpKd3c38/PzALS1tXH27FkuXLgAQGNjI4lEgjNnzgBQVVVFdXU13d3dABQUFBCLxejq6lpfrNne3s7AwMD6OTyamppYWFhYXx9RU1NDRUUFa1epLioqoqWlhc7OTpaXlwE4ePAg/f39TE5OEo/Hue+++5idnWVoaAiAXbt2UVJSQm9vLwDFxcU0Nzdz9OhRnHOYGYcOHaKvr4/p6WkAWlpamJqaYnh4eEvfp/HxcczslrxPsPpJs616nzo7O8nJybkl7xNAc3Pzlr9PIiJhd9PT7ZvZy0DldZ76HeBp59z2q7adds79zLoVM3scyHPOfdPMvgDEnHOP3yysoaHBrX2j8EFHRweHDx8OOiMlPrWCX70+tYJOty8ifrjpnhXn3Cdv9JyZjZvZjuRelR3A9U7f2Qbcb2a/CRQA2WYWd85ttL6FSCRys7RQKSoqCjohZT61gl+9PrWKiPgirQsZmtm3gEnn3BNm9jWgxDn3lQ22/wIp7lmJxWJubXe5iGwN7VkRER+ku8D2CeABMxsEHkg+xsxiZvZkOi+8th7BF2trX3zgUyv41etTq4iIL9JaYOucm+Q6i2adcz3Al67z838G/FmKr51O2m23tkDSBz61gl+9PrWKiPgitGewFREREYE016xsJd/WrKysrJCR4cfs51Mr+NXrUytozYqI+CG0/6qunf/CF/39/UEnpMynVvCr16dWERFfhHZY8e3Y/9pJvHzgUyv41etTq4iIL0I7rIiIiIhAiIeVvLy8oBM2pbm5OeiElPnUCn71+tQqIuKL0A4riUQi6IRNmZ2dDTohZT61gl+9PrWKiPgitMPK2sXhfLF2oTkf+NQKfvX61Coi4ovQDisiIiIiEOLzrJjZLODPZZfhDmAi6IgU+dQKfvX61ArQ4JwrDDpCRGQjaZ1uf4ud8elkVWbW40uvT63gV69PrbDaG3SDiMjN6DCQiIiIhJqGFREREQm1MA8rfxx0wCb51OtTK/jV61Mr+NcrIh9CoV1gKyIiIgLh3rMiIiIiomFFREREwi2Uw4qZfcrMzpjZO2b2taB7NmJmT5nZRTM7GXTLzZhZjZm9YmanzazfzH4r6KaNmFmOmb1mZn3J3q8H3XQzZhYxszfM7LtBt9yMmQ2b2VtmdkIfYRaRMAvdmhUziwADwAPACPA68M+cc6cCDbsBMzsIxIFnnHNNQfdsxMx2ADucc71mVggcBz4b4q+tAfnOubiZZQGdwG85534ccNoNmdm/B2JAkXPuM0H3bMTMhoGYc86nk9iJyIdQGPesHADecc4NOecWgeeAhwJuuiHn3DFgKuiOVDjnxpxzvcn7s8BpoCrYqhtzq+LJh1nJW7im66uYWTXwq8CTQbeIiHyQhHFYqQLOXfV4hBB/Q/WVme0C9gHdwZZsLHlY5QRwEXjJORfm3t8HvgKsBB2SIgf8wMyOm9ljQceIiNxIGIcVu87PhfZ/0z4yswLgO8CXnXMzQfdsxDmXcM7dA1QDB8wslIfazOwzwEXn3PGgWzbhPudcC/Bp4N8kD2mKiIROGIeVEaDmqsfVwPmAWj5wkms/vgP8hXPu/wTdkyrn3HtAB/CpgFNu5D7gHyfXgTwH/CMzezbYpI05584nf7wI/A2rh2BFREInjMPK60CdmdWaWTbwKPBCwE0fCMkFq38CnHbO/V7QPTdjZmVmtj15Pxf4JPB2sFXX55z7bedctXNuF6t/Z3/knPuNgLNuyMzyk4usMbN84JeB0H+iTUQ+nEI3rDjnloHHge+zugD0fznn+oOtujEz+0ugC2gwsxEz+2LQTRu4D/g8q//rP5G8PRh01AZ2AK+Y2ZusDrEvOedC/5FgT1QAnWbWB7wGvOic+7uAm0RErit0H10WERERuVro9qyIiIiIXE3DioiIiISahhUREREJNQ0rIiIiEmoaVkRERCTUNKyIiIhIqGlYERERkVD7/2wqCJ2RuGDdAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.plot(xx,yy,'k',lw=3,label='$f(x)=' + sb.latex(f)+'$')\n",
"plt.plot(xx,yy_Taylor,'--',label=text_Taylor)\n",
"plt.plot(xx,yy_Pade,'--',label=text_Pade)\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.ylim([-0.5,1])\n",
"plt.gca().set_aspect(1/plt.gca().get_data_ratio())\n",
"plt.legend(fontsize=12,bbox_to_anchor=(1.04,0.8), loc=\"upper left\")"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"# SymPyとの比較\n",
"# http://mpmath.org/doc/current/calculus/approximation.html#pade-approximation-pade \n",
"from mpmath import pade,taylor,polyval\n",
"def g(x):\n",
" return sb.exp(-x)\n",
" \n",
"def Pade_Sympy(g,n,m): \n",
" a = taylor(g, 0, n+m)\n",
" p, q = pade(a, n, m)\n",
" R = polyval(p[::-1], x)/polyval(q[::-1], x)\n",
" return R\n",
"\n",
"R_sympy = Pade_Sympy(g,n,m)\n",
"RR_sympy = sb.lambdify(x,R_sympy)\n",
"\n",
"yy_sympy = RR_sympy(xx)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.legend.Legend at 0x11533b6d8>"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAD8CAYAAABZ/vJZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXlclVX++N/nXrYLXBCQRQEFFEigUCSJIrTFbLLNGbe2aRmzmnLax5rKry3Or3WaFptpnaZ92rOmph0NIwpRTCzFxAUVUCAWWeTC+f3BkgvChefCfY6d9+vlS557zz3Pm/sAn/uc5fMRUko0Go1Go+kvFncLaDQajUZNdADRaDQazYDQAUSj0Wg0A0IHEI1Go9EMCB1ANBqNRjMgdADRaDQazYBwSQARQjwnhKgUQqw7zPNCCPGoEGKTEGKtECLNFefVaDQajftw1R3I88DpvTz/GyC+89984B8uOq9Go9Fo3IRLAoiUcgVQ3UuTc4AXZAffAMOEECNccW6NRqPRuAePITpPJLB9v+Oyzsd27d9ICDGfjjsU7D7WiYGBgXj6BeHv54vVaqWxsREADw8PbDYb9fX1Xa/D39+fxsZG2traAPDz86O1tZV9+/YB4OPjgxCCpqYmADw9PfH29qahoQEAi8WCn5/fgPuQUmK1WvHz82Pv3r20t7cD4O/vT0tLC62trQDYbDaklDQ3NwPg5eWFp6cne/fuBcBqteLr6+uSPrq8AOx2O01NTTgcDqSU+Pn50dbWRktLCwDe3t69vscH9wHg6+vbrz4Gep2EEN3PG71O+/cxGNepsbGx+9jodRrIe9zfPjZu3LhHShmKRjMAhiqAiB4eOySHipTyKeApgPSRVvmfeaE8HrKYh6+dO9h+hsnJyWHKlCnu1nAKlVxBLV+VXAGEEFvd7aBRl6EKIGVA9H7HUcDO3l6wUw7n4taF1G+sGlQxV5GamupuBadRyRXU8lXJVaMxylAt410G/L5zNdZxQK2UcldvL6gigDIZxs8+EdTUNQyNpQH2H0IwOyq5glq+KrlqNEZx1TLeV4E8IFEIUSaE+IMQ4kohxJWdTT4ENgObgKeBP/bZaVvHWLTw9OGp95a7QnNQ2bx5s7sVnEYlV1DLVyVXjcYoLhnCklKe18fzEri6P3164mCu9QvOsORzf+GZLLxouiFHjUaj0bgW0+5EH+brzTTLd2Rbv2ecdTsOR5u7lXolJibG3QpOo5IrqOWrkqtGYxTTBpDhw+w80ziZG/ZdySce2bz0v5XuVuqV4OBgdys4jUquoJavSq4ajVFMG0CampoobYvg7fZsavHn5RXF7lbqlcLCQncrOI1KrqCWr0quGo1RTBtAAGakx3Z/XdLo271BS6PRaDTux7QBxMPDgytnnETwvnIWWN7gRl5k5bfm/XQXFBTkbgWnUckV1PJVyVWjMYppA4jNZsPuZ2N8fR43er3DZQHf8vkHb7pb67CotIFMJVdQy1clV43GKKYNIF0bss46YzqLvmxm9ptNvPX2u262OjzLl5t/r0oXKrmCWr4quWo0RhmqVCYDZtq0aVx4oUdn8rgNrF+/nqSkJHdrHUJXMjwVUMkV1PJVyVWjMYpp70C68PX15Ywzzug+fv31191oc3iE6ClfpDlRyRXU8lXJVaMximnvQOx2e/fXM2fNZtvG1VyWMYyVlVW0t7djsZgr9k2ePNndCk6jkiuo5auSq0ZjFHP9Fd6PrloOAKdOO51zfzOZK0aWMDOkhPe/Mt9qrKKiIncrOI1KrqCWr0quGo1RTBtAuorhAIQE2lnZFMtTjuk86TiTp/63yo1mPVNTU+NuBadRyRXU8lXJVaMximmHsA5mUurR/HXjhI6Dhj2mHMbSaDSaXxOm/Qvs6+t7wPG1s05FtnSU8cR/OK9/lu8Gq8OTlpbmbgWnUckV1PJVyVWjMYppA0hXvesuAu1+hLdVcqplFQ97LuXVT79xk1nPVFdXu1vBaVRyBbV8VXLVaIxi2gDS0tJyyGOzJsVyucd/mWFdyVj5E/taHT280j1s2bLF3QpOo5IrqOWrkqtGYxTTBpCeWDDrVJ5tnMy9rXPJ9cjgyXe+cLeSRqPR/GoxbQDx9vY+5DEfby9KZAz/bDubCoJ5eWWJG8x6Ji4uzt0KTqOSK6jlq5KrRmMU0wYQq9Xa4+OXT/0lWd0uaxg1dQ1DpdQr+298NDsquYJaviq5ajRGMW0A6ch9dShzT8vEu34Hc6xfcg9P8Oq7Hw2xWc+otIFMJVdQy1clV43GKC4JIEKI04UQG4QQm4QQt/Tw/CghxJdCiNVCiLVCiDN66scZLBYLU0N+5tb2p7kooJDi9/9hTF6j0Wg0A8JwABFCWIGlwG+AJOA8IcTB6XJvB16XUk4A5gJP9NWvh8fh9zhef/HvuDOnmUvebeKlD1ZQXl4+YH9XERIS4m4Fp1HJFdTyVclVozGKK+5AJgGbpJSbpZT7gNeAcw5qI4GAzq8DgZ19dWqz2Q773JgxY1jlmcG/i1qpa27j5ZdfHpi5C0lOTna3gtOo5Apq+arkqtEYxRWpTCKB7fsdlwEZB7VZDHwihFgA+AGn9tSREGI+MB8gLCyMnJwcoGNli91u7x5fDgkJ4fe//z1fffUVAEufeILrr7+eNWvWUFdXB0B6ejoVFRVs396hFh8fj7e3N+vWraOr/4SEBHJzc4GOVV+ZmZkUFBTQ0NAxMZ+RkUFZWRk7duwAIDExEavVyvr16wGIiIggNjaWvLw8GhoaCA0NJSMjg/z8/O5kkJmZmZSWlnbfJSUlJdHW1saGDRs63rzISKKiosjP79hZ7+/vT3p6Onl5ed17YbKysti4cSOVlZUApKSk0NLSQklJxyq06OhowsPDKSgoACAgIIC0tDRyc3O7c4plZ2dTXFxMVVUVDQ0NnHDCCdTX17N582YAYmJiCA4OprCwI1FlUFAQqampLF++HCklQggmT55MUVFRd76ntLQ0qquru/c+9HSdkpOTWbFiBdBxV5mVlUVhYWG/rlNFRUV3mnSj1wk6PpwM1nX67LPP8PHxccl1go4Kh4N9nTSagSKMFsARQswCpkkp53UeXwRMklIu2K/NDZ3nekgIkQk8C6RIKdsP129iYqLs+uXtifr6elKSx/GHrAjCRo5g2Gl/Ye5pmYa+FyPk5OQwZcoUt52/P6jkCmr5quQKIIRYJaVMd7eHRk1cMYRVBkTvdxzFoUNUfwBeB5BS5gE+wPDeOu2rMI/dbufoGX/ktvhNXOr3NS9/nNtfb5fS25yN2VDJFdTyVclVozGKKwLId0C8ECJWCOFFxyT5soPabANOARBCjKMjgOzurVN/f/8+T3zGiZN4yDGLK1pvoLg9yq17QrKystx27v6ikiuo5auSq0ZjFMMBRErpAK4BPgZ+oGO1VbEQ4i4hxNmdzW4ELhdCFAGvApfIPsbODrcPZH/mn3sST9RmkdM+nnbvAJa88KGh78UIXWPSKqCSK6jlq5KrRmMUl9xvSyk/BD486LFF+329HjihP30enI23Jzw8rBwb0sp3nXkX319fzYP9OYkL6ZoUVgGVXEEtX5VcNRqjmHYnurMsuug0ImQl/+fxb24MXsGn36x1t5JGo9H8KjDtjJ+fn59T7Y4eO4qwpm1cavuYRunNnLcnMfW4YwbZ7lDS09VZyKKSK6jlq5KrRmMU096BtLa2Ot02a1I6d7dewKx9/8falnDq9zYNolnPVFRUDPk5B4pKrqCWr0quGo1RTBtA9u3b53Tb6+ZM45m6TIplDMIWwF3/en8QzXqmayOcCqjkCmr5quSq0RjFtAGkP3h5enBs0C8VDN8rrnKjjUaj0fw6MG0A6UoH4Sx3XXI6o9rLeJT7+T/HoxSuXj1IZj0THx8/pOczgkquoJavSq4ajVFMG0D62ol+MElxUcRXfc10r9XMCNvB8/98dJDMeqanCopmRSVXUMtXJVeNxiimDSBdSe76w2WXzefid5tJfLyBZ196nZ9//nkQzHqmK/mfCqjkCmr5quSq0RjFtAFkIBx//PF8L8ZRVidpbGzkueeec7eSRqPRHLGYNoB4enr2+zVCCBYs6E4CzNMvv0Fzi/OruYwQFhY2JOdxBSq5glq+KrlqNEYxnM59sJg4caJctWpVv1/X1NTE5FNO4eGsemw+XrwXfTt3Xj5jEAwPxOFwKJOJVSVXUMtXJVfQ6dw1xjDtHUhXsaD+YrPZGDXlPJJtVcSLHXxZdPiaIq6kqzCVCqjkCmr5quSq0RjFtAHECHdccjZ/aL6OzJbH2OZ/NK99kuduJY1GozniMG0AsVgGrpaaMJrtTX404AvAA+8PfoptlZZvquQKavmq5KrRGMW0AcTZZIqH4+ZzOoZ1Be0M87Xw9SAPZWVmuq+cbn9RyRXU8lXJVaMximln+5wpKNUbs07J4N63H+aVkGeIFruZ+0IDyx662UV2h1JQUKBMJlaVXEEtX5VcjbJq1aowDw+PZ4AUTPxhVDNg2oF1Dodj3sSJEyt7amDaAOJMQam++P2JSZQXBxMgGrFYYNO2XYwdNcIFdocy0El/d6CSK6jlq5KrUTw8PJ6JiIgYFxoaWmOxWMy5nFMzYNrb28Xu3buTysvLnwHO7qnNEf2pYcHsqSysOZfslr+zxprEjU8OfZZejeYIJiU0NLROB48jE4vFIkNDQ2vpuMPsuc0Q+vQLo3Mg0DERnz1uFK2dN1pr9gZSvqfGcL89kZGRMSj9DgYquYJaviq5ugCLDh5HNp3X97BxwrQBpD8FpXrjrnnnQMMevGhllt8q7vrHSy7p92DKysoGpd/BQCVXUMtXJVeNxiguCSBCiNOFEBuEEJuEELccps1sIcR6IUSxEOKVvvrsT0Gp3vDx9uLMOA+e3HcbD3o+ifdniw1P0PfEjh07XN7nYKGSK6jlq5KrRmMUwwFECGEFlgK/AZKA84QQSQe1iQduBU6QUiYD1xk9b3/42zWzeHtVFesq29i4q45nnnlmKE+v0Wg0RySuuAOZBGySUm6WUu4DXgPOOajN5cBSKWUNgJSyxyVh+9PfglK94eXlReLMv3DMP/byYYmD++67j+bmZpf1D5CYmOjS/gYTlVxBLV+VXH8N7Ny50+P444+PDwgIGD979uzRV199deRdd93lVMbLo48+elxBQYHr/hANgObmZnHMMccc1dDQIF5++eXA888/f5Q7fQ7GFct4I4H9C0GXAQfPJCYACCFWAlZgsZTyfwd3JISYD8wHGDFiBDk5OQDExcVht9spKioCICQkhOTkZFasWNHxTXh4kJWVRWFhIXV1dQCkp6dTUVHRXaP69N+cQWjYEiorK9m5cyeLF9/JPffc3Z27yNvbm8zMTAoKCrqXYmZkZFBWVtY9LJGYmIjVamX9+vUAREREEBsbS15eHg6Hg23btpGRkUF+fn53PZPMzExKS0spLy8HICkpiba2NjZs6NjYGBkZSVRUFPn5+QD4+/uTnp5OXl4eLS0dZXqzsrLYuHEjlZUdcTclJYWWlhZKSkoAiI6OJjw8nIKCAgACAgJIS0sjNzcXh8MBQHZ2NsXFxVRVVeFwOPDx8aG+vp7NmzcDEBMTQ3BwMIWFHbv2g4KCSE1NZfny5UgpEUIwefJkioqKqKnpWIiQlpZGdXU1W7Zscdl1io+Px9vbu7uuRlhYGEFBQd0/C0avE3TkSxus61RWVtb9GqPXCSA1NXXQr9ORzKJFiyLi4uJavv7665KdO3d6jB8/Pmnz5s3fO/Paa6+9tvy2226L/Pjjj38abM/D4ePjIy+55JI9l19++ag9e/Z4fPjhh25z6QnD2XiFELOAaVLKeZ3HFwGTpJQL9mvzAdAKzAaigK+AFCnlYSs+JSYmyq5fRFfx2GOPcdtfFnLTWUlUhE3irjvvIyTQ7pK+c3JymDJlikv6GmxUcgW1fFVyBWPZeIuKirakpqbucbWTK0lISEh66KGHtp911ln1d9xxR3hJSYnPa6+9ttWZ1zY2NooRI0akrl27tnj06NGuWdXTA9u2bfOYMWPGmIMff+edd34aNWqUY9WqVT6TJk1K3rRp09rY2NhB8zgcRUVFw1NTU2N6es4VQ1hlQPR+x1HAzh7avCelbJVSlgIbgCEvHj1v3jwuvuwiFiWUcEvAJ1z/yH+GWkGj0QwBzc3Nwm63jy8pKbHNmTNnbEJCQtKnn34amJ2dXd/V5sorr4yaOnVq9x/uK664IiozMzOhpaVFAPj6+srk5OS9y5YtC3CF00MPPTR8zJgxyXa7fXx2dnb8jh07PAA6g8SGg/+NGjXKUV5ebl2wYMGohQsX7nj22WeDXeHhSlwxhPUdEC+EiAV2AHOB8w9q8y5wHvC8EGI4HUNam3vrdCAFpfrCZrPROvJ43mrYwZtt2Xxd60v5nhoihgcZ7jsiIsIFhkODSq6glq9Krq5ECDFxqM4lpeyzUJCPj4/Mycn58bTTTkusqqoqAggKCkpNTk7unvxcvHjxrvj4+KO//vprW25urv+XX34ZkJeXt8Hb27t7WCYhIaG5qKjIdnD/J5100tiCggL/ns6dnp7e8OWXX27a/7Fbbrkl4v333w967733No0dO3bfxRdfPOqmm26KfPXVVw97N9TU1CTmzp0b+/e//337hAkTmiZOnHjUzTffvDswMLC9r+9/qDAcQKSUDiHENcDHdMxvPCelLBZC3AUUSCmXdT53mhBiPdAG3CylrOqt38HKavrggjmkLBTgF4KwwRUPv8F7S+Yb7jc2NtYFdkODSq6glq9Krkc63333ne9RRx3VvWa/vr7eGhAQ0P3HNyIiou3yyy+vvOSSS2IbGhqsK1as+DEkJOSAHEp2u729vLz8kE+zBweI3tixY4fHI488MiI/P399SkpKC8Dll1++Z8GCBaN7e53NZpM5OTnd51m7du2Pzp5zqHDJPhAp5YdSygQp5Rgp5ZLOxxZ1Bg9kBzdIKZOklEdLKV/rq8/Byink7+vDjLFe3cc/NNn5cYvxtftdE7QqoJIrqOWrkuuRzpo1a3xTUlKauo4DAgLa6urqDvibl5aW1lhSUmJbtGjRjrFjxx4yv1BfX28JDAx0GPH44IMPAlpbWy0nnHDCOLvdPt5ut4+fMWNGvN1uN57wz82YNpniYHLfH2fywbXPcVNwDjO8VzLvsQUse+gmd2tpNMrizLDSULNu3Trb+eef3z3ScdRRRzWtX7/eZ/LkyY0A3377re2GG24Y9dvf/rbqhRdeGH7llVdWH9zHxo0bfc4777xDHs/Ozo7vbQhrxYoVJV3H1dXV1qlTp9Z89NFHvQ7bq4hpU5kYKSjVF16eHlx23EiSxFbCxM/EeVbw1eofDPVpsx0yTGpaVHIFtXxVcj3S+fHHH33T09O770BOO+202hUrVtgBSktLPWfMmDH24Ycf3vr8889v27Bhg+2DDz44YElmU1OTKC4u9jvrrLPqDu57xYoVJY2Njat7+rd/8AA49thjG7/55ht7bm6uL0B1dbXlpZdeGtbebpqpjAFj2gDiimSKvbHwwjP4fzVTObPlHt5lCn969gtD/amURE8lV1DLVyXXI5lt27Z51NXVWVNTU7snza+44oqqL774IrCqqsp6+umnx1911VUVF1xwQa3dbm+/+uqryxctWhS5fx+vvvrqsIyMjPqYmBhDS2dPPfXUvTfeeOOuOXPmjPH19Z0wbty4lI8++ihgMD8kDxWG94EMFuPGjZM//GDsrqAvXvzwK+5Y8cuHi9szfZl3zkkD6is/P1+ZPx4quYJaviq5wpG/D+RgrrnmmsiwsLDWRYsW9ZkN45hjjjnq2Wef3XLssce6Nm2FYvS2D8S0cyBDcXt30Rkn8shHS9njF8MYsYP3PtjAZWdNHtDwWdeOZhVQyRXU8lXJ9dfI448/7vSKGTOuejIb6t9DGeQfV05jcstyPva8mTt4khdf+Le7lTQajUYJTBtA/P17XODgco5NHkuMDXbVt1Ows43/u+0Wamtr+91PZmbmINgNDiq5glq+KrlqNEYxbQDpSlA3FNz8l0Wc/Ladqz9sZuvOSu66665+91FaWjoIZoODSq6glq9KrhqNUUwbQFxVkdAZ/P39ufu+v3UfP/LII+SvXtuvPrqyuKqASq6glq9KrhqNUUwbQIaaOXPmMHnyZNLHRbP8hiSW/OuNIZnI12g0GlUxbQAZ6g1ZQgjufuAR/njusZzgu5UFgSu5/lHns/UmJSX13cgkqOQKavmq5KrRGMW0AcQd+1NOPDaV99tO5B+Os7ih9Sre3SLYtG2XU69ta1MnrY1KrqCWr0quGo1RTBtAXF1y1lmevOUy7qs9jXp8ET52Lnjwbade5+riV4OJSq6glq9Krpq+EUJMXLdu3eCkBj8CMG0AcRehQQH8MSMEAEE7Kf5VPPTSh2620mg0/SEyMvJoHx+fNF9f3wkhISGpM2fOjKmtrdV/71yMad9QLy+vvhsNEn++cDohe0t51vNBnvV6iJ9Wf8rO3Yck5DyAyMjIXp83Eyq5glq+Krke6bz22msljY2NqwsKCtYXFRX53XrrrSPc7XSkYdoAMhgVCfvDqzf9lk+bk6iUw2j0DmHOX1/ttX1UVNQQmRlHJVdQy1cl118LsbGxrSeffHLtDz/8YHvkkUdC4uLikv38/CZERUUd/cADDwzfv+0dd9wRHhoaekxYWNgxf//730P2f66pqUnMnz8/asSIEUeHhISknn/++aMaGhrE0H435sK0ubD27t3r1vMnjB6JJe5kTtl6MvX4gjc8/Or/uP6803tsn5+fz5QpU4ZWcoCo5Apq+ark6mru/mD9yGdzS536lH/WMSP2PHZ+2gHlXBe8Ujj6/bW7hh/uNX/Iit11x5lJO/vrtWnTJs/PP/88cPr06TXh4eGO999/f9O4ceNaPvroI/+ZM2fGZ2ZmNmZlZTW++eabAU888UTE//73vw2JiYn7LrjgggMqBl599dVRW7Zs8V6zZs16Ly8v+bvf/S5u4cKFI5cuXWq8Ip2imPYOxAzcfcXvsNT/krRzaV45pTsq3Gik0Wic5fzzzx9rt9vHZ2dnH3XcccfVL1myZNfcuXNrk5OTWywWC9OnT2844YQT6r788kt/gP/85z/Bc+bM2XPsscc2BwQEtC9ZsqQ7WLW3t/Pqq68Of+yxx7aHh4e3BQUFtf/lL3/Z9e677wa77zt0P6a9A7Fare5WwGKx8J+bz2Xaw8uZ61vAAt9X+Out7/PEC28d0naocne5ApVcQS1flVyPdF555ZVN5557bv3+j73++usBS5YsGbllyxaf9vZ2mpubLcnJyU0A5eXlnmlpad1DH/Hx8fu6vt61a5dHc3Oz5bjjjhu3f39tbW16CMuM+Pr6ulsBgKNiIrlqgh8py18manQzvps+4KWXXuLCCy88oF16+oBKKrgFlVxBLV+VXF3NHWcm7RzIEFMXj52ftvWx89nad8uB0dTUJC6++OIxTzzxxJbzzz//Z29vb3nqqaeO6dpzFh4e3rp9+/bu1TubNm3q/joiIsLh4+PTvnbt2uLY2Nihy7NkclwyhCWEOF0IsUEIsUkIcUsv7WYKIaQQos/fMnfPgezPLZecw2fWU7jg7UYeytvHVVddxU8//XRAm7y8PDfZ9R+VXEEtX5Vcf200NzeLffv2WcLCwlo9PT3l66+/HrBy5cqArudnz55d/frrrw9ftWqVT319veX2228f2fWc1Wpl7ty5e66++uroHTt2eEBHWdy33noroKdz/VowHECEEFZgKfAbIAk4TwhxSD4HIYQd+BOQ70y/ZstDdffDT/JtUwwADQ0NzD3vPPY2/bLZcSizBxtFJVdQy1cl118bQUFB7ffcc8+23//+92MCAwPHv/LKKyGnnHJKd+2G2bNn182fP79i2rRpiWPGjEk56aSTDqiFvnTp0rK4uLiWjIyMcf7+/hNOOeWUhB9++MFn6L8T82C4pK0QIhNYLKWc1nl8K4CU8v8d1O7vwGfATcBNUsqC3vpNTEyUZtvVu2rVKjIzMwkL8uOB85J4zfN3vPfADQDk5OQos/pGJVdQy1clV/j1lbTV9J/BLmkbCWzf77gMOKAotBBiAhAtpfxACHHT4ToSQswH5gOMHDmSnJwcAOLi4rDb7RQVFQEQEhJCcnIyK1as6PgmPDzIysqisLCQurqODw3p6elUVFSwfXuHWnx8PN7e3qxbtw6AsLAwEhISyM3NBcDb25vMzEwKCgpoaGgAICMjg7KyMnbs6Fill5iYyJ9uvZPprcs4yWsdwjGcO54cxSmJHSsPu+ph5+fnd5c2zczMpLS0tDvNd1JSEm1tbd0pLyIjI4mKiiI/v+PGzN/fn/T0dPLy8ro/zWZlZbFx40YqKztWhKWkpNDS0kJJSQkA0dHRhIeHU1DQEZMDAgJIS0sjNzcXh8MBQHZ2NsXFxVRVVQFQU1NDfX09mzdvBiAmJobg4GAKCwsBCAoKIjU1leXLlyOlRAjB5MmTKSoqoqamBoC0tDSqq6vZsmXLoF6n4447rvtnwdnrZLVaWb9+PQARERHExsZ2Dy/ZbLZBu07Dhw/vdnXFdUpNTR3066TRDBRX3IHMAqZJKed1Hl8ETJJSLug8tgBfAJdIKbcIIXJw4g4kJSVFdv0RMRMORxtnXf9Xrh/2Jbe1XkZlizfPX5BEuK9QJhPr+vXrlXEFtXxVcgV9B6Lpm97uQFwxiV4GRO93HAXsvxLDDqQAOUKILcBxwLK+JtKHsqBUf/DwsPLPv1zJ/JqL2U0QwtuXPzz3DRs2b3G3mtN03cmogkq+KrlqNEZxRQD5DogXQsQKIbyAucCyriellLVSyuFSyhgpZQzwDXB2X3cgZmb0iFDuPzcR6ehYJv67oB95/IN80038azQazWBiOIBIKR3ANcDHwA/A61LKYiHEXUKIswfa71AXlOovs089jnOimrnI+gkPeD7F38Pe5bzFT7tbyylSUlLcrdAvVPJVyVWjMYpL9oFIKT+UUiZIKcdIKZd0PrZISrmsh7ZTnLn7cEdBqf7y6HXnUdQczYb2KJ5z/Ib8fVH839PvuFurT1RbaqqSr0quGo1RTJsLy10FpfrLa0sW8Ns9V/GcfcRhAAAgAElEQVRO+4kAPP+j5LVPzL2ZrGtVkCqo5KuSq0ZjFNMGEFXw9fHmg9tmIes6kiyO9Gzgk/8s7V4yqdFoNEcqpg0g7iwo1V9iI8NZMi2K4JadvGtdyJMjl3HH70+murr3IlTuIjo6uu9GJkIlX5VcNX3Tn5K2TU1NYsyYMcnbtm0zbY7B/WlqahKxsbHJXalZBoJpA4i7C0r1l3NPOZ5bTovnlcI6Vm5v4/1vSzn77LNpbGx0t9ohhIeHu1uhX6jkq5LrkYw7Str+7W9/C83IyKgfNWqUA+Cnn37ynDZt2pigoKBUu90+PiEhIenRRx8N6auf/vLBBx/YLRbLRF9f3wl+fn4TYmJiUh555JE+z2Oz2eQFF1ywZ/HixREDPbdpA4iZkik6Q0FBAbN/M4XIS55l2kuN1LbAypUrmTlzJs3N5ppY7doJrQoq+arkeqQz1CVtn3/++dBLLrmkquv4/PPPj42MjNy3devW76urq9f861//Ko2IiBiUDW6hoaGtjY2Nq+vr61fffffdZTfeeGPMqlWr+szTdemll1a/+eabIU1NTQNKS2/aAKIqs+fM5b4HH+4+Hj+shhkLH2Jfq8ONVhrNr5ehKGlbUlLitX37du+TTjqp+5Pv999/7zdv3rw9AQEB7Z6enpxwwglNs2fPrgOYMmXK2CVLloTt339CQkLSiy++OAw6hs7uvffe0NGjR6f4+flNuPbaa0cWFxd7jx8//ih/f/8JZ5xxRlxzc/Mhf/QtFgsXXXTRz3a73bF27VpbX+cZM2ZMa0BAQNsXX3zhN5D31rQBxAwFpfpDQMAvWZ2vu+46brvtdhZeOIW/Jqznb4GvcNJNT+BwtLnR8Bf2d1UBlXxVcnU5iwMnsjhw4gGPPX/mWBYHTuT7NwK7H1v56HAWB07kzct+KRn78zZPFgdO5P64Yw54/dJJ41gcOJGtXw+4QFBXSdtjjjmmsaukbX19/eonn3yydNGiRdG5ubm+AF0lbT/88MONmzZtWvf5558fcDGvvvrqqE2bNvmsWbNm/aZNm74vLy/3Wrhw4UiAVatW2aKiolr2H3ofP3783muuuWb0U089FVRSUnLApO6FF15Y9Z///Ke7mmFeXp6tsrLSa9asWd3ZgT/55JPA1atXr1++fPkP//jHPyLmzZs3+pVXXtm8devWtRs2bLA9/fTTh1RDbGtr44UXXhhWX19vTUtLa3TmPGPGjGlevXr1gN5f0wYQsxSUcpa0tLQDju+6606Kw6ZT1B7H3x0z2WEbQ/ZN/zRFEDnY1eyo5KuS65HOUJa0rampsfr5+R3wy/3ee+/9dNxxx9Xff//9I4866qijjzrqqKTly5f7AlxwwQU/b9261ef777/3BvjXv/4VcuaZZ1b7+Ph0b4C75ZZbdgUHB7enp6c3x8fHN5188sl1SUlJ+0JCQtpOPvnk2v3/6O/evdvTbrePDw4OHv/Xv/515BNPPFGampra4sx5/P39237++ecBfWI37WqBrkyrqpCbm0tWVlb3scVi4Z37rmfyjd7ssMUBsNMnhuyb/snyB6/E08N9d1gHu5odlXxVcnU5i2tXHfLYJR9sOuSxE/60hxP+dGASxmGjWnt8/dXf/jBQnaEsaRscHNy2d+/eA36pQ0ND25544okdwI5du3Z5XHPNNVGzZs0aW15evtZms8np06dXP/fccyEPPvjgznfffTf4lVdeOaBK3ciRI7vHvX18fNrDw8O7509sNlt7RUWF537naq2oqFh78HvgzHkaGhqsw4YNG9AnW9PegaiwE31/utJx74+Hh5XlD/2R8KYtAMSIXdxvf5mzb7yf5pZ9h7QfKnpyNTMq+ark+mujq6TtddddV1FZWVlUX1+/ZvLkybX9LWlbX1+/putfY2PjaoD09PTGsrIy78MlgR0xYoTjlltuKd+9e7dnZWWlFeCyyy6reuutt4KXLVsWYLPZ2k899dRBWTnU13l++uknnwkTJgxouahpA8iRgoeHla8evILwxi0s9niBLGsxF9u/IfPGp6nf2+RuPY3mV8NglrQdM2ZM66hRo5pzcnK6J6OvuuqqyO+++86ntbWVmpoay6OPPho2atSoloiIiDaAU089da/FYmHhwoVRs2bNqmKQ6O08paWlnrW1tR4nn3zygIKXaQOI3W53t0K/yM7OPuxzXp4erPzblTzQcCYvO07hTsfvqfGPYdLNz1O+p2YILTvozdWMqOSrkuuvjcEuaXvppZfufv7557tXbjU2Nlpmzpw5NjAwcEJcXNzRZWVlXu+8884BQ3qzZ8+uKikpsV122WWDFkB6O89zzz0XPHPmzD02m21AQz6GC0oNFsnJybK4uNjdGk7z/fffc/TRR/fapr29nal//ic/eXQsPhG0E1/zNS/ceS0REQPey9NvnHE1Eyr5quQKuqCUK2lqahIpKSlJX3zxxcbRo0c7td/j8ccfD/nXv/41fNWqVYNav7un8zQ1NYmkpKSk3NzcDZGRkYcdex3sglKDgmpjyV0lSHvDYrHw6f1Xku7VsbjjJl7m3eFPcMesNH74YcBzhf3GGVczoZKvSq4a12Kz2eRPP/1U7GzwqK+vtzz99NOhl1566aAG4cOdx2azydLS0uLegkdfmDaAHKlYLBbevOtypg+vxnf7Sryskq07Kzn++OP57LPP3K2n0WiGgLfeeisgLCwsdfjw4a1XXHHFoH3qGOzzmHYIa8KECXL16tXu1nCampoagoKC+vWaD//7X+6+Zg7fbOmYv7Jardz50FJuXXA5FsvgxfaBuLoTlXxVcgU9hKXpGyWHsNra3L/hrj/U19f33eggzpg+naVvrWDkyI7FHqnxIziu8kXOvOF+ausHLxfYQFzdiUq+KrlqNEYxbQBRrbLb5s2bB/S6tLQ0vvvuOyYen80/ZwRziuf3zPAvYuLCl8lfNzjFiQbq6i5U8lXJ1QW0t7e3DygJn0YNOq9v++GeN20A+TUxcuRIPv3ovyzaO4eXHKfwoGM2joBIZj+7mr8+/7679TSaw7Fu9+7dgTqIHJm0t7eL3bt3BwLrDtfGtKlMvL2dquFiGmJiYgy9PijAn/8+fAuX3/8izVUWhBU8vH3w2/Q259y0kRcXX0WAv2vygxl1HWpU8lXJ1SgOh2NeeXn5M+Xl5SnoD6NHIu3AOofDMe9wDVwyiS6EOB14BLACz0gp7z3o+RuAeYAD2A1cJqXc2lufqk2i19XVuSwT67//u4JF/9vKDYGfc63HO3zXnsCcPfN5/IJjOeOE8Yb7d6XrUKCSr0quYGwSXaMx/KlBCGEFlgK/AZKA84QQSQc1Ww2kSymPAd4E7u+rXzNW8uuNwsJCl/V18fRsvlg4lf/VxlHQnsADrXNoDxjJVe9s4YK7/2W4togrXYcClXxVctVojOKK285JwCYp5WYp5T7gNeCc/RtIKb+UUnZFhG+AKBec94hmTFQE/33kdv5pu4b81o5svsLDk9DG9Vw+ayo//vijmw01Gs2vHVfMgUQC2/c7LgMyemn/B+Cjnp4QQswH5kNHbemcnBwA4uLisNvtFBUVARASEkJycjIrVqwAwMPDg6ysLAoLC6mr60hfk56eTkVFBdu3d6jFx8fj7e3NunUd80FhYWEkJCSQm5sLdMy5ZGZmUlBQ0J1KPiMjg7KyMnbs2AFAYmIiVquV9evXAxAREUFsbCx5eXk0NTWRn59PRkYG+fn5NDV1JErMzMyktLSU8vJyAJKSkmhra2PDho6MApGRkURFRZGfnw+Av78/6enp5OXl0dLSwoWZoznH4s+1r6wiNkhwr8czeBzTxtGTx3PCWRcyZ84c4uLiCA8P7y6nGhAQQFpaGrm5ud07+rOzsykuLqaqqoqmpiZqamqor6/vXjUUExNDcHBw9yfooKAgUlNTWb58OVJKhBBMnjyZoqIiamo68nelpaVRXV3Nli1bBvU6BQYGdv8sGL1OADabzeXXCehO497lmpKSQktLCyUlHavpoqOj+3WdAFJTUwf9Omk0A8XwHIgQYhYwTUo5r/P4ImCSlHJBD20vBK4BJkspe12nm56eLnV96V9obG7hqnueYMKq2/H3bOfy95sBGDduHPf8bSm/Pf0kNxtqVETPgWiM4IohrDIger/jKGDnwY2EEKcCtwFn9xU8QL0NWcuXLx/U/n19vPn3Pddz8n15PLkr+ZcnmmuwfvX/mH7tEgp+cG4PwmC7uhqVfFVy1WiM4ooA8h0QL4SIFUJ4AXOBZfs3EEJMAJ6kI3hUuuCcpmOoUsIcc8wx5OV9w8MPP4y/vz9/mzGSczzzOd++it89s4Zzb3+aypraXvswa/qaw6GSr0quGo1RDAcQKaWDjmGpj4EfgNellMVCiLuEEGd3NnsA8AfeEEKsEUIsO0x3yiLE0O2l8vDw4LrrrmPN9+t5tHUGbziy+ZtjJsLTmzWOkUy78zX++OBLNB2m6uFQuroClXxVctVojGLaZIp6DsR5nnrnC+799CfaA0YCklc9lxAiarmm9iImJMZzz/wZbq3BrjEveg5EYwTT7h7tWh2jCu5c0TJ/xsn8+MilnBNeS2jjVqItlYSKWnbZxvKfbTYSrnuBW//5VvdKH9VW36jkq5KrRmMU06YyUa2gVNeSSXfh5enBI9efz56f6/jTwxZqGhqo9+xMfeI/nIk/3sN1025g4gWLiIoe5VbX/uLu97Y/qOSq0RjFtAFEMzCGDwvglTuvYuvOSq557G2KmkM4w7aOmcN+5Nij2xl7+TyGBQ/npptu4pLL/kB46HB3K2s0GkUx7RzIrzkXlisp3VHBoqUvErHybnbWNPLK9x13djYfL+665nesaEvlsnNO5dzJE91senjM+t72hEquoOdANMYw7RyIagWlqqur3a3QI7GR4bz415u4fdk2jrngHkaMGAHApVPHcZPff7nd/y2u+6ic+Kuf5rpHX2dPrfn235j1ve0JlVw1GqOYNoCoVlCqK02EWQkMDGThwoWUlpaycOFCdoUdz4dtk3jFcQoArfaRrNy5j6fuvZ6zb3yA5z/Mpa3tsHVkhhSzv7f7o5KrRmMU0wYQzeDg7e3N6aefzltPL2XXxFv5qnYE0tGxX2Sm9Sv+4v0GV/h+wuIVtYy97mXOueNZvlr9g5utNRqNGTFtAFGtoFRcXJy7FZwmLi4OIQR/OPdkVi39E//7YzrHee9gdX0QH7el80bbZACkXzCtjmZ2PjKNq09P4oEHHnBLyVbV3luN5teCaQOI1arWxje73e5uBac52HVcXDSv3Tmflx9ewvaJf2FjQwDtzR2Zbs8hhzkxtSSKzfz5z39mzJgxpKWlMeu2x/n3R1+zzzH4c1Uqv7cazZGMaQOIagWlVNpAdjhXq9XKvHNP4etH/0TRnWdwUWwTn1aGcvuXrTy/prW7XYx/Ey97LKJ95d+J//PbTLruSf78z3f4cWv5kPqaEZVcNRqj6H0gmh4JCrRz9xUz4YqZ1NQsIX7ZMt544w0+/fRTxo+NwMoufpb+CC9fKvHliy21HPP0H3mgcSzbvROZmjqKC6dmMCJEnSWtGo2mf5g2gHh4mFatR0JCQtyt4DT9dQ0KCuLiiy/m4osvpr6+nsde/5jpazawiyDw62hzguV7LvT4nCi/PVzSOpmNa/axdM0Kple/SlT4CNKypnLccccRGBg46L7uRCVXjcYopt1IqFoyxfb2diwW044IHoCrXB0OB29/nsery7+nvqaC04dtY4uI5IP2TABGiQpWeF9P5d52wh9sQAhBcnIyR59yDj4jk5gyfixnHJfC8ADfIfEdClRyBb2RUGMM0waQxMRE2VVOVAVycnKYMmWKuzWcYrBcd5Tv5oUPv+Kz77ezucGDccEOFlYtZkdNCxe83ZEcUyCoviMcT9HGiS2PUE0A1uafiffYw8iwUCaNG82p6UmMiRjWnRpdv7eDhw4gGiOoNU6kMTWREaHcetlvubXzuHT7Dr79egTfrMxl/PivWLt2LSNGRlCPDW9aqaZjxVKbzzDu8FzKxL0b+f3Xt3Dvyp+htZnR7WWcHljOPrwJCwtj7NixeHl5ue8b1Gg0B2DaAKJaYR6V5myGyjU2OpLYObOZM2c2AA0NDXy58lvu+KaETeVVSMc2CByJ8PDELhrxEa1sl2EdL/b04TxLHlfyEX/5vJnkO+/EarUy6bgMjs3MoNISTlvgaBJGBpE6JpIJCaOIDvHD2811T1T6OdBojGLaISzV5kA0A6Ns5y7+m1vIyvXbKKvYQ2mLPw57BFZbADe3LGWu51dc8UEzb/3QkQTyglNTeOmEbXzSNpH5rTcCYKGdBz3/yfb24Tza9Bv8RSvDbYJxwz2Zmx5JdHQ0I0eOxNPT053fqinRQ1gaI5j245Jq+0AKCwtJS0tzt4ZTmMk1auQIrpg9nSv2e6y+vp6vC9dRWDyHv/50FLuG5TFq1E62bdtGnVcYn7aF8l17Ynf7EVTxW2sulZZhPOw1m1qg1gH37byRmGU7mPpiI0UVkvDwcE4+ZTIxUcMpbo+hynMkIX6ehAf6EhliZ3R4MLEjhxMa4EuInxfDfL3w8ujfhLiZ3luNZrAxbQBRLRtvXV2duxWcxuyudrudaZMzmTa5YzVX18T03r17KSr+kbzizZRsqSDm5++p2NtOrUVyw7Dz8PKxHdBPmLWO4d4WqpskUkrKy8s50baJK30+4p7WC3imbRzb6+DY+h+5r/IhcopT+X3rNd2vv7Hxcfyad/MdqQSHRhASEsJu/KnFDy97KMHD7IQNsxMREkCQvw1/bw9+3PkzEXFN+Pt44OflgdWi1lCsRtMfTBtANJqD8fPz4/hJEzl+0qG1SxoaGij5qZSEDaWs31rO1spa5pfPwLI1H4dvGaKuEiklq/aN5g3HCIplTPdrh4tahom9+PDLbnsrbSwI/pp2Kbn87lW0d470fnR9KqcHlHL5nht4tz0dqGaa5Vuu8PiA/7Zl8GzbdMj9An8auWnfU3i2t7HaIw273U5AQAA79+6jRVrZ5R2Hp48ffj5eBPoI/G0++Pn6Yffzwc/bE29PK9nxoUQE+hzwfW4or8dqAS+rFS8PC55WgaeHBS9rxz+LDliaIcQlAUQIcTrwCGAFnpFS3nvQ897AC8BEoAqYI6Xc0luffn5+rlAbMtLT1RlGVskVnPP19/dnQurRTEg9usfn9+3bx65du9iyrYwNW3cxurwK6571VNQ28V1jOxlt83BIQTPrsNgC8PX15e7SUXhVl3QHD4AmYaNK2qmR/t2PRYvdpFk2sbo9vvuxEFHHJQHf8lN1Oxc+9GX34xv+MpYEz0qm1txPifQGHCz0eJWrPN7nvta53Nt2NgApYjNprffwdgXcuTYCm82GzWbj2uO98LDAza3zaaIjuJxtWck4yzbeb8vkBzkaK5Ixll38OfRbLAEj+Tl4At7e3giLhdxvv6HN4kWBVzpBNv35UWMMwz9BQggrsBSYCpQB3wkhlkkp1+/X7A9AjZRyrBBiLnAfMKe3fltbW3t72nRUVFTg7+/fd0MToJIruMbXy8uL0aNHM3r0aCafePh2bW1tVFdXU1FRwZ49j7F7926e+O0e9uzp+Pd0rTf3VsDe1l1YHLtoabfwmocneX7nUW0ZTqt1O8LLlzpvWLJhODXlZQf0X+yIZq81gDr5y+ZJgcQhLbTwyyR/gGhkrH8zWysdrF27tvvxD7OHYxP7uLl1fvdjp1lXcab1G9a3j2a9jMGB4ChKOaX2DV7Jbe3egxPgb6P2Rk/q22wc3fAsjvpdht5TjcYVH0EmAZuklJsBhBCvAecA+weQc4DFnV+/CTwuhBCylyVg+/btc4Ha0LF9+3bGjBnjbg2nUMkVhtbXarUSGhpKaGio06+RUrJ3715qamr47LPPSEgYQ21tLbW191BbW8v9M+qpq6ujvr6el2o8qG9x4NGST/A+B82t7dzf5ss9bTNxyBZa25eB1ZNcL8Ex633YW7P7gHNdVfk7bJ6Cvfuqabd6I6wevOOdTrFHFOsZ3d1ui4zgiTVWvtva9Iun8OCjtmPZ1xWo2tWaZ9SYD1cEkEhg+37HZUDG4dpIKR1CiFogBNizfyMhxHxgPkBYWBg5OTlAR40Fu93enek0JCSE5ORkVqxY0fFNeHiQlZVFYWFh9wRxeno6FRUVbN/eoRYfH4+3tzfr1q2jq/+EhARyc3OBjvojmZmZFBQU0NDQkco8IyODsrIyduzYAUBiYiJWq5X16ztiY0REBLGxseTl5dHQ0EB+fj4ZGRnk5+fT1NTxi5uZmUlpaSnl5R2ZapOSkmhra6Nrl31kZCRRUVHk5+cDHUMx6enp5OXldVdlzMrKYuPGjVRWVgKQkpJCS0sLJSUlAERHRxMeHk7XsueAgADS0tLIzc3F4ehY/pqdnU1xcTFVVVU0NDRQU1NDfX19d32PmJgYgoODKSwsBDryX6WmprJ8+XKklAghmDx5MkVFRdTU1ACQlpZGdXV1dxW+wbpOUsrunwWj1wnAZrMN2nUaMWIEra2t+Pr6MmnSpH5fpxNPPJE1a9awc+dOmpuziYmJYc+ePWzevJmWlhYCAwPx9PQk6Mcf2bdvH1arlZCQEH78cRcR+1YQ3OpgeGgYO8oreafxWFoDWzjjDD/27t1LbX0DN30hEFaBb/vjNO7teP80moFieB+IEGIWME1KOa/z+CJgkpRywX5tijvblHUe/9TZpupw/aampkqVUmPv2LGDyMhId2s4hUquoJavSq6g94FojOGKrG9lQPR+x1HAzsO1EUJ4AIFAdW+dqrYTXaUKiiq5glq+KrlqNEZxRQD5DogXQsQKIbyAucCyg9osAy7u/Hom8EVv8x9A99CCKnQNuaiASq6glq9KrhqNUQzPgXTOaVwDfEzHMt7npJTFQoi7gAIp5TLgWeBFIcQmOu485ho9r0aj0Wjci0sWgkspPwQ+POixRft93QzM6k+fquUtCgsLc7eC06jkCmr5quSq0RjFtJVvVBtLTkhIcLeC06jkCmr5quSq0RjFtAGka4mmKnQtB1YBlVxBLV+VXDUao5g2gGg0Go3G3Jg2gKhUVxrUGnJTyRXU8lXJVaMxii4opdH8itEbCTVGMO3HfNUKSqkU7FRyBbV8VXLVaIxi2gCiWkEplSb9VXIFtXxVctVojGLaAKLRaDQac2PaOZC0tDTZlRlWBZqamrDZbH03NAEquYJaviq5gp4D0RjDtHcgqhWUKisr67uRSVDJFdTyVclVozGKaQOIagWlumpRqIBKrqCWr0quGo1RTBtANBqNRmNuTBtAfHx83K3QLxITE92t4DQquYJaviq5ajRGMW0AUa2glNVqdbeC06jkCmr5quSq0RjFtAFEtYJSXfW3VUAlV1DLVyVXjcYopg0gGo1GozE3pg0gqhWUioiIcLeC06jkCmr5quSq0RjFtAFEtaymsbGx7lZwGpVcQS1flVw1GqOYNoCollMoLy/P3QpOo5IrqOWrkqtGYxTTBhCNRqPRmBtDAUQIESyE+FQIUdL5f1APbcYLIfKEEMVCiLVCiDlOiSlWUEql/EcquYJaviq5ajRGMZRMUQhxP1AtpbxXCHELECSlXHhQmwRASilLhBAjgVXAOCnlz731rQtKaTSDj06mqDGC0Y/55wD/7vz638C5BzeQUm6UUpZ0fr0TqARC++p47969BtWGlvz8fHcrOI1KrqCWr0quGo1RPAy+PlxKuQtASrlLCBHWW2MhxCTAC/jpMM/PB+YDhIWFkZOTA0BcXBx2u52ioiIAQkJCSE5OZsWKFR3fhIcHWVlZFBYWUldXB0B6ejoVFRVs374dgPj4eLy9vVm3bh1d/SckJJCbmwt0rPrKzMykoKCgewI/IyODsrKy7gR5iYmJWK3W7s1iERERxMbGkpeXR0NDA/n5+WRkZJCfn9+9ETIzM5PS0lLKy8sBSEpKoq2tjQ0bNgAQGRlJVFRU9x8ef39/0tPTycvLo6WlBYCsrCw2btxIZWUlACkpKbS0tFBSUgJAdHQ04eHh3dXwAgICSEtLIzc3F4fDAUB2djbFxcVUVVXR0NBATU0N9fX1bN68GYCYmBiCg4PpSqEfFBREamoqy5cvR0qJEILJkydTVFRETU0NAGlpaVRXV7Nly5ZBvU6NjY3dPwtGrxN0DDMN1nWqqanpdjV6nQBSU1MH/TppNAOlzyEsIcRnQE+L228D/i2lHLZf2xop5SHzIJ3PjQBygIullN/0JZaYmCi7fnlVICcnhylTprhbwylUcgW1fFVyBT2EpTGG0TmQDcCUzruPEUCOlPKQbHJCiAA6gsf/k1K+4UzfEydOlKtWrRqw21DT0tKizN4VlVxBLV+VXEEHEI0xjM6BLAMu7vz6YuC9gxsIIbyAd4AXnA0eQPewgCqUlpa6W8FpVHIFtXxVctVojGI0gNwLTBVClABTO48RQqQLIZ7pbDMbyAYuEUKs6fw3vq+OVatI2DV2rgIquYJaviq5ajRGMTSJLqWsAk7p4fECYF7n1y8BLxk5j0aj0WjMh2l366m2ISspKcndCk6jkiuo5auSq0ZjFNMGECOT++6gra3N3QpOo5IrqOWrkqtGYxTTBpDm5mZ3K/QLlZYcq+QKavmq5KrRGMW0AUSj0Wg05sa0AcTLy8vdCv0iMjLS3QpOo5IrqOWrkqtGYxTTBhDVKhJGRUW5W8FpVHIFtXxVctVojGLaAKKTKQ4eKrmCWr4quWo0RjFtANFoNBqNuTFtALFare5W6Bf+/v7uVnAalVxBLV+VXDUaoxhKpjiY6IJSGs3go5Mpaoxg2jsQ1eZAumpNqIBKrqCWr0quGo1RTBtA2tvb3a3QL1TKHqySK6jlq5KrRmMU0wYQjUaj0Zgb086BqFZQyuFw4OFhtELw0KCSK6jlq5Ir6DkQjTFMewei2lDAxo0b3a3gNCq5glq+KrlqNEYxbTyC2xAAAAdESURBVABRraBUZWWluxWcRiVXUMtXJVeNxiimDSAajUajMTemDSCqFZRKSUlxt4LTqOQKavmq5KrRGMW0AcSsk/uHQ6U5G5VcQS1flVw1GqOYNoCoVlCqpKTE3QpOo5IrqOWrkqtGYxRDAUQIESyE+FQIUdL5f1AvbQOEEDuEEI8bOadGo9FozIHRO5BbgM+llPHA553Hh+NuYLmzHatWUCo6OtrdCk6jkiuo5auSq0ZjFKMB5Bzg351f/xs4t6dGQoiJQDjwibMdq1ZQKjw83N0KTqOSK6jlq5KrRmMUo1tmw6WUuwCklLuEEGEHNxBCWICHgIuAU3rrTAgxH5gPEBYWRk5ODgBxcXHY7XaKiooACAkJITk5mRUrVnR8Ex4eZGVlUVhYSF1dHQDp6elUVFSwfft2AOLj4/H29mbdunV09Z+QkEBubi4A3t7eZGZmUlBQQENDAwAZGRmUlZWxY8cOABITE7Faraxfvx6AiIgIYmNjycvLo6GhgdDQUDIyMsjPz6epqQmAzMxMSktLKS8vByApKYm2tjY2bNgAdJRAjYqK6i5E5O/vT3p6Onl5ed0TsllZWWzcuLF7j0FKSgotLS3d4+3R0dGEh4fTlb04ICCAtLQ0cnNzcTgcAGRnZ1NcXExVVRUNDQ2ccMIJ1NfXs3nzZgBiYmIIDg6msLAQgKCgIFJTU1m+fDlSSoQQTJ48maKiImpqagBIS0ujurqaLVu2DOp1qqioQAjhkusEHSv8Bus65ebm4uPj45LrBJCamjro10mjGSh9pjIRQnwGRPTw1G3Av6WUw/ZrWyOlPGAeRAhxDeArpbxfCHEJkC6lvKYvscTERNn1y6sCOTk5TJkyxd0aTqGSK6jlq5Ir6FQmGmP0eQcipTz1cM8JISqEECM67z5GAD1tw80EThRC/BHwB7yEEA1Syt7mS5QrKBUQEOBuBadRyRXU8lXJVaMxiqFkikKIB4AqKeW9QohbgGAp5Z97aX8JTt6B6IJSGs3go+9ANEYwOol+LzBVCFECTO08RgiRLoR4xkjHXePbqtA1l6ICKrmCWr4quWo0RjE0iS6lrKKHiXEpZQEwr4fHnweed7JvI2pDTtckqAqo5Apq+arkqtEYxbQ70TUajUZjbkxbUEq1OZD29nYsFjXisUquoJavSq6g50A0xjDtT3rX+nxVKC4udreC06jkCmr5quSq0RjFtAFEtbHkro1fKqCSK6jlq5KrRmMU0wYQjUaj0Zgb0wYQX19fdyv0i9TUVHcrOI1KrqCWr0quGo1RTBtA2tra3K3QL+rr692t4DQquYJaviq5ajRGMW0AUa2yW1eyOxVQyRXU8lXJVaMximkDiEaj0WjMjWn3gQgh6gF10vHCcGCPuyWcRCVXUMtXJVeARCml3d0SGjUxWg9kMNmg0gYnIUSBKr4quYJaviq5Qoevux006qKHsDQajUYzIHQA0Wg0Gs2AMHMAecrdAv1EJV+VXEEtX5VcQT1fjYkw7SS6RqPRaMyNme9ANBqNRmNidADRaDQazYAwZQARQpwuhNgghNjUWWvdtAghnhNCVAoh1rnbpS+EENFCiC+FED8IIYqFENe626k3hBA+QohvhRBFnb53utupL4QQViHEaiHEB+526QshxBYhxPdCiDV6Oa9mIJhuDkQIYQU20lFjvQz4DjhPSrnerWKHQQiRDTTw/9u3nxcbowCM499nMUqDbKTJVSxkY4FkMztJfoWlBStlQ5GF8k/IH8BGZFJDKYUppCm/mgklFpIyjZqFxKyEx+K+i1m413UXzrnT86m3+567errd7nPec86Fy7Y3lc7TjaQRYMT2tKTlwBRwqOLPVsCw7XlJQ8AkcMr2k8LROpJ0BtgGrLC9v3SebiR9ALbZHqQ/PkZFanwC2Q68s/3e9ndgDDhYOFNHth8Bn0vn6IXtT7anm/tvwBtgTdlUnbltvhkONVddM54FJLWAfcDF0lki/ocaC2QN8HHBeIaKf+QGlaR1wBbgadkk3TVLQi+AOWDCds15LwBngV+lg/TIwD1JU5KOlw4Tg6fGAtEf3qt21jmIJC0DxoHTtr+WztON7Z+2NwMtYLukKpcJJe0H5mxPlc7yD0ZtbwX2ACea5diIntVYIDPA2gXjFjBbKMui0+wljANXbd8onadXtr8AD4HdhaN0MgocaPYVxoAdkq6UjdSd7dnmdQ64SXv5OKJnNRbIc2CDpPWSlgCHgVuFMy0Kzab0JeCN7fOl8/yNpFWSVjb3S4GdwNuyqf7M9jnbLdvraH9n79s+UjhWR5KGm4MUSBoGdgHVnySMulRXILZ/ACeBu7Q3ea/bfl02VWeSrgGPgY2SZiQdK52pi1HgKO3Z8Yvm2ls6VBcjwANJr2hPLCZsV388dkCsBiYlvQSeAbdt3ymcKQZMdcd4IyJiMFT3BBIREYMhBRIREX1JgURERF9SIBER0ZcUSERE9CUFEhERfUmBREREX34DB4enCGX9fVMAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.figure()\n",
"plt.plot(xx,yy,'k',lw=3,label='$f(x)=' + sb.latex(f)+'$')\n",
"plt.plot(xx,yy_Pade,'--',lw=3,label='Pade')\n",
"plt.plot(xx,yy_sympy,':',label='Pade(SymPy)')\n",
"plt.legend()\n",
"plt.grid()\n",
"plt.ylim([-0.5,1])\n",
"plt.gca().set_aspect(1/plt.gca().get_data_ratio())\n",
"plt.legend(fontsize=12,bbox_to_anchor=(1.04,0.8), loc=\"upper left\")"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$$f_{\\rm{Pade}}(x) =\\frac{0.03 x^{2} - 0.33 x + 1.0}{0.03 x^{3} + 0.2 x^{2} + 0.67 x + 1.0}$$"
],
"text/plain": [
"<IPython.core.display.Math object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/latex": [
"$$f_{\\rm{Pade(SymPy)}}(x) =\\frac{0.03 x^{2} - 0.33 x + 1.0}{0.03 x^{3} + 0.2 x^{2} + 0.67 x + 1.0}$$"
],
"text/plain": [
"<IPython.core.display.Math object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"text_sympy = r\"$f_{\\rm{Pade(SymPy)}}(x) =\"+ sb.latex(RoundFloat(R_sympy.expand().simplify(),2))+'$'\n",
"display(Math(text_Pade))\n",
"display(Math(text_sympy))"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment