Henrywood and Agarwal, Equation (13)

Average Error: 59.4 → 30.3

Time: 15.1s

Precision: binary64

Math

\[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

↓

\[\begin{array}{l} \mathbf{if}\;d \leq -2.1 \cdot 10^{+71}:\\ \;\;\;\;0\\ \mathbf{elif}\;d \leq -1.55 \cdot 10^{-162}:\\ \;\;\;\;0.25 \cdot \left(h \cdot \left(\left(M \cdot M\right) \cdot \frac{D \cdot D}{d \cdot d}\right)\right)\\ \mathbf{elif}\;d \leq 1.75 \cdot 10^{-182}:\\ \;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(\left(\frac{d}{D} \cdot \frac{d}{D}\right) \cdot \left(2 \cdot \frac{c0}{h \cdot w}\right)\right)\\ \mathbf{elif}\;d \leq 4.2 \cdot 10^{-86}:\\ \;\;\;\;\mathsf{fma}\left(0.25, \left(\frac{D}{d} \cdot \frac{D}{d}\right) \cdot \left(M \cdot \left(h \cdot M\right)\right), \frac{0}{w} \cdot \left(c0 \cdot c0\right)\right)\\ \mathbf{elif}\;d \leq 2.3 \cdot 10^{+152}:\\ \;\;\;\;0.25 \cdot \frac{{D}^{2} \cdot \left(h \cdot {M}^{2}\right)}{{d}^{2}}\\ \mathbf{else}:\\ \;\;\;\;0\\ \end{array} \]

FPCore

(FPCore (c0 w h D d M)
 :precision binary64
 (*
  (/ c0 (* 2.0 w))
  (+
   (/ (* c0 (* d d)) (* (* w h) (* D D)))
   (sqrt
    (-
     (*
      (/ (* c0 (* d d)) (* (* w h) (* D D)))
      (/ (* c0 (* d d)) (* (* w h) (* D D))))
     (* M M))))))

↓

(FPCore (c0 w h D d M)
 :precision binary64
 (if (<= d -2.1e+71)
   0.0
   (if (<= d -1.55e-162)
     (* 0.25 (* h (* (* M M) (/ (* D D) (* d d)))))
     (if (<= d 1.75e-182)
       (* (/ c0 (* 2.0 w)) (* (* (/ d D) (/ d D)) (* 2.0 (/ c0 (* h w)))))
       (if (<= d 4.2e-86)
         (fma
          0.25
          (* (* (/ D d) (/ D d)) (* M (* h M)))
          (* (/ 0.0 w) (* c0 c0)))
         (if (<= d 2.3e+152)
           (* 0.25 (/ (* (pow D 2.0) (* h (pow M 2.0))) (pow d 2.0)))
           0.0))))))

C

double code(double c0, double w, double h, double D, double d, double M) {
	return (c0 / (2.0 * w)) * (((c0 * (d * d)) / ((w * h) * (D * D))) + sqrt(((((c0 * (d * d)) / ((w * h) * (D * D))) * ((c0 * (d * d)) / ((w * h) * (D * D)))) - (M * M))));
}

↓

double code(double c0, double w, double h, double D, double d, double M) {
	double tmp;
	if (d <= -2.1e+71) {
		tmp = 0.0;
	} else if (d <= -1.55e-162) {
		tmp = 0.25 * (h * ((M * M) * ((D * D) / (d * d))));
	} else if (d <= 1.75e-182) {
		tmp = (c0 / (2.0 * w)) * (((d / D) * (d / D)) * (2.0 * (c0 / (h * w))));
	} else if (d <= 4.2e-86) {
		tmp = fma(0.25, (((D / d) * (D / d)) * (M * (h * M))), ((0.0 / w) * (c0 * c0)));
	} else if (d <= 2.3e+152) {
		tmp = 0.25 * ((pow(D, 2.0) * (h * pow(M, 2.0))) / pow(d, 2.0));
	} else {
		tmp = 0.0;
	}
	return tmp;
}

Julia

function code(c0, w, h, D, d, M)
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D))) + sqrt(Float64(Float64(Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D))) * Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D)))) - Float64(M * M)))))
end

↓

function code(c0, w, h, D, d, M)
	tmp = 0.0
	if (d <= -2.1e+71)
		tmp = 0.0;
	elseif (d <= -1.55e-162)
		tmp = Float64(0.25 * Float64(h * Float64(Float64(M * M) * Float64(Float64(D * D) / Float64(d * d)))));
	elseif (d <= 1.75e-182)
		tmp = Float64(Float64(c0 / Float64(2.0 * w)) * Float64(Float64(Float64(d / D) * Float64(d / D)) * Float64(2.0 * Float64(c0 / Float64(h * w)))));
	elseif (d <= 4.2e-86)
		tmp = fma(0.25, Float64(Float64(Float64(D / d) * Float64(D / d)) * Float64(M * Float64(h * M))), Float64(Float64(0.0 / w) * Float64(c0 * c0)));
	elseif (d <= 2.3e+152)
		tmp = Float64(0.25 * Float64(Float64((D ^ 2.0) * Float64(h * (M ^ 2.0))) / (d ^ 2.0)));
	else
		tmp = 0.0;
	end
	return tmp
end

Wolfram

code[c0_, w_, h_, D_, d_, M_] := N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[Sqrt[N[(N[(N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] - N[(M * M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

↓

code[c0_, w_, h_, D_, d_, M_] := If[LessEqual[d, -2.1e+71], 0.0, If[LessEqual[d, -1.55e-162], N[(0.25 * N[(h * N[(N[(M * M), $MachinePrecision] * N[(N[(D * D), $MachinePrecision] / N[(d * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 1.75e-182], N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(d / D), $MachinePrecision] * N[(d / D), $MachinePrecision]), $MachinePrecision] * N[(2.0 * N[(c0 / N[(h * w), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 4.2e-86], N[(0.25 * N[(N[(N[(D / d), $MachinePrecision] * N[(D / d), $MachinePrecision]), $MachinePrecision] * N[(M * N[(h * M), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(N[(0.0 / w), $MachinePrecision] * N[(c0 * c0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 2.3e+152], N[(0.25 * N[(N[(N[Power[D, 2.0], $MachinePrecision] * N[(h * N[Power[M, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[Power[d, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 0.0]]]]]

Derivation

1. Split input into 5 regimes

2. if d < -2.09999999999999989e71 or 2.29999999999999985e152 < d

   1. Initial program 62.7

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

   2. Taylor expanded in c0 around -inf 62.6

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(-1 \cdot \left(\left(\frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)} + -1 \cdot \frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)}\right) \cdot c0\right)\right)} \]

   3. Simplified 30.6

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(c0 \cdot 0\right)} \]

   4. Taylor expanded in c0 around 0 25.8

      \[\leadsto \color{blue}{0} \]

   if -2.09999999999999989e71 < d < -1.5499999999999999e-162

   1. Initial program 54.3

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

   2. Taylor expanded in c0 around -inf 55.5

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(0.5 \cdot \frac{{D}^{2} \cdot \left(w \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot c0} + -1 \cdot \left(\left(\frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)} + -1 \cdot \frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)}\right) \cdot c0\right)\right)} \]

   3. Simplified 37.1

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\mathsf{fma}\left(0.5, \left(\frac{D}{d} \cdot \frac{D}{d}\right) \cdot \frac{w}{\frac{c0}{M \cdot \left(M \cdot h\right)}}, c0 \cdot 0\right)} \]

   4. Taylor expanded in c0 around 0 32.3

      \[\leadsto \color{blue}{0.25 \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2}}} \]

   5. Simplified 32.4

      \[\leadsto \color{blue}{0.25 \cdot \left(h \cdot \left(\left(M \cdot M\right) \cdot \frac{D \cdot D}{d \cdot d}\right)\right)} \]

   if -1.5499999999999999e-162 < d < 1.74999999999999992e-182

   1. Initial program 63.1

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

   2. Taylor expanded in c0 around inf 62.9

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(2 \cdot \frac{{d}^{2} \cdot c0}{{D}^{2} \cdot \left(w \cdot h\right)}\right)} \]

   3. Simplified 47.1

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(\left(\frac{d}{D} \cdot \frac{d}{D}\right) \cdot \left(\frac{c0}{w \cdot h} \cdot 2\right)\right)} \]

   if 1.74999999999999992e-182 < d < 4.2e-86

   1. Initial program 57.6

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

   2. Taylor expanded in c0 around -inf 61.1

      \[\leadsto \color{blue}{-0.5 \cdot \frac{\left(\frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)} + -1 \cdot \frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)}\right) \cdot {c0}^{2}}{w} + 0.25 \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2}}} \]

   3. Simplified 37.3

      \[\leadsto \color{blue}{\mathsf{fma}\left(0.25, \left(\frac{D}{d} \cdot \frac{D}{d}\right) \cdot \left(M \cdot \left(M \cdot h\right)\right), \frac{0}{w} \cdot \left(c0 \cdot c0\right)\right)} \]

   if 4.2e-86 < d < 2.29999999999999985e152

   1. Initial program 54.9

      \[\frac{c0}{2 \cdot w} \cdot \left(\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} + \sqrt{\frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} \cdot \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)} - M \cdot M}\right) \]

   2. Taylor expanded in c0 around -inf 54.7

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\left(0.5 \cdot \frac{{D}^{2} \cdot \left(w \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot c0} + -1 \cdot \left(\left(\frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)} + -1 \cdot \frac{{d}^{2}}{{D}^{2} \cdot \left(w \cdot h\right)}\right) \cdot c0\right)\right)} \]

   3. Simplified 35.6

      \[\leadsto \frac{c0}{2 \cdot w} \cdot \color{blue}{\mathsf{fma}\left(0.5, \left(\frac{D}{d} \cdot \frac{D}{d}\right) \cdot \frac{w}{\frac{c0}{M \cdot \left(M \cdot h\right)}}, c0 \cdot 0\right)} \]

   4. Taylor expanded in w around inf 30.7

      \[\leadsto \color{blue}{0.25 \cdot \frac{{D}^{2} \cdot \left(h \cdot {M}^{2}\right)}{{d}^{2}}} \]
3. Recombined 5 regimes into one program.

4. Final simplification 30.3

   \[\leadsto \begin{array}{l} \mathbf{if}\;d \leq -2.1 \cdot 10^{+71}:\\ \;\;\;\;0\\ \mathbf{elif}\;d \leq -1.55 \cdot 10^{-162}:\\ \;\;\;\;0.25 \cdot \left(h \cdot \left(\left(M \cdot M\right) \cdot \frac{D \cdot D}{d \cdot d}\right)\right)\\ \mathbf{elif}\;d \leq 1.75 \cdot 10^{-182}:\\ \;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(\left(\frac{d}{D} \cdot \frac{d}{D}\right) \cdot \left(2 \cdot \frac{c0}{h \cdot w}\right)\right)\\ \mathbf{elif}\;d \leq 4.2 \cdot 10^{-86}:\\ \;\;\;\;\mathsf{fma}\left(0.25, \left(\frac{D}{d} \cdot \frac{D}{d}\right) \cdot \left(M \cdot \left(h \cdot M\right)\right), \frac{0}{w} \cdot \left(c0 \cdot c0\right)\right)\\ \mathbf{elif}\;d \leq 2.3 \cdot 10^{+152}:\\ \;\;\;\;0.25 \cdot \frac{{D}^{2} \cdot \left(h \cdot {M}^{2}\right)}{{d}^{2}}\\ \mathbf{else}:\\ \;\;\;\;0\\ \end{array} \]

Reproduce

herbie shell --seed 2022211 
(FPCore (c0 w h D d M)
  :name "Henrywood and Agarwal, Equation (13)"
  :precision binary64
  (* (/ c0 (* 2.0 w)) (+ (/ (* c0 (* d d)) (* (* w h) (* D D))) (sqrt (- (* (/ (* c0 (* d d)) (* (* w h) (* D D))) (/ (* c0 (* d d)) (* (* w h) (* D D)))) (* M M))))))