{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "DhcM3fYEm0m3" }, "source": [ "# Problem 6\n", "A flea starts at $(0,0)$ on the infinite two-dimensional integer lattice and executes a biased random walk: At each step, it hops north or south with probability $\\frac{1}{4}$, east with probability $\\frac{1}{4} + \\epsilon$, and west with probability $ \\frac{1}{4} - \\epsilon$. The probability that the flea returns to $(0,0)$ sometime during its wanderings is $\\frac{1}{2}$. What is $\\epsilon$?" ] }, { "cell_type": "markdown", "metadata": { "id": "dJTH5KpWrQ2j" }, "source": [ "Here, we mostly follow the notation and exposition given in Chapter 6 of [Bornemann's book](https://epubs.siam.org/doi/book/10.1137/1.9780898717969)\n", "\n", "### Modelling as a Linear System\n", "\n", "Let $p_N, p_S, p_E, p_W$ be the probabilities that the flea goes north, south, east, and west respectively, and also define $q(x,y)$ to be the probability that the flea returns to the origin from this point. \n", "\n", "We therefore have that $q(0,0) = 1$ and that the probability that the flea ever returns to the origin is $$ p = p_E q(1,0) + p_W q(-1,0) + p_N q(0,1) + p_S q(0,-1) $$\n", "\n", "Conditioning on each direction that the flea can travel, we also have that $q$ satisfies the following linear equation:\n", "\n", "$$q(x,y) = p_E q(x+1, y) + p_W q(x-1,y) + p_N q(x,y+1)+ p_S q(x,y-1) $$\n", "\n", "Instead of dealing of $q$, we instead approximate it with $\\hat{q}(x,y) = \\begin{cases} q(x,y) \\; \\; |x|, |y| < N \\\\ 0 \\; \\; \\text{otherwise}\\end{cases}$\n", "\n", "Using the same recurrence, we then get a system of $(2N+1)^2$ linear equations that we can solve for the values of $\\hat{q}$" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 34 }, "executionInfo": { "elapsed": 820, "status": "ok", "timestamp": 1603144959944, "user": { "displayName": "Husnain Raza", "photoUrl": "", "userId": "16146658556003895350" }, "user_tz": 300 }, "id": "6n4_1lnimtlv", "outputId": "87e26f43-04eb-49cd-88dd-711812c8034d", "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABL/klEQVR4nO3deVxU5eIG8GcWZoZ9kV1RxH1DEBVBsSxyS9O0wi0Tl9TckrL0Vtr2SyvLfd/Nct8qTS3KDVEUEU1xBQUXQEAWQbaZ8/vDossVjZGZObM8389nPrcOZw4P59rweM573lciCIIAIiIiIjMhFTsAERERkS6x3BAREZFZYbkhIiIis8JyQ0RERGaF5YaIiIjMCssNERERmRWWGyIiIjIrcrEDGJpGo8Ht27dhb28PiUQidhwiIiKqBkEQUFBQAG9vb0ilT742Y3Hl5vbt2/Dx8RE7BhERET2FtLQ01KlT54n7WFy5sbe3B/Dw5Dg4OIichoiIiKojPz8fPj4+Fb/Hn8Tiys3ft6IcHBxYboiIiExMdYaUcEAxERERmRWWGyIiIjIrLDdERERkVlhuiIiIyKyw3BAREZFZYbkhIiIis8JyQ0RERGaF5YaIiIjMCssNERERmRWWGyIiIjIrLDdERERkVlhuiIiIyKyw3OiIIAiIv5GD0nKN2FGIiIgsmsWtCq4vKVmF6L8kFtZWMrT1dUZIg1oI8auFVrUdIZexQxIRERkKy42O3Mp9gFq2CmQXluLIlSwcuZIFALBTytG+vgtCG9RCB79aaO7lAKn035drJyIioqcjEQRBEDuEIeXn58PR0RF5eXlwcHDQ6bEFQcCVzPs4djULscnZOJ6cg7wHZZX2cbS2Qgc/F4T41UJoQ1c0creDRMKyQ0RE9CTa/P5mudEjtUZA0p18xF7LRmxyNuJScnC/pLzSPq52CnTwq1VxG6u+qy3LDhER0f9guXkCQ5ab/1Wu1uDcrTzEJmcj9lo2Tl7PQXFZ5QHIng6qh0Xnr7Lj42Jj0IxERETGiOXmCcQsN/+rpFyNxLS8v67sZOH0jVyUqiuXnTrO1gitKDuu8HRUiZSWiIhIPCw3T2BM5eZ/FZepcfrGPRz76zZWYlouyjWV/+/pF1gbM3q3gKONlUgpiYiIDI/l5gmMudz8r8KScpy8nlMxZufcrTwIAuDhoMTMfq3wXFMPsSMSEREZhDa/v0WfgGXRokXw9fWFSqVCcHAw4uLinrh/bm4uxo0bBy8vLyiVSjRu3Bh79+41UFrDslXK8WwTd0zr2Qw/ju+E7WND4edmi4z8EgxfewpTtiYiv7js3w9ERERkQUQtN5s3b0ZUVBRmzJiB06dPo3Xr1ujWrRsyMzOr3L+0tBQvvPACrl+/jm3btuHSpUtYsWIFateubeDk4mhT1xl7J4ZhVFh9SCTA1vib6DbnMA5dvit2NCIiIqMh6m2p4OBgtGvXDgsXLgQAaDQa+Pj4YMKECZg6deoj+y9duhRff/01Ll68CCurpxtzYkq3pZ7k1PUcvLs1EdeziwAAA9v74D89m8FexbE4RERkfkzitlRpaSni4+MRHh7+TxipFOHh4YiNja3yPT/++CNCQkIwbtw4eHh4oGXLlvjiiy+gVqsf+31KSkqQn59f6WUO2vq64JdJnRHZ0RcAsDEuDd3nHsHRv2ZGJiIislSilZusrCyo1Wp4eFQeFOvh4YH09PQq35OcnIxt27ZBrVZj7969+Oijj/DNN9/g888/f+z3mTlzJhwdHStePj4+Ov05xGStkGFG7xbY9GYH1HWxwa3cBxiy6gQ+3HUOhf8zWSAREZGlEH1AsTY0Gg3c3d2xfPlyBAUFISIiAh988AGWLl362PdMmzYNeXl5Fa+0tDQDJjaMDn618MukMAwNqQcA2HA8Fd3mHkbstWyRkxERERmeaOXG1dUVMpkMGRkZlbZnZGTA09Ozyvd4eXmhcePGkMlkFduaNWuG9PR0lJaWVvkepVIJBweHSi9zZKuU49M+LfHDqGDUcbbGzXsPMHDFcXz843kUlfIqDhERWQ7Ryo1CoUBQUBCio6Mrtmk0GkRHRyMkJKTK93Ts2BFXr16FRvPPLL6XL1+Gl5cXFAqF3jObgtAGrtj3dmcMDq4LAFh77Dp6zDuCuJQckZMREREZhqi3paKiorBixQqsW7cOSUlJGDt2LAoLCxEZGQkAGDp0KKZNm1ax/9ixY5GTk4NJkybh8uXL2LNnD7744guMGzdOrB/BKNkp5fi/l1vhuxHt4e2owo3sIkQsj8WnP13Ag9LHD74mIiIyB3Ixv3lERATu3r2L6dOnIz09HQEBAdi3b1/FIOPU1FRIpf/0Lx8fH+zfvx+TJ0+Gv78/ateujUmTJuH9998X60cwamGN3LBvcmd8sScJm06mYXVMCv64lInZr/ojqJ6L2PGIiIj0gssvWIiDlzIxdfs5pOcXQyIBRoX5IeqFxlBZyf79zURERCIziXluyLCebeKO/ZM745WgOhAEYPnhZLw4/wgSUu+JHY2IiEinWG4siKO1FWa/2hqr3mgLd3slrt0tRP8lx/DlvosoKedYHCIiMg8sNxbo+WYeODC5M14OrA2NACw5eA295h/F2Zu5YkcjIiKqMZYbC+Vko8CciAAsfz0IrnZKXMm8j5cXH8OmuFSxoxEREdUIy42F69rCE79O7owXW3lBrREwdcc5zP3tMixsnDkREZkRlhuCs60CCwcFYsJzDQEAc3+7gmk7zqFcrfmXdxIRERkflhsCAEgkErzTtQk+79sSUgmw6WQaRn8Xz6UbiIjI5LDcUCVDOtTD0iFBUMqliL6YiUErTiD7fonYsYiIiKqN5YYe0bWFJ34YFQwnGyucScvFK0tjkZpdJHYsIiKiamG5oSoF1XPBtjGhqO1kjZSsQvRbcgx/3soTOxYREdG/Yrmhx2robocdb4WimZcDsu6XIGJZLA5fvit2LCIioidiuaEn8nBQYcvoDujYsBYKS9UYvvYktsffFDsWERHRY7Hc0L+yV1lhzbD26BPgjXKNgHe2JmLxwaucC4eIiIwSyw1Vi0IuxZzXAjC6sx8A4Kt9lzDjx/NQa1hwiIjIuLDcULVJpRJM69kM03s1h0QCrI+9gXHfn0ZxGRfdJCIi48FyQ1ob3qk+FgwMhEImxb7z6Xh91QnkFpWKHYuIiAgAyw09pV7+3lg/oj3sVXKcvH4PryyNxa3cB2LHIiIiYrmhp9fBrxa2jQmFp4MKVzPvo9/iGCTdyRc7FhERWTiWG6qRJp722PFWKBp72CEjvwSvLY3FsWtZYsciIiILxnJDNebtZI2to0PRvr4LCkrKMWz1SfyUeFvsWEREZKFYbkgnHG2ssH54e/Rs5YlStQYTNiZg5ZFksWMREZEFYrkhnVFZybBgYBsMC/UFAHy+Jwmf/3wBGs6FQ0REBsRyQzolk0owo3dzTOvRFACw8mgKJm0+g5JyzoVDRESGwXJDOieRSDD6mQaYE9EacqkEPyXexrDVJ5FfXCZ2NCIisgAsN6Q3LwfWwZrIdrBVyBCbnI3XlsbibkGJ2LGIiMjMsdyQXoU1csPm0SFws1fiYnoBIpbF4k4eJ/sjIiL9YbkhvWtZ2xFbR4egtpM1krMK8erSWKRmF4kdi4iIzBTLDRmEr6stNo/uAN9aNrh57wFeXXYMVzPvix2LiIjMEMsNGUwdZxtsGR2CRu4PZzOOWBaLC7e5XAMREekWyw0ZlLuDCptHh6CFtwOyC0sxYHksElLviR2LiIjMCMsNGZyLrQI/jOqANnWdkF9cjiErT+BEcrbYsYiIyEyw3JAoHK2t8N2IYIT41UJhqRpvrInDoct3xY5FRERmgOWGRGOrlGNNZDt0aeKG4jINRq07hf3n08WORUREJo7lhkSlspJh2ett0aPlwwU33/r+NHafuSV2LCIiMmEsNyQ6hVyKBQMD0S+wNtQaAW9vPoMtJ9PEjkVERCaK5YaMglwmxexXW2NwcF0IAvDe9rNYG5MidiwiIjJBLDdkNKRSCT7v2xIjO9UHAHz80wUsPnhV5FRERGRqWG7IqEgkEnzwYjNMfL4RAOCrfZcwe/8lCIIgcjIiIjIVLDdkdCQSCaJeaIypPZoCABb+cRWf/ZzEgkNERNXCckNGa8wzDfBpnxYAgNUxKfjPznNQa1hwiIjoyVhuyKgNDfHF16/4QyoBNsal4Z0tZ1Cu1ogdi4iIjBjLDRm9V9v6YN6AQMilEuw6cxvjfjiNknK12LGIiMhIsdyQSejd2htLhgRBIZNi//kMvLk+HsVlLDhERPQolhsyGS8098CqYW2hspLi0OW7GLYmDvdLysWORURERoblhkxKWCM3rB8eDDulHMeTc/D6qhPIKyoTOxYRERkRlhsyOe3ru+D7kcFwtLZCQmouBq44juz7JWLHIiIiI8FyQyaptY8TNr3ZAa52Cly4k4+I5ceRkV8sdiwiIjICLDdkspp5OWDz6BB4OqhwNfM+Xll6DKnZRWLHIiIikbHckElr4GaHrWNCUNfFBmk5D/DqsmO4klEgdiwiIhIRyw2ZPB8XG2wdE4JG7nbIyC/Ba8tice5mntixiIhIJCw3ZBY8HFTYMjoE/nUcca+oDANXHMeJ5GyxYxERkQhYbshsONsq8P3IYATXd8H9knIMXR2HPy5lih2LiIgMjOWGzIq9ygrrhrfHc03dUVKuwZvrT2HP2TtixyIiIgNiuSGzo7KSYemQIPTy90KZWsCEjaex5WSa2LGIiMhAWG7ILCnkUswbEIiB7X2gEYD3tp/FyiPJYsciIiIDYLkhsyWTSvDFy63wZmc/AMDne5Iw59fLEARB5GRERKRPLDdk1iQSCab1aIp3XmgMAJgXfQWf/ZzEgkNEZMZYbsjsSSQSTHi+EWb0bg4AWB2Tgve3n4Vaw4JDRGSOWG7IYkR2rI+vX/GHVAJsOXUTEzcmoLRcI3YsIiLSMaMoN4sWLYKvry9UKhWCg4MRFxf32H3Xrl0LiURS6aVSqQyYlkzZq219sHhwG1jJJNhz7g5GrT+FB6VqsWMREZEOiV5uNm/ejKioKMyYMQOnT59G69at0a1bN2RmPn7yNQcHB9y5c6fidePGDQMmJlPXvaUXVr3RDiorKQ5dvos3Vschv7hM7FhERKQjopebb7/9FqNGjUJkZCSaN2+OpUuXwsbGBqtXr37seyQSCTw9PSteHh4eBkxM5qBzYzdsGBEMe6UccddzMGjFceQUloodi4iIdEDUclNaWor4+HiEh4dXbJNKpQgPD0dsbOxj33f//n3Uq1cPPj4+6NOnD86fP//YfUtKSpCfn1/pRQQAbX1dsPHNDnCxVeDPW/l4bVks0vOKxY5FREQ1JGq5ycrKglqtfuTKi4eHB9LT06t8T5MmTbB69Wrs3r0bGzZsgEajQWhoKG7evFnl/jNnzoSjo2PFy8fHR+c/B5mulrUdsWV0CDwdVLiaeR+vLjuG1OwisWMREVENiH5bSlshISEYOnQoAgIC8Mwzz2DHjh1wc3PDsmXLqtx/2rRpyMvLq3ilpXEafqqsobsdto4JgW8tG6TlPMArS4/hUnqB2LGIiOgpiVpuXF1dIZPJkJGRUWl7RkYGPD09q3UMKysrBAYG4urVq1V+XalUwsHBodKL6H/5uNhgy5gQNPGwR2ZBCSKWxyIxLVfsWERE9BRELTcKhQJBQUGIjo6u2KbRaBAdHY2QkJBqHUOtVuPcuXPw8vLSV0yyEO72Kmwe3QGtfZyQW1SGQSuOI/ZattixiIhIS6LfloqKisKKFSuwbt06JCUlYezYsSgsLERkZCQAYOjQoZg2bVrF/p9++ikOHDiA5ORknD59GkOGDMGNGzcwcuRIsX4EMiNONgp8PzIYIX61UFiqxrA1cfj9Ysa/v5GIiIyGXOwAERERuHv3LqZPn4709HQEBARg3759FYOMU1NTIZX+08Hu3buHUaNGIT09Hc7OzggKCsKxY8fQvHlzsX4EMjN2SjnWRLbD+B9O47ekTLy5Ph5zIgLQu7W32NGIiKgaJIKFrSCYn58PR0dH5OXlcfwNPVGZWoN3tyZi95nbkEiAz/q0xJAO9cSORURkkbT5/S36bSkiY2Ulk2LOawEYHFwXggB8uOtPfPvrZa4oTkRk5FhuiJ5AKpXg874tMfH5RgCA+dFX8J+d51Cu5oKbRETGiuWG6F9IJBJEvdAYn/dtCYkE2BiXhjEbTqO4jAtuEhEZI5Ybomoa0qEelgxuA4Vcit+SMjB45QnkFnE9KiIiY8NyQ6SF7i29sGFEMBxUcsTfuIdXl8bidu4DsWMREdF/Ybkh0lL7+i7YOiYUng4qXMm8j36Lj+FyBpdrICIyFiw3RE+hiac9tr8ViobudkjPL8YrS44hLiVH7FhERASWG6KnVtvJGtvGhKBNXSfkF5djyKoT2H++6tXsiYjIcFhuiGrg4XINHRDezAOl5RqM3RCP70/cEDsWEZFFY7khqiFrhQxLh7TBgHY+0AjABzv/xBxO9kdEJBqWGyIdkMukmNmvVcVkf/Oir+A/O//kZH9ERCJguSHSkUcn+0vF2O852R8RkaE91argGo0GV69eRWZmJjSayn8z7dy5s06CEZmqIR3qwdVOgYmbzuDXCw8n+1v1Rls42SjEjkZEZBG0XhX8+PHjGDRoEG7cuPHImAKJRAK12rj/lspVwclQ4lJyMHLdSeQXl6ORux3WDW8PbydrsWMREZkkva4KPmbMGLRt2xZ//vkncnJycO/evYpXTg7n+SD6Gyf7IyISh9ZXbmxtbZGYmIiGDRvqK5Ne8coNGdqt3Ad4Y3Ucrmbeh4NKjlXD2qGdr4vYsYiITIper9wEBwfj6tWrTx2OyNI8MtnfSk72R0SkT1oPKJ4wYQLeeecdpKeno1WrVrCysqr0dX9/f52FIzIXf0/2N2HjafyWlImxG+LxWd+WGBxcT+xoRERmR+vbUlLpoxd7JBIJBEHggGKif1Gu1uDDXX9i08k0AMCk5xvh7fBGkEgkIicjIjJu2vz+1vrKTUpKylMHI7J0f0/2526vxPzfr2Je9BVkFpTgsz4tIJdx2ikiIl3QqtyUlZXhueeew88//4xmzZrpKxORWZNIJIjq2gRuDipM3/0nNsalIut+CeYPCIS1QiZ2PCIik6fVXxWtrKxQXFysryxEFuX1DvWwZHAbKORS/HohAwOWxyKzgP99ERHVlNbXwceNG4cvv/wS5eXl+shDZFG6t/TChhHBcLKxQuLNPLy8iHPhEBHVlNYDil9++WVER0fDzs4OrVq1gq2tbaWv79ixQ6cBdY0DiskYpWQVYvjak0jJKoS9Uo7FQ9ogrJGb2LGIiIyGXue5cXJyQv/+/dGtWzd4e3vD0dGx0ouItFff1RY7xoaiva8LCkrKMWzNSWyMSxU7FhGRSdL6yo2p45UbMmYl5Wq8v+0sdp25DQAY/Ywf3u/WFFIpHxUnIsum1ys3RKQ/SrkMcyIC8HZ4IwDAskPJGPfDaRSXGff8UURExkTreW7q16//xAnHkpOTaxSIyNJJJBK8Hd4Y9WrZ4L1tZ/HLn+m4k3ccK4a2hZu9Uux4RERGT+ty8/bbb1f697KyMiQkJGDfvn2YMmWKrnIRWbyXA+vA29EaozfE40xaLl5eHIM1w9qhkYe92NGIiIyazsbcLFq0CKdOncKaNWt0cTi94ZgbMjXJd+9j+NqTuJ5dBHuVHEsGB6FTI1exYxERGZQoY2569OiB7du36+pwRPQXPzc77HirI9r5OqOguBzD1sRh80k+SUVE9Dg6Kzfbtm2Di4uLrg5HRP/FxVaBDSOD0SfAG+UaAe9vP4cv912ERmNRDzsSEVWL1mNuAgMDKw0oFgQB6enpuHv3LhYvXqzTcET0D6VchrkRAfCtZYt50Vew5OA1pGYX4ZvXWkNlxTWpiIj+pnW56dOnT6VyI5VK4ebmhmeffRZNmzbVaTgiqkwikWDyC41R18UGU3ecxZ5zd3A77wFWDG0LVzs+SUVEBHASP7HjED2148nZGP1dPPIelMHHxRprhrVDQ3c+SUVE5kmvA4plMhkyMzMf2Z6dnQ2ZjJfGiQylg18t7HgrFPVq2SAt5wFeXnwMMVezxI5FRCQ6rcvN4y70lJSUQKFQ1DgQEVVfAzc77HyrI9rWe/gk1Rur47DlZJrYsYiIRFXtMTfz588H8PCe/8qVK2FnZ1fxNbVajcOHD3PMDZEI/n6S6r1tZ/Fj4m28t/0srmcX4t2uTbgmFRFZpGqPualfvz4A4MaNG6hTp06lW1AKhQK+vr749NNPERwcrJ+kOsIxN2SuBEHAnF8vY/7vVwEAL/p74ZtX+SQVEZkHbX5/V/vKTUpKCgCgS5cu2LFjB5ydnWuWkoh0SiKRIKprE9SrZfvwSaqzd3An9+GTVLX4JBURWRCtx9z88ccfcHZ2RmlpKS5duoTy8nJ95CKip9Q/qA7WDw+Gg0qO06m56Ls4BlczC8SORURkMFqXmwcPHmDEiBGwsbFBixYtkJr6cBr4CRMmYNasWToPSETaC2lQCzvHdURdl3+epDp0+a7YsYiIDELrcjN16lQkJibi4MGDUKlUFdvDw8OxefNmnYYjoqf38Emq0IonqSLXxGHlkeTHPvFIRGQutC43u3btwsKFC9GpU6dKMxW3aNEC165d02k4IqqZWnZKfD8qGK+1rQONAHy+JwlTtp1FSbla7GhERHqjdbm5e/cu3N3dH9leWFhYqewQkXFQymX4sr8/ZvRuDqkE2BZ/EwOXH0dmQbHY0YiI9ELrctO2bVvs2bOn4t//LjQrV65ESEiI7pIRkc5IJBJEdqyPdcPbVww07rMwBudu5okdjYhI57ReOPOLL75Ajx49cOHCBZSXl2PevHm4cOECjh07hkOHDukjIxHpSFgjN+we3wkj153EtbuFeGXpMXz9amu81Npb7GhERDqj9ZWbTp06ITExEeXl5WjVqhUOHDgAd3d3xMbGIigoSB8ZiUiH6rvaYue4jujSxA0l5RpM3JiAr/dfhEbDgcZEZB60WhW8rKwMo0ePxkcffVQxY7Gp4QzFRA+pNQK+2n8Ryw4lAwDCm3lg7oAA2Cm1vqBLRKR3elsV3MrKCtu3b69ROCIyDjKpBNN6NMOciNZQyKX4LSkD/RbHIDW7SOxoREQ1ovVtqb59+2LXrl16iEJEYng5sA62jA6Bu70SlzPu46VFR3HsWpbYsYiInprW158bNWqETz/9FDExMQgKCoKtrW2lr0+cOFFn4YjIMAJ8nPDj+E4Y/d0pJN7Mw+ur4vBx7+Z4PcRX7GhERFrTaswNgCeOtZFIJEhOTq5xKH3imBuixysuU+P97Wex+8xtAMDg4LqY0bsFFHKtL/ISEemUNr+/tS43po7lhujJBEHA0kPJ+Gr/RQgCEFzfBYsHt+HK4kQkKr0NKCYi8yeRSDD22QZYObQt7JRynEjJQZ9FMUi6ky92NCKiamG5IaIqPd/MAzvfCkW9Wja4ee8B+i85hn1/posdi4joX7HcENFjNfKwx+5xHdGxYS0UlaoxZkM8FkRf4criRGTUWG6I6ImcbBRYF9kew0J9AQDf/HoZ439IQFFpubjBiIgeo1rlpl+/fsjPf3i/ff369SgpKdFrKCIyLnKZFB+/1AIz+7WClUyCPefu4NWlsbiV+0DsaEREj6hWufn5559RWFgIAIiMjERenm5XEl60aBF8fX2hUqkQHByMuLi4ar1v06ZNkEgk6Nu3r07zEFHVBravi+9HdkAtWwXO385Hn4VHcep6jtixiIgqqdYkfk2bNsW0adPQpUsXCIKALVu2PPYxrKFDh2oVYPPmzYiKisLSpUsRHByMuXPnolu3brh06RLc3d0f+77r16/j3XffRVhYmFbfj4hqpn19F+we3xGj1scj6U4+Bq44js/6tMSA9nXFjkZEBKCa89wcO3YMUVFRuHbtGnJycmBvbw+JRPLowSQS5ORo97e44OBgtGvXDgsXLgQAaDQa+Pj4YMKECZg6dWqV71Gr1ejcuTOGDx+OI0eOIDc397FLQpSUlFS6jZafnw8fHx/Oc0NUQ0Wl5XhnSyJ++esJqoHt6+Ljl5pDKZeJnIyIzJHO57kJDQ3F8ePHcffuXQiCgMuXL+PevXuPvLQtNqWlpYiPj0d4ePg/gaRShIeHIzY29rHv+/TTT+Hu7o4RI0b86/eYOXMmHB0dK14+Pj5aZSSiqtko5Fg0qA3e7doYEgmwMS4VEcuO404ex+EQkbi0floqJSUFbm5uOvnmWVlZUKvV8PDwqLTdw8MD6elVz6dx9OhRrFq1CitWrKjW95g2bRry8vIqXmlpaTXOTUQPSaUSjH+uEdYMawdHayucSctF7wVHcTw5W+xoRGTBtF44s169esjNzcWqVauQlJQEAGjevDlGjBgBR0dHnQf8bwUFBXj99dexYsUKuLq6Vus9SqUSSiWnjSfSp2ebuOOn8Z0wesPDcTiDV57ABz2bIbKjb5W3sImI9EnrKzenTp1CgwYNMGfOHOTk5CAnJwdz5sxBgwYNcPr0aa2O5erqCplMhoyMjErbMzIy4Onp+cj+165dw/Xr19G7d2/I5XLI5XKsX78eP/74I+RyOa5du6btj0NEOlK3lg12jA1FnwBvqDUCPv35At7efAYPStViRyMiC6P1wplhYWFo2LAhVqxYAbn84YWf8vJyjBw5EsnJyTh8+LBWAYKDg9G+fXssWLAAwMMBxXXr1sX48eMfGVBcXFyMq1evVtr24YcfoqCgAPPmzUPjxo2hUCie+P24cCaRfgmCgDUx1/F/e5Og1gho6mmP5a+3Rd1aNmJHIyITptdVwa2trZGQkICmTZtW2n7hwgW0bdsWRUVFWoXdvHkz3njjDSxbtgzt27fH3LlzsWXLFly8eBEeHh4YOnQoateujZkzZ1b5/mHDhj3xaan/xXJDZBjHk7Mx/ofTyLpfCkdrK8wbEIBnmzx+egcioifR66rgDg4OSE1NfWR7Wloa7O3ttT0cIiIiMHv2bEyfPh0BAQE4c+YM9u3bVzHIODU1FXfu3NH6uEQkrg5+tfDThE4I8HFC3oMyRK49iYW/X4FGw3WpiEi/tL5yM3HiROzcuROzZ89GaGgoACAmJgZTpkxB//79MXfuXH3k1BleuSEyrJJyNT756QJ+OPHwL0UvNPfAN6+1hoPKSuRkRGRK9HpbqrS0FFOmTMHSpUtRXv5w4TwrKyuMHTsWs2bNMvonk1huiMSx+WQqPtp1HqVqDfxcbbF8aBAaumt/tZeILJNey83fioqKKp5OatCgAWxsTGOwIMsNkXgS03IxZkM87uQVw1YhwzevtUb3ll5ixyIiE2CQcmOqWG6IxJV1vwTjfziN48kPZzQf+2wDvNu1CWRSzodDRI+n1wHFREQ14WqnxIYRwRgVVh8AsOTgNQxbE4d7haUiJyMic8FyQ0QGJ5dJ8cGLzTF/YCCsrWQ4ciULvRYcxZ+38sSORkRmgOWGiETzUmtv7BwXinq1bHAr9wH6LzmG7fE3xY5FRCZO63JTWFiojxxEZKGaejrgx/Gd0KWJG0rKNXhnayJm7P4TpeUasaMRkYnSutx4eHhg+PDhOHr0qD7yEJEFcrS2wqo32mHS840AAOtib2DwyuPIzC8WORkRmSKty82GDRuQk5OD5557Do0bN8asWbNw+/ZtfWQjIgsilUow+YXGWDm0LeyVcpy8fg+9FhxF/I0csaMRkYnRutz07dsXu3btwq1btzBmzBj88MMPqFevHnr16oUdO3ZUTOxHRPQ0wpt74McJndDYww6ZBSUYsPw41sakwMJmrSCiGtDJPDcLFizAlClTUFpaCldXV4wZMwZTp041yon9OM8NkWkoLCnHe9vPYs/Zh2vL9fL3wqz+/rBTykVORkRiMMg8NxkZGfjqq6/QvHlzTJ06Fa+88gqio6PxzTffYMeOHejbt+/THpqICLZKORYODMRHvZpDLpXg57N38NLCo7icUSB2NCIyclpfudmxYwfWrFmD/fv3o3nz5hg5ciSGDBkCJyenin2uXbuGZs2aobTU+Cbl4pUbItMTfyMH475PQHp+MaytZPiiX0u8HFhH7FhEZEB6vXITGRkJb29vxMTE4MyZMxg/fnylYgMA3t7e+OCDD7Q9NBFRlYLquWDPxE4Ia+SKB2VqTN6ciP/sPIfiMrXY0YjICGl95aaoqMgox9JUF6/cEJkutUbA/OgrmP/7FQgC0Kq2IxYPbgMfF9P9TCKi6tHrlRt7e3tkZmY+sj07OxsymUzbwxERVZvsr8fF10a2h7ONFc7dysOL84/gtwsZYkcjIiOidbl53IWekpISKBSKGgciIvo3zzR2w56JYQis64T84nKMXH8Ks365iHI1ZzUmIqDaz1TOnz8fACCRSLBy5UrY2dlVfE2tVuPw4cNo2rSp7hMSEVXB28kam98MwcxfkrAm5jqWHrqGhNR7WDAoEO72KrHjEZGIqj3mpn79+gCAGzduoE6dOpVuQSkUCvj6+uLTTz9FcHCwfpLqCMfcEJmfn8/exvvbzqKwVA03eyUWDAxEB79aYsciIh3S5ve31gOKu3Tpgh07dsDZ2blGIcXCckNknq7dvY+3NpzGpYwCSCXAlG5NMbqzH6RSidjRiEgH9FpuTB3LDZH5elCqxge7zmHH6VsAgPBm7vjm1QA42liJnIyIakrn5SYqKgqfffYZbG1tERUV9cR9v/32W+3SGhjLDZF5EwQBm06mYcaP51FarkEdZ2ssGRyEVnUcxY5GRDWgze/vag0oTkhIQFlZWcU/P45Ewsu/RCQuiUSCge3rolVtR4z9Ph5pOQ/Qf8kxzHipOQa1r8vPKSILwNtSRGS28orK8M7WRPyW9HAenH6BtfH5yy1ho+Dim0SmxiALZxIRGTtHGyusGBqEaT2aQiaVYEfCLfRdFINrd++LHY2I9KhaV2769etX7QPu2LGjRoH0jVduiCzTieRsjN+YgLsFJbBVyDCrvz96t/YWOxYRVZPOx9w4OnIgHhGZtmC/WtgzsRMmbkzA8eQcTNiYgPgb9/Cfns2gkPMiNpE54ZgbIrIo5WoNvv31MhYfvAYACPBxwsJBgajjzMU3iYwZx9wQET2GXCbFe92bYtUbbeGgkuNMWi56zjuC/efTxY5GRDpSrSs3bdq0QXR0NJydnREYGPjERylPnz6t04C6xis3RPS3tJwiTNiYgDNpuQCAyI6+mNaDt6mIjJHOx9z06dMHSqUSANC3b98aByQiMgY+LjbYMjoEX++/iBVHUrAm5jrib9zDwoFtULcWb1MRmSqOuSEiAvDbhQy8uy0RuUVlsFfK8dUr/ujRykvsWET0F4OsLXXq1CkkJSUBAJo3b46goKCnOYzBsdwQ0ePcyn2AiX89RQUAQ0Pq4T89m0FlJRM5GRHptdzcvHkTAwcORExMDJycnAAAubm5CA0NxaZNm1CnTp2nDm4ILDdE9CRlag2+OXAZSw89fJqqhbcDFg5qg/qutiInI7Jsen1aauTIkSgrK0NSUhJycnKQk5ODpKQkaDQajBw58qlDExEZAyuZFFN7NMWayHZwsVXg/O189F5wFD8m3hY7GhFVk9ZXbqytrXHs2DEEBgZW2h4fH4+wsDAUFRXpNKCu8coNEVVXel4xJm5MQNz1HADAwPZ1MaN3c96mIhKBXq/c+Pj4VKwQ/t/UajW8vTmVORGZD09HFX4YFYzxXRpCIgE2xqVybSoiE6B1ufn6668xYcIEnDp1qmLbqVOnMGnSJMyePVun4YiIxCaXSfFutyZYP7w9atkqcDG9AL0XHMXOhJtiRyOix6jWbSlnZ+dKE/cVFhaivLwccvnDaXL+/mdbW1vk5OToL60O8LYUET2tzPxiTNz0cG0qAHitbR188lJLWCt4m4pI33Q+id/cuXN1kYuIyKS5O6jw/cgOmB99BfN/v4Itp27iTFouFg1qg0Ye9mLHI6K/cBI/IqKncOxqFiZtPoO7BSVQWUnxWZ+WeLWtj9ixiMyWwRbOLC4uRn5+fqUXEZElCG3oir0Tw9CpoSuKyzSYsu0soracQWFJudjRiCye1uWmsLAQ48ePh7u7O2xtbeHs7FzpRURkKdzslVg3vD3eeaExpBJgx+lbeGnhUVxM51/0iMSkdbl577338Pvvv2PJkiVQKpVYuXIlPvnkE3h7e2P9+vX6yEhEZLRkUgkmPN8IP4zqAA8HJa7dLUSfhTHYFJcKC7vrT2Q0tB5zU7duXaxfvx7PPvssHBwccPr0aTRs2BDfffcdNm7ciL179+orq05wzA0R6Uv2/RJEbUnEoct3AQAvtfbGF/1awU5ZrWc3iOgJ9DrmJicnB35+fgAABweHike/O3XqhMOHDz9FXCIi81DLTok1w9rh/e5NIZNK8GPibfSafwTnbuaJHY3Iomhdbvz8/JCSkgIAaNq0KbZs2QIA+OmnnyoW0iQislRSqQRjn22AzW92gJejCtezi9BvSQxWHknmbSoiA9G63ERGRiIxMREAMHXqVCxatAgqlQqTJ0/GlClTdB6QiMgUtfV1wS+TwtC1uQfK1AI+35OE4WtPIut+idjRiMxejee5uX79esW4G39/f13l0huOuSEiQxIEARuO38Bne5JQWq6Bm70ScyMC0LGhq9jRiEyKNr+/OYkfEZEBJN3Jx4SNCbiaeR8SCTD2mQaY/EJjWMlqNN0YkcXQ+yR+0dHR6NWrFxo0aIAGDRqgV69e+O23354qLBGRJWjm5YAfx3fEwPY+EARg8cFreG1ZLNJyisSORmR2tC43ixcvRvfu3WFvb49JkyZh0qRJcHBwQM+ePbFo0SJ9ZCQiMgs2Cjlm9vPHwkGBsFfJkZCai57zj2DP2TtiRyMyK1rflqpTpw6mTp2K8ePHV9q+aNEifPHFF7h165ZOA+oab0sRkTFIyynCxE0JSEjNBQAMaOeDGb1bcIVxosfQ622p3NxcdO/e/ZHtXbt2RV4e53IgIqoOHxcbbBkdgnFdGkAiATadTEPvhUeRdIdLNxDVlNbl5qWXXsLOnTsf2b5792706tVLJ6GIiCyBlUyKKd2aYsOIYLjbK3E18z76LIrB+tjrnBOHqAaqNSf4/PnzK/65efPm+L//+z8cPHgQISEhAIDjx48jJiYG77zzjn5SEhGZsY4NXfHLpDC8uzURf1y6i+m7z+PolSx89Yo/nGwUYscjMjnVGnNTv3796h1MIkFycnKNQ+kTx9wQkbESBAGrY65j1i9JKFML8HJUYW5EAIL9aokdjUh0nOfmCVhuiMjY/XkrDxM2JiAlqxBSCTDx+UaY8FwjyKQSsaMRiUbv89z8TRAE3hcmItKxlrUd8dOETujfpg40AjD3tysYuOI4buc+EDsakUl4qnKzfv16tGrVCtbW1rC2toa/vz++++67pw6xaNEi+Pr6QqVSITg4GHFxcY/dd8eOHWjbti2cnJxga2uLgICAGn1vIiJjZKeU45vXWmNuRABsFTLEpeSg5/wjOHA+XexoREZP63Lz7bffYuzYsejZsye2bNmCLVu2oHv37hgzZgzmzJmjdYDNmzcjKioKM2bMwOnTp9G6dWt069YNmZmZVe7v4uKCDz74ALGxsTh79iwiIyMRGRmJ/fv3a/29iYiMXd/A2tgzMQz+dRyRW1SGN7+Lx/Tdf6K4TC12NCKjpfWYm/r16+OTTz7B0KFDK21ft24dPv74Y6SkpGgVIDg4GO3atcPChQsBABqNBj4+PpgwYQKmTp1arWO0adMGL774Ij777LNHvlZSUoKSkn9W4c3Pz4ePjw/H3BCRSSkt12D2gUtYfvjhQxtNPe2xcFAgGrrbi5yMyDD0Oubmzp07CA0NfWR7aGgo7tzRbgrx0tJSxMfHIzw8/J9AUinCw8MRGxv7r+8XBAHR0dG4dOkSOnfuXOU+M2fOhKOjY8XLx8dHq4xERMZAIZfiPz2bYW1kO7jaKXAxvQC9FhzFxrhUjn0k+h9al5uGDRtiy5Ytj2zfvHkzGjVqpNWxsrKyoFar4eHhUWm7h4cH0tMff185Ly8PdnZ2UCgUePHFF7FgwQK88MILVe47bdo05OXlVbzS0tK0ykhEZEyebeKOvZPCENbIFcVlGkzbcQ5jNsTjXmGp2NGIjEa1JvH7b5988gkiIiJw+PBhdOzYEQAQExOD6OjoKkuPPtjb2+PMmTO4f/8+oqOjERUVBT8/Pzz77LOP7KtUKqFUKg2Si4jIENztVVgX2R4rjybj6/2XsP98BhLTjuDbiNYIbeAqdjwi0Wldbvr374+4uDh8++232LVrFwCgWbNmiIuLQ2BgoFbHcnV1hUwmQ0ZGRqXtGRkZ8PT0fOz7pFIpGjZsCAAICAhAUlISZs6cWWW5ISIyR1KpBG92boDQBq6YuDEByVmFGLzyBEZ3boCoFxpDIa/RTB9EJk2rP/1lZWUYPnw4nJ2dsWHDBsTHxyM+Ph4bNmzQutgAgEKhQFBQEKKjoyu2aTQaREdHVyztUB0ajabSoGEiIkvRsrYjfp7YCQPb+0AQgKWHruGVpceQklUodjQi0WhVbqysrLB9+3adBoiKisKKFSuwbt06JCUlYezYsSgsLERkZCQAYOjQoZg2bVrF/jNnzsSvv/6K5ORkJCUl4ZtvvsF3332HIUOG6DQXEZGpsFHIMbOfP5YOaQNHayucvZmHF+cfwZaTaRxsTBZJ69tSffv2xa5duzB58mSdBIiIiMDdu3cxffp0pKenIyAgAPv27asYZJyamgqp9J8OVlhYiLfeegs3b96EtbU1mjZtig0bNiAiIkIneYiITFX3ll5o7eOEqM2JiE3Oxnvbz+Lg5UzMfNkfjjZWYscjMhit57n5/PPP8c033+D5559HUFAQbG1tK3194sSJOg2oa1xbiojMnVojYPnhZHxz4BLKNQK8HVX4NiIAHbgAJ5kwvS6c+aQVwrkqOBGR8UhMy8WkTQm4nl0EiQR469kGeDu8MaxkHGxMpoergj8Byw0RWZLCknJ88tN5bDl1EwDQ2scJ8wcEoF4t2395J5Fx4argREQEALBVyvHVK62xaFAbOKjkSEzLRc95R7At/iY/v8lsPVW5WbVqFVq2bAmVSgWVSoWWLVti5cqVus5GREQ68qK/F355uzPa13dBYaka725NxISNCch7UCZ2NCKd07rcTJ8+HZMmTULv3r2xdetWbN26Fb1798bkyZMxffp0fWQkIiIdqO1kjY2jOmBKtyaQSSX4+ewd9Jx3BHEpOWJHI9IprcfcuLm5Yf78+Rg4cGCl7Rs3bsSECROQlZWl04C6xjE3RERAQuo9TNp0Bqk5RZBKgPFdGmLi840g52BjMlJ6HXNTVlaGtm3bPrI9KCgI5eXl2h6OiIhEEFjXGXsnhaF/mzrQCMD836/i1WWxSM0uEjsaUY1pXW5ef/11LFmy5JHty5cvx+DBg3USioiI9M9OKcc3r7XGgoGBsFfJkZCai57zj2Bnwk2xoxHViNa3pSZMmID169fDx8cHHTp0AACcOHECqampGDp0KKys/pkF89tvv9VtWh3gbSkiokfdvFeEyZvP4OT1ewCAPgHe+KxvSzioOLMxGQe9znPTpUuXau0nkUjw+++/a3Nog2C5ISKqWrlag8UHr2Fe9BWoNQLqOFtjbkQA2vq6iB2NiJP4PQnLDRHRk8XfuIe3NycgLecBpBJg3F+DjTmzMYnJYJP4ERGR+Qmq54y9E8PQr01taARgwe9X8cqSY0i+e1/saETVwnJDRESPsFdZ4dvXArBwUCAcra2QeDMPL84/ih9OpHJmYzJ6LDdERPRYvfy9se/tMHRsWAsPytT4z85zGLU+Htn3S8SORvRYLDdERPREXo7W+G54MD58sRkUMil+S8pAt7lH8MfFTLGjEVWJ5YaIiP6VVCrByDA/7BrXEY097JB1vwSRa09i+u4/8aBULXY8okpYboiIqNqaezvgx/GdMLxjfQDA+tgb6LXgCP68lSdyMqJ/sNwQEZFWVFYyTO/dHN+NaA93eyWu3S3Ey4tjsOTgNag1HGxM4mO5ISKipxLWyA373+6M7i08UaYW8OW+ixi44jhu3uP6VCQulhsiInpqzrYKLBnSBl+94g9bhQxxKTnoMe8Idp+5JXY0smAsN0REVCMSiQSvtfXB3klhaFPXCQXF5Zi06QwmbkxA3oMyseORBWK5ISIinahXyxZbRodgcnhjyKQS/Jh4Gz3mHkbstWyxo5GFYbkhIiKdkcukmBTeCNvGhMC3lg1u5xVj0MrjmPlLEkrLNWLHIwvBckNERDoXWNcZeyaGYUA7HwgCsOxQMvouisHVzAKxo5EFYLkhIiK9sFXKMau/P5a9HgRnGytcuJOPF+cfxfrY61yfivSK5YaIiPSqWwtP7H+7Mzo3dkNJuQbTd59H5NqTyCwoFjsamSmWGyIi0jt3BxXWRbbDJy+1gFIuxcFLd9F97hHsP58udjQyQyw3RERkEBKJBG+E+uKnCZ3Q3MsBOYWlGP1dPKZsTURBMR8ZJ91huSEiIoNq7GGPneNCMfbZBpBIgK3xN9Fj3hHEpeSIHY3MBMsNEREZnFIuw/vdm2LzmyGo42yNm/ceIGJ5LGb9chEl5VxlnGqG5YaIiETTvr4LfpkUhtfa1oEgAEsPXUPfRcdwKZ2PjNPTY7khIiJR2aus8NUrrbHs9SC42CqQdCcfvRccxcojydBwlXF6Ciw3RERkFP5+ZPz5pu4oVWvw+Z4kDF55ArdyH4gdjUwMyw0RERkNN3slVr7RFjP7tYKNQobY5Gx0n3MYOxNucuI/qjaWGyIiMioSiQQD29fF3olhCKzrhIKSckzenIjxPyTgXmGp2PHIBLDcEBGRUfJ1tcXW0SF4t2tjyKUS7Dl3B93mHsahy3fFjkZGjuWGiIiMllwmxfjnGmHnWx3RwM0WmQUleGN1HKbv/hMPSvnIOFWN5YaIiIxeqzqO2DMxDMNCfQEA62Nv4MX5R5CYlitqLjJOLDdERGQSVFYyfPxSC3w3oj08HJRIzipEvyXHMO+3KyhXa8SOR0aE5YaIiExKWCM37H+7M3r5e0GtETDnt8vovzQWyXfvix2NjATLDRERmRwnGwUWDmqDeQMCYK+SIzEtFy/OP4oNx2/wkXFiuSEiItPVJ6A29r/dGaENauFBmRof7voTkWtPIjO/WOxoJCKWGyIiMmneTtbYMCIYH/VqDoVcioOX7qLb3MP45dwdsaORSFhuiIjI5EmlEozoVB8/T+iEFt4OuFdUhrHfn0bU5jPIe1AmdjwyMJYbIiIyG4097LHzrY4Y16UBpBJgR8ItdJ97GEevZIkdjQyI5YaIiMyKQi7FlG5NsXVMKHxr2eBOXjGGrDqBj388z4n/LATLDRERmaWges7YOykMr3eoBwBYe+w6Xpx/BGc48Z/ZY7khIiKzZaOQ47O+LbFu+D8T//VfcgzfHriE0nJO/GeuWG6IiMjsPdPYDQfefgZ9Aryh1giY//tV9FsSg8sZBWJHIz1guSEiIovgaGOFeQMCsXBQIJxsrPDnrXz0WnAUK48kQ6PhxH/mhOWGiIgsSi9/bxx4uzO6NHFDabkGn+9JwsAVx5GWUyR2NNIRlhsiIrI47g4qrB7WDjP7tYKNQoYTKTnoPvcwtpxM4/INZoDlhoiILJJEIsHA9nXxy6QwtPN1RmGpGu9tP4tR60/hbkGJ2PGoBlhuiIjIotWrZYtNb4ZgWo+mUMik+C0pk8s3mDiWGyIisngyqQSjn2mAHyd0RDMvB+QUlnL5BhPGckNERPSXpp4O2DUuFG89y+UbTBnLDRER0X9RymV4r3tTbB0TgnpcvsEksdwQERFVIaieC36ZFIYhHeoC4PINpsQoys2iRYvg6+sLlUqF4OBgxMXFPXbfFStWICwsDM7OznB2dkZ4ePgT9yciInpaNgo5Pu/biss3mBjRy83mzZsRFRWFGTNm4PTp02jdujW6deuGzMzMKvc/ePAgBg4ciD/++AOxsbHw8fFB165dcevWLQMnJyIiS/FMYzfsf7szXmr9z/INLy+OwaV0Lt9gjCSCyLMVBQcHo127dli4cCEAQKPRwMfHBxMmTMDUqVP/9f1qtRrOzs5YuHAhhg4d+q/75+fnw9HREXl5eXBwcKhxfiIisiw/n72ND3f9idyiMihkUkR1bYxRYX6QSSViRzNr2vz+FvXKTWlpKeLj4xEeHl6xTSqVIjw8HLGxsdU6RlFREcrKyuDi4lLl10tKSpCfn1/pRURE9LT+Xr7huabuKFVrMOuXi3htWSyuZxWKHY3+Imq5ycrKglqthoeHR6XtHh4eSE9Pr9Yx3n//fXh7e1cqSP9t5syZcHR0rHj5+PjUODcREVk2dwcVVr3RFl/194edUo74G/fQY94RrI+9zkU4jYDoY25qYtasWdi0aRN27twJlUpV5T7Tpk1DXl5exSstLc3AKYmIyBxJJBK81s4H+94OQ4hfLTwoU2P67vMYujoOt3IfiB3PoolablxdXSGTyZCRkVFpe0ZGBjw9PZ/43tmzZ2PWrFk4cOAA/P39H7ufUqmEg4NDpRcREZGu1HG2wfcjg/Fx7+ZQWUlx9GoWus85jK2nuAinWEQtNwqFAkFBQYiOjq7YptFoEB0djZCQkMe+76uvvsJnn32Gffv2oW3btoaISkRE9FhSqQTDOtbH3olhCKzrhIKSckzZdhaj1scjs6BY7HgWR/TbUlFRUVixYgXWrVuHpKQkjB07FoWFhYiMjAQADB06FNOmTavY/8svv8RHH32E1atXw9fXF+np6UhPT8f9+/fF+hGIiIgAAH5udtg6OgTvdW8CK5kEvyVloNucw9jLRTgNSvRyExERgdmzZ2P69OkICAjAmTNnsG/fvopBxqmpqbhz558/FEuWLEFpaSleeeUVeHl5Vbxmz54t1o9ARERUQS6T4q1nG+LH8Z3QzMsB94rK8Nb3pzFxYwJyi0rFjmcRRJ/nxtA4zw0RERlKabkGC36/gsUHr0GtEeBur8SX/f3Rpam72NFMjsnMc0NERGTOFHIp3unaBNvHhsLPzRaZBSWIXHsSU7efxf2ScrHjmS2WGyIiIj0L8HHC3olhGN6xPgBg08k0dJ97GLHXskVOZp5YboiIiAxAZSXD9N7NsXFUB9RxtsbNew8wcMVxfPLTeRSXqcWOZ1ZYboiIiAwopEEt7Hu7Mwa2fzhj/pqY6+g5/wgSUu+JnMx8sNwQEREZmJ1Sjpn9/LEmsh3c7ZVIvluI/kuO4ev9F1FarhE7nsljuSEiIhJJlybuODC5M/oEeEMjAIv+uIY+i2KQdIeLPNcEyw0REZGInGwUmDcgEIsHt4GzjRWS7uTjpYVHsfD3KyhX8yrO02C5ISIiMgI9W3nhwORn8EJzD5SpBcw+cBn9lxzD1cwCsaOZHJYbIiIiI+Fmr8Ty14Pw7WutYa+SI/FmHnrOP4oVh5Oh1ljUnLs1wnJDRERkRCQSCfq1qYMDkzujc2M3lJZr8H97kxCxLBbXswrFjmcSWG6IiIiMkJejNdZFtsPMfq1gq5Dh1I176DHvCNYduw4Nr+I8EcsNERGRkZJIJBjYvi72vd0ZIX618KBMjRk/nsfglSeQllMkdjyjxXJDRERk5HxcbPD9yGB82qcFrK1kiE3ORve5h7ExLhUWtv51tbDcEBERmQCpVIKhIb74ZVIY2tZzRmGpGtN2nMOwNSdxJ++B2PGMCssNERGRCfF1tcXm0SH4oGczKORSHLp8F13nHMb2+Ju8ivMXlhsiIiITI5NKMKqzH/ZO7ITWdRxRUFyOd7YmYtT6eGQWFIsdT3QsN0RERCaqobs9to8NxZRuTWAlk+C3pAx0nXMYPyXeFjuaqFhuiIiITJhcJsW4Lg3x04ROaOHtgNyiMkzYmIBx359GTmGp2PFEwXJDRERkBpp6OmDXuI6Y9HwjyKUS7Dl3B13nHML+8+liRzM4lhsiIiIzYSWTYvILjbHzrY5o7GGHrPulGP1dPCZvPoO8ojKx4xkMyw0REZGZaVXHET9N6IQxzzSAVALsTLiFrnMP4Y9LmWJHMwiWGyIiIjOklMswtUdTbBsbCj9XW2TklyByzUm8v+0sCorN+yoOyw0REZEZa1PXGXsmhmFEp/qQSIDNp9LQfe4RHL2SJXY0vWG5ISIiMnPWChk+6tUcm0Z1QF0XG9zKfYAhq07gg53ncL+kXOx4OsdyQ0REZCGC/Wrhl0lhGBpSDwDw/YlUdJ97GMeumddVHJYbIiIiC2KrlOPTPi3xw8hg1Hayxs17DzBoxQnM2P0nikrN4yoOyw0REZEFCm3oiv2TO2NwcF0AwLrYG+g+9whOJGeLnKzmWG6IiIgslJ1Sjv97uRW+G9Ee3o4qpOYUYcCK4/jkp/N4UKoWO95TY7khIiKycGGN3LBvcmcMaOcDQQDWxFxHz/lHcOp6jtjRngrLDREREcFBZYVZ/f2xJrIdPB1USMkqxKvLYvF/ey6guMy0ruKw3BAREVGFLk3csX9yZ7wSVAeCAKw4koKe84/gdOo9saNVG8sNERERVeJobYXZr7bGqjfawt1eieS7hXhlyTHM+uWiSVzFYbkhIiKiKj3fzAMHJnfGy4G1oRGApYeuofeCo0hMyxU72hOx3BAREdFjOdkoMCciAMtfD4KrnRJXMu+j35Jj+Hr/RZSUG+dVHJYbIiIi+lddW3ji18md8VJrb6g1Ahb9cQ0vLYjBn7fyxI72CJYbIiIiqhZnWwXmDwzEksFtUMtWgUsZBeizKAbfHriE0nKN2PEqsNwQERGRVnq08sKByZ3xYisvqDUC5v9+FX0WxeDC7XyxowFguSEiIqKnUMtOiUWD22DhoEA421gh6U4+Xlp4FPN+u4IytbhXcVhuiIiI6Kn18vfGgcnPoFsLD5RrBMz57TJeXhwj6iPjLDdERERUI272SiwdEoR5AwLgaG2FVrWdoLKSiZZHLtp3JiIiIrMhkUjQJ6A2QvxqwVohXrEBWG6IiIhIh9wdVGJH4G0pIiIiMi8sN0RERGRWWG6IiIjIrLDcEBERkVlhuSEiIiKzwnJDREREZoXlhoiIiMwKyw0RERGZFZYbIiIiMissN0RERGRWWG6IiIjIrLDcEBERkVlhuSEiIiKzYnGrgguCAADIz88XOQkRERFV19+/t//+Pf4kFlduCgoKAAA+Pj4iJyEiIiJtFRQUwNHR8Yn7SITqVCAzotFocPv2bdjb20Mikej02Pn5+fDx8UFaWhocHBx0emz6B8+zYfA8GwbPs+HwXBuGvs6zIAgoKCiAt7c3pNInj6qxuCs3UqkUderU0ev3cHBw4H84BsDzbBg8z4bB82w4PNeGoY/z/G9XbP7GAcVERERkVlhuiIiIyKyw3OiQUqnEjBkzoFQqxY5i1nieDYPn2TB4ng2H59owjOE8W9yAYiIiIjJvvHJDREREZoXlhoiIiMwKyw0RERGZFZYbIiIiMissN0+waNEi+Pr6QqVSITg4GHFxcU/cf+vWrWjatClUKhVatWqFvXv3Vvq6IAiYPn06vLy8YG1tjfDwcFy5ckWfP4LJ0PW5HjZsGCQSSaVX9+7d9fkjmARtzvP58+fRv39/+Pr6QiKRYO7cuTU+pqXQ9Xn++OOPH/nz3LRpUz3+BKZBm/O8YsUKhIWFwdnZGc7OzggPD39kf35GV03X59kgn88CVWnTpk2CQqEQVq9eLZw/f14YNWqU4OTkJGRkZFS5f0xMjCCTyYSvvvpKuHDhgvDhhx8KVlZWwrlz5yr2mTVrluDo6Cjs2rVLSExMFF566SWhfv36woMHDwz1YxklfZzrN954Q+jevbtw586dildOTo6hfiSjpO15jouLE959911h48aNgqenpzBnzpwaH9MS6OM8z5gxQ2jRokWlP893797V809i3LQ9z4MGDRIWLVokJCQkCElJScKwYcMER0dH4ebNmxX78DP6Ufo4z4b4fGa5eYz27dsL48aNq/h3tVoteHt7CzNnzqxy/9dee0148cUXK20LDg4WRo8eLQiCIGg0GsHT01P4+uuvK76em5srKJVKYePGjXr4CUyHrs+1IDz8j6dPnz56yWuqtD3P/61evXpV/tKtyTHNlT7O84wZM4TWrVvrMKXpq+mfvfLycsHe3l5Yt26dIAj8jH4cXZ9nQTDM5zNvS1WhtLQU8fHxCA8Pr9gmlUoRHh6O2NjYKt8TGxtbaX8A6NatW8X+KSkpSE9Pr7SPo6MjgoODH3tMS6CPc/23gwcPwt3dHU2aNMHYsWORnZ2t+x/ARDzNeRbjmKZOn+fkypUr8Pb2hp+fHwYPHozU1NSaxjVZujjPRUVFKCsrg4uLCwB+RldFH+f5b/r+fGa5qUJWVhbUajU8PDwqbffw8EB6enqV70lPT3/i/n//rzbHtAT6ONcA0L17d6xfvx7R0dH48ssvcejQIfTo0QNqtVr3P4QJeJrzLMYxTZ2+zklwcDDWrl2Lffv2YcmSJUhJSUFYWBgKCgpqGtkk6eI8v//++/D29q74xc3P6Efp4zwDhvl8trhVwckyDBgwoOKfW7VqBX9/fzRo0AAHDx7E888/L2IyIu316NGj4p/9/f0RHByMevXqYcuWLRgxYoSIyUzTrFmzsGnTJhw8eBAqlUrsOGbrcefZEJ/PvHJTBVdXV8hkMmRkZFTanpGRAU9Pzyrf4+np+cT9//5fbY5pCfRxrqvi5+cHV1dXXL16teahTdDTnGcxjmnqDHVOnJyc0LhxY/55forzPHv2bMyaNQsHDhyAv79/xXZ+Rj9KH+e5Kvr4fGa5qYJCoUBQUBCio6Mrtmk0GkRHRyMkJKTK94SEhFTaHwB+/fXXiv3r168PT0/PSvvk5+fjxIkTjz2mJdDHua7KzZs3kZ2dDS8vL90ENzFPc57FOKapM9Q5uX//Pq5du8Y/z1qe56+++gqfffYZ9u3bh7Zt21b6Gj+jH6WP81wVvXw+63W4sgnbtGmToFQqhbVr1woXLlwQ3nzzTcHJyUlIT08XBEEQXn/9dWHq1KkV+8fExAhyuVyYPXu2kJSUJMyYMaPKR8GdnJyE3bt3C2fPnhX69Olj8Y8ZCoLuz3VBQYHw7rvvCrGxsUJKSorw22+/CW3atBEaNWokFBcXi/IzGgNtz3NJSYmQkJAgJCQkCF5eXsK7774rJCQkCFeuXKn2MS2RPs7zO++8Ixw8eFBISUkRYmJihPDwcMHV1VXIzMw0+M9nLLQ9z7NmzRIUCoWwbdu2So8gFxQUVNqHn9GV6fo8G+rzmeXmCRYsWCDUrVtXUCgUQvv27YXjx49XfO2ZZ54R3njjjUr7b9myRWjcuLGgUCiEFi1aCHv27Kn0dY1GI3z00UeCh4eHoFQqheeff164dOmSIX4Uo6fLc11UVCR07dpVcHNzE6ysrIR69eoJo0aNsuhfuH/T5jynpKQIAB55PfPMM9U+pqXS9XmOiIgQvLy8BIVCIdSuXVuIiIgQrl69asCfyDhpc57r1atX5XmeMWNGxT78jK6aLs+zoT6fJYIgCLq7DkREREQkLo65ISIiIrPCckNERERmheWGiIiIzArLDREREZkVlhsiIiIyKyw3REREZFZYboiIiMissNwQERGRWWG5ISIiIrPCckNERERmheWGiIiIzArLDRGZhbi4OISFhcHe3h62trZo1aoVTp48KXYsIhKBXOwARES6MGDAAISGhmL58uVQqVRITk6Gh4eH2LGISAQsN0RkFsrLy1G3bl00bNgQVlZWqF+/vtiRiEgkEkEQBLFDEBHV1KlTp/Dyyy/jzp07UKlUuHXrFhwdHcWORUQiYLkhIrPwwgsvwN7eHlOnToWLiwv8/PwglXJYIZElYrkhIpOXlZUFNzc3nDlzBq1btxY7DhGJjOWGiMxC3bp1ERgYiOnTp8PV1RUpKSkoLS1F165dxY5GRAbGa7ZEZBZ++eUXaDQadOvWDY0bN8aoUaOQkZEhdiwiEgGv3BAREZFZ4ZUbIiIiMissN0RERGRWWG6IiIjIrLDcEBERkVlhuSEiIiKzwnJDREREZoXlhoiIiMwKyw0RERGZFZYbIiIiMissN0RERGRWWG6IiIjIrPw/UvkjSVQlF5EAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "import numpy as np\n", "from scipy.sparse import coo_matrix\n", "from scipy.sparse.linalg import spsolve\n", "from scipy.linalg import solve\n", "from scipy.optimize import bisect\n", "\n", "def solve(N, eps):\n", " pN, pS, pE, pW = 0.25, 0.25, 0.25+eps, 0.25-eps\n", " pts = [(x,y) for x in range(-N,N+1) for y in range(-N,N+1)]\n", " # we form a linear equation A x = 0, where x = [q(x_i, y_i)]\n", " matrix = [] # entries are tuples (i,j,A_ij)\n", " for i, pt in enumerate(pts):\n", " if pt == (0,0):\n", " matrix.append([i,i,1])\n", " else:\n", " x,y = pt\n", " entries = []\n", " matrix.append((i,i,1))\n", " for neighbor_pt, prob in zip([(x+1, y), (x-1,y), (x,y+1), (x,y-1)], [pE, pW, pN, pS]):\n", " if neighbor_pt in pts: matrix.append([i, pts.index(neighbor_pt), -prob])\n", " matrix = np.array(matrix)\n", " row, col, data = matrix[:, 0].astype(int), matrix[:, 1].astype(int), matrix[:, 2]\n", " A = coo_matrix((data, (row,col))).tocsc()\n", " b = np.array([0]*A.shape[0])\n", " b[pts.index((0,0))] = 1\n", " q = spsolve(A,b)\n", " p = pE*q[pts.index((1,0))] + pW*q[pts.index((-1,0))] + pN*q[pts.index((0,1))] + pS*q[pts.index((0,-1))]\n", " return p\n", "\n", "F = np.vectorize(lambda t: solve(10,t))\n", "X = np.linspace(0, 0.25, 20)\n", "Y = F(X)\n", "\n", "import matplotlib.pyplot as plt\n", "\n", "plt.plot(X, Y)\n", "plt.xlabel(\"ε\"); plt.ylabel(\"probability of return\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With this, we get a rough estimate of $ p \\approx 0.06 $" ] }, { "cell_type": "markdown", "metadata": { "id": "ddnj2lNgxiqo" }, "source": [ "### Attempt 2: Combinatorics\n", "We define $ E $ to be the flea's expected number of visits to the origin. Note that \n", "$$ E = \\sum_{k>0} k P(\\text{flea returns to origin in k steps}) = \\sum_{k=1}^\\infty k p^{k-1} (1-p) = \\frac{1}{1-p}$$\n", "\n", "\n", "Next, note that any walk of the flea that returns to the origin has to take an even number of steps. Note also that any such walk of $2k$ steps must have an equal number of steps going north and south as well as an equal number of steps going east and west. Therefore, defining $p_{2k}$ to be the probability that the flea returns to the origin after $2k$ steps, we get\n", "\n", "$$p_{2k} = \\sum_{j=0}^{k}\\binom{2k}{j,j,k-j,k-j} p_{N}^j p_{S}^j p_{E}^{k-j} p_{W}^{k-j} = \\sum_{j=0}^k \\binom{2k}{k} \\binom{k}{j}^2 (p_{E} p_{W})^{k-j} (p_{N} p_{S} )^{j} $$\n", "\n", "Finally, note that we can write $$\\begin{align}\n", "E &= \\sum_{k>0} E(\\text{flea returns to origin after 2k steps}) \\\\\n", "&= \\sum_{k=1}^\\infty p_{2k} \\\\\n", "&= \\sum_{k=1}^\\infty \\sum_{j=0}^k \\binom{2k}{k} \\binom{k}{j}^2 (p_{E} p_{W})^{k-j} (p_{N} p_{S} )^{j}\n", "\\end{align}$$. \n", "\n", "Note we have that $p = \\frac{1}{2} \\implies E = \\frac{1}{1-\\frac{1}{2}} = 2$, so we can use standard root-finding techniques to find an $\\epsilon$ such that $E = 2$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": false }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAABLLklEQVR4nO3deXxM5+IG8GdmMtmTicgmkoiIiERsQWKnVftWWqqtfS9Kafuji1bvvUUXWqr2pah9CaWlaheCkNiykIUEWUhk3zPn94fKbS4iEzNzZnm+n898rkzOnDxzbmoe57znfSWCIAggIiIiMhBSsQMQERERqRPLDRERERkUlhsiIiIyKCw3REREZFBYboiIiMigsNwQERGRQWG5ISIiIoNiInYAbVMqlbh//z5sbGwgkUjEjkNERETVIAgCcnNz4erqCqm06nMzRldu7t+/D3d3d7FjEBERUQ0kJyfDzc2tym2MrtzY2NgAeHxwbG1tRU5DRERE1ZGTkwN3d/eKz/GqGF25eXIpytbWluWGiIhIz1RnSAkHFBMREZFBYbkhIiIig8JyQ0RERAaF5YaIiIgMCssNERERGRSWGyIiIjIoLDdERERkUFhuiIiIyKCw3BAREZFBYbkhIiIig8JyQ0RERAaF5YaIiIgMCsuNGt1My8Xth/lixyAiIjJqLDdqcv1eNoasPId31pxHSnah2HGIiIiMFsuNmjjbmqOWpSnuZRXi3TXnkZFXLHYkIiIio8RyoyaONmbYPC4IrgpzxD/Ix4h1F5BdWCp2LCIiIqPDcqNGde0ssHlcEBysTXHjfg7GbriIwpJysWMREREZFZYbNfNytMbGMUGwMTdB+J1HmLj5EorLWHCIiIi0heVGA/xcbbFhdGtYyGU4dfMBpm+NRFm5UuxYRERERoHlRkMC69lj9YhWMJVJcehGKj7edRVKpSB2LCIiIoPHcqNBHRo64Ke3W0AmlWBPxD18sf8GBIEFh4iISJNYbjSsu78LFg1pBokE2BR2BwsPxbLgEBERaRDLjRYMaF4X/xkYAABYcTIey47HiZyIiIjIcLHcaMnbQR74rE9jAMB3f97E2jOJIiciIiIyTCw3WjSuoxdmdGsIAPjXgShsvZAkciIiIiLDw3KjZdNfbYiJnbwAAJ/svYaQiHsiJyIiIjIsLDdaJpFIMLuXL4YH14MgALN2XsGh66lixyIiIjIYLDcikEgkmNffH28EuqFcKWDa1ss4HpsudiwiIiKDwHIjEqlUgoWDm6JP0zooLRcwadMlnI17KHYsIiIivcdyIyKZVIIfhjZHt8bOKC5TYuwv4bh4O1PsWERERHqN5UZkcpkUy95pgU4+jigsLcfo9RcRmZwldiwiIiK9xXKjA8xMZFj5biCCveyRV1yGEWvP48b9bLFjERER6SWWGx1hYSrD2pGtEVivFnKKyjB87QXEpuaKHYuIiEjvsNzoECszE6wf3RpN3RTIzC/BO2vOI/5BntixiIiI9ArLjY6xNZdj45g28Ktji4d5xXh7dRhuP8wXOxYREZHeYLnRQXaWptg8LgiNnG2QlvO44CRnFogdi4iISC+IWm7mz5+P1q1bw8bGBk5OThg4cCBiY2OrfM2NGzcwePBgeHp6QiKR4IcfftBOWC2zt3pccBo4WuF+dhGGrQ7D/axCsWMRERHpPFHLzcmTJzFlyhSEhYXhyJEjKC0tRffu3ZGf//zLMAUFBfDy8sKCBQvg4uKixbTa52hjhi3jg+FZ2xJ3HxXi7dVhSM0uEjsWERGRTpMIgiCIHeKJBw8ewMnJCSdPnkSnTp1euL2npydmzJiBGTNmPHeb4uJiFBcXV3ydk5MDd3d3ZGdnw9bWVh2xNe5+ViGGrjqH5MxCeDlYYduEYDjZmosdi4iISGtycnKgUCiq9fmtU2NusrMfz+1ib2+vtn3Onz8fCoWi4uHu7q62fWuLq50Fto4PRl07CyQ8zMfba87jQW7xi19IRERkhHSm3CiVSsyYMQPt27dHkyZN1LbfOXPmIDs7u+KRnJystn1rk1stS2wdH4w6CnPEpefhnTVhyMhjwSEiIvpfOlNupkyZguvXr2Pbtm1q3a+ZmRlsbW0rPfSVR+3HBcfZ1gw30/LwzprzeJRfInYsIiIinaIT5Wbq1Kk4cOAAjh8/Djc3N7Hj6DRPBytsGR8MRxszxKTmsuAQERH9D1HLjSAImDp1Kvbu3Ytjx46hfv36YsbRGw0crbF1fBAcrM0QlZKDd9eeR1YBCw4REREgcrmZMmUKNm/ejC1btsDGxgapqalITU1FYeF/53MZMWIE5syZU/F1SUkJIiMjERkZiZKSEty7dw+RkZGIi4sT4y2IxtvJ5u+CY4ob93MwfO0FZBeUih2LiIhIdKLeCi6RSJ75/Pr16zFq1CgAQJcuXeDp6YkNGzYAAG7fvv3MMzydO3fGiRMnXvgzVbmVTB/Epubi7dVhyMgvQVM3BTaNDYLCQi52LCIiIrVS5fNbp+a50QZDKzcAEJOag7dXn0dmfgmauSmwkQWHiIgMjN7Oc0M14+tii1/HBaGWpRxX7mZjxLoLyC7kJSoiIjJOLDcGonEdW/w6LvhxwUnOYsEhIiKjxXJjQPxcWXCIiIhYbgzMk4Jjx4JDRERGiuXGAPm52mLLPwvO2vMsOEREZDRYbgzUk4LzZJDx8LXnOQ8OEREZBZYbA+bnaost44Nhb2WKq3ez8c7aMM5kTEREBo/lxsA1rmOLreODUdvKFNfvPZ4Ph2tRERGRIWO5MQKNXGywdUIwHKxNEZWSg7fXPJ7wj4iIyBCx3BgJH2cbbB0fDAdrM0Sn5ODt1WF4mFcsdiwiIiK1Y7kxIg2dbbBtQjCcbMwQk5qLYavCkJ5bJHYsIiIitWK5MTLeTtbYPrEtXGzNcSs9D2+tCkNaDgsOEREZDpYbI1TfwQrbJwbDVWGOhAf5eGtVGFKyC8WORUREpBYsN0aqXm0rbJ/YFnXtLJD4MB9DV4bh7qMCsWMRERG9NJYbI+Zub4ntE4PhYW+JpMwCDF0ZhqQMFhwiItJvLDdGzq2WJXZMbAsvByvcyyrE0FXnkPgwX+xYRERENcZyQ3BRmGPbhGA0dLJGSnYRhq48h7j0XLFjERER1QjLDQEAnGzNsXVCMHxdbJCeW4y3VoUhJjVH7FhEREQqY7mhCg7WZtg6PhhN6triYV4J3loVhmt3s8WORUREpBKWG6qklpUpfh0XjObudsgqKMXba8Jw6c4jsWMRERFVG8sNPUVhIcfmcUFo42mP3KIyjFh7HucTMsSORUREVC0sN/RM1mYm2DCmNTp4OyC/pBwj11/A6VsPxI5FRET0Qiw39FyWpiZYM7IVujZyRFGpEmM3hOOvqDSxYxEREVWJ5YaqZC6XYeXwVujp74KSciUmbb6E367cFzsWERHRc7Hc0AuZmkjx09stMLC5K8qUAqZvi8DO8GSxYxERET0Tyw1Vi4lMiu+HNMewNu5QCsBHu65i07nbYsciIiJ6CssNVZtMKsHXrwdgdHtPAMDn+25g5cl4cUMRERH9D5YbUolEIsHcvn6Y0rUBAGD+HzFY9GcsBEEQORkREdFjLDekMolEgo96+OLjno0AAEuOxeFfB6JZcIiISCew3FCNvdfFG/P6+wMA1oUmYs6eayhXsuAQEZG4WG7opYxs54lv32gKqQTYdjEZM7ZHorRcKXYsIiIyYiw39NLebOWOpcNawkQqwW9X7mPipksoKi0XOxYRERkplhtSiz5N62D1iFYwM5HiWEw6Rq2/gLziMrFjERGREWK5IbXp6uuEX8a0gbWZCcISMvHO6jA8yi8ROxYRERkZlhtSq2Cv2tgyPgh2lnJcuZuNt1aFIT2nSOxYRERkRFhuSO2autlhx8S2cLIxQ2xaLt5ceQ7JmQVixyIiIiPBckMa4eNsg12T2sHD3hJ3MgowePlZ3EzLFTsWEREZAZYb0hiP2pbYOaktfJytkZ5bjCErz+FKcpbYsYiIyMCx3JBGOduaY/uEtmjmboesglK8vToMZ+Mfih2LiIgMGMsNaVwtK1P8Oi4I7RrURn5JOUatv4jDN1LFjkVERAaK5Ya0wtrMBOtGtUZ3P2eUlCkxefMl7AxPFjsWEREZIJYb0hpzuQw/v9MSbwa6QSkAH+26ijWnE8SORUREBoblhrTKRCbFN280xfiO9QEA/z4YjW8Px3BFcSIiUhuWG9I6iUSCT3o3xkc9GgEAlh2Pxyd7r3NFcSIiUguWGxKFRCLBlK7e+Pr1AEgkwNYLSZi29TKKy7jgJhERvRyWGxLV20EeWPZ2S5jKpPj9WirGbLjIBTeJiOilsNyQ6HoH1MH60a1hZSpDaFwG3l4dhoy8YrFjERGRnmK5IZ3Q3tsBWycEw97KFFfvZuPNFedw9xHXoyIiItWx3JDOaOpmh52T2qKunQUSHuZj8PKziE3lelRERKQalhvSKQ0crbF7cjv4OFsjLacYb644i4u3M8WORUREeoTlhnSOi8IcOya2RWC9WsgpKsO7a87jr6g0sWMREZGeYLkhnWRnaYrNY4Pwqq8TisuUmLj5EnZc5HINRET0Yiw3pLMsTGVYMTwQbwS6oVwp4OPdV7HseBxnMyYioiqx3JBOk8uk+PaNppjcpQEA4NvDsfhy/w3OZkxERM/FckM6TyKR4P96+mJuXz8AwC/n7uD9rRGczZiIiJ6J5Yb0xpgO9bFkWAvIZRIcvJaCUesuIqeoVOxYRESkY1huSK/0b+aKDaPbwNrMBOcSMjBkxTmk5RSJHYuIiHQIyw3pnfbeDtg2IRgO1maISc3FoJ/PIi49T+xYRESkI1huSC81qavA3vfaob6DFe5lFeKNFWdx6c4jsWMREZEOYLkhveVub4ldk9qimbsdsgpK8fbqMBzhZH9EREaP5Yb0Wm1rM2wdH4SujRwfT/a3KRy/nr8jdiwiIhIRyw3pPUtTE6we0QpDW7lDKQCf7r2O7w7HcrI/IiIjxXJDBsFEJsWCwQGY0a0hAOCn43H4cOdVlJYrRU5GRETaxnJDBkMikWBGNx8sHBwAmVSC3ZfvYsyGi8grLhM7GhERaRHLDRmcoa09sGZEK1jIZTh96yGGrDiHdM6FQ0RkNFhuyCB19XXC9onBcLA2RVRKDl7/+SxupeWKHYuIiLSA5YYMVlM3O+yZ3B5ef8+FM3j5WYQlZIgdi4iINKxG5UapVOLmzZs4c+YMTp06Vemhivnz56N169awsbGBk5MTBg4ciNjY2Be+bufOnfD19YW5uTkCAgLw+++/1+RtkBHwqG2J3ZPbIbBeLeQUlWHE2gvYf+W+2LGIiEiDVC43YWFh8Pb2RuPGjdGpUyd06dKl4tG1a1eV9nXy5ElMmTIFYWFhOHLkCEpLS9G9e3fk5+c/9zVnz57FsGHDMHbsWERERGDgwIEYOHAgrl+/rupbISNRy8oUv44LQq8mLigpV+L9rRFYfiKet4oTERkoiaDi3/DNmzeHj48P5s2bhzp16kAikVT6vkKhqHGYBw8ewMnJCSdPnkSnTp2euc3QoUORn5+PAwcOVDwXHByM5s2bY8WKFS/8GTk5OVAoFMjOzoatrW2Ns5L+KVcK+M/BaKwLTQQAvBPkgXn9/WEi49VZIiJdp8rnt4mqO7916xZ27doFb2/vGgd8nuzsbACAvb39c7c5d+4cZs6cWem5Hj16ICQk5JnbFxcXo7i4uOLrnJyclw9KekkmlWBuPz+421vgqwNR+PV8ElKyi7B0WAtYman8nwIREekolf/JGhQUhLi4OLUHUSqVmDFjBtq3b48mTZo8d7vU1FQ4OztXes7Z2RmpqanP3H7+/PlQKBQVD3d3d7XmJv0zun19LH8nEGYmUhyLScfQVbxVnIjIkKj8z9Vp06Zh1qxZSE1NRUBAAORyeaXvN23atEZBpkyZguvXr+PMmTM1ev3zzJkzp9KZnpycHBYcQs8mLtg6IRjjfgnH9XuPbxVfN6o1GrnYiB2NiIheksrlZvDgwQCAMWPGVDwnkUggCAIkEgnKy8tVDjF16lQcOHAAp06dgpubW5Xburi4IC2t8srPaWlpcHFxeeb2ZmZmMDMzUzkTGb6WHrWw9712GL3+IhIe5uON5Wex/N1AdGjoIHY0IiJ6CSpflkpMTHzqkZCQUPG/qhAEAVOnTsXevXtx7Ngx1K9f/4Wvadu2LY4ePVrpuSNHjqBt27Yq/WwiAKhX2wp73muHNp72yC0uw6j1F7AjPFnsWERE9BJUuluqtLQUvr6+OHDgABo3bvzSP/y9997Dli1bsG/fPjRq1KjieYVCAQsLCwDAiBEjULduXcyfPx/A41vBO3fujAULFqBPnz7Ytm0bvv76a1y+fLnKsTpP8G4pepbisnJ8vOsq9kU+ngNnaldvzOru89TdgEREJA5VPr9VOnMjl8tRVKS+gZfLly9HdnY2unTpgjp16lQ8tm/fXrFNUlISUlJSKr5u164dtmzZglWrVqFZs2bYtWsXQkJCqlVsiJ7HzESGH4Y2x7RXHt8F+NPxOEzfFomiUtUvsxIRkbhUnufm66+/xs2bN7FmzRqYmOjf7bM8c0MvsiM8GZ/suYYypYBW9Wph1YhWsLcyFTsWEZFRU+XzW+Vy8/rrr+Po0aOwtrZGQEAArKysKn1/z549qifWIpYbqo7QuIeYtPkScovKUK+2JdaPag0vR2uxYxERGS2NTuJnZ2dXcccUkaFq7+2APZPbYfSGi7iTUYDXfz6LlcMDEexVW+xoRET0AiqfudF3PHNDqniQW4zxG8MRmZwFuUyC+YOa4o3AqqcrICIi9dPYgGIiY+NoY4ZtE4LRJ6AOSssFfLjzCr47HAul0qj+TUBEpFdUvixVv379Km+PVXWuGyJdZy6XYemwFvB0sMSy4/H46Xgcbmfk47s3m8FcLhM7HhER/Q+Vy82MGTMqfV1aWoqIiAgcOnQIH330kbpyEekUqVSCj3r4wrO2FT7Zew0HrqbgXlYhVg1vBUcbzoBNRKRL1DbmZtmyZQgPD8f69evVsTuN4Zgbelnn4jMwafMlZBeWoq6dBdekIiLSAlHG3PTq1Qu7d+9W1+6IdFbbBrWx9712qO9ghXtZhRi8/CyOx6SLHYuIiP6mtnKza9cu2Nvbq2t3RDrNy9Eae99rh2Ave+QVl2HsLxexPjQRRnbzIRGRTlJ5zE2LFi0qDSgWBAGpqal48OABfv75Z7WGI9Jldpam2DgmCJ+FXMOO8LuY91sU4h/k4Yt+/pDLeCMiEZFYVC43AwYMqFRupFIpHB0d0aVLF/j6+qo1HJGuMzWRYuHgpmjgaI0Fh2KwOSwJtx8WYNnbLaGwlIsdj4jIKHESPyI1+fNGKmZsj0RBSTm8HK2wdmRr1HewevELiYjohTQ6oFgmkyE9/enBkxkZGZDJOOcHGa/u/i7YOaktXBXmSHiQj4HLQnE2/qHYsYiIjI7K5eZ5J3qKi4thasqVk8m4+bsqEDK1PZq52yG7sBQj1l7AlvNJYsciIjIq1R5zs2TJEgCARCLBmjVrYG393xWSy8vLcerUKY65IQLgZGOO7ROC8dGuq/jtyn18svcabqbl4rM+jWHCgcZERBpX7TE39evXBwDcuXMHbm5ulS5BmZqawtPTE1999RWCgoI0k1RNOOaGtEUQBCw9FodFR24CADr5OGLpsBZQWHCgMRGRqlT5/FZ5QHHXrl2xZ88e1KpV66VCioXlhrTtj2spmLnjCgpLOdCYiKimNDqg+Pjx46hVqxZKSkoQGxuLsrKyGgclMga9Aupg56S2cLH970Dj0DgONCYi0hSVy01hYSHGjh0LS0tL+Pv7Iynp8WDJadOmYcGCBWoPSGQImtRVYP/U9mj+ZKDxugvYeO42ZzQmItIAlcvN7NmzceXKFZw4cQLm5uYVz3fr1g3bt29XazgiQ+Jka45tE4IxqEVdlCsFzN13A5+GXEdpuVLsaEREBkXlchMSEoKffvoJHTp0qDRTsb+/P+Lj49UajsjQmMtl+H5IM8zu5QuJBNhyPgnvrjmPzPwSsaMRERkMlcvNgwcP4OTk9NTz+fn5lcoOET2bRCLBpM4NsGZEK1iZynA+MRP9fzqD6JQcsaMRERkElctNq1atcPDgwYqvnxSaNWvWoG3btupLRmTgXm3sjL1T2qNebUvcfVSIwcvP4vCNVLFjERHpPZUXzvz666/Rq1cvREVFoaysDD/++COioqJw9uxZnDx5UhMZiQyWj7MN9k1pjylbLiM0LgMTN13CrNd8MPUVb54JJSKqIZXP3HTo0AFXrlxBWVkZAgIC8Oeff8LJyQnnzp1DYGCgJjISGTQ7S1P8MroNRrXzBAB8f+Qmpmy5jIISTrNARFQTKk3iV1paiokTJ+Lzzz+vmLFY33ASP9Jl2y8m4bOQ6ygtF+DrYoPVI1rB3d5S7FhERKLT2CR+crkcu3fvfqlwRPR8Q1t7YOv4YDhYmyEmNRf9fzqDc/EZYsciItIrKl+WGjhwIEJCQjQQhYgAoJWnPfZPbY+Augo8KijFu2vPc8I/IiIVqDyguGHDhvjqq68QGhqKwMBAWFlVXiPn/fffV1s4ImPlameBnZPaYvbuqwiJvI+5+27gxr0cfDXQH2YmshfvgIjIiKm8cGZVY20kEgkSEhJeOpQmccwN6RNBELD6dAIW/BEDpQC09LDDincD4WRr/uIXExEZEI2uCq7vWG5IH528+QDTtlxGTlEZnGzMsGJ4IFp61BI7FhGR1mh0VXAi0r7OPo7YP7UDfJytkZ5bjLdWhmHHxWSxYxER6SSWGyI94elghT3vtUcPf2eUlCvx8e6r+DzkOkrKuPAmEdE/sdwQ6RFrMxMsfycQM1/zgUQCbAq7g3fWhCE9t0jsaEREOoPlhkjPSKUSvP9qQ6wZ0Qo2Zia4ePsR+i8NRUTSI7GjERHphGqVm0GDBiEn5/GKxRs3bkRxcbFGQxHRi73a2Bn7praHt5M1UnOKMHRlGLZfTBI7FhGR6Kp1t5SpqSnu3LmDOnXqQCaTISUlBU5OTtrIp3a8W4oMTW5RKWbtuII/o9IAAG8HeeCLfn6cD4eIDIoqn9/VmsTP19cXc+bMQdeuXSEIAnbs2PHcHY8YMUL1xERUYzbmcqx4NxDLjsdh0V83seV8EmJScrD83UA4cz4cIjJC1Tpzc/bsWcycORPx8fHIzMyEjY0NJBLJ0zuTSJCZmamRoOrCMzdkyI7HpGP6tgjkFJXB0cYMy99piVae9mLHIiJ6aRqdxE8qlSI1NZWXpYh01O2H+Zi46RJi03JhIpVgbj8/DA+u98x/kBAR6QuNTuKXmJgIR0fHGocjIs16PB9OO/RpWgdlSgFz993ArJ1XUFRaLnY0IiKtqNHyC1lZWVi7di2io6MBAH5+fhg7diwUCoXaA6obz9yQsRAEAWtOJ2LBoRiUKwX4u9pixbuBcLe3FDsaEZHKNHrmJjw8HA0aNMDixYuRmZmJzMxMLF68GA0aNMDly5drHJqI1EsikWB8Jy9sGtsGta1MceN+Dvr9dAYnbz4QOxoRkUapfOamY8eO8Pb2xurVq2Fi8vhmq7KyMowbNw4JCQk4deqURoKqC8/ckDG6n1WIyb9expXkLEgkwMxuPpjS1RtSKcfhEJF+0OiAYgsLC0RERMDX17fS81FRUWjVqhUKCgpUT6xFLDdkrIrLyvHl/ihsvfB4or9ujZ3w/ZDmUFjIRU5GRPRiGr0sZWtri6Skp2dBTU5Oho2Njaq7IyItMTORYf6gAHwzuClMTaT4Kzod/X86g+iUHLGjERGplcrlZujQoRg7diy2b9+O5ORkJCcnY9u2bRg3bhyGDRumiYxEpEZDWrtj96R2qGtngTsZBXj951DsjbgrdiwiIrVR+bJUSUkJPvroI6xYsQJlZWUAALlcjsmTJ2PBggUwMzPTSFB14WUposce5Zfg/W0ROH3rIQDg3WAPfN6XyzYQkW7S6JibJwoKChAfHw8AaNCgASwt9eP2UpYbov8qVwr48egtLDl6CwDQ3N0OP7/TEq52FiInIyKqTCvlRl+x3BA97XhMOmZsj0R2YSnsrUyx5K0W6NDQQexYREQVNDqgmIgMT1dfJxyY1gH+rrbIzC/B8HXn8dOxW1AqjerfPkRkIFhuiAgA4G5vid2T2+Gt1u4QBOC7P29i3MZwZBWUiB2NiEglLDdEVMFcLsOCwU3xzRtNYWYixbGYdPRdegZX72aJHY2IqNpULjf5+fmayEFEOmRIK3fsea8dPOwtcfdRId5Yfg6/nr8DIxuiR0R6SuVy4+zsjDFjxuDMmTOayENEOsLfVYHfpnXAa37OKClX4tO91zFzxxUUlJSJHY2IqEoql5vNmzcjMzMTr7zyCnx8fLBgwQLcv39fE9mISGQKCzlWDQ/EnF6+kEkl2BtxDwOXhSIuPU/saEREz1XjW8EfPHiATZs2YcOGDYiOjkaPHj0wZswY9O/fv2JBTV3EW8GJauZ8QgambY1Aem4xrExlmD+4Kfo3cxU7FhEZCa3Pc7N06VJ89NFHKCkpgYODAyZNmoTZs2fr5MR+LDdENZeeW4TpWyNxLiEDADA8uB4+69uYsxoTkcZpZZ6btLQ0fPPNN/Dz88Ps2bPxxhtv4OjRo/j++++xZ88eDBw4sKa7JiId5WRjjk1j22BqV28AwKawO3hj+TkkZxaInIyI6L9UPnOzZ88erF+/HocPH4afnx/GjRuHd999F3Z2dhXbxMfHo3Hjxigp0b35MXjmhkg9jsemY+b2SDwqKIWtuQm+e7MZuvu7iB2LiAyURs/cjB49Gq6urggNDUVkZCSmTp1aqdgAgKurKz799FNVd01EeqRrIyccfL8jWnjYIaeoDBM2XcJ/DkahtFwpdjQiMnIqn7kpKCjQybE01cUzN0TqVVKmxMJDMVh7JhEA0NLDDj+9zcU3iUi9NHrmxsbGBunp6U89n5GRAZmMgwqJjI2piRSf9/XDyuGBsDE3weWkLPRechrHY57+e4KISBtULjfPO9FTXFwMU1PTlw5ERPqph78LDk7riIC6CmQVlGL0hotY8EcMyniZioi0rNoT0ixZsgQAIJFIsGbNGlhbW1d8r7y8HKdOnYKvr6/6ExKR3vCobYldk9vi64PR+OXcHaw4GY/w25lYMqwFL1MRkdZUe8xN/fr1AQB37tyBm5tbpUtQpqam8PT0xFdffYWgoCDNJFUTjrkh0o7fr6Xg/3ZdRW5xGWpZyrFoSHN09XUSOxYR6SmNTuLXtWtX7NmzB7Vq1XqpkGJhuSHSnjsZ+Zi6JQLX7mUDACZ28sKHPRpBLqvxFFtEZKQ0OqD4+PHjais2p06dQr9+/eDq6gqJRIKQkJAXvmbZsmVo3LgxLCws0KhRI2zcuFEtWYhI/erVtsKuyW0xqp0nAGDlqQQMXXkOdx9x0j8i0pxqjbmZOXMm/vWvf8HKygozZ86scttFixZV+4fn5+ejWbNmGDNmDAYNGvTC7ZcvX445c+Zg9erVaN26NS5cuIDx48ejVq1a6NevX7V/LhFpj5mJDF/290dQfXt8vPsqLidloc+SM/j2jaac9I+INKJa5SYiIgKlpaUVf34eiUSi0g/v1asXevXqVe3tN23ahIkTJ2Lo0KEAAC8vL1y8eBELFy5kuSHScb0C6sDfVYFpWy/jyt1sTNh0CaPbe2J2L1+uTUVEalWtcnP8+PFn/lnbiouLYW5uXuk5CwsLXLhwAaWlpZDL5c98TXFxccXXOTk5Gs9JRM/mUdsSOye1wzeHYrDmTCLWh95G+O1H+OntFqhX20rseERkIPRqVF+PHj2wZs0aXLp0CYIgIDw8HGvWrEFpaSkePnz4zNfMnz8fCoWi4uHu7q7l1ET0T6YmUnzW1w9rRrSCnaUc1+5lo8+SM9h/5b7Y0YjIQFTrbqnqjId5Ys+ePTULIpFg7969Va4mXlhYiClTpmDTpk0QBAHOzs5499138c033yA1NRXOzs5PveZZZ27c3d15txSRDrifVYjp2yJw8fYjAMBbrd3xRT9/WJjyMhURVab2u6X+eebjRQ9NsrCwwLp161BQUIDbt28jKSkJnp6esLGxgaOj4zNfY2ZmBltb20oPItINrnYW2Do+GNNe8YZEAmy7mIwBy87gZlqu2NGISI9Va8zN+vXrNZ1DJXK5HG5ubgCAbdu2oW/fvpBK9eoKGxH9zUQmxazujRBUvzZmbI/EzbQ89P/pDOb29cewNu4q36hARCRqI8jLy0NkZCQiIyMBAImJiYiMjERSUhIAYM6cORgxYkTF9jdv3sTmzZtx69YtXLhwAW+99RauX7+Or7/+Woz4RKRGHRo64I/pHdHJxxFFpUp8svcapm6JQHZhqdjRiEjPVOvMTcuWLXH06FHUqlULLVq0qPJfUpcvX672Dw8PD0fXrl0rvn4yh87IkSOxYcMGpKSkVBQd4PEaVt9//z1iY2Mhl8vRtWtXnD17Fp6entX+mUSkuxxtzLBhVGusPp2Abw/H4uC1FFy5m4Ulw1qgpYd+zopORNpXrXIzYMAAmJmZAUCVA35V1aVLl+euMg4AGzZsqPR148aNq5xnh4j0n1QqwcTODRDkVRvTtl5GcmYhhqw4h5ndfTCpUwNIpbxMRURVU3ltKX3HtaWI9EdOUSk+3Xsdv/19m3h779pYPKQ5nGzNX/BKIjI0Gl0484nw8HBER0cDAPz8/BAYGFiT3Wgdyw2RfhEEATvD7+KL/TdQWFoOeytTfP9mM64wTmRkNFpu7t69i2HDhiE0NBR2dnYAgKysLLRr1w7btm2ruItJV7HcEOmnuPQ8TNsageiUx7OMc+kGIuOi0VXBx40bh9LSUkRHRyMzMxOZmZmIjo6GUqnEuHHjahyaiKgq3k7W2Pteu4oVxteH3sbAZWcRl54nbjAi0jkqn7mxsLDA2bNn0aJFi0rPX7p0CR07dkRBQYFaA6obz9wQ6b9jMWn4cOdVZOaXwEIuwxf9/DC0NefEITJkGj1z4+7uXrFC+D+Vl5fD1dVV1d0REansFV9nHJreER28HVBYWo7Ze65hypbLyC7gnDhEVINy8+2332LatGkIDw+veC48PBzTp0/Hd999p9ZwRETP42Rrjo1j2mBOL1+YSCX4/Voqev14CucTMsSORkQiq9ZlqVq1alU63Zufn4+ysjKYmDyeJufJn62srJCZmam5tGrAy1JEhufq3Sy8vzUCtzMKIJUAU7p64/1XG0Iu47IsRIZClc/vak3i98MPP6gjFxGRRjR1s8PB9zviy/03sPPSXSw9FofTtx5iyVst4FHbUux4RKRlnMSPiAzKgav3MWfPNeQWlcHazATz+vtjUMu6HGxMpOc0OqD4n4qKipCTk1PpQUQkpr5NXfHH9I5o42mPvOIyzNp5BdO2RnCwMZERUbnc5OfnY+rUqXBycoKVlRVq1apV6UFEJDa3WpbYOiEYH3b3gUwqwYGrKRxsTGREVC43H3/8MY4dO4bly5fDzMwMa9aswbx58+Dq6oqNGzdqIiMRkcpkUgmmvtIQuya1Rb3alrifXYS3Vofhm0MxKClTih2PiDRI5TE3Hh4e2LhxI7p06QJbW1tcvnwZ3t7e2LRpE7Zu3Yrff/9dU1nVgmNuiIxPXnEZ5v092BgAAuoq8MNbzdHA0VrkZERUXRodc5OZmQkvLy8AgK2tbcWt3x06dMCpU6dqEJeISLOszUzw7ZvN8PM7LaGwkOPavWz0XXIGW84nwcjuqSAyCiqXGy8vLyQmJgIAfH19sWPHDgDAb7/9VrGQJhGRLuodUAeHZnREuwa1UVhajk/2XsP4jZeQkVcsdjQiUiOVy83o0aNx5coVAMDs2bOxbNkymJub44MPPsBHH32k9oBEROpUR2GBzWOD8GnvxpDLJPgrOg09fjiN47HpYkcjIjV56Xlubt++XTHupmnTpurKpTEcc0NET9y4n40PtkfiZtrjlcWHB9fDJ70bw8JUJnIyIvpfqnx+cxI/IjJqRaXlWHgoButDbwMAGjha4YehLRDgphA3GBFVovFJ/I4ePYq+ffuiQYMGaNCgAfr27Yu//vqrRmGJiMRkLpfhi37+2DimDZxszBD/IB+v/xyKn47dQrnSqP7tR2QwVC43P//8M3r27AkbGxtMnz4d06dPh62tLXr37o1ly5ZpIiMRkcZ18nHE4Rmd0DvABWVKAd/9eRNDVp5DUkaB2NGISEUqX5Zyc3PD7NmzMXXq1ErPL1u2DF9//TXu3bun1oDqxstSRFQVQRCw5/I9fLH/BvKKy2Bl+vjMzput3Lg+FZGINHpZKisrCz179nzq+e7duyM7O1vV3RER6RSJRILBgW4V61Pll5Tj491XMWHTJTzkLeNEekHlctO/f3/s3bv3qef37duHvn37qiUUEZHY3O0fr081u5cv5DIJjkSlocfiUzgSlSZ2NCJ6AZPqbLRkyZKKP/v5+eE///kPTpw4gbZt2wIAwsLCEBoailmzZmkmJRGRCGRSCSZ1boBODR0xc0ckYlJzMX5jOIa0csPcfv6wNqvWX6FEpGXVGnNTv3796u1MIkFCQsJLh9IkjrkhopooLivHoj9vYtXpBAgC4FbLAt+/2QxBXrXFjkZkFDjPTRVYbojoZYQlZGDWjiu4l1UIiQQY39ELM1/zgbmcE/8RaZLG57l5QhAELjpHREYl2Ks2Ds3oiCGt3CAIwKpTCRjwUyhu3OcNFUS6okblZuPGjQgICICFhQUsLCzQtGlTbNq0Sd3ZiIh0ko25HN+80QyrhgeitpUpYtNyMXDZ44n/ysqVYscjMnoql5tFixZh8uTJ6N27N3bs2IEdO3agZ8+emDRpEhYvXqyJjEREOqm7vwsOf9AJ3f2cUVr+eOK/N1acQ8KDPLGjERk1lcfc1K9fH/PmzcOIESMqPf/LL7/gyy+/RGJioloDqhvH3BCRuj2Z+O/L/TeQW1wGc7kUs3v6YkRbT0ilnPiPSB00OuYmJSUF7dq1e+r5du3aISUlRdXdERHpvScT/x3+oBM6eDugqFSJL3+Lwrtrz+PuIy7fQKRtKpcbb29v7Nix46nnt2/fjoYNG6olFBGRPnK1s8DGMW3w1QB/WMhlOBufgZ4/nMaOi8m8+YJIi1S+LLV7924MHToU3bp1Q/v27QEAoaGhOHr0KHbs2IHXX39dI0HVhZeliEgbEh/m48OdV3DpziMAwKu+Tpg/KABOtuYiJyPSTxqf5+by5ctYtGgRoqOjAQCNGzfGrFmz0KJFi5ol1iKWGyLSlnKlgDWnE/D9nzdRUq6EnaUc8/r7o38zVy7CSaQijZWb0tJSTJw4EZ9//nm1Zy3WNSw3RKRtsam5mLUzEtfv5QAAejVxwb8HNkFtazORkxHpD40NKJbL5di9e/dLhSMiMjaNXGyw9732mNGtIUykEvxxPRXdF5/Coeu8CYNIE1QeUDxw4ECEhIRoIAoRkeGSy6SY0c0HIVPao5GzDTLySzBp82VM3xaBR/klYscjMigqL2nbsGFDfPXVVwgNDUVgYCCsrKwqff/9999XWzgiIkPTpK4C+6e1x5Kjt7D8RDz2Rd7H2fgMfP16AF7zcxY7HpFBqNEkfs/dGVcFJyKqtsjkLHy48wri0h/PaDyoRV180c8fCku5yMmIdA9XBa8Cyw0R6ZKi0nIs/usmVp9KgFIAnG3NMH9QAF7x5Vkcon/iquBERHrCXC7DnF6NsWtyO3g5WiEtpxhjNoRj1o4ryC4oFTsekV6qUblZu3YtmjRpAnNzc5ibm6NJkyZYs2aNurMRERmNlh618Pv7HTG+Y31IJMDuy3fx2uKTOBqdJnY0Ir2jcrmZO3cupk+fjn79+mHnzp3YuXMn+vXrhw8++ABz587VREYiIqNgLpfh0z5+2DWpLbwcrJCeW4yxv4Rj5vZIZBXwjiqi6lJ5zI2joyOWLFmCYcOGVXp+69atmDZtGh4+fKjWgOrGMTdEpA+KSsvx/Z+xWHMmEYIAONqY8Y4qMmoaHXNTWlqKVq1aPfV8YGAgysrKVN0dERE9Q6WzOI5WeJBbjPEbwzGD8+IQvZDK5Wb48OFYvnz5U8+vWrUK77zzjlpCERHRY4H17PH7+x0xsbMXpBIgJPI+XuPsxkRVUvmy1LRp07Bx40a4u7sjODgYAHD+/HkkJSVhxIgRkMv/Oz/DokWL1JtWDXhZioj0VUTSI3y86ypu/T0vTp+mdTCvvz8cuEYVGQGNznPTtWvXam0nkUhw7NgxVXatFSw3RKTPisvKseToLaw4mYBypQB7K1N82d8f/ZrW4UrjZNA4iV8VWG6IyBBcv5eND3deQUxqLgDgNT9n/HtgEzjbmoucjEgztDaJHxERiaNJXQX2T+2AD7r5QC6T4EhUGrotOokdF5M5uSoZPZYbIiI9ZWoixfRuDfHbtA5o6qZAblEZPt59FSPWXcDdRwVixyMSDcsNEZGe83WxxZ7J7TCnly/MTKQ4feshui8+hV/O3oZSybM4ZHxYboiIDICJTIqJnRvgj+kd0dqzFgpKyvHF/hsYsvIc4h/kiR2PSKtYboiIDIiXozW2T2iLrwb4w8pUhvA7j9Drx9NYdjwOpeVKseMRaQXLDRGRgZFKJRjR1hOHP+iETj6OKClT4tvDsRjwUyiu38sWOx6RxrHcEBEZKLdalvhldGt8/2YzKCzkiErJwYBloVh4KAZFpeVixyPSGJYbIiIDJpFIMDjQDX/N7Iw+AXVQrhSw/EQ8ev14GucTMsSOR6QRLDdEREbA0cYMy95piVXDA+Fsa4bEh/kYuioMn+y9hpyiUrHjEakVyw0RkRHp7u+CIzM7Y1gbDwDAlvNJeG3RSRy+kSpyMiL1YbkhIjIytuZyzB8UgG0TglHfwQppOcWYuOkS3vv1EtJzisSOR/TSWG6IiIxUsFdt/DG9I97r0gAyqQS/X0vFq4tOYtuFJC7hQHqN5YaIyIiZy2X4uKcv9k9tj4C6j5dwmL3nGoatDkPiw3yx4xHVCMsNERHB31WBve+1w6e9G8NcLkVYQiZ6/HCKk/+RXmK5ISIiAI+XcBjfyQtHPuiMjg0dKib/67f0DCKSHokdj6jaWG6IiKgSd3tLbBzTBouHNkMtSzliUnMxaPlZfLn/BvKKy8SOR/RCLDdERPQUiUSC11u44eisLhjUsi4EAdhw9jZeW3QSR6LSxI5HVCVRy82pU6fQr18/uLq6QiKRICQk5IWv+fXXX9GsWTNYWlqiTp06GDNmDDIyOMsmEZEm2FuZYtGQ5tg8Ngge9pZIyS7C+I3hmLTpElKzeds46SZRy01+fj6aNWuGZcuWVWv70NBQjBgxAmPHjsWNGzewc+dOXLhwAePHj9dwUiIi49ahoQMOz+iESZ0f3zZ+6EYqui06iU3nbkOp5G3jpFskgo5MZiCRSLB3714MHDjwudt89913WL58OeLj4yueW7p0KRYuXIi7d+9W6+fk5ORAoVAgOzsbtra2LxubiMjoRN3PwZy913AlOQsA0MLDDvMHBcDXhX+nkuao8vmtV2Nu2rZti+TkZPz+++8QBAFpaWnYtWsXevfu/dzXFBcXIycnp9KDiIhqzs/VFnsmt8O8/v6wNjNBRFIW+i45gwV/xKCwhKuNk/j0qty0b98ev/76K4YOHQpTU1O4uLhAoVBUeVlr/vz5UCgUFQ93d3ctJiYiMkwyqQQj23niyMxO6OHvjDKlgBUn49H9h5M4EZsudjwycnpVbqKiojB9+nTMnTsXly5dwqFDh3D79m1MmjTpua+ZM2cOsrOzKx7JyclaTExEZNjqKCywcngrrBoeiDoKcyRnFmLU+ouYtjUC6bkccEzi0KsxN8OHD0dRURF27txZ8dyZM2fQsWNH3L9/H3Xq1Hnhz+GYGyIizcgrLsPiIzexPjQRSgGwMTfB//X0xdttPCCVSsSOR3rOYMfcFBQUQCqtHFkmkwEAF3kjIhKZtZkJPu/rh31TOlSsU/VZyHUMXnEW0Skc70jaI2q5ycvLQ2RkJCIjIwEAiYmJiIyMRFJSEoDHl5RGjBhRsX2/fv2wZ88eLF++HAkJCQgNDcX777+PNm3awNXVVYy3QERE/yPATYGQKe3xZT+//w44XnoGX/8ejYISznBMmifqZakTJ06ga9euTz0/cuRIbNiwAaNGjcLt27dx4sSJiu8tXboUK1asQGJiIuzs7PDKK69g4cKFqFu3brV+Ji9LERFpT2p2Eeb9dgN/XE8FANS1s8C8/v7o5ucscjLSN6p8fuvMmBttYbkhItK+o9FpmLvvBu5lFQIAuvs548v+/nC1sxA5GekLgx1zQ0RE+unVxs44MvPxDMcmUgn+jEpDt0UnsfpUAkrLlWLHIwPDckNERFphaWqC2b18cfD9jmhVrxYKSsrxn9+j0W/pGVy6kyl2PDIgLDdERKRVjVxssGNiW3wzuClqWcoRk5qLwcvPYfbuq3iUXyJ2PDIALDdERKR1UqkEQ1q74+isLhjSyg0AsO1iMl75/gR2XEzmYpz0UjigmIiIRHfxdiY+23sdsWm5AIBW9WrhXwOboHEd/j1Nj3FAMRER6ZXWnvY48H4HfNLbF5amMoTfeYS+S8/g3weikFfMuXFINSw3RESkE+QyKSZ0aoC/ZnZGT38XlCsFrDmTiG7fn8TBqymciZ6qjeWGiIh0iqudBVYMD8T6Ua3hYW+J1JwiTNlyGSPWXUDCgzyx45EeYLkhIiKd1NXXCX9+0AnTX20IUxMpTt96iJ4/nMb3f8aiqLRc7Hikw1huiIhIZ5nLZfjgNR/8OaMTOvs4oqRciaXH4tBt0Un8FZUmdjzSUSw3RESk8zwdrLBhdGuseLclXBXmuPuoEOM2hmPcLxeRnFkgdjzSMSw3RESkFyQSCXo2qYO/ZnWuWMbhr+h0dFt0EkuO3uKlKqrAckNERHrlyTIOh2Z0RFuv2iguU2LRkZvo8cMpHI9JFzse6QCWGyIi0kveTjbYMj4IS4a1gLOtGe5kFGD0hosYvzGcl6qMHMsNERHpLYlEgv7NXHF0VhdM6OQFE6kER/5ecfyHv27yUpWR4vILRERkMG6l5WLuvhs4l5ABAPCwt8QX/fzwamNnkZPRy1Ll85vlhoiIDIogCDhwNQX/ORiN1JwiAMCrvk6Y288P9WpbiZyOaoprSxERkdGSSCTo18wVR2d1xsTOXpDLJDgak47XFp/C93/GorCEl6oMHc/cEBGRQYtLz8O8327g9K2HAIC6dhb4rE9j9GziAolEInI6qi5elqoCyw0RkfERBAGHb6ThXweicC+rEADQwdsBX/b3g7eTjcjpqDpYbqrAckNEZLwKS8qx/GQ8VpyMR0mZEiZSCUa188T0bg1hYy4XOx5VgeWmCiw3RESUlFGAfx2MwpG/16dysDbD//VshMEt3SCV8lKVLmK5qQLLDRERPXEiNh1f/RaFhIf5AIDm7naY198fzdztxA1GT2G5qQLLDRER/VNJmRLrQxOx5Ogt5JeUQyIBhgS646OejeBgbSZ2PPoby00VWG6IiOhZ0nKKsPCPGOyJuAcAsDE3wYxuPhjRth7kMs6cIjaWmyqw3BARUVUu3cnEF/tv4Pq9HABAQydrzO3nh44NHUVOZtxYbqrAckNERC9SrhSwIzwZ3x6ORWZ+CQDgNT9nfNanMWc5FgnLTRVYboiIqLqyC0vx41+38Mu52yhXCjCVSTGuY31M6eoNKzMTseMZFZabKrDcEBGRqm6l5eKrA1EVsxw72Zhhdi9fDGxel7eOawnLTRVYboiIqCYEQcCRqDT8+2A0kjILAAAtPOzwZT/eOq4NLDdVYLkhIqKXUVxWjrVnEvHTsTgU/L0I5xuBbvi4RyM42ZqLnM5wsdxUgeWGiIjUIS2nCN8cisXuy3cBAFamMkx5xRtj2teHuVwmcjrDw3JTBZYbIiJSp4ikR5j3WxQik7MAAO72Fvi0tx96+Dtz1XE1YrmpAssNERGpm1IpYN+Ve1jwRwzScooBAG29auPzvn7wc+VnjTqw3FSB5YaIiDQlv7gMK07GY9WpBBSXKSGVAENbe2BWdx8u5fCSWG6qwHJDRESadvdRAeb/EYODV1MAADZmJpj6ijdGtfeEmQnH49QEy00VWG6IiEhbLiRm4qsD/13KoV5tS3zSuzG6+3E8jqpYbqrAckNERNqkVArYdekuvv0zFg9yOR6nplhuqsByQ0REYsgrLsOKE/FYdToBJWVKSCTA0FbumNndB042nB/nRVhuqsByQ0REYrr7qAAL/ojBgb/H41iZyvBeV2+M7cD5carCclMFlhsiItIFl+5k4qsD0bjy9/w4de0sMLuXL/o2rcPxOM/AclMFlhsiItIVT+bHWfhHLFJzigAALT3s8HlfP7TwqCVyOt3CclMFlhsiItI1hSXlWHUqAStOxqOw9PF6Vf2bueLjno3gVstS5HS6geWmCiw3RESkq9JyivDd4VjsunwXggCYmkgxrkN9TO7SADbmcrHjiYrlpgosN0REpOuu38vGvw9GISwhEwDgYG2KD17zwdBW7jCRSUVOJw6Wmyqw3BARkT4QBAF/Rafj69+jkfgwHwDQ0Mkan/RpjC4+jkY36JjlpgosN0REpE9KypT49fwd/Hj0FrIKSgEAHRs64JPejdG4jvF8jrHcVIHlhoiI9FF2QSmWHruFX87dRmm5AIkEGBLojlndfeBka/iTALLcVIHlhoiI9FlSRgEWHorBwWuPJwG0NJVhQicvTOjkBUtTE5HTaQ7LTRVYboiIyBBcupOJfx+MRkRSFgDAycYMs7r74I1Ad8ikhjceh+WmCiw3RERkKARBwO/XUrHwUAySMgsAAI2cbTCnty86G9igY5abKrDcEBGRoSkuK8emc3ew9FgcsgsfDzru4O2AOb194e+qEDmderDcVIHlhoiIDFVWQQl+OhaHjefuoKT88crjg1q44cMePqijsBA73kthuakCyw0RERm65MwCfHM4Fr9duQ8AMDORYmyH+pjUpQFs9XSmY5abKrDcEBGRsYhMzsLXB6Nx4fbjmY7trUzx/iveeDuoHkxN9GumY5abKrDcEBGRMREEAUei0rDwUAziHzye6diztiU+7umLXk1c9GbQMctNFVhuiIjIGJWVK7E9PBmLj9zCw7xiAEBzdzt80rsx2tS3Fzndi7HcVIHlhoiIjFl+cRlWnUrA6tMJKCgpBwB0a+yM2b0awdvJRuR0z8dyUwWWGyIiIiA9twg//nUL2y4mo1wpQCoBhrZ2x4xuPnDWweUcWG6qwHJDRET0X/EP8vDNoRgcvpEGADCXSzGugxcmdvaCjQ7dWcVyUwWWGyIioqeF387E/D9icOnOIwCP76ya2tUb7wR7wMxEJnI6lpsqsdwQERE9myAI+DMqDd/8484qd3sLfNi9Efo1dYVUxDWrWG6qwHJDRERUtbJyJXZeuovFR24iPffxnVV+dWwxu5cvOvk4ipKJ5aYKLDdERETVU1hSjnWhiVhxIh65xWUAgPbetfF/PX3R1M1Oq1lYbqrAckNERKSazPwSLDseh01/r1kFAH2a1sGH3RuhvoOVVjKw3FSB5YaIiKhmkjMLsPivm9gbcQ+CAJhIJRja2h3TX20IJw3fPs5yUwWWGyIiopcTnZKDbw/H4lhMOgDAQi7DmA6emNCpARQWmrl9nOWmCiw3RERE6hGWkIGFh2IQkZQFALCzlOO9Lg0woq0nzOXqvX1clc9vUZcEPXXqFPr16wdXV1dIJBKEhIRUuf2oUaMgkUieevj7+2snMBEREVUI9qqNPZPbYeXwQHg7WSOroBRf/x6Drt+dqFi/Sgyilpv8/Hw0a9YMy5Ytq9b2P/74I1JSUioeycnJsLe3x5tvvqnhpERERPQsEokEPfxdcHhGJ3zzRlO4Kszh7WQNB2sz0TKZiPaTAfTq1Qu9evWq9vYKhQIKhaLi65CQEDx69AijR49+7muKi4tRXPzf9piTk1OzsERERPRcMqkEQ1q5o38zV2QVlIqaRdQzNy9r7dq16NatG+rVq/fcbebPn19RihQKBdzd3bWYkIiIyLiYy2VwUYi78Kbelpv79+/jjz/+wLhx46rcbs6cOcjOzq54JCcnaykhERERiUHUy1Iv45dffoGdnR0GDhxY5XZmZmYwMxPvuh8RERFpl16euREEAevWrcPw4cNhamoqdhwiIiLSIXpZbk6ePIm4uDiMHTtW7ChERESkY0S9LJWXl4e4uLiKrxMTExEZGQl7e3t4eHhgzpw5uHfvHjZu3FjpdWvXrkVQUBCaNGmi7chERESk40QtN+Hh4ejatWvF1zNnzgQAjBw5Ehs2bEBKSgqSkpIqvSY7Oxu7d+/Gjz/+qNWsREREpB+4/AIRERHpPL1ZfoGIiIhI3VhuiIiIyKCw3BAREZFBYbkhIiIig8JyQ0RERAaF5YaIiIgMit6uLVVTT+58z8nJETkJERERVdeTz+3qzGBjdOUmNzcXAODu7i5yEiIiIlJVbm4uFApFldsY3SR+SqUS9+/fh42NDSQSiVr3nZOTA3d3dyQnJ3OCQA3icdYOHmft4HHWHh5r7dDUcRYEAbm5uXB1dYVUWvWoGqM7cyOVSuHm5qbRn2Fra8v/cLSAx1k7eJy1g8dZe3istUMTx/lFZ2ye4IBiIiIiMigsN0RERGRQWG7UyMzMDF988QXMzMzEjmLQeJy1g8dZO3ictYfHWjt04Tgb3YBiIiIiMmw8c0NEREQGheWGiIiIDArLDRERERkUlhsiIiIyKCw3L7Bs2TJ4enrC3NwcQUFBuHDhQpXb79y5E76+vjA3N0dAQAB+//33St8fNWoUJBJJpUfPnj01+Rb0grqPMwBER0ejf//+UCgUsLKyQuvWrZGUlKSpt6AX1H2c//d3+cnj22+/1eTb0HnqPs55eXmYOnUq3NzcYGFhAT8/P6xYsUKTb0EvqPs4p6WlYdSoUXB1dYWlpSV69uyJW7duafIt6AVVjvONGzcwePBgeHp6QiKR4IcffnjpfdaIQM+1bds2wdTUVFi3bp1w48YNYfz48YKdnZ2Qlpb2zO1DQ0MFmUwmfPPNN0JUVJTw2WefCXK5XLh27VrFNiNHjhR69uwppKSkVDwyMzO19ZZ0kiaOc1xcnGBvby989NFHwuXLl4W4uDhh3759z92nMdDEcf7n73FKSoqwbt06QSKRCPHx8dp6WzpHE8d5/PjxQoMGDYTjx48LiYmJwsqVKwWZTCbs27dPW29L56j7OCuVSiE4OFjo2LGjcOHCBSEmJkaYMGGC4OHhIeTl5WnzrekUVY/zhQsXhA8//FDYunWr4OLiIixevPil91kTLDdVaNOmjTBlypSKr8vLywVXV1dh/vz5z9x+yJAhQp8+fSo9FxQUJEycOLHi65EjRwoDBgzQSF59pYnjPHToUOHdd9/VTGA9pYnj/L8GDBggvPLKK+oJrKc0cZz9/f2Fr776qtI2LVu2FD799FM1Jtcv6j7OsbGxAgDh+vXrlfbp6OgorF69WgPvQD+oepz/qV69es8sNy+zz+riZannKCkpwaVLl9CtW7eK56RSKbp164Zz58498zXnzp2rtD0A9OjR46ntT5w4AScnJzRq1AiTJ09GRkaG+t+AntDEcVYqlTh48CB8fHzQo0cPODk5ISgoCCEhIRp7H7pOk7/PT6SlpeHgwYMYO3as+oLrGU0d53bt2mH//v24d+8eBEHA8ePHcfPmTXTv3l0zb0THaeI4FxcXAwDMzc0r7dPMzAxnzpxR91vQCzU5zmLs81lYbp7j4cOHKC8vh7Ozc6XnnZ2dkZqa+szXpKamvnD7nj17YuPGjTh69CgWLlyIkydPolevXigvL1f/m9ADmjjO6enpyMvLw4IFC9CzZ0/8+eefeP311zFo0CCcPHlSM29Ex2nq9/mffvnlF9jY2GDQoEHqCa2HNHWcly5dCj8/P7i5ucHU1BQ9e/bEsmXL0KlTJ/W/CT2giePs6+sLDw8PzJkzB48ePUJJSQkWLlyIu3fvIiUlRTNvRMfV5DiLsc9nMbpVwcX21ltvVfw5ICAATZs2RYMGDXDixAm8+uqrIiYzHEqlEgAwYMAAfPDBBwCA5s2b4+zZs1ixYgU6d+4sZjyDtW7dOrzzzjuV/uVL6rF06VKEhYVh//79qFevHk6dOoUpU6bA1dX1qbMRVDNyuRx79uzB2LFjYW9vD5lMhm7duqFXr14QOJG/3mG5eQ4HBwfIZDKkpaVVej4tLQ0uLi7PfI2Li4tK2wOAl5cXHBwcEBcXZ5TlRhPH2cHBASYmJvDz86u0TePGjY329LKmf59Pnz6N2NhYbN++XX2h9ZAmjnNhYSE++eQT7N27F3369AEANG3aFJGRkfjuu++Mstxo6vc5MDAQkZGRyM7ORklJCRwdHREUFIRWrVqp/03ogZocZzH2+Sy8LPUcpqamCAwMxNGjRyueUyqVOHr0KNq2bfvM17Rt27bS9gBw5MiR524PAHfv3kVGRgbq1KmjnuB6RhPH2dTUFK1bt0ZsbGylbW7evIl69eqp+R3oB03/Pq9duxaBgYFo1qyZeoPrGU0c59LSUpSWlkIqrfzXtUwmqzhLaWw0/fusUCjg6OiIW7duITw8HAMGDFDvG9ATNTnOYuzzmdQ2NNkAbdu2TTAzMxM2bNggREVFCRMmTBDs7OyE1NRUQRAEYfjw4cLs2bMrtg8NDRVMTEyE7777ToiOjha++OKLSrca5ubmCh9++KFw7tw5ITExUfjrr7+Eli1bCg0bNhSKiopEeY+6QN3HWRAEYc+ePYJcLhdWrVol3Lp1S1i6dKkgk8mE06dPa/396QpNHGdBEITs7GzB0tJSWL58uVbfj67SxHHu3Lmz4O/vLxw/flxISEgQ1q9fL5ibmws///yz1t+frtDEcd6xY4dw/PhxIT4+XggJCRHq1asnDBo0SOvvTZeoepyLi4uFiIgIISIiQqhTp47w4YcfChEREcKtW7eqvU91YLl5gaVLlwoeHh6Cqamp0KZNGyEsLKzie507dxZGjhxZafsdO3YIPj4+gqmpqeDv7y8cPHiw4nsFBQVC9+7dBUdHR0Eulwv16tUTxo8fr9b/Q/WVOo/zE2vXrhW8vb0Fc3NzoVmzZkJISIim34bO08RxXrlypWBhYSFkZWVpOr7eUPdxTklJEUaNGiW4uroK5ubmQqNGjYTvv/9eUCqV2ng7Okvdx/nHH38U3NzcBLlcLnh4eAifffaZUFxcrI23otNUOc6JiYkCgKcenTt3rvY+1UEiCBwpRURERIaDY26IiIjIoLDcEBERkUFhuSEiIiKDwnJDREREBoXlhoiIiAwKyw0REREZFJYbIiIiMigsN0RERGRQWG6IiIjIoLDcEBERkUFhuSEiIiKDwnJDRAbjwoUL6NixI2xsbGBlZYWAgABcvHhR7FhEpGUmYgcgIlKXt956C+3atcOqVatgbm6OhIQEODs7ix2LiLSM5YaIDEZZWRk8PDzg7e0NuVyO+vXrix2JiEQgEQRBEDsEEZE6hIeH4/XXX0dKSgrMzc1x7949KBQKsWMRkZax3BCRwXjttddgY2OD2bNnw97eHl5eXpBKObSQyNiw3BCRQXj48CEcHR0RGRmJZs2aiR2HiETEckNEBsPDwwMtWrTA3Llz4eDggMTERJSUlKB79+5iRyMiLeL5WiIyGH/88QeUSiV69OgBHx8fjB8/HmlpaWLHIiIt45kbIiIiMig8c0NEREQGheWGiIiIDArLDRERERkUlhsiIiIyKCw3REREZFBYboiIiMigsNwQERGRQWG5ISIiIoPCckNEREQGheWGiIiIDArLDRERERmU/wftH4DibxtGGAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# from math import factorial\n", "import numpy as np\n", "from math import comb as C\n", "\n", "def solve2(N, eps):\n", " pN, pS, pE, pW = 0.25, 0.25, 0.25+eps, 0.25-eps \n", " return sum([\n", " sum([C(2*n, n) * C(n,k)**2 * (pN*pS)**k * (pE*pW)**(n-k) for k in range(n+1)\n", " ]) for n in range(N)\n", " ])\n", "\n", "F = lambda t: solve2(200, t)\n", "X = np.linspace(0.05, 0.1, 100)\n", "Y = np.vectorize(F)(X)\n", "\n", "plt.plot(X, Y)\n", "plt.xlabel(\"ε\"); plt.ylabel(\"probability of return\")\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "0.061912524105980984" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.optimize import bisect\n", "\n", "bisect(lambda eps: solve2(250, eps) - 2, 0.06, 0.062)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With this approach, we are able to get 3 accurate digits. However, this approach has two main issues\n", "\n", " - the sum as is is very costly to compute" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 1.65 s, sys: 3 µs, total: 1.65 s\n", "Wall time: 1.66 s\n" ] }, { "data": { "text/plain": [ "2.0000000000030997" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%time solve2(250, 0.061912524105980984)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " - attempting to sum the series to more terms results in an overflow error" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "scrolled": true }, "outputs": [ { "ename": "OverflowError", "evalue": "int too large to convert to float", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mOverflowError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[5], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43msolve2\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m500\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m0.061912524105980984\u001b[39;49m\u001b[43m)\u001b[49m\n", "Cell \u001b[0;32mIn[2], line 7\u001b[0m, in \u001b[0;36msolve2\u001b[0;34m(N, eps)\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21msolve2\u001b[39m(N, eps):\n\u001b[1;32m 6\u001b[0m pN, pS, pE, pW \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0.25\u001b[39m, \u001b[38;5;241m0.25\u001b[39m, \u001b[38;5;241m0.25\u001b[39m\u001b[38;5;241m+\u001b[39meps, \u001b[38;5;241m0.25\u001b[39m\u001b[38;5;241m-\u001b[39meps \n\u001b[0;32m----> 7\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28msum\u001b[39m([\n\u001b[1;32m 8\u001b[0m \u001b[38;5;28msum\u001b[39m([C(\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39mn, n) \u001b[38;5;241m*\u001b[39m C(n,k)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m (pN\u001b[38;5;241m*\u001b[39mpS)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mk \u001b[38;5;241m*\u001b[39m (pE\u001b[38;5;241m*\u001b[39mpW)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m(n\u001b[38;5;241m-\u001b[39mk) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(n\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m 9\u001b[0m ]) \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(N)\n\u001b[1;32m 10\u001b[0m ])\n", "Cell \u001b[0;32mIn[2], line 8\u001b[0m, in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21msolve2\u001b[39m(N, eps):\n\u001b[1;32m 6\u001b[0m pN, pS, pE, pW \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0.25\u001b[39m, \u001b[38;5;241m0.25\u001b[39m, \u001b[38;5;241m0.25\u001b[39m\u001b[38;5;241m+\u001b[39meps, \u001b[38;5;241m0.25\u001b[39m\u001b[38;5;241m-\u001b[39meps \n\u001b[1;32m 7\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28msum\u001b[39m([\n\u001b[0;32m----> 8\u001b[0m \u001b[38;5;28msum\u001b[39m([C(\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39mn, n) \u001b[38;5;241m*\u001b[39m C(n,k)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m2\u001b[39m \u001b[38;5;241m*\u001b[39m (pN\u001b[38;5;241m*\u001b[39mpS)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mk \u001b[38;5;241m*\u001b[39m (pE\u001b[38;5;241m*\u001b[39mpW)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m(n\u001b[38;5;241m-\u001b[39mk) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(n\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m 9\u001b[0m ]) \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(N)\n\u001b[1;32m 10\u001b[0m ])\n", "Cell \u001b[0;32mIn[2], line 8\u001b[0m, in \u001b[0;36m\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21msolve2\u001b[39m(N, eps):\n\u001b[1;32m 6\u001b[0m pN, pS, pE, pW \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m0.25\u001b[39m, \u001b[38;5;241m0.25\u001b[39m, \u001b[38;5;241m0.25\u001b[39m\u001b[38;5;241m+\u001b[39meps, \u001b[38;5;241m0.25\u001b[39m\u001b[38;5;241m-\u001b[39meps \n\u001b[1;32m 7\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28msum\u001b[39m([\n\u001b[0;32m----> 8\u001b[0m \u001b[38;5;28msum\u001b[39m([\u001b[43mC\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mn\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mn\u001b[49m\u001b[43m)\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mC\u001b[49m\u001b[43m(\u001b[49m\u001b[43mn\u001b[49m\u001b[43m,\u001b[49m\u001b[43mk\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m2\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43m \u001b[49m\u001b[43m(\u001b[49m\u001b[43mpN\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mpS\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mk\u001b[49m \u001b[38;5;241m*\u001b[39m (pE\u001b[38;5;241m*\u001b[39mpW)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m(n\u001b[38;5;241m-\u001b[39mk) \u001b[38;5;28;01mfor\u001b[39;00m k \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(n\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m 9\u001b[0m ]) \u001b[38;5;28;01mfor\u001b[39;00m n \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(N)\n\u001b[1;32m 10\u001b[0m ])\n", "\u001b[0;31mOverflowError\u001b[0m: int too large to convert to float" ] } ], "source": [ "solve2(500, 0.061912524105980984)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We bypass both of these issues by evaluating the sum symbolically." ] }, { "cell_type": "markdown", "metadata": { "id": "v_ysEil_3mtn" }, "source": [ "### Attempt 3: Symbolic manipulation\n", "\n", "Substituting $x = p_E p_W, y = p_N p_S$, we attempt to evaluate the inner sum symbolically with `sympy`" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 83 }, "executionInfo": { "elapsed": 1115, "status": "ok", "timestamp": 1603767346645, "user": { "displayName": "Husnain Raza", "photoUrl": "", "userId": "16146658556003895350" }, "user_tz": 300 }, "id": "YpKn5vaq35FE", "outputId": "b8d2aa01-486a-4e91-e061-84c1c123d04b" }, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle x^{k} \\left(\\begin{cases} {{}_{2}F_{1}\\left(\\begin{matrix} - k, - k \\\\ 1 \\end{matrix}\\middle| {\\frac{y}{x}} \\right)} & \\text{for}\\: \\left|{\\frac{y}{x}}\\right| \\leq 1 \\\\\\sum_{j=0}^{k} x^{- j} y^{j} {\\binom{k}{j}}^{2} & \\text{otherwise} \\end{cases}\\right)$" ], "text/plain": [ "x**k*Piecewise((hyper((-k, -k), (1,), y/x), Abs(y/x) <= 1), (Sum(y**j*binomial(k, j)**2/x**j, (j, 0, k)), True))" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy\n", "\n", "# sympy.init_printing()\n", "x,y = sympy.symbols(\"x y\")\n", "k,j = sympy.symbols(\"k j\", integer=True, positive=True)\n", "sympy.Sum(sympy.binomial(k, j)**2 * x**(k-j) * y**j, (j,0,k)).doit()" ] }, { "cell_type": "markdown", "metadata": { "id": "53M0wZUo37RQ" }, "source": [ "Therefore, we have that \n", "$$E =\\sum_{k=0}^{\\infty}\\sum_{j=0}^k \\binom{2k}{k} \\binom{k}{j}^2 x^{k-j} y^{j} = \\sum_{k=0}^{\\infty} \\binom{2k}{k} x^k {{}_{2}F_{1}\\left(\\begin{matrix} - k, - k \\\\ 1 \\end{matrix}\\middle| {\\frac{y}{x}} \\right)} $$\n", "\n", "where ${}_2F_1$ is the [hypergeometric function](https://en.wikipedia.org/wiki/Hypergeometric_function)" ] }, { "cell_type": "markdown", "metadata": { "id": "eVfRrzB7lMWk" }, "source": [ "To avoid overflow when calculating the binomial coefficients, we recursively calculate them - defining $a_k = \\binom{2k}{k} x^k $, we have $a_0 = 0$ and that $a_k$ satisfies the following recurrence: \n", "\n", "$$\\begin{align*}\n", "\\frac{a_{k+1}}{a_k} &= \\frac{\\binom{2k+2}{k+1}x^{k+1}}{\\binom{2k}{k}x^k} = \\frac{(2k+1)(2k+2)x}{(k+1)^2} \\\\\n", "&\\implies a_{k+1} = a_k \\cdot \\frac{(2k+1)(2k+2)x}{(k+1)^2} \\\\\n", "&\\iff a_{k} = a_{k-1} \\cdot \\frac{(2k-1)\\cdot2k \\cdot x}{k^2}\n", "\\end{align*}$$\n", "\n", "We thus have the sum is \n", "\n", "$$ E = \\sum_{n=0}^\\infty a_n \\cdot {{}_{2}F_{1}\\left(\\begin{matrix} - n, - n \\\\ 1 \\end{matrix}\\middle| {\\frac{y}{x}} \\right)} $$\n", "\n", "We also use `mpmath` for extended precision" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 229 }, "executionInfo": { "elapsed": 1107, "status": "error", "timestamp": 1603767346653, "user": { "displayName": "Husnain Raza", "photoUrl": "", "userId": "16146658556003895350" }, "user_tz": 300 }, "id": "rLxGibLv36cM", "outputId": "aabdb8dd-d1b9-4a2a-a6de-bb0af394042e" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 55.5 s, sys: 275 ms, total: 55.7 s\n", "Wall time: 1min\n" ] }, { "data": { "text/plain": [ "0.06191395447403192" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from scipy.special import hyp2f1\n", "from functools import lru_cache\n", "import mpmath\n", "\n", "mpmath.mp.dps = 20\n", "\n", "def solve3(N, eps):\n", " pN, pS, pE, pW = 0.25, 0.25, 0.25+eps, 0.25-eps \n", " x,y = mpmath.mpf(pE*pW), mpmath.mpf(pN*pS)\n", " # return x,y\n", " # Automatically creates cache to store previous values\n", " @lru_cache(maxsize=None)\n", " def a(n):\n", " return 1 if n == 0 else a(n-1) * (2*n-1) * 2*n*x / (n*n)\n", " return sum([a(n)*mpmath.hyp2f1(-n,-n,1,y/x) for n in range(N)])\n", "\n", "%time bisect(lambda t: solve3(2000, t) - 2, 0.0619, 0.062)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that this answer is accurate to more than 10 digits - however we can significantly speed up execution time as follows:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**TODO**: Sage and `sympy` have difficulty with the inner binomial sum; however Mathematica / Wolfram Alpha can [prove](https://www.wolframalpha.com/input?i=sum%28C%28k%2Cj%29%5E2+*+z%5Ej%29+from+j+%3D+0+to+k) that the inner sum can be represented with the [Legendre polynomials](https://en.wikipedia.org/wiki/Legendre_polynomials) $P_n(x)$\n", "\n", "$$ \\sum_{j=0}^k \\binom{k}{j}^2 z^j = (1-z)^k P_k\\left(\\frac{1+z}{1-z}\\right)$$\n", "\n", "We have, following [this StackExchange answer](https://math.stackexchange.com/a/766225), that:\n", "\n", "$$\n", "\\begin{align}\n", "P_n(x) &= \\frac{1}{2^n n!}\\frac{d^n}{dx^n}\\left((x+1)^n\\cdot(x-1)^n\\right) \\\\\n", "&= \\frac{1}{2^n n!} \\sum_{k=0}^n \\binom{n}{k} \\cdot \\left( \\frac{d^k}{dx^k}(x+1)^n \\right) \\cdot \\left( \\frac{d^{n-k}}{dx^{n-k}} (x-1)^n \\right) \\\\\n", "&= \\frac{1}{2^n n!} \\sum_{k=0}^n \\binom{n}{k} \\cdot \\frac{n!}{(n-k)!} (x+1)^{n-k} \\cdot \\frac{n!}{k!} (x-1)^k \\\\\n", "&= \\frac{1}{2^n} \\sum_{k=0}^n \\binom{n}{k}^2 (x-1)^k (x+1)^{n-k} \\\\\n", "&= \\left(\\frac{x+1}{2}\\right)^n \\sum_{k=0}^n \\binom{n}{k}^2 \\left(\\frac{x-1}{x+1}\\right)^k \\\\\n", "&\\implies P\\left(\\frac{z+1}{z-1}\\right) = (1-z)^{-n} \\sum_{k=0}^n \\binom{n}{k}^2 z^k & \\; \\left(z \\mapsto \\frac{x-1}{x+1}\\right) \\\\\n", "&\\implies (1-z)^n P\\left(\\frac{z+1}{z-1}\\right) = \\sum_{k=0}^n \\binom{n}{k}^2 z^k\n", "\\end{align}\n", "$$\n", "\n", "as desired. Therefore, \n", "\n", "$$ \\begin{align*}\n", "E &= \\sum_{k=0}^\\infty \\binom{2k}{k} x^{k} \\left[ \\sum_{j=0}^k \\binom{k}{j}^2 \\left(\\frac{y}{x}\\right)^{j} \\right] \\\\\n", "&= \\sum_{k=0}^\\infty \\binom{2k}{k} x^k (1-y/x)^k P_k\\left(\\frac{1+y/x}{1-y/x}\\right)\\\\\n", "&= \\sum_{k=0}^\\infty \\binom{2k}{k} (x-y)^k P_k\\left(\\frac{x+y}{x-y}\\right) \\\\\n", "\\end{align*}\n", "$$\n", "\n", "Expressing this as a sum of Legendre polynomials, we have the following lemma:\n", "\n", "**Lemma (Bornemann p136)**: Let $f(z) = \\sum_{k \\ge 0} a_k z^k $. Using [Laplace's integral representation of the Legendre polynomials](https://dlmf.nist.gov/18.10#E5), we have\n", "\n", "$$\\begin{align}\n", "\\sum_{k\\ge0} a_k P_k(x) z^k &= \\sum_{k\\ge0} a_k \\left( \\frac{1}{\\pi} \\int_0^\\pi (x + \\sqrt{x^2-1} \\cos{\\theta})^k d\\theta \\right) z^k \\\\\n", "&= \\frac{1}{\\pi} \\int_0^\\pi \\left( \\sum_{k\\ge0}a_k (x + \\sqrt{x^2-1} \\cos{\\theta})^k z^k \\right) d\\theta \\\\\n", "&= \\frac{1}{\\pi} \\int_0^\\pi f\\left((x + \\sqrt{x^2-1} \\cos{\\theta})\\right) d\\theta\n", "\\end{align}$$\n", "\n", "Using this lemma on the sum, we get the following:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{1}{\\sqrt{4 \\epsilon^{2} \\sqrt{- \\frac{1}{4 \\epsilon^{2}} + \\frac{1}{64 \\epsilon^{4}}} \\cos{\\left(\\theta \\right)} + 4 \\epsilon^{2} + \\frac{1}{2}}}$" ], "text/plain": [ "1/sqrt(4*epsilon**2*sqrt(-1/(4*epsilon**2) + 1/(64*epsilon**4))*cos(theta) + 4*epsilon**2 + 1/2)" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import sympy\n", "from sympy import pi\n", "from sympy.abc import x,y,theta,phi\n", "\n", "epsilon = sympy.symbols(\"epsilon\", real = True, positive = True)\n", "\n", "u = x+y; Z = x-y\n", "X = u/Z\n", "\n", "f = lambda t: 1/sympy.sqrt(1 - 4*t)\n", "\n", "\n", "I = sympy.collect(f((X + sympy.sqrt(X*X-1) * sympy.cos(theta)) * Z), theta)\n", "\n", "# x = pE pW = 1/16 - eps^2 , y = pN pS = 1/16\n", "I = sympy.expand(I.subs({x: sympy.Rational(1,16) - epsilon**2, y: sympy.Rational(1,16)}))\n", "\n", "I" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This looks like an [elliptic integral](https://en.wikipedia.org/wiki/Elliptic_integral#Incomplete_elliptic_integral_of_the_first_kind), so we help `sympy` along by rewriting $cos \\theta$ in terms of $\\phi = \\frac{\\theta}{2}, d\\phi = d\\theta / 2$ \n", "\n", "$$ \\cos{\\theta} = \\cos{2 \\phi} = 1-2\\sin^2{\\phi}$$\n", "\n", "We thus have:\n", "\n", "$$ \\frac{1}{\\pi} \\int_0^\\pi f(\\cdots) d\\theta = \\frac{1}{\\pi} \\int_0^{2\\pi} f(\\cdots) \\frac{d \\phi}{2} = \\frac{1}{2 \\pi} \\int_0^{2\\pi} f(\\cdots) d \\phi $$" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/latex": [ "$\\displaystyle \\frac{2 \\sqrt{2} \\sqrt{\\epsilon} K\\left(\\frac{2 \\epsilon \\sqrt{1 - 16 \\epsilon^{2}}}{8 \\epsilon^{3} + \\epsilon \\sqrt{1 - 16 \\epsilon^{2}} + \\epsilon}\\right)}{\\pi \\sqrt{8 \\epsilon^{3} + \\epsilon \\sqrt{1 - 16 \\epsilon^{2}} + \\epsilon}}$" ], "text/plain": [ "2*sqrt(2)*sqrt(epsilon)*elliptic_k(2*epsilon*sqrt(1 - 16*epsilon**2)/(8*epsilon**3 + epsilon*sqrt(1 - 16*epsilon**2) + epsilon))/(pi*sqrt(8*epsilon**3 + epsilon*sqrt(1 - 16*epsilon**2) + epsilon))" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/latex": [ "$\\displaystyle 2.0000144264652$" ], "text/plain": [ "2.00001442646520" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from IPython.display import display\n", "integral = sympy.integrate(I.subs({sympy.cos(theta): 1 - 2*sympy.sin(phi)**2}), (phi, 0, 2*pi)) / (2*pi)\n", "display(integral)\n", "integral.subs({epsilon: 0.061912524105980984}).n()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With this, we are easily able to get over 1000 digits in a few seconds." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 294 ms, sys: 12 µs, total: 294 ms\n", "Wall time: 294 ms\n" ] }, { "data": { "text/plain": [ "mpf('0.06191395447399094284817521647321217699963877499836207606146725885993101029759615845907105645752087861371677762164603547703521657619786238920767692652100542462229364610602198310866645011155144480116905438357082333052235585247093496178675886579320261925322981842755398323065953839669902970037920145029333952863791551342841317381955154708160407987604794496625267390407048570909176364536425275876175019368190664354266513041076927589642937945735889870031373692899494980937812909549989657672408260042853941856806894557709937997149148640770482619263308825889098614336343956453277691625205180439236827297229820851545210184243583239147819016954127322193915516777984472434921418355899583976518392024997997265373859014398997622343683013391697219911388642490211986246607441523114844470587182633801852319326368640141921291336160036558158746187604490956929277980914586157049027359267817771067750751040233227627161401518091321963347432305885543319798172695687346868337236137394209763856503927615840572738408317472572053')" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import mpmath\n", "mpmath.mp.dps = 1000\n", "\n", "f = lambda t: sympy.lambdify(epsilon, integral, \"mpmath\")(t) - mpmath.mpf('2')\n", "\n", "%time mpmath.findroot(f, 0.0619, tol = mpmath.mpf(\"1e-2000\")).real" ] } ], "metadata": { "colab": { "authorship_tag": "ABX9TyPYuRb2DNXHchHX3uMFD7VX", "collapsed_sections": [], "name": "Problem 6", "provenance": [] }, "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.10.6" } }, "nbformat": 4, "nbformat_minor": 1 }