Quadratic roots, narrow range

Average Error: 28.0 → 0.3

Time: 17.2s

Precision: binary64

Cost: 13632

\[\left(\left(1.0536712127723509 \cdot 10^{-8} < a \land a < 94906265.62425156\right) \land \left(1.0536712127723509 \cdot 10^{-8} < b \land b < 94906265.62425156\right)\right) \land \left(1.0536712127723509 \cdot 10^{-8} < c \land c < 94906265.62425156\right)\]

Math

\[\frac{\left(-b\right) + \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}}{2 \cdot a} \]

↓

\[\frac{c}{\frac{b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}}{-2}} \]

FPCore

(FPCore (a b c)
 :precision binary64
 (/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a)))

↓

(FPCore (a b c)
 :precision binary64
 (/ c (/ (+ b (sqrt (fma c (* a -4.0) (* b b)))) -2.0)))

C

double code(double a, double b, double c) {
	return (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a);
}

↓

double code(double a, double b, double c) {
	return c / ((b + sqrt(fma(c, (a * -4.0), (b * b)))) / -2.0);
}

Julia

function code(a, b, c)
	return Float64(Float64(Float64(-b) + sqrt(Float64(Float64(b * b) - Float64(Float64(4.0 * a) * c)))) / Float64(2.0 * a))
end

↓

function code(a, b, c)
	return Float64(c / Float64(Float64(b + sqrt(fma(c, Float64(a * -4.0), Float64(b * b)))) / -2.0))
end

Wolfram

code[a_, b_, c_] := N[(N[((-b) + N[Sqrt[N[(N[(b * b), $MachinePrecision] - N[(N[(4.0 * a), $MachinePrecision] * c), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / N[(2.0 * a), $MachinePrecision]), $MachinePrecision]

↓

code[a_, b_, c_] := N[(c / N[(N[(b + N[Sqrt[N[(c * N[(a * -4.0), $MachinePrecision] + N[(b * b), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] / -2.0), $MachinePrecision]), $MachinePrecision]

Derivation

1. Initial program 28.0

   \[\frac{\left(-b\right) + \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}}{2 \cdot a} \]

2. Simplified 28.0

   \[\leadsto \color{blue}{\frac{\left(-b\right) + \sqrt{b \cdot b - \left(4 \cdot a\right) \cdot c}}{a \cdot 2}} \]
   Proof
   (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 4 a) c)))) (*.f64 a 2)): 0 points increase in error, 0 points decrease in error
   (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 4 a) c)))) (Rewrite<= *-commutative_binary64 (*.f64 2 a))): 2 points increase in error, 0 points decrease in error

3. Applied egg-rr 27.3

   \[\leadsto \frac{\color{blue}{\frac{\frac{b \cdot b - \mathsf{fma}\left(b, b, c \cdot \left(a \cdot -4\right)\right)}{\sqrt{b + \sqrt{\mathsf{fma}\left(b, b, c \cdot \left(a \cdot -4\right)\right)}}}}{-\sqrt{b + \sqrt{\mathsf{fma}\left(b, b, c \cdot \left(a \cdot -4\right)\right)}}}}}{a \cdot 2} \]

4. Simplified 27.1

   \[\leadsto \frac{\color{blue}{\frac{b \cdot b - \mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}{-\left(b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}\right)}}}{a \cdot 2} \]
   Proof
   (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 c (*.f64 a -4) (*.f64 b b))) (neg.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))) (*.f64 a 2)): 0 points increase in error, 0 points decrease in error
   (/.f64 (/.f64 (-.f64 (*.f64 b b) (Rewrite<= fma-def_binary64 (+.f64 (*.f64 c (*.f64 a -4)) (*.f64 b b)))) (neg.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))) (*.f64 a 2)): 12 points increase in error, 0 points decrease in error
   (/.f64 (/.f64 (-.f64 (*.f64 b b) (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 b b) (*.f64 c (*.f64 a -4))))) (neg.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))) (*.f64 a 2)): 0 points increase in error, 12 points decrease in error
   (/.f64 (/.f64 (-.f64 (*.f64 b b) (Rewrite=> fma-def_binary64 (fma.f64 b b (*.f64 c (*.f64 a -4))))) (neg.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))) (*.f64 a 2)): 0 points increase in error, 12 points decrease in error
   (/.f64 (/.f64 (Rewrite<= /-rgt-identity_binary64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) 1)) (neg.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))) (*.f64 a 2)): 12 points increase in error, 0 points decrease in error
   (/.f64 (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) 1) (neg.f64 (+.f64 b (sqrt.f64 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 c (*.f64 a -4)) (*.f64 b b))))))) (*.f64 a 2)): 12 points increase in error, 0 points decrease in error
   (/.f64 (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) 1) (neg.f64 (+.f64 b (sqrt.f64 (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 b b) (*.f64 c (*.f64 a -4)))))))) (*.f64 a 2)): 12 points increase in error, 0 points decrease in error
   (/.f64 (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) 1) (neg.f64 (+.f64 b (sqrt.f64 (Rewrite=> fma-def_binary64 (fma.f64 b b (*.f64 c (*.f64 a -4)))))))) (*.f64 a 2)): 0 points increase in error, 12 points decrease in error
   (/.f64 (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) 1) (neg.f64 (Rewrite<= rem-square-sqrt_binary64 (*.f64 (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4)))))) (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4)))))))))) (*.f64 a 2)): 12 points increase in error, 0 points decrease in error
   (/.f64 (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) 1) (Rewrite=> distribute-lft-neg-in_binary64 (*.f64 (neg.f64 (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4))))))) (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4))))))))) (*.f64 a 2)): 0 points increase in error, 12 points decrease in error
   (/.f64 (/.f64 (Rewrite=> /-rgt-identity_binary64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4))))) (*.f64 (neg.f64 (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4))))))) (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4)))))))) (*.f64 a 2)): 12 points increase in error, 0 points decrease in error
   (/.f64 (Rewrite<= associate-/l/_binary64 (/.f64 (/.f64 (-.f64 (*.f64 b b) (fma.f64 b b (*.f64 c (*.f64 a -4)))) (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4))))))) (neg.f64 (sqrt.f64 (+.f64 b (sqrt.f64 (fma.f64 b b (*.f64 c (*.f64 a -4))))))))) (*.f64 a 2)): 0 points increase in error, 12 points decrease in error

5. Taylor expanded in b around 0 0.4

   \[\leadsto \frac{\frac{\color{blue}{4 \cdot \left(c \cdot a\right)}}{-\left(b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}\right)}}{a \cdot 2} \]

6. Applied egg-rr 0.5

   \[\leadsto \color{blue}{\frac{4 \cdot \frac{0.5}{a}}{b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}} \cdot \left(-c \cdot a\right)} \]

7. Simplified 0.5

   \[\leadsto \color{blue}{\left(c \cdot \left(-a\right)\right) \cdot \frac{\frac{2}{a}}{b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}}} \]
   Proof
   (*.f64 (*.f64 c (neg.f64 a)) (/.f64 (/.f64 2 a) (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))): 0 points increase in error, 0 points decrease in error
   (*.f64 (Rewrite<= distribute-rgt-neg-in_binary64 (neg.f64 (*.f64 c a))) (/.f64 (/.f64 2 a) (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))): 0 points increase in error, 0 points decrease in error
   (*.f64 (neg.f64 (*.f64 c a)) (/.f64 (/.f64 (Rewrite<= metadata-eval (*.f64 4 1/2)) a) (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))): 5 points increase in error, 0 points decrease in error
   (*.f64 (neg.f64 (*.f64 c a)) (/.f64 (Rewrite<= associate-*r/_binary64 (*.f64 4 (/.f64 1/2 a))) (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))))): 0 points increase in error, 0 points decrease in error
   (Rewrite<= *-commutative_binary64 (*.f64 (/.f64 (*.f64 4 (/.f64 1/2 a)) (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b))))) (neg.f64 (*.f64 c a)))): 0 points increase in error, 0 points decrease in error

8. Applied egg-rr 0.5

   \[\leadsto \color{blue}{\frac{-1}{\frac{\left(b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}\right) \cdot \left(a \cdot 0.5\right)}{a \cdot c}}} \]

9. Simplified 0.3

   \[\leadsto \color{blue}{\frac{c}{\frac{b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}}{-2}}} \]
   Proof
   (/.f64 c (/.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) -2)): 0 points increase in error, 0 points decrease in error
   (/.f64 c (/.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (Rewrite<= metadata-eval (*.f64 -1 2)))): 1 points increase in error, 0 points decrease in error
   (/.f64 c (/.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 -1 (Rewrite<= metadata-eval (/.f64 1 1/2))))): 0 points increase in error, 1 points decrease in error
   (/.f64 c (/.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 -1 (/.f64 (Rewrite<= *-inverses_binary64 (/.f64 a a)) 1/2)))): 10 points increase in error, 0 points decrease in error
   (/.f64 c (/.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 -1 (Rewrite<= associate-/r*_binary64 (/.f64 a (*.f64 a 1/2)))))): 0 points increase in error, 4 points decrease in error
   (/.f64 c (Rewrite<= associate-/l/_binary64 (/.f64 (/.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (/.f64 a (*.f64 a 1/2))) -1))): 4 points increase in error, 0 points decrease in error
   (/.f64 c (/.f64 (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 a 1/2)) a)) -1)): 1 points increase in error, 9 points decrease in error
   (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 c -1) (/.f64 (*.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 a 1/2)) a))): 9 points increase in error, 0 points decrease in error
   (/.f64 (Rewrite=> *-commutative_binary64 (*.f64 -1 c)) (/.f64 (*.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 a 1/2)) a)): 0 points increase in error, 6 points decrease in error
   (Rewrite=> associate-/l*_binary64 (/.f64 -1 (/.f64 (/.f64 (*.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 a 1/2)) a) c))): 0 points increase in error, 5 points decrease in error
   (/.f64 -1 (Rewrite<= associate-/r*_binary64 (/.f64 (*.f64 (+.f64 b (sqrt.f64 (fma.f64 c (*.f64 a -4) (*.f64 b b)))) (*.f64 a 1/2)) (*.f64 a c)))): 1 points increase in error, 0 points decrease in error

10. Final simplification 0.3

    \[\leadsto \frac{c}{\frac{b + \sqrt{\mathsf{fma}\left(c, a \cdot -4, b \cdot b\right)}}{-2}} \]

Alternatives

Alternative 1
Error	9.4
Cost	14788

\[\begin{array}{l} t_0 := \frac{\sqrt{b \cdot b + c \cdot \left(a \cdot -4\right)} - b}{a \cdot 2}\\ \mathbf{if}\;t_0 \leq -1:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;\frac{\left(c \cdot 4\right) \cdot \frac{a}{-2 \cdot \left(a \cdot \frac{c}{b} - b\right)}}{a \cdot -2}\\ \end{array} \]

Alternative 2
Error	0.4
Cost	7808

\[\frac{\frac{4 \cdot \left(c \cdot a\right)}{\left(-b\right) - \sqrt{b \cdot b + c \cdot \left(a \cdot -4\right)}}}{a \cdot 2} \]

Alternative 3
Error	0.5
Cost	7616

\[\left(c \cdot a\right) \cdot \frac{\frac{-2}{a}}{b + \sqrt{b \cdot b + c \cdot \left(a \cdot -4\right)}} \]

Alternative 4
Error	9.6
Cost	7492

\[\begin{array}{l} \mathbf{if}\;b \leq 1.44:\\ \;\;\;\;\left(\sqrt{b \cdot b + a \cdot \left(c \cdot -4\right)} - b\right) \cdot \frac{0.5}{a}\\ \mathbf{else}:\\ \;\;\;\;\frac{\left(c \cdot 4\right) \cdot \frac{a}{-2 \cdot \left(a \cdot \frac{c}{b} - b\right)}}{a \cdot -2}\\ \end{array} \]

Alternative 5
Error	11.7
Cost	1216

\[\frac{\left(c \cdot 4\right) \cdot \frac{a}{-2 \cdot \left(a \cdot \frac{c}{b} - b\right)}}{a \cdot -2} \]

Alternative 6
Error	11.7
Cost	1088

\[\frac{\left(c \cdot a\right) \cdot 2}{b - a \cdot \frac{c}{b}} \cdot \frac{-0.5}{a} \]

Alternative 7
Error	11.7
Cost	1088

\[\frac{-2 \cdot \frac{c \cdot a}{b - a \cdot \frac{c}{b}}}{a \cdot 2} \]

Alternative 8
Error	12.0
Cost	1024

\[\frac{-c}{b} - \frac{c}{b} \cdot \frac{c \cdot a}{b \cdot b} \]

Alternative 9
Error	23.2
Cost	256

\[\frac{-c}{b} \]

Alternative 10
Error	63.0
Cost	192

\[\frac{b}{a} \]

Reproduce

herbie shell --seed 2022341 
(FPCore (a b c)
  :name "Quadratic roots, narrow range"
  :precision binary64
  :pre (and (and (and (< 1.0536712127723509e-8 a) (< a 94906265.62425156)) (and (< 1.0536712127723509e-8 b) (< b 94906265.62425156))) (and (< 1.0536712127723509e-8 c) (< c 94906265.62425156)))
  (/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a)))