Cubic critical, medium range

Average Error: 43.5 → 0.6

Time: 20.1s

Precision: binary64

Cost: 14144

\[\left(\left(1.1102230246251565 \cdot 10^{-16} < a \land a < 9007199254740992\right) \land \left(1.1102230246251565 \cdot 10^{-16} < b \land b < 9007199254740992\right)\right) \land \left(1.1102230246251565 \cdot 10^{-16} < c \land c < 9007199254740992\right)\]

Math

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

↓

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

FPCore

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

↓

(FPCore (a b c)
 :precision binary64
 (*
  (/ (* 3.0 (/ c a)) (+ (/ b a) (/ (sqrt (fma a (* c -3.0) (* b b))) a)))
  -0.3333333333333333))

C

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

↓

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

Julia

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

↓

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

Wolfram

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

↓

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

Derivation

1. Initial program 43.5

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

2. Simplified 43.5

   \[\leadsto \color{blue}{\frac{b - \sqrt{\mathsf{fma}\left(a, c \cdot -3, b \cdot b\right)}}{a} \cdot -0.3333333333333333} \]
   Proof
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b)))) a) -1/3): 0 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (fma.f64 a (*.f64 c (Rewrite<= metadata-eval (neg.f64 3))) (*.f64 b b)))) a) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (Rewrite<= fma-def_binary64 (+.f64 (*.f64 a (*.f64 c (neg.f64 3))) (*.f64 b b))))) a) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 (*.f64 c (neg.f64 3)) a)) (*.f64 b b)))) a) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (+.f64 (Rewrite<= associate-*r*_binary64 (*.f64 c (*.f64 (neg.f64 3) a))) (*.f64 b b)))) a) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (+.f64 (*.f64 c (Rewrite<= distribute-lft-neg-in_binary64 (neg.f64 (*.f64 3 a)))) (*.f64 b b)))) a) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (+.f64 (Rewrite<= *-commutative_binary64 (*.f64 (neg.f64 (*.f64 3 a)) c)) (*.f64 b b)))) a) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (Rewrite<= +-commutative_binary64 (+.f64 (*.f64 b b) (*.f64 (neg.f64 (*.f64 3 a)) c))))) a) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (-.f64 b (sqrt.f64 (Rewrite<= cancel-sign-sub-inv_binary64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) a) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (Rewrite<= *-lft-identity_binary64 (*.f64 1 (/.f64 (-.f64 b (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) a))) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (*.f64 (Rewrite<= metadata-eval (/.f64 -1 -1)) (/.f64 (-.f64 b (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) a)) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (Rewrite<= times-frac_binary64 (/.f64 (*.f64 -1 (-.f64 b (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) (*.f64 -1 a))) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (/.f64 (*.f64 -1 (-.f64 b (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) (Rewrite<= neg-mul-1_binary64 (neg.f64 a))) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (Rewrite<= neg-mul-1_binary64 (neg.f64 (-.f64 b (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))))) (neg.f64 a)) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (Rewrite<= sub0-neg_binary64 (-.f64 0 (-.f64 b (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))))) (neg.f64 a)) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (/.f64 (Rewrite<= associate-+l-_binary64 (+.f64 (-.f64 0 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) (neg.f64 a)) -1/3): 0 points increase in error, 25 points decrease in error
   (*.f64 (/.f64 (+.f64 (Rewrite<= neg-sub0_binary64 (neg.f64 b)) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) (neg.f64 a)) -1/3): 25 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) (neg.f64 a)) (Rewrite<= metadata-eval (/.f64 -1 3))): 25 points increase in error, 0 points decrease in error
   (Rewrite<= *-commutative_binary64 (*.f64 (/.f64 -1 3) (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) (neg.f64 a)))): 0 points increase in error, 25 points decrease in error
   (Rewrite<= times-frac_binary64 (/.f64 (*.f64 -1 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) (*.f64 3 (neg.f64 a)))): 25 points increase in error, 0 points decrease in error
   (/.f64 (*.f64 -1 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) (Rewrite<= distribute-rgt-neg-in_binary64 (neg.f64 (*.f64 3 a)))): 0 points increase in error, 25 points decrease in error
   (/.f64 (*.f64 -1 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c))))) (Rewrite=> neg-mul-1_binary64 (*.f64 -1 (*.f64 3 a)))): 0 points increase in error, 25 points decrease in error
   (Rewrite=> times-frac_binary64 (*.f64 (/.f64 -1 -1) (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) (*.f64 3 a)))): 25 points increase in error, 0 points decrease in error
   (*.f64 (Rewrite=> metadata-eval 1) (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) (*.f64 3 a))): 25 points increase in error, 0 points decrease in error
   (Rewrite=> *-lft-identity_binary64 (/.f64 (+.f64 (neg.f64 b) (sqrt.f64 (-.f64 (*.f64 b b) (*.f64 (*.f64 3 a) c)))) (*.f64 3 a))): 0 points increase in error, 25 points decrease in error

3. Applied egg-rr 43.8

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

4. Simplified 43.1

   \[\leadsto \color{blue}{\frac{\frac{b}{a} \cdot \frac{b}{a} - \frac{\frac{\mathsf{fma}\left(a, c \cdot -3, b \cdot b\right)}{a}}{a}}{\frac{b}{a} + \frac{\sqrt{\mathsf{fma}\left(a, c \cdot -3, b \cdot b\right)}}{a}} \cdot -0.3333333333333333} \]
   Proof
   (*.f64 (/.f64 (-.f64 (*.f64 (/.f64 b a) (/.f64 b a)) (/.f64 (/.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b)) a) a)) (+.f64 (/.f64 b a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) -1/3): 0 points increase in error, 0 points decrease in error
   (*.f64 (/.f64 (-.f64 (*.f64 (/.f64 b a) (/.f64 b a)) (/.f64 (/.f64 (Rewrite<= rem-square-sqrt_binary64 (*.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))))) a) a)) (+.f64 (/.f64 b a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) -1/3): 4 points increase in error, 2 points decrease in error
   (*.f64 (/.f64 (-.f64 (*.f64 (/.f64 b a) (/.f64 b a)) (/.f64 (Rewrite<= associate-*r/_binary64 (*.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) a)) (+.f64 (/.f64 b a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) -1/3): 0 points increase in error, 6 points decrease in error
   (*.f64 (/.f64 (-.f64 (*.f64 (/.f64 b a) (/.f64 b a)) (Rewrite<= associate-*l/_binary64 (*.f64 (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a)))) (+.f64 (/.f64 b a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) -1/3): 6 points increase in error, 0 points decrease in error
   (Rewrite<= associate-/r/_binary64 (/.f64 (-.f64 (*.f64 (/.f64 b a) (/.f64 b a)) (*.f64 (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) (/.f64 (+.f64 (/.f64 b a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a)) -1/3))): 2 points increase in error, 4 points decrease in error
   (Rewrite<= associate-/l*_binary64 (/.f64 (*.f64 (-.f64 (*.f64 (/.f64 b a) (/.f64 b a)) (*.f64 (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a))) -1/3) (+.f64 (/.f64 b a) (/.f64 (sqrt.f64 (fma.f64 a (*.f64 c -3) (*.f64 b b))) a)))): 0 points increase in error, 6 points decrease in error

5. Taylor expanded in b around 0 0.6

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

6. Final simplification 0.6

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

Alternatives

Alternative 1
Error	6.2
Cost	13764

\[\begin{array}{l} \mathbf{if}\;b \leq 0.00025:\\ \;\;\;\;\left(b - \sqrt{\mathsf{fma}\left(b, b, c \cdot \left(a \cdot -3\right)\right)}\right) \cdot \frac{-0.3333333333333333}{a}\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(-0.375, \frac{c \cdot c}{b \cdot \frac{b \cdot b}{a}}, -0.5 \cdot \frac{c}{b}\right)\\ \end{array} \]

Alternative 2
Error	6.2
Cost	7620

\[\begin{array}{l} \mathbf{if}\;b \leq 0.00025:\\ \;\;\;\;-0.3333333333333333 \cdot \frac{b - \sqrt{b \cdot b + a \cdot \left(c \cdot -3\right)}}{a}\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(-0.375, \frac{c \cdot c}{b \cdot \frac{b \cdot b}{a}}, -0.5 \cdot \frac{c}{b}\right)\\ \end{array} \]

Alternative 3
Error	6.3
Cost	7492

\[\begin{array}{l} \mathbf{if}\;b \leq 0.00025:\\ \;\;\;\;\frac{-0.3333333333333333}{a} \cdot \left(b - \sqrt{b \cdot b + a \cdot \left(c \cdot -3\right)}\right)\\ \mathbf{else}:\\ \;\;\;\;c \cdot \left(\frac{-0.5}{b} + a \cdot \left(-0.375 \cdot \frac{c}{{b}^{3}}\right)\right)\\ \end{array} \]

Alternative 4
Error	6.3
Cost	7492

\[\begin{array}{l} \mathbf{if}\;b \leq 0.00025:\\ \;\;\;\;-0.3333333333333333 \cdot \frac{b - \sqrt{b \cdot b + a \cdot \left(c \cdot -3\right)}}{a}\\ \mathbf{else}:\\ \;\;\;\;c \cdot \left(\frac{-0.5}{b} + a \cdot \left(-0.375 \cdot \frac{c}{{b}^{3}}\right)\right)\\ \end{array} \]

Alternative 5
Error	6.4
Cost	7296

\[c \cdot \left(\frac{-0.5}{b} + a \cdot \left(-0.375 \cdot \frac{c}{{b}^{3}}\right)\right) \]

Alternative 6
Error	12.3
Cost	320

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

Reproduce

herbie shell --seed 2022340 
(FPCore (a b c)
  :name "Cubic critical, medium range"
  :precision binary64
  :pre (and (and (and (< 1.1102230246251565e-16 a) (< a 9007199254740992.0)) (and (< 1.1102230246251565e-16 b) (< b 9007199254740992.0))) (and (< 1.1102230246251565e-16 c) (< c 9007199254740992.0)))
  (/ (+ (- b) (sqrt (- (* b b) (* (* 3.0 a) c)))) (* 3.0 a)))