Skip to content

Instantly share code, notes, and snippets.

@pmineiro
Last active August 17, 2022 18:40
Show Gist options
  • Save pmineiro/3826bae079f6825b6b279c9b2ddd2ba6 to your computer and use it in GitHub Desktop.
Save pmineiro/3826bae079f6825b6b279c9b2ddd2ba6 to your computer and use it in GitHub Desktop.
Approximation of PolyGamma[1, b] for b >= 1 ... for the closed-form pdrop portion of the CS
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 8,
"id": "4bd09d76",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO29d3Rk5ZXu/byVo1IpdCtLrc4BOhMaDAYbcBNsGHtgnC72gLHH13NZXjMftsf2zHjufPZw567veq5NGAfsMYYBjA0NtHEiDKlzjqjVUiu0WjlVrlPv98ept1QqVUkVTlWdU7V/a7HsrpaOXrVU5zl7PzswzjkIgiCI0kRX6AMQBEEQhYNEgCAIooQhESAIgihhSAQIgiBKGBIBgiCIEoZEgCAIooQxFPoA6VBdXc1bW1sLfQyCIAhNceDAgRHOeU2iv9OUCLS2tmL//v2FPgZBEISmYIz1JPs7SgcRBEGUMJoQAcbYbYyxxycnJwt9FIIgiKJCEyLAOd/FOb+/vLy80EchCIIoKjQhAgRBEERuIBEgCIIoYUgECIIgSpiSEIFpXxCvnhjEyIy/0EchCIJQFZoQgWyrg3pGPfjCfxzA250jCp+MIAhC22hCBLKtDlq1xAmLUYfDvRMKn4wgCELbaEIEssWg12FDYwUOXZgrAuEwxz+/cgonB6YKdDKCIIjCUhIiAAAbmytwcmAK/pAUfe304DQef7MLz+zvLeDJCIIgCkfpiEBTJQJSGCdinvrf6hwGAJwYoE5kgiBKk9IRgeYKAMDhmJTQW52jAICTA1MIh3lBzkUQBFFISkYE6sosqC+34FDEHPYFJew9P4oquwnugISeMU+BT0gQBJF/SkYEAGBjcyUOXRgHABy8MA5fMIzPXNkCILOU0JQviEtTPkXPSBAEkU80IQJKTRHd2FyBvnEvhqf9eLtzBHodw2eubIVRz+Z4BakQlMK45/H3cP3/eh1/OHkpq3MRBEEUCk2IgFJTRKO+QO8E3np/BBubKlBlN2F5rRPH+9MTmB+81okTA1OocZpx/3/sx0/fPp/V2QiCIAqBJkRAKdbWl8OgY3j9zBCO9k/i6o7qyOtlODkwBc5TM4dPDEzi//6pEx/b2IDdf30Nblxdh3/YdRLf3X06l8cnCIJQnJISAYtRjzX1ZXj2QB84B3YsnxWBUXcAl6YWny0UCIXx1WeOoNJuwrdvWwObyYBHPrUZd25qwONvnsPwNM0nIghCO5SUCADAxqYKBEJh2E16XN4kp4fWNshpplTM4UdeP4fTg9P4fz+2HhU2EwBAr2P4wrXLEObAb08M5u7wBEEQClN6ItBcCQC4ot0Fo17+9lcvLQNjwPH+xc3hV45dxNUdLty4pm7O6yvqHOiodeDlowPKH5ogCCJHlJwIbG6pBGPAdStroq85zAa0ueyLRgJBKYyukRmsb6iY93eMMexcvxR7zo9haJrKRgmC0AYlJwJNVTa88pVrcM+25jmvr6kvW7RMtHvEjaDEsaLOkfDvd25YCs6B3x6nlBBBENqg5EQAkNM/Bv3cb31tfTn6J7yY8ASSft7ZSzMAgBV1zoR/v6LOiRV1Drx09KJyhyUIgsghmhABpZrFFmJtfRkALDhW+uylaegY0FGbOBIAgJ3r67Gve0yxTuIxdwC/PX4R33npJD7x2Lv40X91pVzKShAEsRiaEAGlmsUWQojA8QV8gbOXptFcZYPFqE/6MTs3LAHnwO5j2UcDnUPTuOZ7f8IDvziIX7zXg9EZP/7p5VP40pMHMeMPZX19giAIQ6EPoBZcDjMcZgMGJ5PX+Z+9NJ00FSToqHVi1RInXj52Ef/t6raMzxOS5H4Ek0GHn967DZc1lcOk1+Hf/6sL3919GmcuTeNHn9mC9prkUQlBEMRiaCISyBcWox6+mKUzsfhDErpHPYuKAAB8aE0d9veMwxPI/Gn90TfO4UjfJL7z0XXY1lYFs0EPxhjuv3YZnvzLKzDuDuDBZ47QCGyCILKCRCAGq0kHXyCxCHQNuyGFOVYsWVwEVi8tA+fy52TCyYEp/J8/vo9bNyzFrRvq5/39lctc+MbONTjSO4Fd1JdAEEQWkAjEYDXq4Q0mFoGzl6YBIGl5aCzLI8bx+0PTaZ8hJIXx1WePoNxqwnfuWJf04+7c2IA1S8vwL789A1+SMxMEQSwGiUAMi4mAXsfQVm1f9DotLjsMOob3IyWl6fBu1yhOXZzCN29djUq7KenH6XQMf7dzNfonvHjine60vw5BEARAIjAHi1EPb5J00NlLM2irtsNsSF4ZJDAZdGittuP9ofRF4E+nh2Ay6PChuLEUibiqoxo3rKrFD/7UidEZGlxHEET6kAjEYDXpk6ZW5Mqg1Ctxltc60JmBCLx+ZhhXtrtgM6VWuPW1j6yGJyjhh6+fS/trEQRBkAjEkCwd5A1IuDCWWmWQYHmtAz2j7rTy9edH3Dg/4sYHV9Wm/DkdtQ7ctLYOLxwegESVQgRBpAmJQAyWJCJwbngGnCcfF5GIjjonwly+safKn04PAQCuX5m6CADAR9YvxciMH/u7x9L6PIIgCBKBGGRPIDzv9TODqVcGCUSFUDopoddOD2FZjR3NLlvKnwPIomE26LCbBtcRBJEm1DEcg9WY2BM4OzQNk16HFtfilUGCtmo7dAwpm8Nufwh7zo/iv13VmvLXENjNBly3sga7j1/Et25dA52OpX2NVDkxMIkfvnYO3qAEo57BYtTjE1uaoqs6CYLQFiQCMVhNuoTpoLOD02ivsUeX0KSCxahHi8uOzhR7Bd7qHEFQ4rg+DT8glo+sX4pXT1zCod5xbG6pyugaC+ELSvj+H9/HY292ocxiQEOlFSGJY3jajxcOD+DOjQ34xs7VcDnMin9tgiByhyZEgDF2G4DbOjo6cvp1rEY9pDBHUArPueH3T3jRmkYUIFhW40i5V+C100Nwmg3Y2prZDfyDq2phMujw8tFBxUXg4qQXn/zRHnQNu/HxzY34xs7V0dWavqCEH7zWiUffOIc/nh7CD/5iU3R3M0EQ6kcTnkA+pogCiE4HjY8G3H4JDkv6erm8zoHzI24Epfk+Qyycc7x2ZgjXrKhOK9qIxWkx4trlckpIyXlCnHP87XNHMTjpw88/tw0Pf/yyqAAA8r/ZVz+8Eq985RrUlZnxxScPoGs4/dJYgiAKgyZEIF9YTbIIxM8P8gRCsJkWbxKLZ3mtA6EwR8/owhVCZy/N4NKUH9elWRUUz0fWL8HFSR+O9E1kdZ1Ynt7Xi/96fwRfu2UVrl1Rk/Tjltc58ePPboVRr8N9P9+PKV9QsTMQBJE7SARisCaLBAIS7Ck2b8WyvFYuKV0sJXR+RP77NUvL0v4asdywug5GPcMrCuwyAIC+cQ/+6aWTuGqZC5/c3rLoxzdV2fDDT25Cz6gH/+Ppw9S3QBAagEQghkQiEJLCCITCKXfwxrKsVvYRFqsQ6hv3AgCaKtMrDY2n3GrEtrYqvN05mtV1ADkN9P/86igA4Ht3bUi54uiKdhe+fdsa/On0EB57k7qYCULtkAjEYImkfGLnB3kigpBJOshmMqCx0pqSCDjMBpRZs/fpNzRW4OylafiT7EVIld8eH8TbnaP4+s7VaKpKT5w+dUULblxdh0deO7fgzmaCIAoPiUAMiSIBjz8iAub0RQCQfYH3Ly1cJto37kVjpRWMZV/fv76hHKEwjza4ZcpzB/qwpMyCu7c2p/25jDH8zU0rMRMI4dE3urI6B0EQuYVEIAYhArENY2I7WCaeACAbpl0j7gXz433jHjRUWDO6fjzrG+QKquP9UxlfY2TGj9fPDuOOjfXQZ9h4tnKJE7dfVo8n3jmPoWlfxmchCCK3kAjEEC0RjRkd4Qlkng4CgIYKKwKhMMYXSIv0T8iRgBI0VlpRbjXiWP9kxtfYdUQeRnfnxsaszvLgjSsQlDh++Bp5AwShVkgEYkiUDnL75UggE2MYAJyR/oIZX+J9w5PeIKZ9ITRmaQoLGGNY11CG41mIwK8P9WNtfRlWprBKcyFaq+34+OZG/HLPBfSNe7K6FkEQuYFEIAaLSf7nmOMJBLLzBBzmiAj4E4tAf6QyqEGhSAAA1tWX48zgNAKhhZvUEtE5NI2jfZP42MYGRc7y329YDgB4jLwBglAlJAIxiEjAn0AEMvUERKfxdJJIQDwhK5UOAoB1DeUISOHoXuR0eP5gP3QMuP3y+QvuM6Ghwopb1i/BrqMDGYkSQRC5RROzg/LFrCcQkw4KiHRQZpFAmcUIAJhO0kHbPxGJBBQyhoFYc3gS6xpSH7URDnO8cHgA166oQa3Toth57ri8Hi8cHsCbZ4dxYwprM5XEG5Cwv2cMb3WOYH/3OCqsRqxY4sTKOieuX1mLcpsxr+chCLVBIhCDUa+DUc/iSkSzE4HF0kF9415YjXpULbBUPl2aq2xwmg04PpCeL7C3ewz9E1787c0rFTsLAFyzvAaVNiNeODKQVxH4zaF+fOPXx+AOyGOv1zeUo2/cizfODiMU5qi0GfHVD6/EPduaM66CIgitQyIQR/x2MdEsZjdnlw5KLgIeNCjUIyDQ6RjWNpThWJplou+cG4WOATeuVvZGbdTrsHPDUjx3oA8z/lBUGHOFLyjhOy+dxJN7LmBrayX+6voObGuripr7gVAYxwcm8b3dp/F3vzmOJ/dcwMN/tiGtqIkgigXyBOKIXyzj8UvQMcBsyOyfStzwknkCSpaHxrK+oRynLk4tOsE0lpMDk1hW48hY8Bbio5c3wBcM4/cnc7v9bMITwCceexdP7rmAL1zbjqfuuwLXraydU91lMuiwqbkST99/Bf7vX2zEuDuAex5/D4d7lRu8RxBagUQgDqtJP88TsJkMGT+pW4x6mPS6BYxhr6J+gGBdQzkCoXBa6y2P90/l7Gl4c0slGiut+M2hgZxcH5DnHT30q2M4dXEKj316M772kdUwLDCamzGGWzfU4/kvXYVKuwmf/vEeHFVwAitBaAESgTiscekgb0DK2A8QOCwGzPjnG8Mz/hAmPEHFegRiETfzVJvGhqf9GJzyYW19dpNMk8EYw+2X1eOtzhGMzPhz8jWe3teL354YxN/ctBI3rV2S8ufVV1jx1P1XoMJmxKd+tAfH+jLvsSAIrUEiEIfsCcymUNwBKev0iMNsSNgsJnoEcpEOanPZYTfpcSJFETgRMZHX1ucuL/7RjQ2QwhwvH1Vm1HUsnUMz+IddJ3DN8mr85Y72tD+/ocKKp+67Ak6LEQ/84kDSai6CKDZIBOKwGHVzlsp4/JktlInFYTYkNIb7J+QeASUbxQQ6HcPa+nKcGEjNHBYftyZHkQAArKiTSzNfPaGsL+APSfjKU4dgMxnwrx+/LOWx1/E0Vtrw/Xs24uKkF//8yilFz0gQaqWgIsAYW8MYe4Yx9ghj7M8KeRZBfDrIneFWsVicFgOmEkQCfTmMBAB5yctApA9hMU4MTKLFZUO5Nbd181cuc+HQhYm0DOvFeGZ/H05enML37tqA2rLs+hs2t1Tivmvb8dTeXrx+ZkihExKEeslYBBhjP2GMDTHGjse9fjNj7AxjrJMx9tAil7kFwL9xzr8I4DOZnkVJrKZEnkB26SCnJXE6qG/cC5NBh2q7OavrJ6PGacbwjB+cL77h63j/FNblMBUk2NpaBW9QwskUI5TFCIc5fvr2eVzWWI4bV2e3nlPw4I0rsLzWgYd+dQyTXkoLEcVNNpHAEwBujn2BMaYH8APIN/c1AO6JPO2vZ4y9FPdfLYD/AHA3Y+xhAK4szqIYFmN8dZAEe4ZzgwRJ00HjXjRWWDNOXyxGrdOMoMQx4Vn4RjbpCeLCmAdrG3KXChJsaa0EAOzrHlPkem+cHUbXsBuf29GmWK+FxajH//7E5Rie8eOfX6a0EFHcZCwCnPM3AcS/k7cB6OScd3HOAwCeBnAH5/wY5/zWuP+GIv/9FYCHAIwk+jqMsfsZY/sZY/uHh4czPW7KWI36OVu5PP4QrMYsjWFLYhEQjWK5osYpRxhD0wtX45y4mHtTWFBXZkFzlQ37u8cVud5P3j6PujIzPrJ+qSLXE6xvLMdnr2zFcwf70DPqVvTaBKEmlPYEGgD0xvy5L/JaQhhjrYyxxwH8HMDDiT6Gc/4453wL53xLTU2NoodNhDUuEvAEs48EnBYjpn3BeWmZXDWKCWojIjC8mAhEOotzVR4az5bWSuzrHkspTbUQZy9N47/eH8FnrmyFcYF+gEx54APt0OsYHnmd9iEQxUtBjWHOeTfn/H7O+Sc5528V8iwC4QmIG5THn70n4DAbEJQ4/DFTNL0BCSMzgZz0CAiESbrYZq8TA5NYWm5BtSM33kQ8W1urMOoO4PxIdk/YP337PCxGHf5iW/orMFOhtsyCu7c24VcH+6KD/gii2FBaBPoBNMX8uTHymmawGPUIcyAghREIhRGQwrArUB0EzJ0flIvpofHUpBgJHB+YyksqSLA14gtkkxIacwfw/MF+3LmpEZUKDt+L5wsfWAYAeOwNigaI4kRpEdgHYDljrI0xZgJwN4AXs70oY+w2xtjjk5O57+SM7hkOhKNpIasCfQLA3O1iURHIYTrIYTbAZtIv6Al4AiGcG57BujyYwoJlNQ5U2oxZmcMvHR2APxTGZ69sVe5gCWiosOKuTY14el8vhqbUsyvZF5Sw9/wY3jw7jM6h6aQDCgliMTLOczDGngJwHYBqxlgfgG9zzn/MGPsygFcB6AH8hHN+IttDcs53Adi1ZcuW+7K91mKIG743KIFDTgll2zHsjO4UmH2jjkRuzCJvnytqnOYFReDUxSlwnh9TWMAYw5bWKuzvyTwSeOv9ETRWWrNegZkKX7quA88e6MNjb3bhm7euyfnXS4Y/JOGnb3fj9ycv4VjfJAJxvRbLauy49+o23LWpMesHF6J0yPjuxjm/J8nrrwB4JeMTFZjYPcNSWBYBJTqGAWA6Zn7QqFu+MbtynIevdZoxvIAnIOr182UKC7a2VuL3Jy9heNofTVulihTmeLdrFDsVrghKRrPLhp3rl+LZ/b3425tXwmzI/w32nXMj+LtfH0fXiBuXNVXg3qtbsaW1CmUWAwanfOif8GL3sUH83W+O43/97gzuu6YdD3xgGe1JIBaF9gnEYTFG9gwHJITC8pOWEs1iwNx00OhMAGaDLmu/YTFqnRacGkzemHVx0geDjmFpuXKbxFJhS2sVAOBAzxhuXpfezfxY/ySmfSFc3VGdi6Ml5M5NDXjxyADeODOMD6cxnC5bQlIY33zhOJ7a24vmKht+/rltuHZF4iq5L35gGfZ1j+PxN8/h4VfP4L2uUfyfuzcqurCIKD40MTson56AJSYSmN0vrJAnEJO3HZ7xo9phVnSZTCJqnGYMTyVPB014g6iwGXN+jnjW1ZfDbNBhXwbm8NudckvJVcvy1194dUc1XHYTXjiSu1HY8XDO8Y1fywJw/7Xt+N2D1yYVAEBOs21rq8KPPrsV371zPfacH8Nt//YWjccmFkQTIsA538U5v7+8PPd566gxHJTgEfuFs/YE5i+WGZ0JwOXI/RNajdOMaX9oTu9DLBOeQM7nBSXCZNDhssYKHLqQmQisXlqW81RaLGI72h9OXsqbCfu9357Bf+7vxVc+2IGvf2R19AElFe7e1oznHrhS/v+Pv4fjKU6TJUoPTYhAPokawwEJbr9CkUCCEtFRtx+uPITpi5WJTniCqLQVJl2wrNaOnlFPWp/jDUjY3z2OHR35nzJyx+X18IfC+J3CU1AT8e9vduHRN87hk9ub8eCHVmR0jQ2NFfj1l65Cpc2Ee5/YR70OREJIBOKINYaVKhE1G+ZvF5Mjgdw/yUa7hmcSm8MTHjkdVAiaqmwYdQfSerLe3zOGgBTOqx8g2NQsb0d74XBuU0JH+ybwz7tPYef6pfjHO9ZllaqrLbPgp/duhS8o4d6f7sUU7Ukg4tCECBTCE/AFJbgj6SB7lsYwMHe7GOc8b+mgWmekaziJLzDhCaCiQJFAS5UdANKazfNW5wiMejn3nW/ysR0tHOb49osn4LKb8d271itS3bOizonHPrUZXcNufOkXBxEOZzeugyguNCECefUETLGegBwJ2LKcHQTIvoCIBKb9IQSkcM5GSMey2BC5CW8QFQXwBACgxSWPzLiQRkronc5RbGyuzLpiK1PuuFzejvbKMeW3owHA84f6cejCBB66ZVW0v0QJruqoxnc+ug5vdY7gyb0XFLsuoX00IQL5JDYd5PaHYNAxmBQYTha7YnJ0JgAAeYkEquwm6HUsoSfgD8lCV6h0ULMQgbHURGDcHcDxgUnsKEAqSLByiROrljhzkhKa8gXx3d2nsbG5AnduTDp3MWPu3tqEqztc+Jfdp3FJRd3PRGEhEYgjWiIaCMMTkGA16RUpn3SYDZj2CxHIT6MYAOh1DC67KeEQucnInoFCpYPKLEZU2IzoSVEE3u0aBecoiB8Qy4fX1OHQhXHF9xB//w/vY9Ttxz/cvjYnOyYYY/ifH12PgBTG37+YdSM/USSQCMSh1zGYDLpIn0BIET8AmLtdbCQSCVTnIRIAgNqyxKMjxqMiUJhIAABaqmwpp4NOXZyCjgHrG/I34iIR29tdCHPgQBZjL+IZmPDiiXe68edbmrChsUKx68bTWm3HX9+4HLuPD+L3Jy/l7OsQ2kETIpBPYxiQU0KyMSwp4gcAkZ0CEWNYjIzI1+jmWqclYTpowiOLUaFKRAGg2WVPOR3UMyov4TEZCvtru7G5AgYdw97zymxHA4BfHehDKMzxpes6FLtmMu67ph2rljjxrReOJ+0fIUoHTYhAPo1hQB4d4Q1Ikf3CyohAIk8gXzffGkfiSGAisj+3EM1iguYqK/onvCktnu8Z80QrigqJzWTA+sZyxUQgHOZ45kAvrmx3RX2SXGLU6/D3t6/FxUkfnj3Qu/gnEEWNJkQg31iN+qgxrFQVilgxKZeH+lFuNebtiba2zIzRGX90IJ4gGgkUcLZMS5UdUphjIIVGpguj7rzcJFNhW1sVjvRNKPIk/V7XKHrHvPjzrU2Lf7BCbG+rwuaWSjz2RldKAkwULyQCCbBERMATkBQb8Ba7XWzEnZ8eAUGN04wwn01DCcQC+kKViAKzFUKLdQ5P+YIY9wTRUqUOEdjeVoWgxHGoN3tf4Jn9vXBaDLh5Xf4G0zHG8MAHlqF/wouXj+am3DUTQlIY+7vH8PCrp3Hfz/fjq88cwT+9dBI/fuv8osuRiMygKaIJsJr00WaxZrMyN52ymPlBozP+vPQICETX8NCUP9o8BsjGsFHPFEt5ZUJLimWiwjxuUUkksLmlCowBe8+P4aplmVcrTXqD2H18EJ/Y0pTWbCAluGFVLZbXOvDI6+dwx+X1eR8iGEsgFMYPXuvEE+90Y9IbhF7H0F5th9sfwrgnCG9Qwvd2n8Ztl9Xjczta87r/otghEUiAWDbvDUiwKfTGjJ0fNDoTQEetQ5HrpkJ0flBcl+ukV+4WLuSbv85pgcmgW1QERKTQrAJPAJB9lNVLyrL2BV483A9/KIxPbMlfKkig08nRwFefPYLXzgzhg6vq8n4GADg9OIWvPnMEJwamcPPaJbjtsnrsWF49x6s6NzyDn73TjecO9OFXB/vwlRuW43/csDwnpbSlhibSQQWpDgrJnkC2W8UEDrP8Cz3jC2E0z+kg8fQfP1J63F24bmGBTsfQVGlddHREz5j892rxBADZFzh4YRyBUOY59Wf292H10rK8rveM5fbL61FfbsEjrxdmh/Iz+3tx+7+9jcFJHx791GY8+unN2Llh6bxihWU1DvzjHevw7tduwMc3N+L7f3wfD/ziAK3VVABNiEDeq4NMciTgUbg6CAAmvAGMewJw5TEdNDs6Ym7D2IQ3UNDyUEGLa/FpohdGPah2mKL/jmpge1sVfMEwjmU4pvnc8AyO9U/iE1saCxaNGfU6/OU17djXPZ73cdN/PHUJD/3qKLa3V+F3D16bkidSbjXiX/5sA7516xr88fQQ7vzh26ra/axFNCEC+cZq1GPKF0IozBUTAbFT4MKYB5znr1EMkI1up8Uwz1ib8ARRXsBGMUFzlQ29Yx5wnnywWc+oB80qMYUFYohdpimh97pGAQDXr6xV7EyZ8LGNDTDoGHYdzd/CnKN9E/jyLw9hXUM5Hvv05rS65xlj+NyONvzs3m3oHfPii08ezCoaK3VIBBJgNeox5pbLJ5UqERUi0D0ipzXyuRAFkM3h+F4BeZeAOkTAHZAwGvk3T8SFMQ9aXOrwAwQuhxkdtQ7sOT+a0efvPT+GGqe54GZ3pd2EHcur8dKRiwsKsVL0jnnwuSf2weUw4cef3Zrxe2zH8mo8/PENONAzjr/fRWMwMoVEIAFWkz5aU29XqGNYpDG6I2mPfCyUiaXGaZ4fCXgLN0Y6lpZFykT9IQkDk17VRQIAsLW1Cge6xzO6ee47P4ZtbVUFNeYFt22oR/+EF4d6c7uKMhzm+MrThxCUOJ64d1s0VZkpt26oxwMfWIZf7rmAp2g6akaQCCQgtlRPyWYxYHZ2fv4jAQsuxXgCvqAEXzBc0G5hgRCB3iQVQr1jXnCunvLQWFYtcWLaH0q7hr1v3IOBSR+2teZ/L0IiPrS2DiaDDrtyvEP5hSPyqOxv3rpGsQq5v7lpJa5dUYNvvXCc1mhmAIlAAizG2X8WpTwBs0EPk0EXfdrNpycAyCmXgQlfNHcqGsXUYAw3Vi4cCVyIVAapUQRaq+UU1fmR1BfjAMC+btlH2KoSESizGHH9yhq8fPTivM5ypfAEQvje7jO4rLFc0VHZeh3D9+++HOVWI/5x18m8pLSKCRKBBFhzEAkAgNNsgD8UhkHHUKbgwpBUaKuWxzP0jss32vHIyIhCThAVWIx6LCmzRMtA41Fbj0As7RmKwN7zY3BaDFi5xJmLY2XEbZfVY2jar+hgvFgeff0cBqd8+NZtaxSv76+wmfDgh1Zgb/cYXj1B01HTQRMiUIg+AYFSngAwmxKqspvy3uQinliFMT2hgjHSsTS7ko+U7hn1wGbS5z16SoX6CitMeh3Op7EiE5BFYEtLpSLrI5Xig6tqYTPp8VIOqoT6xj147M0u3H5ZPTa35Cb6+ao2eOIAAB6BSURBVPMtTVhe68B3d5+iaqE00IQI5LtPIHaxvJIjFYQ5nG8/AJj/xDrpjUQCVnXcWJeUWZLu7b0wJpeHqsFAjUevY2h22XB+OHURGJ3x49ywG9vaXDk8WfrYTAbcsLoOu48PKj5U7l9/dxaMAQ/dskrR68Zi0Ovw9Z2r0T3qwS/e68nZ1yk21NN5oyJyYQwDs2WihXiirbSbUGEzRkVALJSptKsjErAYdfAneXrrGXXndcxGurS67OhOIxLY1y0PndvWVpmrI2XMzvVLsevIAPZ1ZzcTKZYxdwAvHR3Ap65oQX2FVZFrJuO6FTW4Znk1vv+n93HXpsa898H0jLrx2Jtd6B3zoH/CC49fwl/fuBx3b21S5UMMoJFIIN/MSQcpKAJidES+lsnE0+qyR0VgdoKoOiIBs0GfUATCYY7eca/qegRiaa+xo3vUg3CKhure82MwG3RY35C7DWKZcuUyFxgD9p1Xbmvabw71IyjxvIzKZozh6x9ZjQlPED97tzvnXy+WAz3j+NgP38FvDvVjyhfCqiVONFZa8bXnj+HLvzyESa+y60iVgiKBBMSmg6wKpoNEJJDvHgFBe7Ud70a6VCc8AZgMujmVUIXEbNDBH5w/m39wSq5oUmOPgKDVZUcgFMbApDda6bQQ+7rHcHlTRcE3pCWi3GrEyjpntHopWzjneGZ/LzY0lmPVkvzMR1q9tAw7Oqrxn/t68eXrO/Liv7189CIefOYwlpZb8NwDV6K9Ro5cw2GOx97swr/+7gwO907g6fuvQJPKfpfV91uoAkQkYNQzRd+ohfQEALlC6OKkD96AFO0WVkuIak6SDupR2QjpRLSlUSE04w/hxMBkdOSEGhGD8UIK+ALH+6dwenAaH8/zlNS7tzWhf8KLtzpHcv61Xjl2EX/1y4NY31CO5794VVQAAHlA4hevW4ZnH7gSY+4AHn71TM7Pky4kAgkQnoCSfgAQEwkUqMolWiE06sa4J6CaVBAgp4NCYT7vxhPtEVBheaigvWZu5dVCHO+fRJgDm1vU5wcItrZWwROQcGJgKutrPbO/F2aDDrdfVq/AyVLnQ2vqUGkz4ul9ue0invYF8e0XT2B9Qzme/MvtSR/wNjZX4t6rW/HikQGcVODfVUlIBBIgUkBKbRUTOApoDAOzT6zdI25MeIOqKQ8F5HQQAATiRGBgwgfGgPoKS6JPUwW1TjNsJj26UhAB0RXdqmKPQzSwZZsS8gUlvHC4H7esW5L3znSzQY+7NjXi9ycvJa06U4L/7w/vY2TGj3/66LpFlwJ94dplKLMY8L9/r65ogEQgASIdpKQfAMjNYgDyOkY6FiECXSNuTHrUKQL+4FwR8AYlWAx6GPTq/VVljMkVQimIQP+EF4wBS1UsakvKLWiqsmYtAq+eGMSUL1SQhTmAnBIKShy/OtCXk+ufGZzGE+904+6tTbisaXGTv9xmxBc+sAx/ODWEAz3KGe/Zot53VgERZqlSC2UEK+qcKLcaC2Zy2s0G1DrN6B5RXzpIPEXF+wK+oASzSszrhWirtqfkCfSNe1HrNMNsKNxKz1TY2lqF/RkOxhM8d6APjZVWXNFemH6IjlontrRU4j/39So+SoJzjm++cBxOiwF/c1PqvQ/3Xt2KaocJD796WjXjLdT/7kL+O4YtBuEJKPtG3d7uwpFvfxiVBaoOAuSbVZdIB6mkRwBA9EbvD82tEPIHw9Gfh5ppq7ajd9y7aJNV/7gXDTmulVeCba1VGHUHcC6NJrhYfEEJe7rGcPPaJQVdAXn3tmZ0jbgVH4XxyrFB7D0/hr+9aRWq0ng/20wGfPn6DrzXNYb3unIzniNdNCEC+e4Y1ukYzAadoj0CaqG9xo4zg9MIhMKqigTEk7EvLh3kC0mqKWNdiFYxm2mRXcn9E140pFBGWmi2RqqX9meYEjrSO4GAFMb2AkUBgp3rl8Jq1OPlYxcVve4v3utBq8uWUe/D3duaYTLo8MdT6phxpP53V4GwmvSKewJqoNVlj+5lVcNCGUHUE4iLBHxBSfWpEyDGdF+gczgc5rg4qY1IoL3aDpfdhL0ZioB48t7aWtgqKKtJj+3tVYqWivaNe/Bu1yju2tSY0ewni1GPjU0V2JOjQX3pQiKQhKZKm6pr0zNF3KwA9QyPA2YjgXhPwB8KayISiJruC6RPhqb9CEocjZXqFwHGGLa0VmZsDu/tHsOqJU5VLC3a0VGNrmE3Bia8ilzvN4f6AQAfzWIc9va2KpwYmMS0r/BdxOp/dxWIZx+4Eg/euKLQx1CcuSJQ+DeoIOoJxKeDghLMi5TeqYFKmxHlVuOCkUD/hJwqatCACACyOdw75sXgZHqL3INSGAd6xlXTELdjuTwDSYlogHOO5w/2Y3tbVVadv9vbXQhzYL8KqoRIBJJgMaq7LDFTml02iCZhdUUCydJB4ejfqRnGGFoXqRDqG5efRBs1kA4CgC2RfoHDaa6cPDEwBU9AUo0IrKxzotphxtsKiMDh3gl0jbhx16bGrK6zqbkSBh3DHhWYw+p/dxGKYjboo+kINRrDidNB6o8EADmP3j2S3BgWIqCVSCDaCZ32rgR5PpVaRIAxhh0dLrzdOZJ1WebzB/thNuhwy/olWV3HatJjQ2N59N+qkJAIlCCiW1ULkYA/KGlGBFpcNvRPeOd9D4L+CS8qbUbFx5HkijKLES67KboXO1X2nh9De7UdtU71NMRd3VGNkZkATg9OZ3wNf0jCrqMDuGntEjgV2Ay4vd2Fo32T8ARCWV8rG0gESpCVkaY1Nd1cF/QENJAOAoCKyGgEtz+JCIx7NRMFCFpctgWjm3jCYY6958dUEwUIhC+QTUrotdPDmPAEcecmZfYjb2+rQijMcbAnvXSb0mjj3UUoyn//4HI8ff8VhT7GHCwarw4CZjvM3f7ET3b9E140Vmir4qzVZU8rEjhzaRpTvpDqRGBpuRXLauxZmcO/OzmIKrsJOzqUWbazuaUSOoaCp4S08e4iFKXcZsTqpfmZ7Z4qyTqGfZHZQVogKgIJwnvOuUYjATsGJn3wJdj1kAjRH6A2EQDkUtE9XWNJ03WLcejCBDa3VCpWMOK0GLGuoRzvFbhfgESAUAWmyBtrfsdwWBOzg4DZMSOJ0kHjniC8QUkTjWKxtFbLkcuFRTqhBXvPj6GhwprScp18s2N5DbxBCYcupJ9+GXMHcH7EjU3Nyja/bW+rwuHeiZRFNhdo491FFD0GvQ4GHZvzlBaSwpDCXDORgFgalMjo69dYZZBAFBGkMiEVAA5eGMeWAncJJ2N7exX0OoZ3MkgJHbog1/NvalZ2Jei2NhcCoXDaZbhKogkRyPcAOaIwyCsmZyMBX8QfUJOBvRCi6ieRJ9A3HmkU01okEBEBseFtIXxBCRcnfeiI2aylJsosRrRX23EqgwqhgxfGYdAxbGhUVgTEWI0jJAILk+8BckRhMBvnLpsXIbJW0kF2c/J0UH9kZIEWRkbEUm4zosK2cCe0IPo9Vqn3e+yodeDc0Ezan3ewZwKrl5YpPk+swmaC1ajH8HTuFt8shjbeXURJYDbo5qSDhCBoJR1kXyAd1DfuhcNsyPuGLSVocdlTigTEBNUmFfoBgo5aB3rGPGmZwyEpjCN9E4qnggQuhwmj7kBOrp0KJAKEapBFQMORQCQdNJMkEmiosIKxws3Wz5RWly3lhTkAVGkKCzpqHZDCPK3ehzOXpuEJSNiUo73Q1Q5zTldgLoY23l1ESWA26Od6AkIENBIJWIw66FhyY1hrprCg1WXHwGTyTmhB77gHJoMOtc7CrE9NhWURv6IzjZTQwUg1kdKVQYJqhwkjMxQJEATMxiTpII1EAowx2E2GpJ6A1kxhQWu1DZwDvWMLj2LuG/OiscJa0E1ii7GsxgHG0hOBQz3jqHaYc+bnVDvMGKVIgCCSp4O0Uh0EADazfl510LQviElvULORQEu0QmjhlFDfuEf136PVpEdDhRWdw+lEAuPY1FyRs1Se8ATC4cLsHCYRIFSDJa46SKSGtDI7CJB9gfiOYa1WBgmivQKLmMO9496sZuzni45aR8qRwOiMH92jnpz5AQDgspshhTkmvYVZMKOddxdR9JgNujmdkyI1pKVIwG42wBOYmw4SG63qNZoOqrQZ4bQYFmwYc/tDGHMHNCF0HTUOdA3PQErhyftQjv0AAKiOeCiFModJBAjVYDbE9wloLxKwmfTRHc6Ccbf8hOeyq2d/QzowxtDqsi/YKyAqg9RcHiroqHXAHwpHu7gXYrZJLHc9StWR34tCmcPaeXcRRU98n4AWPQE5EpgrAlORPbJlCsygLxSt1Qv3CoiOaE1EArWRCqHhxTuH3x+aQXuNPae/gyISGHVTJECUOGbj3LERfo2NjQAiIhBXHTTtk0XBadHGMplEtLps6Bv3IBA36lsQbRTTiCcApFYhNDjpy3kaT0SIIwXqGiYRIFTD/HSQ6BPQzq+pPUE6aMobhN2k7Z3VLS47wnzW5I6nd9wLq1GviZRXhc2EaocpJRG4OOnD0vLcbkirtJmgYyhY17B2fyuJomN+Okh7kYDNNN8YnvIFUabBcRGxtLrkJ/xk5nDfuAeNldrpiF5Ws3iFUCAUxsiMH0vKchsJ6HQMVfbCdQ2TCBCqQfQJiGXg/pAEo55Br+Lmo3gcZj3cgdCcheZT3pCm/QAAWBpJiVya8iX8+94xbZSHCkSZ6EKL58X3mutIAChs1zCJAKEazEY9OAeCkvzG9AXDmhkZIbCZDeAc8MaUuk75gpr2AwDAFonGvEmWn4hIQCt01Dow5QtheIGn78GICCzJiwhQJEAQ0dy/SAn5QpJmRkYI7Am2ixVDOkiMUE4kApPeIKZ8IU2UhwpSMYcvTuYvEnA5TBilSIAodczGucvm/RqMBBKNk5bTQdqOBMwGHRgDvIH5IiDKQ5tUvEcgnuW1TgCLiEDEBF+ahyY/igQIArORgKgK8oUkzYyRFtii46RnRWC6CCIBxhisRn1CERCD5dQ8QjqeujIzHGbDopGA02yIrg3NJS6HCZ6AlHACba7R1juMKGpm00EiEpA0s1BGILaLiQohzjmmfNo3hgHIIpAgHRSNBDQkAowxLCm3LLjRa3DSlxc/AJAjAQAFSQnlTQQYY+2MsR8zxp6Lec3OGPsZY+zfGWOfzNdZCHUiUj+iYcwfCmvOE4jfM+wJSJDCXPPGMCCX6iZOB3nhNBtQZtXW9+i0GKKNfIm4OJVPERCjI/KfEkrpHcYY+wljbIgxdjzu9ZsZY2cYY52MsYcWugbnvItz/vm4l+8E8Bzn/D4At6d1cqLoEKmfqDEclDTnCYjUgTCGoyMjNJ4OAuS5SIkigd4xDxqrbJrpERCUWYyY9iWf3Dk46c2LKQzIk0QBdUcCTwC4OfYFxpgewA8A3AJgDYB7GGNrGGPrGWMvxf1Xm+S6jQB6I/8/9aWfRFESnw7yBbUYCUSqgyK53Smv/L9FkQ5KIgJ9415NlYcKFooEglIYQ9N+LCnPz/eVaJLos/t7cbh3IudfO6V3GOf8TQBjcS9vA9AZecIPAHgawB2c82Oc81vj/htKcuk+yEKQ8lmI4iWaDhKeQEjSVLcwEFMdFEkHzUYC2kqVJMJi1M/rhga01yMgcFqM0Z9PPMPTfnCen/JQYHZ+kBgd4Q1I+Pqvj+Fn73Tn/Gtnc+NtwOxTPCDf0BuSfTBjzMUYexTARsbY1yIvPw/gLsbYIwB2Jfm8+xlj+xlj+4eHh7M4LqF2opGAqA4KhjU1NwiYNYbdkZvldBFMEBXYTPo5+x4AICSF4Q5IqLKpf2ZQPGUWA6aSRAKiRyBfnoDFqIfTbIga1YcujCMocYzlYZ5Q3h5POOejAB6Ie80N4N5FPu9xAI8DwJYtWwqzf43ICxZjfDpIe5GASa+DQceixnA0HVQEnkCiElFPRBREM5mWcFoMCITC8Ifme08XJyM9AnkSAWB2zSQAvHdeTrzkQwSyeczqB9AU8+fGyGsEkRHz00FhzYkAYww202zaRKQbiqE6yGqanw4SY7PteailVxpnJDpL5AsMRruF85fmqnaYo+Ok3+saBaB+EdgHYDljrI0xZgJwN4AXlTnWXBhjtzHGHp+cnMzF5QmVIKqDos1iQUlz6SBArhCaiUYCRSQCxvnpIGGA2zQaCQCJReDipA82kz6vnd5yJOCHLyjhcO8EGMvPoplUS0SfAvAugJWMsT7G2Oc55yEAXwbwKoBTAJ7hnJ/IxSE557s45/eXl+duxRtReGIjAc45/KFwdJSElrDFbBeb8oVgMeo0V+qaCGsCYzgaCZi0J3KzkcB8c1g0iuWz7NXlMGN0JoDDvRMIhMLY0lIJXzCcsDdDSVL6yXHO70ny+isAXlH0RETJEjtAbnarmPYiAbtJP9sn4A0WhSkMzPYJcM6jN8doJGDWnsiVLRgJ5K9HQFDtMGPME8A7nSNgDLhp7RLs6x7HqNuPRlPuurG19w4jipbZ6qBwtGtYi0/QdrMhagxP+0JFYQoDgMU017MBZgfKFWUkkONlMvFUO0zgHNh9fBBrlpah1WUHkHtfQBMiQJ5AacAYgymyWEZ0DWsxErCZDNES0SlfUPMTRAVW49y5SMBsJGDXYCQgPAFRwSWQwhyXpv0FiQQAebn99jYXKuN6B3KFJt5h5AmUDmLFZHS1pCYjAf2sJ+ANRp84tY4twU4B4QlYNRgJiDRdfMPYyIwfUpjnrUdAELuf+Yr2quifx0kEiFJCLJv3RSIBrY2SBuamg6aKKR0ktosligQ0WB3kSOIJDEzkv0cAmB0dwRiwra0KVZGhcrlOB2lPvomixmzQzfEENBkJzDOGi+NtZk0gAiI1ZNNgJKDXMdhN+nkiMJjnbmFBdWSI3Mo6JypsJnDOYdQzSgcB5AmUEmZjJB0U9QS0JwI2kwHeoDxCuhhWSwrEjT42HeT2h2DUy16OFnEmmCR6sQCNYoA8X8ppMWBHRzUA2SOrtJkoHQSQJ1BKRNNBQe2mg8Q46TF3AEGJF02JqNUk/yzmeAIBSZNRgCDRJNHBKR/MBh0qbfn9uTHG8Ju/uhoPfmhF9LUquynnkYB2f3pEUWIx6uALSppOB4maeZFWKIZuYQCwGiORQMwKRLc/pEk/QFBmNWLaPz8SWJrnRjHBshrHnD9X2U1UIkqUFuZIiahPwyWiomZeDCErlnSQNVF1UFDS5PA4QaJI4NKUD7XO/PoBySARIEqO2XSQtpvFgNnccrEZw7F9Ah5/SJPD4wROizE630kwMu1HTZm5QCeai4tEgCg15OogSdPNYiI9EhWBYosE5pSISpocHidIFAkMT/tR41CHCFTaTZj0BhGUwot/cIZo4h1G1UGlg9moRyA2EtBidVDkyXhQpIOKxRiO/Cx8c4zhkCZHRgjiRcAbkDDtD6HGqQ4RiDaMeXIXDWhCBKg6qHSIegKiOkiDpYcOc3wkoN2bZCxGPYNex+LSQVJU9LRImcWIgDT7+yZ2/KpFBKoivQO5TAlp7x1GFDVibIQ/FAZj2hQBUTI5OCU8geKIBBhjsBnnLpt3B7RdHRS/U2BoWm0ikPuuYe29w4iixmzQRzqG5YUyhSjTy5bZ6iAfTHqdJoUsGZa4PcMev/b7BIDZSaJix69aPAESAaLkkDuGw5GtYtp8whR9AoFQGGVWgyaFLBmxi2U45/AEtW0Ml8WtmByOpINqKRIgiMJgNugQkMLwBiVNVgYBgFGvi45RKJZUkMBmml027w+FIYW5JhfKCJxxk0RHpv1gbPbmW2hE13LJiwBVB5UOYlbQlDekyblBApEndxZJeajAEuMJeDS8UEYQ7wkMz/jhsptg0Kvj1mjQ61BuNZIIUHVQ6SDy5xPegKZz6SJPXiyNYoLYSECMy9ZyOiiRJ1CtEj9A4Mrx/CDtvsuIokT4AFqPBMQQuWJLB1kTRQIaLhF1xnsC037VVAYJquwmjM2QCBAlgnj6n/QGNTk8TiDy5MXSIyCwmGZFILpkXsORgBDrqVgRUFkkUGk3UbMYUTqI0dFT3qAmx0gL7KbijARsxtl0kFgtqeUSUb2OwWE2YNoXBOccwzPqiwQoHUSUFCIdNO0PabZEFJhdvF4sc4ME1iKLBIDZ0RFTvhACobDqRKDKLi+W4Zzn5PokAoSqiDWDtVoiCsRGAtp9Sk5EbJ+Atwg8AUCO1qZ9wdlGMRWKQCjMMeUNLf7BGaCJdxmViJYOsSKg5UhAeALOIksHWU3ygD8pzDW9ZD4Wp8WAKW9Idd3CgmjDWI58AU2IAJWIlg6xU0M1HQmI6qAiM4ZjJ4lGPQGNRwJOiwHT/mC0W1iNkQAAjLn9Obm+dt9lRFEyNx2k3SfMojWGTbOLZUQkYNXwzwkQy+ZDqk0HuSKTREdzVCZKIkCoitgbv7abxYrTGLbERgIBCVajHnqdtmcjCWN4eNoPo56hXGU/s0q7fJ5clYlq911GFCXFEgk0VdlgNuhUM4hMKWL3DLv9Ic1XBgEiEghGewTUNvAvGgnkqEyURIBQFcVSHfThNXXY8/UbUGFTxyAypYhNB3kCkqaHxwmcFgOCEkffuEd1qSBAFl6rUZ+zrmHtvsuIosQ8Jx2k3RsMY6zoBACYjc68AUnzqyUFImXXNeJWpQgAkdERFAkQpUCxRALFijXOEyiGdJDo5VDj8DhBe409Z/u2tS/jRFFh0DHoGBDm2o4EihUxIsITkD0BrTeKAbOTRAH1VQYJ/uPz23N2bXrUIlQFYyx686dIQH2ISMBbRJFAbEOfWkUgl2jiXUYdw6WFGByXq/CXyByLSf7ZeAMhuAMhTQ+PE8yJBFSaDsolmhAB6hguLYQvoOU+gWJF3PS9kY5higS0D73LCNUxmw7S/g2m2IimgwJhuAPkCRQDJAKE6hBegJaXyhQreh2DyaCDOxCCLxguikjAYTJA9IeptTool2hfxomiQ0QCWl4qU8xYjfroHJti6BPQRRbLhMO8KCKbdCm975hQPcILoHSQOrEa9RiNTLQsho5hQB70Z9Sra1xEviARIFSHOZoOokhAjdhMeoxExi4XQyQAyL6As8gWAKVKaX7XhKqZTQcVx1NmsWGJSQdZi8ATAIBPX9lSNIKWLqX5XROqJpoOokhAlVhNenQOFY8nAACf3N5S6CMUDHqXEarDbNDBoGMw6OnXU43YTHoEpLD8/4vEEyhl6F1GqA6zQU+NYiom1rAvlkiglKF3GqE6ml02NLvshT4GkYTYdZLF0CdQ6pCME6rjix9Yhi9c217oYxBJiL3xl2JdfbFBP0FCdeh0DDqUZs22FrBQJFBUaCIdRFNECUI9iBu/jtGQv2JAEz9BmiJKEOpBeAJ2k0F1S9mJ9NGECBAEoR5EgxiVhxYHJAIEQaSFEAEqDy0OSAQIgkgLkQ6iSKA4IBEgCCItoiJAkUBRQCJAEERaRD0BKg8tCkgECIJIi9jqIEL7kAgQBJEWFAkUFyQCBEGkhbj508iI4oBEgCCItLAYKRIoJkgECIJIC1EVRJFAcUAiQBBEWpRbjfjwmjpsb6sq9FEIBSApJwgiLfQ6hsc/s6XQxyAUgiIBgiCIEoZEgCAIooQhESAIgihhSAQIgiBKGBIBgiCIEoZEgCAIooQhESAIgihhSAQIgiBKGMY5L/QZUoYxNgygB0A5gMkMLlENYETRQxHJyPRnpHbU+n0V4ly5/pq5uL4S18z2GoW4f7VwzmsS/YWmREDAGHucc35/Bp+3n3NOrY55INOfkdpR6/dViHPl+mvm4vpKXDPba6jt/qXVdNCuQh+AWJRi/Rmp9fsqxLly/TVzcX0lrpntNVT1O6TJSCBTKBIgCEKrUCSgDI8X+gAEQRAZkpP7V0lFAgRBEMRcSi0SIAiCIGIgESAIgihhSAQIgiBKmJIVAcaYnTH2M8bYvzPGPlno8xAEQaQDY6ydMfZjxthz2VynqESAMfYTxtgQY+x43Os3M8bOMMY6GWMPRV6+E8BznPP7ANye98MSBEHEkc49jHPexTn/fLZfs6hEAMATAG6OfYExpgfwAwC3AFgD4B7G2BoAjQB6Ix8m5fGMBEEQyXgCqd/DFKGoRIBz/iaAsbiXtwHojKhmAMDTAO4A0AdZCIAi+3cgCEKbpHkPU4RSuPk1YPaJH5Bv/g0AngdwF2PsEaisjZsgCCKGhPcwxpiLMfYogI2Msa9lenFDtqfTKpxzN4B7C30OgiCITOCcjwJ4INvrlEIk0A+gKebPjZHXCIIgtEBO72GlIAL7ACxnjLUxxkwA7gbwYoHPRBAEkSo5vYcVlQgwxp4C8C6AlYyxPsbY5znnIQBfBvAqgFMAnuGcnyjkOQmCIBJRiHsYDZAjCIIoYYoqEiAIgiDSg0SAIAiihCERIAiCKGFIBAiCIEoYEgGCIIgShkSAIAiihCERIAiCKGFIBAiCIEoYEgGCIIgS5v8HU/+QHt+TzhMAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAD8CAYAAACRkhiPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAdhElEQVR4nO3deXjU5b028PubnUAgZAUSQvaEiFAUQUX2LZxasZuV2p5WEYsFjlZbxJ7F9z2nvtSqXRAUsXCwWqW81FMR2fdFUALKJgRCQiBhSQgQCJD9e/4gYgiZZJJZnpn53Z/rmuvqPDPzmzs+JXd+yzwjqgoiIrImP9MBiIjIHJYAEZGFsQSIiCyMJUBEZGEsASIiC2MJEBFZWIDpAG0RFRWliYmJpmMQEXmV3bt3n1PV6OYe86oSSExMRE5OjukYREReRUQKbT3Gw0FERBbGEiAisjCWABGRhbEEiIgsjCVARGRhLAEiIguzRAnU1Sv+8XkxKmvqTEchIvIoliiBPScu4Om/fYF3d9q8VJaIyJIsUQJ3JUZgaHo05mzMw6XKGtNxiIg8hiVKAABmjMvAxas1eHPzMdNRiIg8hmVKoE9cF0z4Rg8s2FaAs5cqTcchIvIIlikBAHh2TAbq6hV/XHfUdBQiIo9gqRJIiAzFI4N6YUnOSRwrrTAdh4jIOEuVAABMG5mKkAA/vLI613QUIiLjLFcCUZ2CMXloMlYeOIPPT1wwHYeIyCi3lYCIJIvIAhFZ2tKYOzw+JBlRnYLw0qrDUFV3vjURkUexqwREZKGIlIjIgSbj2SKSKyJ5IjKzpW2oar6qTmptzB06BQdg+sg07Mw/j81HSt399kREHsPePYFFALIbD4iIP4C5AMYDyAIwUUSyROR2EVne5Bbj1NROMHFgAhIiQvHSqlzU13NvgIisya4SUNUtAM43GR4IIK/hr/lqAIsBTFDV/ap6f5NbSXsDisgTIpIjIjmlpc77qz0owA/Pjk3HodOXsGzvKadtl4jImzhyTiAOwMlG94saxpolIpEiMg9AfxF53tZYU6o6X1UHqOqA6Ohmvye53b7Vtwdu69EZr6zJRVUtF5cjIutx24lhVS1T1SmqmqKqs2yNuZOfn+C57EwUXbiG9z494e63JyIyzpESKAbQs9H9+IYxrzIkLQqDUyPx2oY8XObickRkMY6UwC4AaSKSJCJBAB4GsMw5sdxH5PrewPkr1Xhra4HpOEREbmXvJaLvA9gBIENEikRkkqrWApgGYDWAQwCWqOpB10V1nb7x4fjm7d3x5635KL1cZToOEZHb2Ht10ERV7a6qgaoar6oLGsZXqGp6wzH9F10b1bV+OS4D1bX1eG0DF5cjIuuw3LIRtiRFdcQP7uqJ9z49gcKyK6bjEBG5BUugkadGpSHQ3w+vrjliOgoRkVuwBBqJ6RyCSfclYdneUzhQXG46DhGRy7EEmnhiWDK6hgbipVWHTUchInI5lkATnUMCMXVEKrYePYdtR8+ZjkNE5FIsgWb86O5eiAvvgN+t5lLTROTbWALNCAn0xy/GpGNfUTlW7D9jOg4RkcuwBGz4dv84ZMSG4eXVh1FTV286DhGRS7AEbPD3E8zIzsDxsqv4266Trb+AiMgLsQRaMDIzBncldsWf1h/F1epa03GIiJyOJdACEcHM8ZkovVyFhdu4uBwR+R6WQCvu7BWBMVmxeHNzPi5cqTYdh4jIqVgCdpgxLgNXqmsxd2Oe6ShERE7FErBDWmwYvntHPP6yoxBFF66ajkNE5DQsATv9Ykw6IMAf1nKpaSLyHSwBO/UI74Cf3puIDz4vwuEzl0zHISJyCpZAG/x8eAo6BQfg5VW5pqMQETkFS6ANwkODMGVYCtYfLsGu4+dNxyEichhLoI0eG5yEmLBg/HYlF5cjIu/HEmijDkH+eHp0OnYXXsDaL8+ajkNE5BCWQDs8NCAeyVEd8fLqXNTVc2+AiLwXS6AdAvz98MtxGThaUoG/7ykyHYeIqN1YAu00vk839Ivvgj+sPYLKmjrTcYiI2oUl0E4igufGZ+J0eSXe2VFoOg4RUbuwBBxwb0oUhqZHY87GPJRfqzEdh4iozVgCDpoxLgPl12rw5uZjpqMQEbUZS8BBfeK64IF+PbBwewHOXqo0HYeIqE1YAk7w7Nh01NYp/riOi8sRkXdhCThBr8iO+OGgBCzJOYljpRWm4xAR2Y0l4CTTR6YhOMAPr67h4nJE5D1YAk4SHRaMx4ckY8X+M9h78qLpOEREdmEJONHkIUmI6BiEl1ZxcTki8g5uKwERSRaRBSKytNFYbxGZJyJLReRJd2VxlbCQQEwbkYpPjpVh69FzpuMQEbXKrhIQkYUiUiIiB5qMZ4tIrojkicjMlrahqvmqOqnJ2CFVnQLgIQCD2xreEz1ydwLiu3bAS6sOo56LyxGRh7N3T2ARgOzGAyLiD2AugPEAsgBMFJEsEbldRJY3ucXY2rCIPADgYwAr2vUTeJjgAH88MyYdB09dwkf7TpmOQ0TUIrtKQFW3AGj6VVoDAeQ1/IVfDWAxgAmqul9V729yK2lh28tUdTyAR5p7XESeEJEcEckpLS2176cybMI34pDZLQyvrjmC6tp603GIiGxy5JxAHICTje4XNYw1S0QiRWQegP4i8nzD2HARmS0ib8LGnoCqzlfVAao6IDo62oG47uPvJ3guOxMnzl/F4l0nTMchIrIpwF1vpKplAKY0GdsEYJO7MrjT8IxoDEyKwOz1R/HdO+LRMdht/6mJiOzmyJ5AMYCeje7HN4wRri81PXN8Js5VVGPBtgLTcYiImuVICewCkCYiSSISBOBhAMucE8s33JHQFWOzYjF/Sz7KKqpMxyEiuoW9l4i+D2AHgAwRKRKRSapaC2AagNUADgFYoqoHXRfVO83IzsDV6lrM3cilponI89h1oFpVJ9oYXwEfubTTVVJjwvC9O+Px7s5CPHZfIuK7hpqORER0A5eNcIOnR6cDAvx+7RHTUYiIbsIScIMe4R3w6L2J+J/Pi3H4zCXTcYiIbmAJuMmTw1MQFhyA363iUtNE5DlYAm4SHhqEKcNTsOFwCT4raPrhayIiM1gCbvTovUmI7RyM3648xKWmicgjsATcqEOQP54alY49Jy5i7ZdnTcchImIJuNtDA+KRHNURL6/ORW0dF5cjIrNYAm4W4O+HX43LwNGSCnywh6tsEJFZLAEDsvt0Q7+e4fjDuiOorKkzHYeILIwlYICI4LnsDJwur8Rfdhw3HYeILIwlYMi9KVEYmh6NuRuPofxajek4RGRRLAGDZozLQPm1Gry5mYvLEZEZLAGD+sR1wQP9emDh9gKcvVRpOg4RWRBLwLBnx6ajtk4xe/1R01GIyIJYAob1iuyIHw5KwOJdJ1Fw7orpOERkMSwBDzB9ZBqCA/zwyhouLkdE7sUS8ADRYcF4/L4kfLzvNPYVXTQdh4gshCXgISYPTUZExyC8tOqw6ShEZCEsAQ8RFhKIqSNSsT2vDFuPlpqOQ0QWwRLwID+6OwFx4R3w0qrDqK/nUtNE5HosAQ8SHOCPZ8ak40DxJSzff9p0HCKyAJaAh3mwfxwyu4Xh1TW5qK7lUtNE5FosAQ/j7yd4LjsThWVXsXjXCdNxiMjHsQQ80PCMaAxKisDs9UdxparWdBwi8mEsAQ8kInhufCbOVVTjz1sLTMchIh/GEvBQdyR0RfZt3TB/yzGcq6gyHYeIfBRLwIP9KjsDlbX1mLMhz3QUIvJRLAEPlhLdCd+7Ix7vfXoCRReumo5DRD6IJeDhnhqdBgjwp3VcapqInI8l4OF6hHfAj+/uhb/vKUJeyWXTcYjIx7AEvMDPh6egQ6A/Xl1zxHQUIvIxLAEvENkpGI8PScbKA2e41DQROZXbSkBEkkVkgYgsbTQ2XES2isg8ERnurize6PEhSYjoGITfrjwMVS4uR0TOYVcJiMhCESkRkQNNxrNFJFdE8kRkZkvbUNV8VZ3UdBhABYAQAEVtCW41YSGBmD4yFZ8cK8PmI1xqmoicw949gUUAshsPiIg/gLkAxgPIAjBRRLJE5HYRWd7kFmNju1tVdTyA5wD83/b9CNbxyKBe6BUZit+uPIw6LjVNRE5gVwmo6hYA55sMDwSQ1/AXfjWAxQAmqOp+Vb2/ya3Exna/WibzAoDgdv4MlhEU4IcZ4zJx+MxlfLCHO05E5DhHzgnEATjZ6H5Rw1izRCRSROYB6C8izzeMfUdE3gTwDoA5Nl73hIjkiEhOaSkPg/zT7d3Qr2c4Xl1zBJU1dabjEJGXc9uJYVUtU9UpqpqiqrMaxj5Q1Z+p6g9UdZON181X1QGqOiA6OtpdcT2WiODX4zNx5lIlFm7n4nJE5BhHSqAYQM9G9+MbxsjFBiVHYlRmDOZtOobyazWm4xCRF3OkBHYBSBORJBEJAvAwgGXOiUWteWZsOi5V1uLPW/NNRyEiL2bvJaLvA9gBIENEikRkkqrWApgGYDWAQwCWqOpB10Wlxm7r0QXf7NsdC7YVcKlpImo3e68Omqiq3VU1UFXjVXVBw/gKVU1vOM7/omujUlPPjElHZU0d3th0zHQUIvJSXDbCi6VEd8J374jHOzsLcbr8muk4ROSFWAJe7qnRaVBVzF7PpaaJqO1YAl4uvmsoHhnUC0tyipBfWmE6DhF5GZaAD5g6IhXBAX74/VouNU1EbcMS8AHRYcGYdF8Slu87jQPF5abjEJEXYQn4iMlDkxEeGoiXV+eajkJEXoQl4CM6hwTiyWEp2HykFDvzy0zHISIvwRLwIT+5NxGxnYPx/1Yc4lLTRGQXloAPCQn0x6//qTf2FZXjvU8LTcchIi/AEvAxD/TrgftSo/C7VbkouVxpOg4ReTiWgI8REfzXg31QVVeP/1p+yHQcIvJwLAEflBTVET8fnoKP9p7CFn4fMRG1gCXgo6YMS0FSVEf8n2UHUVtX3/oLiMiSWAI+KiTQHzPHZyL/3BX8nd9HTEQ2sAR82NisWPTrGY4/rTvK7yMmomaxBHyYiGDGuAycKq/EXz89YToOEXkgloCPG5wahcGpkXh9Yx4qqmpNxyEiD8MSsIBfjs1A2ZVqLNxWYDoKEXkYloAF9E/oinG3xeKNTceQV8LvHCCir7EELOI/J/RBhyB/TH//c54kJqIbWAIWEds5BK98vy8Onb6EWSv4SWIiuo4lYCEjM2Mx6b4kvL2jEGsOnjEdh4g8AEvAYmZkZ6BPXGfM/GA/rxYiIpaA1QQH+OM3D96O81eq8dedXG6ayOpYAhb0jZ7hGJIWhbe25vMkMZHFsQQsatqIVJyrqMbiz/hJYiIrYwlY1KDkSAxMjMCbW/JRVcu9ASKrYglY2LSRqThdXokP9hSbjkJEhrAELGxIWhT6xXfB65vyUM8vpieyJJaAhYkIHh2chJPnr+HzkxdNxyEiA1gCFjciIwYBfoJ1h86ajkJEBrAELK5LaCAGJkVg3ZcsASIrclsJiEiyiCwQkaWNxoaIyDwR+bOIfOKuLHSzMVmxOFpSgePnrpiOQkRuZlcJiMhCESkRkQNNxrNFJFdE8kRkZkvbUNV8VZ3UZGyrqk4BsBzA220NT84xuncsAPCQEJEF2bsnsAhAduMBEfEHMBfAeABZACaKSJaI3C4iy5vcYlrZ/g8BvNfG7OQkPSNCkdktDGt5SIjIcgLseZKqbhGRxCbDAwHkqWo+AIjIYgATVHUWgPvtDSAiCQDKVfWyva8h5xvdOxZvbD6GC1eq0bVjkOk4ROQmjpwTiANwstH9ooaxZolIpIjMA9BfRJ5v9NAkAP/dwuueEJEcEckpLS11IC61ZHRWLOrqFZuOlJiOQkRu5LYTw6papqpTVDWlYW/hq/EXVNXmSWFVna+qA1R1QHR0tHvCWlDfuC6ICQvmISEii3GkBIoB9Gx0P75hjLyQn59gVO9YbM4tbXFl0UuVNZixdC9OXbzmxnRE5CqOlMAuAGkikiQiQQAeBrDMObHIhAf69cCV6jr8/91FNp+zaPtxLMkpwrK9p9yYjIhcxd5LRN8HsANAhogUicgkVa0FMA3AagCHACxR1YOui0qudndyBO5ICMe8TcdQU1d/y+MVVbVYuL0AAPBZwXl3xyMiF7CrBFR1oqp2V9VAVY1X1QUN4ytUNb3hOP+Lro1KriYimD4qDcUXr+F/Pr/1yN5fdxbi4tUa9I3vgl3Hz6OOi84ReT0uG0E3GZ4ejT5xnfH6xrybfslX1tThra35GJIWhccGJ+FyZS0On7lkMCkROQNLgG4iIpg2Ig3Hy65i+b6vj/sv/uwEzlVUY9qIVAxMigDAQ0JEvsCuD4uRtYzNikV6bCf8af1RVNVePzcwb3M+BiZGYFByJAAgvmsHfFZwHo8OTjIZlYgcxBKgW/j5CZ4alY6p7+3BjKX7boz//qF+N/73wKQIbM4thapCREzEJCInYAlQs77ZtzvuShyF6oarhIID/BEdFnzj8UFJEfhgTzGOlV5BakwnUzGJyEEsAbIppnOIzccGJl0/LPRZwXmWAJEX44lhapfEyFBEhwXjs4Iy01GIyAEsAWoXEcHApAh8WnAeqvy8AJG34uEgardBSRH4eN9p5J+7grjwDi0+NyTQ/5YxVb1x9VFTQf5+8PPjCWciV2MJULt99XmBUa9ubvW5z45Jx/RRaTeN/es/DuC9T080+/xv9AzHP6YOdjykDYkzP8a/jEzFM2MzXPYeRN6AJUDtlhEbhle+3w+ll6tafN7rm/JQeP7qLeOFZdf3IH50d6+bxtd8ecYt33c8e0MeS4AsjyVA7SYi+N6d8a0+792dhTYf694lBE8OT7lp7HT5NX7pPZGb8MQwEZGFsQSIiCyMJUBEZGE8J0BeI6+kAq+uyUVN3a2fSwgLCcBvHuyDjsHO+7/0h18Uo14V3+7f+nkPIm/FEiCvseVIKVYeOIOM2DD4N/oMQUVVLU6cv4of39MLdyR0ddr7vbOjEAdPXcLQtGhEdgpu/QVEXoglQF5nyc/uQZfQwBv3N+WW4Kf/vcsl73Wtpg4LtxfgV+MyXbJ9ItN4ToCoFW9/UojyqzWmYxC5BEuAqAWxnYNRUVWLt3ccNx2FyCVYAkQtSI3phNG9Y7BwewEqqmpNxyFyOpYAUSumjkjFxas1+Gszn3wuv1qDj/aewodfFOPDL4qRc7z937u8v6gcVbV1jkQlajOeGCZqRf+ErrizV1cs33caPxt28xIXC7blY/aGvBv3A/wEX/5nNoIC2vb3VdGFq/jWnG2YOLAnZn2nr1NyE9mDewJEdojoGITa+ls/n1BZW4+gAD+sf3YYJg9JQm29or4d369wseHE896T5Q5nJWoLlgCRg/xFkBLdCREd+VkC8j4sASIiC2MJEBFZGE8ME3m4pt/hLOL8r910x3uQZ2IJEHm4B1//BHtPXrxx/1fjMjB1RKrTtr+78AImvrUT1Q3f99wh0B9Ln7wHt/Xo0uprr1XX4duvb8fkIcn4rh1fMESeh4eDiDxc3tnL6J8QjqdHpyEsJADHSiucuv2iC1dRXVuPn9zTC0+PTkOAv+C19XmtvxDA4l0ncPjMZadnIvfhngCRF7gzoSueHp2OpbuLXPYeP7k3EcnRnVBXr3htQx6OnL2M9Ngwm8+vqq3Dm5vzXZaH3IN7AkR0k0cHJyE0yB9zN7a8N/D33cU4c6nSTanIVdxWAiKSLCILRGRpo7EsEVkiIm+IyPfclYWIbIvoGIQf3d0LH+09hYJzV5p9Tk1dPV7flId+PcMR4MeTyN7MrhIQkYUiUiIiB5qMZ4tIrojkicjMlrahqvmqOqnJ8HgAr6nqkwD+uU3JichlHh+ShEB/P7yxqfm9gWVfnELRhWuYPiIVvJDIu9l7TmARgDkA/vLVgIj4A5gLYAyAIgC7RGQZAH8As5q8/jFVLWlmu+8AeEFEHgAQ2bboROQqMWEhmDgwAe/uLERFVS0EN/+m33X8PHp374xRvWMMJbStpq4ef1x3BA/flYCeEaGm43g8u0pAVbeISGKT4YEA8lQ1HwBEZDGACao6C8D9dm63BMDUhkL5wN7QROR6U4al4POTF3Hk7K1X/nQNDcLz4zM98vMEy744hbkbj6Fb5xD8+J5E03E8niNXB8UBONnofhGAQbaeLCKRAF4E0F9EnlfVWQ3F8msAHQG8bON1TwB4AgASEhIciEsmNbemWkvrrLV9CTZytm5dQvDh1MGmY7RJXb1iro1DWNQ8t10iqqplAKY0GTuOhl/wLbxuPoD5ADBgwAD+biAim1YeOI380uZPZlPzHLk6qBhAz0b34xvGiOzW3NEEzzvAQN6gvl4xZ0MeYsK4mmtbOFICuwCkiUiSiAQBeBjAMufEIiJqm3WHzuLwmct4Ymiy6Shexa7DQSLyPoDhAKJEpAjAC6q6QESmAViN61cELVTVgy5LSpbn7mOBPPboGapq67Dt6DnU1F1f26hzSCDuSYm86aS0qmLOxjz0igzFA/164DcfH2p1u8dKK3D07OUb9/vGh6NHeAfn/wAezt6rgybaGF8BYIVTExG1hseLLGXVgTN4avEXN42te2YYUmM63bhfWHYV+4rK8cK3suBv54fXfv7uHuQ2KoGM2DCsfGoI/Cz24TcuG0HUgqbXx5P7VdbUAQDefmwg/uP+rJvGvlJbf30vIaqT/ecDrtXUYURGNFY+NQT/fn8Wcs9extpDZ52U2nuwBIjs1HTNfdvPc+A92v9Sn5cW0wnxXZ17uCY8NAi9u3fGT+7phV6RoZizIc/uefYVLAEiO9izP+DI56Y88DNXlhLg74epw1Oxv7gcm4+Umo7jViwBIiIAD/aPQ1x4B7xmsb0BlgCRt3Hy7ycL/b5rUVCAH6YMS8buwgvYmX/edBy34ZfKEHkRVx428sR1gNzt+wN64rUNeViwLR/3pNy8puX7n53ACx8ehLbSwuGhQVj7i6EIDw1yZVSnYQkQETUICfRHVo/OKL1cdctjuWeuX046eYjtD6MdLanA2i/P4lxFFUuAiMgkVxzlCgn0w4zsTJuPf7T3FNZ+6V2XmfKcABGRhbEEiIgsjCVARGRhLAEiIg/38urDWLitwCXbZgkQEXm49YdK8GlBmUu2zRIgIrIwlgARkYWxBIiILIwlQERkYSwBIiILYwkQEVkYS4CIyMJYAkREFsYSILdobQ32W57fjiUgPeXLUdr6s970Wjt+CA/5MW/ikZk8MZQHEm/6GjURKQVQ2MxDXQCU2zEWBeCcC6K1prks7tqOPa9x5DltGfekOQGcMy+umhN7ntfS447Mi7fPSXu344x/K66aE8CxeemlqtHNPqKqXn8DMN/OsRxPyeeu7djzGkee05ZxT5oTZ82Lq+bEnue19Lgj8+Ltc+LKeTE1J66cF185HPSRnWOmOCtLe7Zjz2sceU5bxj1pTgDn5HHVnNjzvJYe99Z58fZ/K143J151OMhRIpKjqgNM56CvcU48D+fEM7lqXnxlT8Be800HoFtwTjwP58QzuWReLLUnQEREN7PangARETXCEiAisjCWABGRhVm2BEQkWUQWiMhS01noayLyoIi8JSJ/E5GxpvMQICK9RWSeiCwVkSdN56GviUhHEckRkfvbuw2fKgERWSgiJSJyoMl4tojkikieiMwEAFXNV9VJZpJaSxvn5R+qOhnAFAA/MJHXCto4J4dUdQqAhwAMNpHXKtoyLw2eA7DEkff0qRIAsAhAduMBEfEHMBfAeABZACaKSJb7o1naIrR9Xv6t4XFyjUVow5yIyAMAPgawwr0xLWcR7JwXERkD4EsAJY68oU+VgKpuAXC+yfBAAHkNf/lXA1gMYILbw1lYW+ZFrnsJwEpV3ePurFbR1n8rqrpMVccDeMS9Sa2ljfMyHMDdAH4IYLKItOv3eUD743qNOAAnG90vAjBIRCIBvAigv4g8r6qzjKSzrmbnBcB0AKMBdBGRVFWdZyKcRdn6tzIcwHcABIN7AiY0Oy+qOg0AROSnAM6pan17Nm6FEmiWqpbh+nFn8iCqOhvAbNM56GuqugnAJsMxyAZVXeTI633qcJANxQB6Nrof3zBGZnFePA/nxDO5dF6sUAK7AKSJSJKIBAF4GMAyw5mI8+KJOCeeyaXz4lMlICLvA9gBIENEikRkkqrWApgGYDWAQwCWqOpBkzmthvPieTgnnsnEvHABOSIiC/OpPQEiImoblgARkYWxBIiILIwlQERkYSwBIiILYwkQEVkYS4CIyMJYAkREFsYSICKysP8FnTjMS2gSVpoAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def approxpolygammaone(b):\n",
" assert b >= 1\n",
" \n",
" if b >= 10:\n",
" z = 1 / b\n",
" # Assuming[z > 0, FullSimplify@Series[PolyGamma[1, 1/z], { z, 0, 8 }]]\n",
" return z * (1 + z * (1/2 + z * (1/6 + z**2 * (-1/30 + z**2 * (1/42)))))\n",
" else:\n",
" # MiniMaxApproximation[PolyGamma[1, z], { z, { 1, 10 }, 4, 4 }] \n",
" numerator = 52.654073150313565+b*(176.21984028201697+b*(230.4407865623894+b*(232.30138091080408+b*(0.00006342628796017858))))\n",
" denominator = 1+b*(-7.709336904239235+b*(80.59885744577618+b*(114.25915616534267+b*(232.30352278139097))))\n",
" \n",
" return numerator / denominator\n",
" \n",
"def flass():\n",
" from matplotlib import pyplot as plt \n",
" import numpy as np\n",
" import scipy.special as sc\n",
" \n",
" b = np.geomspace(1, 10, 100)\n",
" exact = sc.polygamma(1, b)\n",
" approx = [ approxpolygammaone(x) for x in b ]\n",
" plt.plot(b, np.abs(exact - approx))\n",
" plt.xscale('log')\n",
" plt.yscale('log')\n",
" plt.show()\n",
" \n",
" plt.figure()\n",
" \n",
" b = np.geomspace(10, 10000, 100)\n",
" exact = sc.polygamma(1, b)\n",
" approx = [ approxpolygammaone(x) for x in b ]\n",
" plt.plot(b, np.abs(exact - approx))\n",
" plt.xscale('log')\n",
" plt.yscale('log')\n",
" plt.show()\n",
" \n",
"flass()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "db6902b5",
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"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.8.13"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment