Floating point errors in the solution of the quadratic equation
===============================================================
In this notebook I attempt to derive formula for the errors in
the solution of the quadratic equation $A x^2 + B x + C = 0$. 
I respect two kind of errors: Propagated uncertainty in $A, B,$
and $C$. And roundoff errors in the solution formula due to
floating point computations.


```
import scipy
import numpy as np
import sympy
#from sympy import *
from IPython.display import display, Math, Latex
def prettylatex(x):
    display(Math(sympy.latex(x)))
import sys
import os
import interval
import math
```


```
def itostring(i):
    assert len(i) == 1
    return u'{}+/-{} ⊂ {}'.format(i.midpoint[0].inf,0.5*(i[0].sup-i[0].inf),i)

def eps_half_interval_around(x, eps):
    return interval.interval([x*(1-0.5*eps), x*(1+0.5*eps)])
```


```
def gamma_t(dtype):
    # Get machine epsilon. (Isn't the machine epsilon actually half of that? ... Well not that it's a big deal ...)
    eps = np.finfo(np.float16).eps.astype(np.float64)
    # Then make the gamma function as defined in PBRT pg.216
    def gamma(n):
        return n*eps/(1-n*eps)
    return gamma

def typedlike(x, t):
    """ return x cast to same type as t. t.dtype must give the numpy type."""
    return np.array(x, dtype=t.dtype)
```


```
def errorformula1(A, B, C, rA, rB, rC):
    xi = typedlike(1 if B < 0 else -1,A)
    one = typedlike(1,A)
    four = typedlike(4, A)
    two  = typedlike(2, A)
    three = typedlike(3, A)
    eps  = np.finfo(A.dtype).eps
    Q = B*B - A*C
    sQ = np.sqrt(Q)
    E1 = eps*(np.abs((B - sQ*xi)/A) + np.abs(sQ/A/two) + np.abs(C/sQ) + np.abs(B*B/four/A/sQ) + np.abs(sQ/A/four))
    E2 = rA*np.abs(A*(xi*C/A/sQ + (B-xi*sQ)/(two*A*A))) + rB/two*np.abs(B/A*(B*xi/sQ-one)) + rC*np.abs(C/sQ)
    return E1+E2

def errorformula2(A, B, C, rA, rB, rC):
    xi = typedlike(1 if B < 0 else -1,A)
    one = typedlike(1,A)
    four = typedlike(4, A)
    two  = typedlike(2, A)
    three = typedlike(3, A)
    eps  = np.finfo(A.dtype).eps
    Q = B*B - A*C
    sQ = np.sqrt(Q)
    E1_den1 = B - xi*sQ
    E1_den2 = E1_den1*E1_den1
    E1 = eps*(four*np.abs(C/E1_den1) + two*np.abs(C*sQ/E1_den2) + np.abs(four*A*C*C/sQ/E1_den2) + np.abs(B*B*C/sQ/E1_den2) + np.abs(C*sQ/E1_den2))
    E2 = four*rA*np.abs(A*C*C/sQ/E1_den2) + two*rB*np.abs((B*C*(B*xi-sQ))/sQ/E1_den2) + rC*np.abs(C*(four*A*C*xi/sQ/E1_den2 - two/E1_den1))
    return E1+E2
```


```
sqrt = np.sqrt

def quadratic(A, B, C):
    Q = B*B - typedlike(4, A)*A*C
    sQ = sqrt(Q)
    sol1 = (-B + sQ)/(typedlike(2, A)*A)
    sol2 = (-B - sQ)/(typedlike(2, A)*A)
    return sol1,sol2

def better_quadratic(A, B, C, rA, rB, rC):
    # from PBRT pg. 1080
    D = B*B - typedlike(4, A)*A*C
    sD = sqrt(D)
    Q = -B+sD if B<0 else -B-sD
    t0 = Q/2/A
    t1 = 2*C/Q
    e0 = errorformula1(A, B, C, rA, rB, rC)
    e1 = errorformula2(A, B, C, rA, rB, rC)
    if t0<t1:
        t0, t1, e0, e1 = t1, t0, e1, e0
    return t0, t1, e0, e1

# TODO: implement improved solution precision method
#def quadratic_interval(A, B, C):
#     Q = B*B - 4*A*C
#     sQ = interval.imath.sqrt(Q)
#     sol1 = (-B + sQ)/(2*A)
#     sol2 = (-B - sQ)/(2*A)
    # return sol1, sol2
```


```
def print_quadratic_solution_and_float_errors(A, B, C, rA, rB, rC):
    print "A={}, B={}, C={}".format(A, B, C)
    print "========================"
    eps = float(np.finfo(np.float16).eps)
    #print "Float16 Eps: ", eps
    print "Float64 solution:"
    A64 = np.float64(A)
    B64 = np.float64(B)
    C64 = np.float64(C)
    A16, B16, C16 = np.float16(A),np.float16(B),np.float16(C)
    t0_64,t1_64 = quadratic(A64,B64,C64)
    print t0_64, t1_64
    print "Residual:"
    print (A*t0_64*t0_64 + B*t0_64 + C), (A*t1_64*t1_64 + B*t1_64 + C)
    print "---------------------"
    print "Float16 solution:"
    t0_16, t1_16 = quadratic(A16, B16, C16)
    print t0_16, t1_16
    print "Residual:"
    print (A*t0_16*t0_16 + B*t0_16 + C), (A*t1_16*t1_16 + B*t1_16 + C)
    print "---------------------"
    print "Improved Float16 Solution"
    t0_better_16, t1_better_16, e0_better_16, e1_better_16 = better_quadratic(A16, B16, C16, np.float16(eps*rA), np.float16(eps*rB), np.float16(eps*rC))
    print t0_better_16, t1_better_16
    print "Residual:"
    print (A*t0_better_16*t0_better_16 + B*t0_better_16 + C), (A*t1_better_16*t1_better_16 + B*t1_better_16 + C)
    print "---------------------"
    print "Error Estimates"
    rA = rA*eps
    rB = rB*eps
    rC = rC*eps
#     isol1, isol2 = quadratic_interval(eps_half_interval_around(A, rA), eps_half_interval_around(B, rB), eps_half_interval_around(C, rC))
#     def fmt(i, x):
#         return "[{:0.5f}, {:0.5f}]".format(i[0].inf-x, i[0].sup-x)
#     print "Solution 1 bounds around true solution:", fmt(isol1, t0_64)
#     print "Solution 2 bounds around true solution:", fmt(isol2, t1_64)
    assert type(e0_better_16) == np.float16
    assert type(e1_better_16) == np.float16
    print "Analytic bound of solution 1: +/-", e0_better_16
    print "Analytic bound of solution 2: +/-", e1_better_16
    print "Actual errors of improved solution:", (t0_better_16-t0_64), (t1_better_16-t1_64)
    print ""
```


```
print_quadratic_solution_and_float_errors(1, 100, 1, 1, 1, 1)
#print_quadratic_solution_and_float_errors(1, -100, 1, 1, 1, 1)
#print_quadratic_solution_and_float_errors(3, 210, 3002 , 1, 1, 1)
#print_quadratic_solution_and_float_errors(3, -210, 3002 , 1, 1, 1)
#print_quadratic_solution_and_float_errors(1, -100, 0.001 , 1, 1, 1)
#print_quadratic_solution_and_float_errors(3, 100, 3, 100, 0, 0)
#print_quadratic_solution_and_float_errors(3, 100, 3, 0, 100, 0)
#print_quadratic_solution_and_float_errors(3, 100, 3, 0, 0, 100)
#print_quadratic_solution_and_float_errors(3, 100, 3*(1+np.finfo(np.float16).eps*100), 0, 0, 100)
```

    A=1, B=100, C=1
    ========================
    Float64 solution:
    -0.010001000200048793 -99.98999899979995
    Residual:
    1.2212453270876722e-13 0.0
    ---------------------
    Float16 solution:
    0.0 -100.0
    Residual:
    1.0 1.0
    ---------------------
    Improved Float16 Solution
    -0.01 -100.0
    Residual:
    9.999999999998899e-05 1.0
    ---------------------
    Error Estimates
    Analytic bound of solution 1: +/- 4.89e-05
    Analytic bound of solution 2: +/- 0.4883
    Actual errors of improved solution: 1.0002000487927681e-06 -0.010001000200048793
    



```
A, B, C, rA, rB, rC, Q_, e1, e2, e3, e4, e5, e6, S_, em, xi = sympy.symbols('A B C rA rB rC Q e1 e2 e3 e4 e5 e6 S e_m xi')
```


```
roundoff_rel_errs = [e1, e2, e3, e4, e5, e6]
```


```
Q = B*B - 4*A*C
# There are two types of solutions. xi switches the sign, as in the improved solution formula.
S1_exact = (-B + xi*sympy.sqrt(Q))/(2*A)
S2_exact = 2*C/(-B + xi*sympy.sqrt(Q))
prettylatex(S1_exact)
prettylatex(S2_exact)
```


$$\frac{1}{2 A} \left(- B + \xi \sqrt{- 4 A C + B^{2}}\right)$$



$$\frac{2 C}{- B + \xi \sqrt{- 4 A C + B^{2}}}$$



```
# Represent solutions with floating point errors.
# e1  ... e6 are relative roundoff errors after arithmetic operations. 
# Their values depend on the next closest floating point value, and is generally unkown.
Q = B*B - 4*A*C
EQ = (B*B*(1+e3) - 4*A*C*(1+e2))*(1+e1)
S1 = (-B + xi*sympy.sqrt(EQ)*(1+e4))/(2*A)*(1+e5)*(1+e6)
S2 = 2*C/(-B + xi*sympy.sqrt(EQ)*(1+e4))*(1+e5)*(1+e6)
prettylatex(S1)
prettylatex(S2)
```


$$\frac{1}{2 A} \left(- B + \xi \sqrt{\left(e_{1} + 1\right) \left(- 4 A C \left(e_{2} + 1\right) + B^{2} \left(e_{3} + 1\right)\right)} \left(e_{4} + 1\right)\right) \left(e_{5} + 1\right) \left(e_{6} + 1\right)$$



$$\frac{2 C \left(e_{5} + 1\right) \left(e_{6} + 1\right)}{- B + \xi \sqrt{\left(e_{1} + 1\right) \left(- 4 A C \left(e_{2} + 1\right) + B^{2} \left(e_{3} + 1\right)\right)} \left(e_{4} + 1\right)}$$



```
# Do simple error propagation by linearization:
# Taylor series expansion to first order around true coefficients and e1..e6=0.
# That is  F(A + A*rA, ..., e1, ....) \approx F(A, ..., 0, ...) + |dF/dA|*A*rA + .... + |dF/de1|*e1 + ....
derivs = {}
for s in [A, B, C] + roundoff_rel_errs:
    d = sympy.diff(S1,s).subs({e1:0, e2:0, e3:0, e4:0, e5:0, e6:0}).subs(Q , Q_)
    derivs[s] = d
    prettylatex(sympy.Eq(sympy.Derivative(S_, s) ,d))
```


$$\frac{d}{d A} S = - \frac{C \xi}{A \sqrt{Q}} - \frac{1}{2 A^{2}} \left(- B + \sqrt{Q} \xi\right)$$



$$\frac{d}{d B} S = \frac{\frac{B \xi}{\sqrt{Q}} - 1}{2 A}$$



$$\frac{d}{d C} S = - \frac{\xi}{\sqrt{Q}}$$



$$\frac{d}{d e_{1}} S = \frac{\xi}{2 A \sqrt{Q}} \left(- 2 A C + \frac{B^{2}}{2}\right)$$



$$\frac{d}{d e_{2}} S = - \frac{C \xi}{\sqrt{Q}}$$



$$\frac{d}{d e_{3}} S = \frac{B^{2} \xi}{4 A \sqrt{Q}}$$



$$\frac{d}{d e_{4}} S = \frac{\sqrt{Q} \xi}{2 A}$$



$$\frac{d}{d e_{5}} S = \frac{1}{2 A} \left(- B + \sqrt{Q} \xi\right)$$



$$\frac{d}{d e_{6}} S = \frac{1}{2 A} \left(- B + \sqrt{Q} \xi\right)$$



```
# Then estimate each error term by derivative times uncertainty magnitude. 
# EStimate e1...e6 as machine epsilon.
roundoff_total_err = 0
for s in roundoff_rel_errs:
    term = derivs[s]
    roundoff_total_err += sympy.Abs(sympy.ratsimp(term)).subs(Q , Q_)
roundoff_total_err = (sympy.simplify(roundoff_total_err))*em
prettylatex(roundoff_total_err)
```


$$e_{m} \left(\left|{\frac{1}{A} \left(B - \sqrt{Q} \xi\right)}\right| + \frac{\left|{\frac{\sqrt{Q} \xi}{A}}\right|}{2} + \left|{\frac{C \xi}{\sqrt{Q}}}\right| + \frac{1}{4} \left|{\frac{B^{2} \xi}{A \sqrt{Q}}}\right| + \frac{1}{4} \left|{\frac{\xi}{A \sqrt{Q}} \left(4 A C - B^{2}\right)}\right|\right)$$



```
# Take supplied errors for A,B, and C.
propagated_error = 0
for s, rs in [(A, rA), (B, rB), (C, rC)]:
    term = derivs[s]
    err = sympy.Abs((term))*rs
    propagated_error += err
prettylatex(propagated_error)
```


$$rA \left|{- \frac{C \xi}{A \sqrt{Q}} - \frac{1}{2 A^{2}} \left(- B + \sqrt{Q} \xi\right)}\right| + \frac{rB}{2} \left|{\frac{1}{A} \left(\frac{B \xi}{\sqrt{Q}} - 1\right)}\right| + rC \left|{\frac{\xi}{\sqrt{Q}}}\right|$$



```
# In total ...
total_error1 = roundoff_total_err + propagated_error
prettylatex(sympy.Eq(sympy.symbols('ERR1'), total_error1))
```


$$ERR_{1} = e_{m} \left(\left|{\frac{1}{A} \left(B - \sqrt{Q} \xi\right)}\right| + \frac{\left|{\frac{\sqrt{Q} \xi}{A}}\right|}{2} + \left|{\frac{C \xi}{\sqrt{Q}}}\right| + \frac{1}{4} \left|{\frac{B^{2} \xi}{A \sqrt{Q}}}\right| + \frac{1}{4} \left|{\frac{\xi}{A \sqrt{Q}} \left(4 A C - B^{2}\right)}\right|\right) + rA \left|{- \frac{C \xi}{A \sqrt{Q}} - \frac{1}{2 A^{2}} \left(- B + \sqrt{Q} \xi\right)}\right| + \frac{rB}{2} \left|{\frac{1}{A} \left(\frac{B \xi}{\sqrt{Q}} - 1\right)}\right| + rC \left|{\frac{\xi}{\sqrt{Q}}}\right|$$



```
prettylatex(total_error1.subs(B-xi*sympy.sqrt(Q_), sympy.symbols('G')).subs(Q, Q_))
```


$$e_{m} \left(\left|{\frac{G}{A}}\right| + \frac{3}{4} \left|{\frac{\sqrt{Q} \xi}{A}}\right| + \left|{\frac{C \xi}{\sqrt{Q}}}\right| + \frac{1}{4} \left|{\frac{B^{2} \xi}{A \sqrt{Q}}}\right|\right) + rA \left|{\frac{C \xi}{A \sqrt{Q}} - \frac{G}{2 A^{2}}}\right| + \frac{rB}{2} \left|{\frac{1}{A} \left(\frac{B \xi}{\sqrt{Q}} - 1\right)}\right| + rC \left|{\frac{\xi}{\sqrt{Q}}}\right|$$



```
# Now the same for solution type 2
derivs = {}
for s in [A, B, C] + roundoff_rel_errs:
    d = sympy.diff(S2,s).subs({e1:0, e2:0, e3:0, e4:0, e5:0, e6:0}).subs(Q , Q_)
    derivs[s] = d
roundoff_total_err = 0
for s in roundoff_rel_errs:
    term = derivs[s]
    roundoff_total_err += sympy.Abs((term)).subs(Q , Q_)
roundoff_total_err = (sympy.simplify(roundoff_total_err))*em
propagated_error = 0
for s, rs in [(A, rA), (B, rB), (C, rC)]:
    term = derivs[s]
    err = sympy.Abs((term))*rs
    propagated_error += err
total_error2 = roundoff_total_err + propagated_error
prettylatex(sympy.Eq(sympy.symbols('ERR2'),total_error2))
```


$$ERR_{2} = e_{m} \left(4 \left|{\frac{C}{B - \sqrt{Q} \xi}}\right| + 2 \left|{\frac{C \sqrt{Q} \xi}{\left(B - \sqrt{Q} \xi\right)^{2}}}\right| + 4 \left|{\frac{A C^{2} \xi}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + \left|{\frac{B^{2} C \xi}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + \left|{\frac{C \xi \left(4 A C - B^{2}\right)}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right|\right) + 4 rA \left|{\frac{C^{2} \xi}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + 2 rB \left|{\frac{C \left(\frac{B \xi}{\sqrt{Q}} - 1\right)}{\left(B - \sqrt{Q} \xi\right)^{2}}}\right| + rC \left|{\frac{4 A C \xi}{\sqrt{Q} \left(- B + \sqrt{Q} \xi\right)^{2}} + \frac{2}{- B + \sqrt{Q} \xi}}\right|$$



```
prettylatex(total_error2.subs(Q, Q_).subs(B-xi*sympy.sqrt(Q_), sympy.symbols('G')))
```


$$e_{m} \left(4 \left|{\frac{C}{G}}\right| + 3 \left|{\frac{C \xi}{G^{2}} \sqrt{Q}}\right| + 4 \left|{\frac{A C^{2} \xi}{G^{2} \sqrt{Q}}}\right| + \left|{\frac{B^{2} C \xi}{G^{2} \sqrt{Q}}}\right|\right) + 4 rA \left|{\frac{C^{2} \xi}{G^{2} \sqrt{Q}}}\right| + 2 rB \left|{\frac{C}{G^{2}} \left(\frac{B \xi}{\sqrt{Q}} - 1\right)}\right| + rC \left|{\frac{4 A C \xi}{G^{2} \sqrt{Q}} - \frac{2}{G}}\right|$$



```
# Let's try. However, note that all calculations are done in default floating point format.
# But I take the low precision eps of float16. Normally the error calculations themselves would be subject to roundoff error.
def print_absolute_errors(A_val, B_val, C_val, eA, eB, eC, computation_eps, val_type = None):
    if val_type is not None:
        A_val = val_type(A_val)
        B_val = val_type(B_val)
        C_val = val_type(C_val)
        eA = val_type(eA)
        eB = val_type(eB)
        eC = val_type(eC)
    #eps = np.finfo(np.float16).eps
    d = {A : A_val, B : B_val, C : C_val, em : computation_eps}
    d[xi] = 1 if d[B]<0 else -1
    de = d.copy()
    de.update({ rA: eA, rB : eB, rC: eC})
    S1_val = S1_exact.subs(d).evalf()
    S2_val = S2_exact.subs(d).evalf()
    e1_val = total_error1.subs(Q_ , Q).subs(d).evalf()
    e2_val = total_error2.subs(Q_ , Q).subs(d).evalf()
    #if S1_val < S2_val:
    #    S1_val, S2_val, e1_val, e2_val = S2_val, S1_val, e2_val, e1_val
    print "Sol = ", S1_val, S2_val
    print "Err1 = ", e1_val, "(", total_error1.subs(Q_ , Q).subs(de).evalf() ,")" 
    print "Err2 = ", e2_val, "(", total_error2.subs(Q_ , Q).subs(de).evalf() ,")" 

print_absolute_errors(1, 100, 1, 0, 0, 0.1, np.finfo(np.float32).eps, val_type = np.float32)
```

    Sol =  -99.9899902343750 -0.0100010018795729
    Err1 =  100.0*rA + 1.0001*rB + 0.010002*rC + 3.57604003511369e-5 ( 0.00103596039116383 )
    Err2 =  1.0004e-6*rA + 0.00010003*rB + 0.010002*rC + 3.57675622453257e-9 ( 0.00100020365789533 )



```
# Compute the relative error.
rel_err1 = (total_error1 / sympy.Abs(S1_exact.subs(Q,Q_))).subs(Q, Q_)
prettylatex(rel_err1)
```


$$\frac{2}{\left|{\frac{1}{A} \left(B - \sqrt{Q} \xi\right)}\right|} \left(e_{m} \left(\left|{\frac{1}{A} \left(B - \sqrt{Q} \xi\right)}\right| + \frac{3}{4} \left|{\frac{\sqrt{Q} \xi}{A}}\right| + \left|{\frac{C \xi}{\sqrt{Q}}}\right| + \frac{1}{4} \left|{\frac{B^{2} \xi}{A \sqrt{Q}}}\right|\right) + rA \left|{- \frac{C \xi}{A \sqrt{Q}} - \frac{1}{2 A^{2}} \left(- B + \sqrt{Q} \xi\right)}\right| + \frac{rB}{2} \left|{\frac{1}{A} \left(\frac{B \xi}{\sqrt{Q}} - 1\right)}\right| + rC \left|{\frac{\xi}{\sqrt{Q}}}\right|\right)$$



```
# Relative error for type 2 solution.
rel_err2 = (total_error2 / sympy.Abs(S2_exact.subs(Q,Q_)))
rel_err2 = sympy.simplify(rel_err2)
# Sympy appears to have problems canceling absolute terms :-(
prettylatex(rel_err2)
```


$$\frac{1}{\left|{\frac{C}{B - \sqrt{Q} \xi}}\right|} \left(\frac{e_{m}}{2} \left(4 \left|{\frac{C}{B - \sqrt{Q} \xi}}\right| + 2 \left|{\frac{C \sqrt{Q} \xi}{\left(B - \sqrt{Q} \xi\right)^{2}}}\right| + 4 \left|{\frac{A C^{2} \xi}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + \left|{\frac{B^{2} C \xi}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + \left|{\frac{C \xi \left(4 A C - B^{2}\right)}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right|\right) + 2 rA \left|{\frac{C^{2} \xi}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + rB \left|{\frac{C \left(B \xi - \sqrt{Q}\right)}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right| + rC \left|{\frac{2 A C \xi - \sqrt{Q} \left(B - \sqrt{Q} \xi\right)}{\sqrt{Q} \left(B - \sqrt{Q} \xi\right)^{2}}}\right|\right)$$



```
# Study relative errors
eps64 = np.finfo(np.float64).eps
eps32 = np.finfo(np.float32).eps
def compute_error_propagation_coefficients(A_val, B_val, C_val, computation_eps = eps64):
    """ Compute linear factors in the propagation of errors when computation is done using 'computation_eps'. Relative input errors rA,rB,rC are 
        measured in machine epsilons (float32). The printed formula approximates the relative error w.r.t the true solution in units of
        epsilon of float32. """
    dABC = {A : A_val, B : B_val, C : C_val, em : computation_eps}
    dABC[xi] = 1 if d[B]<0 else -1
    print "Solutions: ", S1_exact.subs(Q_, Q).subs(dABC).evalf(), S2_exact.subs(Q_, Q).subs(dABC).evalf()
    rel_error1_in_eps32 = (rel_err1/eps32).subs(Q_ , Q).subs(dABC).subs({rA: A*rA*eps32, rB : B*rB*eps32, rC : C*rC*eps32}).evalf()
    rel_error2_in_eps32 = (rel_err2/eps32).subs(Q_ , Q).subs(dABC).subs({rA: A*rA*eps32, rB : B*rB*eps32, rC : C*rC*eps32}).evalf()
    print "Solution Type 1: {},\nSolution Type 2: {}".format(rel_error1_in_eps32, rel_error2_in_eps32)
    
# Worst catastrophic cancelation in the BB - 4AC term!
compute_error_propagation_coefficients(1, 2.*(1+eps32),  1)
```

    Solutions:  -1.00048840047384 -0.999511837944738
    Solution Type 1: 1024.50009155273*A*rA + 1023.99996948242*B*rB + 1023.50009155273*C*rC + 3.81656206993623e-6,
    Solution Type 2: 1023.50009155273*A*rA + 1023.99996948242*B*rB + 1024.50009155273*C*rC + 3.81656206993623e-6



```
# Same here. Float32 calculations.
compute_error_propagation_coefficients(1, 2.*(1+eps32),  1, computation_eps = eps32)
```

    Solutions:  -1.00048840047384 -0.999511837944738
    Solution Type 1: 1024.50009155273*A*rA + 1023.99996948242*B*rB + 1023.50009155273*C*rC + 2049.00115919127,
    Solution Type 2: 1023.50009155273*A*rA + 1023.99996948242*B*rB + 1024.50009155273*C*rC + 2049.00115919127



```
# Benign cases
compute_error_propagation_coefficients(1, 3,  1, computation_eps = eps64)
compute_error_propagation_coefficients(1, 3, 1, computation_eps = eps32)
```

    Solutions:  -2.61803398874989 -0.381966011250105
    Solution Type 1: 1.17082043397669*A*rA + 0.447213605681647*B*rB + 0.170820383068248*C*rC + 5.95253473659895e-9,
    Solution Type 2: 0.170820383068248*A*rA + 0.44721356884395*B*rB + 1.1708203234636*C*rC + 5.95253473659895e-9
    Solutions:  -2.61803398874989 -0.381966011250105
    Solution Type 1: 1.17082043397669*A*rA + 0.447213605681647*B*rB + 0.170820383068248*C*rC + 3.19574283420307,
    Solution Type 2: 0.170820383068248*A*rA + 0.44721356884395*B*rB + 1.1708203234636*C*rC + 3.19574283420307



```
# Catastrophic cancelation in -B + sqrt(....) prevented by improved solution formula.
compute_error_propagation_coefficients(1, 100, 1, computation_eps = eps64)
compute_error_propagation_coefficients(1, 100, 1, computation_eps = eps32)
```

    Solutions:  -99.9899989998000 -0.0100010002000500
    Solution Type 1: 1.00010005767486*A*rA + 0.0100020003020874*B*rB + 0.000100030004041248*C*rC + 5.58812178674155e-9,
    Solution Type 2: 0.000100030004041248*A*rA + 0.0100020000040344*B*rB + 1.0000999703994*C*rC + 5.58812178674155e-9
    Solutions:  -99.9899989998000 -0.0100010002000500
    Solution Type 1: 1.00010005767486*A*rA + 0.0100020003020874*B*rB + 0.000100030004041248*C*rC + 3.00009995911915,
    Solution Type 2: 0.000100030004041248*A*rA + 0.0100020000040344*B*rB + 1.0000999703994*C*rC + 3.00010009534472



```

```

fe65fb8523be42bdcd418cb27c010d67df57f01b
ipynb
Jupyter Notebook
misc/FloatingPointAdventures.ipynb
DaWelter/NaiveTrace
a904785a0e13c394b2c221bc918cddb41bc8b175
FSFAP
2018-04-25T08:14:14.000Z
2022-01-29T06:19:16.000Z
misc/FloatingPointAdventures.ipynb
DaWelter/NaiveTrace
a904785a0e13c394b2c221bc918cddb41bc8b175
FSFAP
misc/FloatingPointAdventures.ipynb
DaWelter/NaiveTrace
a904785a0e13c394b2c221bc918cddb41bc8b175
FSFAP
Qwen/Qwen-72B
1. YES
2. YES