Skip to content

Instantly share code, notes, and snippets.

@timvieira
Last active May 29, 2023 13:46
Show Gist options
  • Star 3 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save timvieira/656d9c74ac5f82f596921aa20ecb6cc8 to your computer and use it in GitHub Desktop.
Save timvieira/656d9c74ac5f82f596921aa20ecb6cc8 to your computer and use it in GitHub Desktop.
Don't use scipy.special.digamma, if you care about speed.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Inefficiency of scipy.special.digamma\n",
"\n",
"In this notebook, we look at the inefficiency the digamma function in scipy ([`scipy.special.digamma`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.digamma.html)).\n",
"\n",
"**UPDATE:** A fix will appear in scipy - see [issue \\#8127](https://github.com/scipy/scipy/issues/8127).\n",
"\n",
"**Warning:** All of this discussion is going to be specific to the version of the relevant libraries as well as the Python version. Here is what I'm using:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Python version: 2.7.14 |Anaconda custom (64-bit)| (default, Oct 27 2017, 18:21:12) \n",
"[GCC 7.2.0]\n",
"scipy version: 1.0.0\n",
"numpy version: 1.14.1\n"
]
}
],
"source": [
"import scipy\n",
"import scipy.special as sp\n",
"import pylab as pl\n",
"import numpy as np\n",
"import sys\n",
"\n",
"print 'Python version:', sys.version\n",
"print 'scipy version:', scipy.__version__\n",
"print 'numpy version:', np.__version__"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XuUHOV55/Hv05cZ3TWSRoDQSBoJ\nyYAw95Gw4wAysEbYJ4jEKCs768VZO9gGNtllL4bsHvssCdkl2Q17souC8cIJcdYRt6yRExyCzc03\nQCMQEEFkRhfQSFiMNJKQQBpNdz/7R781arW61T1Sz3RX6/c5R0fdb1e9VdU91U+/9bzvW+buiIiI\n1Eqi3jsgIiLNRYFFRERqSoFFRERqSoFFRERqSoFFRERqSoFFRERqSoFFRERqSoFFRERqSoFFRERq\nKlXvHaiH9vZ27+zsrPduiIjEytq1a3e6+/RKy52UgaWzs5Pu7u5674aISKyY2dvVLKdLYSIiUlMK\nLCIiUlMKLCIiUlMKLCIiUlMKLCIiUlNVBRYzW2pmG8ysx8xuK/F6q5k9FF5/0cw6C167PZRvMLOr\nK9VpZnNDHW+FOltC+VfN7HUzW2dmPzGzhZW2ISIio69iYDGzJHAPcA2wEPhc4Zd68CVgt7vPB+4G\n7grrLgRWAOcAS4GVZpasUOddwN3uvgDYHeoG+K67n+vuFwB/DPzpsbYx7HdCRERqopoWy2Kgx903\nufshYBWwrGiZZcCD4fGjwJVmZqF8lbsPuPtmoCfUV7LOsM4VoQ5CndcBuPv7BdsbD0T3VC63DZEh\new8M8vi6bfXeDZGTQjWBZSawteB5bygruYy7Z4C9wLRjrFuufBqwJ9Rx1LbM7GYz20i+xfK7w9g/\nOcn9v5d7+b1V63h374F674pI06smsFiJMq9ymVqV5x+43+PuZwBfB/7zMPYPM7vRzLrNrLuvr6/E\nKtLMduwbAKD/g0N13hOR5ldNYOkFZhU87wC2l1vGzFLAZKD/GOuWK98JtIU6ym0L8pfOrhvG/uHu\n97l7l7t3TZ9ecaobaTI7Q2DZ8+FgnfdEpPlVE1jWAAtCb60W8ony1UXLrAZuCI+vB552dw/lK0Kv\nsbnAAuClcnWGdZ4JdRDqfBzAzBYUbO8zwFsF2y61DZEhO/fnA8vuD9ViERlpFSehdPeMmd0CPAkk\ngQfcfb2Z3QF0u/tq4H7gO2bWQ76lsiKsu97MHgbeADLAze6eBShVZ9jk14FVZvaHwCuhboBbzOwq\nYJB8b7EbKm1DJLJzfz6gqMUiMvIs30g4uXR1dblmNz65fPy//oh39x7k33/qI9xyxYLKK4jIUcxs\nrbt3VVpOI++l6bk7u0KLZbdaLCIjToFFmt77BzIcyuYAXQoTGQ0KLNL0+kLiHmCPkvciI06BRZpe\n1CMsnTT1ChMZBQos0vT6whiWzmnj2XNAl8JERpoCizS9qMWy4NQJyrGIjAIFFml6O/cPkExYvsXy\n4SFyuZOvi73IaFJgkaa3c98hpo5vYer4FnIO+wYylVcSkeOmwCJNb+f+AaZPaKVtXAugnmEiI02B\nRZrezv0DtE9sZcq4NKBBkiIjTYFFmt7O/Ydon9CiFovIKFFgkabm7vQNXQrLt1jUM0xkZCmwSFPb\nN5DhUCZH+4RWpoQWiwZJiowsBRZpatHgyPaJLUweqxaLyGhQYJGmFt05sn1CK8mEMWlMSjkWkRGm\nwCJNLbrBV/uEVgCmjG9RrzCREabAIk0tms5l+sR8YGkb16L5wkRGmAKLNLWd+wdIGEOJ+ynj0roU\nJjLCFFikqe3cP8DU8fn8CkDb2LR6hYmMMAUWaWp9+/KDIyNt41rUK0xkhCmwSFPbuX9gKL8C+Uti\n+w5myIRbFYtI7aXqvQMix+O/P7mB//OTTUPPF3VO5S//1WLM7Ijldu4fYG77+KHnQ6PvDwwO9RST\n49e9pZ//+OhrPPa1X2HK+JbKK8hJQS0WiaXXtu1l0pg0N3y8k8sWTOfHb+1k/fb3j1jG3enbN1B0\nKUyDJGvpf/zDL9i08wO27TlQ712RBqLAIrGUyeaYPXUct3/6bP7k+vNpSSV4pHvrEcvsH8gwEKZz\niUzRRJQ18/I7u/n5pl0AZHTzNCmgwCKxlMk6qWT+stfkcWmuPuc0Hn91OwOZ7NAyxYMj4XCLRYMk\nT9zKZzYOPVbOSgpVFVjMbKmZbTCzHjO7rcTrrWb2UHj9RTPrLHjt9lC+wcyurlSnmc0NdbwV6mwJ\n5bea2Rtm9pqZ/cjM5hSskzWzdeHf6uN7KyROBnM5UonDf77LL+5gz4eD/PCN94bKigdHglostbLh\nl/v44Zs7uGTuVAAGs2qxyGEVA4uZJYF7gGuAhcDnzGxh0WJfAna7+3zgbuCusO5CYAVwDrAUWGlm\nyQp13gXc7e4LgN2hboBXgC53Pw94FPjjgu0fcPcLwr9rh/UOSCwVtlgAPjG/nRmTx/DI2sOXwwrn\nCYsox1Ibf/5sD+NaknzpV+cCkMmpxSKHVdNiWQz0uPsmdz8ErAKWFS2zDHgwPH4UuNLy3XOWAavc\nfcDdNwM9ob6SdYZ1rgh1EOq8DsDdn3H3D0P5C0DH8A9XmkUm50e0WJIJ47MXdfD8L/r45d6DAPTu\nzieU2yceTt5PaE2RSpgGSZ6Arf0f8v3X3uW3LplNe2gNKscihaoJLDOBwqxobygruYy7Z4C9wLRj\nrFuufBqwJ9RRbluQb8X8oOD5GDPrNrMXzOy6UgdhZjeGZbr7+vrKHavERCabI508smvx9Rd3kHNY\n+WwPv/vXr/BHP3iT2VPHMW384RaLmdE2Lq0cywn4qxfexoAvXzqPdAjuGV0KkwLVjGOxEmXFf0Xl\nlilXXiqgHWv5wxsy+xdAF3B5QfFsd99uZvOAp83sdXffeEQl7vcB9wF0dXXpLIi5TM5JJY/8M+ps\nH8/izqn85c/fZmw6ydcuP4OvXHbG0HQukbZxLew9oBbL8chkc/zNK9v45FmncOqkMfR/cGioXCRS\nTWDpBWYVPO8AtpdZptfMUsBkoL/CuqXKdwJtZpYKrZYjtmVmVwH/Cbjc3QeicnffHv7fZGbPAhcC\nRwQWaS6D2RypxNG/Q37/M2fzozd38IWPz+GUiWNKrts2Ns3uD9RiOR7P/aKPvn0DLL84fyU6ajUO\n6lKYFKjmUtgaYEHordVCPhlf3PNqNXBDeHw98LS7eyhfEXqNzQUWAC+VqzOs80yog1Dn4wBmdiHw\nLeBadx/q+mNmU8ysNTxuBz4BvDGcN0HiJ5P1koHlgllt/LtPnVk2qEC+xaIcy/F5pLuX9gktfPKs\nUwBIDl0KU4tFDqvYYnH3jJndAjwJJIEH3H29md0BdLv7auB+4Dtm1kO+pbIirLvezB4m/0WfAW52\n9yxAqTrDJr8OrDKzPyTfE+z+UP4nwATgkTBtxzuhB9jZwLfMLEc+UP43d1dgaXKlLoVVa8q4NOu3\nq8UyXLv2D/DDN3fwxV/pJB3e+yi4K3kvhaqaK8zdnwCeKCr7RsHjg8DyMuveCdxZTZ2hfBP5XmPF\n5VeVqf9nwLnHPgJpNpnc0cn7auWT92qxDNf31m0nk3OWdx2+ih0FGCXvpZBG3kss5S+FHd+fb9u4\nFg4O5jg4mK28sAD5edce6d7K+R2TOfO0iUPl0VgijWORQprdWGJpMJs7YoDkcBwefT/IaZOTtdyt\nptK3b4BnN7yHA7s/OMQ//XIff3DdR49YJupurJH3UkiBRWIpP0Dy+ALLjLZ8Yn/d1j0snXxaLXer\nqfz+/3udp97YMfR8YmuKa887/YhlklGLRcl7KaDAIrHj7mRPIHl/6fx2Zk8dx58/t5Grzzn1qHu4\nSH4usKfe2MFXLp/HFz6Wn5Zv0tg0k8akj1hOyXspRTkWiZ3oSyx9nC2WVDLBVy6fx6tb9/Dzjbtq\nuWtN497nNjKuJclXLzuDjinj6Jgy7qigAkreS2kKLBI70ZfY8bZYAD57UQfTJ7Zyz7M9tdqtprG1\n/0NWv7qdzy+eXfGukMmEYabkvRxJgUViZzB8iR1vd2OAMekkv3PpXH7as4tXt+6p1a41hW89v5Gk\nGV++dF5Vy6cTCSXv5QgKLBI7UYuleA6w4fr8JXOYPDbNSrVahry37yAPd/fy2Ytnctrk8rMXFEol\nTcl7OYKS9xI70ZfYiVwKg/wU+jf8Sid/9qO3eGvHPhacOrHySk3of/zDBh5ak59s/OBglkw2x1cu\nO6Pq9ZMJU/JejqAWi8TOiSbvC/32r3QyNp3kz587OecsfXfvAe59biMzp4zlyrNP4TPnzeDOXz+X\nzvbxVdeRTiaUY5EjqMUisVOL5H1kyvgWPn/JbP7iZ1v4t1d9hFlTx51wnXHy7ec3k3P4sxUXHvex\npxKmXmFyBLVYJHZqkbwv9OVL55Iw+PaPN9Wkvrjo/+AQf/3SOyy74PQTCqjppJL3ciQFFomdWiXv\nIzMmj+U3LuzgoTVb6ds3UHmFJvEXP93MgcEsX7u8+nxKKamk6VKYHEGBRWJnMEreH+cklKV8dckZ\nDGZzPPDTzTWrs5HtH8jw4M/f5upzTj3hTgtJXQqTIgosEjvZKHlfo0thAHPbx3PNuTP4zs/fZu+B\n5r9Xy3dfzB/nTUvmn3Bd6YSS93IkBRaJnehLrBbJ+0I3LTmD/QMZ/uqFt2tab6M5OJjl2z/ezK/O\nb+f8WW0nXF9+HItaLHKYAovETpQorkV340LnnD6ZJWdO54GfbObAoea9V8tjL/fSt2+Am5acWG4l\nkkomdM97OYICi8ROrZP3hW7+5Hx2fXCIh9a8U/O6G0Emm+Pe5zZy/qw2Pn7GtJrUmU5o5L0cSYFF\nYmdwhC6FASzqnMqizinc9/wmDmWa78vy715/l639B7h5yRk1u12AkvdSTIFFYiebrX3yvtBNS+az\nfe9BHl+3bUTqr5dczln5zEYWnDKBq84+tWb1auS9FFNgkdgZSt7XsLtxoSVnTufsGZO497mN5Joo\nd/D0P73Hhh37uOmTZ5Co4WXE/DiW5nmf5MRpSheJncERbrGYGTctOYN//dev8NSbO7j6nHjevtjd\n+a3/8yKvb9sLwMBgjo4pY/m1otsLn6iUps2XIgosEjsj1d240KfPncGdf/cmD63ZGtvAsvbt3fxs\n4y6uPudUTm8bC8A1H51R8/ctrWnzpYgCi8RO9Os4NQK9wiLJhPEbF83k3uc2suP9g5w6qbp7kzSS\nh7u3Mq4lyZ/+5gWMbx25U13T5ksx5VgkdqKR96kRuhQWuf7iDnIOf/Ny/JL4Hx7K8Hevvctnzp0x\nokEFlLyXo1UVWMxsqZltMLMeM7utxOutZvZQeP1FM+sseO32UL7BzK6uVKeZzQ11vBXqbAnlt5rZ\nG2b2mpn9yMzmFKxzQ1j+LTO74fjeComLzAjMFVbKvOkT6JozhUfWbsU9Xr/In3j9l3xwKMvyrlkj\nvi1Nmy/FKp6ZZpYE7gGuARYCnzOzhUWLfQnY7e7zgbuBu8K6C4EVwDnAUmClmSUr1HkXcLe7LwB2\nh7oBXgG63P084FHgj8M2pgLfBC4BFgPfNLMpw30jJD5GOnlf6De7ZrGp7wNefmfPiG+rlh7p3krn\ntHEs6hz5UyGlafOlSDU/+RYDPe6+yd0PAauAZUXLLAMeDI8fBa60/OirZcAqdx9w981AT6ivZJ1h\nnStCHYQ6rwNw92fc/cNQ/gLQER5fDTzl7v3uvht4inwQkyY1Gsn7yKfPm8HYdJJH124d8W3Vyju7\nPuTFzf0s75pVs0GQx5LWtPlSpJozcyZQeFb1hrKSy7h7BtgLTDvGuuXKpwF7Qh3ltgX5VswPhrF/\nmNmNZtZtZt19fX0lD1TiYTSS95EJrSk+fe4Mvv/qu7GZP+zRtVtJGPzGRaVOndpLJRK6FCZHqCar\nV+rsLf4rKrdMufJSAe1Yyx/ekNm/ALqAy4exf7j7fcB9AF1dXToLYmwoeT8KgQVgeVcHj73cyz+7\n+znGpJMY8O8+dSZLP9p43ZAz2RyPvbyNX10wnRmTx47KNnWjLylWTYulFyjMAHYA28stY2YpYDLQ\nf4x1y5XvBNpCHUdty8yuAv4TcK27R7f6q2b/pIlEyfuRmISylEvmTuXGy+ZxfkcbZ546kZ37B3i4\nuzEvjf3d6++ybc8BfuuS2aO2TSXvpVg1LZY1wAIzmwtsI5+M/3zRMquBG4CfA9cDT7u7m9lq4Ltm\n9qfA6cAC4CXyrYyj6gzrPBPqWBXqfBzAzC4EvgUsdff3Crb9JPBHBQn7TwG3D+M9kJgZzDnppI1K\n/gDyI/F//9NnDz2/7bHXeOL1d8nlvKZTo5yowrnA/lkN5wKrJJVMkMk57j5qn4k0tootlpDvuIX8\nF/ibwMPuvt7M7jCza8Ni9wPTzKwHuBW4Lay7HngYeAP4e+Bmd8+WqzPU9XXg1lDXtFA3wJ8AE4BH\nzGxdCFq4ez/wB+QD4BrgjlAmTSqTzY14V+NjWdQ5lfcPZvjFe/vqtg+lRHOBfW1JbecCqyS6L44G\nSUqkqpFT7v4E8ERR2TcKHh8ElpdZ907gzmrqDOWbyPcaKy6/6hj79wDwQPkjkGYymPVRy6+Usqhz\nKgBrNvdz1mmT6rYfhdydlc/2MLNtLL92fm3nAqsk6p2XyTrp5KhuWhqURt5L7GRzPuKj7o9l1tSx\nnDqplZe27K7bPhR7cXM/L7+zh69ePo/0KHTDLpQaarEogS95CiwSO5lcblTGsJRjZizqnMqazf0N\nMyL/nmd6aJ/QMioj7YtFQV4JfIkosEjsDGa95ve7H65FnVP55fsH6d19oK77AfB6715+/NZO/tWv\nzmVMHa5FRUF+UC0WCRRYJHYy2fq2WKAgz7Kl/v1EVj7bw8QxKb7wsTmVFx4BQ8l7tVgkUGCR2Bms\nc44F4MzTJjJxTIo1dc6z9Ly3n79f/0tu+HgnE8ek67IPhcl7EVBgkRjKdzeub2BJJoyuOVPq3mK5\n97mNtKYS/PYnOuu2D9FnoUthElFgkdjJ5ryu41giXZ1T6XlvP/0fHKrL9rftOcD3XtnGikWzmTah\ntS77AIeT91mNY5Gg/menyDANZn1UpsyvZPHcfJ6lu06tlm8/vwmA37lsXl22H4mC/KBuTyyBAovE\nTr27G0fO65hMSyrBzzbuGvVt7/1wkFVr3uHXL5zJzLbRmWyynLS6G0uR+p+dIsNU75H3kdZUkqvO\nPoXHXu5l38HBUd32C5t3cXAwx28uGv1xK8WGkvfKsUigwCKxk+9uXP/AAvC1y+ez72CGv3rhnVHd\n7prN/bSkEpzXMXlUt1tK1N1Yd5GUiAKLxE6jJO8Bzu2YzKUL2rn/J5s5ODh6NwJbs6WfCzraaE3V\nf3Ku6PYFSt5LpDHOTpFhaJTkfeSmJfPZuX+AR0bpHi0fHsrwj9vfZ9Hckb+ffTWGRt4reS+BAovE\nTiZX32nzi31s3lQumt3Gt57fNCpfrq+8s4dszukKo//rTcl7KdY4Z6dIlTLZ+o+8L2Rm3LRkPr27\nD/D9V0f+5qUvbe7HDC6e0yAtloSS93IkBRaJncFc/UfeF7virFM489SJrHx2I7kRzjV0v93P2adN\nYlKdpnApFrVYlLyXiAKLxE426w0xjqVQImHc9Mkz6HlvP0+9uWPEtjOYzfHy23uGBmc2AiXvpVhj\nnZ0iVYjued9oPnPuDGZPHcfKZzeO2H1a1m9/nwODWbo6G+MyGDB0YzEl7yWiwCKxU+973peTSib4\nyuXzeHXrnhEbjb9mc376mMUNkriHght9qcUiQeOdnSIVNFryvtBnL+pg+sRWVj7bMyL1r9nSz5xp\n4zhl0pgRqf94DCXv1WKRQIFFYmcwlxv1+7pXa0w6ye9cOpef9uxi3dY9Na3b3el+ezddcxqntQJK\n3svRGvPsFDmGbM6HEsaN6POXzGHy2DT3PruxpvW+FaboX9wgAyMjSt5LMQUWiRV3b4h73h/LhNYU\nKxbN4odv7qBv30DN6v3eK9tIGCw585Sa1VkLad3zXooosEisRL+KG627cbHlXR1kcs73XtlWk/qy\nOedvXt7G5R+ZzqkNlF+Bw3eQ1Mh7iTT22SlSJDMUWBq3xQIw/5SJXDCrjUfWbq1J1+Mfv9XHL98/\nyPKu+k+TXyw5FFjUYpG8qgKLmS01sw1m1mNmt5V4vdXMHgqvv2hmnQWv3R7KN5jZ1ZXqNLO5oY63\nQp0tofwyM3vZzDJmdn3R9rNmti78Wz38t0HiIhorkW7A7sbFlnd18Isd+3mtd+8J1/XI2l7axqW5\n8uzGugwG+Slt0kljUDkWCSqenWaWBO4BrgEWAp8zs4VFi30J2O3u84G7gbvCuguBFcA5wFJgpZkl\nK9R5F3C3uy8Adoe6Ad4Bvgh8t8RuHnD3C8K/a6s6coml6FJYIyfvI792/um0phI8svbEZj3e8+Eh\nnlq/g+sumNkQ0+SXkkoklLyXIdX87FsM9Lj7Jnc/BKwClhUtswx4MDx+FLjSzCyUr3L3AXffDPSE\n+krWGda5ItRBqPM6AHff4u6vAWpvn8SiLq2NOPK+2KQxaa756GmsXrf9hO7VsvrV7RzK5rj+4o4a\n7l1tpRKmkfcypJrAMhMo/MnVG8pKLuPuGWAvMO0Y65YrnwbsCXWU21YpY8ys28xeMLPrSi1gZjeG\nZbr7+vqqqFIaUTSDbqMn7yPLu2bx/sEM//DG8c8f9kh3LwtnTOKjM+t/t8hyUklT8l6GVHN2lvpp\nWPwXVG6ZWpVXMtvdu4DPA//TzM44qhL3+9y9y927pk+fXkWV0oiiL69Gm924nI/Pm8bMtrE8urb3\nuNZ/Y/v7vL5tL8u7Gre1AvlAr2nzJVJNYOkFCruidADFN50YWsbMUsBkoP8Y65Yr3wm0hTrKbeso\n7r49/L8JeBa4sPJhSRwNJe9j0mJJJIylHz2NFzbtYiAz/Mth3/7xJsa1JPn1C6tpuNdPOmEaeS9D\nqjk71wALQm+tFvLJ+OKeV6uBG8Lj64GnPd/HcjWwIvQamwssAF4qV2dY55lQB6HOx4+1c2Y2xcxa\nw+N24BPAG1Ucl8RQnJL3kUWdUzmUyfH6MHuHbe3/kNWvbufzi2fTNq5lhPauNlJJJe/lsIqBJeQ7\nbgGeBN4EHnb39WZ2h5lFPbDuB6aZWQ9wK3BbWHc98DD5L/q/B25292y5OkNdXwduDXVNC3VjZovM\nrBdYDnzLzKLlzwa6zexV8kHpv7m7AkuTilPyPrIoTHH/0pb+Ya33rec3kjD48qXzRmK3akrJeymU\nqrwIuPsTwBNFZd8oeHyQ/Bd+qXXvBO6sps5Qvol8r7Hi8jXkL40Vl/8MOLfiQUhTGErex2AcS2Ta\nhFbOmD6e7i27q17nvX0Hebi7l+sv7uC0yY010r4UJe+lUHzOThEOt1gafeR9scVzp9K9pb/q2xbf\n/5PNZLI5vnLZUf1QGlIqoeS9HKbAIrGSiVnyPtI1ZyrvH8ywYce+isvuPTDI/33hHT5z3ul0to8f\nhb07cemkkvdyWLzOTjnpDc0VFqPkPTB0j/o1VeRZ/va17ewfyPCVyxo/txJRd2MppMAisRKXSSiL\ndUwZy2mTxrCmijzLS5v7mT6xlXNOnzQKe1YbyYRyLHKYAovESnQpLE7Je8hP1NjVOYU1m/srznbc\nvWU3izunkp/hKB7SSdM972VIvM5OOenFNXkP+cthv3z/IL27D5RdZtueA2zbc2Coi3JcpBIJTZsv\nQxRYJFai6/hxS95DfqAk5C91lbMmvLZobmPd174SJe+lUPzOTjmpxW2usEIfOXUiE8ek6H77GIFl\nSz8TW1OcdVp88iug7sZyJAUWiZXDvcLi96ebTBhdc6Ycu8WypZ+L5kyJ1ZQ1AEnlWKRA/M5OOakN\nJe9jmGMB+OjMyWze+QGHMkf/ut/9wSF+sWP/UNfkOEmrV5gUUGCRWBmMaXfjyNz28eQc3un/8KjX\n1r6d74oc5WLiJJVU8l4OU2CRWMnE6J73pUQj6bfs/OCo19Zs6aclmeC8jsa9oVc5uue9FIrn2Skn\nrUyMuxsDzJ0WAsuuowPLS1v6Oa9jMmPSjXlf+2NRd2MppMAisRLn5D3AlPEttI1Ls6moxXLgUJbX\ne/fGrptxJKXkvRSI59kpJ624J+8BOqeNP+pS2Lqte8jknMUxzK9Avvu3kvcSUWCRWBmM6SSUhea1\nHx1Y1oaxLRfNjteI+4gmoZRCCiwSK5lsjlTCYjWPVrHO9vFs33uQA4eyQ2Xrtu7ljOnjmTwuXcc9\nO37RPe8rzYMmJwcFFomVTM5jfRkMDvcMe7s/32pxd9Zt3cP5s9rquVsnJBWm2NF97wUUWCRmMlmP\nbVfjyLyiLsfv7j3Izv0DXBDrwJIP9krgCyiwSMxkcjmSTdJiiXqGvbp1DwDnd8Q4sCQUWOQwBRaJ\nlcGsx7arcWRCa4r2Ca1DLZZ1vXtoSSY4a8bEOu/Z8Ys+E41lEVBgkZjJZHOkY95igahnWH5al1e3\n7uHs0yfRmorfwMhI9Jlo6nwBBRaJmWZI3gN0to9j084PyOac13v3cn4Mp3EpFCXv1eVYQIFFYiaT\ni3/yHvJ5lp37B1i3dQ8fHMrGOr8CBTkWtVgEBRaJmUw2F7t7lZQS9QxbvW4bQKy7GoN6hcmRqgos\nZrbUzDaYWY+Z3Vbi9VYzeyi8/qKZdRa8dnso32BmV1eq08zmhjreCnW2hPLLzOxlM8uY2fVF278h\nLP+Wmd0w/LdB4mIw60OXXeIs6hn2/dfeZWJraijQxJWS91Ko4hlqZkngHuAaYCHwOTNbWLTYl4Dd\n7j4fuBu4K6y7EFgBnAMsBVaaWbJCnXcBd7v7AmB3qBvgHeCLwHeL9m8q8E3gEmAx8E0zi+e8GFJR\nJtccyfvOMMtx/weHOG/WZBIxb4UpeS+FqvnptxjocfdN7n4IWAUsK1pmGfBgePwocKXl59xYBqxy\n9wF33wz0hPpK1hnWuSLUQajzOgB33+LurwHFP4muBp5y93533w08RT6ISRPKZD3W84RFxqSTnD55\nDBDv8SuRoRaLkvdCdYFlJrC14HlvKCu5jLtngL3AtGOsW658GrAn1FFuW8ezf5jZjWbWbWbdfX19\nFaqURpXJ5ZriUhgcvhwW9/yv8tUDAAAM80lEQVQKKMciR6rmDC3187D4r6fcMrUqP5aq1nH3+9y9\ny927pk+fXqFKaVTN0mKB/G2KgVhP5RJJR92NdSlMgFQVy/QCswqedwDbyyzTa2YpYDLQX2HdUuU7\ngTYzS4VWS6ltldq/JUV1PVthHYmpwZwzrklaLMu7ZtE2Ls2pk8bUe1dOWHKou7EuhUl1LZY1wILQ\nW6uFfDJ+ddEyq4GoN9b1wNOenz97NbAi9BqbCywAXipXZ1jnmVAHoc7HK+zfk8CnzGxKSNp/KpRJ\nE8pkc6SbpMVywaw2/sPVZ9V7N2piKHmvS2FCFYEltBxuIf9l/SbwsLuvN7M7zOzasNj9wDQz6wFu\nBW4L664HHgbeAP4euNnds+XqDHV9Hbg11DUt1I2ZLTKzXmA58C0zWx+20Q/8AflgtQa4I5RJE8pk\nm2PkfbNRd2MpVM2lMNz9CeCJorJvFDw+SP4Lv9S6dwJ3VlNnKN9EvtdYcfka8pe5Sm3jAeCBYx6E\nNIXBJkreN5OUuhtLAZ2hEivZnDfNpbBmktaNvqSAAovESibrJJtgrrBmM5S81zgWQYFFYmawSabN\nbzbRxKC6FCagwCIx0yzT5jeboQGSSt4LCiwSM4PZXOzvINmMUupuLAV0hkqsZHOuS2ENKLoUllWL\nRVBgkZhR8r4xJTVXmBTQGSqxMtgk0+Y3GyXvpZACi8RGNue4oxxLA1LyXgrpDJXYGAxfWuoV1nii\nGaeVvBdQYJEYiUZ161JY4zEzUgkjqwGSggKLxEh0rw8l7xtTKmm6H4sACiwSI4Ph17BaLI0plUgo\neS+AAovESPRrWMn7xpRKmuYKE0CBRWJEyfvGphaLRBRYJDaUvG9s6aSS95KnwCKxEV1m0aWwxqTk\nvUR0hkpsDA7lWNRiaUSpRELjWARQYJEYGUre69bEDSmVMI28F0CBRWIk6m6s5H1jSiWVvJc8BRaJ\njaHkvXIsDUnJe4noDJXYUHfjxpZKmKbNF0CBRWIko+R9Q8uPY1GLRRRYJEaGuhsred+Q1N1YIjpD\nJTbU3bixpZLqbix5VQUWM1tqZhvMrMfMbivxequZPRRef9HMOgteuz2UbzCzqyvVaWZzQx1vhTpb\njrUNM+s0swNmti78u/d43wxpbIdH3uv3UCNKa9p8CSqeoWaWBO4BrgEWAp8zs4VFi30J2O3u84G7\ngbvCuguBFcA5wFJgpZklK9R5F3C3uy8Adoe6y24j2OjuF4R/Xx3WOyCxoeR9Y9OlMIlU89NvMdDj\n7pvc/RCwClhWtMwy4MHw+FHgSjOzUL7K3QfcfTPQE+orWWdY54pQB6HO6ypsQ04SSt43NiXvJVJN\nYJkJbC143hvKSi7j7hlgLzDtGOuWK58G7Al1FG+r3DYA5prZK2b2nJldWuogzOxGM+s2s+6+vr4q\nDlsajZL3jS0/bb5aLFJdYCn187D4r6fcMrUqP9Y23gVmu/uFwK3Ad81s0lELut/n7l3u3jV9+vQS\nVUmji5L3abVYGlIqkdClMAGqCyy9wKyC5x3A9nLLmFkKmAz0H2PdcuU7gbZQR/G2Sm4jXGbbBeDu\na4GNwEeqOC6JmUxWLZZGlk6aLoUJUF1gWQMsCL21Wsgn41cXLbMauCE8vh542t09lK8IPbrmAguA\nl8rVGdZ5JtRBqPPxY23DzKaHzgCY2bywjU3VvwUSF9FlFiXvG1MqaUM99+Tklqq0gLtnzOwW4Ekg\nCTzg7uvN7A6g291XA/cD3zGzHvItlRVh3fVm9jDwBpABbnb3LECpOsMmvw6sMrM/BF4JdVNuG8Bl\nwB1mlgGywFfdvf/43xJpVBnNFdbQlLyXSMXAAuDuTwBPFJV9o+DxQWB5mXXvBO6sps5Qvol8r7Hi\n8pLbcPfHgMcqHoTEXnQpLKkcS0PSXGES0U8/iY2h5L0uhTWkVFLJe8lTYJHYyORyJBOGhi81pnTS\nhu6ZIyc3BRaJjUzONTiygaUSCdwhp8thJz0FFomNTNY1T1gDi3rrqdUiOkslNjLZnBL3DSxqTSrP\nIgosEhuDOVfivoFFA1cVWESBRWIjk82R0hiWhpXWpTAJdJZKbGRyrlH3DSwK+hp9LwosEhtK3je2\noeS9Rt+f9HSWSmxE41ikMSl5LxEFFomNwazGsTSyoeS9ciwnPQUWiY1MNqdLYQ0suk/OoFosJz2d\npRIbSt43tqjFouS9KLBIbGSyrinzG5iS9xKpatp8yRvIZHnmn/rqvRsnrb79A5wysbXeuyFlREH/\nZxt3seP9gTrvjRQ6f9ZkZkweO2rbU2AZhv0HM3z1r9bWezdOameeNrHeuyBlTB3fAsCfPLmhznsi\nxc46bSI/+L1LR21mcAWWYZg8Ns0Tv3tpvXfjpDa3fXy9d0HKWHj6JJ7/D59k/0Cm3rsiBX7as5M7\nn3iTZza8xxVnnToq21RgGYZUMsHC0yfVezdEGtbsaePqvQtSZMGpE/iLn23hnmc28skzTxmVVosy\noSIiTSydTHDjZfNY+/ZuXtrcPyrbVGAREWly/3zRLNontLDy2Y2jsj0FFhGRJjcmneS3PzGX537R\nxz9u2zvi21NgERE5CXzh43OY2Jriz0eh1aLkvYjISWDSmDRfXXIGBwezuPuIJvEVWEREThI3f3L+\nqGxHl8JERKSmqgosZrbUzDaYWY+Z3Vbi9VYzeyi8/qKZdRa8dnso32BmV1eq08zmhjreCnW2HO82\nRERk9FUMLGaWBO4BrgEWAp8zs4VFi30J2O3u84G7gbvCuguBFcA5wFJgpZklK9R5F3C3uy8Adoe6\nh72N4b4RIiJSG9W0WBYDPe6+yd0PAauAZUXLLAMeDI8fBa60fGZoGbDK3QfcfTPQE+orWWdY54pQ\nB6HO645zGyIiUgfVBJaZwNaC572hrOQy7p4B9gLTjrFuufJpwJ5QR/G2hruNI5jZjWbWbWbdfX2a\noVhEZKRUE1hK9UkrvpNPuWVqVX482ziywP0+d+9y967p06eXWEVERGqhmsDSC8wqeN4BbC+3jJml\ngMlA/zHWLVe+E2gLdRRva7jbEBGROqgmsKwBFoTeWi3kE+Wri5ZZDdwQHl8PPO3uHspXhB5dc4EF\nwEvl6gzrPBPqINT5+HFuQ0RE6qDiAEl3z5jZLcCTQBJ4wN3Xm9kdQLe7rwbuB75jZj3kWxErwrrr\nzexh4A0gA9zs7lmAUnWGTX4dWGVmfwi8EurmeLZRztq1a3ea2dtVvD/ltJNvXcVdsxwH6FgaVbMc\nS7McB5zYscypZiHL/+iX4TCzbnfvqvd+nKhmOQ7QsTSqZjmWZjkOGJ1j0ch7ERGpKQUWERGpKQWW\n43NfvXegRprlOEDH0qia5Via5ThgFI5FORYREakptVhERKSmFFjKOJEZnRtNFcfyRTPrM7N14d+X\n67GflZjZA2b2npn9Y5nXzcz+LBzna2Z20WjvY7WqOJYlZra34DP5xmjvYzXMbJaZPWNmb5rZejP7\nvRLLxOJzqfJY4vK5jDGzl8zs1XAs/6XEMiP3Hebu+lf0j/zYmo3APKAFeBVYWLTMTcC94fEK4KF6\n7/cJHMsXgf9d732t4lguAy4C/rHM658GfkB+mp+PAS/We59P4FiWAH9b7/2s4jhmABeFxxOBX5T4\n+4rF51LlscTlczFgQnicBl4EPla0zIh9h6nFUtqJzOjcaKo5llhw9+fJD44tZxnwl573AvnpgWaM\nzt4NTxXHEgvu/q67vxwe7wPe5OhJYGPxuVR5LLEQ3uv94Wk6/CtOqI/Yd5gCS2knMqNzo6lq9mfg\ns+EyxaNmNqvE63FQ7bHGxcfDpYwfmNk59d6ZSsKllAvJ/zouFLvP5RjHAjH5XCx/76t1wHvAU+5e\n9nOp9XeYAktpJzKjc6OpZj+/D3S6+3nADzn8KyZu4vKZVONlYI67nw/8L+B7dd6fYzKzCcBjwL9x\n9/eLXy6xSsN+LhWOJTafi7tn3f0C8hPzLjazjxYtMmKfiwJLaScyo3OjqXgs7r7L3QfC028DF4/S\nvtVa08x07e7vR5cy3P0JIG1m7XXerZLMLE3+i/j/uvvflFgkNp9LpWOJ0+cScfc9wLPk77BbaMS+\nwxRYSjuRGZ0bTcVjKbrefS35a8txtBr4l6EX0seAve7+br136niY2WnR9W4zW0z+XN1V3706WtjH\n+4E33f1PyywWi8+lmmOJ0ecy3czawuOxwFXAPxUtNmLfYRVnNz4Z+QnM6NxoqjyW3zWza8nPDt1P\nvpdYwzGzvybfK6fdzHqBb5JPSuLu9wJPkO+B1AN8CPx2ffa0siqO5Xrga2aWAQ4AKxr0h8sngC8A\nr4fr+QC/D8yG2H0u1RxLXD6XGcCDZpYkH/wedve/Ha3vMI28FxGRmtKlMBERqSkFFhERqSkFFhER\nqSkFFhERqSkFFhERqSkFFhERqSkFFhERqSkFFhERqan/Dx55GBNEZXkAAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f12d01588d0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"from time import time\n",
"\n",
"def timefn(f, x):\n",
" \"Runtime of `f(x)`.\"\n",
" b4 = time()\n",
" f(x)\n",
" return time() - b4\n",
"\n",
"def bestofk_time(f, x, K):\n",
" \"Runtime as the best of `K` evaluations of `f(x)`.\"\n",
" return min(timefn(f, x) for _ in xrange(K))\n",
"\n",
"\n",
"# Use best-of-100 runtime for comparing \"idealized\" runtimes\n",
"xs = np.linspace(0.001, 3, 100)\n",
"R = 100\n",
"pl.plot(xs, [bestofk_time(sp.digamma, x, R) for x in xs]);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Wow! What happened between 1 and 2?"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A faster implementation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here is an infamous fast implementation of digamma due to [Tom Minka](https://tminka.github.io/) from his [`lightspeed`](https://github.com/tminka/lightspeed/) package."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"def py_fast_digamma(x):\n",
" \"Faster digamma function assumes x > 0.\"\n",
" r = 0\n",
" while (x<=5):\n",
" r -= 1/x\n",
" x += 1\n",
" f = 1/(x*x)\n",
" t = f*(-1/12.0 + f*(1/120.0 + f*(-1/252.0 + f*(1/240.0 + f*(-1/132.0\n",
" + f*(691/32760.0 + f*(-1/12.0 + f*3617/8160.0)))))))\n",
" return r + np.log(x) - 0.5/x + t"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's go ahead Cythonize it, since that'll give us a nice speedup for very little effort."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%load_ext Cython"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%%cython\n",
"cimport cython\n",
"from libc.math cimport log\n",
"\n",
"@cython.cdivision(True)\n",
"cpdef double fast_digamma(double x):\n",
" \"Faster digamma function assumes x > 0.\"\n",
" cdef double r = 0\n",
" while (x<=5):\n",
" r -= 1./x\n",
" x += 1\n",
" cdef double f = 1./(x*x)\n",
" cdef double t = f*(-1/12.0 + f*(1/120.0 \n",
" + f*(-1/252.0 + f*(1/240.0 + f*(-1/132.0\n",
" + f*(691/32760.0 + f*(-1/12.0 + f*3617/8160.0)))))))\n",
" cdef double y = r + log(x) - 0.5/x + t\n",
" return y"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"First, test that we got the implemention correct"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Test that the implementations match on the domain `xs`\n",
"assert np.allclose(sp.digamma(xs), map(fast_digamma, xs))\n",
"assert np.allclose(sp.digamma(xs), map(py_fast_digamma, xs))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, check that the implementations are faster!"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZYAAAD8CAYAAABU4IIeAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xl8VNXd+PHPmZkEwpZAAAUCCasQ\nSMISwqJFwQUUFESriFXcHuvv0VqXVvFxobX6qN3sY6ttte4LoFYFEUtFobWyiyyyBwgkhB0S1kBm\n5vv7Y+4MYZgkk2QmMzd8368XL2bucs65M5n5zrnfc881IoJSSikVKY5YN0AppVTDooFFKaVURGlg\nUUopFVEaWJRSSkWUBhallFIRpYFFKaVURGlgUUopFVEaWJRSSkWUBhallFIR5Yp1A2KhdevWkpGR\nEetmKKWUrXz77bf7RKRNddudlYElIyODZcuWxboZSillK8aYbeFsp6fClFJKRZQGFqWUUhGlgUUp\npVREaWBRSikVURpYlFJKRZQGFqWUUhGlgUUppVREnZXXsSj7W1lYwpfrdgeeZ6WlcGnmOTFs0dlj\n35ETvLtoOx6vF4AWSQncdn5nHA4T45apeKGBRdnSC19u4sv1ezAGRKBxgoPlj19Kk0T9k462P365\niTcXbgu89gAXdG9Nz3NbxLZhKm7oqTBlSyc9Xvp3SmHrM6OZdudgysq9fLluT6yb1eAdLivnw2+L\nuLpfB7Y+M5q/3ZwLQLlbYtwyFU80sChbcnsEl8P35zswoxVtmzdi1qriGLeq4fto+Q6OnvQwaWgG\nAE6n7/SX2zotphRoYFE25fZ6cVrn9J0OwxVZ7Zi3YS9HTrhj3LKGS0R4c2EBOWnJ9O2YAoDL4Q8s\n2mNRp4QVWIwxo4wxG4wx+caYySHWNzLGTLfWLzbGZFRY94i1fIMxZmR1ZRpjOltlbLLKTLSWP2CM\nWWuMWWWM+dIYk15hH48xZoX1b2btXgplJ26v4HKeShZfmdOOk24vc9furmIvVRff5O9ny96jgd4K\nEAjubo8GFnVKtYHFGOMEXgQuBzKBG4wxmUGb3Q4cFJFuwPPAc9a+mcAEoDcwCnjJGOOspszngOdF\npDtw0Cob4DsgV0SygQ+BX1eo/7iI9LX+XVWjV0DZkscrgV/LAP06tqRdcmM9HRZFbywoILVpIqOz\n2wWWJTh9XyEe7bGoCsLpseQB+SKyRUROAtOAsUHbjAXetB5/CFxsjDHW8mkickJEtgL5Vnkhy7T2\nGWGVgVXmOAARmScix6zli4C0mh+uaijKPYLTcerP1+EwjM5qx7827qX0eHkMW9YwFR44xpfrd3ND\nXicauZyB5f4eS7nmWFQF4QSWDkBhhedF1rKQ24iIGygFUqvYt7LlqUCJVUZldYGvF/N5heeNjTHL\njDGLjDHjwjgmZXMer/e0HgvAmJz2lHuEf67ZFaNWNVzvLyvEABMHdTptuf898OipMFVBOIP+Q131\nFPxXVNk2lS0PFdCq2v5URcb8CMgFLqywuJOIFBtjugBfGWNWi8jmoP3uBO4E6NTp9A+Hsp/gHAtA\nTloyHVsl8dbCbRwq8/026ZCSxKg+58aiiQ2GiPDpymKGdm1N+5Sk09b5R+Zp8l5VFE5gKQI6Vnie\nBgSfyPZvU2SMcQHJwIFq9g21fB+QYoxxWb2W0+oyxlwCPApcKCIn/MtFpNj6f4sxZj7QDzgtsIjI\ny8DLALm5ufopsLngHAuAMYbrBnTkd19sZPWO0sDyf94/jB7nNK/vJjYYa4oPUbD/GHdd2PWMdf7g\nrjkWVVE4p8KWAt2t0VqJ+JLxwSOvZgKTrMfXAl+JiFjLJ1ijxjoD3YEllZVp7TPPKgOrzBkAxph+\nwF+Bq0QkcCWcMaalMaaR9bg1cD6wtiYvgrIfd1COxe8nF3dn1S8uY+WUy/j6oeEkuhy8tbCg3tvX\nkHy6qhiXwzCy95k9v8CoMM2xqAqqDSxWz+EeYA6wDnhfRNYYY540xvhHYL0KpBpj8oEHgMnWvmuA\n9/F90f8DuFtEPJWVaZX1MPCAVVaqVTbAb4BmwAdBw4p7AcuMMSvxBaVnRUQDSwPn9npJcIaem6pF\n4wSSkxLo2KoJV+W056PlOzhUpgn92hARPlu1k/O7taZl08Qz1if4T4VpjkVVENbESiIyG5gdtOyJ\nCo/LgB9Wsu/TwNPhlGkt34Jv1Fjw8ksqKX8BkFX1EaiGxuOVwK/lqtwyNIMPvy3iw2VF3HZB53po\nWcOysqiUooPH+enF3UOud+qpMBWCXnmvbKncc2aOJZQ+HZLp3ymFtxYW4NUvvxqbtbKYRKeDy0Kc\nBoNTo8J0uLGqSAOLsiVfjyW8P99JQzMo2H+Mf2/aG+VWNSxer/DZ6p0M69Ga5KSEkNv4e43aY1EV\naWBRtlRVjiXY5X3a0aZ5I95cUBDdRjUAZeUe1hYfYm3xIWat3snO0jLGZLevdHvNsahQ9OYVypbC\nzbEAJLocTMzrxAtfbaJg31EyWjeNcuvsa/LfV/HJilNXEzROcHBJFTdQ09mNVSgaWJTtiEjYORa/\nGwd14sV5+by9aBuPjwme6k4B7D5UxqxVO7kqpz1XZPnmA+vUqgnNGlX+NaGzG6tQNLAo2/F/h7mc\n4Z/JbduiMZdnteP9ZYU8eFkPvdNkCO8u3o5HhAcv60F6ani9Op3SRYWiORZlO/7TLuGeCvO7ZWg6\nh8vcfPzdjmg0y9ZOur28t3g7w89rG3ZQgYoXSGpgUadoYFG2408U1+RUGED/Ti3p3b4Fby3Yhoh+\nEVb0+fc72XfkxGn3WgmHMQanw2iORZ1GA4uyHf+v45r2WIwxTBqawYbdh1m05UA0mmZbbywooHPr\npvygW+sa7+sLLBqo1SkaWJTt+K+ZSKhBjsXvqpz2tGySoPOHVbCqqITvtpdw85B0HDUM1gAJDqM5\nFnUaDSzKdtye2uVYABonOLl+YCf+uXY3xSXHI900W3pzwTaaJjq5dkDt7p2nPRYVTAOLsh3/l1hN\ncyx+PxrcCRHh3cXbItksW9p/5ASfripmfP80mjcOfXV9dVxOh+ZY1Gk0sCjb8Z8Kq8lw44rSWjbh\nkl7nMHVJIWXlnkg2zXamLS3kpNvLpKHptS7D5TA6pYs6jQYWZTt17bGAb/6wA0dP8tmqnZFqlu24\nPV7eXbSN87ul0q1t7W+E5nIYyjXHoirQwKJspy45Fr+hXVPp1rYZby4sOGuHHs9dt5vi0jImDcmo\nUzlOp/ZY1Ok0sCjbcQdGhdU+sBhjmDQknVVFpawoLIlU02zlzQXb6JCSxMW9Kp8LLBwJDocm79Vp\nNLAo2/EErmOp25/v+P5pNG/kOitnPd6w6zALt+znpiHpder5ga/n6NHkvapAA4uynXLrVFhdciwA\nTRu5uGZAGp+t3snewyci0TTbeGthAY1cDq7P7VjnspyaY1FBNLAo2/HU8sr7UG4ekk65R5i6ZHud\ny7KL0uPlfLR8h+9i0RD3sa8pl+ZYVBANLMp2AqPC6pBj8evSphnDerTh3cXbAj2hhu7Db4s4Xu6p\n8bxglXFpjkUF0cCibCdwHUsdcyx+twxNZ/ehE8xZsysi5cUzr1d4e2EBuekt6dMhOSJluhwmMFJP\nKdDAomyoPALDjSu6sEdbOrVqwlsLGv6V+P/auJeC/ce4OUK9FdApXdSZNLAo2/FEYLhxRU6H4eYh\n6SwpOMDa4kMRKTNevbmwgLbNGzGq97kRKzPB6dAcizqNBhZlO7WdNr8qPxzQkaQEZ4Meerx131Hm\nb9jLxEGdSHRF7qPv1FNhKojen1XZzqkbfUXuyzG5SQLj+nXgo+VFPHJFT1Ka1H20VDxYsHkfm3Yf\nAeDrTXtJcBomDuoU0TpceipMBdHAomzHP5NuJEaFVXTT4HSmLtnOzJXF3FzHaU7iweGycm59fSkn\n3Kd6Ez8ckEbb5o0jWo8ON1bBNLAo2/FEYBLKUDLbt6DHOc2YtXJngwgsc9ft5oTby5u35ZFljQBr\n2aR2U+NXRYcbq2CaY1G2E40ci9+Y7PYs3XaAXaVlES+7vn22aiftkxvzg26tadU0kVZNEzEm8q+Z\n5lhUMA0synaikWPxG53dDhGYvdre0+mXHi/nXxv3ckVWu1rdbrgmNMeigoX1yTTGjDLGbDDG5Btj\nJodY38gYM91av9gYk1Fh3SPW8g3GmJHVlWmM6WyVsckqM9Fa/oAxZq0xZpUx5ktjTHqFfSZZ228y\nxkyq3Uuh7MITpRwLQNc2zejVrgWzVhVHvOz69M81uyj3CGNy2ke9Ls2xqGDVBhZjjBN4EbgcyARu\nMMZkBm12O3BQRLoBzwPPWftmAhOA3sAo4CVjjLOaMp8DnheR7sBBq2yA74BcEckGPgR+bdXRCpgC\nDALygCnGmJY1fSGUfUTiRl9VGZPdjuXbS9hRcjwq5deHWat20rFVEjlpkbm6vipOh0MnoVSnCafH\nkgfki8gWETkJTAPGBm0zFnjTevwhcLHxncwdC0wTkRMishXIt8oLWaa1zwirDKwyxwGIyDwROWYt\nXwSkWY9HAl+IyAEROQh8gS+IqQbKfyosGjkWgCuzfb/yP7Npr+Xg0ZN8k7+P0Vnto5JTCebSafNV\nkHACSwegsMLzImtZyG1ExA2UAqlV7FvZ8lSgxCqjsrrA14v5vAbtUw3IqRt9RSdF2Cm1Cdlpycyy\n6W2L/7FmF26vMCa7Xb3U53JqjkWdLpzhxqF+8gT/FVW2TWXLQ30jVLX9qYqM+RGQC1xYg/ZhjLkT\nuBOgU6fIXiCm6pf/13G0eiwAo7Pa8czn63ln0TaaNTr9Y2IMDOveJiJTzkfDZ6t2kpHahN7tW9RL\nfb4eiwYWdUo4gaUIqHg3oDQg+ByBf5siY4wLSAYOVLNvqOX7gBRjjMvqtZxWlzHmEuBR4EIR8d+Z\nqQi4KKis+cEHISIvAy8D5Obm6qfAxso90c2xAIzJac/v/rmRxz75PuT6SUPS+eXYPlGrv7Z2lBxn\nweZ93D28W72cBgNfjsWtORZVQTiBZSnQ3RjTGdiBLxk/MWibmcAkYCFwLfCViIgxZibwnjHm90B7\noDuwBF8v44wyrX3mWWVMs8qcAWCM6Qf8FRglInsq1D0H+N8KCfvLgEdq8Boom/F4BafDRPWLs0NK\nEgseGcHhMvcZ6372wUq+KyyJWt118e4i3wzN1w+s+50hw5XgNIHZEJSCMAKLiLiNMffg+wJ3Aq+J\nyBpjzJPAMhGZCbwKvG2MycfXU5lg7bvGGPM+sBZwA3eLiAcgVJlWlQ8D04wxT+EbCfaqtfw3QDPg\nA+sLZbuIXCUiB4wxv8IXAAGeFJEDdXhNVJxzW4El2lo3a0TrZo3OWD4woxWv/mcLZeUeGic4o96O\ncJWVe5i2tJBLM88hrWWTeqvX6TB4xXevl2hfM6PsIawpXURkNjA7aNkTFR6XAT+sZN+ngafDKdNa\nvgXfqLHg5ZdU0b7XgNcqPwLVkHi83qieBqtO347JlHuEtTsP0b9T/Ixsn7VqJweOnmRSPU9H438v\n3F4hUQOLQq+8VzZU7qmfHktlcjqmALAyjk6HiQhvLiige9tmDOmaWq91O60ZEDSBr/w0sCjb8Xgl\nakONw9EuOYlzWjSKq8DyXWEJq3eUcvPQjHpL2vv5b7imeRblp4FF2U595ViqkpOWwsqi0pi2oaI3\nFxTQvJGL8f3q/xIu/3uhI8OUnwYWZTtuT2xzLOA7HbZ131FKjp2MaTsA9hwuY/bqnVybm0bTRvV/\nJ4yKORalQAOLsiGPV6IyAWVN9PXnWeKg1zJ1cSHlHonZPWRcTs2xqNNpYFG24/ZKVKbMr4mstGSM\niX0C/6Tby7uLt3FhjzZ0bt00Jm0InArTHIuyaGBRtuOJgxxLi8YJdG3TLOaBZc6aXew5fIJbhmbE\nrA0uzbGoIBpYlO2Ux0GOBfwJ/BJEYveF+uaCAtJTm3BhjzYxa4NTcywqiAYWZTvxkGMB6NsphX1H\nTlJ0MDb3bfl+RynLth3kpsHpMb3iPUFzLCqIBhZlO77hxrH/0+2b5k/gx+Z02FsLC0hKcPLD3Pqb\nFywUf4+lXO97ryyx/3QqVUPuGE/p4nfeuc1JdDn4dtvBeq+79Hg5M1YUc3X/DiQnJdR7/RX53wvt\nsSg/DSzKdtweiYvAkuhyMKx7az5dWcwJt6de616+7SAn3N7A3S5jyT/cWHMsyk8Di7KdeMmxANw8\nJIN9R04ye3X93m1yRWEJDgPZ9XBP++qcGhWmp8KUjwYWZTvxkmMBuKBba7q0bsqbC7bVa70ri0ro\n3rZ5TK60D+bUU2EqSHx8OpWqAbfXS0IcnAoDcDgMNw9JZ0VhSb1d0yIirCwsIadj7HsrUHESSg0s\nykcDi7Idd4ynzQ92zYA0miY6eXNhQb3UV3jgOAePlQem7481nTZfBdPAomwnnnIsAM0bJ3DNgDRm\nrdzJviMnol7fd4W+UWh94ySwuHS4sQqigUXZTjzlWPxuHpLOSY+XaUu2R72ulYWlNE5w0OOc5lGv\nKxyaY1HB4uvTqVQY4inH4tetbXMu6Naadxdvj/roqJVFJfRpnxzTm51VpDkWFSw+/jKVqgFPnOVY\n/G4eks7O0jK+WLs7anWUe7x8v6M0bvIrcCrHorMbKz8NLMp2yuMsx+J3ca9z6JCSFNUk/oZdhznh\n9sZNfgV0dmN1Jg0synY8cXA/llCcDsNNQ9JZtOUA63cdikodK6whzXEVWJyaY1Gni79Pp1LVcHu8\ncXkqDOD63I40cjl4a2F0LphcWVhCq6aJpLVMikr5tRGYhFIDi7JoYFG24+uxxGdgadk0kbF92/Px\n8h2UHi+PePkri0rISUvGmPg5fn/v0aPDjZUl9vNBKFVDvhxL/P4munlIBu8vK+KDZYXc8YMuESv3\ncFk5m/YcYXTWqYkny8vLKSoqoqysLGL11JRXhFeuakdy0hHWrVsXs3aoyGncuDFpaWkkJNRu5mwN\nLMp24rnHAtCnQzID0lsybWkht1/QOWK9i/9s2ocI9E8/lV8pKiqiefPmZGRkxKwX4/EKnuJS2iU3\npk3zxjFpg4ocEWH//v0UFRXRuXPnWpURvz/7lApBROLinvfVGde3Pfl7jrBh9+GIlfnWwm10SEli\naNfWgWVlZWWkpqbG9NSYv+oY3qFZRZAxhtTU1Dr1gjWwKFvxX4QXzz0WgFF92uEwMGtlZKbT37Dr\nMAu37OdHg9PPCKqxzrf4a9e40nDU9W9KA4uyFf+Q1njOsQC0ad6IIV1T+Wz1TiQCP+XfWlhAI5eD\nCQNjexviUPxfQlUd5rJly7j33nvrqUW1N3LkSA4frrqXmZaWRklJ1TNZz507l3HjxgHw8ccf85vf\n/CZibbSDsD6dxphRxpgNxph8Y8zkEOsbGWOmW+sXG2MyKqx7xFq+wRgzsroyjTGdrTI2WWUmWsuH\nGWOWG2Pcxphrg+r3GGNWWP9m1vxlUHZhlx4LwJjs9mzdd5Q1xXW7pqX0eDkfLd/BVTntadk0MUKt\niyxjDFJFnyU3N5cXXnihHltUO3PmzKF588jOwXb11Vfz85//PKJlxrtqA4sxxgm8CFwOZAI3GGMy\ngza7HTgoIt2A54HnrH0zgQlAb2AU8JIxxllNmc8Bz4tId+CgVTbAduAW4L0QzTwuIn2tf1eFdeTK\nljzW1d3xnmMBGNX7XJwOw6xVdTsd9uG3RRwv9zBpaEZkGhZBR48eZfTo0Vx76fkMH5LL9OnTWbp0\nKUOHDiUnJ4e8vDwOHz7M/PnzGTNmDAC/+MUvuOmmmxgxYgTdu3fnlVdeAeCmm25ixowZgbJvvPFG\nZs48/Xfi4cOHufzyy8nJyaFPnz58+OGHgK8XMXnyZPLy8hg0aBBbtmwBYPfu3YwfP57c3Fzy8vJY\ntGhRoJxJkyaRlZVFdnY2n3zySaAcf2/kyiuvZMCAAfTu3Zu//e1v1b4Wn332Geeddx4XXHDBacfx\nt7/9jfvuuw+ATZs2MWjQIPLy8nj88cdJSfENxDh06BAjRoygf//+ZGdnM2vWLADy8/Pp06cPt912\nG7179+bmm29mzpw5DB06lB49erBs2TIAHnvsMW655RYuu+wyMjIy+OSTT3jwwQfp06cPo0ePxu12\nAzBlyhQGDhxInz59uOuuuyLSmw5JRKr8BwwB5lR4/gjwSNA2c4Ah1mMXsA/fqdfTtvVvV1mZ1j77\nAFeouq1lbwDXBi07Ut1xVPw3YMAAUfa093CZpD88S95asDXWTQnLTa8ulvOf/VK8Xm+t9vd4vHLh\nr7+S8S99E3L92rVr69K8Ovvwww/ljjvukO+LSmTHwWNSUlIinTt3liVLloiISGlpqZSXl8u8efNk\n9OjRIiIyZcoUyc7OlmPHjsnevXslLS1NduzYIfPnz5exY8eKiEhJSYlkZGRIeXn5afVNmzZN7rrr\nrsDzkpISERHp0KGDPPvssyIi8uqrrwbKue6662ThwoUiIrJ161bp3bu3iIg88MAD8uCDD4qIiNfr\nlQMHDgTKOXjwoIiI7N+/X0REjh49Kr169Qq5jd/Ro0elQ4cOkp+fL16vV8aPHx9owyuvvCI//elP\nRURk5MiR8v7774uIyB//+EdJTk4WEZGTJ0/KoUOHRERk9+7d0q1bNxER2bRpk7hcLlmzZo14PB7J\nycmRO+64I/DaX3PNNSIi8uijj8qwYcOkvLxcli1bJklJSfLPf/5TRETGjBkjn3766WnH5PV6ZcKE\nCTJ79uxK39tQf1vAMgnjOzac4cYdgMIKz4uAQZVtIyJuY0wpkGotXxS0bwfrcagyU4ESEXGH2L4q\njY0xywA38KyIfBLGPsqG/DmWeJs2vzJjstvx0IerWFlUWqtpWL7ZvI+C/ce4/9Ie1W77y0/XsLaO\np92CZbZvwZQre1e6Pisri5/97Gd4E6cwevQYunc8h3bt2jFw4EAAWrRoEXK/sWPHkpSURFJSEsOH\nD2fJkiWMGzeOu+++mz179vDRRx9xzTXX4HKd/hWVnZ3N5MmTmTx5MldeeSXnn39+YN0NN9wA+Ho6\nkyf7zq7PnTuXDRs2BLY5ePAgx48fZ+7cuYFeijGGli1bntHG559/PtBjKioqYvPmzeTm5oY8nrVr\n19KjRw+6du0aaMNbb711xnaLFy9m9uzZAEycOJHHHnsM8P3Af/jhh/nPf/6Dw+GgsLCQffv2AdCt\nWzcyM30ndDIzM7nkkksCr/0zzzwTKPuKK67A5XKRlZUFwKWXXhrYrqCgAIAvv/yS3/zmN5SVlbFv\n3z4GDBjA5ZdfHvKY6iKcwBLqnENw/6mybSpbHupboartq9NJRIqNMV2Ar4wxq0Vk82kNNOZO4E6A\nTp06hVGkikd2yrEAjMw8l0edq5m1srhWgWXmimKaN3Ixsve5UWhd3fXo0YNvv/2WV979O88++QRj\nrhgV1oii4G38z2+66Sbeffddpk2bxmuvvXbGfr169WLZsmXMnj2bn//854wZM4b/+Z//CVkm+L6w\nlyxZQmJi4hnLq2rn3Llz+fe//82iRYtISkriggsuqHb4bV1GUr311luUlpayfPlyXC4XaWlpgfoa\nNWoU2M7hcASeOxyOwCmuits5HI7Tjte/3bFjx7jnnntYvnw5HTp04LHHHovahbXhBJYioOJQlDSg\nuJJtiowxLiAZOFDNvqGW7wNSjDEuq9cSqq4ziEix9f8WY8x8oB+wOWibl4GXAXJzc3VkpE3573US\nj7Mbh5LcJIELurVm7rrdPDYmODVZtZNuL3PW7OLSzHNonOCsdvuqehbRUlxcTKtWrbjq2utJTWnB\nR1PfpLi4mKVLlzJw4EAOHz5MUtKZ85rNmDGDRx55hKNHjzJ//nyeffZZAG655Rby8vI499xz6d37\nzOPZsWMHrVu35qabbiIpKYlp06YF1k2fPp2f/exnTJ06NdCTueSSS3jxxRe5//77AVixYgV9+/bl\nsssu409/+hO//e1vERFKSkpO67WUlpbSqlUrkpKSWLNmDUuXLq3ydcjMzGTjxo1s3bqVjIwMpk6d\nGnK7vLw8Pv74Y6655prT2l5aWkrbtm1xuVx88cUX7Nixo8r6auP48eM4HA5at27N4cOH+fvf/86N\nN94Y8XogvFFhS4Hu1mitRHzJ+OCRVzOBSdbja4GvrPNxM4EJ1qixzkB3YEllZVr7zLPKwCpzBlUw\nxrQ0xjSyHrcGzgfWhnFcyobcXvsk7/3yOqdSsP8YB4+erNF+/8nfy6EyN2Ny2kWpZXW3evVq8vLy\nGH/JBfzx97/mySefZPr06fzkJz8hJyeHSy+9NOSv4ry8PEaPHs3gwYN5/PHHad/eN03NOeecQ69e\nvbj11lsD2xYWFnLVVb4xOStXrmTgwIH07duXX//614HeCsCxY8fIy8vjz3/+M7/73e8AePHFF/nm\nm2/Izs4mMzMzMFBgypQp7N69mz59+tC3b1++/vrr09o3evRojh07Rk5ODk8++SSDBgWf/fcZOXIk\ne/bsoUmTJvzlL3/h8ssv5wc/+AFduoSeyueFF17gueeeIy8vjz179pCcnAz4emoLFiwgNzeXDz74\ngO7du4f1+tdEamoqkyZNok+fPlx99dWVHlNEhJOIAa4ANuLrBTxqLXsSuMp63Bj4AMjHFzi6VNj3\nUWu/DcDlVZVpLe9ilZFvldnIWj4QXw/oKLAfWGMtHwqsBlZa/99e3fFo8t6+Nuw6JOkPz5JZK4tj\n3ZSwfZO/V9IfniVfrd9do/3un/adZP9ijpwo91S6TayT937rdpbKtv1Hw9p2ypQp8pvf/CbkuqNH\nj0qXLl0CSflwhUqox6MjR44EBnK8/fbbMn78+Bi3qHLRTt4jIrOB2UHLnqjwuAz4YSX7Pg08HU6Z\n1vItQF6I5UvxnRoLXr4AyKr2IFSDUG6dCrNTjyU7LQVjfFPeDz+vbVj7lJV7+Ofa3VyRdS6Jrvgf\nqGAwdR66OnfuXG677TYeeOCBwC/5hmbp0qXcd999eL1eWrZsyeuvvx7rJkWFTkKpbMU/KizBJjkW\ngGaNXHRv24yVhVVfrV3Rvzbu5cgJN6Oz21e/cRyoSd76F7/4Rcjll1xyCdu3b69V/UVFRbXar75d\ndNFFrFixItbNiLr4/ymkVAV3LXt6AAAgAElEQVR2zLEA5KSlsLKoNOxf9bNW7aRlkwSGdk2Ncssi\nw6CTUKpTNLAoWwnMFWaT61j8+nZK4cDRkxQdPF7ttsdPevhy3W5G9WlHQpzPiRZgdBJKdYpN/mqV\n8im32XBjv5w03zUs34VxOuxfG/dy7KSHMdnxOxosWCRyLKrh0MCibMVjswsk/c47tzmNXI6w8izf\nbjtAosvBwIxW9dCyyLDXu6GiTQOLshW75lgSnA6yOiSHFVhWFpbSu30LW4wGA9+1GSMvGMD9d91e\n/cYVlJSU8NJLL4W17S9+8Qt++9vfAvDEE08wd+7cGrczlq699trAxJhHjhzhxz/+MV27dqV3794M\nGzaMxYsXV7pvQUEB7713au7dN954g3vuuSci7frTn/4UlZFp9vjLVcri9vhHhdnvTzenYwrfF5cG\nTueF4vZ4Wb2jNHDqzA5eeuklXn3v7/zuz6/WaL+aBJaKnnzyycB8WXawZs0aPB5P4KLJO+64g1at\nWrFp0ybWrFnDG2+8EZgXLJTgwBJJt912W1RuZ2C/T6c6q3m89ruOxS+nYwpl5V42VnG74k17jnC8\n3FOrecVi4a677mLLli3cefP1vPaXP7JkyRKGDh1Kv379GDp0aGACyDVr1pCXl0ffvn3Jzs5m06ZN\nTJ48mc2bN9O3b9+Q9yt5+umnOe+887jkkktOm0jylltuCUyXP3v2bHr27MkFF1zAvffeG5iav7J2\nvPHGG4wbN44rr7ySzp0786c//Ynf//739OvXj8GDB3PgwAHANyz4/vvvZ9iwYfTq1YulS5cyfvx4\nunfvHpg4EmDcuHGBqfVffvnlkK/Ru+++y9ixYwHYvHkzixcv5qmnnsJhDUDp0qULo0eP5vHHH+f/\n/u//Avs9+uijvPDCC0yePJmvv/6avn378vzzzwO+qXRGjRpF9+7deeihhwL7TJ06laysLPr06cPD\nDz8cWN6sWTMeffRRcnJyGDx4MLt37wagSZMmZGRksGTJkurf7JoI5yrKhvZPr7y3r09X7pD0h2fJ\nxl2HYt2UGtu276ikPzxL3llUUOk2Uxdvk/SHZ8nWvUfCKjMerrxPT0+Xb9cXyIZdhwLT5IuIfPHF\nF4Ery++55x555513RETkxIkTcuzYsdOmsQ+2bNky6dOnjxw9elRKS0ula9eugav1J02aJB988IEc\nP35c0tLSZMuWLSIiMmHChMDU/JW14/XXX5euXbvKoUOHZM+ePdKiRQv585//LCIi9913nzz//PMi\nInLhhRfKQw89JCIif/jDH6Rdu3ZSXFwsZWVl0qFDB9m3b5+InJqG/tixY9K7d+/A8oqGDRsmq1at\nEhGRGTNmyLhx40Ie89atW6Vfv34iIuLxeKRLly6yb9++02454D+Gzp07S0lJiRw/flw6deok27dv\nlx07dkjHjh1lz549Ul5eLsOHD5ePP/5YREQAmTlzpoiI/PznP5df/epXgfKeeuop+e1vf3tGe6J+\n5b1S8cJj0xwLQMdWSbRqmsjKwhJuHJQecpsVhSUkJyWQntqk5hV8Phl2ra5jK4OcmwWXP1vtZv4L\nJEtLS5k0aRKbNm3CGEN5eTkAQ4YM4emnn6aoqCjwy78qX3/9NVdffTVNmvheB/9cYRWtX7+eLl26\n0LlzZ8A3bb6/11BZOwCGDx9O8+bNad68OcnJyVx55ZWAb3r5VatWBbbz15mVlUXv3r1p1843Sq9L\nly4UFhaSmprKCy+8wMcffwz45jTbtGkTqamnX3u0c+dO2rRpU+1rmJGRQWpqKt999x27d++mX79+\nZ5Tld/HFFwdmJ8jMzGTbtm3s37+fiy66KFDXjTfeyL///W/GjRtHYmJioDc3YMAAvvjii0BZbdu2\nZf369dW2ryb0VJiylXIb51iMMeSkJbOiigT+isIScjqm1GkK9ljwXyD5+OOPM3z4cL7//ns+/fTT\nwASUEydOZObMmSQlJTFy5Ei++uqr6sus5jWQKoY3V9YOqN009MH7uN1u5s+fz9y5c1m4cCErV66k\nX79+ISfcTEpKCizv3bs3K1euxOsNnWe74447eOONN3j99de57bbbKj2+iu1xOp243e4qX4+EhITA\n6+nf3q+srCzkDNR1oT0WZSt2zrEA9GzXgv/k78PjlTOO4dhJNxt3H+ayzHNqV3gYPYuose55X1pa\nSocOvnvzvfHGG4HVW7ZsoUuXLtx7771s2bKFVatWkZOTw+HDofNNw4YN45ZbbmHy5Mm43W4+/fRT\nfvzjH5+2Tc+ePdmyZQsFBQVkZGQwffr0wLrK2hFJpaWltGzZkiZNmrB+/frAbY+D9erVi/z8fDIy\nMujatSu5ublMmTKFJ598EmMMmzZtYu3atYwdO5arr76aJ554gvLy8kDCvnnz5pW+ThUNGjSIn/70\np+zbt4+WLVsydepUfvKTn1S738aNG0+7YVok2O9nnzqr2e1GX8E6tmxCuUfYdejMX7bf7ziEV3xX\n6duNARB46KGHeOSRRzj//PPxeDyB9dOnTw9MUb9+/XpuvvlmUlNTOf/88+nTp88Zyfv+/ftz/fXX\n07dvX6655hp+8IMfnFFnUlISL730EqNGjeKCCy7gnHPOCZweqqwdkTRq1CjcbjfZ2dk8/vjjDB48\nOOR2o0ePZv78+YHnf/vb39i1axfdunUjKyuL//qv/wrcNiAxMZHhw4dz3XXX4XT67sGTnZ2Ny+Ui\nJycnkLwPpV27djzzzDMMHz6cnJwc+vfvHxg0UJVvvvkm8qPswknENLR/mry3rze+2SrpD8+S/UdO\nxLoptfL1Rt8U+gs3n5nk/eu/8iX94Vmy93BZ2OXFQ/JeRKTwwFFZW1xa7/UePnxYRHz3cP9//+//\nye9///t6b0N1jh07JoMGDRK3213ttv772m/cuLEeWiayfPly+dGPfhRyXV2S99pjUbZi1wsk/Tq2\n8p3LLjxw7Ix1KwtLSWuZROtmjc5YF+9iNQnlK6+8Qt++fenduzelpaVnnC6LB0lJSfzyl7+s9q6Q\na9eupVu3blx88cVRudFXKPv27eNXv/pVxMvVHIuyFX+Oxa6nwtqnJOEwoQPLisIS+tnwNBgQyLHU\nt/vvvz9w2+F4NnLkyGq3yczMDFydX18uvfTSqJSrPRZlK3bvsSQ4HbRLTqIwaJbjvYdPsKPkuG0u\njAzmz7EoBRpYlM3YeUoXv46tktge1GNZVeQbgpxj18Ci0+arCuz76VRnJX+PxaYdFgA6tWpyxqmw\nVUWlOAz0aW/PW/Lqjb5URRpYlK14vF5cDmO7Cwgr6tiyCXsOn6Cs/NQw2LU7D9GlTTOSEp0xbFld\n+HIsotFFoYFF2YzbI7a7yVewTtZ0LUUHT/Va1hYfIrNdi1g1qU5eeOEFhuX15ZGf/FeN9jtbp83P\nyMggKyuLnJwcLrvsMnbt2lXpfsGv0fz58wNTs9TVrFmzmDJlSkTKCqaBRdmK2yu2uy1xsLSWvsBS\neMCXwC85dpIdJcfJbG/PwPLSSy/x7gef8MwfX6nR6bCzddp8gHnz5rFy5Upyc3P53//930r3re1r\nFI7Ro0czc+ZMjh07c4RiXdn7E6rOOqGmQrEb/7Us/gT+2p2HAOhtw8DinzZ/0g3X8vYrL7FYp80/\nQ8Vp84MNGzaM/Px8Xn311dOGTb/yyis88MADIV+jI0eOcO2119KzZ09uvPHGwOnHL7/8kn79+pGV\nlcVtt93GiRMnAF8PacqUKfTv35+srKzAhJPGGC666CJmzZpV5XtcK+FcRdnQ/umV9/b1Px+tkgG/\n+mesm1EnXq9Xzntstvzq0zUiIvLKvzfX+Ip7v3i48j49PV3WbS2SlYUHZf/BgzptfpCK0+b7X6+9\ne/eKiMjdd98tDz30kBw5ckS6dOkiJ0+eFBGRIUOGyKpVq854jebNmyctWrSQwsJC8Xg8MnjwYPn6\n668Dr8WGDRtEROSmm24KHEt6erq88MILIiLy4osvyu233x4o75133pF77rkn5Hug0+ars0ZD6LEY\nY+jYsslpPZZzWjSq8xX3zy15jvUHIjv9ec9WPXk47+Fqt/O/I6Ulpdx+6606bX4FoabNHz58OE6n\nk+zsbJ566imaNm3KiBEjmDVrFr169aK8vJysrCwKCgrOOO68vDzS0tIA6Nu3LwUFBTRv3pzOnTvT\no0cPACZNmsSLL77IfffdB8D48eMB35T5H330UaCstm3bUlxcXMm7UHt6KkzZSkPIsQB0bNUkcJGk\nnRP3fv7AMmXKEzptfpCK0+b7zZs3jxUrVvDWW2+RkuK7dqnilPm33nprpcdW0ynzK+5TH1Pmg07p\nomymIfRYwHcty5KtBygr95C/5wgX92pb5zLD6VlES8Ubfem0+aerOG1+VQYNGkRhYSHLly8P9JzC\nnTK/Z8+eFBQUkJ+fT7du3Xj77be58MILq91v48aN9OnTp9rtasr+P/3UWaXc47X9cGOAtJZJHDnh\nZmnBAdxeIbOdPS+MPMX3njz44M902vwgwdPmV+W6667j/PPPp2XLlgBVvkYVNW7cmNdff50f/vCH\nZGVl4XA4uOuuu6qtb968eYwePTqsttVIOImYhvZPk/f2ddfby+TS38+PdTPqbM73OyX94Vny6Mer\nanSP+2DxkLwXETlw9ISsLDwoZSernxo+khratPmjR4+WuXPn1kOrRHbt2iUjRoyodH3Up803xowy\nxmwwxuQbYyaHWN/IGDPdWr/YGJNRYd0j1vINxpiR1ZVpjOlslbHJKjPRWj7MGLPcGOM2xlwbVP8k\na/tNxphJNY6uyjbcXsHZQHIsAHPW7KZpopNOrWpxj/s44u9D1vd19w1l2vySkhJ69OhBUlISF198\ncb20a/v27fzud7+LStnV5liMMU7gReBSoAhYaoyZKSJrK2x2O3BQRLoZYyYAzwHXG2MygQlAb6A9\nMNcY08Pap7IynwOeF5Fpxpi/WGX/GdgO3AL8LKh9rYApQC6+v+tvrbIO1vzlUPHO7fGS0ABOhfkD\ny97DJ8hNb4nD5nkjf46lmhxyxDWUafNTUlLYuHFjPbXGZ+DAgVErO5yffnlAvohsEZGTwDQg+Gqf\nscCb1uMPgYuNb0jHWGCaiJwQka1AvlVeyDKtfUZYZWCVOQ5ARApEZBXgDap7JPCFiBywgskXwKgw\nj1/ZjLuBJO+bNXLRqmkigG2vuK/IWH0W0TmOFeEFlg5AYYXnRdaykNuIiBsoBVKr2Ley5alAiVVG\nZXXVpn2qgfB4xbY3+QrWsaVvmGddhxpLfXcTqhJHTVG1V9e/qXACS6hPcXCtlW0TqeVVCWsfY8yd\nxphlxphle/furaZIFa8aynUsAGnW6bC69FgaN27M/v37Yx5cAqfCYtoKFQkiwv79+2ncuHGtywjn\nOpYioGOF52lA8KWa/m2KjDEuIBk4UM2+oZbvA1KMMS6r1xKqrlDtuyiorPnBG4nIy8DLALm5ufr3\nb1Nuj5dGjRrG5Vfd2zYjKcFJj3Oa17qMtLQ0ioqKiPWPpRPlHvYeOYn3YCKNXHad+l/5NW7cOHB1\nf22E8wldCnQ3xnQGduBLxk8M2mYmMAlYCFwLfCUiYoyZCbxnjPk9vuR9d2AJvl7GGWVa+8yzyphm\nlTmjmvbNAf7XGNPSen4Z8EgYx6VsqKFcIAnwXz/owpU57WmcUPsv4oSEhMCUJrG0ZOsB/uu9hbxz\n+yD6dm8d6+aoGKv2nILVc7gH3xf4OuB9EVljjHnSGOOfwOdVINUYkw88AEy29l0DvA+sBf4B3C0i\nnsrKtMp6GHjAKivVKhtjzEBjTBHwQ+Cvxpg1Vh0HgF/hC4BLgSetZaoBcjegHEvTRi66tmkW62ZE\nhD/Yu73BY2vU2SiscwoiMhuYHbTsiQqPy/B94Yfa92ng6XDKtJZvwTdqLHj5UnynuULV8RrwWpUH\noRoEt6fh5FgaEv8QcI9XzzIrndJF2Yzb68XZAK5jaWj8PZZyjwYWpYFF2UxDGm7ckPh7kdpjUaCB\nRdlMQxpu3JD4JwbVHIsCDSzKZnw5FsNzS57jsf88Vv0Oql74e5HaY1GggUXZjNsrHGIz76x7hxmb\nZ7CldEusm6SoMCpMcywKDSzKZtxeD+vK3qVV41YkOBKYtn5arJukgASn76vErT0WhQYWZTOexivZ\n79nAT/r9hFEZo5iRP4MjJ4/EullnPWfgVJjmWJQGFmUjJzwnkFazSHGmc3W3q5nYayLH3MeYsbm6\nyRlUtLl0uLGqQAOLsgUR4ZVVr2ASDjKg+SScDid9Wvchu3U209ZPwyv6SzmWnJq8VxVoYFFxb+PB\njdz+z9v566q/Un4oi05J2YF1E3pOoOBQAYuKF8WwhUpzLKqihjFNbD05fPIwj3/zeKybcVY54TnB\nwuKFNEtsxqODHmPyG01w9jr1e2hkxkh+u+y3PLX4KXq07FFFSZVr7GrM5IGTSWmcElhW7innmSXP\ncKCsZtPOjeo8ilEZ8XOfuVdXv0qf1n0Y1G5QVOs5NSpMe45KA0uNeMXL9sPbY92Ms851513Hf+f8\nN01cLZjM54GL8QASnYnc2+9e3l3/bq3eGxEhvySfjBYZ3JVzV2D5nG1z+GDjB3RO7ozLEd7HZP/x\n/Xy35ztGdBxBojOxxm2JtPyD+fxh+R/oltKNj676CGOiN2OBKzAJpfZYlAaWGklulMxHV30U62ac\ntY6f9ACcMaXLNT2u4Zoe19S63Lu+uIsPNnzA7Vm3k+BIAGDquqlktMjgk7Gf4DDhnTFeULyAH3/x\nY+YUzOHKrlfWuj2RMnX9VADyS/JZumspee3OmNs1YowxOB1GcywK0ByLspFyayhrpO/HMrHXRPYc\n38OX274E4Pt937Nq3yom9JwQdlABGNxuMBktMgJf6LF06OQhPt3yKZd3vpyURin10ianw2iPRQEa\nWJSNeKyhrJGehPL89ueT1iwt8OU7df1UmriaMLbr2BqV4zAObuh5A6v3rWb13tURbWNNzcifwXH3\ncW7pfQvju4/nq8Kv2HlkZ1TrdDmM5lgUoIFF2Yj/17DTGdk/W6fDyYSeE1i+ZzkLihfwj63/4Kqu\nV9EsseY34bqq61U0cTWJaa/FK16mrZ9G3zZ9yUzN5Przrgdg+obpUa1XeyzKTwOLsg3/zLkJUZg2\nf1y3cSS5kvjZv37GSe9Jbuh5Q63KaZbYjLHdxvKPgn+w//j+CLcyPN/s+Ibth7czsZfvDuLtm7Xn\norSL+Pumv3PCcyJq9SY4HZpjUYAm75WN+Cc4jMY975MbJTO6y2g+3Pghg9sNpktKl1qXNaHnBKau\nn8pzS54jp21OBFsZntlbZtM6qTWXdLoksGxir4l8VfgVv17y6zodW1Wk+QY2lTXn3XUrz1jnMi6u\n6HIFzRObR6VuVbl1+9fRyNkoau97KBpYlG34fw27onQHyR/1+hGzNs/i1t631qmcLslduCjtIj4v\n+JzPCz6PUOtq5r7+95HgTAg8zzs3j8zUTN7f+H70Km0Ja07AmiWhVxcdKeLB3AejV786w0nPSe6a\nexctElswY9yMGg1GqQsNLMo2/Ofvo3Wjr64pXVk0cRFOh7POZf3fiP/j0IlDEWhVzRljSG6UfMay\nd654h6Mnj0at3tF//Jp+HVN4alzWGet+ufCXfLTpI/6773+T5EqKWhvU6eYUzOFA2QEOlB1g0c5F\nDG0/tF7q1cCibMOfY4nmrYkjEVTAN0Ks4pX88SDBkRDVNiWa5hhpFrKOG3vdyNztc/l86+eM7z4+\nam1Qp5u2fhrpLdI5fPIwU9dNrbfAosl7ZRvRzLGouvNdIBl6uPGAcwbQvWV33lv3HiKa4K8P/uux\nbuh5A9d0v4Z/Ff2LosNF9VK3BhZlG/4cS0KEhxuryPBdxxI6aBhjmNhzIhsObuC7Pd/Vc8vOThWv\nx7ruvOtwGEfUh5z76SdU2YY7Slfeq8hwOaue0uWKzr5RYe+tf68eW3V22n98P59v/TxwPda5Tc9l\nRKcRfLTpI467j0e9fs2xKNtwR+nKexUZToeD8ioCS5OEJozvNp531r3D8t3LdehxFH26+VPKveXc\n0OvU9VgTe07ki21fMHvL7DrNrRcODSzKNvy/hrXHEp9cVeRY/K7veT1vr3ubSf+YVE+tOnsNbjeY\nLsmnrl0J5LnWv8f47uOjO9t11EpWKsICw401xxKXqsqx+HVs3pG3Ln+L3Ud311Orzl79z+l/2nNj\nDI/kPULThKZRDSqggUXZSH0MN1a153Iaysqrn4Qyp00OtKmHBqkzDDx3YL3Uoz/9lG3ocOP45nQ4\ndBJKBYQZWIwxo4wxG4wx+caYySHWNzLGTLfWLzbGZFRY94i1fIMxZmR1ZRpjOltlbLLKTKyqDmNM\nhjHmuDFmhfXvL7V9MVR80+HG8S0hjByLOjtU+wk1xjiBF4HLgUzgBmNMZtBmtwMHRaQb8DzwnLVv\nJjAB6A2MAl4yxjirKfM54HkR6Q4ctMqutA7LZhHpa/27C9UglWvyPq45w8ixqLNDOD/98oB8Edki\nIieBaUDwHZDGAm9ajz8ELja+7NBYYJqInBCRrUC+VV7IMq19RlhlYJU5rpo61FnCozmWuOZy6v1Y\nlE84gaUDUFjheZG1LOQ2IuIGSoHUKvatbHkqUGKVEVxXZXUAdDbGfGeM+Zcx5gdhHJOyocB1LFGa\n3VjVjcuh92NRPuGMCgv1KQ7+66lsm8qWhwpoVW1fVR07gU4ist8YMwD4xBjTW0ROm1rWGHMncCdA\np06dQhSl4l20ZzdWdeNymMDIPXV2C+cTWgR0rPA8DSiubBtjjAtIBg5UsW9ly/cBKVYZwXWFrMM6\nzbYfQES+BTYDPYIPQkReFpFcEclt00bHOtqRW3MscU1zLMovnMCyFOhujdZKxJeMnxm0zUzAfynt\ntcBX4pvCdCYwwRrR1RnoDiyprExrn3lWGVhlzqiqDmNMG2swAMaYLlYdW8J/CZRdeDyaY4lnmmNR\nftWeChMRtzHmHmAO4AReE5E1xpgngWUiMhN4FXjbGJOPr6cywdp3jTHmfWAt4AbuFhEPQKgyrSof\nBqYZY54CvrPKprI6gGHAk8YYN+AB7hKRA7V/SVS8ckf5DpKqbjTHovzCuvJeRGYDs4OWPVHhcRnw\nw0r2fRp4OpwyreVb8I0aC14esg4R+Tvw92oPQtme5ljim9NhKPdojkXplffKRnQSyvjmm4RSeyxK\nA4uyEZ02P765nDqli/LRwKJsw+314jDg0MASl7THovw0sCjbcHtF8ytxzGkFFr2nvdJPqbINj1d0\nRFgcS7DeGz0dpjSwKNtwe0QT93HMafUm9XSY0sCibMPt9WriPo753xsdcqw0sCjbcHsl8KtYxR9/\nb1J7LEo/pco2PB4JnMdX8UdzLMpPA4uyjXKvV3Mscczfm9SJKJUGFmUbHq9ojiWO+d8bnTpfaWBR\ntuH2Ci69333c8g8F1xyL0k+psg2PR3ss8czp0ByL8tHAomzDrTmWuObSHIuyaGBRtuHWHEtcc2qO\nRVk0sCjb8GiOJa4laI5FWfRTqmyj3KOnwuKZM3DlvQaWs11Yd5BUPmXlHr5YuzvWzThr7Tl8gnNb\nNI51M1Ql/DmWrzftpbjkeIxboyrq2zGFjq2a1Ft9Glhq4OgJNz+Z+l2sm3FW69M+OdZNUJVo3TwR\ngD/M3RTjlqhgme1a8Nm9F2BM/fT4NbDUQHJSAnMfGBbrZpzV6vNXl6qZnue2YMHkERw76Y51U1QF\nX6zdw3P/WM+ybQcZmNGqXurUwFIDLqeDbm2bx7oZSsWt9ilJsW6CCtI+JYk/z8/nzQUF9RZYNHmv\nlFINWJNEF9flduQf3+9i96GyeqlTA4tSSjVwNw1JxyPCu4u21Ut9GliUUqqBS09tyvDz2vLeku2c\ncHuiXp8GFqWUOgtMGprBviMn+Xz1rqjXpYFFKaXOAj/o1prOrZvy5sKCqNelo8KUUuos4HAYHhvd\nC5fTgYhE9ZoWDSxKKXWWuLjXOfVSj54KU0opFVFhBRZjzChjzAZjTL4xZnKI9Y2MMdOt9YuNMRkV\n1j1iLd9gjBlZXZnGmM5WGZusMhNrW4dSSqn6V21gMcY4gReBy4FM4AZjTGbQZrcDB0WkG/A88Jy1\nbyYwAegNjAJeMsY4qynzOeB5EekOHLTKrnEdNX0hlFJKRUY4PZY8IF9EtojISWAaMDZom7HAm9bj\nD4GLjS8zNBaYJiInRGQrkG+VF7JMa58RVhlYZY6rZR1KKaViIJzkfQegsMLzImBQZduIiNsYUwqk\nWssXBe3bwXocqsxUoERE3CG2r00dkff5ZNi1OmrFK6VUVJ2bBZc/G9UqwumxhBqTFnwnn8q2idTy\n2tRxegONudMYs8wYs2zv3r0hdlFKKRUJ4fRYioCOFZ6nAcWVbFNkjHEBycCBavYNtXwfkGKMcVm9\nlorb16aOABF5GXgZIDc3t/a3uItypFdKKbsLp8eyFOhujdZKxJconxm0zUxgkvX4WuArERFr+QRr\nRFdnoDuwpLIyrX3mWWVglTmjlnUopZSKgWp7LFY+4x5gDuAEXhORNcaYJ4FlIjITeBV42xiTj68X\nMcHad40x5n1gLeAG7hYRD0CoMq0qHwamGWOeAr6zyqY2dSillKp/xvej/+ySm5sry5Yti3UzlFLK\nVowx34pIbnXb6ZX3SimlIkoDi1JKqYjSwKKUUiqiNLAopZSKKA0sSimlIuqsHBVmjNkLbKtDEa3x\nXcxpdw3lOECPJV41lGNpKMcBdTuWdBFpU91GZ2VgqStjzLJwhtzFu4ZyHKDHEq8ayrE0lOOA+jkW\nPRWmlFIqojSwKKWUiigNLLXzcqwbECEN5ThAjyVeNZRjaSjHAfVwLJpjUUopFVHaY1FKKRVRGlgq\nYYwZZYzZYIzJN8ZMDrG+kTFmurV+sTEmo/5bGZ4wjuUWY8xeY8wK698dsWhndYwxrxlj9hhjvq9k\nvTHGvGAd5ypjTP/6bpUEk04AAANnSURBVGO4wjiWi4wxpRXekyfqu43hMMZ0NMbMM8asM8asMcb8\nNMQ2tnhfwjwWu7wvjY0xS4wxK61j+WWIbaL3HSYi+i/oH76p/DcDXYBEYCWQGbTNfwN/sR5PAKbH\nut11OJZbgD/Fuq1hHMswoD/wfSXrrwA+x3dX0cHA4li3uQ7HchEwK9btDOM42gH9rcfNgY0h/r5s\n8b6EeSx2eV8M0Mx6nAAsBgYHbRO17zDtsYSWB+SLyBYROQlMA8YGbTMWeNN6/CFwsTEm1G2SYy2c\nY7EFEfk3vnvxVGYs8Jb4LMJ3N9J29dO6mgnjWGxBRHaKyHLr8WFgHdAhaDNbvC9hHostWK/1Eetp\ngvUvOKEete8wDSyhdQAKKzwv4sw/sMA24ruNcimQWi+tq5lwjgXgGus0xYfGmI4h1ttBuMdqF0Os\nUxmfG2N6x7ox1bFOpfTD9+u4Itu9L1UcC9jkfTHGOI0xK4A9wBciUun7EunvMA0soYWK2sHRPpxt\n4kE47fwUyBCRbGAup37F2I1d3pNwLMc3fUYO8Efgkxi3p0rGmGbA34H7RORQ8OoQu8Tt+1LNsdjm\nfRERj4j0BdKAPGNMn6BNova+aGAJrQio+Ks9DSiubBtjjAtIJj5PbVR7LCKyX0ROWE9fAQbUU9si\nLZz3zRZE5JD/VIaIzAYSjDGtY9yskIwxCfi+iN8VkY9CbGKb96W6Y7HT++InIiXAfGBU0KqofYdp\nYAltKdDdGNPZGJOIL7E1M2ibmcAk6/G1wFdiZcHiTLXHEnS++yp855btaCZwszUKaTBQKiI7Y92o\n2jDGnOs/322MycP3Wd0f21adyWrjq8A6Efl9JZvZ4n0J51hs9L60McakWI+TgEuA9UGbRe07zBWJ\nQhoaEXEbY+4B5uAbVfWaiKwxxjwJLBORmfj+AN82xuTji/ITYtfiyoV5LPcaY64C3PiO5ZaYNbgK\nxpip+EbltDbGFAFT8CUlEZG/ALPxjUDKB44Bt8ampdUL41iuBf6fMcYNHAcmxOkPl/OBm4DV1vl8\ngP8BOoHt3pdwjsUu70s74E1jjBNf8HtfRGbV13eYXnmvlFIqovRUmFJKqYjSwKKUUiqiNLAopZSK\nKA0sSimlIkoDi1JKqYjSwKKUUiqiNLAopZSKKA0sSimlIur/AxBaUYX9rTVWAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f12c67109d0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"pl.plot(xs, [bestofk_time(sp.digamma, x, R) for x in xs], \n",
" label='scipy.special.digamma')\n",
"pl.plot(xs, [bestofk_time(fast_digamma, x, R) for x in xs], \n",
" label='fast digamma (Cython)')\n",
"pl.plot(xs, [bestofk_time(py_fast_digamma, x, R) for x in xs], \n",
" label='fast digamma (Python)')\n",
"pl.legend(loc='best');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Holy cow! That's a lot faster!"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Accuracy\n",
"\n",
"Let's compare the accuracy of the two implementations against the very high-precision implementation in `mpmath`, which we'll treat as the gold standard."
]
},
{
"cell_type": "code",
"execution_count": 44,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"9.85083473332755\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEKCAYAAAA8QgPpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzsnXd4HNW9sN+jsuq2JEuWbHV3G/de\nsY3pEDCmYwMGAqGEkBvSuAkfuUlIclO4CQSS0JtNtenVgMG9V9yNVS3JsizLalY/3x9nR7ta7a52\nV9sknfd59Ozu7MycI2lmfufXhZQSjUaj0WjcJSTQE9BoNBpN90QLEI1Go9F4hBYgGo1Go/EILUA0\nGo1G4xFagGg0Go3GI7QA0Wg0Go1HaAGi0Wg0Go/QAkSj0Wg0HqEFiEaj0Wg8IizQE/AlSUlJMjs7\nO9DT0Gg0mm7F9u3by6WUyZ3tF9QCRAgRA6wBHgFWA08BjcDXUsplnR2fnZ3Ntm3bfDtJjUaj6WEI\nIfJd2c9vJiwhRIYQYrUQ4oAQYp8Q4gEXDvsF8Kb5/SLgbSnlncAVPpuoRqPRaFzCnxpIM/CglHKH\nECIO2C6EWAWEAn+02fd2YCywH4g0b0sH9prft/hhvhqNRqNxgt8EiJSyBCgxv68WQhwA0qSUq4DL\nbfcXQswHYoBRwFlgJUqI7EI7/zUajSbgBMQHIoTIBiYAmx3tI6X8lXnfpUA5ygfyTyHEZcAHTs59\nF3AXQGZmpremrNFoNBobhL/7gQghYoFvgEellCt9OdbkyZOldqJrNBqNewghtkspJ3e2n19NQUKI\ncGAFsMzXwkOj0Wg0vsWfUVgCeA44IKV8zF/jajQajcY3+FMDmQXcDJwnhNhl/rnUj+Nrujkf7C5m\nf3FVoKehseLIiWre2VlEfZMOjOyN+DMKax0g/DWeN2ltlYSEdMup9xhKz9Rz/2s7AZgzNIkfnDuY\nWUP6oRRbTaB44qujvL+7mN9/eIBbZ2Zz8/QsEmJMgZ6Wxk8EdSZ6oNlbdIbff7SfnQWVXDQ6lRum\nZDBjUD8tTAJAeU0DABeOSmFnYSVLntvMOQP7cNe5g7h0zADCQ3VkdyA4VdtAVr9oBifH8tiqw/zr\n6++4bnI6d8weRGa/6EBPT+NjtACxQ3HlWf762SFW7jxOYoyJK8YPZNX+E3ywu5isftHcN38Iiyak\nEaYfWn6jsq4JgNtn5zAhM553dx7n6TXHeOD1Xfz500PcNiubG6ZmEhuhL2l/crq2iSHJsTy3dAqH\nSqt5Zu0xlm8p4JVN+VwyegB3njuI8RnxgZ6mxkfou80OT319lA/3lnD33MHcO38wfSLDqW9q4bN9\npTyz9hg/f3sP//r6O360YAiXjx2oV79+4HRdIwAJ0SYiwkK5fkom107KYPWhMv6z5hi//+gA//jy\nCIunZXHbrGxS+kR2ckaNN6isa2TkgD4ADE+N46/XjuNnFw3nhfV5LNucz0d7S5iancid5w5iwYj+\nWnvvYfg9D8SfeJoHUl7TwNnGFjISO6rgUko+33+C/1t1mIOl1QzoG8ktM7K5cWoG8dHa9usrXtmY\nx8Pv7WPLfy+gvx3hsKuwkmfWHuOTvSWEhgi+N24gd84Z1PZw0/iGkQ9/yuJpmfz68lEdvqtpaOaN\nrYU8vy6X45VnGZQUw/fnDGLRxDQiw0MDMFuNq7iaB6IFiIe0tkq+PFjGC+tz2fDdKSLDQ7h2UgZ3\nzM4hOynGJ2P2Zh7/8giPrTrM4d9fginMscZXWFHHc+tyeXNbIXWNLcweksT35+Qwd1iydrh7mfqm\nFkY8/Ck/u2g4980f4nC/5pZWPv62lKfXfMe3x6voF2Pi5hlZ3Dw9i36xEX6cscZVtADBf5noB0ur\neH5dLu/uLKaptZWLRqXyg7mDmJCZ4POxewv/88E+3txayL7fXuzS/mfqmli2JZ+XNuRxoqqBof1j\n+f6cHK4cr1e/3qL0TD3T//glj141msXTsjrdX0rJpmMVPLP2GF8dLCMiLIRFE9O5Y3Y2Q/rH+WHG\nGlfRAgT/lzIpq6rnpY15vLIxn6r6ZqblJHL3vMHM06vfLvNfb+xiS24F6395nlvHNTa38tHeYp5Z\nk8v+Esvqd8n0LJL06rdLHCip4pJ/rOWpxRO5dMwAt449WlbNc+tyWbHjOI3Nrcwbnsz3Zw/SodlB\nghYgBK4WVk1DM69vKeDZtbmUVtUzLCWWO2br1W9XWPrCFsprGvjw/jkeHS+lZOOxUzy3NpcvD5Zh\nCgth0YQ0vj8nR69+PWTjd6e48ZlNLL9zGjMHJ3l0jvKaBpZtKuCVTXmU1zQyIjWOO2bncMX4gUSE\n6XslUGgBQuCLKTY2t/L+7mKeW5fLAavV7y0zsknUyVZuceWT6+kTGcYrd0zr8rmOltXw/PpcVmwv\noqG5lfnDk7lDr37d5pO9JdyzbAefPDCny8EK9U0tvL+7mOfX5XKwtJqk2AhumZHFTdMytaYYALQA\nIfACxMB29Ws43G+fnUOOdri7xNy/rGZsejxP3DjBa+c8VdPAq1ar3+Epcdw+O1trii6yfHMB//3O\nXjY9tIDUvt4Jm5ZSsv7oKZ5dd4yvD53EFBbCwvEDuX12DiNSdUSdv3BVgOg8ED8ghGDm4CRmDk7i\nyAmVbPXG1kJe2ZTPnKFJ3DIjm/NG9CdUx8g75HRtIwnR4V49Z7/YCB44fyg/mDuID8ya4i9W7OV/\nPz3ETVMzuXlGls4ncYKRmxPvxf+LEILZQ5OYPTSJo2U1vLghlxXbj/PmtiJmD0nijtkqok7nkwQH\nWgMJEGXV9by+pZBlm/M5UdVAekIUt87I5ropGfSN8u6DsrvT3NLKkF99wo8WDOUnFwzz2TiGpvjC\n+jy+OHCCUCG4dMwAbpuVrSPq7PDoR/t5dVMBB37nWmScp1TWNbJ8S0FbRN3g5Bhum5XDoolpRJv0\nGtgXaBMWwS1ADJpaWvl83wle2pDHlrwKosJDWTQxjdtm6dBGg1M1DUz6/Rc88r1R3DYrxy9j5p+q\n5aUN+by1rZDqhmbGZ8Rz++wcLhmdqisPmPnpW7vZcLScDQ8t8Mt4TS2tfLy3hOfW5bKn6Ax9o8K5\nYUoGS6Zn2U361XiOFiB0DwFizbfHz/Dihjze311MY3Mrc4ZaVPbe7Nw9WlbD+Y99w9+vH8/CCWl+\nHbumoZkV24t4YX0ueafqGNg3kltnZnPDlEz6etmk1t34/ktbKa6s5+MHPIuM8xQpJdvyT/P8ulw+\n338CKSULRqawdGY2MwfrQAhvoAUI3U+AGJyqaWD5ZlWQrqxaJ8Fty6vgmn9v5MXbpjBveP+AzKG1\nVfLVwTKeW5fLxmOniAoPZeGENJbOzGZ4au/UFK/+1wYiw0NY9v3pAZtDceVZlm3O57UthVTUqkCI\n22Zls3BC77xXvIUWIHRfAWLQ2NzKh3uKeWatCgNOijWxeJpKgkuO6z2hjV/sP8H3X97Ge/fNYlwQ\nVHbdX1zFyxvzeGfncRqaW5k+KJGbp2dz4Tkpvcq8dd7fvmbkgD48edPEQE+F+qYWPthdzAvr89hf\nUkV8dDg3Ts3k5ulZDIyPCvT0uh1agND9BYiBlJKN353i2XW5fHWwDFNoCFeMH8iNUzOZmBnf41X2\nt7YV8rO397DmZ/ODqsfE6dpG3thWyKub8ik6fZaUPhEsmaZyF3pDjaeJv1vFpWNS+f3CMYGeShtS\nSjbnVvDi+jw+31+KEIKLzknh5unZTB+U2OPvFW+hw3h7EEIIZg5JYuYQFdr4wvpc3tl5nLe3FzEs\nJZYbp2ZyzaR04iJ7pk3e6AUSHxNcv19CjIm75w7mzjmD+PpQGS9uyONvqw7zxOqjXDU+jaWzsnts\nNeDWVkllXSMJQVaBWgjB9EH9mD6oH4UVdbyyKZ83thby8d5ShqXEcvOMbBZNSCNG943xCkGtgQgh\nYoA1wCPAd8ADQBLwpZTyX50d31M0EHvUNDTzwe5iXt9SwO6iM8RGhHHt5HSWzswmq1/PSk7886cH\neXrNMY48eknQryCPnKjmhQ15rNxRRH1TK1NzElk6M5sLR6X0qAZkZ+qaGPfbz/n1ZSP5/pxBgZ6O\nU842KvPWy5vy+PZ4FXGRYVw3OYNbZmT1uHvFWwSdCUsIkQG8DKQCrcDTUsp/dHLMb4FaYJ+U8kPz\nthDgGSnlHZ2N2ZMFiDW7Cyt5YX0uH+0toalFcu6wZBZPy2TBiP494qH10Mq9rNpfyrZfXxDoqbhM\nZV0jb24r5OWNyryV2ieSG6dmcuPUDLv9TLobeeW1zPvr1/zt2nFcPSk90NNxCSklOwoqeXFDHp/s\nLaFFSuYP78+tM7OZMyRJJydaEYwCZAAwQEq5QwgRB2wHFgKhwB9tdr8dGIvSNiKBcinlh0KIK4Bf\nAv+UUi7vbMzeIkAMyqrqWb6lgDe2FlJypp7UPpEsmZ7JjVO7t03+nle3c6Sshi9+MjfQU3GbFnP0\n1ssb81h7pJywEMHFo1O5dWY2k7MSgl6jcsTOgtNc9dQGnl86mfNGpAR6Om5zoqqeZZsLWL65gPKa\nBnKSYlgyPYtrJqXrRF6CUIB0GFiI91CCYJWD7x8FYoBRwFngKillq/m7j6SUlzk47i7gLoDMzMxJ\n+fn5vph+UNPc0srqQyfbHlpGPaGbp2czJr1voKfnNjc8vZGWVslbd88M9FS6xLGTNby6qYC3thdS\nXd/MyAF9uHVGVrcMOV19sIzbXtzKyntnMrEbZ+k3NLfw6belvLQhjx0Flebw7IHcMqPn+q9cIagF\niBAiG+XbGC2lrOpk36VAOVADLAIigD1Syic7G6e3aSD2sLXJj0nry+JpmVw5Po0oU/d4aF389zVk\nJEbzzC2dXs/dgrrGZt7dWczLG/M4WFpNgjnkdEk3CjlduaOIn7y5m9U/nddjCoJ+e/wMr2zM591d\nKjx7anYit8zM4qJzel/1gaAVIEKIWOAb4FEp5UpfjuWxANn9Ony3Ghb9x/uTChBnzjbxzo4ilm8p\n4PCJGlUGYmoGN0/PIj0heEJj7THtD18wd1gyf75mXKCn4lWMkNMX1ueyav8JAOYP78+NUzOZNzw5\nqP1Xz63L5Xcf7mfX/7uA+CCLxOoqlXWNvLWtiFc25VNQUUdKnwgWT8vihqkZ9I/r/v4rVwhKASKE\nCAc+BD6TUj7m6/E8FiCrHoFNT8HDJ70/qQAjpWRLbgUvbczj029LAbhgVApLZ+YEZZy8lJLhD3/K\nbTOzeejSkYGejs8orKjjja2FvLGtkJPVDaT2ieSGqRncMCXTa6XSvcnfPj/Ek6uPcvTRS3us87m1\nVfL14TJe2pDPN4dPEh4quHBUKjdMzWDW4J7tdA+6PBChnkzPAQf8ITy6hCkWWhqhuRHCetbqSgjB\ntEH9mDaoH8crz/LKxnxe31rAZ/tOMCI1jusmZ3D5uAFBs9I629RCY3Nrj1vl2pKRGM1PLxrOA+cP\n5csDZSzfUsDfvzjCE18d5fyRSis5d2jwlDE/XddI36jwoJmPLwgJEZw3IoXzRqRw7GQNyzYXsHJH\nER/tLSEjMYol07K4fkpGj782neHPKKzZwFpgLyqMF+C/pZQf+2pMjzWQTf+CT38JP8+F6ETvTyzI\nqG9q4f1dljj5EAGzhiRx09RMLjwnNaB9So5XnmXWn77iT4vGcMPUzIDNIxDkn6pl+ZYC3t5WxKna\nRtITorhxaibXTc4IeCmbHy7fwf7iKr766byAzsPfNDS38Nm+EyzblM/m3Aoiw0NYNDGdJdOyGDWw\n5zjdg04DkVKuA7rHcsVkdgo21vYKARIZHsp1UzK4bkoGR05U896uYt7ZeZx7lu0gPSGK22blcN3k\nwGS6n641mhb1vlVeVr8YHrpkJD+5YBir9p9g+eYC/vLZIf7+xWEuHTOAJdOzAhYKXFnX5NVGUt2F\niLBQrhg3kCvGDeRASRUvbchjxfYilm8uYFx6X26YmskV4wb2mkz3oM5E7yoeayDfroC3b4d7N0P/\nEd6fWDegpVWyav8Jnlt3jK15p4kxhXLVxDRumZHNsBT/VZ9dd6ScJc9t5o27pjNtUD+/jRusHC2r\nYdnmfN7eXkR1fTMjUuNYPC2ThRPS/CrgL3t8Lal9Inlu6RS/jRmsVNY18s7O47xmDlCJjQjjqglp\nLJ6e2W3b8AadBtKtMJkfkI01gZ1HAAk1J7xdPDqV3YWVvLwxnze3FfHqpgKm5SRyywz/VJ812qYm\nxPQ+DcQeQ/rH8sj3zuFnFw3nvV3FvLopn4ff28efPjnIFePTWDwtk9Fpvs/1qaxr6rYPR28TH23i\ntlk5LJ2ZzY6C0yzbXMAb21TL6slZCSyZnsUlY1KJCOseYfPuoAWIPdpMWL1XgFgzLiOev2XE86vL\nRvLGVlV99r7lO0jpE8H1UzK5bnK6z0KBK33Qd7snEG0K48apmdwwJYPdRWd4dVM+7+ws4rUtBYzL\niGfxtEy+N3agz3J9Ttd5v0d9d0cIwaSsRCZlJfLwZaNYsaOIZZsL+PEbu/jdhyaumZzOTVMze1T9\nLW3CskfxLnh6LtywHEbYTXjv1bS0Sr4+VMbLG/NZc0SFOs8eksT1UzK4YFSKV1daj395hMdWHebw\n7y/BFBa8eRHBwJm6JlbuVA+to2U19IkM45pJGdw0LZMh/WO9Nk5DcwvDf/0pP7toOPfNH+K18/ZE\nWlsl678r59VN+XxxoIyWVsnsIUncNC2T80emBO01rU1YXcFkvtkaawM7jyAlNESwYGQKC0amUHS6\njre2FfHWtkJ+uHwn8dHhLByfxnWTM7wSlXK6rpHYiLCgvdGCib7R4W2mlC25Fby6uYBXNuXx/Ppc\npg9KZPG0LC48p+sCvq28vtZAOiUkRDBnaDJzhiZzoqqeN7cW8tqWAu5dtoOkWBNXT0rnximZZHfT\nbH4tQOwRYRYgDdWBnUc3ID0hmv+6YBg/WjCU9UfLeXNbIcs3F/DihjxGDejD1ZPSWTh+oMfFHHtr\ntE9XsM71OVk9ire2q//J/a/tJDHGxKIJadww1XOtpM0v1Qsj47pCSp9I7l8wlHvnD2HNkZO8trmA\nZ9fm8vSaY8wblszSWTndriqwFiD2sA7j1bhEaIjg3GHJnDssmdO1jXywp5i3txfxuw/387+fHOSy\nsQO4dWY2491sSXs6CJsWdSeS4yK4d94Q7j53MGuOnOSNrYW8uCGPZ9flMi0nkZtnZHHhqFS3NLzT\ntVoD6QqhIYL5w/szf3h/TlTVs3xzAcs2F3Dr81vI6hfNVRPSWDg+rVtoJVqA2CPc7BDWAsQjEmJM\n3DIjm1tmZHP4RDXLNxfw9vYi3tl5nPEZ8dw5ZxAXj3YtQfG01kC8QkiIYN7w/swb3p+T1Q28vb2I\n5Vvy+eHynSTFRnDd5HSunZzhUmHESq2BeI2UPpH81wXDuHf+YD7eW8KbW4v4x5dH+PsXR5iYGW+u\nDDGQ2CDNK9FOdEc8OgAm3w4XPerdSfVSquubWLnjOC+szyXvVB2ZidHcMTuHhRPSnPZfmPuX1YxL\nj+fxGyf4cba9g9ZWyTdHTrJsUz6rD52kpVUyJTuBaydlcOnYAQ4fWss3F/Df7+xl00MLgrJOV3en\nuPIs7+0qZsWOIo6W1RAVHsplY1Xi6Lj0vn5JHNVO9K5iitVhvF4kLjKcW2dms2R6Fqv2l/KfNcd4\n5P19PPrxAS46J5VrJqUze0hSB63kdK0OF/UVIVamlLKqelbuPM6bWwv5+Yo9PPL+Pi4Zncp1UzKY\nltO+yOZpHVrtUwbGR3HPvMHcPXcQOwsreXNrIR/sVibh0Wl9WDIti+8FSbZ74GcQrJhitAnLB6gE\nxQFcdE4q3x6vYsWOIt7ddZwPdheTFh/FtWZTSlp8FM0trVTVN/fKMib+pn+fSO6eO5gfnDuIHQWV\nvL29iA/3FLNy53GG9o9lyfQsrpqYRp/IcCrrGokMD+l2TbC6G0IIJmYmMDEzgV9dNpJ3dx7n1U0F\n/HLlXn774X4uHzuA66dkMDEzcJ0ttQnLEf+aDfEZcONr3p2UpgMNzS18sb+M17cWsO5oOQCTMhM4\nd1gyj606zG++N4qls3ICPMvex9nGFj7Yo7Ld9xSdaTOlFJ2uI/9UHRsfWhDoKfY6VF/307y5tYgP\n9hRT19iiytlMz+KqCWle85UEZT8Qf9MlAfL8xRASBks/9O6kNE4prKhj5Y7jfPJtCQdLVRj1P2+a\nwOVjBwZ4Zr2b3YWVvL61gPd3FVPb2MKoAX34+IE5gZ5Wr6a2oZkPdhfz6uZ8vj1eRYwplEvHDGDR\nxHSm5SR2KRxYCxC6KEBevRrqKuCu1d6dlMZlcstr2ZF/msvGDtDmkiChtqGZj/eWkJYQxczBSYGe\njgallewuOsNrmwv4aG8JNQ3NpMVH8cdFYzh3WLJH59RO9K5iioHKwkDPoleTkxTTY/pt9xRiIsK4\ndnJGoKehsUIIwfiMeMZnxPObK87h8/2lrNxxnAF+iJDTAsQRpljtRNdoNN2KKFMoV45P48rxaX4Z\nTxcYcoQpFhp1KRONRqNxRFALECFEjBBiuxDiciFEiBDiUSHEE0KIW30+uBHG24N9RBqNRtMV/CZA\nhBAZQojVQogDQoh9QogHXDjsF8Cb5vdXAmlAE1Dkq3m2YYqB1mZoafT5UBqNRtMd8acPpBl4UEq5\nQwgRB2wXQqwCQoE/2ux7OzAW2A8YnqDhwEYp5X+EEG8DX/p0tkZJ94YaCPOskqxGo9H0ZPwmQKSU\nJUCJ+X21EOIAkCalXAVcbru/EGI+EAOMAs4CK4B689ctPp+wUdK9sQZidC9ujUajsSUgUVhCiGxg\nArDZ0T5Syl+Z910KlANfAU8IIeYAa5yc+y7gLoDMzEzPJ6lLums0Go1T/C5AhBCxKG3ix1LKqs72\nl1K+aPXxDhf2fxp4GlQioYfT1F0JNRqNphP8GoUlhAhHCY9lUsqV/hzbbdoEiA7l1Wg0Gnv4MwpL\nAM8BB6SUj/lrXI/RJiyNRqNxij81kFnAzcB5Qohd5p9L/Ti+e2gBotFoNE7xZxTWOqD7dIs3WUVh\naTQajaYDQZ2JHlAirPJANBqNRtMBLUAcERYFCG3C0mg0GgdoAeKIkBDd1laj0WicoAWIM0wxOoxX\no9FoHKAFiDN0TxCNRqNxiBYgztAmLI1Go3GIFiDO0BqIRqPROEQLEGdExEKD9oFoNBqNPbQAcYY2\nYWk0Go1DOs1EF0K4WhO90pXqut0KLUA0Go3GIa6UMnkJkDgvQyKBF4GXvTCn4EH7QDQajcYhnQoQ\nKeV8f0wkKDHFqjwQKUF0nzJeGo1G4w/c9oEIIWKEEKG+mEzQYYoB2QrN9Z3vq9FoNL2MTgWIECJE\nCHGTEOIjIUQZcAgoFULsE0L8RQgx1PfTDBCedCVsrIP374eqEt/MSaPRaIIEVzSQ1cBg4CEgVUqZ\nLqVMBuYAm4A/CSGW+HCOgaOtIq8bobwFG2DHy3DkM9/MSaPRaIIEV5zo50spm2w3SikrUO1pV5hb\n1fY8PGkqVfqteq0s8P58NBqNJojoVAMxhIcQ4u/mtrQO9+lxeCRA9qrX0/nen49Go9EEEe440WuA\n94UQMQBCiAuFEOt9M60gwZOuhCcMDUQLEI1G07NxuaWtlPLXQoibgK+FEA1ALfBLn80MFfEFrAEe\nAY4BvwFOAV9KKd/25diA+wKk6SyUH1HvtQai0Wh6OC5rIEKIBcCdKMGRDPxISrnWjeMzhBCrhRAH\nzBFcD7hw2C+AN83vLwGekFLeA9zi6rhdwl0TVtkBkC2QOhZqy1REliNqyqC5oetz1Gg0mgDhjgnr\nV8DDUsp5wDXAG0KI89w4vhl4UEo5EpgO3CeEGCWEGCOE+NDmp78Q4nxgP3DCfPwrwA1CiL8A/dwY\n13PcDeM1zFcjLlevZwrt79faCv+aCav+X9fmp9FoNAHEHRPWeVbv9wohLkFFYc108fgSoMT8vloI\ncQBIk1KuAi633V8IMR+IAUYBZ4GPpZT3mZMYV7o67y7RpoG4aMIq3auEzqC58PUflBkreXjH/aqO\nQ+1J2P0anP8/EB7pvTlrNBqNn3ClmKKQUkrb7VLKErNZy+E+Ts6ZDUwANjvaR0r5K/O+S4FyIFMI\n8d8oofIXJ+e+C7gLIDPT1TqQDgiPAhECDXYEiJRw9jREJ1q2lX4LKedAQrb67MiRXn5YvdafgcOf\nwDlXdW2eGo1GEwBcMWF9JYS437YqrxDCBMwQQrwE3OrqgEKIWJTm8mNXqvdKKV+UUn4opcyTUt4l\npVwspVznZP+npZSTpZSTk5OTXZ2Wo8k6LqiYuwb+MgSO7zAGViaslNEQmwJhkXA6z/55DUd7ZDzs\neq1rc9RoNJoA4YoAOQK0AO8IIYqFEPuFEMfM228E/k9K+aIrg5kTDlcAy6SU/jFDdRVTjH0TVvlh\n5TDf9JT6XJkPDVWQOkYJnr4ZjpMJyw9BZF+YtBSOfqEc6hqNRtPNcEWAzJRSPoUq554JLAAmSimz\npJR3Sil3uTKQOQnxOeCAlPIxj2fsb0yx9gVIjdm3v+8dVffKyEBPHaNeE7KcmLCOQNJwGH+TEkJ7\n3rS/n7vUn4G9b0PtKe+cT6PRaJzgigD5TAixEUhBhc8OBDwpTzsLuBk4Twixy/xzqQfn8S+OmkrV\nnICwKGhtga3PKvOVCIH+o9T38VmOc0HKD0PSMOVgHzhROdO9weanYcUd8LdhsOw62P+eMq1pNBqN\nD3CllMmDwGKUGSsHeBjYa87leMPVgaSU66SUQko5Vko53vzzsccz9xeOfCA1ZdBvCAy/FLY9D0Vb\nIXEwmKLV9wlZUF+ptAJrzlYq4ZM8TH0ef5MSPiV7uj7Xgo2QOAhm3Kciwt68xbEfRqPRaLqIS3kg\nUspjqKKKD0spF0ophwLTgP/z6eyCAUc+kJoTEJcC0++BsxXKl5E62vJ9vDnmwFYLMRzoSWYBMvpq\nCAmHvW91bZ6trUqIDZoHF/wWvvcPtb22vGvn1Wg0Gge4nEgopTxs87lGSrnJ+1MKMiJi7Yfx1pSp\naKvs2SryCiz+D1AmLOjoSDdCeA0BEp2ohM2ZItfnVLIHdi1vv+3kQeXEz5imPkfFq1dbDUij0Wi8\nhNsdCXsd9nwgra1mAdJfRVwkqA+/AAAgAElEQVRNv0dtHzDOso+jXJDyQxBqsggYUGavJidlT2zZ\n8jS8d1977aLQnFKTMVW9RvZVr/WVrp9Xo9H4l+0vwdruE1NkixYgnWHPB3L2NLQ2KQ0EYNxNsHgF\nDLKq7BKVAKY4+yasxMEQapXDGe7AUe+IugrVavfQJ5ZthZshJhkSctTnyHjLXDUaTXCy500lQFqa\nAz0Tj/CkJ/r3fDGRoMXwgVhHMxkhvLH91WtICAw9X70aCKFMU7YayMlDFgd62xhuaiB15jDdgx9a\nthVuVuYro2VLmwbioQmrqV5pWv6muRHW/Z8aX6Pp6dSdgsZqKN0d6Jl4hCcayKNen0UwE5UAyPYr\n+TYBkuL82ISs9j6Q5gYVFZVkK0BinFfuteVshXr9brVqt1tzEiqOWcxXoOprhUV6ZsKSEh4fb0mS\ndDgPO1FmXaVgA3zxGzjyuXfPq9EEI8ZiMM9hcY2gxhMBYrcrYY/F8FVYh8MameOdCRAjF8TQXipy\nVeKgrQAJj4Emd0xYpyB5BLQ0qOivoi1qu+FAN4iM9+wB31AN1SXq3M54Ywm86XIVG9cw5lvxnXfP\nq+nZbHwS9nQxktHftLZaBEiuy50xggpPBEjvykxLNPsUTudatrmjgTTVWi6S8kPqtYMGEu26D6S1\nVWlDwy+B6CQ48KEyX4WaYMD49vtGxSstwV2M+R7f7tiM1Vin8k7y1rnnv+kMQ4Cc0gJE4yItTfDV\no/Dlb7ueOCul4xJE3qbhjFpQhprUvdQN/SDaid4ZdjUQcxZ6RJyLx5r9IG0hvEPb7xce7boJq75S\nOdBj+ishcuRztXoZML5jWfjIvp6ZsOrMJrKGKjh1xP4+RVuhtVkFExRsbP/doU/gG4cFk53jTICU\nfqt8JBqNNcU71ULtTIF63xU+/zX8YxxUn+h8365i3GeDFyg/a0n384NoAdIZEbEqusnWhGWE8DrD\nSCbMN7eOLz+iiiwafUYMTDHQfNY1p7Xhi4lOhJHfUw/54h3t/R8Gnpqw6qxqaRVttb9PwUZAqCTI\nY9+0/+6rR2H17z1r6+vIhFVVDP+ZA3ted/+cmp5N7hr1KkLhwPuen2f/e7Dxn2qBVl3i2TkqCyB/\noxJkJw8514iM+2zUFeo1r/uZsTwRIH4QzUFGQo7yXxjUlHZuvgJlqho4EVY9DG8tVReVrfYBSgMB\n1yKxjFVLdD/ImWvpmmjr/wClgXTFhIVwLEDy16vEyYypkGslQE59Byf2qvee1PiqN1f4rzlheQ9q\ndSZbtWlL05G8daoGXc65sP99z8xY5Ufh3fuUWRg8C38v2QOPT4QXLoan58GTU+Hdex3vb9xnySNU\ncdWuCJCTh+DPg9s/p/yA2wJESnmBLyYS1CRkt19NGxpIZ4SZ4I7P4bxfw8GPLEUUbTE0ElsBsvVZ\n2PDP9tuMiy4qUZmshl6oPtvTQKLiPTRhmcdImwRF2zp+39wIhVsha6YSYiV7LIJt3zvqtf85Klve\n3VBga42p4pjlfalZKFUVu3c+TWBoqve8DltTPbx2I5zY3/m+zY3KB5g9R63kK76DMheOs6axTtWN\nCw2Hq/6jtrkrQJob4J0fKMvA4rfhhtdUu4bdy1W0pD2MROCYJMiZAwWblD/HwB1BmLsG6srh1FH3\n5t1FtAnLFRKyoarIYn+vOeGaBgLqojz3Z3D3ehh7PYy5ruM+hgZi64ze+zbsfKX9NiOE1+iEOPcX\ncMlfIC6143kj49Uq3t2H+NkKZQ4YskDdjA3V7b8v2a1MblkzVftepGX1tP9dSJ8Cs3+scmAKNrg3\ndv0ZFX4M7W8GLUC6F+v/AU9MgmKXuj205+QBOPSxa9fO8e1q4ZU9G0ZcDgilhbjD9hegbB9c/ayl\nHJG7AmT1H9S9csU/YegFMOJSuPh/VXHTjx60n9dkLNSi+6n5G36Q+ip4fTH8e47rQqTUXIzV1fbb\nXsKTRMIYc1/y3kNijjKfnClUK42zp+0/sJ2RPAwWPQ3pkzp+50gDaajpaIttu+jMAqT/CJh2l/0x\nI/sCUvlJ3KHulLqo06eq39vWMWnc2JkzlJZiilUroFPfqQf9qIXqZo7oAzuX2R+j1YE5qv6MpbaY\ntQZywtxvpeq4e7+LJjAc+1oFWbx7r/uBD0YUlCvRfXnrAKEewLH91aJm/3vujVe0Vfkrhyyw1JBz\nR4AUbIYNj8PEW2DYhZbt4ZFw2WNKK1pnp+5s3Sm1WAqPhqzZatuOl+HZ81WS8Im9yjTlCsYCy17d\nPh/SqQARQoQIIW4SQnwkhCgDDgIl5nLufxFC2DHq9zCMulanc6H2pHrvignLVQwBYhuJ1VitHqhN\nZy3b6iogJEw9nDujraCim2YsQ4CkTVSfbf0g+Rug31D1NwgNVzftsW8sN+6oK1Vo8jlXqW32Lurd\ny5WNuOZk++31Z5R21zfDooE01CjbbkiY0kB0j5PgpqleaQYDxquV/Zo/u3e8YS52JTIxb41acBgL\nqpFXKA2m3EH0oD2Kd1lC4MMiVF6WqwKkpQnevQf6psNFf+j4/eD5MOZaWPdYxznVVaj7TAiITYbk\nkbDjJWWKusJsus63STBsaYaThztuM8x9QaiBrAYGAw8BqVLKDCllf2AOsAn4kxBiiQ/nGHjaBEie\n6zkg7tBmwrL55xsP3upSy7a6U8r/0VkEGFjqYbkbiWVc2NGJSlAUbbd819qqIrCyZli25cxV4b7b\nX1AaSXyG2j5+sQqv3P9uxzHyN6oVarWNSaqhSmlOiYMsGkrZfkAqQdXSYPG3aIKT4h3q/3Tuz9Q1\nsPYxOL7D9ePbNJBOHobNDVC4RWkfBiPNlZZc1ULOVqqF4UCrHKqoBNeDT0r2KA3jvIcdh/Vf+CiE\nRsC6v7ffXnfKIvgAJt6sfDl3roYJSyBuIOStb3/M1mfgqentc1VOHVV/bwg+DQTVB+R3Uso9Uso2\nY7qUskJKuUJKeTXgcmOpbklsqlI1T+dZZaF7UwNxEIVl+B6sBcjZivYXnTOMeljuRmLVnYLoBPU+\nfYrSQIxVf9l+JZCyZln2HzRXvVYWKPOVQcZU1XTLnhmr2PxAsRUG9WfUvPsNsWgghno+7GL1qs1Y\nwU2+lYnzoj+oxda797h+HRr14zqLSizaBs31ygFt0DdNLWIOudirzsi9GGArQFzUQAzt3Pp+sCUu\nRVXqtg1Nryu3RH2BagS39EOVgCyEWjDlb2ivce97RyUfWjvmjfsDlNXCj7jSkbAJQAjxd3Nfc4f7\n9FhCQlRSYEWu5WHuVQ3Ejgmrpcmyqqix1kDM2oEreNoTxDBhAaRPhtoyy4rH+uFg0P8cy/6jrrRs\nF0Kp7wUb2guKxlrVvwQsQQGg2gMbGki/wcr0VlehbpDIvsonA9qRHiy0NMP6x2GNTdJo/gZljonp\np67BhU8pbfKFS9svhhzhqg8kby1gftBaM+R85bdzRQiUmJ38AydYtkXFuyFAtkCfNCW4nNE3vWPP\nH+v7zB7Zs9S9b/gCq08ojQvah86X7lEaTkTfoNRADGqA94UQMQBCiAuFEOs7OaZLCCHmCSHWCiH+\nbX6/UAjxjBDiPSHEhZ2fwYsYobyGBhKT7L1zt2kgVjeMdeRTtY0AiUpw7bye9ASRsr2QSp+sXg9+\npByj+9+DPumWJElQAnbUQhg0X62erMkxayf5VpdKyR7lnDd+HwPD2W9oIKC0kBPfQsoYy01qq4Hs\nXKay1DX+o+IYvHCJynH66vcWc2Nri3rIWZs4B8+HxW8qDf65C53n8liXEulUgKxTUVO298Og+er6\nMhIMnVG8C/pmttfq3dFACrcqLb0z+qaphY91uZLOBIjhWDcKLR76GJCQOlb5HI3oytK90H8kRPUN\nSh8IAFLKXwOvAV8LIdYBDwK/dGcwIUSGEGK1EOKA2Qn/QGfDogRXJFAkpXxXSnknsBS43p2xu0xi\njtmEVaousLAI7527zQdipYFYXwjWkVidXXTWtPUEcUOA1Jvr8xhj9D9HaUifPQQvX6mcesMu6uiD\nufwxuMWOryNtkir7Yl1ttNjKHm4tQIzEwcg+qmcKKMfjif2qXXBsigovttZAmurh/fvhU7cuxd5H\nYx28fYd3Es0OfAj/mq0ihC75i6pGsPVZ9V3pXmVGsTXpDD4Pln6gFkbPX+zYj1VbbjFdORMgLc3K\nhGXPdJQ+WUUGHvu689+lZBcMHNd+m6sCpLpUlU+xl4NlS990dV8Z1oSWJnWvObuXk4aqhaqh9R/8\nSC1kp92tzF9l+5TALd2rBKkpzu8aSFjnuyiEEAuAO4FaYABwh5TSxRizNpqBB6WUO4QQccB2IcQq\nIBT4o82+twNrpZTfCCFSgMeAxebvfg086ebYXSMhW90YZQeVT8SbtIXxWmsg1gLEfNFJ6Z4PJCIO\nRIh7Jizr2HRQja9uXA5njiutIz5DrdhcJcwEmdNsBMhO5SBstCo0CZZ5RvY124FDVUXgploVaRMS\nqsKnrQXIqSPqxsxbq1a2/Qa7PrfexPHt8O3bkDIK5jzo+XmkVFpHQrbSKvqmq0S+nctg/q/smzgN\n0ibB9a/Ci5fCd1/BmGs67tPmHBbOBcjJgyoXyYgUtCY0XDnWHSXwGdSfUZrU+MXtt0clqPtMSufB\nKoY5Kd0VAWIOLDlzXP3NrEsSOaLND7JeLa5yv4Gpd8Ggeer7Y98oH0pdudJKyg8HrwYC/Ap4WEo5\nD7gGeEMIcZ7zQ9ojpSyRUu4wv68GDgBpUsq9UsrLbX7KrJz2p4EIofhf4BPjPH7DiMQq3uFdBzqo\napwhYe01EHsmrIYqFbnkqgYihPsFFa1LpRgMmgcTFitnZUJ2+8ZZrpA9W5mhjHMf36Fu/OjE9j4Q\nawESGq6EyOHP1LZUc25In4HtTVhlBy3vbZMuNRaMQp6FDkrTuHyeI+qhO+V29SAE9VBrOAN731T+\nrvgsxz6BzOnKVm9tw7emMk+9JmQ7d6IbWuxAOwIElBnrdK7zbHjDgW4dgQVKgLQ0uuDE32Kugj3W\n+X5g+VudKVSv1lnozsiarY7Z9pya04jL1d+231ClYRkO9NQx5u6pQSpApJTnSSnXmd/vBS4Bfu/p\nwEKIbGACsNnJPouEEP8BXgH+CdwPnA9cI4S428Exdwkhtgkhtp08edLeLp5hCJDmeu860EE96MNj\n2l+wRjRFbKpFgBgP4CgXNRBQZix3TFi2iYreINscJZO3Ts2l4jt100Yn2piwrAQIKD9IU63SRJJH\nqm22AuTkQfX94AVqFdzSs+M5PMaIaCva0rU8GiO6yYiIA2XCSR0Lm59W4dm2Tm1rQkLVgsKRf8LQ\nQPqPcv4wLN6pcqESB9n/fvB89erMjGVkyQ+Y0H674VPpzIxVuFVFb7lizu5jFqiGI91W03dEttlE\nt/YxZc4yzGWD5ilt77g5xD7lHFX4Ndic6E4ir0qABc72cXLOWGAF8GMppcM0aSnlSinlD6SU10sp\nv5ZSPi6lnCSlvFtK+W8HxzwtpZwspZycnOxFR3e8lXPY2xoIdOwJYlwI/YZ0FCCuaiBg1kA8MGG5\nI6Q6Y+BE5efJW2cV9TJRjWHPhGUkSRp+kKRhllL1fdLbJxOePKjMVlPvVNFihsaiaY+hgdSdat/b\nxl0OfaKEhbGiBrUAmvYDlcBXV+5cgIAqeng6z3615tP56rqI7e/chFW8U4XGOtKGk4ZB3ADnZqyS\nXcq0FGNzP7kiQJob1Rxc8X+A8utF9nVfgCSPVPNpqFLtG0LMRUAGzVWLq52vqGKvkX2UDyQINZDV\nQoj7hRDtDN9CCBMwQwjxEuByWzohRDhKeCyTUq50a7aBxBRt8X14WwMBc08QqxvGuBCShijzQGNd\nxzpYruBuQUVXL2x3CDOpasF56ywJZQMnqDGsTVjWUVhg8WcY5itQGkhTneV3KjugqpkOuUA9MHa8\n5L159yTKD1u0OE/NWLXlSoMZfmnH70ZfbXnwOsuJAEvekD0tpLJAmS6dtXlublBRd/b8HwZCKDNW\n7jcqMswexbuUELLFFQFSuleF2bsSgWXQN8N9ARISAplmgTzie5bt2bOVf7PquKV+lykm+DQQ4GKg\nBXhNCFEihNgvhMgFjgA3Av8npXzRlcHMmspzwAEp5WMezjlwGGYsXwgQU3R7E1abBmKuFFNT6tnD\n3V0T1tkKFVXTWbMsd8meraJGjn6p/o7RiY5NWIYGYoTyptgIEFBaSFO9Wk33H6mc/eMXK6d7yR7Y\n9G9VUvv9+3Xpk8Y6qCxU1WpNce1L07S2wnv3uRbyeuRzFR47/OKO34VHwYwfqoeZI7OSQfII1RDN\nnh+kMl8Fa5jMbZ7tFQI9sU81Mhs4oeN31gyer4SAUWjQmvozypRq28UTXBMgbW2kXdRAQJmxqgwB\n4oY5evQiZdIzBK8xR2PuqWYfTESsMn378Xp3JZGwXkr5lJRyFpCJMltNkFJmSSnvlFK6U25zFnAz\ncJ4QYpf5x85yJkgx2tv6woQVHmNjwjL7QIz+IdWlVhedi3kg4NyE1dyoHh7WjkYjTNg9q2TnGH6Q\n/HUWx2d0otK0jGJ79WfUAy7UHByYNlGtZodfYjmPYUuuKlaratmqHkigSkHIVtV46tNfqAfAjpft\nF7LrTVR8B0hIHq7+psbDDyD3a9j5KnzxP52f59DHSsuz99AFOPencPe6zq8dIZQZK3dN+4dda6sS\ndIYAARVpZYtR3LMzATJonnq1Z8YqMQsVWwc6uCZACrcoc6qxoHEF62TCunIVTBBm6vy4MdfAvRs7\n+loGzVOvbRpIrLr+m+z8zXyEy050IcRy4CfAeCDck8GklOuklEJKOVZKOd7842LNgSDA0EDcrcTr\nCrYaSGO1yi41wv+qzRqICLHkd7iCYcKytyopP6QeHvus8jfcyXR3h7SJlnwXw/RgrL4MM5ZRxsQg\nsi/c9rF68Bm0aSDHLdns/c2mmYRsFUo6/T5VPv9Hu2D0NapX9uHPvf87dReMIn5Jw5TJpfRby2Jl\nx8vq9fg2i0PWHk31cPQrJcy9sbjIOVfVlbOuNltbpsxC8VlWBUbt+EGKd6prx9ovaY/Y/iqP6eCH\nHc1YRvsBjzWQrZDhhvkKLOG7DTUd62B5wtjrVLh0prmZnGE18KMfxJ14zP8AdcDVwCohxDvmXI7e\nw+AFypZvCBJvYtsXvaFGqaSGsKouVQ/aqAT3wmgj41X4X7OdfgRGYciTVqGw3riw7REarkI4wbJy\nNMYxTHP1Z5Qz0BlxqYBQGkjZARX+nGiV+zH353DxH5TfRAi44gm1QltxB+x+HT7+meoat9y/eagB\npfwIINTfKWOqypsp3gW1p1Ry2vglavW65RnH58hbp0xKwy5xvI87tPlBrMxYhlM9PsuqvI+dh2Hx\nLnUNuSLIpv1ACcaPHrQsog59Cmv+CkMvUlVwbQmPUos3awFSU6ZMon8fA38ZqkJrXcn/sMZYDFYd\ndy8h2BH9R8Ltn1oEntGdNEgFSCbwDXC3lHI88Dbggt7bg8iYojoMhkd5/9ym2PaJhI01altUgrqY\nq0s8u+icFVQ0SqmXHbBs85UAAXXDRvSxOC6N36XOgQZij9Bw5YMyNJB+Q5ybAUzRcMNyFa//zg+U\nxiVbVbRWb6nqW35YJYCaoiHNXJqmaAvseUMtLmbcC+NuhG9XWPITbDn0kVrk5JzrnTklZCtTlbXv\nxQjhTchy0uKgThX07Mx8ZTDpVpj1Y1Up+ps/K9PTW0vVouKa5+wfI0THbPTCLUrzSRmtmkXNvF81\niHMH61wQbwgQWyLMAsSPjnSXM9GBYcC1wCghRDWwB5grhPgI2COl9GLSRS/EZE8DiVMXc1yK0hbq\nKtwPr7XuCdJnQPvvDA2k/LCyP4eE+ObCNph6p1K7DVU7yo4G4opNuc9ApYFU5NqPorElPgNu/0wJ\nnczpakX6wiWqLP2Iyzz7XboT1q2UY/opTaRwq0oITJukcgim3qlKhW9/UfkyrDn2Nex4RfV3McKp\nvUHOXDjwvjIvhYRakgj7ZlgS7mxNWCe+VRqUqwIE4PzfqGv96z+oTol9Bqi2s84CRWwFiOEnvPJJ\nzxdYbQKkSN3L/c/x7DyOCHIN5HEp5RVSyiHAecCrQDSqvEj38WMEK+G2UVhVlgs8boBZA/HAP9FW\nUNGOI90oDNlUp6JfWlvVTeMrARIS2v7mM8Zx5ANxRJ+BKjHudJ7F/9EZSUOU2SQswlyfK7J9eZWe\nSmur+lsZ0XygzFhHV6m8jQk3q23Jw5VTdtvz7Qv+lX4Lb9ysgjkutam621Vy5qr/ueHkrixQ0Vmm\naMvDsMlGgBgOdGchvLYYpsxhFysT6ZKV9k1X1tj2BKnMV9qzOwEstsQNUD7MM8eVpudtTd94XvhR\nA3FHgHwuhMgTQnwCPIIqpvialPJ2KaWb3iRNB0zmTHQjbNEwYYGy+xs+kGg3L+BIwyFox4RVW2Z5\nX3bA7Gxv9Z0AscWuD8QFAdI33WzukJYILHcIi1APUcOR2pOpLlbXVZKVAEmfrExX4dEqf8Ng6l1K\nS/vyN5C7VoXLLrtWXYeL37Jos95ixGVKG/rwx+p/fzrfUuXZkRP9+A5lwoyz0aY7IzQcbnwdHtht\niaZ0hj0NxOjT4SmhYWre5YdVdFlnZUzcpe1v5r+eIO6UMpmA6kz4U2ANqrRIF6qyadphRCgZWojh\nRAezBlLqmXnJWVvbmhMqvhzMWcQeZLp3BaN9aN1p5eBsqHKtVa+1mctVDcSWrNlqde1O7+vuiJGB\nbpiwwOL8PWdR+6CFYRerIJENT8BLl8O/Zqpw8sVvtc889xamaLjqP0poffqQWhQYAqStQrUdDcRV\nB7otQrheRduuAMl2f0xb+qZbanB5+z4zBaEPxDYDHagGtpnfD7SqYlLprCyJphPaKvLWmROCrDSQ\n2BRLlra7PpDOTFj9hqhKn2UHLf0HfOVEt4eRjd5Yo7Qfl0xY5lyQkPDOk9YckT0bkKp204juk4qE\nlOoaqS1X14dtGQ5brEN4DVJGq3BnWydwSKjyFVWXqqTPk4eU09y6EoC3yZiiqgMbTanOMXe0bLPn\nWwkQKZU5zh9+K+umUq2tSjsadlHXz2tULwbfOdH96ANxxYn+EqovhzORL4EXgZe9MKfeie2Ky3Ci\nQ3t13atRWGUqUa//CLMG4oMyJp0RnaDGtS2k6AxDA0kaqkwTnmD4QfLXdw8B0nQWPv+1KhhpJNeF\nRaqoQGeBBOWHVcKadfJrSIgKd7aHEMrJ3GeA6uznD879uYqKK91jZcKyo4E0VJt71fhhgROVoAR1\nU71a4LQ0eE8DMfC6BuJ/H0inAkRKOd8fE+n1WPdFb21Vdsw2AWKVuOjuzRMarlZztias5kZ1Y8Sm\nqLDk3LUWn4g3Cyl2RpS5nIknAsQT/4dBeKRKqusOfpAT++Ht25WQH79YaRNRCbD6D7DiTrjra8v1\nY0v5ERVA4O3KAt4kzASLnlG/o9FHxNakCxaNwJ1EWk8xnOX1lZYILK8IkAzLe28LkNAwtagIMg1E\n4w/aVPY6S+SJtRPdwJOLzl45k1pz1HVssuov0NKgOrx5OoanRPdTN2hbN0IXBEjcALWqNtrtekr2\nbPj6T0o787aD2Fsc/gzevEUtJpashCELLN/FZ8ArV8EXjziOkCo/0r6GUrDSfwTcu8HyOSTUnFxr\n9TA0FkFdiYRyFetsdCPBMcEF53tn9LHqk+KL+8zPPUHc6Uj4EzubzwDb3ayHpbFH24qr1qKCRtgR\nIJ5oB/YKKhraRmyKpcpw/gaVtGj4Y/yB0VSqTQNxwYkeFgH3b+u6pmT4QQo22S8QGAxseVqFtn7/\nC5UPZM3g82D6vbDpKRh6IQy9oP33DdUqCss6Aqs7YVuh2riG/SHs2wmQPEB4J5DAOIe7JYlcxc89\nQdwJ450M3A2kmX/uAuYBzwghHBhUNS7TZvOts6wgDJtmZLxSTcF7GkiNlQAxak1VfOebQorOiO6n\n5mb4X1y9qWL7W4ouekraZCUwg9WMJaWKOso5t6PwMFjwiIqke/fejv9jo1tdv24qQGxLugfChGUI\nkD5prkdwOcMQIFGJ7nf2dAU/9wRx5zfoB0yUUj4opXwQJVCSgXOBpT6YW+8i3CoKy6jEa2ggQlhq\nQHmy+rLXE8TIQo9JVuMYzkt/mq/AokUYdmZXTFjeos0PEqQJhUbJC3sVYw3CI2HhU0qj3PSv9t9t\nekqZ+rxVfsTfmGKCxISV5736d1EJ6l731X1mimnfDtvHuFsLq9HqcxOQJaU8CzR4dVa9EeuokzYB\nYlVqITZVCQKjI5k72DNhtWkg5ugco9mQP0N4rcczuuS5kgfiTQaOt+RKBBtG863Osq4HTlC9sjc+\naVmln9gHBz6A6XcHr3+nM4zkWoNAmrC8JUCE2RTmKwES4V8fiDsCZDmwSQjxiBDiN8B6VJOpGGC/\nLybXq7AO420zYcVavk8c5LkN1pEJK6KvpTBkf3NEk781EEOAVOSqv4Er/RG8iRGu2RyEa6DinSrX\nJcWFPIx5D6lcoY1Pqc/f/FmZM6bd7ds5+hKTTY+cs6fV3yPcQcSZN4mIAxEKVSWqmZs3K3DP/TnM\nuM9757PG5F8fiMtGZCnl74QQHwOzUTkhd0spjYTCxb6YXK+iLZHQ2olupYFc9Gj71Zg7RMWrsOCW\nZovfoLasfT2gQGkgUVYaiL+1D2i/0vRFn5euULwTUka5ZntPHQ2jrlRmrMHzYf97MOcn/v9/epPw\nmPbVgevN0XL+8NEZFXmNrHFvCpAx13jvXLYEsQYC0Ay0ml+bvD+dXkyoSa14GusstWysNZDoxK5p\nINBeC6kpa9+aN2AaiFHS/ZR//R8GrjQPCgStrea+F24UDZz7S/XwWHadWqXP+KHv5ucPOmgglf7x\nfxhEJUCJOcDUFz2AfIEpzn4TLh/hTkfCB4BlQBLQH3hVCHG/rybW6xDC3BOkrmMYb1cxolasHek1\nJ9pnJyePgKTh7j2wvDoaiGoAABnASURBVIH1ClkLEAunc6HhjHtly1NGqZLrjdWqPHt31j7A3OLA\nxoTljwgsg6gEy2o+oZPuh8GCoYH4qS+6O3GQdwDTpJS1AEKI/wU2Ak/4YmJCiHnA74B9wOvmzW2f\npZRf+2LcgGLcMA3VKk7cW7Ze4yFZVwH9zN37ak6q/AKD8Cj44ZaOx/qa8GgVStvSEBgBYjxknQmQ\nsoOqx/p1r7iWp+INXO37bcuC/wdImPkjr0/J7xgLKoP6yvZas68x7pvwaBWt2B1o64te55d8LndM\nWAKwbizcgvP6WO0PFiJDCLFaCHFACLHPrNE4QwI1QCRQZOdzz8PoCdJYo1RRb9l6jYq1hjredFat\nbq01kEAhhMWMFawayMEPVVMlZz3DvU3xTpX742614cQcuPbFzossdgcME5bR4iAQJixQ5qtgLgVj\njSE0/ORId0cDeQHYLIR4ByU4FgLPu3F8M/CglHKHuZf6diHEKiAU+KPNvrcDa6WU3wghUoDHgJtt\nPvc8x73RlTCsxnvmK1A5HrGpqi3n1DvbJxEGA9GJKmPaX6t7a1wRICe+Va8nDykHtT84vkO1XfW0\nWGRPIDwakKp4pClGaSD+NmFB9/F/gCXwprEG8P397U4U1mNCiK+BWSgBcqs7JUyklCVAifl9tRDi\nAJAmpVwFXO7k0NNAhJSy1fqzo52FEHehsuTJzLStRB/khMeoKKxGU3sHelcRQjVQMspIt9XBCgIN\nBCxmpEBoIKZYCAlzLkCMjO6TB/0zp9YWFf0zoeetkdzCui96WKQKAvFnTkt3FCBtPUH8k0zoSj+Q\napT5qG2T1XdSSun2slEIkQ1MADY72WcRcBEQD/zT9rOj46SUTwNPA0yePNk/niRvYYpWRQVDwr2r\ngYBqFHTgfdXrwchCDxYBEhVAAWKEazoSII21cOo79d5fAqT8iFpI+DugIdiw7vFtJND604RlLGy6\nkwDxc08QV8q5O+k87z5CiFhgBfBjZw2opJQrgZU2m20/9yzCo1XikghpnwPiDTKmqdfCLVBnjq0P\nJhMWBEaAgHMBUnYAkMoEWHZARbf42h5ebM5Ad9eB3tOwbnFgrGG1Ccs5fu4J4oNqXo4RQoSjhMcy\ns4DQWGMyTFg13jVhAQwYq6KdCjerCCwInsgSw4keiERCcC5ASveo19GLlA3e8B/5kuKdypzZXavo\negvrvujG/8efJqyBE8w/3UgT9LMG4jcBIlTv2+eAA1LKx/w1brfCqD5q3Y3QW4RFqJuhcIsyYUUl\nBo+Dts2EFaCaTVEJln7wtpR+qwSbUSrdH2asimNKeHhS96wnYRQYbayx1MHy5zXSb7Bq1hUbJAst\nVzD1UAGCcr7fDJwnhNhl/ukGvUT9SFsYb7X3NRBQjvSSXarKa7CYryCwYbxg1kDstPwFFYGVMlqV\nTAcVieVrasuDRzsMJNZOdH9W4u3OGBpITzNhSSnXSSmFlHKslHK8+edjf43fLTCqjzZUe9+JDsoP\n0tKo2tcG06oqezaMWuh+zoO3cGTCam1VGkjqaCVwI/uqtrK+pu6U/0vKBCOBNmF1R6y1Nj+gW9oG\nE0bmeWuz901YoDQQUHH1waSB9E2D614K3PhRCeZik03tzXqnc5VPKnWMcpwnj/CPBlJ3CmKSfD9O\nsGNdYDQQJqzuSGgYhEX5LYzXr050TSdYlx4w+UCAxPa39HUOJgESaNqSCW3MWEb+h1FOPXmEJRLL\nVzTWKS1UayDtNZD6SvVgDI8M7Jy6A36syKsFSDBhXfvKFyYssITzahu7BUfZ6Ce+VRWSDdNa8gjV\nv926xLi3MVr7agHSvkfO2dPafOUqfuwJogVIMGGyEiC+cKKDxYylNRALxoPJVoCUfquioYymW0bv\neG/5QeoqOt7oRo6ONmGpKLSwKLMA8XMZk+6M1kB6KeFWJixfaSCDz1Or2wFjfXP+7ogjDaR0b/tu\ngIYm4i0/yLJr4JNftN9Wa2ggWoAAloKK/i5j0p0xxQVlMUWNr7H2gfgqqS4xB35+zDfn7q7YEyB1\nFVBVpBzoBnED1P/FG7kgLc1QskeZyKzRJqz2mMyh7WcrIT4j0LPpHkTE+ifhFa2BBBf+MGFpOhJl\npyfIiX3qNdVKAzEiscq8IEAq86G1SfXbtqbNhKUFCKDug8Ya/zeT6s6YtAmrd+IPE5amIxF9VP0x\newLE2oQFyg/iDQ2k/LB6rS5tH9VVW660Ev2wVIRHW6KwtAnLNUwx2oneK9EaSGAICVEP7LNW5UxO\nHVWCxTbYIHmE0hK6GollCJCWRhvTmTmJsLs0MPI1phjl/2is0VnorhIRpzWQXkm4FiABwzYbveI7\nSBzU8UGeNEy9GiXePcUQIKC0EAOdRNgeUyxUFav3WitzDcOEZXRy9CFagAQThhM9PFpllGr8h60A\nOfWdpX+8NYZzu95B7SxXKT+iQlQBqkss22vLtQPdGlO0RcBqE5ZrGObvplqfD6UFSDARalL2b619\n+B9rAdLcoApO9hvScT+j4KOj4ouuUn4YMqer97YaiBYgFkwxtPUC0SYs1zD5r6CiFiDBhBDqn68d\n6P7HWoCczgPZCol2NBBDgNSf8Xys2lNqrJw56rN1JFZduTZhWWMdWKJNWK7Rri+6b9ECJNgwRWsN\nJBBYCxDDv2HPhOUNAWL4P1LHqfMZGkhLs5qDTiK0YJ0bpU1YruHHniBagAQb4dG+qcSrcU5UghIK\nrS3KgQ7KiW5LmEn9j+z5QOoqXHNcGgIkaYhKTjR8IIYA0yYsCyatgbiNH3uCaAESbMT2Vz8a/2LY\n1+vPqBDeqERLr3ZbIvt2FCC15fC3EfDUdNj9uioN74jywxAWCX0zVJhw9Qm1XScRdkRrIO5j8l9P\nEC1Ago2rn4VL/hzoWfQ+rMuZOIrAMoiM72jCOlMELQ1KCLzzA3hioiUZ0ZbyI8pBHxJq1kDMJiwj\nt0SbsCwYD0NTbPC0YA524gbAhJvVq48JWgEihJgnhFgrhPi3+X2mEOJ9IcTzQohfBnp+PqNvutZA\nAoG1AKk4Zt+BbhDZt6MAMZIQr38VbnxDOcq3Pmf/+PLDqsovQFyKcqJLqetg2cPIjdLmK9fpMxCu\n/CcMHO/zofwmQIQQGUKI1UKIA0KIfUKIBzo5RAI1QCRQBAwDPpJS3g6M8vF0Nb0NQ4CcKYKq451o\nIPYEiNFyNRGGXwxZMyFvXcdjm+pVHSwjITFugCUbXZdy74jhENYhvEGJP7PVmoEHpZQ7hBBxwHYh\nxCogFPijzb63A2ullN8IIVKAx4AfAb8SQlwPvOLHeWt6A8YD6vh29dqZACm3KeleZ9ZADL9J9iz4\nYpWqimqtUVYcUyHCbQIkVb1Wl1idQ2sgbRgmLO3/CEr8JkCklCVAifl9tRDiAJAmpVwFXO7k0NNA\nBHAb8IiUco0Q4m3gBV/PWdOLsBUgzkxYUfEdEwlte3Znm3M88tfDOVdZ9muLwDKbsGKtBEhtOUT0\n1bZ+a4z6cEb4tCaoCEi9DCFENjAB2Oxkn0XARUA88E+gHPiNEOImIM/JcXcBdwFkZmZ6a8qano7x\ngCreqV4700AaqlTIbojZCny2QjXyCTOpzwPGKfNL3jobAXLEfH5zlnubBlJqTiLU2kc72kxYWgMJ\nRvwuQIQQscAK4MdSyipH+0kpVwIrbTZf09n5pZRPA08DTJ48WXayu0ajCA2z+DZi+jvPxYnsq8xQ\njTUQaW78VVcB0VZ2+tBw1X8+b337Y8sPQ590i2mmnQA5pSOwbDGc6NoHEpT4NQpLCBGOEh7LzAJC\nowkejIeUvRpY1tjLRj97uuNDLnu26p9uXfr91BGL+QpUv3UjG71W18HqgCHI9d8lKPFnFJYAngMO\nSCkf89e4Go3LtAkQOxno1hh+DutkwrMVls6GBtZ+EFA91kt2Q9rE9vsZ2ejahNWRiFi4fhlMuCXQ\nM9HYwZ8ayCzgZuA8IcQu88+lfhxfo3GOIUCcOdDBdQ1k4HhVDDBvncrz+Oinap+Z97ffLy7VyoSl\nBUgHRl6uBWuQ4s8orHWAbrOmCV7aNBAPBEhdRcfSJ6HhkDlNCZDdr0PhJrjiiY6CJjYVij9V+SDa\nB6LpRgRtJrpG43c89YG0tpp7dttx9GbPhrL98PmvIG0yjF/ScZ+4VIs5TCcRaroRuu2dRmMQmwIh\n4ZCQ43w/I6TUyP1oOKOismx9IABZs9VrXQUsWWEJ+7XGumaRNmFpuhFagGg0BlPvhEHzLclrjogw\nh+4aGohtFro1Aycos9ToReq9PeJSLO+1CUvTjeh1AqSpqYmioiLq6+sDPZWgITIykvT0dMLDe3kG\ndFQCZEzpfL+QUCVEDAFiaCL2TFhhJvjRDpVk6AhrDUQ7izXdiF4nQIqKioiLiyM7OxsVWdy7kVJy\n6tQpioqKyMnpxHSjsWBdUNGoxGvPhGXs64xYaw1ECxBN96HXOdHr6+vp16+fFh5mhBD069dPa2Tu\nYt0TxDBheZotbWSjh0bodsaabkWv00AALTxs0H8PD7DuStjWitaBBtIZ4VFKIJliQP8vNN2IXqeB\ndCdefPFFiouL2z5nZ2dTXl7u5AiN3+hgwhJdqxgbl+q5ANJoAoQWIEGMrQDRBBHtBMhp9Tkk1PPz\nDb0QBi/wztw0Gj+hBUgAePnllxk7dizjxo3jqquuIicnh6amJgCqqqrIzs7mrbfeYtu2bSxevJjx\n48dz9uxZAJ544gkmTpzImDFjOHjwIAAVFRUsXLiQsWPHMn36dPbs2QPAb37zG26//XbmzZvHoEGD\nePzxxwPzC/dEomx8IF3VHi78HVzwP12fl0bjR3qlD8Tgfz7Yx/5ihxXlPWLUwD488r1zHH6/b98+\nHn30UdavX09SUhIVFRU8+OCDfPTRRyxcuJDXX3+dq6++mmuvvZYnn3ySv/71r0yePLnt+KSkJHbs\n2MFTTz3FX//6V5599lkeeeQRJkyYwLvvvstXX33FLbfcwq5duwA4ePAgq1evprq6muHDh3PPPffo\ncF1v0NYT5P+3d++xVdZ3HMffX7DSOi5O64WLm2yByQBrBToHiWA1oNZYo/whwbK6EcMWi8ZLxBm3\nbjNxybJpxhIdRggjDErQrJ2Aix3WS1yFYlQ03SZz2DWg4JGVjstcyXd/nKe1nN4OB/o85/T5vJKG\nc3p+z+n321853/P7/Z7ze04EGylqu3GJH41AQrZ9+3YWLlxIYWHyA2PnnXceS5cuZc2a5AUW16xZ\nw5133tnn8bfeeisAM2bMYO/evQC8/vrrVFRUAFBaWkoikaCtLfnuuKysjBEjRlBYWMiFF17IJ598\nMlipxUv37UyOHer7FF6RISzWI5D+RgqDxd17nPU0Z84c9u7dyyuvvMKJEyeYNm1an8ePGDECgOHD\nh9PR0dH1nKk6f0Zn+9Rj5DR1LyBHP4PzJ/XfXmQI0ggkZNdeey2bNm0ikUgAyfULgCVLlrBo0aKT\nRh+jRo2ivb19wOe8+uqrWb9+PQANDQ0UFhYyevToQYheunRdE6Qt+Ul0nUElMaQCErKpU6fyyCOP\nMHfuXIqKirjvvvsAWLx4MYcOHWLRokVdbSsrK1m2bNlJi+i9qa6upqmpicsvv5wVK1awdu3aQc8j\n9jpHIEc/TW6mqDUQiSHrbfpjqJg5c6Y3NTWd9L3m5mamTJkSUUR927x5M7W1taxbty6Sn5+tv5es\n9fF78PQcKPslbLkfbvgFfOuuqKMSOSPMbJe7zxyoXazXQLJFVVUV27ZtY+vWrVGHIunqHIEc2pv8\nV1NYEkMqIFlg5cqVUYcgp6rzmiCf/fPk+yIxkrUFxMyGAT8DRgNNwGvAI8AYd18YZWwinD0SbBgc\n+ih5X6fxSgyFtohuZpeY2ctm1mxm75vZPQMcUg6MB/4HtLr7h+7+vcGPVCQNFux9dahzBKJFdImf\nMEcgHcD97v6WmY0CdpnZS8Bw4PGUtt8FvgH8xd1/a2abgT+HGKvIwPLHaA1EYi20AuLu+4H9we12\nM2sGxrv7S8BNqe3NrBX4PLh7Iqw4RdLW+VkQG/7FZW5FYiSSz4GY2aVAMfBmP82eBxaY2UrgVTM7\n38yeBorN7OF+nvsuM2sys6aDBw+eybBD19TUxPLly6MOQ/rSeSZWwZd1HQ+JpdAX0c1sJPAccK+7\n97mTobsfBVLXPJYN9PzuvgpYBcnPgZxGqJGbOXPmSRspSpbpLCCavpKYCnUEYmZ5JIvHend/Psyf\nnU2OHDlCWVkZRUVFTJs2jZqaGnbu3Mns2bMpKiqipKSE9vZ2GhoauOmm5OxedXU1FRUVlJaWMmnS\nJJ555hkAKioqqK2t7XruxYsXU1dXF0lesdN9BCISQ6GNQCy5u9+zQLO7/yqsn9uvbSvg491n9jkv\nng43/LzfJi+++CLjxo1jy5YtALS1tVFcXExNTQ2zZs3i8OHDFBQU9Dju3XffpbGxkSNHjlBcXExZ\nWRlLly7liSeeoLy8nLa2Nt544w1tZRKWzs9+6BReiakwRyBzgAqg1MzeDr5uDPHnZ43p06dTX1/P\nQw89xGuvvUZLSwtjx45l1qxZAIwePZqzzupZ28vLyykoKKCwsJBrrrmGHTt2MHfuXPbs2cOBAwfY\nsGEDt912W6/HyiDQCERiLsyzsF4HsmulcYCRwmCZPHkyu3btYuvWrTz88MPMnz+/xxbvvUlt03m/\noqKC9evXs3HjRlavXj0oMUsvOs/C0hqIxJR2443Avn37OOecc7jjjjt44IEHaGxsZN++fezcuROA\n9vb2Xq/bUVtby/Hjx0kkEjQ0NHSNWCorK3nyySeB5G6/EpKuEYi2MZF40lxHBHbv3s2DDz7IsGHD\nyMvL46mnnsLdqaqq4tixYxQUFFBfX9/juJKSEsrKymhpaeHRRx9l3LhxAFx00UVMmTKFW265JexU\n4i1fayASbyogEViwYAELFizo8f3GxsaT7s+bN4958+Z13Z88eTKrVq3qcdzRo0f54IMPTrqWiIRA\nayASc5rCynH19fVcdtllVFVVMWbMmKjDiZdxV8DsKvh6adSRiERCI5AcUV1d3ev3r7vuOlpaWsIN\nRpLOGgHzH4s6CpHIaAQiIiIZiWUBGcqX8c2Efh8ikonYFZD8/HwSiYReNAPuTiKRID8/P+pQRCTH\nxG4NZMKECbS2tpLrO/WeSfn5+UyYMCHqMEQkx8SugOTl5TFx4sSowxARyXmxm8ISEZEzQwVEREQy\nogIiIiIZsaF8NpKZHQQ+yvDwQuDTMxhOlJRL9hkqeYByyVank8tX3f2CgRoN6QJyOsysyd2HxPVk\nlUv2GSp5gHLJVmHkoiksERHJiAqIiIhkRAWkbz33Tc9dyiX7DJU8QLlkq0HPRWsgIiKSEY1AREQk\nI7EvIGZ2vZn9zcz2mNmKXh4fYWY1weNvmtml4UeZnjRyqTSzg2b2dvC1NIo4B2Jmq83sgJm918fj\nZma/DvJ818yuDDvGdKSRxzwza+vWHz8KO8Z0mdklZvaymTWb2ftmdk8vbXKlX9LJJev7xszyzWyH\nmb0T5PGTXtoM7uuXu8f2CxgO/AP4GnA28A7wzZQ2PwCeDm7fDtREHfdp5FIJ/CbqWNPI5WrgSuC9\nPh6/EdgGGHAV8GbUMWeYxzzghajjTDOXscCVwe1RwN97+fvKlX5JJ5es75vg9zwyuJ0HvAlcldJm\nUF+/4j4CKQH2uPuH7v45sBEoT2lTDqwNbm8GrjUzCzHGdKWTS05w91eBz/ppUg78zpMagXPNbGw4\n0aUvjTxyhrvvd/e3gtvtQDMwPqVZrvRLOrlkveD3/J/gbl7wlbqoPaivX3EvIOOBf3W730rPP6Su\nNu7eAbQB54cS3alJJxeA24Lphc1mdkk4oZ1x6eaaC74dTEFsM7OpUQeTjmAapJjkO97ucq5f+skF\ncqBvzGy4mb0NHABecvc++2QwXr/iXkB6q8SpFTydNtkgnTj/CFzq7pcD9XzxziTX5EqfDOQtkltG\nFAErgT9EHM+AzGwk8Bxwr7sfTn24l0Oytl8GyCUn+sbdT7j7FcAEoMTMpqU0GdQ+iXsBaQW6vwuf\nAOzrq42ZnQWMITunJQbMxd0T7v7f4O4zwIyQYjvT0um3rOfuhzunINx9K5BnZoURh9UnM8sj+YK7\n3t2f76VJzvTLQLnkWt+4+7+BBuD6lIcG9fUr7gVkJzDJzCaa2dkkF5nqUtrUAd8Jbi8EtnuwIpVl\nBswlZT76ZpJzv7moDlgSnPVzFdDm7vujDupUmdnFnfPRZlZC8v9jItqoehfE+SzQ7O6/6qNZTvRL\nOrnkQt+Y2QVmdm5wuwC4DvhrSrNBff2K3RUJu3P3DjO7G/gTybOYVrv7+2b2U6DJ3etI/qGtM7M9\nJCv37dFF3Lc0c1luZjcDHSRzqYws4H6Y2QaSZ8EUmlkr8GOSC4S4+9PAVpJn/OwBjgJ3RhNp/9LI\nYyHwfTPrAI4Bt2fpmxOAOUAFsDuYcwf4IfAVyK1+Ib1ccqFvxgJrzWw4yQK3yd1fCPP1S59EFxGR\njMR9CktERDKkAiIiIhlRARERkYyogIiISEZUQEREJCMqICIikhEVEBERyYgKiEiIzGxWsJllvpl9\nKbiOQ+r+RSI5QR8kFAmZmT0G5AMFQKu7Px5xSCIZUQERCVmwV9lO4Dgw291PRBySSEY0hSUSvvOA\nkSSvhpcfcSwiGdMIRCRkZlZH8oqRE4Gx7n53xCGJZCTWu/GKhM3MlgAd7v77YBfVN8ys1N23Rx2b\nyKnSCERERDKiNRAREcmICoiIiGREBURERDKiAiIiIhlRARERkYyogIiISEZUQEREJCMqICIikpH/\nAydeHxO6pyqYAAAAAElFTkSuQmCC\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x7f12c4ddfad0>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"import mpmath; mpmath.prec = 128\n",
"def err(f, x): return float(abs(mpmath.digamma(x) - mpmath.mpf(f(x))))\n",
"pl.plot(xs, [err(fast_digamma, x) for x in xs], label='cython')\n",
"pl.plot(xs, [err(scipy.special.digamma, x) for x in xs], label='scipy')\n",
"\n",
"print np.mean(np.log2(np.array([err(fast_digamma, x) for x in xs]))\n",
" - np.log2([err(scipy.special.digamma, x) for x in xs]))\n",
"\n",
"pl.xlabel('x')\n",
"pl.ylabel('$\\log | \\psi(x) - f(x) |$')\n",
"pl.yscale('log', basey=2)\n",
"pl.legend(loc='best');"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Looks like the scipy version is more accurate by about 9 bits than our Cython version. The Cython version is still accurate to over 40 bits, which on a 64-bit float is pretty darn good."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python [default]",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.14"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment