{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Daudet-Flajolet-Vallée for Gauss map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$[(z-x_0)^j]G_{s,\\mathcal A}[(z-x_0)^i] = (-1)^i \\sum_{u=0}^j \\binom{j}{u} (-x_0)^{j-u} \\binom{\\color{red}{2}s+u+i-1}{i} \\zeta_{\\mathcal A}(\\color{red}{2}s+u+i,x_0)$$\n",
    "\n",
    "où\n",
    "\n",
    "$$\\zeta_{\\mathcal A}(s,x) := \\sum_{m \\in \\mathcal A} \\frac{1}{(m+x)^s}.$$\n",
    "\n",
    "Ici on considère $\\mathcal A = \\mathbb N^*$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### La fonction zeta de Hurwitz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hurwitz_zeta?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "def h_zeta(s,x):\n",
    "    return hurwitz_zeta(s,x)-x^(-s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fonction de définition de la matrice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = lambda i,j : (-1)^i * sum(binomial(j,u) * \\\n",
    "                              (-x0)^(j-u) * \\\n",
    "                              binomial(2*s+u+i-1,i) * \\\n",
    "                              h_zeta(2*s+u+i,x0) \\\n",
    "                              for u in range(j+1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Définition de quelques constantes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "s = 1\n",
    "x0= 1\n",
    "prec = 256"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Complex Ball Field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 464 ms, sys: 1 ms, total: 465 ms\n",
      "Wall time: 463 ms\n",
      "CPU times: user 11.9 ms, sys: 1 ms, total: 12.9 ms\n",
      "Wall time: 12.7 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[-0.3031723850818134338319718557923656049497010097838671109597337762715 +/- 7.74e-68]"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10\n",
    "CBF = ComplexBallField(prec)\n",
    "%time M = matrix(N, N, lambda i,j : CBF(f(i,j)))\n",
    "%time sorted(map(real_part,M.eigenvalues()))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 11.5 s, sys: 0 ns, total: 11.5 s\n",
      "Wall time: 11.5 s\n",
      "CPU times: user 229 ms, sys: 0 ns, total: 229 ms\n",
      "Wall time: 228 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[-0.30366300271208809189566779332484408721221959665473021 +/- 6.41e-54]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 30\n",
    "%time M = matrix(N, N, lambda i,j : CBF(f(i,j)))\n",
    "%time sorted(map(real_part,M.eigenvalues()))[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 26.6 s, sys: 3.97 ms, total: 26.6 s\n",
      "Wall time: 26.6 s\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "unable to certify the eigenvalues",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<timed eval>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.7/site-packages/sage/misc/superseded.py\u001b[0m in \u001b[0;36mwrapper\u001b[0;34m(*args, **kwds)\u001b[0m\n\u001b[1;32m    281\u001b[0m                             self.stacklevel)\n\u001b[1;32m    282\u001b[0m                 \u001b[0mwrapper\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_already_issued\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 283\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwds\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    284\u001b[0m         \u001b[0mwrapper\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_already_issued\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    285\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/usr/lib/python3.7/site-packages/sage/matrix/matrix_complex_ball_dense.pyx\u001b[0m in \u001b[0;36msage.matrix.matrix_complex_ball_dense.Matrix_complex_ball_dense.eigenvalues (/var/tmp/portage/sci-mathematics/sage-9.2/work/sage-9.2/src-python3_7/build/cythonized/sage/matrix/matrix_complex_ball_dense.c:8988)\u001b[0;34m()\u001b[0m\n\u001b[1;32m    719\u001b[0m             \u001b[0meigval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_acb_vec_init\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mn\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    720\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0macb_mat_eig_multiple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0meigval\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mvalue\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meigval_approx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0meigvec_approx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprec\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 721\u001b[0;31m                 \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"unable to certify the eigenvalues\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    722\u001b[0m             \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_acb_vec_to_list\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0meigval\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_parent\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_base\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    723\u001b[0m         \u001b[0;32mfinally\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: unable to certify the eigenvalues"
     ]
    }
   ],
   "source": [
    "N = 40\n",
    "%time M = matrix(N, N, lambda i,j : CBF(f(i,j)))\n",
    "%time sorted(map(real_part,M.eigenvalues()))[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Flottants\n",
    "\n",
    "Le corps CBF ne parvient pas à certifier les valeurs propres de trop grandes matrices.\n",
    "\n",
    "Nous changeons pour des flottants de précision arbitraires"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min 32s, sys: 163 ms, total: 2min 33s\n",
      "Wall time: 2min 33s\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib/python3.7/site-packages/sage/repl/ipython_kernel/__main__.py:1: UserWarning: Using generic algorithm for an inexact ring, which will probably give incorrect results due to numerical precision issues.\n",
      "  from ipykernel.kernelapp import IPKernelApp\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6.76 s, sys: 42 ms, total: 6.81 s\n",
      "Wall time: 6.76 s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.3036630028987326585078485894353481019684291685283653062335525685201680545065"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 60\n",
    "%time M = matrix(N, N, lambda i,j : f(i,j).n(prec))\n",
    "%time sorted(M.eigenvalues())[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "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.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}