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