{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Scaled Bernoulli numbers\n",
    "\n",
    "This SageMath notebook is a supplement to: Marc Mezzarobba, *Rounding Error Analysis of Linear Recurrences Using Generating Series*.\n",
    "\n",
    "It uses F. Johansson's interval arithmetic library Arb (http://arblib.org/) to derive some numerical estimates used in the section on scaled Bernoulli numbers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "umax = 2**(-16)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rigorous check that Rouché's theorem applies in the proof of Lemma 19"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = RBF(1 + umax)\n",
    "rho = RBF(62/10)\n",
    "(a**2 - 1)*(1 + (a**2*rho).cosh()), (sin(rho)/rho).abs()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Estimates used at the end of the proof of Lemma 20"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = RBF(umax).union(0)\n",
    "K = 2*(RBF(pi).cosh() - 1)\n",
    "phi = K*u # interval containing all possible values of K*u\n",
    "alpha = RBF(pi)/(1 + phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code for computing R' exported from the Maple expression.\n",
    "#\n",
    "# One could likely get tighter bounds by computing the derivative (and\n",
    "# the second derivative) by composition of power series, but this\n",
    "# would require some work since alpha depends on u.\n",
    "\n",
    "def R_u():\n",
    "    t1 = alpha\n",
    "    t3 = 1+u\n",
    "    t4 = t1*t3\n",
    "    t5 = t3^2\n",
    "    t6 = t5*t1\n",
    "    t7 = t6.cosh()\n",
    "    t9 = - t1*u+t4*t7-t1\n",
    "    t10 = u^2\n",
    "    t11 = 2*u\n",
    "    t13 = t1.cosh()\n",
    "    t14 = t1.cos()\n",
    "    t15 = -t5* t7+t10+t11+t13+t14\n",
    "    t16 = 1/t15\n",
    "    t17 = 2*t16*t9\n",
    "    t20 = t4.sinh()\n",
    "    t23 = t4.cosh() \n",
    "    t27 = 2*t9*t1\n",
    "    t28 = t1.sinh()\n",
    "    t32 = t1.sin()\n",
    "    t36 = 1/t1\n",
    "    t38 = -1/t15\n",
    "    t44 = -t1*t13-t1*t14+t1*t23+2*t32\n",
    "    t45 = t1^2\n",
    "    t52 = t15^2\n",
    "    t60 = t6.sinh()\n",
    "    t67 = t38 *t36*(t1*t20*(t17*t3+t1)+2*t16*t9*t23-t13*t17-t28*t16*t27+t14*t17+t32*t16*t27)-2*t16*t9*t38/t45*t44-(2*t7*t3+t60*(t17*t5+2*t4)*t5-t11-2+t32*t17-t28*t17)/t52* t36*t44\n",
    "    return t67.above_abs()\n",
    "\n",
    "R_u().above_abs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note that the numerator is 2*K*(2+phi), not 2*phi*(2+phi),\n",
    "# since we want a bound of the form C*u.\n",
    "\n",
    "((2*K*(2+phi)/(2*alpha)**2)).above_abs()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code exported from the Maple expression, as above.\n",
    "\n",
    "def g_u(w):\n",
    "    t1 = 1+u\n",
    "    t2 = t1^2\n",
    "    t3 = w*t2\n",
    "    t4 = t3.cosh()\n",
    "    t6 = t1*t4-u-1\n",
    "    t7 = w.tan()\n",
    "    t10 = 2-1/t7*w\n",
    "    t12 = w*t1\n",
    "    t13 = t12.sinh()\n",
    "    t16 = 1/w\n",
    "    t17 = w.sin()\n",
    "    t19 = t3.sinh()\n",
    "    t20 = w.sinh()\n",
    "    t22 = (t19-t20)*t16\n",
    "    t23 = t16*t17+t2-t22-1\n",
    "    t28 = t12.cosh()\n",
    "    t29 = w.cosh()\n",
    "    t31 = t23^2\n",
    "    t35 = alpha\n",
    "    t37 = t35*t1\n",
    "    t38 = t2*t35\n",
    "    t39 = cosh (t38)\n",
    "    t41 = -t35*u+t37*t39-t35\n",
    "    t42 = u^2\n",
    "    t43 = 2*u\n",
    "    t45 = t35.cosh()\n",
    "    t46 = t35.cos()\n",
    "    t47 = -t2*t39+t42+t43+t45+t46\n",
    "    t48 = 1/t47\n",
    "    t49 = 2*t48*t41\n",
    "    t52 = t37.sinh()\n",
    "    t55 = t37.cosh()\n",
    "    t59 = 2*t41*t35\n",
    "    t60 = t35.sinh()\n",
    "    t64 = t35.sin()\n",
    "    t68 = 1/t35\n",
    "    t70 = -1/t47\n",
    "    t71 = w^2\n",
    "    t72 = t35^2\n",
    "    t73 = 1/t72\n",
    "    t75 = -t71*t73+ 1\n",
    "    t76 = 1/t75\n",
    "    t84 = -t35*t45-t35*t46+t35*t55+2*t64\n",
    "    t92 = t47^2\n",
    "    t101 =  t38.sinh()\n",
    "    t109 = t72^2\n",
    "    t113 = t75^2\n",
    "    t119 = 1/t23*(2*t10*t6+t13*w)+2*t6/t31*(t10*( t22-t2+1)+t28-t29)-2*t76*t70*t68*(t35*t52*(t1*t49+t35)+2*t48*t41*t55-t45*t49- t60*t48*t59+t46*t49+t64*t48*t59)+4*t48*t41*t76*t70*t73*t84+2*(2*t39*t1+t101*(t2 *t49+2*t37)*t2-t43-2+t64*t49-t60*t49)*t76/t92*t68*t84+4*t49*t71/t113*t70/t109* t84\n",
    "    return t119\n",
    "\n",
    "# To bound |∂g/∂u| on the circle |w| = λπ, we cover the unit circle by\n",
    "# n squares centered at exp(2*i*π/n), scale everything by λπ, and run\n",
    "# the above implementation of ∂g/∂u in complex interval arithmetic on\n",
    "# each of these squares.\n",
    "\n",
    "def bound_g_u(n):\n",
    "    k = 1\n",
    "    lambda_pi = RBF(3/2).sqrt()*RBF(pi)\n",
    "    itheta0 = CBF(0, 2*pi)/n\n",
    "    halfside = (RBF(pi)/n).sin()\n",
    "    bound = RBF.zero()\n",
    "    for k in range(n):\n",
    "        mid = (k*itheta0).exp()\n",
    "        # mid.add_error(r) adds r to the radii of both the real and\n",
    "        # imaginary part of mid\n",
    "        w = lambda_pi*mid.add_error(halfside)\n",
    "        val = g_u(w)\n",
    "        bound = bound.max(val.above_abs())\n",
    "    return bound\n",
    "\n",
    "bound_g_u(65536)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The constants in the statement"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "2*72+747"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "892/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The constants in the corollary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "RBF(1).exp()*11/10, RBF(1).exp()*446"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}