{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wc7L7dzG7D19"
   },
   "source": [
    "# Problem 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "L6RGyGJH7D1-"
   },
   "source": [
    "What is $ \\lim_{\\epsilon\\to0} \\int_\\epsilon^1 x^{-1} \\cos(x^{-1} \\log x) dx $ ? "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "Y0rADprz7D1_"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import scipy\n",
    "from math import pi\n",
    "from numpy import cos, log\n",
    "from tqdm import tqdm_notebook as tqdm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "q0YFG9H37D2D"
   },
   "source": [
    "Let $f(x) =  x^{-1}  \\cos(x^{-1} \\log x) $ be our integrand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "id": "lExnZkKb7D2E",
    "outputId": "f0f07ee6-cb96-4151-c7f0-fde3413568e3"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAGFCAYAAADTktgVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6D0lEQVR4nO3deXhU9d3+8Xv27AkhkLCEXRZBoIIiWNxKwdrSWtHHn/pYsbgWfCxYF9xQa9Vaxa0qbmhroVrtqihKcW3FDUWWSigisiZs2beZzJzfH5OZzCSTZZLMzEnyfl3XXGROzpx85hDIne9qMQzDEAAAgAlZE10AAABAcwgqAADAtAgqAADAtAgqAADAtAgqAADAtAgqAADAtAgqAADAtLp8UDEMQ2VlZWI5GAAAup8uH1TKy8uVmZmp8vLyRJcCAAA6WZcPKgAAoPsiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqAAAANMiqDSjqKxGFz7zkd7cUpjoUgAA6LEIKs1Y8vctev+/h3TZ8+sTXQoAAD0WQaUZhytrE10CAAA9HkEFAACYFkEFAACYFkEFAACYFkEFAACYFkGlGZ/sLE50CQAA9HgEFQAAYFoElTb4bBetKwAAJAJBpQ3+9d9DiS4BAIAeiaACAABMi6ACAABMi6ACAABMi6ACAABMi6ACAABMi6DSBoaR6AoAAOiZCCoAAMC0CCoAAMC0CCoAAMC0CCptYLEkugIAAHomggoAADAtggoAADAtgkobMD0ZAIDEIKi0wXMffJ3oEgAA6JEIKm1QXOWRz0ezCgAA8Ra3oHLPPffIYrHo5z//efBYTU2N5s+fr969eystLU1z5sxRUVFRvEqKytqtBxJdAgAAPU5cgsonn3yiJ554QuPHjw87vnDhQr3yyit66aWX9O6772rfvn0666yz4lFS1Iqr3IkuAQCAHifmQaWiokIXXHCBnnrqKfXq1St4vLS0VM8884yWLl2q0047TZMmTdKzzz6rDz74QB9++GGsywIAAF1AzIPK/Pnz9f3vf18zZswIO75+/Xp5PJ6w46NHj9agQYO0bt26Zq9XW1ursrKysAcAAOie7LG8+AsvvKDPPvtMn3zySZPPFRYWyul0KisrK+x4bm6uCgsLm73m3Xffrdtvv72zSwUAACYUsxaV3bt36+qrr9aKFSuUlJTUadddvHixSktLg4/du3d32rVb8sz7TFEGACDeYhZU1q9frwMHDujYY4+V3W6X3W7Xu+++q4cfflh2u125ublyu90qKSkJe11RUZHy8vKava7L5VJGRkbYIx4Kisrj8nUAAECDmHX9fOc739GmTZvCjl188cUaPXq0rr/+euXn58vhcGjt2rWaM2eOJKmgoEC7du3S1KlTY1UWAADoQmIWVNLT0zVu3LiwY6mpqerdu3fw+Lx587Ro0SJlZ2crIyNDV111laZOnaoTTjghVmUBAIAuJKaDaVvzwAMPyGq1as6cOaqtrdWsWbP02GOPJbIkAABgIhbD6Npb7pWVlSkzM1OlpaWdOl5lyA2rmhzbec/3O+36AACgdez1AwAATIugAgAATIugAgAATIugAgAATIugEoXaOm+iSwAAoEchqEThaZbRBwAgrggqUSgoZBl9AADiiaAShW+OVCW6BAAAehSCShS+2F2S6BIAAOhRCCoAAMC0CCpReqfgQKJLAACgxyCoROnDHUcSXQIAAD0GQSVKhrr0Ho4AAHQpBJUo7ThYmegSAADoMQgqUVrzn6JElwAAQI9BUAEAAKZFUAEAAKZFUAEAAKZFUGmHKnddoksAAKBHIKi0w72rCxJdAgAAPQJBpR0+31Wc6BIAAOgRCCrtwJJvAADEB0EFAACYFkGlHfYUVye6BAAAegSCSjscqXTL4/UlugwAALo9gko71Xi8iS4BAIBuj6ACAABMi6ACAABMi6DSTp/uZC0VAABijaDSThc/90miSwAAoNsjqAAAANMiqAAAANMiqAAAANMiqAAAANMiqHSAYbA9IQAAsURQ6YCla7YlugQAALo1gkoHPPLW9kSXAABAt0ZQAQAApkVQ6SB3HbsoAwAQKwSVDvrtW/9NdAkAAHRbBJUOenXT/kSXAABAt0VQ6SCvjynKAADECkGlgwgqAADEDkEFAACYFkGlg/YUV2vX4apElwEAQLdEUOkE1778RaJLAACgWyKodILiKneiSwAAoFsiqHQCBtQCABAbBJVOQFABACA2CCqdYOfhKr26cV+iywAAoNshqHSSBSs/T3QJAAB0OwQVAABgWgQVAABgWgQVAABgWgSVTvT0+zsSXQIAAN0KQaUT3bnqy0SXAABAt0JQAQAApkVQAQAApkVQAQAAphXToHL33XfruOOOU3p6uvr27aszzzxTBQUFYefU1NRo/vz56t27t9LS0jRnzhwVFRXFsqyYemNLoTxeX6LLAACgW4hpUHn33Xc1f/58ffjhh1qzZo08Ho9mzpypysrK4DkLFy7UK6+8opdeeknvvvuu9u3bp7POOiuWZcXU5c+v12/f2p7oMgAA6BYshmHEbUe9gwcPqm/fvnr33Xd10kknqbS0VH369NHKlSt19tlnS5K2bt2qMWPGaN26dTrhhBNavWZZWZkyMzNVWlqqjIyMTqt1yA2r2v3awb1T9O61p3ZaLQAA9FRxHaNSWloqScrOzpYkrV+/Xh6PRzNmzAieM3r0aA0aNEjr1q2LeI3a2lqVlZWFPczGU0fXDwAAnSFuQcXn8+nnP/+5TjzxRI0bN06SVFhYKKfTqaysrLBzc3NzVVhYGPE6d999tzIzM4OP/Pz8WJcetTpf3BqpAADo1uIWVObPn6/NmzfrhRde6NB1Fi9erNLS0uBj9+7dnVRh5/ESVAAA6BRxCSoLFizQq6++qrffflsDBw4MHs/Ly5Pb7VZJSUnY+UVFRcrLy4t4LZfLpYyMjLCH2RyudGv+is9U7fYmuhQAALq0mAYVwzC0YMEC/fWvf9Vbb72loUOHhn1+0qRJcjgcWrt2bfBYQUGBdu3apalTp8aytJhbtWm/Vn68K9FlAADQpdljefH58+dr5cqV+vvf/6709PTguJPMzEwlJycrMzNT8+bN06JFi5Sdna2MjAxdddVVmjp1aptm/JhdeY0n0SUAANClxTSoPP7445KkU045Jez4s88+q7lz50qSHnjgAVmtVs2ZM0e1tbWaNWuWHnvssViWBQAAuoiYBpW2LNGSlJSkRx99VI8++mgsS0kIBtUCANAx7PUTQ+46nwoKy9sU2AAAQFMElRh64r0dmvXge/rDh98kuhQAALokgkocLF2zLdElAADQJRFU4oCVagEAaB+CShzUeQkqAAC0B0ElDup8bFIIAEB7EFTioM5naNOeUr2+aX+iSwEAoEshqMSBYUjnP/2hrlzxmbbsK010OQAAdBkElTgpr6mTJG0rKk9wJQAAdB0ElThjuAoAAG1HUIkzm9WS6BIAAOgyCCpxZiGnAADQZgSVBCiv8WjpmwXaW1Kd6FIAADA1gkqcuet8uvPVL/XwW9s1d/nHiS4HAABTI6jEmcdr6NWN+yRJ/z1QkeBqAAAwN4JKnNX5fHLaue0AALQFPzHjzF3H/GQAANqKoBJndT5DbFEIAEDbEFTizOsz5PU1jSr/2VemdwoOJKAiAADMy57oAnoan8+Q0SinGIahMx5+X5L0jwUnavzArPgXBgCACdGiEmdew5DRKKkcrnQHP/76UGW8SwIAwLQIKnHmM9RkjMqRkKBSUVsX34IAADAxgkqc+XyGbCHr6BuGoSq3N/g8sMsyAAAgqMSd1zBkDdmYsM5nqDokqISGFgAAejqCSpz5fIZCN1D2+gxVexpaUard4S0qhhEeZAAA6EkIKnHWeGqyv0WlYRG4ak94KHnmX1/r6CWr9Y8v9sWlPgAAzISgEmdew5AlZIyK12uoKqQVpXHXz52rvpRhSPe9URC3GgEAMAuCSpz5GrWoeHw+1YS0ojTXzbO3pDqmdQEAYEYElTjzGuEr03p94bN+Qrt+QltaQmcKAQDQUxBU4szrC29VqfMZYeEkdNPCwxUN66u4vT55vGxoCADoWQgqcebzGaoLbVHxGmEBJDSoNF78jTVWAAA9DUElzryGIW/IEvoen08eb8hzb/NBpaza0+x1v9hdoq2FZZ1YKQAAiUdQiTNfhDEqoeGkNrRFpaZtLSob95TozMf+rdmP/Iu9ggAA3QpBJc58vvCgUuc1VBfSouIOCS3ljVtUaiK3qKzauF+GIXm8hl7fvL+TKwYAIHEIKnHmaTQ92eszVOdrZoxKkxaVyEFl097S4McbdpV0QpUAAJgDQSXOQoOIJNX5fHLXGRE/X9mkRSVy1883h6uCH++g6wcA0I0QVOKs8RTjukYtKp4Wun4ijVGprfNqX2nDYnC7Dlc1WaYfAICuiqASZ02CSuMxKi20qNR4mq5au7e4WoYhuexWOWwWub0+7S9lFVsAQPdAUImzxl0/jWf9uMNmAIUHkyp30xaVg+W1kqT+Wcnql5ksSSosrYm6Lo/Xp0MVtVG/DgCAWCKoxJnb23j35PAVZz1eI7hyba3Hfzywen7jDQsl6Uilf/Xa7FSncjNckqSisugCx/7Sas164D1NvvOfuv7ljTIMuo4AAOZAUImzJoNpveEr1UoNrSqBNVV6pTglRd6w8FBIUOmbkSRJKiqLrkXlV6u+DA7CffHT3Xp1I1OcAQDmQFCJs0iDaRsfawgq/mCSleKQ1EyLSv1+QL1TncpNrw8q5W0PKgfLa7Vqkz+YzDw6V5L0xHtftfn1AADEEkElziKNUalr1B3kqYvcohK568ffzdM7zam+9V0/B6Po+nljS6EMQ5owMFP3zBkvh82izXvL9OV+luMHACQeQSXOmraoNN0VOdiiUj9GJSvZ36JS7Wk6mPZwsOvH1TBGJYoWlfe2HZQkzRybp+xUp04d1VeS9OaWojZfAwCAWCGoxFmk6cmeRi0q7rrGXT8ttag0dP30SfN3/RxoY4uKYRj6bFexJGnK0GxJ0nfG+IPK2wUH2nQNAABiiaASZxG7fny+iOc0dP3Ut6i0MusnO9UfaIqrmt9lOdTuI9U6VOGWw2bRuAGZkqRT6ltUvthTErw2AACJQlCJM3eklWkbt6g0nvVTH0CqIyz4VlrtDyVZKQ71TgsEFXdwinNLPt/tb00Z2z9TSQ6bJCk3I0kjc9NkGNInO4+0+X3FgmEYTJUGgB7OnugCeprG3Txen0+eRi0qgXPaMusnsHFhepIjeJ7XZ6i8pk6Z9c+bU1BYLkka2z8j7PjkIdnaVlShT3ce0ayxeW16X52luNKt5z/8Rqs3F+qrgxWq8xnKTXfpuKHZOvNbA3TSUX1ks1riWhMAIHEIKnHWeB8ej9eQp66ZMSrBwbSR11Hx+QxV1K9Wm+ayy2W3Kc1lV0VtnQ5X1rYaVLYVVUiSRuamhx2fPLiXVn60S59+UxzNW+uwv3y2R7f9Y0uTzRf3ldbo7xv26e8b9mlUbroWnzE62EUFAOjeCCoJ5jOajlHxNFnwLdCiUifDMGSpX6q20l2nQM9IepL/rzI71amK2joVV7U+vmT7AX+LylF908KOHzfEP7B2895S1Xi8wW6hWDEMQ79eXaBl7/rXbxmdl65Lpg/TcUN6yWW3acehCq35T5H+vH6PCorKNffZT3TWtwboth+NVUZSy2EMANC1MUYlwTwhs36S6wNB0wXf/C0qPqMhvEgNuyk7bBa57P6/ysB4lsMVLQeVGo9Xu45USZKOatSiMrBXsnIzXPJ4DX2xu6Td762tHvznf4Mh5f++c5RW/d90nT1poAb3TlVeZpKmDc/Rktlj9f51p+mSbw+V1SL95fO9OuOh97VlX2nM6wMAJA5BJcG8Pp/q6oNJqqs+qNT5ZBhGyGDahlaD0O6f8pDxKYFWluz61pfWWlS+Olghn+Ef/5JTPwg3wGKxaMLALEnSpr2xDQKvbdqvh9b+V5J0x4/GatF3RzY7BiUzxaGbf3C0XrpiqvKzk7WnuFpnP75OqzcXxrRGAEDiEFQSzL+Evr9FJcXp777xeH3yeI1gt06K0y6nzf9XFTrzp6LWP+Mn0O0j+Rd+k6QjlS1PUd5x0L+3z1F904IhJ9SE/CxJ0hd7YhdUdh+p0vUvb5QkXX7yMP1k6pA2vW7S4Gy9umC6ph+Vo2qPV1f8Yb2eem9HzOoEACQOQSXBvD4jOOsnxdnQohLo9pEkl92qJIf/ryp05k9g0GmaKzSo+FtUAkvrN2dPcbUkKT87JeLnj6lfV2XTnpI2v5doGIahm/62WeW1dZo0uJd+MXNUVK/PTHHo2bnH6aKpgyVJv3rtS927eivTmQGgmyGoJJi/m8f/cSCoeLy+sLEoLrs12NpSE9qiEuz6aQgqgTEqrbWo7Cn2j08Z2CtyUBk/0B9Udh6uUmkbF5CLxurNhXpv20E5bVbdd84EOWzRfyvabVbd9sOxuu50f8h57J2vdONfNzeZWQUA6LoIKgkWGjwCYcTtbRif4rRbZbFYlFwfYqoijFFJczWMYekdDCpta1EZmJUc8fNZKU4Nqm9t6exxKjUer3756n8kSVecMlxDc1LbfS2LxaKfnTJCd/34GFks0h8/3qX/++PnTVYABgB0TQSVBKsOCyohXT/1xwOzeQIzgiKNUckIbVGpnyF0pJVWkIYWlchBRWpoVfmik7t//vjxLu0rrVG/zCT97JThnXLN86cM0qPnHyuHzaJVm/brkt9/qip3000cAQBdC+uoJFiNp+E3/1RXw2DaQIuKy+4PKIEWleqQH77lEbp+Asvot9SiYhhGQ4tKM10/kj+ovLpxvzZ14oDaardXj77tn4q84LQRnbpGyxnH9FOay67Ln1+v97Yd1IXPfKzlFx3X6sJ3QKL4fIa8hiGvz5Av8KdPTY6Ffd4w5PUp+LEkGYZkqOFjSTKk4JitQGeo/2mk8wKfN5ocC7BYJIv8rZgWi+SfnGcJPx48z3888DqrxRJ2PHBepNcHzg1+zQjXDVzT//CfYw051vD5kK8dYdIAugaCSoKFdv0EBsx66kKDiv9YoLUltEUl2PUToUWluIUxKocq3Kqt88lqkfIyk5o9b3z9FOWNndii8ocPv9GhiloN7JWscybld9p1A04a2Ud/uGSKLn72Y63/pljnPrlOv593vPqmN/8+0fN4fYYq3XWqdnv9D4//UVP/qHb7Gj0PPcdXPzPP/3DXGeHPvYY8daHnGHIHz/WpztsQThA/oeElEGwsCg8zVmt04afJ56yBz4W+ttHXa+Wc5r9eyOetbTs/ECib/3qR3o//PjQ+/4RhvTW4d/u76TvCFEHl0Ucf1W9+8xsVFhZqwoQJeuSRR3T88ccnuqy4qKkPJA6bJTgF2e0N6fpxhHf9RBqjkp4UOkbFPz25orZOtXXeYItMqEC3T15Gkpz25nv/xvbPkMXiX8L+YHmt+qS72vcm63l9hp77YKck6arTRrT4tTti0uBeevHyqbrwmY+1tbBc5yxbpz/Mm9LsDCd0PR6vT8VVbpVUeVRc6VZxlUclVW6V19SpvMaj8to6VdTUqaLW/yivCfzpUUVNnSoj7JtlRjarRbb6H0z+Py0hx/x/Blsf1LTVIND64P+8peFY8POW4GsVcp3A5wIfB1pnAg0yvvqPAy05hhHeKmOo/hwjtHUm9LkhnxF+zdCvEbiuzwh8rulro2UY9S1VwXeEaDx83rd6blB58cUXtWjRIi1btkxTpkzRgw8+qFmzZqmgoEB9+3b//Vxq6v/DdNiswZkv7ha7fkKDir/VJHR6ckayXTarRV6foSOVbvXLbDoGpS3dPpI/AA3LSdVXByu1aW+JThud2673GLD2yyLtLalWrxSHfjRxQIeu1Zox/TL08hVT9b/PfKRvDlfp7GUf6Pl5U5rsawRz8PkMFVe5dbCiVgfKanWwvDb48aGK2oZQUv9nRW3njD+yWvyD2JMcNiU7rUqy25TstPmf1z+SHNawYy67TU671f/Lhd3/79Zps8pht8ppswT/LTtsVjntjZ7brLLb/GHDarGEB5Emx+iqaIlRH4R89cHFF/a8IQj5Qo4FAlXo8/DXh5zva/s1I/0Z8ZrB10e6ftPzo75m6Pm+Zq5Zf+8CXz/S6xt/TcMwlJeRuFbphAeVpUuX6tJLL9XFF18sSVq2bJlWrVql5cuX64YbbkhwdbEX6MqxWy3BFgZ3S10/7tDBtE3HqFgsFvVKcepQRa0OV7QWVJofSBswIT9LXx2s1IZdHQ8qv1/3jSTp3OMGxXz/IEkakpOql6+Ypguf+Uj/PVCh/3linZ656DhNGtwr5l8bDbw+QwfKa7SvpFp7iqu1r8T/8b6S6rAwUhflr8kWi5SZ7FCvFKeyUvx/ZiY7lJ5kV5rLrrQku9Lr/0xzOZTmsod9zr+Rp5WxC11UsGtD/P11dwkNKm63W+vXr9fixYuDx6xWq2bMmKF169ZFfE1tba1qaxsGipaW+gd6lpWVdWptvtqqTr1ecyrKLfLVVslqc8hbUyVfbZUqystUXGL3H/e4/O/NUy1fbZWOlJQG3+uR4lL/OXXVYe8/0+bRgdoq7TlwWIPSm/4j3rHvgHy1Vert9LZ630Zn2+SrrdJH2/aq7IR+7X6fOw5V6L0tu2SxSD8ck9Xpf1/NSbFIT583Vleu+Eyb95bqf377lm6ffbR+GOMWnZ7EMAyVVHn0zZFK7TpcpV1HqrWnuEr7S2u0v7RaB8raHkKyku3qk56knDSXctKcykl3qXeqU1n1YSQz2aGsFIeykp3KSHY0u91Cy7ySzyt3da1a37oTQKylp6e3/AuDkUB79+41JBkffPBB2PFrr73WOP744yO+ZsmSJfW9ljx48ODBgwePrv4oLS1tMSskvOsnWosXL9aiRYuCz0tLSzVo0CDt3r1bGRkZbbpGWVmZ8vPzW3zNuCVvdEq9remT5tTBCrf6ZSbpgimDdN+b2/SD8f107OBeuuOV/+jU0X30yHnH6sn3vtLDa7frrG8N0B1njpMknXzvWzpc6dFffjZVI3Mb3sc1f9qgN7YU6frTR+nCCPvn/PC3/9KOg5V6+ieTdcLw3i3W5/H6NOWutXLX+fSPq07UsJy0qN9j0eFinfqbt2V1peqJCyfpxBE5UV+jM/h8hn771nY9+b5/X6Cpw3rrrh+PU58E9r3GQlu+v1tS5a7TjoMV2n6gUl8drND2A+XafqBS+0trWnxdboZLg7JTNLh3igb0Slb/zGT1z0pSv8xk9UlPamfrh/l19H4jOtzv+IrH/U5Pb3nsYEKDSk5Ojmw2m4qKisKOFxUVKS8vL+JrXC6XXK6ms08yMjKivoktvcbqis8MEY/NLqvLrqTUFGVkZMjqSpHFlSKbK0VWV4rS0/019srKktWVIp8jOVhzpVyyuhzKy8lWRkZDvXk52bK6ylVjcTV5f4ZhqLDaIqsrRaMG9VVGRuujuCcO66dPvynWV8U+TRwW/Tfqnz7bL3tGHw3OTtasiUMTOkjw5rMmacLwfrr25S/00d5qnfPsF7rvnPEdHn9jRq39mzDq19PZsq9UW/aV6T/7ylRQVB4cwxTOKqsrRTlpLg3vk6qhOakakpOqIb1TNSQnRYOzU4MDvnuq9vwfhPbjfsdXIu93QoOK0+nUpEmTtHbtWp155pmSJJ/Pp7Vr12rBggWJLC1uAtOT/YNp6/f6qWs66yel0RL67jpfcJn4dFf4gmbZ9cvoH65s2gN/uNKtGo9PFosiDrSNZGJ+lj79plif7y7WnEkDo3p/hmHohfX7JEn/c2w/U8xkmD2hv8b0S9dVf9ygL/eX6afPfapzJ+frhu+NDu6V1N3UeX3acajSH0r2lmnLvjJt2Vca3NiysZw0p47qm66jctN0VG66RvZN08jc9G57fwCYV8K7fhYtWqSLLrpIkydP1vHHH68HH3xQlZWVwVlA3ZXDZvEvBBVcR8U/3VEKrKNSH1QaraMSWCAudHpm6IJvUsjqtBVNg0rgt+XW1lAJNXFQliRpw+6SNp0f6uOvj2j7wSr53DX60XjztFqM6Juuv/5smn69eque/fdOvfjpbq35skiLvzdac44daIpA1V4er0/O3OF6+fP9+urIN9qyr0xbC8vCVkEOcNgsGpmbrrH9MzS2f6ZG5aVrZG56MOwCQKIlPKice+65OnjwoG699VYVFhZq4sSJWr16tXJzY/dDzeVyacmSJRG7kOLFabPK422YamwLmZ7sX0K/0V4/wRYVf0AJ7Jyc4rQ16fvPDm5MGCmotL7HT2MT87MkSVv3l6vG441qanFgSvLo5DL1yUzMYkHNSXLYtGT2WJ1xTD/d9NdN2lZUoWtf3qhn/vW1rpk5SjPG9DX91FV3nU/bisq1aW+pNu4p1ea9pdpaWKZ+cx/SHa9vDzs3xWnT0f0ygqHk6P4ZGpmbHrOF93oKM/x/0pNwv+PLDPc74UFFkhYsWBDXrh6Xy6Xbbrstbl8vYg0OW9jqmPb6xaCkxuuo1C/41mhl2vLapou9BWQHNyZsvkVlQDO7JkcyICtZfdJdOlheqy92l2jKsJYH4AYUldXojS2FkqSHFswx7X8sxw3J1qtXTdfyf3+tR9/erq2F5br0959qTL8MzZ02WD+aOCAu6760JhBKNu8t1ca99aFkf7nc3qYtJRlJdh0zMFPjBmRqbP9Mje2foSG9U7vtgNZEMsP/Jz0J9zu+zHC/TRFUeqJAKAlwWC0hK9MaTVpUAmNUgl0/Efb5CchOa0uLStsHC1ssFk0Zmq1XN+7XB18dbnNQWfnRLtX5DB0/JFtj+pl70JvTbtUVJw/X/zsuX0++t0PP/nunvtxfpuv/vEl3v75V3xuXp+8f018nDMuW3Rb7FohDFbXaur9cWwvLtLXQ/+e2wopWQ8n4AVk6ZkCm8rOTTd8aBABtQVBJkMbN7WFdP3URxqg0GkwbGKMSsUWlvuunuMotr88I+y06mlVpQ504IkevbtyvdV8d1sLvtn6+u86nlR/vkiRdOHVwVF8rkbJSnLru9NG67KRhevGT3fr9um+0t6Raf/x4t/748W5lJNl1/NBsnTCst44ZkNmhAabVbq92F1dp56FKfXO4SjsPV2rn4UoVFFboUEXk3a/Tk+w6ZkCm/zHQ/+eg7BRCCYBui6CSIIGdkgNa3evHEb57cktBJbCDsmFIpdWesIGRwa6fKIPKtPr1Vj7fXawqd51SnC1/67y+eX9wI8NZYyNPNTezrBSnLj95uC6ZPkzrvjqsVZv2640thTpS6dY/vzygf355IHhuTpqzfq0Q/2qqyQ5bcA8Yr2Gopn7n3Uq3V4fq97A5WF4b3FQyEotFGpydotF5GRqVl64x/dI1Oi9Dg7JTuvRAXwCIFkElQRq3qNhtFjnt/h9AoYNpncGuH/9fVXUbWlQcNqsykuwqq6nTkcraYFDxr5vh7/rJj6LrR5IGZadoQFay9pZU65OdxTp5ZJ8Wzw/skvy/UwZ36cGaNqtF3z4qR98+Kke//NFY/Wd/mT7acUQffX1EWwvLtKe4Wocq3DoUYYZVW6S77BqSk6rBvVM0pHeqBvVO0cjcdI3MTWs1DAJAT9B1f4K04tFHH9WQIUOUlJSkKVOm6OOPP27x/JdeekmjR49WUlKSjjnmmJjXF2gpCbBbrXLamq6jkmQPn55c5zPk8fpaHKMiSb3T/ANXD4f8AD1U0bCGSv8oBtNK/nEqJ47wt6q8vfVAk8+H3u9jZ/xYn+8qkdNm1flTBjU596mnntL06dPVq1cv9erVSzNmzGj178cM7Darxg/M0qUnDdPTF03Wv64/TVtun6VXFnxbz1w0WfecdYx+MXOkrjpthC4/eZguPnGILp0+VP932ghdd/ooLZl9tB4+71v646Un6J+LTtYXt87Uxttm6pWrvq3fnn+sfjFrlP5ncr4m5me1GlKi/f4OeOGFF2SxWILrFqFtor3fJSUlmj9/vvr16yeXy6WRI0fqtddei1O1XV+09/vBBx/UqFGjlJycrPz8fC1cuFA1NS2vpAy/9957T7Nnz1b//v1lsVj0t7/9rdXXvPPOOzr22GPlcrk0YsQIPffcczGtsVv+yvbiiy9q0aJFWrZsmaZMmaIHH3xQs2bNUkFBgfr27dvk/A8++EDnnXee7r77bv3gBz/QypUr9bwntjU2Hkxrt1rksEdaR6W+6ydk1c8qt7dh5+QILSqSf5zK14cqwwbUBlpTollDJdR3j87Tnz7doze2FOrWHxwd7IJofL/nPfWeJGnGqF7qk950ps8777yj8847T9OmTVNSUpJ+/etfa+bMmdqyZYsGDOhamwWmuvwDWaXMuH3NaL+/A3bu3Klf/OIXmj59etxq7Q6ivd9ut1vf/e531bdvX7388ssaMGCAvvnmG2VlZcW/+C4o2vu9cuVK3XDDDVq+fLmmTZumbdu2ae7cubJYLFq6dGkC3kHXUllZqQkTJuinP/2pzjrrrFbP//rrr/X9739fV1xxhVasWKG1a9fqkksuUb9+/TRr1qyY1NgtW1SWLl2qSy+9VBdffLGOPvpoLVu2TCkpKVq+fHnE8x966CGdfvrpuvbaazVmzBj98pe/jHmNkbp+HGHTk8Nn/TjtVtnrg0G12xsc39Bci0qgu+dQSFDZXT8+Jdpun4DpR+Uo1WnT/tIafbGnJHg89H6n5Q7RgWR/K0ryrg8jXmfFihX62c9+pokTJ2r06NF6+umngysSo3XRfn9Lktfr1QUXXKDbb79dw4YNi2O1XV+093v58uU6cuSI/va3v+nEE0/UkCFDdPLJJ2vChAlxrrxrivZ+f/DBBzrxxBN1/vnna8iQIZo5c6bOO++8LtFKawbf+973dOedd+rHP/5xm85ftmyZhg4dqvvvv19jxozRggULdPbZZ+uBBx6IWY3dLqi43W6tX79eM2bMCB6zWq2aMWOG1q1bF/E169atCzs/HhoHFUfoOiphg2kbzgsdUNswRiV8+fyA3Ax/S8aBsobmz91Hol/sLVSSw6ZTR/t/o/n7Bv+y+I3v92/f/q+8PkO9avZr+ydvt+m6VVVV8ng8ys7OblddPUl7vr8l6Y477lDfvn01b968eJTZbbTnfv/jH//Q1KlTNX/+fOXm5mrcuHG666675A1Z4BGRted+T5s2TevXrw8Gkx07dui1117TGWecEZeae5pIPy9nzZrV4v8/HdXtgsqhQ4fk9XqbrGybm5urwsLCiK8pLCyM6Uq4kbhamp7sNZrM+pHCV6dtbYxKYB+ffSUNQaW9U5NDnTM5X5L05/V7VFlbF3a/tx8o158/2ytJOta5v9n73dj111+v/v37xz0sdkXt+f7+17/+pWeeeUZPPfVUPErsVtpzv3fs2KGXX35ZXq9Xr732mm655Rbdf//9uvPOO+NRcpfWnvt9/vnn64477tC3v/1tORwODR8+XKeccopuvPHGeJTc40T6eZmbm6uysjJVV0fa0LTjul1Q6SqatqhYgi0qXp8RnN3jCpnGnByy6FtrY1T6ZSZJkgrLGr5xgou9Zbd/Z+jpI3I0NCdV5bV1WvHRN8HjPsPQzX/bLK/P0Iwxueqt8jZd75577tELL7ygv/71r0pKSmp3XYisvLxcF154oZ566inl5OQkupwewefzqW/fvnryySc1adIknXvuubrpppu0bNmyRJfWLb3zzju666679Nhjj+mzzz7TX/7yF61atSouXfiIj243mDYnJ0c2m01FRUVhx4uKipSXF3k9j7y8vCbnx1rjFhW71SpHyLHyGk+T80KX0S+vDyqpzQSVvPqgsr+k87p+JMlqtejKU4brupc36uG12zVj1Amy2Wx6/rND+nCPTUkOq5bMPlo3v9X8/Q647777dM899+if//ynxo8f3+6aepJov7+/+uor7dy5U7Nnzw4e8/nqd+y221VQUKDhw4fHtugurD3/n/Tr108Oh0M2W0Nr6JgxY1RYWCi32y2nkw0fm9Oe+33LLbfowgsv1CWXXCJJOuaYY1RZWanLLrtMN910k6xWfh/vTJF+XhYVFSkjI0PJye3/2dKSbvc36HQ6NWnSpLCBmYGBmlOnTo34mqlTp8Z9IGfjWT82qyW4e7LUsE5KaNdPYBn9ardXFTXN7/UjSf3ru372l9bIMPy7NAcG0w7p3bHNAc8+dqAm5meporZOFz67XsMvuler9/hru+NH4zQgK6nF+y1J9957r375y19q9erVmjx5cofq6Umi/f4ePXq0Nm3apA0bNgQfP/zhD3Xqqadqw4YNys/Pj2f5XU57/j858cQTtX379mAglKRt27apX79+hJRWtOd+V1VVNQkjgZBoGEbsiu2hIv28XLNmTYv/33dUtwsqkrRo0SI99dRT+t3vfqcvv/xSV155pSorK3XxxRdLkn7yk59o8eLFwfOvvvpqrV69Wvfff7+2bt0alw2YXI02uQvt+pEkX/2/r7AWFWfTwbTpzYxRCbSoVHu8Kq32aNeRSnl9hlKctmC3UHtZrRYt+99JGpSdor0l1artM0qSdFqfSh2TWtHq/f71r3+tW265RcuXL9eQIUNUWFiowsJCVVRUdKiuniKa7++kpCSNGzcu7JGVlaX09HSNGzeOH5xtEO3/J1deeaWOHDmiq6++Wtu2bdOqVat01113af78+Yl6C11KtPd79uzZevzxx/XCCy/o66+/1po1a3TLLbdo9uzZYa1aiKyioiL4S4zkn368YcMG7drl3wJl8eLF+slPfhI8/4orrtCOHTt03XXXaevWrXrsscf0pz/9SQsXLoxZjd2u60eSzj33XB08eFC33nqrCgsLNXHiRK1evTo4AGjXrl1hCXzatGlauXKlbr75Zt1444066qijpB8cF9Mam3T92KyyWPytKh5vw28Bobv2Jjv8f11Vbm9wMG1zXT9JDpuyU506UunWvpIa7arv9hneJ61T9oXJy0zSa1dP118+26NDFW4d2fSOnn/oV1pxY+v3+/HHH5fb7dbZZ58dds0lS5YkfJfOriDa7290TLT3Oz8/X2+88YYWLlyo8ePHa8CAAbr66qt1/fXXJ+otdCnR3u+bb75ZFotFN998s/bu3as+ffpo9uzZ+tWvfpWot9ClfPrppzr11FODzxctWiRJuuiii/Tcc89p//79wdAiSUOHDtWqVau0cOFCPfTQQxo4cKCefvrpmK2hIkkWo4u3jZWVlSkzM1OlpaXKyOi8HXqH3LCq064VyTXfHan712wLPr/qtBG6ZuYoHX3r6uDGg5K08baZykjyT0G+6o+f65Uv9unGM0brrte2SpI+u+W7YXv5hPrRo//WF7tL9Oj5x2rn4Ur95o0CnTmxvx78f9+K4TsDAKDz8GtXgjRZ8M1qjXg8tOUlpb51ZX9pwwDZjGa6fiTpqL5pkqT/HijXtiL/LJwR9ccAAOgKumXXT1fQtOvH3x3jaDTINnTcSqCbJ7AeSnqSXXZb81lzZG59UCmq0H/2l0mSjhmY1bHCAQCII4JKgjjtTQfTSuHBxGW3ho0nyUrxdwHtOlwV9rw5R+WmS5I+3HFYh+uX0p8wMH570gAA0FF0/SRI4y4eW4Sun8atLoFgsvNwpf95csszNr6VnyWrRcGQMqJvmrJSmOUBAOg6CCoJErpmSujz0OONpzBnJvuDSmB5/dZaVLJSnJo0uFfw+eljW16ADQAAsyGoxJE1JJvYrQ27IQeeS621qIS3hgSCS0sWzhipZIdNA7KSddG0Ie2oGgCAxGGMShzZbVa561tDbFb/arR19Su7BUKLw9ZCUGkUTFprUZGkaSNy9PFN31Gq0y6rtePrpwAAEE+0qMRR6EBZq8US3qISYdaPq9GA28bBpG9621aYTU9yEFIAAF0SQSWO7LbwYBI6tTjwcWgrSuN9fBoPns3LYLdhAED3RlCJI7u1hRaVCF0/Ka7wFpWMZHtwY0JJ6pvhilWpAACYAkEljkJn9NisFtlCgkqgWyi0eyg0lEiSxWLRwF4N22gP7JUSq1IBADAFgkochXb92Bq1qARm+yQ5QoNK07HOuSHdPUN6E1QAAN0bQSWOHCFdPzarRTZb06CSntQwYDbV2XSL8nnfHiqb1aILTxjc4vL5AAB0B0xPjqPQ8Sc2qyVszEogqKSFbDKY4mr613PKqL7acvssJTmahhgAALobfiWPo9CuH2szY1TSQ4JKpBYVSYQUAECPQVCJo9CumsZjVFwRun5CPwYAoCciqMSRs9Gsn9AWlsDibukh3T3ZqWwgCADo2QgqcWRvNJg2KWTl2cAYldDVZ3sTVAAAPRxBJY4cIavO2q0WuRxNB9OOyksPHhvaJzV+xQEAYELM+okjR8iYFIfNKluEWT/9MpP1yzPHyW61qF9mcpNrAADQkxBU4sjWaBPCkCEqYSvSXnjC4HiWBQCAadH1E0ehXT/ORou1hS6vDwAA/AgqcWS1hHf91PmM4HOLhaACAEBjBJUEsdss8oYEFQAA0BRBJY5Ce3ccNqvyMpOaPxkAADCYNp5C9/px2Ky68YwxOlBWq3OPy09gVQAAmBdBJY7sjVamzUlz6Q+XTElgRQAAmBtdP3FkZcAsAABRIajEETkFAIDoEFTiyFPHLB8AAKJBUImjQb1TEl0CAABdCoNp4+j/HZevtV8W6cffGpDoUgAA6BIIKnHUO82lv/zsxESXAQBAl0HXDwAAMC2CCgAAMC2CCgAAMC2CCgAAMC0G08bYaaP7qneqUz+c2D/RpQAA0OUQVGIsPcmu35wzIdFlAADQJdH1AwAATIugEmMOG7cYAID24qdojDnt3GIAANqLn6Ix9t0xuYkuAQCALovBtDEyZWi2rjhluE4Z2SfRpQAA0GXFpEVl586dmjdvnoYOHark5GQNHz5cS5YskdvtDjtv48aNmj59upKSkpSfn6977703FuUkRP+sZJ06qq8sFkuiSwEAoMuKSYvK1q1b5fP59MQTT2jEiBHavHmzLr30UlVWVuq+++6TJJWVlWnmzJmaMWOGli1bpk2bNumnP/2psrKydNlll8WirLiyElAAAOiwmASV008/Xaeffnrw+bBhw1RQUKDHH388GFRWrFght9ut5cuXy+l0auzYsdqwYYOWLl3aLYIKk30AAOi4uP04LS0tVXZ2dvD5unXrdNJJJ8npdAaPzZo1SwUFBSouLm72OrW1tSorKwt7mFHvNFeiSwAAoMuLS1DZvn27HnnkEV1++eXBY4WFhcrNDZ8RE3heWFjY7LXuvvtuZWZmBh/5+fmxKboDTh+bp8tPGpboMgAA6PKiCio33HCDLBZLi4+tW7eGvWbv3r06/fTTdc455+jSSy/tcMGLFy9WaWlp8LF79+4OX7OzLbtwkrJSnK2fCAAAWhTVGJVrrrlGc+fObfGcYcMaWhL27dunU089VdOmTdOTTz4Zdl5eXp6KiorCjgWe5+XlNXt9l8sll4tuFQAAeoKogkqfPn3Up0/b1gXZu3evTj31VE2aNEnPPvusrNbwxpupU6fqpptuksfjkcPhkCStWbNGo0aNUq9evaIpCwAAdFMxGaOyd+9enXLKKRo0aJDuu+8+HTx4UIWFhWFjT84//3w5nU7NmzdPW7Zs0YsvvqiHHnpIixYtikVJAACgC4rJ9OQ1a9Zo+/bt2r59uwYOHBj2OcMwJEmZmZl68803NX/+fE2aNEk5OTm69dZbu/zU5DnHDmz9JAAA0CYWI5AcuqiysjJlZmaqtLRUGRkZnXbdITesivo1/3vCIC2ZPZYdkwEA6CT8RO1EY/plEFIAAOhE/FTtRCybDwBA5yKodKKu3YkGAID5EFQ6UXpSTMYmAwDQYxFUOsmPJvbX98Y1v1AdAACIHkGlk9x/zgTZGUgLAECn4idrJ2EgLQAAnY+g0knIKQAAdD6CSiexkFQAAOh0BJVOcO2sUYkuAQCAbomg0gnmnzoi0SUAANAtEVQAAIBpEVQAAIBpEVQAAIBpEVQAAIBpEVQ66IslMxNdAgAA3RZBpYMykx2JLgEAgG6LoAIAAEyLoAIAAEyLoAIAAEyLoNIBt80+OtElAADQrRFUOmBkbnqiSwAAoFsjqHQEGyYDABBTBJUOmDw4O9ElAADQrRFU2mnCwEw57dw+AABiiZ+07WSx0O8DAECsEVTaiZwCAEDsEVQAAIBpEVTaiQYVAABij6DSTnMmDUx0CQAAdHsElXY677hBiS4BAIBuj6DSTlYrnT8AAMQaQQUAAJgWQQUAAJgWQaUdbHT7AAAQFwSVdvjLldMSXQIAAD0CQaUdUl32RJcAAECPQFBphyQHtw0AgHjgJ247DOyVkugSAADoEQgqAADAtAgqAADAtAgqAADAtAgqAADAtAgqAADAtAgqUbpo6uBElwAAQI9BUInSVd85KtElAADQYxBUosQuPwAAxA9BJUrZqc5ElwAAQI9BUImSxUKbCgAA8UJQAQAApkVQAQAApkVQicLg3mxGCABAPBFUorBk9tGJLgEAgB6FoBIFm5XbBQBAPPGTFwAAmBZBBQAAmBZBJQrD+6QmugQAAHoUgkoUBvZi1g8AAPFEUAEAAKYV86BSW1uriRMnymKxaMOGDWGf27hxo6ZPn66kpCTl5+fr3nvvjXU5AACgC4l5ULnuuuvUv3//JsfLyso0c+ZMDR48WOvXr9dvfvMb3XbbbXryySdjXRIAAOgi7LG8+Ouvv64333xTf/7zn/X666+HfW7FihVyu91avny5nE6nxo4dqw0bNmjp0qW67LLLYlkWAADoImLWolJUVKRLL71Uzz//vFJSmg5CXbdunU466SQ5nc7gsVmzZqmgoEDFxcXNXre2tlZlZWVhDwAA0D3FJKgYhqG5c+fqiiuu0OTJkyOeU1hYqNzc3LBjgeeFhYXNXvvuu+9WZmZm8JGfn995hbfgT5dPjcvXAQAADaIKKjfccIMsFkuLj61bt+qRRx5ReXm5Fi9e3OkFL168WKWlpcHH7t27O/1rRDIxPysuXwcAADSIaozKNddco7lz57Z4zrBhw/TWW29p3bp1crlcYZ+bPHmyLrjgAv3ud79TXl6eioqKwj4feJ6Xl9fs9V0uV5PrAgCA7imqoNKnTx/16dOn1fMefvhh3XnnncHn+/bt06xZs/Tiiy9qypQpkqSpU6fqpptuksfjkcPhkCStWbNGo0aNUq9evaIpCwAAdFMxmfUzaNCgsOdpaWmSpOHDh2vgwIGSpPPPP1+333675s2bp+uvv16bN2/WQw89pAceeCAWJQEAgC4optOTW5KZmak333xT8+fP16RJk5STk6Nbb72VqckAACAoLkFlyJAhMgyjyfHx48fr/fffj0cJHWaxJLoCAAB6Hvb6AQAApkVQAQAApkVQaSMrfT8AAMQdQaWNbFaCCgAA8UZQAQAApkVQAQAApkVQaQM73T4AACQEQaUNGEcLAEBiEFQAAIBpEVQAAIBpEVTaYHiftESXAABAj0RQaYPFZ4xJdAkAAPRIBJVmXP2do4IfZyQlbJNpAAB6NIJKM0JXomX5fAAAEoOg0gbHDMhMdAkAAPRIBJU2sLLgGwAACUFQAQAApkVQAQAApkVQAQAApkVQAQAApkVQacbp4/IkSQN7JSe4EgAAei6LYRhGoovoiLKyMmVmZqq0tFQZGRmdeu19JdXKTnUqyWHr1OsCAIC2YcnVFvTPojUFAIBEousHAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYFkEFAACYlsUwDCPRRXSEYRgqLy9Xenq6LBZLossBAACdqMsHFQAA0H3R9QMAAEyLoAIAAEyLoAIAAEyLoAIAAEyLoAIAAEyLoAIAAEyLoAIAAEzr/wMCQK9ctjTCRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = lambda x: cos(log(x)/x)/x\n",
    "x = np.linspace(0.02,1,10000)\n",
    "y = f(x)\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(1, 1, 1)\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.xaxis.set_ticks_position('bottom')\n",
    "ax.yaxis.set_ticks_position('left')\n",
    "plt.plot(x,y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lb9DXWvg7D2H"
   },
   "source": [
    "As we can see, the integrand is very ill-behaved and oscillates a lot - therefore typical numerical integration algorithms may not work that well. Nonetheless, we will try that as a first attempt."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zRrlwv9D7D2H"
   },
   "source": [
    "## Method 1: Integrating as-is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "pbu8fKkM7D2I",
    "outputId": "d2aa2c4b-a518-4f40-f7e8-506fd2dfeb4e"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_88714/4081173682.py:2: IntegrationWarning: The integral is probably divergent, or slowly convergent.\n",
      "  integrate.quad(f,0,1)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(1.2716741383488144, 1.9276166879370535)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy import integrate\n",
    "integrate.quad(f,0,1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DAAYNIfy7D2L"
   },
   "source": [
    "As expected, this does not work as the integral is not well-behaved (note the given warning and that the error estimate is larger than the value it got for the integral itself).\n",
    "\n",
    "However, since the integral is well behaved between its zeros, we can have points $ 0 = x_0 < x_1 < \\cdots < x_{N-1}  < x_N = 1 $ where $f(x_i) = 0 $ for $0 < i < N$ .\n",
    "\n",
    "We therefore have:\n",
    "\n",
    "$$ \\int_{0}^1 f(x) dx = \\sum_{i=0}^{N-1} \\left( \\int_{x_i}^{x_{i+1}} f(x) dx \\right) $$\n",
    " \n",
    "where the integrand over each $[x_{i}, x_{i+1}]$ is well-behaved. We have:\n",
    "\n",
    "\n",
    "$$\\begin{align*}\n",
    "f(x) = 0 &\\implies \\cos(x^{-1} \\log{x}) = 0 \\\\\n",
    "&\\implies x^{-1} \\log{x} = \\pi\\left(n-\\frac{1}{2}\\right), n \\in \\mathbb{Z}_{<0} \\\\\n",
    "&\\implies x = e^{-W\\left(\\frac{\\pi}{2} - n\\pi\\right)}\n",
    "\\end{align*}$$\n",
    "\n",
    "where $ W(x) $ is the [Lambert W function](https://en.wikipedia.org/wiki/Lambert_W_function). An implementation of this is below (note that we only use floating point precision):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "id": "zBdIqgIv7D2L",
    "outputId": "298df7f8-0acf-4f7c-8b45-c28908938d9d"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_88714/502815136.py:8: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0\n",
      "Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`\n",
      "  intervals = [integrate.quad(f, nth_zero(i), nth_zero(i+1)) for i in tqdm(range(0, -100000, -1)) ]\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "326250febd69459eb4d5926239c7ca24",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "  0%|          | 0/100000 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_88714/502815136.py:5: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return exp(-lambertw(pi/2-n*pi)).real\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.32336452967614154, 9.410465197082566e-13)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.special import lambertw\n",
    "from math import pi, exp\n",
    "\n",
    "def nth_zero(n):\n",
    "    return exp(-lambertw(pi/2-n*pi)).real\n",
    "\n",
    "\n",
    "intervals = [integrate.quad(f, nth_zero(i), nth_zero(i+1)) for i in tqdm(range(0, -100000, -1)) ]\n",
    "integral_val = sum([i[0] for i in intervals])\n",
    "error = sum([i[1] for i in intervals])\n",
    "\n",
    "integral_val, error"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mgItcC3V7D2R"
   },
   "source": [
    "Note that there is a problem with this approach: namely, the closer that we get to 0, the more the integrand diverges - therefore, the more accurate value we want, the more subintervals we need to integrate over. We also have to worry about the error terms of calculating the functions of each of the values, so it would be nicer if we calculated just one integral instead.\n",
    "\n",
    "**EXERCISE**: Is there any way we can bound the error with this method? We have:\n",
    "\n",
    "$$ \\int_0^1 f(x) dx \\approx \\int_\\epsilon^1 f(x) dx $$ however, what bounds can we place on the error term $\\int_0^\\epsilon f(x) dx  $ ?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ZsJeKC1-7D2S"
   },
   "source": [
    "## Method 2: Complex integration"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "LK2OrJzi7D2T"
   },
   "source": [
    "We have that\n",
    "\n",
    "$$ \\begin{align*} \n",
    "\\int_0^1 \\frac{\\cos(x^{-1} \\log{x})}{x} dx &= \\int_0^1 \\Re\\left(\\frac{e^{i(\\log{z}/z)}}{z}\\right) \\\\\n",
    "& = \\int_0^1 \\Re(z^{-1+i/z}) \\\\\n",
    "&= \\Re\\left(\\int_0^1 z^{-1+i/z}\\right) \\\\\n",
    "&= \\Re\\left(\\int_C z^{-1+i/z}\\right) \\\\\n",
    "\\\\ \n",
    "\\end{align*}\n",
    "$$ \n",
    "\n",
    "where $C$ is a contour from $z = 0 $ to $z = 1 $ such that the integrand is analytic on the whole contour. \n",
    "\n",
    "We use a semicircle oriented clockwise centered at $z = \\frac{1}{2}  $ on the complex plane with a radius of $\\frac{1}{2} $ . This can be parametrized with the map \n",
    "$ r: [0, \\pi] \\to \\mathbb{C} \\;,  r(\\theta) = \\frac{1}{2} - \\frac{1}{2}e^{-i\\theta}$\n",
    "\n",
    "Substituting this in, we have our integral is\n",
    "\n",
    "$$ \\Re\\left(\\int_C z^{-1+i/z}\\right) = \\Re \\left( \\int_0^\\pi   \\frac{ie^{-it}}{1-e^{-it}} \\left(\\frac{1-e^{-it}}{2}\\right)^{\\frac{2i}{1-e^{-it}}} dt  \\right)   $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "id": "GcACwRMw7D2U",
    "outputId": "5a0ac6f5-9b90-4ffd-fe6e-320878cf0684",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.86 ms, sys: 240 µs, total: 4.1 ms\n",
      "Wall time: 5.6 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/husnain/Desktop/RandomProjects/100-Digit-Challenge/venv/lib/python3.10/site-packages/scipy/integrate/_quadpack_py.py:523: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return _quadpack._qagse(func,a,b,args,full_output,epsabs,epsrel,limit)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.3233674316777749, 5.940167106104359e-09)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from numpy import exp\n",
    "\n",
    "ex = lambda t: exp(1j*t)\n",
    "integrand = lambda t: 1j*ex(-t)/(1-ex(-t)) * (0.5*(1-ex(-t)))**(2j/(1-ex(-t)))\n",
    "\n",
    "%time integrate.quad(integrand,0,pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zZL0ec247D2X"
   },
   "source": [
    "Note that this is much faster than the previous method as there is only one integral to evaluate - we also get the 10 digits that we need. We can get further digits by using an extended precision library like `mpmath`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "id": "KjqEdruo7D2Y",
    "outputId": "70b4520e-1ecf-40a4-e217-a19a08f2ab4d",
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 9min 21s, sys: 1.21 s, total: 9min 22s\n",
      "Wall time: 9min 33s\n",
      "0.3233674316777787613993700879521704466510466257254696616810364434317903372106728944319303704641024513802805273312171512843329551003343518853814345933849959100496203187230243755305850326867826084814143621008284544722325433197423289562208383036051806138170228875927796223889414194808897786277656720091636325480930062612368830709515617259626227109122796181848626492520354425805878303512878922713506176972504363427324278062698195459935937934948703377003993887895844494193853252470105126349141166869280842967050545084228387406907420047518022017247554109230583360447095045817596946489788630167768426091736144508596750494011596931816971587028572107158052888912310396674019547711228705979109046809875599383142853899314897723348536769454368361557868526285246345537463616787033605541989910530337300290559520490910138131952127532519967705188100783651118877283568747500960917691186552933459664209305926699509607663216838251160999002569152976749531082493203607852959693352367456372093026716239528268235600132382543 1.0e-1761\n"
     ]
    }
   ],
   "source": [
    "from mpmath import mp, mpf, mpc, exp, quad, pi\n",
    "mp.dps = 1000\n",
    "ex = lambda t: exp(mpc(0,1)*t)\n",
    "integrand = lambda t: mpc(0,1)*ex(-t)/(1-ex(-t)) * (mpf(\"0.5\")*(1-ex(-t)))**(mpc(0,2)/(1-ex(-t)))\n",
    "%time result, error = quad(integrand, [0,pi], error=True)\n",
    "print(result.real, error)"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Problem 1.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}