Quadratic roots, wide range

Percentage Accurate: 17.8% → 97.5%

Time: 13.1s

Alternatives: 7

Speedup: 3.6×

Specification

\[\left(\left(4.930380657631324 \cdot 10^{-32} < a \land a < 2.028240960365167 \cdot 10^{+31}\right) \land \left(4.930380657631324 \cdot 10^{-32} < b \land b < 2.028240960365167 \cdot 10^{+31}\right)\right) \land \left(4.930380657631324 \cdot 10^{-32} < c \land c < 2.028240960365167 \cdot 10^{+31}\right)\]

Math

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

FPCore

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

C

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

Fortran

real(8) function code(a, b, c)
    real(8), intent (in) :: a
    real(8), intent (in) :: b
    real(8), intent (in) :: c
    code = (-b + sqrt(((b * b) - ((4.0d0 * a) * c)))) / (2.0d0 * a)
end function

Java

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

Python

def code(a, b, c):
	return (-b + math.sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a)

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

MATLAB

function tmp = code(a, b, c)
	tmp = (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a);
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]

Initial Program: 17.8% accurate, 1.0× speedup

Math

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

FPCore

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

C

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

Fortran

real(8) function code(a, b, c)
    real(8), intent (in) :: a
    real(8), intent (in) :: b
    real(8), intent (in) :: c
    code = (-b + sqrt(((b * b) - ((4.0d0 * a) * c)))) / (2.0d0 * a)
end function

Java

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

Python

def code(a, b, c):
	return (-b + math.sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a)

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

MATLAB

function tmp = code(a, b, c)
	tmp = (-b + sqrt(((b * b) - ((4.0 * a) * c)))) / (2.0 * a);
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]

Alternative 1: 97.5% accurate, 0.2× speedup

Math

\[\begin{array}{l} \\ \begin{array}{l} t_0 := b \cdot \left(b \cdot b\right)\\ \mathsf{fma}\left(\mathsf{fma}\left(\frac{\left(\left(c \cdot c\right) \cdot \left(c \cdot c\right)\right) \cdot \left(a \cdot 20\right)}{b \cdot \left(t\_0 \cdot t\_0\right)}, -0.25, \left(c \cdot \left(c \cdot c\right)\right) \cdot \left(-2 \cdot {b}^{-5}\right)\right), a \cdot a, \frac{\mathsf{fma}\left(c, c \cdot \frac{a}{b \cdot b}, c\right)}{-b}\right) \end{array} \end{array} \]

FPCore

(FPCore (a b c)
 :precision binary64
 (let* ((t_0 (* b (* b b))))
   (fma
    (fma
     (/ (* (* (* c c) (* c c)) (* a 20.0)) (* b (* t_0 t_0)))
     -0.25
     (* (* c (* c c)) (* -2.0 (pow b -5.0))))
    (* a a)
    (/ (fma c (* c (/ a (* b b))) c) (- b)))))

C

double code(double a, double b, double c) {
	double t_0 = b * (b * b);
	return fma(fma(((((c * c) * (c * c)) * (a * 20.0)) / (b * (t_0 * t_0))), -0.25, ((c * (c * c)) * (-2.0 * pow(b, -5.0)))), (a * a), (fma(c, (c * (a / (b * b))), c) / -b));
}

Julia

function code(a, b, c)
	t_0 = Float64(b * Float64(b * b))
	return fma(fma(Float64(Float64(Float64(Float64(c * c) * Float64(c * c)) * Float64(a * 20.0)) / Float64(b * Float64(t_0 * t_0))), -0.25, Float64(Float64(c * Float64(c * c)) * Float64(-2.0 * (b ^ -5.0)))), Float64(a * a), Float64(fma(c, Float64(c * Float64(a / Float64(b * b))), c) / Float64(-b)))
end

Wolfram

code[a_, b_, c_] := Block[{t$95$0 = N[(b * N[(b * b), $MachinePrecision]), $MachinePrecision]}, N[(N[(N[(N[(N[(N[(c * c), $MachinePrecision] * N[(c * c), $MachinePrecision]), $MachinePrecision] * N[(a * 20.0), $MachinePrecision]), $MachinePrecision] / N[(b * N[(t$95$0 * t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * -0.25 + N[(N[(c * N[(c * c), $MachinePrecision]), $MachinePrecision] * N[(-2.0 * N[Power[b, -5.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(a * a), $MachinePrecision] + N[(N[(c * N[(c * N[(a / N[(b * b), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + c), $MachinePrecision] / (-b)), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 15.7%

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

3. Taylor expanded in a around 0

   \[\leadsto \color{blue}{-1 \cdot \frac{c}{b} + a \cdot \left(-1 \cdot \frac{{c}^{2}}{{b}^{3}} + a \cdot \left(-2 \cdot \frac{{c}^{3}}{{b}^{5}} + \frac{-1}{4} \cdot \frac{a \cdot \left(4 \cdot \frac{{c}^{4}}{{b}^{6}} + 16 \cdot \frac{{c}^{4}}{{b}^{6}}\right)}{b}\right)\right)} \]

5. Applied rewrites 98.4%

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

7. Applied rewrites 98.4%

   \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\frac{\left(\left(c \cdot c\right) \cdot \left(c \cdot c\right)\right) \cdot \left(20 \cdot a\right)}{b \cdot \left(\left(b \cdot \left(b \cdot b\right)\right) \cdot \left(b \cdot \left(b \cdot b\right)\right)\right)}, -0.25, \left(c \cdot \left(c \cdot c\right)\right) \cdot \left(-2 \cdot {b}^{-5}\right)\right), \color{blue}{a \cdot a}, \frac{\mathsf{fma}\left(c, c \cdot \frac{a}{b \cdot b}, c\right)}{-b}\right) \]

8. Final simplification 98.4%

   \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\frac{\left(\left(c \cdot c\right) \cdot \left(c \cdot c\right)\right) \cdot \left(a \cdot 20\right)}{b \cdot \left(\left(b \cdot \left(b \cdot b\right)\right) \cdot \left(b \cdot \left(b \cdot b\right)\right)\right)}, -0.25, \left(c \cdot \left(c \cdot c\right)\right) \cdot \left(-2 \cdot {b}^{-5}\right)\right), a \cdot a, \frac{\mathsf{fma}\left(c, c \cdot \frac{a}{b \cdot b}, c\right)}{-b}\right) \]
9. Add Preprocessing

Alternative 2: 97.2% accurate, 0.2× speedup

Math

\[\begin{array}{l} \\ \begin{array}{l} t_0 := b \cdot \left(b \cdot b\right)\\ \mathsf{fma}\left(\mathsf{fma}\left(c \cdot c, \frac{a}{b \cdot b}, c\right), \frac{-1}{b}, \left(a \cdot a\right) \cdot \mathsf{fma}\left(\frac{-5 \cdot \left(a \cdot \left(c \cdot c\right)\right)}{b \cdot \left(t\_0 \cdot t\_0\right)}, c \cdot c, \left(\left(c \cdot c\right) \cdot -2\right) \cdot \left(c \cdot {b}^{-5}\right)\right)\right) \end{array} \end{array} \]

FPCore

(FPCore (a b c)
 :precision binary64
 (let* ((t_0 (* b (* b b))))
   (fma
    (fma (* c c) (/ a (* b b)) c)
    (/ -1.0 b)
    (*
     (* a a)
     (fma
      (/ (* -5.0 (* a (* c c))) (* b (* t_0 t_0)))
      (* c c)
      (* (* (* c c) -2.0) (* c (pow b -5.0))))))))

C

double code(double a, double b, double c) {
	double t_0 = b * (b * b);
	return fma(fma((c * c), (a / (b * b)), c), (-1.0 / b), ((a * a) * fma(((-5.0 * (a * (c * c))) / (b * (t_0 * t_0))), (c * c), (((c * c) * -2.0) * (c * pow(b, -5.0))))));
}

Julia

function code(a, b, c)
	t_0 = Float64(b * Float64(b * b))
	return fma(fma(Float64(c * c), Float64(a / Float64(b * b)), c), Float64(-1.0 / b), Float64(Float64(a * a) * fma(Float64(Float64(-5.0 * Float64(a * Float64(c * c))) / Float64(b * Float64(t_0 * t_0))), Float64(c * c), Float64(Float64(Float64(c * c) * -2.0) * Float64(c * (b ^ -5.0))))))
end

Wolfram

code[a_, b_, c_] := Block[{t$95$0 = N[(b * N[(b * b), $MachinePrecision]), $MachinePrecision]}, N[(N[(N[(c * c), $MachinePrecision] * N[(a / N[(b * b), $MachinePrecision]), $MachinePrecision] + c), $MachinePrecision] * N[(-1.0 / b), $MachinePrecision] + N[(N[(a * a), $MachinePrecision] * N[(N[(N[(-5.0 * N[(a * N[(c * c), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(b * N[(t$95$0 * t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(c * c), $MachinePrecision] + N[(N[(N[(c * c), $MachinePrecision] * -2.0), $MachinePrecision] * N[(c * N[Power[b, -5.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 17.8%

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

3. Taylor expanded in a around 0

   \[\leadsto \color{blue}{-1 \cdot \frac{c}{b} + a \cdot \left(-1 \cdot \frac{{c}^{2}}{{b}^{3}} + a \cdot \left(-2 \cdot \frac{{c}^{3}}{{b}^{5}} + \frac{-1}{4} \cdot \frac{a \cdot \left(4 \cdot \frac{{c}^{4}}{{b}^{6}} + 16 \cdot \frac{{c}^{4}}{{b}^{6}}\right)}{b}\right)\right)} \]

5. Applied rewrites 97.5%

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

7. Applied rewrites 97.5%

   \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(\frac{\left(\left(c \cdot c\right) \cdot \left(c \cdot c\right)\right) \cdot \left(20 \cdot a\right)}{b \cdot \left(\left(b \cdot \left(b \cdot b\right)\right) \cdot \left(b \cdot \left(b \cdot b\right)\right)\right)}, -0.25, \left(c \cdot \left(c \cdot c\right)\right) \cdot \left(-2 \cdot {b}^{-5}\right)\right), \color{blue}{a \cdot a}, \frac{\mathsf{fma}\left(c, c \cdot \frac{a}{b \cdot b}, c\right)}{-b}\right) \]

9. Applied rewrites 97.2%

   \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(c \cdot c, \frac{a}{b \cdot b}, c\right), \color{blue}{\frac{-1}{b}}, \mathsf{fma}\left(c, \left(c \cdot c\right) \cdot \left(-2 \cdot {b}^{-5}\right), \frac{\left(c \cdot c\right) \cdot \left(\left(\left(c \cdot \left(c \cdot a\right)\right) \cdot 20\right) \cdot -0.25\right)}{b \cdot \left(\left(b \cdot b\right) \cdot \left(\left(b \cdot b\right) \cdot \left(b \cdot b\right)\right)\right)}\right) \cdot \left(a \cdot a\right)\right) \]
10. Step-by-step derivation

11. Applied rewrites 97.2%

    \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(c \cdot c, \frac{a}{b \cdot b}, c\right), \frac{-1}{b}, \mathsf{fma}\left(\frac{-5 \cdot \left(\left(c \cdot c\right) \cdot a\right)}{b \cdot \left(\left(b \cdot \left(b \cdot b\right)\right) \cdot \left(b \cdot \left(b \cdot b\right)\right)\right)}, c \cdot c, \left(\left(c \cdot c\right) \cdot -2\right) \cdot \left({b}^{-5} \cdot c\right)\right) \cdot \left(a \cdot a\right)\right) \]

12. Final simplification 97.2%

    \[\leadsto \mathsf{fma}\left(\mathsf{fma}\left(c \cdot c, \frac{a}{b \cdot b}, c\right), \frac{-1}{b}, \left(a \cdot a\right) \cdot \mathsf{fma}\left(\frac{-5 \cdot \left(a \cdot \left(c \cdot c\right)\right)}{b \cdot \left(\left(b \cdot \left(b \cdot b\right)\right) \cdot \left(b \cdot \left(b \cdot b\right)\right)\right)}, c \cdot c, \left(\left(c \cdot c\right) \cdot -2\right) \cdot \left(c \cdot {b}^{-5}\right)\right)\right) \]
13. Add Preprocessing

Reproduce

herbie shell --seed 2024222 
(FPCore (a b c)
  :name "Quadratic roots, wide range"
  :precision binary64
  :pre (and (and (and (< 4.930380657631324e-32 a) (< a 2.028240960365167e+31)) (and (< 4.930380657631324e-32 b) (< b 2.028240960365167e+31))) (and (< 4.930380657631324e-32 c) (< c 2.028240960365167e+31)))
  (/ (+ (- b) (sqrt (- (* b b) (* (* 4.0 a) c)))) (* 2.0 a)))