fredrik-johansson/Arbdemo.ipynb

## Arbdemo.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Arb\n",
    "\n",
    "A C library for arbitrary-precision ball arithmetic\n",
    "\n",
    "http://arblib.org/\n",
    "\n",
    "### The painful way (C)\n",
    "\n",
    "    #include \"arb.h\"\n",
    "\n",
    "    int main()\n",
    "    {\n",
    "        arb_t x;\n",
    "        arb_init(x);\n",
    "        arb_const_pi(x, 333);\n",
    "        arb_printn(x, 100, 0);\n",
    "        printf(\"\\n\");\n",
    "        arb_clear(x);\n",
    "    }\n",
    "\n",
    "### Wrappers\n",
    "\n",
    "* SageMath (RealBallField, ComplexBallField)\n",
    "* Nemo.jl (ArbField, AcbField)\n",
    "* Python-FLINT\n",
    "* ... and others"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Python-FLINT\n",
    "\n",
    "Very quick installation:\n",
    "\n",
    "    pip install flint-py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Usage:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.00000000000000 +/- 3.61e-16]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from flint import *\n",
    "\n",
    "arb(\"0.1\") * 10"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068 +/- 3.51e-101]\n",
      "[+/- 1.72e-101]\n",
      "[-1.000000000e-90 +/- 6.09e-101]\n"
     ]
    }
   ],
   "source": [
    "ctx.dps = 100\n",
    "x = arb.pi()\n",
    "print(x)\n",
    "print(x.sin())\n",
    "print((x + arb(\"1e-90\")).sin())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEKCAYAAAD5MJl4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJztnXd4FNX6x78nPSS0QIBQQgcJIk0QVJpeBVFRQcWrF8UG1qvXjtgVfooNC6jXhpUrAioqIKKEIj10iEgNhAQSIIQkpO/5/fFmPLObLTO7Mzs7m/N5nn0yuzsz+2Z29rznvJVxziGRSCQSiZoIqwWQSCQSSeghlYNEIpFIaiGVg0QikUhqIZWDRCKRSGohlYNEIpFIaiGVg0QikUhqIZWDRCKRSGohlYNEIpFIaiGVg0QikUhqEWW1AHpo1KgR79Spk9Vi+KSkpAQJCQlWi+ETKadx2EFGQMppNHaRMyMj4zjnPFnPMbZSDs2bN8fGjRutFsMn6enpGDp0qNVi+ETKaRx2kBGQchqNXeRkjGXpPUaalSQSiURSC6kcJBKJRFILy5QDYyyOMbaeMbaVMbaTMfa8VbJIJBKJxBkrfQ7lAC7inBczxqIBrGKMLeKcr7VQJolEIpHAQuXAqZFEcc3T6JqHbC4hkUgkIYClPgfGWCRjbAuAPAC/cs7XWSmPRCKRSAgWCp3gGGONAHwH4H7O+Q6X9yYAmAAAycnJfefMmWOBhPooLi5GYmKi1WL4RMppHHaQEZByGo1d5Bw2bFgG5/xcXQdxzkPiAeAZAI9426dLly7cDixbtsxqETQh5TQOO8jIuZTTaOwiJ4CNXOeYbGW0UnLNigGMsXgAlwD40yp5JBKJxB+WLQP+DMORy0qfQwqAZYyxbQA2gHwOP1koj0QikeiCc2DCBKBbN2DQIGD/fqslMg4ro5W2Aeht1edLJBJJoGzd2hB799L29u1AixbWymMktqqtJJFIQouTJ4E776QZc/36wPjxwG23WS1V8Fi4MOXv7auvBpYvB5Tyb08/bZFQBiGVg0Qi8ZtnngHmzxfPV60i80rnztbJFCw4B1avbvr380svBUaOpO1mzYCnngIYs0g4A5C1lSQSiV+UlQFffeX8GufAF19YI0+wOXIEKCmh+XWjRsDYsYAS1ZqXB2RnWyicAUjlIJFI/OL774FTp2q//sUXgMMRfHmCzc6dYrt7dyAyEujbV7xmg+4CXpHKQSKR+MWAAcDkyUCrVvS3cWN6/eBBMi+FO7t2ie20NPrbr594TSoHiURSJ2nXDnjpJSArC3jySeCGG8R7n39umVhBw3XlAADnqnKQN2wIrjxGI5WDRCIJiMhIoF494F//Eq8tXWqdPMHC3cpBrRw2biQfjF2R0UoSicQQ+vUDunalWXTv3uR3iAjT6SfnzspBWTl06EDO6VOngIICckq3aWONjIEilYNEItEN57XDNKOjw7OMhDvKyoBrrgHWrj2NkpIGSKlJd2CMwngVk1JWln2VQ5jqdYlEYhacAx07AkOGAA89BJw5Y7VEwSc+Hvj0U+C99zYhK8tZUbZtK7YPHgy6aIYhlYNEItHFoUPAgQPAihXAJ5/QQFmXcV1BtWsntrOygiqKoUizkkQi0UVGhtju08feWcBmMHQomZ3atqVtuyKVg0Qi0YVaOaiTvgAyOb39Nu2zbRuwZk3dW1lcfjk97I5UDhKJRBebNontPn2c32MMmDED2LOHnu/c6RzeGS5MnkyRWGfOpKBPH6BBA6slMh6pHCQSiS62bhXbrsoBoDBWRTls3hyeymHGDKCwEAC64oknrJbGHKRDWiKRaKa4GMjNpe3oaIpacuWcc8R2OIa2FhUpigGIialG06be97crUjlIJBLN7Nsnttu3B6Lc2B7UCkO9f7hw+LDYTk4ud+uQf/llKt+dlmbfMhpSOUgkEs0oXc8AoFMn9/uolUM4tc1UcFUO7li3Dli0CMjMtK+ClMpBIpFoRoty6NBBbO/fb+/6Qu7QohzUuQ52TYSTykEikWhGrRw8dXtLSgIaNqTtkhJqfBNOqJVDs2bulYM6S9quiXBSOUgkEs1Mn075C/PmAZdd5n4fxsLb76BFOciVg0RiU44fB/79b+DOO4GPP6YoHIlvEhKAHj2A0aPdRyopqE1L4awckpPL3O6jLrZ35IjJApmEzHOQ1Dlyc4F//EOUXP7oI4ouWbsWaNLEWtnChXB2SmtZOShVWgHg6FGTBTIJqRwkdYqSEmDYMGD3bufX9+4F7r4b+OYbWSvICEaNogGyY0egZ0+rpTEOzrUph2bN6D7inFaplZWUF2InpFlJUqeIjwemTqUs3shI4I47xHvffgvMnm2dbKHO8eNUkdXh8L3v+ecDDzwAXHGFffsZuKO6GpgyhUqV33QTkJhY5Xa/qCggOZm2ObenU14qB0mdIiKC7OUZGRSD/uGH5HdQePbZ8Au9NIpZsygKp1494MUXrZbGGqKigAcfBF5/HfjyS++rzBYtxLYdTUtSOUjqJErHLoB+6Ero5d69wKpV1skVyii+g/JyIDHRWlnsgFo5KCVH7IRUDpI6T/36wA03iOeffGKdLKHMoUNiWx2q6QvOyRxT13j4YfJhLV8OXHCB1dLoRyoHSZ3gwAGKRvJkL7/1VrG9fr02u3pdQ60cUlN97/+f/5BDOi4OWLrUPLlClUsvBa6/Hhg8GGjc2Gpp9COVg6RO8MEHwMCBZDP/+uva7/fvT07GX38Ftm8n34TEGb3K4eRJMkVVVDhH+NiZN96gQX/cOGDZMqulMRcZyioJezgH5s6l7exs941ZGCPfg8Q9hYWiTHVcHDSVqW7dWmxnZ5sjV7DJyKAJBEC5MuoyGeGGnB9Jwp6tW0WWbv36wCWXWCuPHVHP/Nu00ZYLog5hDZeVgzoktXlz7cdVVhovi9lI5SAJe+bPF9tXXgnExloni13Ra1ICwnPlcOyY2G7WzPu+hw8D3boBjRoBZ51lrlxmIM1KkrBH7Qy9+mrf+xcUAL//TuYnucogAlUO4bJyUCuH5s2B06c979uggeiEV1lJ5k07Zd9btnJgjLVhjC1jjO1ijO1kjD1glSyS8KWoiKKPFC66yPv+339P9ZWuvRZ47TVzZbMT/igHtVkpHFYO1dWUJa6gZEB7okEDysgHgDNn6F60E1aalaoAPMw5TwMwAMC9jLE0C+WRhCErV4oY+549fRfW69NHZEivWAGUuS+6WeeoqBADnVblkJREzmvAue+yXTlxQoQ4N24MxMR4358xe2dJW6YcOOe5nPNNNdtFADIBtLJKHkl4og439LVqAGjgUzqclZUBmzebI5fdeO01KlqYnw9cd522YxgLr9WDq0lJC3ZWDiHhc2CMtQPQG8A6N+9NADABAJKTk5Genh5M0fyiuLhYymkggcj5ww99AdQHACQnb0d6+gmfx7Rvfxb27qVf9Vdf7UF5ue+C/HXhWvpDYmJPAJQB9vPP25Cff1LTcaF4PTMyGgOgErOxsaeQnr7Fp5xRUd0BkP1p6dKdcDjyzRfUKDjnlj4AJALIADDa175dunThdmDZsmVWi6CJcJfzxAnOGeMc4DwykvPCQm3HvfMOHQNw/s9/mitjsAm2nOPGiWv54YfajwvF6/nVV+J/uf56es2XnPfcI46ZPt18GT0BYCPXOTZbunJgjEUDmAfgK875fF/7SyR6qKgAHnmEymZERrpPfnPHeeeJbbUzW6KfSZOoimnr1r4duKGOnjBWd/vl22jRAFhoVmKMMQAfA8jknL9hlRyS8KVFC2DaNNrWU4b7nHPI2VhRQclzJ07U7Q5xOTmkYFu1ooJ7epK/unUzTaygc9VVpOTy8oCzz9Z2jFo52K2ng5UrhwsAjAOwnTG2pea1JznnCy2USRKm6Ikvj42lZkDrajxg69cDl11mjlx24I8/qIAcQAPk999bK49VdOjg3BtbC+rVkt1WDlZGK63inDPO+Tmc8141D6kYJCFB//5iu66blo6o/PGtZDyhLuy8cpDlMyQSN6j9DutqxdDVLdQhqP4oB4eDBsatW42TyS707k2Ti4MHRcE+uxASoawSidF89hnw3/9Sme4xY+ivHvr3B3r1Iv/DoEHmyGgX1CsHdUkMLXBOxQ7PnKHnJSXUZtRuFBaS3NHR+o6rXx/o188cmcxGKgdJWLJiBbB6NT2aN9evHDp3lglwCoGYlRij8t5K+Y3cXGoAZCfeegt4/HGqj5SaStdg2TL9isJuSLOSJCxZs0ZsDxhgnRzhQKA+h5YtxXZOTuDyBJOqKmDOHOqb7XCQeWjr1tBTDEePAh9+SIrMKOTKQRJ2FBWJapgREUDfvtbKY2c4r9vKISoK+Pe/aQWqoDXHQYFzMkvl5wMpKUBiorEy7twJ9OhBn9OkCXDvvSR3oMiVgyTs2LJF5DWkpdnTxh0qnDxJs2aAkgjr19d/DjsrB4BqST35pHiuJ88DAIYPp0J9Xbo4K5lAqKoCFi6kXJxu3UjpAJSTY9RnSOUgCTsyMsR2IKuGkyeB558nh7aWPhDhSKCRSoD9lUNEBEUdKehVDklJYtuoXIdly4DLL6dVzKRJwKhR4j2j8lCkcpCEHZs2ie1AlENEBPDcc9RJbtEie7Z6DBQjchzsrhwA/0pnuNvfqFwHpSd6YSGt7NSTlx9+0FcRwBNSOUjCDqNWDo0aiQbyFRXA7t2ByWVH6tUDhg0jk0jnzv6dw47KYcsW4J13gB07yBHtT7luBTOypNWFYMeMAYYOFSa//fuFzy0QpHKQhBUlJc7O6J49AzvfOeeI7bqYxDV0KLVM3b0bmDnTv3PYUTnMmUOO6B49KIxVrRwyMoAlS7Sfy+iVQ2kpsHcvbTMGnHsulXwZMkTsY0QYtlQOkrBiyxbRreuss4CEhMDOp1YO27cHdq66ilo55OZaJ4ce1E2iBg50HtQXLiTTjVaMXjlkZop7vFMn0aGvRw+xz44dgX+OVA6SsMIof4NCmqpxbWZm4OerizRqRO1CY2MpIa601GqJfHPZZdRHvHt3mpGrVw4AhY9qxWjloB741dVh1dt65POEzHOoI1RXUz/lefNoBpydDTDWF927U+KMYlu3O/fcQ+1AMzKM+Z/UJafrknJITwcuuMCYZC/GaMXQsKG+6rhW8swzzs/VyuHWW4ERI7Sfy2izknoFq14tqJWDESsHqRzCnOpq4H//o6gbxU4pqI+8PP0OtlAmMpJme927G3O+rl1pQOOcejuUl9MMOJzZvRv4xz+oPPU555CSaNOGImL8Ta5q1MhYGYNNejplIR87Blx8sTDlaCFYK4euXen+r64mp3RJSWBmVakcwpg9e4Dx470nxYwZQ0t+iXvq1aMGNwcOkJ13zx7tjV7syqRJNMDs2UOPefNoBVFWZrVk1tGmDT38oVEjUqpVVcDp04FPMDwph9hY4I476PO6dw98lSaVQ5jy9dfAnXeKapgA3TQ33UT21PbtgfT0DIwcWdswv2QJObr0NjYJV7p1I+UAALt2WasclPh1s8wznAPnn08O2VOnxOstWwIFBXW7I56/RESQryUvj67fqVP+r9YLCkRiYkxM7fDi998PTFY10iEdhpw5Q+n+imKIiqLZ4IEDwLvvUmZlWhqQllaEdu2cj121irp9XXoppeLbiUOHaDmtRHIYRag4pbdtozIMSUnA4sXmfAZj1Hd73z7nkhEnT1LbVX9t5oWFwIYNFOVjVHkHs3juOWDqVEp+LC425pyZmZQrE6gZt6AAGDyY7oNu3YypoeQJuXIIQ+rVoxT688+nEsOzZzun/3uitJTaQZaV0eBw/fU0CIVaBUpPTJ8OvPkmJQO9+SZw++3GnHf0aFpppaU5OwCDTadOpPArK+m72bCB7MxmkJTkrBSLiujv/PnAXXfpP9+CBcDNN9P2DTfQvRmKcE73UWEhPc/JIYVZUUEr76NHgU8/BQ4fpufDh2s7r1E+lw4dgOXLSU7lOzELuXIIU3r1opIP69ZpUwwAOdnefVc8//13mkXZBSVJrajIWPPHwIEUBTV0qLVmlfh4UcKjqAi45hpzw0LVpTMUli/371x2SYQ7elQohgYNaLU0ezYpy7g44IkngMmTyXwzZ451cjJG8pmJVA5hQFUVkJVV+/UhQyh8UA+jRwMvviiev/IKsHFjYPIFA84pAU6hVy/rZDELxkhhKwEEmZnkWzILddG9O++kgXP2bP/OZRflEB9PA/8DD1AwB2MijLWiwrmIXna28SbMQOEcuOUWijBr1UpU1PUHaVayOQ4HmU8WLaIetYGWiwDI1rx0Kc0Sq6vpR7JpEznAQpXsbLKLA6QQwyVvw5Vhw2g198QT9Pyzz4wzn82dS+HOw4ZRAqF65XDRRYHZyl2VA+ehmfPQqBEwcaLza+ochzZt6PfRrBltc67tnygsJH9Yfj6FtmpdzeuFMQq7VTrvZWVRXSx/kCsHG8M5Nfb4/HO66YYNE1E1gRARAXz8seiDsHMn8MEHgZ/XTNR1j3r2NG/gOX1amB3M5uhR6mL32mvOK8Px4ymeHaDExv37jfm8jz6iwIUBA4AvvnBeOejtHe1Kgwbifjpzhq6jXVA74Zs3B6ZMoZXF6NFAZKS28qfffAP06UM+CrXpVi/TppHP4/ffyWLgjvbtxXYg94ZUDjbmhRecQ9euvRa1oo/8pWNHZ3/D8887hzaGGmqTkhGrJ1fefZcGyIYNg6co580jn9Gjj5KpQKF5c+cM3c8/D/yzOHdurTpkiDHluhUYs49pyZVAynW7O87fiK+yMioCeNttFE3oCXUIeiCTRakcbMonnzgP3v/6FykKI2fM//63UDYnTlB4X6hitr+BMTFYBiucVe3wHDvW+T21svjii8Dr91dW0v10661kQmrTxrlInjKwFxT4//+Hg3Lw17RmRJa0YioCSFl7CmNVrxwCUQ7S52BDFi8GJkwQz4cPJ2URYbCqj40FXn6ZQg8Bmj0/8oj/syczUZuVzFAO6hpLu3YZf35XcnLIZATQ9zp6tPP7V15JIbtFRWQ62LePQl39JSYG+M9/xPPSUjIxHTlCn5GVRSU1Dh+marf+KIhQVw6cU0vQZs0oBPzhhymM29Ws5A9GrBzUysGbT80os5JUDjZj0yYyH1VX0/PevYFvvzUvF+H66yli6dQpWtLqjX4KBhs3irpRkZHO8flG4VqAz2yH6o8/itXA0KG1B6W4OGDQICofDVBGcyDKwZX4eOCll8TzkhLhg/jrL1IeeuoLAaGvHAoLyZQH0P/2+OO0olKSQRmjTOft26m3xeHDQMOG7TF0qO9zG7FyUPudvCkHo8xKUjnYiIMHKbu5pISep6YCP//sX9N3rTBGCXUtW5qbjekvpaXAhReK5y1amFMrqkULimQ5dYpm0jk5gdvhvaHu9HXFFe73ue46+l6GDaNZvZkkJFCphoMHqW7PsWP6/VuhrhxcB1/GnAfypk1p8pGfL3x9PXpomy01bEgTuMpKyrr2R7lqVQ7SrFTHOHWKaiIdPUrPGzUi81JKivmfnZpq/mf4S3w80L+/MMF07GjO5zBGqwfFaZuZaZ5y4NxZOQwb5n6/8ePpESyWLiUl6e8qtV07Cqts2TI063a1bUsZ4FlZYoLhzt+gjtzKy9M2E2GMVg+KUszP1/+7UisHb8cqE6SyMvIR+RtIIpWDTUhIoJIDf/5J9uEffnA2ddRlJk4k5TBmDIVjmoWrcjBrtv7XX2IS0Lixczc6Mygvp2zrs88mf82NN7rfz9+qpApjxtAjVGnUiK6DmoIC8vk4HEI5pKZSiY02bYD8/B0AztV0frVyyMvTrxy0+hwYo9WD4hfyd/UglYNNiI6mga9NGxqkBg+2Tpb16+nHMmCAdTKouekmMrGcOWNu34BgOaXVq4bBg40PNHAlM5OSKBctopXXjTdSlvyOHbQ6uv124/pj2I2LLiJT0MmTIts4Lo7yHAAgPV17ZT61U9ofv4NWsxIglUOdgzFrax1t3Uqz9HXryM6vmHKCSUUFOQJdzUcxMeZncAerK5xaOaibxnuDcxpw/Ikk27ZNbCuFBX/7TdRRGjGi7ioHQJTcDhS1U1pvxFJ1Nd33Cr5WHY89Btx9NykJf02tMs8hhFm1KvD4dSNp1ozabwIkmzp8NFh88QXZrW++mTqWecNop2ewSndv2CC2fUXCHDlC5q3GjZ0d83q45BLK4J08maLTlPMqqG3sZ87QqkmJkpJop21bGqz799cfRJKbKzKik5NFtrknhgyhQIbu3f0P0JArhxDlf/8D/vlPYNw4MieFQl2jlBQKo/3f/+j5xx8Db78dvM93OCjvwuEgJdG9O9CvHzlKe/em6qmtWwOzZtE1W7OGluKBln5QSE0Vjr569SjqJDHRmHOr2bWLwnNXrvTtb2jShGb4VVUUinnqlH7TWkoKKQVFMXDuPjuacxqYlD4hJ0+SUtLDkiXkN8vNpZltKAU7XHwxfa+pqRS+bfR3O3Wq/4mkkZGUlHrokDn3nFs455Y9AHwCIA/ADi37d+nShduBZcuWBXT8L79wHhPDOf0cOX/sMWPkcsUfOZcuFXI1bsx5aanxcrmiyJmby3nv3vTZjRpxfvo05//5j5DnmWdo/4suEq9NmWKsLLt2cV5c7FlGK1CuCcD5b79531eLnCdOiPMlJnLucIj30tLEexkZ+mUdNkwc/8svgclpJKWlQq7ISM4rK+n1jRs537SJ8+xs8RrnnH/3HecjRnDeufNp/uqrQRXVLwBs5DrHZ6vNSrMAjPC1U11i5Upq5F5RQc+7daNknFBh2DARR11QAHz3XfA+u0ULSgLcto1WLfXrA5s3i/f79KG/t91Gf6OinEMRjaBbt8CatpvBuapgGSPKq7uuGtTJfl27kg37kkv8O3eo5jqo7fnq0hR33EH3VevWzvdabi6Fku/ZUx9//RVcWYOFpWYlzvkKxlg7K2UIJTZsoCQ3pYFLairdgOoa8lYTEUGD79NP0/OPPiLzVzDp0YMenDv/YJUyyKNHA6+/TlFMgZSZtgvnngt8+CFtG6EcvFVjnTcvsMzwUFUOqan0+8vKEg2VAM91ldTbRk9AQgXpcwgRtm2jqBCl9V+LFhQxEko2WYXx44FnnyXb/7JlNNM0M1vYEwcOiPLZSUkiDj8+HnjooeDLEyicA7/8QjNVPVFHgawcfv+dIlu6dKFKn+PHe6/GGmjJkFBVDrGxdB3V19LhcA45VX8n559PrU+zszNw1VV9NX1GeTmFCx87RlUOQv0eDXnlwBibAGACACQnJyNdHecXohQXF+uSc/fuRDz2WE+cPk2ppw0aVGLq1C3Izi5xmsUZjV451fTq1RObNjUG58CUKXtx/fXmCepJzuXLmwI4GwDQrl0Bli83P3yKc+D48VhkZdVDXl4sRo486lVGPRw5Eod//YuSRzp3LsJ//5uh6bjKSobo6EGorIzAgQPADz+sQsOG7ov9u8r53XetkJHRGRkZwOnTOWjX7i/88UdbAGQ7rK7OQnq6AU1CaigoSAZAcbHbtuUjPX2nJjmtoLAwClVVFAKWkFCFtWtXOb1fvz7Qpk0x/vorXZNpqbQ0EtdcMwgAEB3tQO/eKzQr29deo449zZuX4+qrj6B+fQ/NHIxEr5PC6AeAdqjDDunVqzlv0EA4wxo0ICdYMAjE6ffxx0LmPn2Mk8kdv/yynN9zD+c//cT5mTPi9cmThQyPPmquDArl5eSwBDhnTDinjXCg/u9/4v+59FJ9x/br55+j9777xHHTptFrd9whXpsxQ58cvli5Upz7vPO0y6mF0lLO8/M5LyvzXz41O3cKWTt3dr+PXjkTE8U5Cwq0HeNwcJ6QII47cULXR3LO7emQrtNUV5PDS+mKlZREy/y+2lapljJ6tAiv3bRJVEU1g02bGmHmTIrbHjhQ/brY9tR2MT+f+iyrfROBEBMjkoo4951roYcM1UJB7z2gOOMB594Wvnj2WWDFCnLwjxxJr3nKcVD46y/yhc2c6bkbmSfMMCtxTsmZCQkUapuURCbZQDGij4Mr6vMoJVJ8UVAgim0mJuoPH/YXS5UDY2w2gDUAujLGshljBnXDtQeRkVTxNCWFbuply+yhGACKpb/zTupStnmzeQXvAGDNmiZ/b6vrGblzRqt54w36Md50E7VWNAqzMqXV/oJztZXr+Rt1PoQ649kXTZtS6e/bbhNZ0JMmUbe7p58WGdNqhgyhIpD33uusSLSgLhSZm0t2/UD59Vfgv/8V5zpzBrjrLmfHsi/69qUOgqNGiexl9eBtlHJo0UJsa3VkuxbcC1bvbaujlYIc5xJ6dO5MSVxK1U87EUgvXD0MH34UXbq0woIF9OMF6Ier/Hjr1aPr6IoS0QQ4l6QIlLQ0KnwIGKccHA7nlZDeSYLSGjUmRt+g6I5Bg+jhiTZtxLU/csR3nR818fE08y0ooFXH8eOBNY/iHHjyydqv791LEVz33OP7HA4H9WiorCTFqoQqqzvhuat+fP/9wMKF56K4mCZ56lWtJ/xZOeipqWQk0qwURMrLgdWra7+elmY/xRBM0tKK8OqrZMJRBi31qqFXL1qFuTJwIP3QL7yQQoQVRREoZqwc9u0TkVdNm+qPUuvblwrlFRdTKQwzOe88Kkh3883+9RIx0rQ0b54wx8XFAQ8+KN57/nm6Hr7IzRUKtWlT7cohKwvYvz8ReXnaB3r1yiHUlUPIRyuFCwcOUOXQ7dtppeBtZibxjLKk7tGDTB+bNlEYpjsSE6mchNFNisxQDq4mJb2mg7g4/cXxuJ/d7N55R/8xaq64gpRZy5aBV9FVr17vvx944QVSGIcPk3lowQLPJcgVWrakgTorSyhowH0PbTXqFY9WE5E/+RFSOYQxCxZQQ3il6cZVV5GSsCI3wCyOHqVGKRMmBKdjXOvWzn20PWGGLGedJbb37AncjAME5oz2l7ffphpCHTqQ/+iWW4LzuS+/bNy5/vlPmumnp1MP7Lg4+l++/pr8IloUJmM0aLv6FZo2pWuTm+t+5fD448B5523EyJHnajaNBbpyCGbekzQrmUhxMTntrrpKKIboaFruupuJ2JVbbiFFd++95FQ3CqPMQEaTmCgS7qqqjInUCsTdPx3xAAAgAElEQVQZ7S9799LA98cfItnryy/pu7zggsBXCMFg4kRqlXvypBjAn3ySVnRvvCF8Mf4wfTqZ+0pK3JcL6dwZ6Ny5GK1aae+OZ6eVg1QOJvH772T6mDlTvJaaSrWT7r8/eBEHwaBxYxEp8u23xp136VKaub3+ehcsWWLceY3ASNOSqzPaX+XgcNBgNn8+Vab1hboJjNK2c/9+8gOsXq0/EslKYmPFtjv/UyAwZlzDJX9WDlo7wBmNpkU3Y6wZgAsAtARQCmAHKKnCgEC08OLYMeDVV7vWqnd/1VUUS96kifvj7Mz11wNvvUXb8+cDM2b432dYzdKlNIAdONAS7dtTeQeAihJGR2tXsLt3k1JevZpMAV27Bi5bt274W2FlZtJM21+Ki+n+2LiRVpj+mhsPHgQ6daLtpk1pReftGn3/Pdnm9+8XobAHD4r327Vzf1xFBU1+srPJRv/ww/7JGypUVASvJH7r1tSjo0ULZ/OkJ86cEau6qKjg9IxX8KocGGPDADwBIAnAZlB57TgAVwPoyBibC+B1zvlpswW1A4sX00BZVCS+wcaNaXl+443htVpQM2AA3fTZ2cCJE2RaUgbyQPjjD7Gtzm+YPh2YNo2aptx9N3Dlld7PM2mSqB47cKAxyqFXLxpQu3VzbgLkDw0aUH8KgCLa/L1P2rUjk1dxMYWIHj3qfTCJiqIKu0qVXUCbcqiuJnu+co4HH9Q3W8/NpcKIOTkU7fTBB9qPVaispM/Wcq1KSymE1hN9+pBMbdsCc+eKVZQeqqq0+bdatdJneo2IoP4pWVlk8ps2jfo6BKUysLf0aQCvAkj18F4USEmM0ZuW7e8j1MtnHDninOZ+9dXUgyBUMbJmvrqvwu23G3PO8nLOV6zgfNy4A7ywULw+erT4rA8+8H2el18W+991lzGyuRLItXQ4OP/8c85feIHzn3/mvKjIfznOP1/8r4sX65ezfXtxfGam5/2aNBH75eTokzErSxzbvLn7fXzJ+frrnLdty/ndd7svN+NwcP7aa9Q/Ij6e85IS9+dxOJxLWuTn0+sHD3L+7rucz5vH+bZt7o8tLOS8bdtinpTEecOGXsUNmO3bnccVvcCP8hlBGdSNeoSScigv5/zkydqvT5nCeWpqMV+0KPgy6cVI5bBmjbh5GzfmvKLCsFPXkrNVK/FZW7b4Pv6PPzgfO5aUxNq1xsnlTUY9TJsm/h+A844d/a8PdNdd4jyvvKJPzqoqzqOixPHqOlauTJjA+S23cP7kk5zn5emTsbJS1KcC3DeM8nU9x44Vx8+c6X6fbt3EPitXut/n1CnRWKtePdHY6NtvxbFXXun+2OpqziMiHF7/D6P497+FPFFR1OhKD/4oB01uFsbYF4yxhqrn7RhjBlQvsR95ecCLL9IS1J2t9ZFHgI8/3ogRdayF0XnniTC7ggJjatu448gR4SitV09bqOL559PSnEIPzZErEG66STihhw6lss5qB6se1NE53np8FxeLhlIKOTmiVlKzZt5NMR98QE7vKVOo9IseoqKc/SrqRjtaUffZ7tfP/T79+4vt9evd79OwIZmdcnLIjKmYqXzlOABk8mnUSFzEEyc0CO4n06ZR2ZIpU4C1a4NjVtLqg18FYB1jbCRj7E4ASwBMN0+s0MLhIAfcuHE0AD7zDNlzP/+cokPUxMQAUVEhGoNpIoxRkp/CnDnmfM66dWL73HODk1NhJnv2kA35iivI1j1/vvtSIFpRKwdvNZamTaPBv21batgEaPM3GIU66kYdqqmFEyfIiQ7Q781Tn231RMCTcgBokE9JIT+Sgq/saIW3396MvDzygWgNJFi4kPxgt96qvaxLbCzt++STlAdjVPSUNzR9BOf8AwB3APgBwAsABnPOfzRTsFDg4EHgueeoqNzFF5NDqLxcvN+8OcK2RaA/KA3qAXIAu85MtVJWRteVu9GxauVg9Spgxw5qGD92LPDjj/6FkaxdSxm9zz1Hju1AK26efbbY/vNP5/tVzb59NOk5dEgMNMFUDupkLr3KQZ0T0rOn50gjLSsHT6jDTL0ph1atypCcrG+SsmgRJQLOmuUcwuzK2rXkLE9Lo/Ivf/6p/TOMQKtZaRyATwDcDOr7vJAxFkB6SWjzzTekDNq3p4Q19Y8GoGXsV19RmKUStSGh66LMCE+dolBUf/jjD4ooSkmhqq9qQkk5rFsHTJ5Mq6SMDP9GdaOT3+rXFxVyq6qAXbvc71dcLEwoSnSOXVYOWkxKAOUZKea5Awecu7r5QuvKwR/UKwxvzbzmzKEaYpmZtNoIxmpBjVZ9NwbAhZzzPACzGWPfgZSEhyr69mbuXDIjqWncmGzDt95K5aHDNSw1EBgDbr+dMm/HjnUOP9XDypX0V2mnqFBV5TyY6lEOx44B779PP7a4OPJBBIranLF/f6KuYw8cIFu3umyGWjlUV1OfhaNH9ffo7tlTmDu3bXNfzvyHH2hVkZUlejboUQ7Hj9PqMDubfD+PP65PRrVyUCd5aUGrcoiJof997VpxnNKzQiEnh/wOrjZ8M5WDkl0PePe37Njh/Ny1dEZ1tfEJf2q0mpWurlEMyvP1AELQtaePwkL3Mce313SViIigvs7ffEM30Tvv0DJPKgbPPP008Nln9CP0N7GouloUZFMXKNy8WSiLNm3cN6LxRGUlmW5++IFmYUb0EejeXczmsrPjceaM9mOnTaPj16wRryk1lfbvp4Hgoosoh0BvQx210vLmlI6NpaKF9erR87feogHpp598r4iPHaPaVi+8QMmdevF35cC5s4nIm3IAnE1L6lWnwtVXU25IcrKziUevcigrE/3ffaEe5L0pB/X93aQJTWpWrKBSHlpriwWCV+XAGHuKMZbk7j3OeQVj7CLG2BXmiGYe+flkDkhNpY5mZWXO719yCdVlOXSI7IPXX09fjCQ4vPgiOR23bnUepFasENtDhug7Z6tWIqqmqKi2qdAf1H0kOGfY6b4dci04B36s8dgpSqpVK1FaoW1b4W/Jy9NvnuvXjwbFO++k6CetJCaSwrr8ct8rB/WAqbUMhBp/fQ5HjojPS0z0nWWsVg7qVaeCsmo5fpyyygFaUSkNfyIivDf6Wby4OerXJ+e+qwnUE+qVg7dV0623im3F9FdRQffDkSOeTYZG4custB3Aj4yxMgCbAOSDMqQ7A+gFYCmAqaZKaCDFxcBLL9EKQD3L++4756V7ZCRVeJRYR0RE7SiUEydotlteDgwerO98jNF3Hx9PTls9qw5v9OwpWoVu3ep7JgvQirVjRxrkqqvpNfUgFhlJ9+P8+eT70uukHjmytvnEaBo3JuXTrBkpCodDn03cdfas1USiNin17ev7GHX7VNdWsVVVpGBOniT5lZBV9YDdqpV3Z3NMjOPvnhHHj/uWXzknYzQBOHrUc/kOtdJUlLU6G3/fPv/LrmvBl3K4lnN+AWPsMVDpjBQApwF8CWAC57zUHLGMZ/Fi6tfsWkysa1eyOUqMp7SUrvtllxmz8po6lcKIN270L9zTjGX4OeeIsF2t7TkbNQKWL6eZoVIgz9V/8vLLVGLCDJTBuHVr/0OBGaPWnP6SkEBmvmbNnFdKvlCvzrQoYsVsduYMTS5OnqQe0wD973v30rXIzxfXQk8V1EaNKv8+l6LofREdTQo1J4f+7yNHnEuYKLiTIyWFamJ160arCTNN3L5ujb6MsZYAbgIwzOW9eFARvpCmuhp46qnaNeR79SLT0jXXmOvUqas8/jhVpC0uphnwNdcYc964OOrsFipoTTpzh9rGrV45AMYULvTEyy/TdxMVRX/vvJN8CEVFNEMNVu7Is8/qP+app6jf9YYN2mogRUZS8EGbNjTrdjdDj4x0rpaalERd7rKyfPfW6NGjEKdOUX0sPQN1aqrognf4sHvloC4Fr/yvjFGRxmDgayH4PoDfAJwFYKPqkVHzN6SprAT+9S9nxdCsGSWvbdoEXHut9YqBc7Idvvsu1aV3Ze1aKmx3552UZ6F1hmU10dGiReNXX2k7proaePNNchwa0UAnGKgTpzZt0j57LC4W0SiMBa+HAyASyKqqhB9m1ixajdWrR+a3UKZlSxoge/TQtv+VV9L3pDVAok8fCqpIT/e9eouO5mjYUP8M3pff4fPPnYMVAkmM9BevyoFz/jbnvBuATzjnHVSP9pxzP2oXBg+HgyqhqkMWhw+nH+S4caETcfTmm+QEvP9+4JNPar9/9CgNlh99BPzf/4WO3L646Sax/eOPotmRN3bsAB56iJSh2T21HQ5jFG2rVsI5W1ysvbfDpk3CGZ2W5l8vZm/s20fZtJdf7txXGSB/gTJTVmakSoJVZWXgiXgS3/iKWHrySeHLAkQp9mCiaQHJOb/bbEGMxtXJc8895IgOdiKJL9S5AGvW1HYwqQdVX6WpQ4lu3WhJnpFB38XcueTz8YaS3wA4OxKrqhjuu4+u1cUXBzaQPvoo2ft37SLZAi3fzRiZhH74gZ6vW+ecpezKkiVkwunRg3JpNm4kp6g7qqspMuW330iZLFmi/f7Ny6PJBFC7/tTXX9PfkhLhC1IPRFqvSUYG2b9zc6mfhTq6Ri9aS14byY4dtMJNTfVeR8oM/vEP+i5TU2v3ky8sdPaNRkc7rzQUKivp2pvWOlRvpT4rH3qrsjocnM+YwfnDD4tqi8HAtaJkZSXnH37I+eWXUyVHVxl796Yy1G++SfuqKSjgfMkSzl96ifPNm2t/1vz5nL/3nn//n5FVWd3xxhuikuTgwb73X7qU83HjOG/XjvN33hGvT5++yaliaSBccomQad68wM6lMH8+5xMm7OXp6ZwXF3vf9+qrxed/8on3faurOW/aVF8FWoWiIs4Zo+MiI0XFUHffucNBlXSVzzl8WNtnfPyxOGbcOO2yKezaxflFF1GVXdf7w52cq1Zxvn69fxVrHQ4qw/3dd6Ik+oUXCvl//13/ORU5HQ4q3713r/dKtlo5epTze+7h/NJL6bp8+63z+3l5VHE2KorzZs20nROyZHdooL6xKyo4P+cccRP+/HPt/f1VXNnZ4kd94436SwabrRxycjiPiBD/+5492o+tqhLbN9548O9zTJwYmEzqvhMvvxzYudRovZYtW4rP377d9/7XXSf291Sa2hOdOoljMzI8y5mXJ/ZLSNB+Py5cKI77xz/0ycY55wcOiONdBzl3cvbpQ/vGxHC+fLm+zxoyRHxWejq9lpIiXtu3j16rrKRy4I8+ShNLX9di2bJlfNAgcR5PpcGNpLqaelS49qDwhj/KIcSMLOFHdLRzwpa78D9//QgvvUTlsQHg119F4k6okJLinMTmzqfiCXWgwLp1ordqoLWs7riDzDRHjwKPPRbYufRSVQVMnEh+gPbttflVrruOysDPm0cBFHpQ54l4K/DmalLSej+mpVFo8XvvkY1cL6mpwqyVlyfuZXdUVQFbttB2RYU6+ZDMsb4y3tWRTUrgQPfuFJ0VGytMMzk5VBHh1Vcp+1vLtVCy+QF99Zv8JSJCJP+1aCGingz/HHNOaw0FBVQCwFMlSqt46ilyXL70EkUcGcUbb1DsflwcsGCBibbHALjzTrH96af6S0Hk5AD79pFRPjqaSkoEQloanaN58+A796OiaDB96y2KQtMSKXfddTRQjR6tv2+COjxWiXzZsycRK1eSTVsZUNXVPvX4YNq2pcKUd90FDHMNdNdARIRzFI5aSbmyfr2zAsjLo0G8Tx/q17F4sffPUteX2ryZrv2vv1KNq6Ii/3IcFJo2pSgvPfkagfLTT5S7kZvruWR5oISVcpgxgyIzzjpLlCcIJgsX0qywutp51GnWjG7CyZM9Ox/9IT6eislt3UoRPqHIyJEiMuboUffhut74/nuxPWiQ8VE9RpObC59lNCZNIuXUvbu+fsJ6Of98sa0oh9mzUzF4MCXAff45vaYelLU0vTcStTLyphzUJUQaNqT8kg0bxGri7be9f446wMF1FaXOKfFHOXz4ITn3Dx4kJa6VSZNIqXbsWDuB8s8/vUf4tWxpvhM9bJQD5xSbDNCX5G2JasZn33MPKYaFC4H582t3/TArqYkxygJ1paJCe8y9mURHA+PHi+ezZ7vf74YbaIX1yy/OOQ7qHIkxY0wR0RB2766PXr3oR/vAA57341w0eNm1iwqqmYW6GVJmJmUHZ2eLEUUJj1RXhtXSWc9ItCoHdWnrG2+kv/feS/e/0rTI2/3es6dYKWZmwmORRHXNLa3Kwd9cqcxMuhf273dOoHz9dVr1NW5MqxKzzEa+CBvlsGaNyChs0MC5K5nZMOZcnOv771vpNp8YycmTVDzwiSesk0HNhAm0rJ81SyhwNdnZZCKYMoUyqZWl+f79wOrVtB0V5dxMyAgKCoxr7ZiUVPH3D3zlSudS42oyM4VdukkT72Gv7tBzX8XHOyfprVsHdO1ahH79yE7euTMNqOoqp8HukaFVOahrKimThPbtgW+/pTyBDz7wPkgnJopJlMMBbN/ufj939YzMQm0OUsvz8suiwmthIVkerCBslIN60Ln++uDHLSvt+669Fpg5c5Nl7Svz8sicsGIF8NprVFTQatq3p6X8Lbe47428fLnYHjhQ5Kco8fgAJTAqVTMD5cMPyQeUlET2fyNITi7/e6CvqKjd/vHgQSrJMGmSeG3IEG15CwcPkiO9e3e6DnpwNS09/PBfWL+eJhDNmtHfgQNpQtWqlf6ChF98QSvD4cNr90DRgtqM5Uk5lJY6D57qkhZjxmhffan9Dp99Rj6H/fudVxzqSqdmJ56pM7yV/+/MGecCfh06eM7/qK6miZXeLndaCQvlUFZGM0+FW24x9/MOHao9g4uOph/HnDlAw4bW1X5o2tTZyedphhRKXHYZzQDvu4/MSwDN7tTOe3XGdaDExIilutYy21oYMUJsuzpIMzJoIF2wQLymdaCPiqKeCbt20exfj7lw4ECxrZ59M0aP5GQaJAsK3Pc78MXKlTTQLlmiPTtcjXrlsGeP+/9tyxbxepcuztFBelArhy+/BC69lOz9itycO/9etJbncDhoQM/MrF351Rvq8ys+B4fDOdnVU1mVsjJygrdpQwmIZlgqwkI5rF5Nyy+ANO0FF5j3WT/+SPbLp56q/Z7e4ltmEBFBjsZevWjV8Mwz1srjCXVUR1ISrbjeeUdEN+3fL2ZQ8fFVGDXKuM9WZvixscZGl6gH+19+cX7PVUlHRmp3XrZuLcpJl5aKLm9aGDqUVkcbNjgrJlciIpzbV2pF3ddB3SBHKw0aiICFigoK3HDl3XfFdiA1qNROaXVjHiXMNStLvN6kifYOcAcPkpJNS9NXYLJTJ7GSzsmhVVxionPQhadAk7g4Ueakqqp2tWkjsMj4YSzqJjDDh5s3QC9bhr8HqVdeoSW7kYOWUTRuTDPVUCsVovDrr2RemTXLs829Uyf6sf73v0Bm5iEkJBhXyqtHD4oG6djR2JINF14oykPv2UMzSSWX4brrqFzDvHn0/KKL9JnJZsygAatPn9otLb3RogXw739r318vo0aR8kpJ8d+ZffbZFMkWF0fKQW3O4ZxKryh4KkC3bx9FwlVUUF6IO9Qrh4gIGnhLS0UnPHXEUI8e2scR9feoJ88hKooUirLa2L6dTI1KW1PAuw+ofXtaUbVv79nHFQhhoRzUNXlc65QYyZAhFJq5cCEt57x1iLKaUFUMzz5LyUUANVT65RfPsiYk0D7p6YcAGKccYmICr6nkjrg4mpwofp6ZM2k1BNDAqY460etcv/rqwOWbOxfYsqUZ6tWjwc8Iv1zfvr7LWvvi+eep+1+fPrUrp+bk0ICvoJgd1WzbJkqnN21KxRvd3VNNmpCPIzaWfsvPPit6OwDOqzs9uQP169NY0KABrSAqK7VHJ/bo4awc0tJE1dzYWOeAAleWL/e/Fa8WLB1CGGMjGGO7GWN7GWN+xdZUVDiXtjVTOUREkN34ttvoCw12ZEcgOByUOBOsJB1PjBkjZmRLl4ricOHCPfeI7VmzhLlz40Zxn0ZGGtffQg9TpwJTpqThvPMoQGDZMvLzfPGFMW1T/eX882kW726g45xylwYNIt+UO6V+9tkiQfD4ce9Nl7ZsocdbbzkrBoAUz6xZpFz0ZOIzRiufv/4C/vhDX9i6q1Na7fdxpyzVmKkYAFhXWwlAJIB9oClhDICtANK8HeOuttLq1aLGSPv2vmuM6EFPLSA1Ztcs0su2bZxfcAFdo6+/Fq9bJecjj4jvDOA8NZXze+/lPDeX5FPXVbJSTj0oMjocnKelif/t9dfp/VGjxGtjxgRfPoeD6iYpMhw75izT448HXyZv6P3OJ07k/KqrOH/7bSpcFywCvTd/+UV8B927OxcDDLSOmBrYrLZSfwB7Oef7OecVAP4HQHePo9atKS748suBK64wTrjp02kJqo6CsiuffkozGoBmRcps1iqmTnWOojl0iOrzXHIJJThdeql/Tev1UlFB0UpKLoURMOZs41cihJQZaUJCYKslzmmWr7VfscJzzwm7dIsWZPpQZ6vffrv/MoUC779P2fT33x/a5l5X+vcX5r2EBOeVgxUNftQwbpGdgTF2LYARnPM7ap6PA3Ae5/w+l/0mAJgAAPHxnfs+//y36NfP3PTn779vibfeooyZmJhqvPXWFpx1VpGPowTFxcVINLJORoCUlERi/Pj+KCiIxtixhzFuXBbi4hyWypmfH4uJE/ugoMBN4gOAsWMP4a67yPhqhpwHDiTgjjvOhcPBkJpags8+2+D7IC+oZSwtjcAtt/RH27ZnMGXKDsye3QbbtjVChw7FuPDC4+jZ0z/tPGdOa8yenYpTp2Jw7717cO212kNUVq9ugsmTyYYRFeXABRccx/LllF3Vs+cpTJ++xS+ZAODttzshKysBJ07E4PXXt6JJkwrfB7nh5Mlo7NjREPHx1ejW7TBiYuqjqioC9eoZl+pfXQ28/35HtGtXgo4dS9C1a1FAASxG3Js//piC1NQz6NGjEJddNggVFZTNN2vWOrRt67kTM+fA8eMxyMuLQ0FBNC680HNG57BhwzI45/pivfQuNYx6ALgWwEeq5+MAvOv9mL78qacMW2l5JC+P865daWk3cCAtwfUQimaQ337jfOdO59eslvPECc6vuMLZxBQfz/nUqc59LcyQs7hYfGZkJOfl5YGdz1XGPXs4Lymh7WHDxGd9/73/n/HOO+I8N96o79iqKs4HDHC+1srjyy/9l4lzznv0EOfasMG/c8yfL85x8cV0PWfN4rxBA84feMB/E68rmZnic2Jjhdlvxw7qyaAX5Xs/eZL6bfz6K+d//umfbNXVnH/2GefPPku9X3yVCy8tFf9LRAS1B/AEbGZWOgJA3d+odc1rXgmGSSQ5mZJ67ryTyjtblb5uJBddRJEQoURSEsXeb9xINZQ++YScepMmmd8VLCGByiO0a0dmLKPvq06dKESSc/+jYFzp14/+NmzoPtPcG5GRZNqKiXGehaem6isW5w4lBwPwL9cBcI54+v13ICurHl55BTh9mpzHejPZPZXwVjury8vpnuOckiybNwfGjnXf09kX77xDkUWXXCIKGuolIoKy6J97joJHfK1o4uJEjojDYXwNJitDWTcA6MwYaw9SCjcAuNHbAU2blgctyiM11X3vBYmxMGZMOKQ/7N4dhIgPUNz69u2U56C1mJs7evem/IzOnf0LVe7aFbj33r2YPp1Cfq65Bpg2LfCQ1smTqdhgSor/dvLUVIoQWrSIButHHunp5FfR4oMqLKTBfuFCyl9Ytar2PmefTfK+8golj+3cSRVdlbpYCxZQeRW9qHMd9PqDAqF3b2o7m5pqfKFNy5QD57yKMXYfgF9AkUufcM69FjNISqpwqhs/dy59yV27Aldd5V+xPYcDePxxWiW4q24arpSXAwsXtsCgQf5XlbQ7wVAMjFGyXceOgecqBJKfcd99lAHdpEkFdu6k1Yd6xh8IRoWPP/ccKQcAOH7ceWmkdcX16KNikDx2rLZzOi2N+qqcPi1yUB58ULx/5ZWUr6CX1FRSPK7la8xm4UITT67XDmXlwzWU9emnhc3Nn1A8h4PCxQDOmzfX1rZRC1bb8n3x7bfUpxngfNYsa2RIT9duRw7168l5aMtYUCB+J9HR1bX6lIcSrj4ogF5bulTb8ep2oN56hO/ZI3psqx96fUKh/L2rgc18DgHj2t5QLxs3UkEzgGYZ7spJhyM7doikp6eeoiV4MOGcKnl27kwzaqXUusQc1JVG27Q5Y1nFYC289hrN9iMiOIYOJevAjz8CF1+s7fgHHyRz8OHD3n0pnTrVLn3TooVz8cRgM2MG+Tz+8x/zKq3qIYRvE98Eqhz69SOH85gxdFO88opxsoUyDz9MeQVlZZV45JHooJuV9u0Tyikvz9r2pocP0+D5559kejH6WugppaCV8nKykW/ZQuZQX47L9u3p+965EygtzQMQOmHWrnTtSt/J77+vxPDhg3Uf7810x7nztXr/fTIhFRZSvtSECfod/UayciVVdQbIB6du82oFtlIOeXlxuPtuutEdDipupuCvLXbwYPqR0WzFGDlDnfr1gR9+AI4fX4srrjCx3ogHSkspYXHZMmqTGAzbvyf69aNVI0Aydexo3LmrqqgIYmoqlUn48svAFQXnNJApTs/LLqM6X95ISaE+z4DxdaoAkuW66yhSKSqKVqaBEB0NxMZ6CDcKgBtvpMqnI0fSdosW/kcWmYE60ktrRdjKSiqFkpVFfpQ77jBOHlsph1OnovHxx1TQ7MgR0eqvSZPA2i36U6rY7gwYAKSnW9NHtEcPMhVUVBjXic1f0tKEcsjMNFY57NlDWcmZmVQK2ogVBGN0/ZTe02vW+FYOZpOYKJobRUXRxC3UJlrl5RQeWlxMYeojRoh6TEaxeTMN8Pn5ZI3Qmxv3yitkYs3J0R52XlwsynrXq0eZ7kZVpbaVcgBIU5aVOZuUtON+07cAABhTSURBVEYZVVRQpMLjj+sreywxh5gY7TMksxgwgO6LtDQRM24Ue/fSD5Vz7Y1jtDBwIJW2HjRIf+c2M4iLowY8p07RaunEicAH3sWLm+Pdd8kiMGqUf0UuT54kP+LYsWTGLC6m1zt0MCcycexYYc3o31+Ua9fKgAGe+zd4olEjUkLFxTRZPnHCuI6JtlIOycnlePVVmpX442945hnSzj/9RHVYrLR1hxrHj1NcuBGloe3E1KnmnfvKK2nFsHOnsT1Gnn+e+m2HEj/9ROGxLVoEtopX2Lat0d9hrc2b61cOr70GPP00TSQLCqhM/KFDFCqrdMEzmuRkoRzy8/UrB39gjMrEV1ZSDo2n5D9/sJVyaNy44u8WoErNc0BbXPHmzZTwo2z/8AMV6arrVFZSIbjXXqMf0l9/md9YvS6RkGC8Y1FPtNFPP9Fq+bzzaAZuVvCB0d0Xs7NFZp4//sQ2beh+BqjA4PPP02sTJhgkoBsGDKBZfHKy/61M/UHdDMlIbKUc1KidN1qW1r16USbkI49Qo4977zVPNjsRFUUNd5T2iM89R2UWzIBzqpnfvz8NVFZXnawLrFpFlT7XraNInJEjrZZIG3ffvQ8JCX2xf79o5KOHa66hKK3CQho8g9G+9/XXzf+MYBIWykGL3ZoxClUcNoxq+oSaw8wqGKOVw5AhtAy+SnfRdO1s306henPmUOev/HzjwzwlzqhLQNupOVW3bkUYOtT/42NigMWLaZxo394wsUzjxx+p9UDLltR+YPx4qyWysXKYMYOccrm5+oqZ+dvnNpwZPJhssZdcYm4pDXWD+8suCx3FsGMHJR3t2kUF2NS9hv3l1ClajbVubc6stayM2kSuWEGz43ffdb/ft9/S/7Z2LZl+/vzTeFlcMSO3wx+6dLFPSZzdu0VfkVCJnrSVcjh6NA5jxlDiz4gRoVdl1M4EIzN0wgSaGS1YAFx7rfmfp5Xp00WmfGqqMcph/nwKK2zQgJr/vPhi4OdUU1QkvrPoaHJQN2xYe7+mTcmUpJiTzFIOS5cCd99Nk7URI8yzg4cr/uQ4KJSUkGn4wAGKvJs0yRiZbKUcTp+Oxvz5wIUXahvMFiyg2dVTT1EyksRamjWj/tu33Wa1JM6oo0rUpSYCQSnTffq0OUl+ycmkxDZvppn64sUUSmkVUVGiDIq/ZbvtzsmTtJLLz6dE0xtu0H7sQw9RdFtOjn4fS3Ex5VUANBl54gljVqu2Ug4KWmrvV1ZSmYi9e8nBOn8+2dUl2vjzT0oc8scZaDfOO48G1m7djLtHGBOx/2efbcw5XRk/Hjj3XPITaa09ZBbq2e7Jk4Gd64UXgDlzeqNXL2DiROOqvprNvn34u6VAnz76lEOrVv6bk5o1owS4M2doMlJQQH7VQLGVcmjevAxvv03mJF9ZmHPnipmMw2HeDzTcyM2lldasWRRVtHp1cCI9rOTCC+lhJG+8QdErOTnuzT1GoO5V7UpFBflSevcOzvfXoQOtulJSAv9/N20Cdu5siJ07aTZtF9TJZ/n5wftcxqhJUGQkOd+N8hvaSjk0bFiJ66+naJfevSnh5p//FPkLaq6/nkInp0wB/vUvYxJz6gKVlVQDyOEgJ+aCBYFHMB05Qorc6mzoYMOYdc7FpUsp6qVDBwrbfughcz8vOtq4pC91DlMHY8tAmUpyMtXnSk4Ofub6e+8Zf05bKQeF3FxK08/OJmeMOyIjqbjWDTfQgCfRRmoqORbfeosyL/0taKjm6aepjMHQocCzz1J0lMRclOqe+/fbzwfw00/AvHlb0KBBL0Puv2CRmEghqeGCLaP99Xj2IyKsLcNrR556in6gixYBZ50V2LlOnABmz6aVyO+/h0aIYzjicADz5on+x/Xri45m/nRItBKKGDv1d7RXuLNnD1UlGDiQcrFChbBXDhL9NG1KJgkjbNWnTpGzlDFy0uktLBYs9u2jMhNjx1INrkBYuZJMckrWudmsXEn1/6+9VoTMvvMO9cpYuJDKkgcLzskhevp08D7T7hw5QiW3166l9gGhgq2UQ25uHEaOpJheBVfl8PnndKEloUHHjrQK2buX7KKh6tw+dIjMX3PmiD7G/vLYYzQLbNCAQqnNprJSDCoZGWRyBWjFfNllwbvmL7xAUTNJSdSNTaKNnByx7W9f79OnyRT84IPGlQaylc+hqCgaixY5Fx5TK4dDh4Bbb6XtUaOAb76xtpFMuMA5DZoNG/qfLNehQ2g7F9XRbDt3+t+TgHPnZjfBsJkPG0YhuMuXk0/HqjagsbGi2J2/fg7Xbm12Y+1aUtD5+ZR4qKXo4ujRZFrKzdXfA0KBc1IMABAfTxnzgV5HWykHBWVmBDgrh48/FiVri4qkYjCC3FzKSP/5Z4oOU7rmhRvJyZRZ2rFjYGHPJSUUsbJ9O/lbgnGtGKMIs4ceos+2CuW3mJAAVPvZR2riRFq5JSX1whtvWJ+/oZfZs6nAJ0ArRy3KIS6Oelp36uT/5zZsSIm+BQXUafHYscD7k9hKObRoUYaZM0Xj8IgI56YiffsC//gH9YVW2iJKAiMmBti4kbaPHiWbtqc6PmoyM2nWbKcCh0b0dkhMpAECoAEyWLPg1q1FhJJVXHcdJYHVr+//OfbupSjE7OxGqKgwTrZgoR6PgpnrAFDSb1QU5Tr4uwJRYyvl0KBBJfr2Fc+bN3dO+Bg1ih6HD1PWoCRwmjQhP86IERTi+uqrvo/ZuZMyd7t2pTyTkSPtbSrwFzOLGIYi8fG+9/HFwYNiO5TNkJ7o359WP8nJwa/IMHmyseezlXIAaPaq4GnZZHVP3XDj0kvJTKKloi3nwD33kO1561bgySep2qs08Um0kJlJvsPvv9+K9u3tV7vl0kvpoYdQ9bPYaNFPHD8uto1uEC7xjDvFUFVVO8GQMTLPxMXRY/ZseyoGI9stSrQTG0tNoPr1K7DlfeMPTZrQhLZ/f89JvVZgu5XD8OEUv52fL+zZ/kaWSALjqaeAr76i0tTPPSdev+ACihQ7dsx+ZdUfe4zqSW3dSiYOPWVXzpyhXIk+fcisVhfbrXJOhfeUJjsJCVZLFNoUFZETuaCAxrR69ayWSGCrIbW0NBJXXgmcfz5FBCjZu2PHktL49NPgJR7VdebPB155hZyHigNWzahRFOVkN5YtA/74g8ogb96s79gtW6ir3nXX2acdp9EMHkxJlD166L9+dZFjx8R2y5aBm5deeIECdrp1CzwR0VYrB4eD/Z2gpITNFRZScbiKCmDJEpq1BhItIdFGXh6Z9fLzrYurN4PevUV01q5dFP2mlYwMsa0OnKhLqCuTqpO7tFBYSCsNu99P06cL68Z773n/fzp1Iv9cTo4xWeVz5lBACEBl97WE0nrCViuHmBhhCN6zh/6uXIm/Q95697ZPW0C7c9ddFBU2Z45IvgkHJk6kjO6cHO8lsd1x3nlklrr44rrbOyQlheL7/QljnjiR/FQdOgDr19u3O9eLL9IK8qOPKNfFF7GxZIIzoneKkY2rbKWjo6MdeO89clj16UOvXXEFlcuYMyfwpA+JPmJj7VfUzReBzPj79w9sphYOvPMOMHOmf8fu30+5IQcOAHFx9o0ISEkRDY9yc4ObNDpxIpl009ICL6FuK+UAkH17xQraTk+nGVpqKvDII5aKJZFIEFhuR2mp2E5JKQtcGIu44w4yEbVs6X+tJH/RYwb1he2UgzqU1YhWeBKJJDTYvp0UxMGDQG5uudXi+I0eM+vJk2SGC0U/i618DoCzcpDd3SRmcuKEaDUrCQ7x8WQOqSuh6cOHk3m2VStg2zarpXHGkq+AMXYdY2wnY8zBGDtXz7FqB8+zz1LIoURiJGvWUJh006bA/fdrO2bGjI64+WZyQqonMHUNzik8c8sW4NdfrZYm9MnOpjytnBwqnGckxcWBtS+wSj/vADAawAo9BzkczKna40cfkU3P6oJjkvCieXNg927aXrPGd7Y058CyZc3wxReU23HkiPkyhioOB/0me/emMhLl9rUOmY56LIuIMC6gZvNmiphLSgqst4Mlli7OeSYAMJ0ZH9XVtffv0qXuhg1KzKF9e/qhnjhBJo7jx70Xcty3DzhxgnrRNmwYWMlvuxMZScpV6edw9CjQtq3v4/btozIrrVrZ36R07Bjw5pu0GkhIoFwHd0RG0nUqL6frZFQL3QYNqCUvQEmd/iroEHSDeKaqSiiHtDTKcZBOaYnRMEY/rnbttFUabdcOmDkzAyUlfVFaWveqsbrSrRuZSFJSnHuveOOBB6hnSEwM8P33xlR4tYrycqoeANA18KQcFGJjtSlQrXTsSMl1e/fSdna2f+dhnHPjpFKfmLGlANwtlCZzzn+o2ScdwCOc841ezjMBwAQAaNgwtW9hIRnRzjvvBF5+ebvRYhtCcXExEo0oqG4yUk7jsIOMQOjKOX58P2RlUSGmDz7YiJYtj4aknK64u56VlQyXXkrmjIgIjiVLViAy0pxx1hO7djVAcnIZkpMpQ3jYsGEZnHNd/l1wzi17AEgHcK7W/Zs3787Jwsv5zTfzkGXZsmVWi6AJKadx2EFGzkNXzmHDOE9Opt92QUHoyumKJzlffZXzjz/mfNEizisrgyuTOwBs5DrHZ1uZldQ+B3UNF4lEYm8UG3lxsTFdzKxGS1LuoUPCES3zHGpgjF3DGMsGMBDAz4yxX7Qeq5S0lTkOkmBw6BCVg/jqK/fvHzxIhdMkxhAOikErjzxCfRxiY4G5c62WpjZWRSt9B+A7vcclJVUgP5+yKE1ylUgkf7N4MXDZZbTduzdw00219xk3jmL6+/fvhtRUe7a2NJrKSqoIeuAA9bi44QarJQpNlJBnhyM0LSEhuJjxjZ0jGST24fzzKbywspJixw8fdm5Bm5tLvR84B9LTm9WpWa83cnOBc86h7eRkqRw8kZJCJqWjR4HWra2Wpja2VA4SSTBo0IASuSorgSuvrG3yyMqiUMG9e4FzzjmFZs3sW2baSFq1Eko1P9+3H2HtWpo9d+hAORKh2E9ZL2vXAu++S7kOF1xAZbxdUUxJRUWh1QFOQSoHicQLP/7oebAaMAD46y8qGLdq1X4AdbTDjwuRkcDQobTdoQPF/XtTDpMmUYVlAFi0CBgxwmwJzScvT/ipYmO97xuqzclspRxOn47GqlWU4CF7N0iCga9ZLGNkQjl5UvanVbNkifZ9lVIlAK3EwgG1mSiQ+kZWYivlcPRoHAYNoqxoLR2WJBKj4ZzMJKE627MbVVW0yvjzT4r8at/eaomMoWtXYMYMakxm1+6UtlIOCp07Wy2BpK7y5JPAyy9TJNPw4VZLY3+iooCvv6ZtzsPD3wBQTaV77vH8/po1pBg7diTHdCj+37YscdWpk9USSOoic+eSYgCAqVNlOLXRhOIAaRbPPgsMHkzO+59+sloa98iVg0SikcpKci6Wl1Nuw6FDxhZMCydKS8nJrOQ6yDa+zqibSIXqZNeWyiFUL6YkvPnnP4H+/YH168mkJCsCe6akBBg5krYTEoCHH65bKwNXHA5RipxzCm9t2pRKlYeqn8WWykGuHCRW0bFj+ETUmEmTJuS0LyoiRXH8OCXEuTJ9OkUedutGfTDCqdz5li3Ao48Ce/bQ/6aYjxgDvvjCWtm0YEvlIFcOEklowxhwxRXUn6FnT/eNbMrLydykdEQrKQnNZDB/iYwEli6lbaMa+QQT2ymHpCS5nJdI7IASheSJvXuFYmjbNrwUA+C8wjx0iKKTQrH6qidsJCohVw0SSXiQmAg89hiQmUllM8KNevUo47tdO/Ir2EkxADZUDtLfIJGEB23binaa4Yq7UiCPP05tjgcPJsURqo56WymHmBgH0tKslkIikUj8IzcXmDaNtmNjgcJC37WXrMJWSXDt2pXgySetlkIikWhlwQLgiScorHX/fqulsZ4VK8T2gAGhqxgAm60cJBKJvZg5E/ilps/jpk11txlSaSnw229kSps6lZTEsGFWS+UdW60cJBKJvVCa/gDAtm1ie948YOxY4M03ySEdzjz9NOV9XHklsHMnlShftIic8aGMXDlIJBLTGDmSYvzPOYfMKAqLFwNz5tCjrIyS4MKVRo1o5QCQme32262VRytSOUgkEtMYOlQ0/lGzdq3YHjgwWNJYw6hRlOx31lnUi9wuSOUgkUiCzuefU9nqNWuAc8+1Whpz6dyZaijZzd8ilYNEIgk6vXvTw1vPg3DCbooBkA5piUQSJDgnp3RFhdWSSLQglYNEIjGdd9+lrOCePYHu3amEtSS0kcpBIpGYzqFD1CcaoDpDEXLkCXnkVySRSEzn4YeBCy+k7WuusVYWiTakQ1oikZhO8+bAypVAVhbQsKHV0ki0IJWDRCIJGrLntn2QZiWJRCKR1EIqB4lEIpHUQioHiUQikdRCKgeJRCKR1EIqB4lEIpHUQioHiUQikdRCKgeJRCKR1IJxzq2WQTOMsSIAu62WQwNNARy3WggNSDmNww4yAlJOo7GLnF055/X1HGC3JLjdnPOQr/7OGNso5TQOO8hpBxkBKafR2ElOvcdIs5JEIpFIaiGVg0QikUhqYTfl8F+rBdCIlNNY7CCnHWQEpJxGE7Zy2sohLZFIJJLgYLeVg0QikUiCgC2UA2NsBGNsN2NsL2PsCavl8QRj7CBjbDtjbIs/0QFmwRj7hDGWxxjboXotiTH2K2NsT83fxlbKWCOTOzmfY4wdqbmmWxhjI62UsUamNoyxZYyxXYyxnYyxB2peD6lr6kXOkLqmjLE4xth6xtjWGjmfr3m9PWNsXc3v/hvGWEwIyjiLMXZAdS17WSWjGsZYJGNsM2Psp5rn+q8l5zykHwAiAewD0AFADICtANKslsuDrAcBNLVaDjdyDQbQB8AO1WvTADxRs/0EgFdCVM7nADxitWwucqYA6FOzXR/AXwDSQu2aepEzpK4pAAYgsWY7GsA6AAMAzAFwQ83r7wO4OwRlnAXgWquvoRt5HwLwNYCfap7rvpZ2WDn0B7CXc76fc14B4H8ArrJYJlvBOV8B4KTLy1cB+Kxm+zMAVwdVKDd4kDPk4Jzncs431WwXAcgE0Aohdk29yBlScKK45ml0zYMDuAjA3JrXLb2eXmQMORhjrQFcDuCjmucMflxLOyiHVgAOq55nIwRv8Bo4gCWMsQzG2ASrhfFBc855bs32UQDNrRTGB/cxxrbVmJ0sN3+pYYy1A9AbNJMM2WvqIicQYte0xgyyBUAegF9B1oJTnPOqml0s/927ysg5V67llJpr+SZjLNZCERWmA3gMgKPmeRP4cS3toBzsxIWc8z4ALgNwL2NssNUCaYHTWjMkZ0EA3gPQEUAvALkAXrdWHAFjLBHAPAAPcs5Pq98LpWvqRs6Qu6ac82rOeS8ArUHWgrMsFqkWrjIyxs4GMAkkaz8ASQAet1BEMMauAJDHOc8I9Fx2UA5HALRRPW9d81rIwTk/UvM3D8B3oJs8VDnGGEsBgJq/eRbL4xbO+bGaH6UDwIcIkWvKGIsGDbhfcc7n17wcctfUnZyhek0BgHN+CsAyAAMBNGKMKSV+QuZ3r5JxRI3pjnPOywF8Cuuv5QUARjHGDoJM8BcBeAt+XEs7KIcNADrXeNtjANwAYIHFMtWCMZbAGKuvbAO4FMAO70dZygIAt9Rs3wLgBwtl8Ygy2NZwDULgmtbYcD8GkMk5f0P1VkhdU09yhto1ZYwlM8Ya1WzHA7gE5B9ZBuDamt0svZ4eZPxTNRlgIDu+pdeScz6Jc96ac94ONFb+zjm/Cf5cS6u96ho97yNBkRb7AEy2Wh4PMnYARVJtBbAzlOQEMBtkPqgE2RtvB9khfwOwB8BSAEkhKucXALYD2AYafFNCQM4LQSajbQC21DxGhto19SJnSF1TAOcA2Fwjzw4Az9S83gHAegB7AXwLIDYEZfy95lruAPAlaiKaQuEBYChEtJLuaykzpCUSiURSCzuYlSQSiUQSZKRykEgkEkktpHKQSCQSSS2kcpBIJBJJLaRykEgkEkktpHKQSCQSSS2kcpBIJBJJLaRykEgCgDHWr6boWlxNlvzOmpo7EomtkUlwEkmAMMZeAhAHIB5ANuf8/ywWSSIJGKkcJJIAqan5tQFAGYDzOefVFoskkQSMNCtJJIHTBEAiqNtanMWySCSGIFcOEkmAMMYWgMojtwcVsbvPYpEkkoCJ8r2LRCLxBGPsZgCVnPOvGWORAFYzxi7inP9utWwSSSDIlYNEIpFIaiF9DhKJRCKphVQOEolEIqmFVA4SiUQiqYVUDhKJRCKphVQOEolEIqmFVA4SiUQiqYVUDhKJRCKphVQOEolEIqnF/wNq4C1fq0pLJwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "# We evaluate the Riemann zeta function using Arb and plot it using mpmath\n",
    "from mpmath import plot\n",
    "ctx.dps = 15\n",
    "plot(lambda x: acb(0.5 + 1j*x).zeta(), [0,40])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: Hilbert matrices\n",
    "\n",
    "Textbook example of ill-conditioned matrices:\n",
    "\n",
    "$$A_{i,j} = \\frac{1}{i+j+1}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 1.0000, 0.50000, 0.33333, 0.25000, 0.20000]\n",
      "[0.50000, 0.33333, 0.25000, 0.20000, 0.16667]\n",
      "[0.33333, 0.25000, 0.20000, 0.16667, 0.14286]\n",
      "[0.25000, 0.20000, 0.16667, 0.14286, 0.12500]\n",
      "[0.20000, 0.16667, 0.14286, 0.12500, 0.11111]\n"
     ]
    }
   ],
   "source": [
    "A = arb_mat.hilbert(5,5)\n",
    "\n",
    "print(A.str(5, radius=False))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.74929513e-12 +/- 3.68e-21]\n"
     ]
    }
   ],
   "source": [
    "print(A.det())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why ball arithmetic might be useful"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.1491893717231122e-120\n"
     ]
    }
   ],
   "source": [
    "from scipy.linalg import hilbert, det\n",
    "print(det(hilbert(15)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[+/- 1.06e-89]\n"
     ]
    }
   ],
   "source": [
    "print(arb_mat.hilbert(15,15).det())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Automatic precision"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "15 [+/- 3.20e-417]\n",
      "30 [+/- 1.69e-767]\n",
      "60 [+/- 3.53e-1575]\n",
      "120 [+/- 3.98e-3186]\n",
      "240 [3.3700336774911741861999225672508298305760992725682801800204310940198298318687928463528636e-5942 +/- 1.77e-6031]\n"
     ]
    }
   ],
   "source": [
    "ctx.dps = 15\n",
    "while 1:\n",
    "    H = arb_mat.hilbert(100,100)\n",
    "    d = H.det()\n",
    "    print(ctx.dps, d)\n",
    "    if d.rel_accuracy_bits() > 53:\n",
    "        break\n",
    "    ctx.dps *= 2\n",
    "\n",
    "ctx.dps = 15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### New feature in Arb 2.16: eigenvalues"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[1.56705069109823 +/- 8.92e-15] + [+/- 5.27e-15]j,\n",
       " [0.2085342186110 +/- 2.01e-14] + [+/- 5.27e-15]j,\n",
       " [3.28792877e-6 +/- 7.64e-15] + [+/- 5.27e-15]j,\n",
       " [0.00030589804015 +/- 6.67e-15] + [+/- 5.27e-15]j,\n",
       " [0.01140749162342 +/- 5.68e-15] + [+/- 5.27e-15]j]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arb_mat.hilbert(5,5).eig()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.00000000000000*x^5 + ([-1.78730158730159 +/- 4.07e-15])*x^4 + ([0.34759117535903 +/- 3.44e-15])*x^3 + ([-0.00383508350430 +/- 2.74e-15])*x^2 + ([1.15292700e-6 +/- 2.19e-15])*x + ([-3.75e-12 +/- 2.51e-15])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[[+/- 3.78e-5] + [+/- 3.64e-5]j,\n",
       " [0.00031 +/- 7.48e-6] + [+/- 3.01e-6]j,\n",
       " [0.01140749 +/- 8.57e-9] + [+/- 8.89e-9]j,\n",
       " [0.20853422 +/- 3.10e-9] + [+/- 1.82e-9]j,\n",
       " [1.56705069 +/- 5.66e-9] + [+/- 5.37e-9]j]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pol = arb_mat.hilbert(5,5).charpoly()\n",
    "print(pol)\n",
    "pol.roots()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example: numerical integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function that both mpmath and SciPy have trouble integrating:\n",
    "\n",
    "$$\\int_0^1 \\operatorname{sech}^2(10(x-0.2)) + \\operatorname{sech}^4(100(x-0.4)) + \\operatorname{sech}^6(1000(x-0.6)) \\,\\, dx \\approx 0.210803$$\n",
    "\n",
    "SciPy gives an error estimate of 3e-9 although the actual error is 0.001"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.209819784443225\n",
      "(0.20973606883387982, 3.0281456740012046e-09)\n"
     ]
    }
   ],
   "source": [
    "from mpmath import quad, sech, plot\n",
    "from scipy.integrate import quad as scipy_quad\n",
    "\n",
    "f = lambda x: sech(10*x-2)**2 + sech(100*x-40)**4 + sech(1000*x-600)**6\n",
    "\n",
    "print(quad(f, [0,1]))\n",
    "print(scipy_quad(f, 0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEKCAYAAAAMzhLIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcVNWZ//HP091sArKjKEYIAkLUICEIro3LuIyjRhOjE5ck/nSM0Rjjy4mZmEyGUSf+kpmMZpzfxExMfmoU1whRDGbEwhUCihrBQFBcQOOCLLaAbGf+OFXeW9Xd1UV3nbq3qr7v16tennvrdtXDtbueOs+551xzziEiItKehqQDEBGRdFOiEBGRopQoRESkKCUKEREpSolCRESKUqIQEZGilChERKQoJQoRESlKiUJERIpqSjqAndW/f3+3zz77JB1GKnz44Yf07t076TBSQecionMR0bmIPPPMM+8554Z05merLlHstttuLFq0KOkwUiGTydDc3Jx0GKmgcxHRuYjoXETM7LXO/qxKTyIiUpQShYiIFKVEISIiRSlRiIhIUUoUIiJSlBKFiIgUpUQhVe+dd2Ddum5JhyFSs6puHoVI3LPPwpQpAFP5wx9gwoSkIxKpPepRSFWbORO2boWtWxuYOTPpaERqkxKFVLVNm6L2Rx8lF4dILVOikKq2ZUvbbREpHyUKqWpKFCLhKVFIVVOiEAlPiUKqmhKFSHhKFFLV4gPYShQiYShRSFVTj0IkPCUKqWpKFCLhKVFIVVOiEAkvWKIws5vN7B0ze7Gd583MbjCzFWb2gplNDBWL1C4lCpHwQvYofgUcV+T544HR2ccFwP8LGEvNWb0a3n67B84lHUmyNJgd2bEDHn4YXnihX93/Xkh5BUsUzrnHgPeLHHIycIvz5gP9zWxYqHhqgXNw/fUwZgwMHw5nnDGVUaPg2mth+/ako0uGehSR++6DY4+FSy89kKeeSjoaqSVJrh67J/BGbHtVdt9bhQea2QX4XgdDhgwhk8lUIr5U2bSpgeuu25d584bm7V+5Er77XZg5cw3f+95L9OmzLaEIk7F27WeB3gC8//4HZDLPJBtQgu68cxSwFwC33voyW7e+UfwH6kBLS0tdfl6UW1UsM+6cuwm4CWDs2LGuubk52YAqbMcOOPFEmDcv2terFzi3nc2bGwH4wx8G8S//ciiZDHSro1szxP+t3bv3pd5+N+Luvz9qjxw5iubmUckFkxKZTKaufyfKJcmrnlaT+/rjDc/ukwI//jE89FC0fckl8P77MGvWE/zDP0T7n3oKvv/9yseXJI1RRHbsaLst0lVJJopZwDnZq5+mAOudc63KTvVu0SLyksEVV8ANN0DPntCtm+Oaa/wYRc4Pfwhz51Y+zqRojCISH6eq1zErCSPk5bF3AE8DY81slZmdZ2YXmtmF2UNmA68AK4CfAxeFiqVaOQeXXRb90U+dCtdc0/q4b3/bD2LmxH+m1ilRRNSjkFCCjVE4587s4HkHfD3U+9eCBx6AJ57w7W7d4JZb2h5/aGiAX/4SRo+GDz+EF16A22+Hs8+ubLxJUKKIKFFIKJqZnVLbt8N3vhNtX3gh7LNP+8cPGwaXXx5tX3VVfdzxTYkiotKThKJEkVIPPghLlvh2nz7+g78jl18OQ4b49uuvw4wZ4eJLA+eUKOLUo5BQlChS6vrro/ZFF8HQoe0fm7PrrvCtb+W/Ri3P0N26NX97y5ba/vd2RIlCQlGiSKEXX4yuXGpogK/vxEjO+ef7K6IAFi+GJ58sf3xpUdiDcK6+Sy4qPUkoShQp9NOfRu3PfQ4+8YnSf3bQIDjrrGg73jOpNW2Vmuq5/KQehYSiRJEyGzf6K5ZyvvGNnX+N+M/MnAnvvdf1uNJIiSKfEoWEokSRMvffDy0tvj1mDBx22M6/xv77w0EH+fbWrXDnneWLL03auqqrnhOFSk8SihJFytx6a9Q+5xww69zrnHNO269ZS9SjyKcehYSiRJEif/mLv59Azpe+1PnX+uIXo8l5CxbA8uVdiy2NlCjyKVFIKEoUKTJjRvQHfvjhMGJE519r0CA44YRo+7bbuhRaKilR5FPpSUJRokiRe++N2l3pTeTEr366776uv17aKFHkU49CQlGiSIm//CWa89DQAKec0vXXPP74aE7FkiWwbFnXXzNN2hrMLpyEV0+UKCQUJYqUmDkzmlV82GGlzcTuSO/ecFzsruW/+U3XXzNN1KPIF08OKj1JOSlRpES8NHTaaeV73VNPbfs9aoESRb54clCPQspJiSIF1q7Nv9lQOcpOOSeeCE3ZxeQXLvSLBdYKJYp8Kj1JKEoUKfDww7Btm29/9rOw117Fj98ZAwbAtGnR9uzZ5XvtpClR5FPpSUJRokiBBx+M2ieeWP7Xj79mLSUKzczOp9KThKJEkbAdO+Chh6Lt+NyHcom/5iOPwObN5X+PJKhHkU+lJwlFiSJhCxdGi/btthtMnFj+99hnH79uFPhFB+fNK/97JEGJIp9KTxKKEkXC4qWg44/3cyhCiPcqaqX8pESRT6UnCUWJImHxD+0QZaecv/7rtt+zmilR5FPpSUJRokjQ++/DM8/4dkMDHHNMuPc67DDYZRffXrECXnst3HtVigaz86n0JKEoUSQok4lmY0+aBP37h3uvHj3y723xyCPh3qtS1KPIp9KThKJEkaD4h/VRR4V/v/h7KFHUHpWeJBQligQlnShyvZlqpUSRT6UnCUWJIiGrVkWrufboAQcfHP49J0yAgQN9++23/Yqy1UxjFPlUepJQlCgSEu9NHHII9OoV/j0bGvKX86j28pN6FPlUepJQlCgSUumyU1vvpURRW1R6klCCJgozO87MlpnZCjO7so3nP2Fmj5rZYjN7wcwCziRID+fSkSjmzYsWI6xGShT5VHqSUIIlCjNrBG4EjgfGA2ea2fiCw64C7nLOHQicAfxnqHjSZNkyePNN3+7XDz7zmcq99+jR0eq0GzbAokWVe+9yU6LIp9KThBKyRzEZWOGce8U5twWYAZxccIwDds22+wFvBownNeK9iebm6H4RlWBWO+UnDWbnU6KQUEJ+RO0JvBHbXgUcVHDMD4CHzewSoDdwdFsvZGYXABcADBkyhEwmU+5YK+rOOz8FDAHgE5/4M5nM6k69TktLS6fOxR577AaMA+Cee9ZyyCHPd+r9k/b22wcAA/P2rV79DpnM0mQCStjGjVOBHgCsW7eBTObZZANKgc7+jUi+Cn6XbdOZwK+cc/9qZlOBW81sP+dc3vch59xNwE0AY8eOdc3NzZWPtEy2b4c//jHavvDC0YwfP7pTr5XJZOjMuRgzBq691reXLh3A5MnNHy/vUU169269r1+/oTQ3l+GG41WoW7eo3bv3rp363ag1nf0bkXwhS0+rgfi92oZn98WdB9wF4Jx7GugJDA4YU+Keew7WrfPtYcNg3LjKx7DHHtH7btkCTz1V+RjKQWMU+VR6klBCJoqFwGgzG2lm3fGD1bMKjnkdOArAzMbhE8W7AWNK3GOPRe0jjvBjBkmIz6d4/PFkYugqJYp8ujxWQgmWKJxz24CLgTnAS/irm5aY2XQzOyl72OXA+Wb2PHAH8GXnqn1hieIKE0VS4gsExmOqJhrMzqfLYyWUoGMUzrnZwOyCfd+PtZcCh4SMIU127Mj/9n744cnFEk8U8+f7D9ju3ZOLpzPUo8in0pOEopnZFfTSS7BmjW8PHpzM+ETOnnvCJz/p25s3V+d8CiWKfCo9SShKFBUUL/Ecdlhy4xM58R5NNY5TKFHkU+lJQlGiqKB4okiy7JRT7eMUShT5VHqSUJQoKsS59CWKeAxPPll95QoNZudT6UlCUaKokFdeidZ36tsXPv3pZOMBGDUKdt/dt9evz58IWA3Uo8in0pOEokRRIfHexKGHQmNjcrHkmFXvOIVzShSFVHqSUJQoKiRtZaec+DhFNSWK7dvbvpWrEoWn0pOUkxJFhaQ1UcRjeeyx6rmPdjwhNDS0vb+eFPYg1KOQclKiqIBVq/wYBUDPnjBpUrLxxO23H/Tv79tvvw0rViQbT6niA9nxxQG3bKmeZFdOShQSkhJFBcRLOlOnpmsGdEODHzPJqZbLZOM9h549oaEhyg7VfNe+zipMDCo9STkpUVRA4US7tKnGcYp4oujeHbp129Hmc/WiMDGoRyHlpERRAWlZ36k9heMU1SCeDHr0gKYm1+Zz9UKlJwlJiSKwNWtgyRLfbmqCKVOSjactEydCr16+vXKlH1NJO/Uo8qn0JCEpUQT25JNRe+LEtu/KlrTu3f3YSU41lJ/ig9nduuX3KLZuTSCghKn0JCEpUQT2xBNRO43jEznVNk4RTwbqUaj0JGEpUQQW/9BNc6KotnGKeKIo7FG0tQZUrVPpSUJSogho48b8+zzEL0NNmylT/BgK+DGV3H0z0qpwjKKpKfqkVOlJPQopLyWKgBYsiK7pHz8eBg1KNp5idtklfyJgvGSWRsV6FPWYKAoTg3P1OfFQwlCiCChedkpzbyKnmsYpChNFY6Mujy1ln0hnKFEEVC3jEznVlCiKXR5bjz2KtsYklCikXJQoAtm2DZ5+OtquhkRxyCFR+9ln4cMPk4ulIyo95VOPQkJSoghk8eLog3avvWDvvZONpxQDB/pFAsEnuvnzk42nmMLLY+OlJyUKT1c+SbkoUQRSLfMnClVL+SleelKPQj0KCUuJIpBqG5/IqZZE0br0VN8T7jRGISEpUQTgXG30KObPT++388LSk3oUrfep9CTlokQRwLJl8O67vj1wIIwbl2w8O2P4cBgxwrc3bvSD2mmk0lM+lZ4kJCWKAArnTzRU2VmuhvJTsdJTPSYKlZ4kpKAfYWZ2nJktM7MVZnZlO8ecbmZLzWyJmd0eMp5KqbaJdoWqIVEUzqPQVU+t96n0JOXSFOqFzawRuBE4BlgFLDSzWc65pbFjRgPfAQ5xzq01s6Gh4qmkah3IzonH/MQT/kMobb2iYvMo6nEwW6UnCSnkn/9kYIVz7hXn3BZgBnBywTHnAzc659YCOOfeCRhPRaxaBa++6tu9evl7UFSbsWNhyBDffv99eOmlZONpS2Gi0Mzs1vuUKKRcQiaKPYE3YtursvvixgBjzOxJM5tvZscFjKci4lc7TZniyyLVxiy/ZJbG8pNKT/lUepKQgpWeduL9RwPNwHDgMTPb3zm3Ln6QmV0AXAAwZMgQMplMhcMs3YwZo8nlw732epVM5tVg79XS0hLsXAwbNhzYB4B77nmbffdNV7di5cp98L8y8NprK9ixI0oUy5evJJN5LaHIkrFkya5Afvf1qafms3Ll5mQCSomQfyP1JGSiWA3sFdsent0XtwpY4JzbCqw0s+X4xLEwfpBz7ibgJoCxY8e65ubmUDF32Te+EbXPPnsEzc0jgr1XJpMh1Lno0wf+8z99e/ny3Whu3i3I+3TW3XdH7XHj9uHDD1/9eHv48JE0N4+sfFAJamrjL3ny5CmMGlX5WNIk5N9IPQlZeloIjDazkWbWHTgDmFVwzP343gRmNhhfinolYExBrV0LL77o242NvvRUrSZM8MkC4I034LWUfUHXjYvyqfQkIQVLFM65bcDFwBzgJeAu59wSM5tuZidlD5sDrDGzpcCjwBXOuZTfW619Tz4Z3Sxm4sTog7YaNTXB1KnRdtrGKXQ/iny66klCCnrRo3NutnNujHNulHPumuy+7zvnZmXbzjn3LefceOfc/s65GSHjCa3aL4stlOb5FK2veqrvwWxd9SQhpezq+OpW7RPtCqU5UbS+6kmlp1L2iXRGSYPZ2YlwhwB7AJuAF4FFzjn9KmZt2gSLFkXbtZAoDjrIf1vfutXPpXj33Wh+RdJ046J8GqOQkIr2KMxsmpnNAR4EjgeGAeOBq4A/mtk/mdmu4cNMvwULog+offdNzwdqV/TqBZMmRdtp6lUoUeRT6UlC6qhHcQJwvnPu9cInzKwJOBG/RMe9AWKrKvFLtQ8/PLEwyu6II6Jbus6bB6eemmw8Oa2vetJgdin7RDqjaI/COXdFW0ki+9w259z9zrm6TxKQnyimTUssjLKLX4KepnlLWj02n0pPElJJg9lmdquZ9YttjzCzR8KFVV02bYq+dUP+h2u1O+SQaDLXCy/Ae+8lG0+OSk/5VHqSkEq96ukJYIGZnWBm5wMPA/8eLqzqMn9+VO7Yd1/Yffdk4ymnPn3gs5+Nth97LLlY4oqVnuoxUaj0JCGVlCiccz8D/g8wE5gOHO6c+23IwKpJvCRTS72JnDSWn1pPuFPpqZBKT1IupZaezgZuBs4BfgXMNrNPB4yrqjz6aNSu9UQR/7cmqdg9szWY3f4+kc4odVHA04BDs/eLuMPMfoNPGAeGCqxabNzoL43NqcVEkRun2LbNr2WVhvkUumd2Po1RSEillp5Oid9UyDn3B+CgYFFVkfj4xLhxsFu6Flkti969YfLkaDsN4xS66imfSk8SUkcT7q4ys4FtPeec22JmR5rZiWFCqw61XnbKSVv5qVjpSYmi/X0indFR6emPwG/NbDPwLPAu0BN/z4gJwP8A1waNMOVqdf5EoWnT4Nrs/+k0DGir9JRPpScJqaPS0+edc4fglwNfAjQCG4DbgMnOucucc+8GjjG1CscnjjgiuVhCmzrVfyADLFkC7yR8d3PNo8in0pOE1FGi+IyZ7QF8CX/ToZ8Bt+BvStQrcGyp99RT0YfS+PEwdGiy8YRUOE4xb15ysUDr0lP88lhd9dT+PpHO6ChR/BfwCLAvsCj2eCb737pWL2WnnPi/MelxisLSk+5H0XqfEoWUS0drPd3gnBsH3Oyc+2TsMdI598kKxZhav/991K7lgeyceKL4n/9JLg7n/KW6OSo9qfQkYZV6eezXQgdSbdauje4/0dAARx6ZbDyVcPDBfulxgD//Obn7aMeTRFMTmGlmtkpPEpLucNdJc+dGf4iTJsHANi8iri09e+YvoR7vUVVSYdkJ1KNQ6UlCUqLopIcfjtrHHJNcHJUW/7fGz0ElFV7xBFrCQ6UnCUmJopPi36b/6q+Si6PS4v/WRx5J5sOo8Ion0MxslZ4kJCWKTnj5ZVi50rd794YpU5KNp5L22y9aRv3992Hx4srH0FbpqbHRj1WAH+yut2/TKj1JSEoUnRAvuUybFn2rrQdmcPTR0XYS5ae2Sk+F7XrrVaj0JCEpUXRCvOxUT+MTOfF/cxID2m2VnkCJopR9Ip2hRLGTtm3ztfmcehqfyIkniiefhJaWyr5/W6Wnwna9DWir9CQhKVHspIULYcMG3x4+HMaOTTaeJAwb5scqwH9zr/Sy4yo9tabSk4SkRLGTCstOuQHUepNk+am90lO8rUShHoWUjxLFTooP3tZj2Skn/m+v9IB2KaWneksUKj1JSEEThZkdZ2bLzGyFmV1Z5LjTzMyZ2aSQ8XTV2rX+jnbgexJHHZVsPEk6/PDoG/zSpZVdzkOlp9ZUepKQgiUKM2sEbgSOB8YDZ5rZ+DaO6wtcCiwofC5tfve76I9v8uTk7xudpF12yV8k8MEHK/feuuqpNZWeJKSQPYrJwArn3CvOuS3ADODkNo77Z+A6YHPAWMrigQei9ol1fQNYL34O4ucmNF311JoShYQUMlHsCbwR216V3fcxM5sI7OWcq+D30c7Ztg0eeijaVqLIPwdz58KHH1bmfTWY3VpbZSaVnqRcOrpndjBm1gD8G/DlEo69ALgAYMiQIWQSuGnz88/3Y+3aA/ExbGbt2vmJ3zu6paUlkXMRN3LkJFau7MNHH8FPfvJHDj10TfD3fO65IcCnAFi37h0ymaW0tLSwceN6oB8Af/jDYjZtWh88lrR4/fVRwF55+1aseJlM5o22f6BOpOFvpBaETBSryf/NHZ7dl9MX2A/ImL/GdHdglpmd5JzLu3uec+4m4CaAsWPHuuYE7hIUr8GfdlpPpk2rfAyFMpkMSZyLuC9+EX74Q99+9dX9ueqq8O+5OvZbtMceQ2luHkomk2Hw4H4f799vvwPr4mZSOb/5Tet9I0aMorl5VOWDSZE0/I3UgpClp4XAaDMbaWbdgTPw990GwDm33jk32Dk3wjk3ApgPtEoSaaHxibb9zd9E7QcfrExdXIPZran0JCEFSxTOuW3AxcAc4CXgLufcEjObbmYnhXrfEFasgD/9ybd79aqPu9mV6qCDYPBg3/7LX+DZZ8O/ZymXx2owW4PZUj5B51E452Y758Y450Y5567J7vu+c25WG8c2V0Nv4uijo9uBil/e+4QTou3f/jb8e2rCXWtKFBKSZmaXIJ4o4qUW8Sp9mayuempNpScJSYmiAxs2wLx50Xb827N4xx4LTdnLIp59Nn+wOQTNzG5NPQoJSYmiA7Nn+zkUABMnwp57Fj++Hu26KxxxRLQ9c2bY91PpqTUlCglJiaIDd90VtU85Jbk40u5zn4vad98d9r1Kueqp3gazVXqSkJQoivjgA9+jyDn99ORiSbvTTouWXJ83z18BFYpKT62pRyEhKVEU8cAD8NFHvn3AAfV5k6JS7b67X1EWwDm4775w79Ve6ameB7OVKCQkJYoi4mUn9SY6Fj9H8XNXbppw15ruRyEhKVG0Y8OG/EUAv/CF5GKpFqeeGpWfHnsM3norzPtoMLs13Y9CQlKiaEe87PTpT8OYMcnGUw123z26+ilk+UljFK2p9CQhKVG0Q2WnzqlE+UlXPbWm0pOEpETRhg0b/N3sclR2Kt2pp0JD9rfq8cfDlJ9UempNpScJSYmiDb/9bVR2mjABRo9ONp5qsttu+eWne+8t/3u0V3rSVU8d7xPpDCWKNvz611FbvYmdFz9nt91W/tfXVU+tqfQkISlRFFi9GubMibbPPDO5WKrVF74QfWgvWAAvvVTe11fpqTWVniQkJYoCt9wS/dFNmwYjRyYbTzUaPDh/ld1f/rK8r6/7UbSm0pOEpEQR4xzcfHO0/dWvJhdLtYufu1tuKe83fJWeWlOikJCUKGKefNLfzQ78iqinnppsPNXs2GNh2DDffvvt/KvIukpLeLSmRQElJCWKmHhv4owzYJddkoul2jU1wTnnRNvxc9tV6lG0ph6FhKREkdXSkj9BTGWnrvvKV6L2Aw/AO++U53U1M7s1JQoJSYki6+674cMPfXv8eJg8Odl4asHYsXDwwb69bVv5LpUt5aqnehvMVulJQlKiyIqXRr7ylWhxO+maeM/s5pv9BQNdFU8CKj156lFISEoUwIsvwhNP+HZjI5x1VrLx1JLTT4/GepYs8avKdtW6dVG7X7+oXc+D2UoUEpISBfAf/xG1P/c5vwqqlEffvnD22dH2T37StdfbsQPWro22BwyI2vXco1DpSUKq+0Sxbh3cemu0fcklycVSqy69NGrPmgUvv9z51/rgg+ibcp8+GszOUY9CQqr7RPHLX8LGjb69//5w2GHJxlOLxo2D447zbefgpz/t/Gu115sAJYpS9ol0Rl0nim3b4Prro+2LL9YgdiiXXRa1f/ELWL++c69TLFHExyhyq//WC5WeJKS6ThT33AOvvebbgwdrEDukY47xlx2Dn7PS2Ql4xRJF375Re8OGzr1+tVKPQkKq20ThHPzoR9H217+umdghmcE3vxlt33BD577xvv9+1C5MFP37R+3O9liqlRKFhBQ0UZjZcWa2zMxWmNmVbTz/LTNbamYvmNkjZrZ3yHji5s6FZ5/17Z49faKQsM46CwYN8u1XX+3cPbXjPYqBA/Of693bX94MsGlTfU26U+lJQgqWKMysEbgROB4YD5xpZuMLDlsMTHLOHQDcA/zfUPEUmj49an/5yzBkSKXeuX716gUXXhhtX331zn/rLVZ6MsufV1FPvQr1KCSkkD2KycAK59wrzrktwAzg5PgBzrlHnXPZa46YDwwPGM/H5s2LJn41NcG3v12JdxXw5afevX37hRf85bI7o1iigPzyU3xiXq1TopCQQiaKPYE3Ytursvvacx7wUMB4PvZP/xS1zz0XRoyoxLsK+IsGLroo2p4+feeW9egoUdRrj0KlJwmpKekAAMzsLGAScEQ7z18AXAAwZMgQMplMp9/rmWcG8OijnwagocFx5JELyGQ2d/r1ktTS0tKlc5GUKVO60aPHFD76qJHFi2H69CUcccS7Jf3ssmXjgaEAvPXWUjIZvyRtdC4+DfgMMm/e87S0rG3zdWrN5s1TgR55+9av/4BM5plkAkqJav0bSR3nXJAHMBWYE9v+DvCdNo47GngJGFrK644ZM8Z11o4dzk2a5Jz/Duvceed1+qVS4dFHH006hE674oro/8OYMc5t3Vrazx19dPRzDz0U7c+di1NOiZ6/557yx51Wu+8e/btzjwkTko4qedX8N1JuwCLXyc/zkKWnhcBoMxtpZt2BM4C8irSZHQj8DDjJOVemuxW07957YdEi3+7RA37wg9DvKO258sqoTLR8een31S521RPU7xiFSk8SUrBE4ZzbBlwMzMH3GO5yzi0xs+lmdlL2sB8BfYC7zew5M9vJoc3SffRR/qD1JZfA8IoMnUtbBg6Ev//7aPt73yttktzOjFHUU6LQYLaEFHQehXNutnNujHNulHPumuy+7zvnZmXbRzvndnPOTcg+Tir+ip337/8Or7zi2wMG+G+0kqxLL4U9s5c3vP02XHttxz+zM1c91dNgthKFhFQXM7Pfestfs58zfXo08UuS07s3XHddtP2Tn8CKFe0fv2NHfi8hnhTa2lfvPQqVnqRc6iJRfPObfn0hgE99Kn/SlyTrb/8Wpk717S1b/KWz7V0uu2FD9Fzfvn4OTCFdHhtRj0LKpeYTxezZcNdd0fYNN7T9ASPJMPPLjjdkfxN//3u4/fa2jy22zlOOehTF94l0Rk0nig0b4Gtfi7bPPReOPDK5eKRtn/kMfOMb0fZll8G7bUyr6OiKJ6jfHoVKTxJSTSeKyy+H11/37YED4cc/TjYead/06dFVaO++68uDhSWojgayoX57FCo9SUg1mygefBD++7+j7Rtv9MtHSDr17Qs//3m0fd998Otf5x9TSqJQj6L4PpHOqMlE8eabfkXYnNNPhzPOSCwcKdFxx8Hf/V20/bWvwZ//HG2rR9E+lZ7rfYicAAALb0lEQVQkpJpLFNu3w5e+BO+957eHDfO9CakOP/4x7LOPb7e0+CS/ObsUVymD2bvuGrU3bKiPb9W5RTsK1cO/XSqj5hLFlVdCbg2whga44w6VnKpJnz7+KrXc/a+few6++lX/Qfj449Fx7d0/pFu3aBnzHTuiy6JrWXsJQYlCyqWmEsWvf50/YP2P/whHtLkeraTZgQf6mfQ5d9wBRx3lL3UGf0ntSUXm8Ndb+am9hKDSk5RLzSSKxx+H886Ltk86Ca66Krl4pGsuvDD/0uZHH43ap58O++7b/s/W24B2ewlBPQopl5pIFEuXwskn+4X/AMaNg1tvjSZxSfUx85MjTzih9f7vfa/4z6pHUXy/yM6q+o/SpUv9JLrcFTG77eYvjY0Pakp1amqC++/3cyy6dfP7zj3XL8NSTL31KJQoJLSqXsxiyRKfJN7J3smid2944AEYOTLZuKR8unXzPYizz/b/v48+uuOfqbceRbz01NQE27a13i/SFVWbKF580SeJ3FIPffrAQw/BpEnJxiVhjBhR+r3N67lHEU8U6lFIuVRl6en3v4fDDouSRN++MGcOHHposnFJOsR7FGvr4JbZ8YSQK9EV7hfpiqrrUaxb153jj4+61bvu6pPElCnJxiXpkbsZEsDixcnFUSnxElM8Uaj0JOVSdT2Kd97p8fEfwB57wNy5ShKSL75C8COPRKWYWlVYesppb8a2yM6qukSRM3kyLFrkl6gWiRs3LlqJdv16WLAg2XhCiyeKxsb855QopByqMlGcdRbMm+fXcRIpZOYXGMz53e+Si6US4iWmhgZoaHBtPifSWVWXKAYP/ohbboGePZOORNLs2GOj9pw5ycVRCYU9inii0IC2lEPVJYqBA7dglnQUknZHHRXNzF+0KH9BwVoTTwYNDeT9fShRSDlUXaIQKcWAAdDc7NvO+VLUbbfV5gdnYaJobFTpScpLiUJq1o03wu67+/bGjX529wEH+P21NGM7ngwaG8FMpScpLyUKqVn77usvethrr2jfkiVw8cX+QohzzvGz+XOLSVarwh5FfDFMJQopByUKqWljxsDzz8N3vxvd0Aj8XfNuvdWvTjt0qL8r4owZ0bph1aR1olDpScpLiUJq3oABcPXVsGqVLztNmJD//IYNcPvtcOaZfvXhAw6ASy/1N8Javjz938pbl56i7bTHLtWh6pbwEOms/v3hoov8DZGeecb3IO67D1auzD/uj3/0j5x+/fzEzokT/WS+sWN9WWvQoMrG355iPQolCimHoInCzI4Drgcagf92zv2w4PkewC3AZ4A1wBedc6+GjEnEzK8yPGkS/OhHvjQ1c6Zf7uPpp1sv+bF+vV8qZu7c/P2DBsHo0X4WeFuPoUOhV6/w/x6VniS0YInCzBqBG4FjgFXAQjOb5ZxbGjvsPGCtc24fMzsDuA74YqiYRAqZ+VLUhAn+HustLX7Oxfz5sHChf7z3Xts/u2aNfxTTq5dPKIMGwcCBUbtvX780fnuPXr2gRw/o3t3/N97u1i2/vNR6Zna0rR6FlEPIHsVkYIVz7hUAM5sBnAzEE8XJwA+y7XuA/zAzc04r1Egy+vSB44/3D/BzMF5/3SeMF1+EZcvgT3/yYxcbN3b8eps2+bGRVavKG2c8gcQTWeHM7BNO8KsY5BJIVx7x5LQz7Ur/XLz91lujufNOEpf0JOH99+/az4dMFHsCb8S2VwEHtXeMc26bma0HBgHtfIcTqSwz2Htv//j856P9O3b4D/9XX40SQfyxerW/X8rWrWHi2rLFPz74IH9/z57Q1BQlivhYS33as+ND6sCJJ3bt56tiMNvMLgAuABgyZAiZTCbZgFKipaVF5yIryXOxxx7+MXly/n7nYNOmRjZsaGLDhm7ZRxMffNCNjRsb2bQp/7F5s//vxo2NbNnSwNatuYfltbdvb/tixW7ddnDIIcvYe2/jttv2rcC/XKrFmjVd++4dMlGsBmJTnRie3dfWMavMrAnohx/UzuOcuwm4CWDs2LGuObc2Q53LZDLoXHj1dC527PC9iY8+ih7btsGAAQ307z+OTCbD1Vfvy5o1/thyPOLjIPHCcEftnTm2HD9X+BrLly9nzJgxJCkNhfS99x7cpV5FyESxEBhtZiPxCeEM4G8LjpkFnAs8DXwemKvxCZHiGhp8ianYCsq5clm9y2TepLk52URRC4IliuyYw8XAHPzlsTc755aY2XRgkXNuFvAL4FYzWwG8j08mIiKSIkHHKJxzs4HZBfu+H2tvBr4QMgYREekaLeEhIiJFKVGIiEhRShQiIlKUEoWIiBSlRCEiIkVZtU1bMLMPgGVJx5ESg9FyJzk6FxGdi4jORWSsc65vZ36wKpbwKLDMOTcp6SDSwMwW6Vx4OhcRnYuIzkXEzBZ19mdVehIRkaKUKEREpKhqTBQ3JR1AiuhcRHQuIjoXEZ2LSKfPRdUNZouISGVVY49CREQqKLWJwsyOM7NlZrbCzK5s4/keZnZn9vkFZjai8lFWRgnn4ltmttTMXjCzR8ysZheY7uhcxI47zcycmdXsFS+lnAszOz37u7HEzG6vdIyVUsLfyCfM7FEzW5z9OzkhiThDM7ObzewdM3uxnefNzG7InqcXzGxiSS/snEvdA78s+cvAJ4HuwPPA+IJjLgL+K9s+A7gz6bgTPBfTgF2y7a/V87nIHtcXeAyYD0xKOu4Efy9GA4uBAdntoUnHneC5uAn4WrY9Hng16bgDnYvDgYnAi+08fwLwEGDAFGBBKa+b1h7FZGCFc+4V59wWYAZwcsExJwP/P9u+BzjKLOlbmAfR4blwzj3qnNuY3ZyPv5tgLSrl9wLgn4HrgM2VDK7CSjkX5wM3OufWAjjn3qlwjJVSyrlwwK7Zdj/gzQrGVzHOucfw9/Zpz8nALc6bD/Q3s2EdvW5aE8WewBux7VW0vkv6x8c457YB64FBFYmusko5F3Hn4b8x1KIOz0W2K72Xc+7BSgaWgFJ+L8YAY8zsSTObb2bHVSy6yirlXPwAOMvMVuHvkXNJZUJLnZ39PAGqc2a2tMPMzgImAUckHUsSzKwB+DfgywmHkhZN+PJTM76X+ZiZ7e+cW5doVMk4E/iVc+5fzWwq/s6a+znndiQdWDVIa49iNbBXbHt4dl+bx5hZE747uaYi0VVWKecCMzsa+C5wknPuowrFVmkdnYu+wH5AxsxexddgZ9XogHYpvxergFnOua3OuZXAcnziqDWlnIvzgLsAnHNPAz3x60DVm5I+TwqlNVEsBEab2Ugz644frJ5VcMws4Nxs+/PAXJcdrakxHZ4LMzsQ+Bk+SdRqHRo6OBfOufXOucHOuRHOuRH48ZqTnHOdXuMmxUr5G7kf35vAzAbjS1GvVDLICinlXLwOHAVgZuPwieLdikaZDrOAc7JXP00B1jvn3uroh1JZenLObTOzi4E5+CsabnbOLTGz6cAi59ws4Bf47uMK/ODNGclFHE6J5+JHQB/g7ux4/uvOuZMSCzqQEs9FXSjxXMwB/srMlgLbgSucczXX6y7xXFwO/NzMLsMPbH+5Fr9Ymtkd+C8Hg7PjMf8IdANwzv0XfnzmBGAFsBH4SkmvW4PnSkREyiitpScREUkJJQoRESlKiUJERIpSohARkaKUKEREpCglChERKUqJQkREilKiEOkiM/tsdm3/nmbWO3vvh/2SjkukXDThTqQMzOxq/LIQvYBVzrl/STgkkbJRohApg+waQwvx98A42Dm3PeGQRMpGpSeR8hiEX2+rL75nIVIz1KMQKQMzm4W/s9pIYJhz7uKEQxIpm1SuHitSTczsHGCrc+52M2sEnjKzI51zc5OOTaQc1KMQEZGiNEYhIiJFKVGIiEhRShQiIlKUEoWIiBSlRCEiIkUpUYiISFFKFCIiUpQShYiIFPW/sXahB6muu2oAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(f, [0,1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculating the integral with Arb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda x, _: (10*x-2).sech()**2 + (100*x-40).sech()**4 + (1000*x-600).sech()**6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.21080273550055 +/- 4.44e-15]\n"
     ]
    }
   ],
   "source": [
    "ctx.dps = 15\n",
    "print(acb.integral(f, 0, 1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.210802735500549277375643255705729154360909186436781190347850505878720613128145500205058689261557642 +/- 5.64e-100]\n"
     ]
    }
   ],
   "source": [
    "ctx.dps = 100\n",
    "print(acb.integral(f, 0, 1))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}