Henrywood and Agarwal, Equation (13)

Percentage Accurate: 24.7% → 42.8%
Time: 31.3s
Alternatives: 8
Speedup: 16.8×

Specification

?
\[\begin{array}{l} \\ \begin{array}{l} t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\ \frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right) \end{array} \end{array} \]
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (/ (* c0 (* d d)) (* (* w h) (* D D)))))
   (* (/ c0 (* 2.0 w)) (+ t_0 (sqrt (- (* t_0 t_0) (* M M)))))))
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
}
real(8) function code(c0, w, h, d, d_1, m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d
    real(8), intent (in) :: d_1
    real(8), intent (in) :: m
    real(8) :: t_0
    t_0 = (c0 * (d_1 * d_1)) / ((w * h) * (d * d))
    code = (c0 / (2.0d0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (m * m))))
end function
public static double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + Math.sqrt(((t_0 * t_0) - (M * M))));
}
def code(c0, w, h, D, d, M):
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D))
	return (c0 / (2.0 * w)) * (t_0 + math.sqrt(((t_0 * t_0) - (M * M))))
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D)))
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_0 + sqrt(Float64(Float64(t_0 * t_0) - Float64(M * M)))))
end
function tmp = code(c0, w, h, D, d, M)
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	tmp = (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
end
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$0 + N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] - N[(M * M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}

\\
\begin{array}{l}
t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\
\frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right)
\end{array}
\end{array}

Sampling outcomes in binary64 precision:

Local Percentage Accuracy vs ?

The average percentage accuracy by input value. Horizontal axis shows value of an input variable; the variable is choosen in the title. Vertical axis is accuracy; higher is better. Red represent the original program, while blue represents Herbie's suggestion. These can be toggled with buttons below the plot. The line is an average while dots represent individual samples.

Accuracy vs Speed?

Herbie found 8 alternatives:

AlternativeAccuracySpeedup
The accuracy (vertical axis) and speed (horizontal axis) of each alternatives. Up and to the right is better. The red square shows the initial program, and each blue circle shows an alternative.The line shows the best available speed-accuracy tradeoffs.

Initial Program: 24.7% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \begin{array}{l} t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\ \frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right) \end{array} \end{array} \]
(FPCore (c0 w h D d M)
 :precision binary64
 (let* ((t_0 (/ (* c0 (* d d)) (* (* w h) (* D D)))))
   (* (/ c0 (* 2.0 w)) (+ t_0 (sqrt (- (* t_0 t_0) (* M M)))))))
double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
}
real(8) function code(c0, w, h, d, d_1, m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d
    real(8), intent (in) :: d_1
    real(8), intent (in) :: m
    real(8) :: t_0
    t_0 = (c0 * (d_1 * d_1)) / ((w * h) * (d * d))
    code = (c0 / (2.0d0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (m * m))))
end function
public static double code(double c0, double w, double h, double D, double d, double M) {
	double t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	return (c0 / (2.0 * w)) * (t_0 + Math.sqrt(((t_0 * t_0) - (M * M))));
}
def code(c0, w, h, D, d, M):
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D))
	return (c0 / (2.0 * w)) * (t_0 + math.sqrt(((t_0 * t_0) - (M * M))))
function code(c0, w, h, D, d, M)
	t_0 = Float64(Float64(c0 * Float64(d * d)) / Float64(Float64(w * h) * Float64(D * D)))
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_0 + sqrt(Float64(Float64(t_0 * t_0) - Float64(M * M)))))
end
function tmp = code(c0, w, h, D, d, M)
	t_0 = (c0 * (d * d)) / ((w * h) * (D * D));
	tmp = (c0 / (2.0 * w)) * (t_0 + sqrt(((t_0 * t_0) - (M * M))));
end
code[c0_, w_, h_, D_, d_, M_] := Block[{t$95$0 = N[(N[(c0 * N[(d * d), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D * D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$0 + N[Sqrt[N[(N[(t$95$0 * t$95$0), $MachinePrecision] - N[(M * M), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}

\\
\begin{array}{l}
t_0 := \frac{c0 \cdot \left(d \cdot d\right)}{\left(w \cdot h\right) \cdot \left(D \cdot D\right)}\\
\frac{c0}{2 \cdot w} \cdot \left(t_0 + \sqrt{t_0 \cdot t_0 - M \cdot M}\right)
\end{array}
\end{array}

Alternative 1: 42.8% accurate, 0.2× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\\ t_1 := \frac{c0}{2 \cdot w}\\ t_2 := \frac{c0}{w \cdot h}\\ t_3 := \frac{\frac{c0}{w}}{h}\\ t_4 := {\left(\frac{d_m}{D_m}\right)}^{2}\\ t_5 := \frac{t_4}{\frac{h}{\frac{c0}{w}}}\\ \mathbf{if}\;c0 \leq -6.5 \cdot 10^{+153}:\\ \;\;\;\;\frac{\frac{c0}{w}}{2} \cdot \mathsf{fma}\left(t_3, t_0, \sqrt{\mathsf{fma}\left(t_3, t_0, M_m\right) \cdot \left(t_3 \cdot t_0 - M_m\right)}\right)\\ \mathbf{elif}\;c0 \leq -2.4 \cdot 10^{+135}:\\ \;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\ \mathbf{elif}\;c0 \leq -7500:\\ \;\;\;\;t_1 \cdot \mathsf{fma}\left(\sqrt{\mathsf{fma}\left(t_2, t_4, M_m\right)}, \frac{d_m}{D_m} \cdot \sqrt{t_2}, t_2 \cdot t_4\right)\\ \mathbf{elif}\;c0 \leq 45:\\ \;\;\;\;\log \left({\left(\sqrt{e^{\frac{c0}{w}}}\right)}^{\left(2 \cdot \left({d_m}^{2} \cdot \frac{c0}{\left(w \cdot h\right) \cdot {D_m}^{2}}\right)\right)}\right)\\ \mathbf{else}:\\ \;\;\;\;t_1 \cdot \mathsf{fma}\left(\sqrt{\mathsf{fma}\left(t_3, t_4, M_m\right)}, \sqrt{t_5 - M_m}, t_5\right)\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (* (/ d_m D_m) (/ d_m D_m)))
        (t_1 (/ c0 (* 2.0 w)))
        (t_2 (/ c0 (* w h)))
        (t_3 (/ (/ c0 w) h))
        (t_4 (pow (/ d_m D_m) 2.0))
        (t_5 (/ t_4 (/ h (/ c0 w)))))
   (if (<= c0 -6.5e+153)
     (*
      (/ (/ c0 w) 2.0)
      (fma t_3 t_0 (sqrt (* (fma t_3 t_0 M_m) (- (* t_3 t_0) M_m)))))
     (if (<= c0 -2.4e+135)
       (* -0.5 (/ (pow c0 2.0) (/ w 0.0)))
       (if (<= c0 -7500.0)
         (*
          t_1
          (fma
           (sqrt (fma t_2 t_4 M_m))
           (* (/ d_m D_m) (sqrt t_2))
           (* t_2 t_4)))
         (if (<= c0 45.0)
           (log
            (pow
             (sqrt (exp (/ c0 w)))
             (* 2.0 (* (pow d_m 2.0) (/ c0 (* (* w h) (pow D_m 2.0)))))))
           (* t_1 (fma (sqrt (fma t_3 t_4 M_m)) (sqrt (- t_5 M_m)) t_5))))))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = (d_m / D_m) * (d_m / D_m);
	double t_1 = c0 / (2.0 * w);
	double t_2 = c0 / (w * h);
	double t_3 = (c0 / w) / h;
	double t_4 = pow((d_m / D_m), 2.0);
	double t_5 = t_4 / (h / (c0 / w));
	double tmp;
	if (c0 <= -6.5e+153) {
		tmp = ((c0 / w) / 2.0) * fma(t_3, t_0, sqrt((fma(t_3, t_0, M_m) * ((t_3 * t_0) - M_m))));
	} else if (c0 <= -2.4e+135) {
		tmp = -0.5 * (pow(c0, 2.0) / (w / 0.0));
	} else if (c0 <= -7500.0) {
		tmp = t_1 * fma(sqrt(fma(t_2, t_4, M_m)), ((d_m / D_m) * sqrt(t_2)), (t_2 * t_4));
	} else if (c0 <= 45.0) {
		tmp = log(pow(sqrt(exp((c0 / w))), (2.0 * (pow(d_m, 2.0) * (c0 / ((w * h) * pow(D_m, 2.0)))))));
	} else {
		tmp = t_1 * fma(sqrt(fma(t_3, t_4, M_m)), sqrt((t_5 - M_m)), t_5);
	}
	return tmp;
}
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(Float64(d_m / D_m) * Float64(d_m / D_m))
	t_1 = Float64(c0 / Float64(2.0 * w))
	t_2 = Float64(c0 / Float64(w * h))
	t_3 = Float64(Float64(c0 / w) / h)
	t_4 = Float64(d_m / D_m) ^ 2.0
	t_5 = Float64(t_4 / Float64(h / Float64(c0 / w)))
	tmp = 0.0
	if (c0 <= -6.5e+153)
		tmp = Float64(Float64(Float64(c0 / w) / 2.0) * fma(t_3, t_0, sqrt(Float64(fma(t_3, t_0, M_m) * Float64(Float64(t_3 * t_0) - M_m)))));
	elseif (c0 <= -2.4e+135)
		tmp = Float64(-0.5 * Float64((c0 ^ 2.0) / Float64(w / 0.0)));
	elseif (c0 <= -7500.0)
		tmp = Float64(t_1 * fma(sqrt(fma(t_2, t_4, M_m)), Float64(Float64(d_m / D_m) * sqrt(t_2)), Float64(t_2 * t_4)));
	elseif (c0 <= 45.0)
		tmp = log((sqrt(exp(Float64(c0 / w))) ^ Float64(2.0 * Float64((d_m ^ 2.0) * Float64(c0 / Float64(Float64(w * h) * (D_m ^ 2.0)))))));
	else
		tmp = Float64(t_1 * fma(sqrt(fma(t_3, t_4, M_m)), sqrt(Float64(t_5 - M_m)), t_5));
	end
	return tmp
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(N[(d$95$m / D$95$m), $MachinePrecision] * N[(d$95$m / D$95$m), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(c0 / N[(w * h), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$3 = N[(N[(c0 / w), $MachinePrecision] / h), $MachinePrecision]}, Block[{t$95$4 = N[Power[N[(d$95$m / D$95$m), $MachinePrecision], 2.0], $MachinePrecision]}, Block[{t$95$5 = N[(t$95$4 / N[(h / N[(c0 / w), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[c0, -6.5e+153], N[(N[(N[(c0 / w), $MachinePrecision] / 2.0), $MachinePrecision] * N[(t$95$3 * t$95$0 + N[Sqrt[N[(N[(t$95$3 * t$95$0 + M$95$m), $MachinePrecision] * N[(N[(t$95$3 * t$95$0), $MachinePrecision] - M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[c0, -2.4e+135], N[(-0.5 * N[(N[Power[c0, 2.0], $MachinePrecision] / N[(w / 0.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[c0, -7500.0], N[(t$95$1 * N[(N[Sqrt[N[(t$95$2 * t$95$4 + M$95$m), $MachinePrecision]], $MachinePrecision] * N[(N[(d$95$m / D$95$m), $MachinePrecision] * N[Sqrt[t$95$2], $MachinePrecision]), $MachinePrecision] + N[(t$95$2 * t$95$4), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[c0, 45.0], N[Log[N[Power[N[Sqrt[N[Exp[N[(c0 / w), $MachinePrecision]], $MachinePrecision]], $MachinePrecision], N[(2.0 * N[(N[Power[d$95$m, 2.0], $MachinePrecision] * N[(c0 / N[(N[(w * h), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision], N[(t$95$1 * N[(N[Sqrt[N[(t$95$3 * t$95$4 + M$95$m), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(t$95$5 - M$95$m), $MachinePrecision]], $MachinePrecision] + t$95$5), $MachinePrecision]), $MachinePrecision]]]]]]]]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\\
t_1 := \frac{c0}{2 \cdot w}\\
t_2 := \frac{c0}{w \cdot h}\\
t_3 := \frac{\frac{c0}{w}}{h}\\
t_4 := {\left(\frac{d_m}{D_m}\right)}^{2}\\
t_5 := \frac{t_4}{\frac{h}{\frac{c0}{w}}}\\
\mathbf{if}\;c0 \leq -6.5 \cdot 10^{+153}:\\
\;\;\;\;\frac{\frac{c0}{w}}{2} \cdot \mathsf{fma}\left(t_3, t_0, \sqrt{\mathsf{fma}\left(t_3, t_0, M_m\right) \cdot \left(t_3 \cdot t_0 - M_m\right)}\right)\\

\mathbf{elif}\;c0 \leq -2.4 \cdot 10^{+135}:\\
\;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\

\mathbf{elif}\;c0 \leq -7500:\\
\;\;\;\;t_1 \cdot \mathsf{fma}\left(\sqrt{\mathsf{fma}\left(t_2, t_4, M_m\right)}, \frac{d_m}{D_m} \cdot \sqrt{t_2}, t_2 \cdot t_4\right)\\

\mathbf{elif}\;c0 \leq 45:\\
\;\;\;\;\log \left({\left(\sqrt{e^{\frac{c0}{w}}}\right)}^{\left(2 \cdot \left({d_m}^{2} \cdot \frac{c0}{\left(w \cdot h\right) \cdot {D_m}^{2}}\right)\right)}\right)\\

\mathbf{else}:\\
\;\;\;\;t_1 \cdot \mathsf{fma}\left(\sqrt{\mathsf{fma}\left(t_3, t_4, M_m\right)}, \sqrt{t_5 - M_m}, t_5\right)\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 2: 44.4% accurate, 0.1× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{c0}{2 \cdot w}\\ t_1 := \frac{c0}{w \cdot h}\\ t_2 := {\left(\frac{d_m}{D_m}\right)}^{2}\\ t_3 := \frac{c0 \cdot \left(d_m \cdot d_m\right)}{\left(w \cdot h\right) \cdot \left(D_m \cdot D_m\right)}\\ \mathbf{if}\;t_0 \cdot \left(t_3 + \sqrt{t_3 \cdot t_3 - M_m \cdot M_m}\right) \leq \infty:\\ \;\;\;\;t_0 \cdot \frac{2 \cdot \left(c0 \cdot {d_m}^{2}\right)}{\left(w \cdot h\right) \cdot {D_m}^{2}}\\ \mathbf{else}:\\ \;\;\;\;\log \left({\left(\sqrt{e^{\frac{c0}{w}}}\right)}^{\left(\mathsf{fma}\left(t_1, t_2, \sqrt{\mathsf{fma}\left(t_1, t_2, M_m\right) \cdot \mathsf{fma}\left(\frac{\frac{{d_m}^{2}}{w}}{h \cdot D_m}, \frac{c0}{D_m}, -M_m\right)}\right)\right)}\right)\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (/ c0 (* 2.0 w)))
        (t_1 (/ c0 (* w h)))
        (t_2 (pow (/ d_m D_m) 2.0))
        (t_3 (/ (* c0 (* d_m d_m)) (* (* w h) (* D_m D_m)))))
   (if (<= (* t_0 (+ t_3 (sqrt (- (* t_3 t_3) (* M_m M_m))))) INFINITY)
     (* t_0 (/ (* 2.0 (* c0 (pow d_m 2.0))) (* (* w h) (pow D_m 2.0))))
     (log
      (pow
       (sqrt (exp (/ c0 w)))
       (fma
        t_1
        t_2
        (sqrt
         (*
          (fma t_1 t_2 M_m)
          (fma (/ (/ (pow d_m 2.0) w) (* h D_m)) (/ c0 D_m) (- M_m))))))))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (2.0 * w);
	double t_1 = c0 / (w * h);
	double t_2 = pow((d_m / D_m), 2.0);
	double t_3 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	double tmp;
	if ((t_0 * (t_3 + sqrt(((t_3 * t_3) - (M_m * M_m))))) <= ((double) INFINITY)) {
		tmp = t_0 * ((2.0 * (c0 * pow(d_m, 2.0))) / ((w * h) * pow(D_m, 2.0)));
	} else {
		tmp = log(pow(sqrt(exp((c0 / w))), fma(t_1, t_2, sqrt((fma(t_1, t_2, M_m) * fma(((pow(d_m, 2.0) / w) / (h * D_m)), (c0 / D_m), -M_m))))));
	}
	return tmp;
}
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(c0 / Float64(2.0 * w))
	t_1 = Float64(c0 / Float64(w * h))
	t_2 = Float64(d_m / D_m) ^ 2.0
	t_3 = Float64(Float64(c0 * Float64(d_m * d_m)) / Float64(Float64(w * h) * Float64(D_m * D_m)))
	tmp = 0.0
	if (Float64(t_0 * Float64(t_3 + sqrt(Float64(Float64(t_3 * t_3) - Float64(M_m * M_m))))) <= Inf)
		tmp = Float64(t_0 * Float64(Float64(2.0 * Float64(c0 * (d_m ^ 2.0))) / Float64(Float64(w * h) * (D_m ^ 2.0))));
	else
		tmp = log((sqrt(exp(Float64(c0 / w))) ^ fma(t_1, t_2, sqrt(Float64(fma(t_1, t_2, M_m) * fma(Float64(Float64((d_m ^ 2.0) / w) / Float64(h * D_m)), Float64(c0 / D_m), Float64(-M_m)))))));
	end
	return tmp
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(c0 / N[(w * h), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[Power[N[(d$95$m / D$95$m), $MachinePrecision], 2.0], $MachinePrecision]}, Block[{t$95$3 = N[(N[(c0 * N[(d$95$m * d$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D$95$m * D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(t$95$0 * N[(t$95$3 + N[Sqrt[N[(N[(t$95$3 * t$95$3), $MachinePrecision] - N[(M$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], Infinity], N[(t$95$0 * N[(N[(2.0 * N[(c0 * N[Power[d$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[Log[N[Power[N[Sqrt[N[Exp[N[(c0 / w), $MachinePrecision]], $MachinePrecision]], $MachinePrecision], N[(t$95$1 * t$95$2 + N[Sqrt[N[(N[(t$95$1 * t$95$2 + M$95$m), $MachinePrecision] * N[(N[(N[(N[Power[d$95$m, 2.0], $MachinePrecision] / w), $MachinePrecision] / N[(h * D$95$m), $MachinePrecision]), $MachinePrecision] * N[(c0 / D$95$m), $MachinePrecision] + (-M$95$m)), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision]]]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{c0}{2 \cdot w}\\
t_1 := \frac{c0}{w \cdot h}\\
t_2 := {\left(\frac{d_m}{D_m}\right)}^{2}\\
t_3 := \frac{c0 \cdot \left(d_m \cdot d_m\right)}{\left(w \cdot h\right) \cdot \left(D_m \cdot D_m\right)}\\
\mathbf{if}\;t_0 \cdot \left(t_3 + \sqrt{t_3 \cdot t_3 - M_m \cdot M_m}\right) \leq \infty:\\
\;\;\;\;t_0 \cdot \frac{2 \cdot \left(c0 \cdot {d_m}^{2}\right)}{\left(w \cdot h\right) \cdot {D_m}^{2}}\\

\mathbf{else}:\\
\;\;\;\;\log \left({\left(\sqrt{e^{\frac{c0}{w}}}\right)}^{\left(\mathsf{fma}\left(t_1, t_2, \sqrt{\mathsf{fma}\left(t_1, t_2, M_m\right) \cdot \mathsf{fma}\left(\frac{\frac{{d_m}^{2}}{w}}{h \cdot D_m}, \frac{c0}{D_m}, -M_m\right)}\right)\right)}\right)\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 3: 42.5% accurate, 0.2× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\\ t_1 := \frac{c0}{w \cdot h}\\ t_2 := \frac{\frac{c0}{w}}{h}\\ t_3 := \frac{\frac{c0}{w}}{2} \cdot \mathsf{fma}\left(t_2, t_0, \sqrt{\mathsf{fma}\left(t_2, t_0, M_m\right) \cdot \left(t_2 \cdot t_0 - M_m\right)}\right)\\ t_4 := {\left(\frac{d_m}{D_m}\right)}^{2}\\ \mathbf{if}\;c0 \leq -2.3 \cdot 10^{+151}:\\ \;\;\;\;t_3\\ \mathbf{elif}\;c0 \leq -1.2 \cdot 10^{+135}:\\ \;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\ \mathbf{elif}\;c0 \leq -120000:\\ \;\;\;\;\frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(\sqrt{\mathsf{fma}\left(t_1, t_4, M_m\right)}, \frac{d_m}{D_m} \cdot \sqrt{t_1}, t_1 \cdot t_4\right)\\ \mathbf{elif}\;c0 \leq 8.2 \cdot 10^{-50}:\\ \;\;\;\;\log \left({\left(\sqrt{e^{\frac{c0}{w}}}\right)}^{\left(2 \cdot \left({d_m}^{2} \cdot \frac{c0}{\left(w \cdot h\right) \cdot {D_m}^{2}}\right)\right)}\right)\\ \mathbf{else}:\\ \;\;\;\;t_3\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (* (/ d_m D_m) (/ d_m D_m)))
        (t_1 (/ c0 (* w h)))
        (t_2 (/ (/ c0 w) h))
        (t_3
         (*
          (/ (/ c0 w) 2.0)
          (fma t_2 t_0 (sqrt (* (fma t_2 t_0 M_m) (- (* t_2 t_0) M_m))))))
        (t_4 (pow (/ d_m D_m) 2.0)))
   (if (<= c0 -2.3e+151)
     t_3
     (if (<= c0 -1.2e+135)
       (* -0.5 (/ (pow c0 2.0) (/ w 0.0)))
       (if (<= c0 -120000.0)
         (*
          (/ c0 (* 2.0 w))
          (fma
           (sqrt (fma t_1 t_4 M_m))
           (* (/ d_m D_m) (sqrt t_1))
           (* t_1 t_4)))
         (if (<= c0 8.2e-50)
           (log
            (pow
             (sqrt (exp (/ c0 w)))
             (* 2.0 (* (pow d_m 2.0) (/ c0 (* (* w h) (pow D_m 2.0)))))))
           t_3))))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = (d_m / D_m) * (d_m / D_m);
	double t_1 = c0 / (w * h);
	double t_2 = (c0 / w) / h;
	double t_3 = ((c0 / w) / 2.0) * fma(t_2, t_0, sqrt((fma(t_2, t_0, M_m) * ((t_2 * t_0) - M_m))));
	double t_4 = pow((d_m / D_m), 2.0);
	double tmp;
	if (c0 <= -2.3e+151) {
		tmp = t_3;
	} else if (c0 <= -1.2e+135) {
		tmp = -0.5 * (pow(c0, 2.0) / (w / 0.0));
	} else if (c0 <= -120000.0) {
		tmp = (c0 / (2.0 * w)) * fma(sqrt(fma(t_1, t_4, M_m)), ((d_m / D_m) * sqrt(t_1)), (t_1 * t_4));
	} else if (c0 <= 8.2e-50) {
		tmp = log(pow(sqrt(exp((c0 / w))), (2.0 * (pow(d_m, 2.0) * (c0 / ((w * h) * pow(D_m, 2.0)))))));
	} else {
		tmp = t_3;
	}
	return tmp;
}
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(Float64(d_m / D_m) * Float64(d_m / D_m))
	t_1 = Float64(c0 / Float64(w * h))
	t_2 = Float64(Float64(c0 / w) / h)
	t_3 = Float64(Float64(Float64(c0 / w) / 2.0) * fma(t_2, t_0, sqrt(Float64(fma(t_2, t_0, M_m) * Float64(Float64(t_2 * t_0) - M_m)))))
	t_4 = Float64(d_m / D_m) ^ 2.0
	tmp = 0.0
	if (c0 <= -2.3e+151)
		tmp = t_3;
	elseif (c0 <= -1.2e+135)
		tmp = Float64(-0.5 * Float64((c0 ^ 2.0) / Float64(w / 0.0)));
	elseif (c0 <= -120000.0)
		tmp = Float64(Float64(c0 / Float64(2.0 * w)) * fma(sqrt(fma(t_1, t_4, M_m)), Float64(Float64(d_m / D_m) * sqrt(t_1)), Float64(t_1 * t_4)));
	elseif (c0 <= 8.2e-50)
		tmp = log((sqrt(exp(Float64(c0 / w))) ^ Float64(2.0 * Float64((d_m ^ 2.0) * Float64(c0 / Float64(Float64(w * h) * (D_m ^ 2.0)))))));
	else
		tmp = t_3;
	end
	return tmp
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(N[(d$95$m / D$95$m), $MachinePrecision] * N[(d$95$m / D$95$m), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(c0 / N[(w * h), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(N[(c0 / w), $MachinePrecision] / h), $MachinePrecision]}, Block[{t$95$3 = N[(N[(N[(c0 / w), $MachinePrecision] / 2.0), $MachinePrecision] * N[(t$95$2 * t$95$0 + N[Sqrt[N[(N[(t$95$2 * t$95$0 + M$95$m), $MachinePrecision] * N[(N[(t$95$2 * t$95$0), $MachinePrecision] - M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$4 = N[Power[N[(d$95$m / D$95$m), $MachinePrecision], 2.0], $MachinePrecision]}, If[LessEqual[c0, -2.3e+151], t$95$3, If[LessEqual[c0, -1.2e+135], N[(-0.5 * N[(N[Power[c0, 2.0], $MachinePrecision] / N[(w / 0.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[c0, -120000.0], N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(N[Sqrt[N[(t$95$1 * t$95$4 + M$95$m), $MachinePrecision]], $MachinePrecision] * N[(N[(d$95$m / D$95$m), $MachinePrecision] * N[Sqrt[t$95$1], $MachinePrecision]), $MachinePrecision] + N[(t$95$1 * t$95$4), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[c0, 8.2e-50], N[Log[N[Power[N[Sqrt[N[Exp[N[(c0 / w), $MachinePrecision]], $MachinePrecision]], $MachinePrecision], N[(2.0 * N[(N[Power[d$95$m, 2.0], $MachinePrecision] * N[(c0 / N[(N[(w * h), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]], $MachinePrecision], t$95$3]]]]]]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\\
t_1 := \frac{c0}{w \cdot h}\\
t_2 := \frac{\frac{c0}{w}}{h}\\
t_3 := \frac{\frac{c0}{w}}{2} \cdot \mathsf{fma}\left(t_2, t_0, \sqrt{\mathsf{fma}\left(t_2, t_0, M_m\right) \cdot \left(t_2 \cdot t_0 - M_m\right)}\right)\\
t_4 := {\left(\frac{d_m}{D_m}\right)}^{2}\\
\mathbf{if}\;c0 \leq -2.3 \cdot 10^{+151}:\\
\;\;\;\;t_3\\

\mathbf{elif}\;c0 \leq -1.2 \cdot 10^{+135}:\\
\;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\

\mathbf{elif}\;c0 \leq -120000:\\
\;\;\;\;\frac{c0}{2 \cdot w} \cdot \mathsf{fma}\left(\sqrt{\mathsf{fma}\left(t_1, t_4, M_m\right)}, \frac{d_m}{D_m} \cdot \sqrt{t_1}, t_1 \cdot t_4\right)\\

\mathbf{elif}\;c0 \leq 8.2 \cdot 10^{-50}:\\
\;\;\;\;\log \left({\left(\sqrt{e^{\frac{c0}{w}}}\right)}^{\left(2 \cdot \left({d_m}^{2} \cdot \frac{c0}{\left(w \cdot h\right) \cdot {D_m}^{2}}\right)\right)}\right)\\

\mathbf{else}:\\
\;\;\;\;t_3\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 4: 50.2% accurate, 0.4× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{c0}{2 \cdot w}\\ t_1 := \frac{c0 \cdot \left(d_m \cdot d_m\right)}{\left(w \cdot h\right) \cdot \left(D_m \cdot D_m\right)}\\ \mathbf{if}\;t_0 \cdot \left(t_1 + \sqrt{t_1 \cdot t_1 - M_m \cdot M_m}\right) \leq \infty:\\ \;\;\;\;t_0 \cdot \frac{2 \cdot \left(c0 \cdot {d_m}^{2}\right)}{\left(w \cdot h\right) \cdot {D_m}^{2}}\\ \mathbf{else}:\\ \;\;\;\;t_0 \cdot \left(c0 \cdot 0\right)\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (/ c0 (* 2.0 w)))
        (t_1 (/ (* c0 (* d_m d_m)) (* (* w h) (* D_m D_m)))))
   (if (<= (* t_0 (+ t_1 (sqrt (- (* t_1 t_1) (* M_m M_m))))) INFINITY)
     (* t_0 (/ (* 2.0 (* c0 (pow d_m 2.0))) (* (* w h) (pow D_m 2.0))))
     (* t_0 (* c0 0.0)))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (2.0 * w);
	double t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	double tmp;
	if ((t_0 * (t_1 + sqrt(((t_1 * t_1) - (M_m * M_m))))) <= ((double) INFINITY)) {
		tmp = t_0 * ((2.0 * (c0 * pow(d_m, 2.0))) / ((w * h) * pow(D_m, 2.0)));
	} else {
		tmp = t_0 * (c0 * 0.0);
	}
	return tmp;
}
D_m = Math.abs(D);
d_m = Math.abs(d);
M_m = Math.abs(M);
public static double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (2.0 * w);
	double t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	double tmp;
	if ((t_0 * (t_1 + Math.sqrt(((t_1 * t_1) - (M_m * M_m))))) <= Double.POSITIVE_INFINITY) {
		tmp = t_0 * ((2.0 * (c0 * Math.pow(d_m, 2.0))) / ((w * h) * Math.pow(D_m, 2.0)));
	} else {
		tmp = t_0 * (c0 * 0.0);
	}
	return tmp;
}
D_m = math.fabs(D)
d_m = math.fabs(d)
M_m = math.fabs(M)
def code(c0, w, h, D_m, d_m, M_m):
	t_0 = c0 / (2.0 * w)
	t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m))
	tmp = 0
	if (t_0 * (t_1 + math.sqrt(((t_1 * t_1) - (M_m * M_m))))) <= math.inf:
		tmp = t_0 * ((2.0 * (c0 * math.pow(d_m, 2.0))) / ((w * h) * math.pow(D_m, 2.0)))
	else:
		tmp = t_0 * (c0 * 0.0)
	return tmp
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(c0 / Float64(2.0 * w))
	t_1 = Float64(Float64(c0 * Float64(d_m * d_m)) / Float64(Float64(w * h) * Float64(D_m * D_m)))
	tmp = 0.0
	if (Float64(t_0 * Float64(t_1 + sqrt(Float64(Float64(t_1 * t_1) - Float64(M_m * M_m))))) <= Inf)
		tmp = Float64(t_0 * Float64(Float64(2.0 * Float64(c0 * (d_m ^ 2.0))) / Float64(Float64(w * h) * (D_m ^ 2.0))));
	else
		tmp = Float64(t_0 * Float64(c0 * 0.0));
	end
	return tmp
end
D_m = abs(D);
d_m = abs(d);
M_m = abs(M);
function tmp_2 = code(c0, w, h, D_m, d_m, M_m)
	t_0 = c0 / (2.0 * w);
	t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	tmp = 0.0;
	if ((t_0 * (t_1 + sqrt(((t_1 * t_1) - (M_m * M_m))))) <= Inf)
		tmp = t_0 * ((2.0 * (c0 * (d_m ^ 2.0))) / ((w * h) * (D_m ^ 2.0)));
	else
		tmp = t_0 * (c0 * 0.0);
	end
	tmp_2 = tmp;
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(c0 * N[(d$95$m * d$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D$95$m * D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[N[(t$95$0 * N[(t$95$1 + N[Sqrt[N[(N[(t$95$1 * t$95$1), $MachinePrecision] - N[(M$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], Infinity], N[(t$95$0 * N[(N[(2.0 * N[(c0 * N[Power[d$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(t$95$0 * N[(c0 * 0.0), $MachinePrecision]), $MachinePrecision]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{c0}{2 \cdot w}\\
t_1 := \frac{c0 \cdot \left(d_m \cdot d_m\right)}{\left(w \cdot h\right) \cdot \left(D_m \cdot D_m\right)}\\
\mathbf{if}\;t_0 \cdot \left(t_1 + \sqrt{t_1 \cdot t_1 - M_m \cdot M_m}\right) \leq \infty:\\
\;\;\;\;t_0 \cdot \frac{2 \cdot \left(c0 \cdot {d_m}^{2}\right)}{\left(w \cdot h\right) \cdot {D_m}^{2}}\\

\mathbf{else}:\\
\;\;\;\;t_0 \cdot \left(c0 \cdot 0\right)\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 5: 50.4% accurate, 0.5× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{c0}{2 \cdot w}\\ t_1 := \frac{c0 \cdot \left(d_m \cdot d_m\right)}{\left(w \cdot h\right) \cdot \left(D_m \cdot D_m\right)}\\ t_2 := t_0 \cdot \left(t_1 + \sqrt{t_1 \cdot t_1 - M_m \cdot M_m}\right)\\ \mathbf{if}\;t_2 \leq \infty:\\ \;\;\;\;t_2\\ \mathbf{else}:\\ \;\;\;\;t_0 \cdot \left(c0 \cdot 0\right)\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (/ c0 (* 2.0 w)))
        (t_1 (/ (* c0 (* d_m d_m)) (* (* w h) (* D_m D_m))))
        (t_2 (* t_0 (+ t_1 (sqrt (- (* t_1 t_1) (* M_m M_m)))))))
   (if (<= t_2 INFINITY) t_2 (* t_0 (* c0 0.0)))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (2.0 * w);
	double t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	double t_2 = t_0 * (t_1 + sqrt(((t_1 * t_1) - (M_m * M_m))));
	double tmp;
	if (t_2 <= ((double) INFINITY)) {
		tmp = t_2;
	} else {
		tmp = t_0 * (c0 * 0.0);
	}
	return tmp;
}
D_m = Math.abs(D);
d_m = Math.abs(d);
M_m = Math.abs(M);
public static double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (2.0 * w);
	double t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	double t_2 = t_0 * (t_1 + Math.sqrt(((t_1 * t_1) - (M_m * M_m))));
	double tmp;
	if (t_2 <= Double.POSITIVE_INFINITY) {
		tmp = t_2;
	} else {
		tmp = t_0 * (c0 * 0.0);
	}
	return tmp;
}
D_m = math.fabs(D)
d_m = math.fabs(d)
M_m = math.fabs(M)
def code(c0, w, h, D_m, d_m, M_m):
	t_0 = c0 / (2.0 * w)
	t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m))
	t_2 = t_0 * (t_1 + math.sqrt(((t_1 * t_1) - (M_m * M_m))))
	tmp = 0
	if t_2 <= math.inf:
		tmp = t_2
	else:
		tmp = t_0 * (c0 * 0.0)
	return tmp
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(c0 / Float64(2.0 * w))
	t_1 = Float64(Float64(c0 * Float64(d_m * d_m)) / Float64(Float64(w * h) * Float64(D_m * D_m)))
	t_2 = Float64(t_0 * Float64(t_1 + sqrt(Float64(Float64(t_1 * t_1) - Float64(M_m * M_m)))))
	tmp = 0.0
	if (t_2 <= Inf)
		tmp = t_2;
	else
		tmp = Float64(t_0 * Float64(c0 * 0.0));
	end
	return tmp
end
D_m = abs(D);
d_m = abs(d);
M_m = abs(M);
function tmp_2 = code(c0, w, h, D_m, d_m, M_m)
	t_0 = c0 / (2.0 * w);
	t_1 = (c0 * (d_m * d_m)) / ((w * h) * (D_m * D_m));
	t_2 = t_0 * (t_1 + sqrt(((t_1 * t_1) - (M_m * M_m))));
	tmp = 0.0;
	if (t_2 <= Inf)
		tmp = t_2;
	else
		tmp = t_0 * (c0 * 0.0);
	end
	tmp_2 = tmp;
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(c0 * N[(d$95$m * d$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(w * h), $MachinePrecision] * N[(D$95$m * D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(t$95$0 * N[(t$95$1 + N[Sqrt[N[(N[(t$95$1 * t$95$1), $MachinePrecision] - N[(M$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$2, Infinity], t$95$2, N[(t$95$0 * N[(c0 * 0.0), $MachinePrecision]), $MachinePrecision]]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{c0}{2 \cdot w}\\
t_1 := \frac{c0 \cdot \left(d_m \cdot d_m\right)}{\left(w \cdot h\right) \cdot \left(D_m \cdot D_m\right)}\\
t_2 := t_0 \cdot \left(t_1 + \sqrt{t_1 \cdot t_1 - M_m \cdot M_m}\right)\\
\mathbf{if}\;t_2 \leq \infty:\\
\;\;\;\;t_2\\

\mathbf{else}:\\
\;\;\;\;t_0 \cdot \left(c0 \cdot 0\right)\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 6: 32.9% accurate, 1.0× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{c0}{w \cdot h}\\ t_1 := t_0 \cdot \frac{d_m \cdot d_m}{D_m \cdot D_m}\\ \mathbf{if}\;c0 \leq -5 \cdot 10^{+154} \lor \neg \left(c0 \leq 4 \cdot 10^{-17}\right):\\ \;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(t_1 + \sqrt{t_1 \cdot \left(t_0 \cdot \left(\frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\right)\right) - M_m \cdot M_m}\right)\\ \mathbf{else}:\\ \;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (/ c0 (* w h))) (t_1 (* t_0 (/ (* d_m d_m) (* D_m D_m)))))
   (if (or (<= c0 -5e+154) (not (<= c0 4e-17)))
     (*
      (/ c0 (* 2.0 w))
      (+
       t_1
       (sqrt (- (* t_1 (* t_0 (* (/ d_m D_m) (/ d_m D_m)))) (* M_m M_m)))))
     (* -0.5 (/ (pow c0 2.0) (/ w 0.0))))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (w * h);
	double t_1 = t_0 * ((d_m * d_m) / (D_m * D_m));
	double tmp;
	if ((c0 <= -5e+154) || !(c0 <= 4e-17)) {
		tmp = (c0 / (2.0 * w)) * (t_1 + sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))));
	} else {
		tmp = -0.5 * (pow(c0, 2.0) / (w / 0.0));
	}
	return tmp;
}
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
real(8) function code(c0, w, h, d_m, d_m_1, m_m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d_m
    real(8), intent (in) :: d_m_1
    real(8), intent (in) :: m_m
    real(8) :: t_0
    real(8) :: t_1
    real(8) :: tmp
    t_0 = c0 / (w * h)
    t_1 = t_0 * ((d_m_1 * d_m_1) / (d_m * d_m))
    if ((c0 <= (-5d+154)) .or. (.not. (c0 <= 4d-17))) then
        tmp = (c0 / (2.0d0 * w)) * (t_1 + sqrt(((t_1 * (t_0 * ((d_m_1 / d_m) * (d_m_1 / d_m)))) - (m_m * m_m))))
    else
        tmp = (-0.5d0) * ((c0 ** 2.0d0) / (w / 0.0d0))
    end if
    code = tmp
end function
D_m = Math.abs(D);
d_m = Math.abs(d);
M_m = Math.abs(M);
public static double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (w * h);
	double t_1 = t_0 * ((d_m * d_m) / (D_m * D_m));
	double tmp;
	if ((c0 <= -5e+154) || !(c0 <= 4e-17)) {
		tmp = (c0 / (2.0 * w)) * (t_1 + Math.sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))));
	} else {
		tmp = -0.5 * (Math.pow(c0, 2.0) / (w / 0.0));
	}
	return tmp;
}
D_m = math.fabs(D)
d_m = math.fabs(d)
M_m = math.fabs(M)
def code(c0, w, h, D_m, d_m, M_m):
	t_0 = c0 / (w * h)
	t_1 = t_0 * ((d_m * d_m) / (D_m * D_m))
	tmp = 0
	if (c0 <= -5e+154) or not (c0 <= 4e-17):
		tmp = (c0 / (2.0 * w)) * (t_1 + math.sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))))
	else:
		tmp = -0.5 * (math.pow(c0, 2.0) / (w / 0.0))
	return tmp
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(c0 / Float64(w * h))
	t_1 = Float64(t_0 * Float64(Float64(d_m * d_m) / Float64(D_m * D_m)))
	tmp = 0.0
	if ((c0 <= -5e+154) || !(c0 <= 4e-17))
		tmp = Float64(Float64(c0 / Float64(2.0 * w)) * Float64(t_1 + sqrt(Float64(Float64(t_1 * Float64(t_0 * Float64(Float64(d_m / D_m) * Float64(d_m / D_m)))) - Float64(M_m * M_m)))));
	else
		tmp = Float64(-0.5 * Float64((c0 ^ 2.0) / Float64(w / 0.0)));
	end
	return tmp
end
D_m = abs(D);
d_m = abs(d);
M_m = abs(M);
function tmp_2 = code(c0, w, h, D_m, d_m, M_m)
	t_0 = c0 / (w * h);
	t_1 = t_0 * ((d_m * d_m) / (D_m * D_m));
	tmp = 0.0;
	if ((c0 <= -5e+154) || ~((c0 <= 4e-17)))
		tmp = (c0 / (2.0 * w)) * (t_1 + sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))));
	else
		tmp = -0.5 * ((c0 ^ 2.0) / (w / 0.0));
	end
	tmp_2 = tmp;
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(c0 / N[(w * h), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(t$95$0 * N[(N[(d$95$m * d$95$m), $MachinePrecision] / N[(D$95$m * D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[Or[LessEqual[c0, -5e+154], N[Not[LessEqual[c0, 4e-17]], $MachinePrecision]], N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(t$95$1 + N[Sqrt[N[(N[(t$95$1 * N[(t$95$0 * N[(N[(d$95$m / D$95$m), $MachinePrecision] * N[(d$95$m / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] - N[(M$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(-0.5 * N[(N[Power[c0, 2.0], $MachinePrecision] / N[(w / 0.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{c0}{w \cdot h}\\
t_1 := t_0 \cdot \frac{d_m \cdot d_m}{D_m \cdot D_m}\\
\mathbf{if}\;c0 \leq -5 \cdot 10^{+154} \lor \neg \left(c0 \leq 4 \cdot 10^{-17}\right):\\
\;\;\;\;\frac{c0}{2 \cdot w} \cdot \left(t_1 + \sqrt{t_1 \cdot \left(t_0 \cdot \left(\frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\right)\right) - M_m \cdot M_m}\right)\\

\mathbf{else}:\\
\;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 7: 32.9% accurate, 1.0× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \begin{array}{l} t_0 := \frac{c0}{w \cdot h}\\ t_1 := t_0 \cdot \frac{d_m \cdot d_m}{D_m \cdot D_m}\\ t_2 := \frac{c0}{2 \cdot w}\\ \mathbf{if}\;c0 \leq -5 \cdot 10^{+154}:\\ \;\;\;\;t_2 \cdot \left(t_1 + \sqrt{t_1 \cdot \left(t_0 \cdot \left(\frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\right)\right) - M_m \cdot M_m}\right)\\ \mathbf{elif}\;c0 \leq 1.02 \cdot 10^{-16}:\\ \;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\ \mathbf{else}:\\ \;\;\;\;t_2 \cdot \left(t_1 + \sqrt{t_1 \cdot \frac{d_m \cdot \left(t_0 \cdot \frac{d_m}{D_m}\right)}{D_m} - M_m \cdot M_m}\right)\\ \end{array} \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (let* ((t_0 (/ c0 (* w h)))
        (t_1 (* t_0 (/ (* d_m d_m) (* D_m D_m))))
        (t_2 (/ c0 (* 2.0 w))))
   (if (<= c0 -5e+154)
     (*
      t_2
      (+
       t_1
       (sqrt (- (* t_1 (* t_0 (* (/ d_m D_m) (/ d_m D_m)))) (* M_m M_m)))))
     (if (<= c0 1.02e-16)
       (* -0.5 (/ (pow c0 2.0) (/ w 0.0)))
       (*
        t_2
        (+
         t_1
         (sqrt
          (- (* t_1 (/ (* d_m (* t_0 (/ d_m D_m))) D_m)) (* M_m M_m)))))))))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (w * h);
	double t_1 = t_0 * ((d_m * d_m) / (D_m * D_m));
	double t_2 = c0 / (2.0 * w);
	double tmp;
	if (c0 <= -5e+154) {
		tmp = t_2 * (t_1 + sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))));
	} else if (c0 <= 1.02e-16) {
		tmp = -0.5 * (pow(c0, 2.0) / (w / 0.0));
	} else {
		tmp = t_2 * (t_1 + sqrt(((t_1 * ((d_m * (t_0 * (d_m / D_m))) / D_m)) - (M_m * M_m))));
	}
	return tmp;
}
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
real(8) function code(c0, w, h, d_m, d_m_1, m_m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d_m
    real(8), intent (in) :: d_m_1
    real(8), intent (in) :: m_m
    real(8) :: t_0
    real(8) :: t_1
    real(8) :: t_2
    real(8) :: tmp
    t_0 = c0 / (w * h)
    t_1 = t_0 * ((d_m_1 * d_m_1) / (d_m * d_m))
    t_2 = c0 / (2.0d0 * w)
    if (c0 <= (-5d+154)) then
        tmp = t_2 * (t_1 + sqrt(((t_1 * (t_0 * ((d_m_1 / d_m) * (d_m_1 / d_m)))) - (m_m * m_m))))
    else if (c0 <= 1.02d-16) then
        tmp = (-0.5d0) * ((c0 ** 2.0d0) / (w / 0.0d0))
    else
        tmp = t_2 * (t_1 + sqrt(((t_1 * ((d_m_1 * (t_0 * (d_m_1 / d_m))) / d_m)) - (m_m * m_m))))
    end if
    code = tmp
end function
D_m = Math.abs(D);
d_m = Math.abs(d);
M_m = Math.abs(M);
public static double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	double t_0 = c0 / (w * h);
	double t_1 = t_0 * ((d_m * d_m) / (D_m * D_m));
	double t_2 = c0 / (2.0 * w);
	double tmp;
	if (c0 <= -5e+154) {
		tmp = t_2 * (t_1 + Math.sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))));
	} else if (c0 <= 1.02e-16) {
		tmp = -0.5 * (Math.pow(c0, 2.0) / (w / 0.0));
	} else {
		tmp = t_2 * (t_1 + Math.sqrt(((t_1 * ((d_m * (t_0 * (d_m / D_m))) / D_m)) - (M_m * M_m))));
	}
	return tmp;
}
D_m = math.fabs(D)
d_m = math.fabs(d)
M_m = math.fabs(M)
def code(c0, w, h, D_m, d_m, M_m):
	t_0 = c0 / (w * h)
	t_1 = t_0 * ((d_m * d_m) / (D_m * D_m))
	t_2 = c0 / (2.0 * w)
	tmp = 0
	if c0 <= -5e+154:
		tmp = t_2 * (t_1 + math.sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))))
	elif c0 <= 1.02e-16:
		tmp = -0.5 * (math.pow(c0, 2.0) / (w / 0.0))
	else:
		tmp = t_2 * (t_1 + math.sqrt(((t_1 * ((d_m * (t_0 * (d_m / D_m))) / D_m)) - (M_m * M_m))))
	return tmp
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	t_0 = Float64(c0 / Float64(w * h))
	t_1 = Float64(t_0 * Float64(Float64(d_m * d_m) / Float64(D_m * D_m)))
	t_2 = Float64(c0 / Float64(2.0 * w))
	tmp = 0.0
	if (c0 <= -5e+154)
		tmp = Float64(t_2 * Float64(t_1 + sqrt(Float64(Float64(t_1 * Float64(t_0 * Float64(Float64(d_m / D_m) * Float64(d_m / D_m)))) - Float64(M_m * M_m)))));
	elseif (c0 <= 1.02e-16)
		tmp = Float64(-0.5 * Float64((c0 ^ 2.0) / Float64(w / 0.0)));
	else
		tmp = Float64(t_2 * Float64(t_1 + sqrt(Float64(Float64(t_1 * Float64(Float64(d_m * Float64(t_0 * Float64(d_m / D_m))) / D_m)) - Float64(M_m * M_m)))));
	end
	return tmp
end
D_m = abs(D);
d_m = abs(d);
M_m = abs(M);
function tmp_2 = code(c0, w, h, D_m, d_m, M_m)
	t_0 = c0 / (w * h);
	t_1 = t_0 * ((d_m * d_m) / (D_m * D_m));
	t_2 = c0 / (2.0 * w);
	tmp = 0.0;
	if (c0 <= -5e+154)
		tmp = t_2 * (t_1 + sqrt(((t_1 * (t_0 * ((d_m / D_m) * (d_m / D_m)))) - (M_m * M_m))));
	elseif (c0 <= 1.02e-16)
		tmp = -0.5 * ((c0 ^ 2.0) / (w / 0.0));
	else
		tmp = t_2 * (t_1 + sqrt(((t_1 * ((d_m * (t_0 * (d_m / D_m))) / D_m)) - (M_m * M_m))));
	end
	tmp_2 = tmp;
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := Block[{t$95$0 = N[(c0 / N[(w * h), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(t$95$0 * N[(N[(d$95$m * d$95$m), $MachinePrecision] / N[(D$95$m * D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[c0, -5e+154], N[(t$95$2 * N[(t$95$1 + N[Sqrt[N[(N[(t$95$1 * N[(t$95$0 * N[(N[(d$95$m / D$95$m), $MachinePrecision] * N[(d$95$m / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] - N[(M$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[c0, 1.02e-16], N[(-0.5 * N[(N[Power[c0, 2.0], $MachinePrecision] / N[(w / 0.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(t$95$2 * N[(t$95$1 + N[Sqrt[N[(N[(t$95$1 * N[(N[(d$95$m * N[(t$95$0 * N[(d$95$m / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / D$95$m), $MachinePrecision]), $MachinePrecision] - N[(M$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\begin{array}{l}
t_0 := \frac{c0}{w \cdot h}\\
t_1 := t_0 \cdot \frac{d_m \cdot d_m}{D_m \cdot D_m}\\
t_2 := \frac{c0}{2 \cdot w}\\
\mathbf{if}\;c0 \leq -5 \cdot 10^{+154}:\\
\;\;\;\;t_2 \cdot \left(t_1 + \sqrt{t_1 \cdot \left(t_0 \cdot \left(\frac{d_m}{D_m} \cdot \frac{d_m}{D_m}\right)\right) - M_m \cdot M_m}\right)\\

\mathbf{elif}\;c0 \leq 1.02 \cdot 10^{-16}:\\
\;\;\;\;-0.5 \cdot \frac{{c0}^{2}}{\frac{w}{0}}\\

\mathbf{else}:\\
\;\;\;\;t_2 \cdot \left(t_1 + \sqrt{t_1 \cdot \frac{d_m \cdot \left(t_0 \cdot \frac{d_m}{D_m}\right)}{D_m} - M_m \cdot M_m}\right)\\


\end{array}
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Alternative 8: 29.5% accurate, 16.8× speedup?

\[\begin{array}{l} D_m = \left|D\right| \\ d_m = \left|d\right| \\ M_m = \left|M\right| \\ \frac{c0}{2 \cdot w} \cdot \left(c0 \cdot 0\right) \end{array} \]
D_m = (fabs.f64 D)
d_m = (fabs.f64 d)
M_m = (fabs.f64 M)
(FPCore (c0 w h D_m d_m M_m)
 :precision binary64
 (* (/ c0 (* 2.0 w)) (* c0 0.0)))
D_m = fabs(D);
d_m = fabs(d);
M_m = fabs(M);
double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	return (c0 / (2.0 * w)) * (c0 * 0.0);
}
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
real(8) function code(c0, w, h, d_m, d_m_1, m_m)
    real(8), intent (in) :: c0
    real(8), intent (in) :: w
    real(8), intent (in) :: h
    real(8), intent (in) :: d_m
    real(8), intent (in) :: d_m_1
    real(8), intent (in) :: m_m
    code = (c0 / (2.0d0 * w)) * (c0 * 0.0d0)
end function
D_m = Math.abs(D);
d_m = Math.abs(d);
M_m = Math.abs(M);
public static double code(double c0, double w, double h, double D_m, double d_m, double M_m) {
	return (c0 / (2.0 * w)) * (c0 * 0.0);
}
D_m = math.fabs(D)
d_m = math.fabs(d)
M_m = math.fabs(M)
def code(c0, w, h, D_m, d_m, M_m):
	return (c0 / (2.0 * w)) * (c0 * 0.0)
D_m = abs(D)
d_m = abs(d)
M_m = abs(M)
function code(c0, w, h, D_m, d_m, M_m)
	return Float64(Float64(c0 / Float64(2.0 * w)) * Float64(c0 * 0.0))
end
D_m = abs(D);
d_m = abs(d);
M_m = abs(M);
function tmp = code(c0, w, h, D_m, d_m, M_m)
	tmp = (c0 / (2.0 * w)) * (c0 * 0.0);
end
D_m = N[Abs[D], $MachinePrecision]
d_m = N[Abs[d], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
code[c0_, w_, h_, D$95$m_, d$95$m_, M$95$m_] := N[(N[(c0 / N[(2.0 * w), $MachinePrecision]), $MachinePrecision] * N[(c0 * 0.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
d_m = \left|d\right|
\\
M_m = \left|M\right|

\\
\frac{c0}{2 \cdot w} \cdot \left(c0 \cdot 0\right)
\end{array}
Derivation
    &prev;&pcontext;&pcontext2;&ctx;
  1. Add Preprocessing

Reproduce

?
herbie shell --seed 2024010 
(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))))))