Henrywood and Agarwal, Equation (9a)

Percentage Accurate: 81.5% → 89.7%
Time: 11.1s
Alternatives: 12
Speedup: 2.1×

Specification

?
\[\begin{array}{l} \\ w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \end{array} \]
(FPCore (w0 M D h l d)
 :precision binary64
 (* w0 (sqrt (- 1.0 (* (pow (/ (* M D) (* 2.0 d)) 2.0) (/ h l))))))
double code(double w0, double M, double D, double h, double l, double d) {
	return w0 * sqrt((1.0 - (pow(((M * D) / (2.0 * d)), 2.0) * (h / l))));
}
real(8) function code(w0, m, d, h, l, d_1)
    real(8), intent (in) :: w0
    real(8), intent (in) :: m
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: d_1
    code = w0 * sqrt((1.0d0 - ((((m * d) / (2.0d0 * d_1)) ** 2.0d0) * (h / l))))
end function
public static double code(double w0, double M, double D, double h, double l, double d) {
	return w0 * Math.sqrt((1.0 - (Math.pow(((M * D) / (2.0 * d)), 2.0) * (h / l))));
}
def code(w0, M, D, h, l, d):
	return w0 * math.sqrt((1.0 - (math.pow(((M * D) / (2.0 * d)), 2.0) * (h / l))))
function code(w0, M, D, h, l, d)
	return Float64(w0 * sqrt(Float64(1.0 - Float64((Float64(Float64(M * D) / Float64(2.0 * d)) ^ 2.0) * Float64(h / l)))))
end
function tmp = code(w0, M, D, h, l, d)
	tmp = w0 * sqrt((1.0 - ((((M * D) / (2.0 * d)) ^ 2.0) * (h / l))));
end
code[w0_, M_, D_, h_, l_, d_] := N[(w0 * N[Sqrt[N[(1.0 - N[(N[Power[N[(N[(M * D), $MachinePrecision] / N[(2.0 * d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]
\begin{array}{l}

\\
w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}
\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 12 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: 81.5% accurate, 1.0× speedup?

\[\begin{array}{l} \\ w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \end{array} \]
(FPCore (w0 M D h l d)
 :precision binary64
 (* w0 (sqrt (- 1.0 (* (pow (/ (* M D) (* 2.0 d)) 2.0) (/ h l))))))
double code(double w0, double M, double D, double h, double l, double d) {
	return w0 * sqrt((1.0 - (pow(((M * D) / (2.0 * d)), 2.0) * (h / l))));
}
real(8) function code(w0, m, d, h, l, d_1)
    real(8), intent (in) :: w0
    real(8), intent (in) :: m
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: d_1
    code = w0 * sqrt((1.0d0 - ((((m * d) / (2.0d0 * d_1)) ** 2.0d0) * (h / l))))
end function
public static double code(double w0, double M, double D, double h, double l, double d) {
	return w0 * Math.sqrt((1.0 - (Math.pow(((M * D) / (2.0 * d)), 2.0) * (h / l))));
}
def code(w0, M, D, h, l, d):
	return w0 * math.sqrt((1.0 - (math.pow(((M * D) / (2.0 * d)), 2.0) * (h / l))))
function code(w0, M, D, h, l, d)
	return Float64(w0 * sqrt(Float64(1.0 - Float64((Float64(Float64(M * D) / Float64(2.0 * d)) ^ 2.0) * Float64(h / l)))))
end
function tmp = code(w0, M, D, h, l, d)
	tmp = w0 * sqrt((1.0 - ((((M * D) / (2.0 * d)) ^ 2.0) * (h / l))));
end
code[w0_, M_, D_, h_, l_, d_] := N[(w0 * N[Sqrt[N[(1.0 - N[(N[Power[N[(N[(M * D), $MachinePrecision] / N[(2.0 * d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]
\begin{array}{l}

\\
w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}
\end{array}

Alternative 1: 89.7% accurate, 0.8× speedup?

\[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;\frac{D\_m \cdot M\_m}{d\_m \cdot 2} \leq 2 \cdot 10^{+104}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{{\left(\frac{\frac{d\_m}{M\_m}}{D\_m} \cdot 2\right)}^{-2}}{\ell}, -h, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d\_m} \cdot M\_m\right) \cdot D\_m, \frac{\left(\left(0.5 \cdot D\_m\right) \cdot h\right) \cdot \frac{M\_m}{d\_m}}{-\ell}, 1\right)} \cdot w0\\ \end{array} \end{array} \]
d_m = (fabs.f64 d)
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
(FPCore (w0 M_m D_m h l d_m)
 :precision binary64
 (if (<= (/ (* D_m M_m) (* d_m 2.0)) 2e+104)
   (* (sqrt (fma (/ (pow (* (/ (/ d_m M_m) D_m) 2.0) -2.0) l) (- h) 1.0)) w0)
   (*
    (sqrt
     (fma
      (* (* (/ 0.5 d_m) M_m) D_m)
      (/ (* (* (* 0.5 D_m) h) (/ M_m d_m)) (- l))
      1.0))
    w0)))
d_m = fabs(d);
D_m = fabs(D);
M_m = fabs(M);
assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
	double tmp;
	if (((D_m * M_m) / (d_m * 2.0)) <= 2e+104) {
		tmp = sqrt(fma((pow((((d_m / M_m) / D_m) * 2.0), -2.0) / l), -h, 1.0)) * w0;
	} else {
		tmp = sqrt(fma((((0.5 / d_m) * M_m) * D_m), ((((0.5 * D_m) * h) * (M_m / d_m)) / -l), 1.0)) * w0;
	}
	return tmp;
}
d_m = abs(d)
D_m = abs(D)
M_m = abs(M)
w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
function code(w0, M_m, D_m, h, l, d_m)
	tmp = 0.0
	if (Float64(Float64(D_m * M_m) / Float64(d_m * 2.0)) <= 2e+104)
		tmp = Float64(sqrt(fma(Float64((Float64(Float64(Float64(d_m / M_m) / D_m) * 2.0) ^ -2.0) / l), Float64(-h), 1.0)) * w0);
	else
		tmp = Float64(sqrt(fma(Float64(Float64(Float64(0.5 / d_m) * M_m) * D_m), Float64(Float64(Float64(Float64(0.5 * D_m) * h) * Float64(M_m / d_m)) / Float64(-l)), 1.0)) * w0);
	end
	return tmp
end
d_m = N[Abs[d], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision], 2e+104], N[(N[Sqrt[N[(N[(N[Power[N[(N[(N[(d$95$m / M$95$m), $MachinePrecision] / D$95$m), $MachinePrecision] * 2.0), $MachinePrecision], -2.0], $MachinePrecision] / l), $MachinePrecision] * (-h) + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision], N[(N[Sqrt[N[(N[(N[(N[(0.5 / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] * D$95$m), $MachinePrecision] * N[(N[(N[(N[(0.5 * D$95$m), $MachinePrecision] * h), $MachinePrecision] * N[(M$95$m / d$95$m), $MachinePrecision]), $MachinePrecision] / (-l)), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision]]
\begin{array}{l}
d_m = \left|d\right|
\\
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
[w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
\\
\begin{array}{l}
\mathbf{if}\;\frac{D\_m \cdot M\_m}{d\_m \cdot 2} \leq 2 \cdot 10^{+104}:\\
\;\;\;\;\sqrt{\mathsf{fma}\left(\frac{{\left(\frac{\frac{d\_m}{M\_m}}{D\_m} \cdot 2\right)}^{-2}}{\ell}, -h, 1\right)} \cdot w0\\

\mathbf{else}:\\
\;\;\;\;\sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d\_m} \cdot M\_m\right) \cdot D\_m, \frac{\left(\left(0.5 \cdot D\_m\right) \cdot h\right) \cdot \frac{M\_m}{d\_m}}{-\ell}, 1\right)} \cdot w0\\


\end{array}
\end{array}
Derivation
  1. Split input into 2 regimes
  2. if (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) < 2e104

    1. Initial program 81.2%

      \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
    2. Add Preprocessing
    3. Step-by-step derivation
      1. lift--.f64N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}} \]
      2. sub-negN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right)}} \]
      3. +-commutativeN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) + 1}} \]
      4. lift-*.f64N/A

        \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}\right)\right) + 1} \]
      5. distribute-lft-neg-inN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)\right) \cdot \frac{h}{\ell}} + 1} \]
      6. lift-/.f64N/A

        \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)\right) \cdot \color{blue}{\frac{h}{\ell}} + 1} \]
      7. clear-numN/A

        \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)\right) \cdot \color{blue}{\frac{1}{\frac{\ell}{h}}} + 1} \]
      8. un-div-invN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)}{\frac{\ell}{h}}} + 1} \]
      9. frac-2negN/A

        \[\leadsto w0 \cdot \sqrt{\frac{\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)}{\color{blue}{\frac{\mathsf{neg}\left(\ell\right)}{\mathsf{neg}\left(h\right)}}} + 1} \]
      10. associate-/r/N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)}{\mathsf{neg}\left(\ell\right)} \cdot \left(\mathsf{neg}\left(h\right)\right)} + 1} \]
      11. neg-mul-1N/A

        \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{-1 \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}}}{\mathsf{neg}\left(\ell\right)} \cdot \left(\mathsf{neg}\left(h\right)\right) + 1} \]
      12. associate-*l/N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{\mathsf{neg}\left(\ell\right)} \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)} \cdot \left(\mathsf{neg}\left(h\right)\right) + 1} \]
      13. metadata-evalN/A

        \[\leadsto w0 \cdot \sqrt{\left(\frac{\color{blue}{\mathsf{neg}\left(1\right)}}{\mathsf{neg}\left(\ell\right)} \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right) \cdot \left(\mathsf{neg}\left(h\right)\right) + 1} \]
      14. frac-2negN/A

        \[\leadsto w0 \cdot \sqrt{\left(\color{blue}{\frac{1}{\ell}} \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right) \cdot \left(\mathsf{neg}\left(h\right)\right) + 1} \]
    4. Applied rewrites91.3%

      \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\frac{{\left(\frac{\frac{d}{M}}{D} \cdot 2\right)}^{-2}}{\ell}, -h, 1\right)}} \]

    if 2e104 < (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d))

    1. Initial program 49.8%

      \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
    2. Add Preprocessing
    3. Step-by-step derivation
      1. lift--.f64N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}} \]
      2. sub-negN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right)}} \]
      3. +-commutativeN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) + 1}} \]
      4. lift-*.f64N/A

        \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}\right)\right) + 1} \]
      5. lift-/.f64N/A

        \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \color{blue}{\frac{h}{\ell}}\right)\right) + 1} \]
      6. associate-*r/N/A

        \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{\frac{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot h}{\ell}}\right)\right) + 1} \]
      7. distribute-neg-frac2N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot h}{\mathsf{neg}\left(\ell\right)}} + 1} \]
      8. lift-pow.f64N/A

        \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}} \cdot h}{\mathsf{neg}\left(\ell\right)} + 1} \]
      9. unpow2N/A

        \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{M \cdot D}{2 \cdot d} \cdot \frac{M \cdot D}{2 \cdot d}\right)} \cdot h}{\mathsf{neg}\left(\ell\right)} + 1} \]
      10. associate-*l*N/A

        \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\frac{M \cdot D}{2 \cdot d} \cdot \left(\frac{M \cdot D}{2 \cdot d} \cdot h\right)}}{\mathsf{neg}\left(\ell\right)} + 1} \]
      11. associate-/l*N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{M \cdot D}{2 \cdot d} \cdot \frac{\frac{M \cdot D}{2 \cdot d} \cdot h}{\mathsf{neg}\left(\ell\right)}} + 1} \]
      12. lower-fma.f64N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\frac{M \cdot D}{2 \cdot d}, \frac{\frac{M \cdot D}{2 \cdot d} \cdot h}{\mathsf{neg}\left(\ell\right)}, 1\right)}} \]
    4. Applied rewrites60.3%

      \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\left(\frac{0.5}{d} \cdot M\right) \cdot D, \frac{\left(D \cdot 0.5\right) \cdot \left(\frac{M}{d} \cdot h\right)}{-\ell}, 1\right)}} \]
    5. Step-by-step derivation
      1. lift-*.f64N/A

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{\frac{1}{2}}{d} \cdot M\right) \cdot D, \frac{\color{blue}{\left(D \cdot \frac{1}{2}\right) \cdot \left(\frac{M}{d} \cdot h\right)}}{-\ell}, 1\right)} \]
      2. *-commutativeN/A

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{\frac{1}{2}}{d} \cdot M\right) \cdot D, \frac{\color{blue}{\left(\frac{M}{d} \cdot h\right) \cdot \left(D \cdot \frac{1}{2}\right)}}{-\ell}, 1\right)} \]
      3. lift-*.f64N/A

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{\frac{1}{2}}{d} \cdot M\right) \cdot D, \frac{\color{blue}{\left(\frac{M}{d} \cdot h\right)} \cdot \left(D \cdot \frac{1}{2}\right)}{-\ell}, 1\right)} \]
      4. associate-*l*N/A

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{\frac{1}{2}}{d} \cdot M\right) \cdot D, \frac{\color{blue}{\frac{M}{d} \cdot \left(h \cdot \left(D \cdot \frac{1}{2}\right)\right)}}{-\ell}, 1\right)} \]
      5. lower-*.f64N/A

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{\frac{1}{2}}{d} \cdot M\right) \cdot D, \frac{\color{blue}{\frac{M}{d} \cdot \left(h \cdot \left(D \cdot \frac{1}{2}\right)\right)}}{-\ell}, 1\right)} \]
      6. lower-*.f6460.2

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d} \cdot M\right) \cdot D, \frac{\frac{M}{d} \cdot \color{blue}{\left(h \cdot \left(D \cdot 0.5\right)\right)}}{-\ell}, 1\right)} \]
    6. Applied rewrites60.2%

      \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d} \cdot M\right) \cdot D, \frac{\color{blue}{\frac{M}{d} \cdot \left(h \cdot \left(D \cdot 0.5\right)\right)}}{-\ell}, 1\right)} \]
  3. Recombined 2 regimes into one program.
  4. Final simplification86.7%

    \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{D \cdot M}{d \cdot 2} \leq 2 \cdot 10^{+104}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{{\left(\frac{\frac{d}{M}}{D} \cdot 2\right)}^{-2}}{\ell}, -h, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d} \cdot M\right) \cdot D, \frac{\left(\left(0.5 \cdot D\right) \cdot h\right) \cdot \frac{M}{d}}{-\ell}, 1\right)} \cdot w0\\ \end{array} \]
  5. Add Preprocessing

Alternative 2: 83.6% accurate, 0.7× speedup?

\[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -2 \cdot 10^{-7}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{M\_m}{\ell \cdot d\_m} \cdot \frac{\left(D\_m \cdot D\_m\right) \cdot M\_m}{d\_m}, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;1 \cdot w0\\ \end{array} \end{array} \]
d_m = (fabs.f64 d)
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
(FPCore (w0 M_m D_m h l d_m)
 :precision binary64
 (if (<= (* (pow (/ (* D_m M_m) (* d_m 2.0)) 2.0) (/ h l)) -2e-7)
   (*
    (sqrt
     (fma (* -0.25 h) (* (/ M_m (* l d_m)) (/ (* (* D_m D_m) M_m) d_m)) 1.0))
    w0)
   (* 1.0 w0)))
d_m = fabs(d);
D_m = fabs(D);
M_m = fabs(M);
assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
	double tmp;
	if ((pow(((D_m * M_m) / (d_m * 2.0)), 2.0) * (h / l)) <= -2e-7) {
		tmp = sqrt(fma((-0.25 * h), ((M_m / (l * d_m)) * (((D_m * D_m) * M_m) / d_m)), 1.0)) * w0;
	} else {
		tmp = 1.0 * w0;
	}
	return tmp;
}
d_m = abs(d)
D_m = abs(D)
M_m = abs(M)
w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
function code(w0, M_m, D_m, h, l, d_m)
	tmp = 0.0
	if (Float64((Float64(Float64(D_m * M_m) / Float64(d_m * 2.0)) ^ 2.0) * Float64(h / l)) <= -2e-7)
		tmp = Float64(sqrt(fma(Float64(-0.25 * h), Float64(Float64(M_m / Float64(l * d_m)) * Float64(Float64(Float64(D_m * D_m) * M_m) / d_m)), 1.0)) * w0);
	else
		tmp = Float64(1.0 * w0);
	end
	return tmp
end
d_m = N[Abs[d], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[(N[Power[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision], -2e-7], N[(N[Sqrt[N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(M$95$m / N[(l * d$95$m), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(D$95$m * D$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] / d$95$m), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision], N[(1.0 * w0), $MachinePrecision]]
\begin{array}{l}
d_m = \left|d\right|
\\
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
[w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
\\
\begin{array}{l}
\mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -2 \cdot 10^{-7}:\\
\;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{M\_m}{\ell \cdot d\_m} \cdot \frac{\left(D\_m \cdot D\_m\right) \cdot M\_m}{d\_m}, 1\right)} \cdot w0\\

\mathbf{else}:\\
\;\;\;\;1 \cdot w0\\


\end{array}
\end{array}
Derivation
  1. Split input into 2 regimes
  2. if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -1.9999999999999999e-7

    1. Initial program 57.5%

      \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
    2. Add Preprocessing
    3. Taylor expanded in M around 0

      \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
    4. Step-by-step derivation
      1. +-commutativeN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
      2. associate-*r/N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
      3. associate-*r*N/A

        \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
      4. associate-*r*N/A

        \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
      5. associate-*l/N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
      6. associate-*r/N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
      7. lft-mult-inverseN/A

        \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
      8. distribute-rgt-inN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
      9. distribute-lft-inN/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
      10. associate-*r*N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
      11. rgt-mult-inverseN/A

        \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
      12. lower-fma.f64N/A

        \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
    5. Applied rewrites44.0%

      \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
    6. Step-by-step derivation
      1. Applied rewrites53.9%

        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right) \cdot D}{\color{blue}{\ell}}, 1\right)} \]
      2. Step-by-step derivation
        1. Applied rewrites50.8%

          \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(D \cdot D\right) \cdot M}{d} \cdot \color{blue}{\frac{M}{\ell \cdot d}}, 1\right)} \]

        if -1.9999999999999999e-7 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l))

        1. Initial program 84.8%

          \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
        2. Add Preprocessing
        3. Taylor expanded in M around 0

          \[\leadsto w0 \cdot \color{blue}{1} \]
        4. Step-by-step derivation
          1. Applied rewrites94.1%

            \[\leadsto w0 \cdot \color{blue}{1} \]
        5. Recombined 2 regimes into one program.
        6. Final simplification80.9%

          \[\leadsto \begin{array}{l} \mathbf{if}\;{\left(\frac{D \cdot M}{d \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -2 \cdot 10^{-7}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{M}{\ell \cdot d} \cdot \frac{\left(D \cdot D\right) \cdot M}{d}, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;1 \cdot w0\\ \end{array} \]
        7. Add Preprocessing

        Alternative 3: 82.4% accurate, 0.8× speedup?

        \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -2 \cdot 10^{-7}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \left(\frac{D\_m}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot \left(D\_m \cdot M\_m\right)\right) \cdot M\_m, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;1 \cdot w0\\ \end{array} \end{array} \]
        d_m = (fabs.f64 d)
        D_m = (fabs.f64 D)
        M_m = (fabs.f64 M)
        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
        (FPCore (w0 M_m D_m h l d_m)
         :precision binary64
         (if (<= (* (pow (/ (* D_m M_m) (* d_m 2.0)) 2.0) (/ h l)) -2e-7)
           (*
            (sqrt
             (fma (* -0.25 h) (* (* (/ D_m (* (* d_m d_m) l)) (* D_m M_m)) M_m) 1.0))
            w0)
           (* 1.0 w0)))
        d_m = fabs(d);
        D_m = fabs(D);
        M_m = fabs(M);
        assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
        double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
        	double tmp;
        	if ((pow(((D_m * M_m) / (d_m * 2.0)), 2.0) * (h / l)) <= -2e-7) {
        		tmp = sqrt(fma((-0.25 * h), (((D_m / ((d_m * d_m) * l)) * (D_m * M_m)) * M_m), 1.0)) * w0;
        	} else {
        		tmp = 1.0 * w0;
        	}
        	return tmp;
        }
        
        d_m = abs(d)
        D_m = abs(D)
        M_m = abs(M)
        w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
        function code(w0, M_m, D_m, h, l, d_m)
        	tmp = 0.0
        	if (Float64((Float64(Float64(D_m * M_m) / Float64(d_m * 2.0)) ^ 2.0) * Float64(h / l)) <= -2e-7)
        		tmp = Float64(sqrt(fma(Float64(-0.25 * h), Float64(Float64(Float64(D_m / Float64(Float64(d_m * d_m) * l)) * Float64(D_m * M_m)) * M_m), 1.0)) * w0);
        	else
        		tmp = Float64(1.0 * w0);
        	end
        	return tmp
        end
        
        d_m = N[Abs[d], $MachinePrecision]
        D_m = N[Abs[D], $MachinePrecision]
        M_m = N[Abs[M], $MachinePrecision]
        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
        code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[(N[Power[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision], -2e-7], N[(N[Sqrt[N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(D$95$m / N[(N[(d$95$m * d$95$m), $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision] * N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision] * M$95$m), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision], N[(1.0 * w0), $MachinePrecision]]
        
        \begin{array}{l}
        d_m = \left|d\right|
        \\
        D_m = \left|D\right|
        \\
        M_m = \left|M\right|
        \\
        [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
        \\
        \begin{array}{l}
        \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -2 \cdot 10^{-7}:\\
        \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \left(\frac{D\_m}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot \left(D\_m \cdot M\_m\right)\right) \cdot M\_m, 1\right)} \cdot w0\\
        
        \mathbf{else}:\\
        \;\;\;\;1 \cdot w0\\
        
        
        \end{array}
        \end{array}
        
        Derivation
        1. Split input into 2 regimes
        2. if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -1.9999999999999999e-7

          1. Initial program 57.5%

            \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
          2. Add Preprocessing
          3. Taylor expanded in M around 0

            \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
          4. Step-by-step derivation
            1. +-commutativeN/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
            2. associate-*r/N/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
            3. associate-*r*N/A

              \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
            4. associate-*r*N/A

              \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
            5. associate-*l/N/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
            6. associate-*r/N/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
            7. lft-mult-inverseN/A

              \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
            8. distribute-rgt-inN/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
            9. distribute-lft-inN/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
            10. associate-*r*N/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
            11. rgt-mult-inverseN/A

              \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
            12. lower-fma.f64N/A

              \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
          5. Applied rewrites44.0%

            \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
          6. Step-by-step derivation
            1. Applied rewrites53.7%

              \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, M \cdot \color{blue}{\left(\left(D \cdot M\right) \cdot \frac{D}{\left(d \cdot d\right) \cdot \ell}\right)}, 1\right)} \]

            if -1.9999999999999999e-7 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l))

            1. Initial program 84.8%

              \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
            2. Add Preprocessing
            3. Taylor expanded in M around 0

              \[\leadsto w0 \cdot \color{blue}{1} \]
            4. Step-by-step derivation
              1. Applied rewrites94.1%

                \[\leadsto w0 \cdot \color{blue}{1} \]
            5. Recombined 2 regimes into one program.
            6. Final simplification81.8%

              \[\leadsto \begin{array}{l} \mathbf{if}\;{\left(\frac{D \cdot M}{d \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -2 \cdot 10^{-7}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \left(\frac{D}{\left(d \cdot d\right) \cdot \ell} \cdot \left(D \cdot M\right)\right) \cdot M, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;1 \cdot w0\\ \end{array} \]
            7. Add Preprocessing

            Alternative 4: 82.9% accurate, 0.8× speedup?

            \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \leq 2 \cdot 10^{-285}:\\ \;\;\;\;1 \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\left(\frac{M\_m}{d\_m} \cdot M\_m\right) \cdot D\_m\right) \cdot D\_m}{\ell \cdot d\_m}, 1\right)} \cdot w0\\ \end{array} \end{array} \]
            d_m = (fabs.f64 d)
            D_m = (fabs.f64 D)
            M_m = (fabs.f64 M)
            NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
            (FPCore (w0 M_m D_m h l d_m)
             :precision binary64
             (if (<= (pow (/ (* D_m M_m) (* d_m 2.0)) 2.0) 2e-285)
               (* 1.0 w0)
               (*
                (sqrt
                 (fma (* -0.25 h) (/ (* (* (* (/ M_m d_m) M_m) D_m) D_m) (* l d_m)) 1.0))
                w0)))
            d_m = fabs(d);
            D_m = fabs(D);
            M_m = fabs(M);
            assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
            double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
            	double tmp;
            	if (pow(((D_m * M_m) / (d_m * 2.0)), 2.0) <= 2e-285) {
            		tmp = 1.0 * w0;
            	} else {
            		tmp = sqrt(fma((-0.25 * h), (((((M_m / d_m) * M_m) * D_m) * D_m) / (l * d_m)), 1.0)) * w0;
            	}
            	return tmp;
            }
            
            d_m = abs(d)
            D_m = abs(D)
            M_m = abs(M)
            w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
            function code(w0, M_m, D_m, h, l, d_m)
            	tmp = 0.0
            	if ((Float64(Float64(D_m * M_m) / Float64(d_m * 2.0)) ^ 2.0) <= 2e-285)
            		tmp = Float64(1.0 * w0);
            	else
            		tmp = Float64(sqrt(fma(Float64(-0.25 * h), Float64(Float64(Float64(Float64(Float64(M_m / d_m) * M_m) * D_m) * D_m) / Float64(l * d_m)), 1.0)) * w0);
            	end
            	return tmp
            end
            
            d_m = N[Abs[d], $MachinePrecision]
            D_m = N[Abs[D], $MachinePrecision]
            M_m = N[Abs[M], $MachinePrecision]
            NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
            code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[Power[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision], 2e-285], N[(1.0 * w0), $MachinePrecision], N[(N[Sqrt[N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(N[(N[(M$95$m / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] * D$95$m), $MachinePrecision] * D$95$m), $MachinePrecision] / N[(l * d$95$m), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision]]
            
            \begin{array}{l}
            d_m = \left|d\right|
            \\
            D_m = \left|D\right|
            \\
            M_m = \left|M\right|
            \\
            [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
            \\
            \begin{array}{l}
            \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \leq 2 \cdot 10^{-285}:\\
            \;\;\;\;1 \cdot w0\\
            
            \mathbf{else}:\\
            \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\left(\frac{M\_m}{d\_m} \cdot M\_m\right) \cdot D\_m\right) \cdot D\_m}{\ell \cdot d\_m}, 1\right)} \cdot w0\\
            
            
            \end{array}
            \end{array}
            
            Derivation
            1. Split input into 2 regimes
            2. if (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) < 2.00000000000000015e-285

              1. Initial program 83.3%

                \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
              2. Add Preprocessing
              3. Taylor expanded in M around 0

                \[\leadsto w0 \cdot \color{blue}{1} \]
              4. Step-by-step derivation
                1. Applied rewrites97.4%

                  \[\leadsto w0 \cdot \color{blue}{1} \]

                if 2.00000000000000015e-285 < (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64))

                1. Initial program 71.5%

                  \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                2. Add Preprocessing
                3. Taylor expanded in M around 0

                  \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
                4. Step-by-step derivation
                  1. +-commutativeN/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
                  2. associate-*r/N/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
                  3. associate-*r*N/A

                    \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
                  4. associate-*r*N/A

                    \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
                  5. associate-*l/N/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
                  6. associate-*r/N/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
                  7. lft-mult-inverseN/A

                    \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
                  8. distribute-rgt-inN/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                  9. distribute-lft-inN/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
                  10. associate-*r*N/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
                  11. rgt-mult-inverseN/A

                    \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
                  12. lower-fma.f64N/A

                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
                5. Applied rewrites45.8%

                  \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
                6. Step-by-step derivation
                  1. Applied rewrites66.5%

                    \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right) \cdot D}{\color{blue}{\ell}}, 1\right)} \]
                  2. Step-by-step derivation
                    1. Applied rewrites63.8%

                      \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot D\right) \cdot \left(-D\right)}{\color{blue}{\left(-\ell\right) \cdot d}}, 1\right)} \]
                  3. Recombined 2 regimes into one program.
                  4. Final simplification78.0%

                    \[\leadsto \begin{array}{l} \mathbf{if}\;{\left(\frac{D \cdot M}{d \cdot 2}\right)}^{2} \leq 2 \cdot 10^{-285}:\\ \;\;\;\;1 \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot D\right) \cdot D}{\ell \cdot d}, 1\right)} \cdot w0\\ \end{array} \]
                  5. Add Preprocessing

                  Alternative 5: 79.8% accurate, 0.8× speedup?

                  \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -4 \cdot 10^{+37}:\\ \;\;\;\;\mathsf{fma}\left(-0.125 \cdot \left(D\_m \cdot D\_m\right), \left(\left(\frac{w0}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot h\right) \cdot M\_m\right) \cdot M\_m, w0\right)\\ \mathbf{else}:\\ \;\;\;\;1 \cdot w0\\ \end{array} \end{array} \]
                  d_m = (fabs.f64 d)
                  D_m = (fabs.f64 D)
                  M_m = (fabs.f64 M)
                  NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                  (FPCore (w0 M_m D_m h l d_m)
                   :precision binary64
                   (if (<= (* (pow (/ (* D_m M_m) (* d_m 2.0)) 2.0) (/ h l)) -4e+37)
                     (fma
                      (* -0.125 (* D_m D_m))
                      (* (* (* (/ w0 (* (* d_m d_m) l)) h) M_m) M_m)
                      w0)
                     (* 1.0 w0)))
                  d_m = fabs(d);
                  D_m = fabs(D);
                  M_m = fabs(M);
                  assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                  double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                  	double tmp;
                  	if ((pow(((D_m * M_m) / (d_m * 2.0)), 2.0) * (h / l)) <= -4e+37) {
                  		tmp = fma((-0.125 * (D_m * D_m)), ((((w0 / ((d_m * d_m) * l)) * h) * M_m) * M_m), w0);
                  	} else {
                  		tmp = 1.0 * w0;
                  	}
                  	return tmp;
                  }
                  
                  d_m = abs(d)
                  D_m = abs(D)
                  M_m = abs(M)
                  w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                  function code(w0, M_m, D_m, h, l, d_m)
                  	tmp = 0.0
                  	if (Float64((Float64(Float64(D_m * M_m) / Float64(d_m * 2.0)) ^ 2.0) * Float64(h / l)) <= -4e+37)
                  		tmp = fma(Float64(-0.125 * Float64(D_m * D_m)), Float64(Float64(Float64(Float64(w0 / Float64(Float64(d_m * d_m) * l)) * h) * M_m) * M_m), w0);
                  	else
                  		tmp = Float64(1.0 * w0);
                  	end
                  	return tmp
                  end
                  
                  d_m = N[Abs[d], $MachinePrecision]
                  D_m = N[Abs[D], $MachinePrecision]
                  M_m = N[Abs[M], $MachinePrecision]
                  NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                  code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[(N[Power[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision], -4e+37], N[(N[(-0.125 * N[(D$95$m * D$95$m), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(N[(w0 / N[(N[(d$95$m * d$95$m), $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision] * h), $MachinePrecision] * M$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] + w0), $MachinePrecision], N[(1.0 * w0), $MachinePrecision]]
                  
                  \begin{array}{l}
                  d_m = \left|d\right|
                  \\
                  D_m = \left|D\right|
                  \\
                  M_m = \left|M\right|
                  \\
                  [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                  \\
                  \begin{array}{l}
                  \mathbf{if}\;{\left(\frac{D\_m \cdot M\_m}{d\_m \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -4 \cdot 10^{+37}:\\
                  \;\;\;\;\mathsf{fma}\left(-0.125 \cdot \left(D\_m \cdot D\_m\right), \left(\left(\frac{w0}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot h\right) \cdot M\_m\right) \cdot M\_m, w0\right)\\
                  
                  \mathbf{else}:\\
                  \;\;\;\;1 \cdot w0\\
                  
                  
                  \end{array}
                  \end{array}
                  
                  Derivation
                  1. Split input into 2 regimes
                  2. if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -3.99999999999999982e37

                    1. Initial program 55.3%

                      \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                    2. Add Preprocessing
                    3. Taylor expanded in M around 0

                      \[\leadsto \color{blue}{w0 + \frac{-1}{8} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell}} \]
                    4. Step-by-step derivation
                      1. +-commutativeN/A

                        \[\leadsto \color{blue}{\frac{-1}{8} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell} + w0} \]
                      2. *-commutativeN/A

                        \[\leadsto \color{blue}{\frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell} \cdot \frac{-1}{8}} + w0 \]
                      3. associate-/l*N/A

                        \[\leadsto \color{blue}{\left({D}^{2} \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}\right)} \cdot \frac{-1}{8} + w0 \]
                      4. associate-*r*N/A

                        \[\leadsto \color{blue}{{D}^{2} \cdot \left(\frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell} \cdot \frac{-1}{8}\right)} + w0 \]
                      5. *-commutativeN/A

                        \[\leadsto {D}^{2} \cdot \color{blue}{\left(\frac{-1}{8} \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}\right)} + w0 \]
                      6. associate-*r*N/A

                        \[\leadsto \color{blue}{\left({D}^{2} \cdot \frac{-1}{8}\right) \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}} + w0 \]
                      7. lower-fma.f64N/A

                        \[\leadsto \color{blue}{\mathsf{fma}\left({D}^{2} \cdot \frac{-1}{8}, \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}, w0\right)} \]
                    5. Applied rewrites35.9%

                      \[\leadsto \color{blue}{\mathsf{fma}\left(\left(D \cdot D\right) \cdot -0.125, \frac{\left(M \cdot M\right) \cdot h}{\ell} \cdot \frac{w0}{d \cdot d}, w0\right)} \]
                    6. Step-by-step derivation
                      1. Applied rewrites35.8%

                        \[\leadsto \mathsf{fma}\left(\left(D \cdot D\right) \cdot -0.125, \left(\left(M \cdot M\right) \cdot h\right) \cdot \color{blue}{\frac{w0}{\left(d \cdot d\right) \cdot \ell}}, w0\right) \]
                      2. Step-by-step derivation
                        1. Applied rewrites42.1%

                          \[\leadsto \mathsf{fma}\left(\left(D \cdot D\right) \cdot -0.125, M \cdot \color{blue}{\left(M \cdot \left(\frac{w0}{\left(d \cdot d\right) \cdot \ell} \cdot h\right)\right)}, w0\right) \]

                        if -3.99999999999999982e37 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l))

                        1. Initial program 85.1%

                          \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                        2. Add Preprocessing
                        3. Taylor expanded in M around 0

                          \[\leadsto w0 \cdot \color{blue}{1} \]
                        4. Step-by-step derivation
                          1. Applied rewrites92.5%

                            \[\leadsto w0 \cdot \color{blue}{1} \]
                        5. Recombined 2 regimes into one program.
                        6. Final simplification78.0%

                          \[\leadsto \begin{array}{l} \mathbf{if}\;{\left(\frac{D \cdot M}{d \cdot 2}\right)}^{2} \cdot \frac{h}{\ell} \leq -4 \cdot 10^{+37}:\\ \;\;\;\;\mathsf{fma}\left(-0.125 \cdot \left(D \cdot D\right), \left(\left(\frac{w0}{\left(d \cdot d\right) \cdot \ell} \cdot h\right) \cdot M\right) \cdot M, w0\right)\\ \mathbf{else}:\\ \;\;\;\;1 \cdot w0\\ \end{array} \]
                        7. Add Preprocessing

                        Alternative 6: 89.3% accurate, 1.2× speedup?

                        \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} t_0 := \frac{D\_m \cdot M\_m}{d\_m \cdot 2}\\ t_1 := \frac{M\_m}{d\_m} \cdot M\_m\\ \mathbf{if}\;t\_0 \leq 2 \cdot 10^{-87}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D\_m}{d\_m} \cdot t\_1\right) \cdot D\_m}{\ell}, 1\right)} \cdot w0\\ \mathbf{elif}\;t\_0 \leq 2 \cdot 10^{+211}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{-0.5 \cdot \left(D\_m \cdot M\_m\right)}{d\_m} \cdot \frac{h}{\ell}, \left(\frac{0.5}{d\_m} \cdot M\_m\right) \cdot D\_m, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{\left(\left(-0.25 \cdot h\right) \cdot D\_m\right) \cdot t\_1}{d\_m}, \frac{D\_m}{\ell}, 1\right)} \cdot w0\\ \end{array} \end{array} \]
                        d_m = (fabs.f64 d)
                        D_m = (fabs.f64 D)
                        M_m = (fabs.f64 M)
                        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                        (FPCore (w0 M_m D_m h l d_m)
                         :precision binary64
                         (let* ((t_0 (/ (* D_m M_m) (* d_m 2.0))) (t_1 (* (/ M_m d_m) M_m)))
                           (if (<= t_0 2e-87)
                             (* (sqrt (fma (* -0.25 h) (/ (* (* (/ D_m d_m) t_1) D_m) l) 1.0)) w0)
                             (if (<= t_0 2e+211)
                               (*
                                (sqrt
                                 (fma
                                  (* (/ (* -0.5 (* D_m M_m)) d_m) (/ h l))
                                  (* (* (/ 0.5 d_m) M_m) D_m)
                                  1.0))
                                w0)
                               (*
                                (sqrt (fma (/ (* (* (* -0.25 h) D_m) t_1) d_m) (/ D_m l) 1.0))
                                w0)))))
                        d_m = fabs(d);
                        D_m = fabs(D);
                        M_m = fabs(M);
                        assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                        double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                        	double t_0 = (D_m * M_m) / (d_m * 2.0);
                        	double t_1 = (M_m / d_m) * M_m;
                        	double tmp;
                        	if (t_0 <= 2e-87) {
                        		tmp = sqrt(fma((-0.25 * h), ((((D_m / d_m) * t_1) * D_m) / l), 1.0)) * w0;
                        	} else if (t_0 <= 2e+211) {
                        		tmp = sqrt(fma((((-0.5 * (D_m * M_m)) / d_m) * (h / l)), (((0.5 / d_m) * M_m) * D_m), 1.0)) * w0;
                        	} else {
                        		tmp = sqrt(fma(((((-0.25 * h) * D_m) * t_1) / d_m), (D_m / l), 1.0)) * w0;
                        	}
                        	return tmp;
                        }
                        
                        d_m = abs(d)
                        D_m = abs(D)
                        M_m = abs(M)
                        w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                        function code(w0, M_m, D_m, h, l, d_m)
                        	t_0 = Float64(Float64(D_m * M_m) / Float64(d_m * 2.0))
                        	t_1 = Float64(Float64(M_m / d_m) * M_m)
                        	tmp = 0.0
                        	if (t_0 <= 2e-87)
                        		tmp = Float64(sqrt(fma(Float64(-0.25 * h), Float64(Float64(Float64(Float64(D_m / d_m) * t_1) * D_m) / l), 1.0)) * w0);
                        	elseif (t_0 <= 2e+211)
                        		tmp = Float64(sqrt(fma(Float64(Float64(Float64(-0.5 * Float64(D_m * M_m)) / d_m) * Float64(h / l)), Float64(Float64(Float64(0.5 / d_m) * M_m) * D_m), 1.0)) * w0);
                        	else
                        		tmp = Float64(sqrt(fma(Float64(Float64(Float64(Float64(-0.25 * h) * D_m) * t_1) / d_m), Float64(D_m / l), 1.0)) * w0);
                        	end
                        	return tmp
                        end
                        
                        d_m = N[Abs[d], $MachinePrecision]
                        D_m = N[Abs[D], $MachinePrecision]
                        M_m = N[Abs[M], $MachinePrecision]
                        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                        code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := Block[{t$95$0 = N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(M$95$m / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision]}, If[LessEqual[t$95$0, 2e-87], N[(N[Sqrt[N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(N[(D$95$m / d$95$m), $MachinePrecision] * t$95$1), $MachinePrecision] * D$95$m), $MachinePrecision] / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision], If[LessEqual[t$95$0, 2e+211], N[(N[Sqrt[N[(N[(N[(N[(-0.5 * N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision] / d$95$m), $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(0.5 / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] * D$95$m), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision], N[(N[Sqrt[N[(N[(N[(N[(N[(-0.25 * h), $MachinePrecision] * D$95$m), $MachinePrecision] * t$95$1), $MachinePrecision] / d$95$m), $MachinePrecision] * N[(D$95$m / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision]]]]]
                        
                        \begin{array}{l}
                        d_m = \left|d\right|
                        \\
                        D_m = \left|D\right|
                        \\
                        M_m = \left|M\right|
                        \\
                        [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                        \\
                        \begin{array}{l}
                        t_0 := \frac{D\_m \cdot M\_m}{d\_m \cdot 2}\\
                        t_1 := \frac{M\_m}{d\_m} \cdot M\_m\\
                        \mathbf{if}\;t\_0 \leq 2 \cdot 10^{-87}:\\
                        \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D\_m}{d\_m} \cdot t\_1\right) \cdot D\_m}{\ell}, 1\right)} \cdot w0\\
                        
                        \mathbf{elif}\;t\_0 \leq 2 \cdot 10^{+211}:\\
                        \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{-0.5 \cdot \left(D\_m \cdot M\_m\right)}{d\_m} \cdot \frac{h}{\ell}, \left(\frac{0.5}{d\_m} \cdot M\_m\right) \cdot D\_m, 1\right)} \cdot w0\\
                        
                        \mathbf{else}:\\
                        \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{\left(\left(-0.25 \cdot h\right) \cdot D\_m\right) \cdot t\_1}{d\_m}, \frac{D\_m}{\ell}, 1\right)} \cdot w0\\
                        
                        
                        \end{array}
                        \end{array}
                        
                        Derivation
                        1. Split input into 3 regimes
                        2. if (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) < 2.00000000000000004e-87

                          1. Initial program 78.9%

                            \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                          2. Add Preprocessing
                          3. Taylor expanded in M around 0

                            \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
                          4. Step-by-step derivation
                            1. +-commutativeN/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
                            2. associate-*r/N/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
                            3. associate-*r*N/A

                              \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
                            4. associate-*r*N/A

                              \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
                            5. associate-*l/N/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
                            6. associate-*r/N/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
                            7. lft-mult-inverseN/A

                              \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
                            8. distribute-rgt-inN/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                            9. distribute-lft-inN/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
                            10. associate-*r*N/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
                            11. rgt-mult-inverseN/A

                              \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
                            12. lower-fma.f64N/A

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
                          5. Applied rewrites66.7%

                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
                          6. Step-by-step derivation
                            1. Applied rewrites85.4%

                              \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right) \cdot D}{\color{blue}{\ell}}, 1\right)} \]

                            if 2.00000000000000004e-87 < (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) < 1.9999999999999999e211

                            1. Initial program 79.8%

                              \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                            2. Add Preprocessing
                            3. Step-by-step derivation
                              1. lift--.f64N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}} \]
                              2. sub-negN/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right)}} \]
                              3. +-commutativeN/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) + 1}} \]
                              4. lift-*.f64N/A

                                \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}\right)\right) + 1} \]
                              5. *-commutativeN/A

                                \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{\frac{h}{\ell} \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}}\right)\right) + 1} \]
                              6. distribute-rgt-neg-inN/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{h}{\ell} \cdot \left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right)\right)} + 1} \]
                              7. lift-pow.f64N/A

                                \[\leadsto w0 \cdot \sqrt{\frac{h}{\ell} \cdot \left(\mathsf{neg}\left(\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}}\right)\right) + 1} \]
                              8. unpow2N/A

                                \[\leadsto w0 \cdot \sqrt{\frac{h}{\ell} \cdot \left(\mathsf{neg}\left(\color{blue}{\frac{M \cdot D}{2 \cdot d} \cdot \frac{M \cdot D}{2 \cdot d}}\right)\right) + 1} \]
                              9. distribute-lft-neg-inN/A

                                \[\leadsto w0 \cdot \sqrt{\frac{h}{\ell} \cdot \color{blue}{\left(\left(\mathsf{neg}\left(\frac{M \cdot D}{2 \cdot d}\right)\right) \cdot \frac{M \cdot D}{2 \cdot d}\right)} + 1} \]
                              10. associate-*r*N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{h}{\ell} \cdot \left(\mathsf{neg}\left(\frac{M \cdot D}{2 \cdot d}\right)\right)\right) \cdot \frac{M \cdot D}{2 \cdot d}} + 1} \]
                              11. lower-fma.f64N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\frac{h}{\ell} \cdot \left(\mathsf{neg}\left(\frac{M \cdot D}{2 \cdot d}\right)\right), \frac{M \cdot D}{2 \cdot d}, 1\right)}} \]
                            4. Applied rewrites79.8%

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\frac{h}{\ell} \cdot \frac{\left(D \cdot M\right) \cdot -0.5}{d}, \left(\frac{0.5}{d} \cdot M\right) \cdot D, 1\right)}} \]

                            if 1.9999999999999999e211 < (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d))

                            1. Initial program 55.8%

                              \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                            2. Add Preprocessing
                            3. Taylor expanded in M around 0

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
                            4. Step-by-step derivation
                              1. +-commutativeN/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
                              2. associate-*r/N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
                              3. associate-*r*N/A

                                \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
                              4. associate-*r*N/A

                                \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
                              5. associate-*l/N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
                              6. associate-*r/N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
                              7. lft-mult-inverseN/A

                                \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
                              8. distribute-rgt-inN/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                              9. distribute-lft-inN/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
                              10. associate-*r*N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
                              11. rgt-mult-inverseN/A

                                \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
                              12. lower-fma.f64N/A

                                \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
                            5. Applied rewrites59.1%

                              \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
                            6. Step-by-step derivation
                              1. Applied rewrites66.7%

                                \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(-0.25 \cdot h\right) \cdot \left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right), \color{blue}{\frac{D}{\ell}}, 1\right)} \]
                              2. Step-by-step derivation
                                1. Applied rewrites62.7%

                                  \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\frac{\left(\left(\left(-0.25 \cdot h\right) \cdot M\right) \cdot \frac{M}{d}\right) \cdot D}{d}, \frac{\color{blue}{D}}{\ell}, 1\right)} \]
                                2. Step-by-step derivation
                                  1. Applied rewrites66.6%

                                    \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\frac{\left(D \cdot \left(-0.25 \cdot h\right)\right) \cdot \left(\frac{M}{d} \cdot M\right)}{d}, \frac{D}{\ell}, 1\right)} \]
                                3. Recombined 3 regimes into one program.
                                4. Final simplification82.5%

                                  \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{D \cdot M}{d \cdot 2} \leq 2 \cdot 10^{-87}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D}{d} \cdot \left(\frac{M}{d} \cdot M\right)\right) \cdot D}{\ell}, 1\right)} \cdot w0\\ \mathbf{elif}\;\frac{D \cdot M}{d \cdot 2} \leq 2 \cdot 10^{+211}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{-0.5 \cdot \left(D \cdot M\right)}{d} \cdot \frac{h}{\ell}, \left(\frac{0.5}{d} \cdot M\right) \cdot D, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{\left(\left(-0.25 \cdot h\right) \cdot D\right) \cdot \left(\frac{M}{d} \cdot M\right)}{d}, \frac{D}{\ell}, 1\right)} \cdot w0\\ \end{array} \]
                                5. Add Preprocessing

                                Alternative 7: 86.7% accurate, 1.5× speedup?

                                \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} t_0 := \frac{M\_m}{d\_m} \cdot M\_m\\ \mathbf{if}\;\frac{D\_m \cdot M\_m}{d\_m \cdot 2} \leq 2 \cdot 10^{+63}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D\_m}{d\_m} \cdot t\_0\right) \cdot D\_m}{\ell}, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{\left(\left(-0.25 \cdot h\right) \cdot D\_m\right) \cdot t\_0}{d\_m}, \frac{D\_m}{\ell}, 1\right)} \cdot w0\\ \end{array} \end{array} \]
                                d_m = (fabs.f64 d)
                                D_m = (fabs.f64 D)
                                M_m = (fabs.f64 M)
                                NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                (FPCore (w0 M_m D_m h l d_m)
                                 :precision binary64
                                 (let* ((t_0 (* (/ M_m d_m) M_m)))
                                   (if (<= (/ (* D_m M_m) (* d_m 2.0)) 2e+63)
                                     (* (sqrt (fma (* -0.25 h) (/ (* (* (/ D_m d_m) t_0) D_m) l) 1.0)) w0)
                                     (* (sqrt (fma (/ (* (* (* -0.25 h) D_m) t_0) d_m) (/ D_m l) 1.0)) w0))))
                                d_m = fabs(d);
                                D_m = fabs(D);
                                M_m = fabs(M);
                                assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                                double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                	double t_0 = (M_m / d_m) * M_m;
                                	double tmp;
                                	if (((D_m * M_m) / (d_m * 2.0)) <= 2e+63) {
                                		tmp = sqrt(fma((-0.25 * h), ((((D_m / d_m) * t_0) * D_m) / l), 1.0)) * w0;
                                	} else {
                                		tmp = sqrt(fma(((((-0.25 * h) * D_m) * t_0) / d_m), (D_m / l), 1.0)) * w0;
                                	}
                                	return tmp;
                                }
                                
                                d_m = abs(d)
                                D_m = abs(D)
                                M_m = abs(M)
                                w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                                function code(w0, M_m, D_m, h, l, d_m)
                                	t_0 = Float64(Float64(M_m / d_m) * M_m)
                                	tmp = 0.0
                                	if (Float64(Float64(D_m * M_m) / Float64(d_m * 2.0)) <= 2e+63)
                                		tmp = Float64(sqrt(fma(Float64(-0.25 * h), Float64(Float64(Float64(Float64(D_m / d_m) * t_0) * D_m) / l), 1.0)) * w0);
                                	else
                                		tmp = Float64(sqrt(fma(Float64(Float64(Float64(Float64(-0.25 * h) * D_m) * t_0) / d_m), Float64(D_m / l), 1.0)) * w0);
                                	end
                                	return tmp
                                end
                                
                                d_m = N[Abs[d], $MachinePrecision]
                                D_m = N[Abs[D], $MachinePrecision]
                                M_m = N[Abs[M], $MachinePrecision]
                                NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := Block[{t$95$0 = N[(N[(M$95$m / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision]}, If[LessEqual[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d$95$m * 2.0), $MachinePrecision]), $MachinePrecision], 2e+63], N[(N[Sqrt[N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(N[(D$95$m / d$95$m), $MachinePrecision] * t$95$0), $MachinePrecision] * D$95$m), $MachinePrecision] / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision], N[(N[Sqrt[N[(N[(N[(N[(N[(-0.25 * h), $MachinePrecision] * D$95$m), $MachinePrecision] * t$95$0), $MachinePrecision] / d$95$m), $MachinePrecision] * N[(D$95$m / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision]]]
                                
                                \begin{array}{l}
                                d_m = \left|d\right|
                                \\
                                D_m = \left|D\right|
                                \\
                                M_m = \left|M\right|
                                \\
                                [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                                \\
                                \begin{array}{l}
                                t_0 := \frac{M\_m}{d\_m} \cdot M\_m\\
                                \mathbf{if}\;\frac{D\_m \cdot M\_m}{d\_m \cdot 2} \leq 2 \cdot 10^{+63}:\\
                                \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D\_m}{d\_m} \cdot t\_0\right) \cdot D\_m}{\ell}, 1\right)} \cdot w0\\
                                
                                \mathbf{else}:\\
                                \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{\left(\left(-0.25 \cdot h\right) \cdot D\_m\right) \cdot t\_0}{d\_m}, \frac{D\_m}{\ell}, 1\right)} \cdot w0\\
                                
                                
                                \end{array}
                                \end{array}
                                
                                Derivation
                                1. Split input into 2 regimes
                                2. if (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) < 2.00000000000000012e63

                                  1. Initial program 81.5%

                                    \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                  2. Add Preprocessing
                                  3. Taylor expanded in M around 0

                                    \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
                                  4. Step-by-step derivation
                                    1. +-commutativeN/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
                                    2. associate-*r/N/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
                                    3. associate-*r*N/A

                                      \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
                                    4. associate-*r*N/A

                                      \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
                                    5. associate-*l/N/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
                                    6. associate-*r/N/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
                                    7. lft-mult-inverseN/A

                                      \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
                                    8. distribute-rgt-inN/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                    9. distribute-lft-inN/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
                                    10. associate-*r*N/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
                                    11. rgt-mult-inverseN/A

                                      \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
                                    12. lower-fma.f64N/A

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
                                  5. Applied rewrites64.8%

                                    \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
                                  6. Step-by-step derivation
                                    1. Applied rewrites85.4%

                                      \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right) \cdot D}{\color{blue}{\ell}}, 1\right)} \]

                                    if 2.00000000000000012e63 < (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d))

                                    1. Initial program 49.8%

                                      \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                    2. Add Preprocessing
                                    3. Taylor expanded in M around 0

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
                                    4. Step-by-step derivation
                                      1. +-commutativeN/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
                                      2. associate-*r/N/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
                                      3. associate-*r*N/A

                                        \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
                                      4. associate-*r*N/A

                                        \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
                                      5. associate-*l/N/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
                                      6. associate-*r/N/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
                                      7. lft-mult-inverseN/A

                                        \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
                                      8. distribute-rgt-inN/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                      9. distribute-lft-inN/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
                                      10. associate-*r*N/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
                                      11. rgt-mult-inverseN/A

                                        \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
                                      12. lower-fma.f64N/A

                                        \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
                                    5. Applied rewrites44.7%

                                      \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
                                    6. Step-by-step derivation
                                      1. Applied rewrites57.5%

                                        \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\left(-0.25 \cdot h\right) \cdot \left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right), \color{blue}{\frac{D}{\ell}}, 1\right)} \]
                                      2. Step-by-step derivation
                                        1. Applied rewrites54.8%

                                          \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\frac{\left(\left(\left(-0.25 \cdot h\right) \cdot M\right) \cdot \frac{M}{d}\right) \cdot D}{d}, \frac{\color{blue}{D}}{\ell}, 1\right)} \]
                                        2. Step-by-step derivation
                                          1. Applied rewrites57.4%

                                            \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(\frac{\left(D \cdot \left(-0.25 \cdot h\right)\right) \cdot \left(\frac{M}{d} \cdot M\right)}{d}, \frac{D}{\ell}, 1\right)} \]
                                        3. Recombined 2 regimes into one program.
                                        4. Final simplification81.0%

                                          \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{D \cdot M}{d \cdot 2} \leq 2 \cdot 10^{+63}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D}{d} \cdot \left(\frac{M}{d} \cdot M\right)\right) \cdot D}{\ell}, 1\right)} \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\mathsf{fma}\left(\frac{\left(\left(-0.25 \cdot h\right) \cdot D\right) \cdot \left(\frac{M}{d} \cdot M\right)}{d}, \frac{D}{\ell}, 1\right)} \cdot w0\\ \end{array} \]
                                        5. Add Preprocessing

                                        Alternative 8: 89.2% accurate, 1.9× speedup?

                                        \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d\_m} \cdot M\_m\right) \cdot D\_m, \frac{\left(\frac{M\_m}{d\_m} \cdot h\right) \cdot \left(0.5 \cdot D\_m\right)}{-\ell}, 1\right)} \cdot w0 \end{array} \]
                                        d_m = (fabs.f64 d)
                                        D_m = (fabs.f64 D)
                                        M_m = (fabs.f64 M)
                                        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                        (FPCore (w0 M_m D_m h l d_m)
                                         :precision binary64
                                         (*
                                          (sqrt
                                           (fma
                                            (* (* (/ 0.5 d_m) M_m) D_m)
                                            (/ (* (* (/ M_m d_m) h) (* 0.5 D_m)) (- l))
                                            1.0))
                                          w0))
                                        d_m = fabs(d);
                                        D_m = fabs(D);
                                        M_m = fabs(M);
                                        assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                                        double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                        	return sqrt(fma((((0.5 / d_m) * M_m) * D_m), ((((M_m / d_m) * h) * (0.5 * D_m)) / -l), 1.0)) * w0;
                                        }
                                        
                                        d_m = abs(d)
                                        D_m = abs(D)
                                        M_m = abs(M)
                                        w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                                        function code(w0, M_m, D_m, h, l, d_m)
                                        	return Float64(sqrt(fma(Float64(Float64(Float64(0.5 / d_m) * M_m) * D_m), Float64(Float64(Float64(Float64(M_m / d_m) * h) * Float64(0.5 * D_m)) / Float64(-l)), 1.0)) * w0)
                                        end
                                        
                                        d_m = N[Abs[d], $MachinePrecision]
                                        D_m = N[Abs[D], $MachinePrecision]
                                        M_m = N[Abs[M], $MachinePrecision]
                                        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                        code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := N[(N[Sqrt[N[(N[(N[(N[(0.5 / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] * D$95$m), $MachinePrecision] * N[(N[(N[(N[(M$95$m / d$95$m), $MachinePrecision] * h), $MachinePrecision] * N[(0.5 * D$95$m), $MachinePrecision]), $MachinePrecision] / (-l)), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision]
                                        
                                        \begin{array}{l}
                                        d_m = \left|d\right|
                                        \\
                                        D_m = \left|D\right|
                                        \\
                                        M_m = \left|M\right|
                                        \\
                                        [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                                        \\
                                        \sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d\_m} \cdot M\_m\right) \cdot D\_m, \frac{\left(\frac{M\_m}{d\_m} \cdot h\right) \cdot \left(0.5 \cdot D\_m\right)}{-\ell}, 1\right)} \cdot w0
                                        \end{array}
                                        
                                        Derivation
                                        1. Initial program 76.5%

                                          \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                        2. Add Preprocessing
                                        3. Step-by-step derivation
                                          1. lift--.f64N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}} \]
                                          2. sub-negN/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right)}} \]
                                          3. +-commutativeN/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) + 1}} \]
                                          4. lift-*.f64N/A

                                            \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}\right)\right) + 1} \]
                                          5. lift-/.f64N/A

                                            \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left({\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \color{blue}{\frac{h}{\ell}}\right)\right) + 1} \]
                                          6. associate-*r/N/A

                                            \[\leadsto w0 \cdot \sqrt{\left(\mathsf{neg}\left(\color{blue}{\frac{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot h}{\ell}}\right)\right) + 1} \]
                                          7. distribute-neg-frac2N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot h}{\mathsf{neg}\left(\ell\right)}} + 1} \]
                                          8. lift-pow.f64N/A

                                            \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{{\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}} \cdot h}{\mathsf{neg}\left(\ell\right)} + 1} \]
                                          9. unpow2N/A

                                            \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{M \cdot D}{2 \cdot d} \cdot \frac{M \cdot D}{2 \cdot d}\right)} \cdot h}{\mathsf{neg}\left(\ell\right)} + 1} \]
                                          10. associate-*l*N/A

                                            \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\frac{M \cdot D}{2 \cdot d} \cdot \left(\frac{M \cdot D}{2 \cdot d} \cdot h\right)}}{\mathsf{neg}\left(\ell\right)} + 1} \]
                                          11. associate-/l*N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{M \cdot D}{2 \cdot d} \cdot \frac{\frac{M \cdot D}{2 \cdot d} \cdot h}{\mathsf{neg}\left(\ell\right)}} + 1} \]
                                          12. lower-fma.f64N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\frac{M \cdot D}{2 \cdot d}, \frac{\frac{M \cdot D}{2 \cdot d} \cdot h}{\mathsf{neg}\left(\ell\right)}, 1\right)}} \]
                                        4. Applied rewrites83.7%

                                          \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(\left(\frac{0.5}{d} \cdot M\right) \cdot D, \frac{\left(D \cdot 0.5\right) \cdot \left(\frac{M}{d} \cdot h\right)}{-\ell}, 1\right)}} \]
                                        5. Final simplification83.7%

                                          \[\leadsto \sqrt{\mathsf{fma}\left(\left(\frac{0.5}{d} \cdot M\right) \cdot D, \frac{\left(\frac{M}{d} \cdot h\right) \cdot \left(0.5 \cdot D\right)}{-\ell}, 1\right)} \cdot w0 \]
                                        6. Add Preprocessing

                                        Alternative 9: 83.9% accurate, 2.1× speedup?

                                        \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D\_m}{d\_m} \cdot \left(\frac{M\_m}{d\_m} \cdot M\_m\right)\right) \cdot D\_m}{\ell}, 1\right)} \cdot w0 \end{array} \]
                                        d_m = (fabs.f64 d)
                                        D_m = (fabs.f64 D)
                                        M_m = (fabs.f64 M)
                                        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                        (FPCore (w0 M_m D_m h l d_m)
                                         :precision binary64
                                         (*
                                          (sqrt
                                           (fma (* -0.25 h) (/ (* (* (/ D_m d_m) (* (/ M_m d_m) M_m)) D_m) l) 1.0))
                                          w0))
                                        d_m = fabs(d);
                                        D_m = fabs(D);
                                        M_m = fabs(M);
                                        assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                                        double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                        	return sqrt(fma((-0.25 * h), ((((D_m / d_m) * ((M_m / d_m) * M_m)) * D_m) / l), 1.0)) * w0;
                                        }
                                        
                                        d_m = abs(d)
                                        D_m = abs(D)
                                        M_m = abs(M)
                                        w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                                        function code(w0, M_m, D_m, h, l, d_m)
                                        	return Float64(sqrt(fma(Float64(-0.25 * h), Float64(Float64(Float64(Float64(D_m / d_m) * Float64(Float64(M_m / d_m) * M_m)) * D_m) / l), 1.0)) * w0)
                                        end
                                        
                                        d_m = N[Abs[d], $MachinePrecision]
                                        D_m = N[Abs[D], $MachinePrecision]
                                        M_m = N[Abs[M], $MachinePrecision]
                                        NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                        code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := N[(N[Sqrt[N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(N[(D$95$m / d$95$m), $MachinePrecision] * N[(N[(M$95$m / d$95$m), $MachinePrecision] * M$95$m), $MachinePrecision]), $MachinePrecision] * D$95$m), $MachinePrecision] / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision] * w0), $MachinePrecision]
                                        
                                        \begin{array}{l}
                                        d_m = \left|d\right|
                                        \\
                                        D_m = \left|D\right|
                                        \\
                                        M_m = \left|M\right|
                                        \\
                                        [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                                        \\
                                        \sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D\_m}{d\_m} \cdot \left(\frac{M\_m}{d\_m} \cdot M\_m\right)\right) \cdot D\_m}{\ell}, 1\right)} \cdot w0
                                        \end{array}
                                        
                                        Derivation
                                        1. Initial program 76.5%

                                          \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                        2. Add Preprocessing
                                        3. Taylor expanded in M around 0

                                          \[\leadsto w0 \cdot \sqrt{\color{blue}{1 + \frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell}}} \]
                                        4. Step-by-step derivation
                                          1. +-commutativeN/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{-1}{4} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot h\right)}{{d}^{2} \cdot \ell} + 1}} \]
                                          2. associate-*r/N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot \left({M}^{2} \cdot h\right)\right)}{{d}^{2} \cdot \ell}} + 1} \]
                                          3. associate-*r*N/A

                                            \[\leadsto w0 \cdot \sqrt{\frac{\frac{-1}{4} \cdot \color{blue}{\left(\left({D}^{2} \cdot {M}^{2}\right) \cdot h\right)}}{{d}^{2} \cdot \ell} + 1} \]
                                          4. associate-*r*N/A

                                            \[\leadsto w0 \cdot \sqrt{\frac{\color{blue}{\left(\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)\right) \cdot h}}{{d}^{2} \cdot \ell} + 1} \]
                                          5. associate-*l/N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\frac{\frac{-1}{4} \cdot \left({D}^{2} \cdot {M}^{2}\right)}{{d}^{2} \cdot \ell} \cdot h} + 1} \]
                                          6. associate-*r/N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right)} \cdot h + 1} \]
                                          7. lft-mult-inverseN/A

                                            \[\leadsto w0 \cdot \sqrt{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) \cdot h + \color{blue}{\frac{1}{h} \cdot h}} \]
                                          8. distribute-rgt-inN/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                          9. distribute-lft-inN/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{h \cdot \left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}\right) + h \cdot \frac{1}{h}}} \]
                                          10. associate-*r*N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}} + h \cdot \frac{1}{h}} \]
                                          11. rgt-mult-inverseN/A

                                            \[\leadsto w0 \cdot \sqrt{\left(h \cdot \frac{-1}{4}\right) \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \color{blue}{1}} \]
                                          12. lower-fma.f64N/A

                                            \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot \frac{-1}{4}, \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell}, 1\right)}} \]
                                        5. Applied rewrites61.7%

                                          \[\leadsto w0 \cdot \sqrt{\color{blue}{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(M \cdot M\right) \cdot D\right) \cdot D}{\left(d \cdot d\right) \cdot \ell}, 1\right)}} \]
                                        6. Step-by-step derivation
                                          1. Applied rewrites79.9%

                                            \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(\frac{M}{d} \cdot M\right) \cdot \frac{D}{d}\right) \cdot D}{\color{blue}{\ell}}, 1\right)} \]
                                          2. Final simplification79.9%

                                            \[\leadsto \sqrt{\mathsf{fma}\left(-0.25 \cdot h, \frac{\left(\frac{D}{d} \cdot \left(\frac{M}{d} \cdot M\right)\right) \cdot D}{\ell}, 1\right)} \cdot w0 \]
                                          3. Add Preprocessing

                                          Alternative 10: 78.7% accurate, 2.5× speedup?

                                          \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;D\_m \cdot M\_m \leq 4 \cdot 10^{-109}:\\ \;\;\;\;1 \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\left(\left(h \cdot M\_m\right) \cdot \left(\frac{w0}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot M\_m\right)\right) \cdot D\_m, -0.125 \cdot D\_m, w0\right)\\ \end{array} \end{array} \]
                                          d_m = (fabs.f64 d)
                                          D_m = (fabs.f64 D)
                                          M_m = (fabs.f64 M)
                                          NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                          (FPCore (w0 M_m D_m h l d_m)
                                           :precision binary64
                                           (if (<= (* D_m M_m) 4e-109)
                                             (* 1.0 w0)
                                             (fma
                                              (* (* (* h M_m) (* (/ w0 (* (* d_m d_m) l)) M_m)) D_m)
                                              (* -0.125 D_m)
                                              w0)))
                                          d_m = fabs(d);
                                          D_m = fabs(D);
                                          M_m = fabs(M);
                                          assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                                          double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                          	double tmp;
                                          	if ((D_m * M_m) <= 4e-109) {
                                          		tmp = 1.0 * w0;
                                          	} else {
                                          		tmp = fma((((h * M_m) * ((w0 / ((d_m * d_m) * l)) * M_m)) * D_m), (-0.125 * D_m), w0);
                                          	}
                                          	return tmp;
                                          }
                                          
                                          d_m = abs(d)
                                          D_m = abs(D)
                                          M_m = abs(M)
                                          w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                                          function code(w0, M_m, D_m, h, l, d_m)
                                          	tmp = 0.0
                                          	if (Float64(D_m * M_m) <= 4e-109)
                                          		tmp = Float64(1.0 * w0);
                                          	else
                                          		tmp = fma(Float64(Float64(Float64(h * M_m) * Float64(Float64(w0 / Float64(Float64(d_m * d_m) * l)) * M_m)) * D_m), Float64(-0.125 * D_m), w0);
                                          	end
                                          	return tmp
                                          end
                                          
                                          d_m = N[Abs[d], $MachinePrecision]
                                          D_m = N[Abs[D], $MachinePrecision]
                                          M_m = N[Abs[M], $MachinePrecision]
                                          NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                          code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[(D$95$m * M$95$m), $MachinePrecision], 4e-109], N[(1.0 * w0), $MachinePrecision], N[(N[(N[(N[(h * M$95$m), $MachinePrecision] * N[(N[(w0 / N[(N[(d$95$m * d$95$m), $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision] * M$95$m), $MachinePrecision]), $MachinePrecision] * D$95$m), $MachinePrecision] * N[(-0.125 * D$95$m), $MachinePrecision] + w0), $MachinePrecision]]
                                          
                                          \begin{array}{l}
                                          d_m = \left|d\right|
                                          \\
                                          D_m = \left|D\right|
                                          \\
                                          M_m = \left|M\right|
                                          \\
                                          [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                                          \\
                                          \begin{array}{l}
                                          \mathbf{if}\;D\_m \cdot M\_m \leq 4 \cdot 10^{-109}:\\
                                          \;\;\;\;1 \cdot w0\\
                                          
                                          \mathbf{else}:\\
                                          \;\;\;\;\mathsf{fma}\left(\left(\left(h \cdot M\_m\right) \cdot \left(\frac{w0}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot M\_m\right)\right) \cdot D\_m, -0.125 \cdot D\_m, w0\right)\\
                                          
                                          
                                          \end{array}
                                          \end{array}
                                          
                                          Derivation
                                          1. Split input into 2 regimes
                                          2. if (*.f64 M D) < 4e-109

                                            1. Initial program 79.9%

                                              \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                            2. Add Preprocessing
                                            3. Taylor expanded in M around 0

                                              \[\leadsto w0 \cdot \color{blue}{1} \]
                                            4. Step-by-step derivation
                                              1. Applied rewrites77.8%

                                                \[\leadsto w0 \cdot \color{blue}{1} \]

                                              if 4e-109 < (*.f64 M D)

                                              1. Initial program 71.3%

                                                \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                              2. Add Preprocessing
                                              3. Taylor expanded in M around 0

                                                \[\leadsto \color{blue}{w0 + \frac{-1}{8} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell}} \]
                                              4. Step-by-step derivation
                                                1. +-commutativeN/A

                                                  \[\leadsto \color{blue}{\frac{-1}{8} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell} + w0} \]
                                                2. *-commutativeN/A

                                                  \[\leadsto \color{blue}{\frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell} \cdot \frac{-1}{8}} + w0 \]
                                                3. associate-/l*N/A

                                                  \[\leadsto \color{blue}{\left({D}^{2} \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}\right)} \cdot \frac{-1}{8} + w0 \]
                                                4. associate-*r*N/A

                                                  \[\leadsto \color{blue}{{D}^{2} \cdot \left(\frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell} \cdot \frac{-1}{8}\right)} + w0 \]
                                                5. *-commutativeN/A

                                                  \[\leadsto {D}^{2} \cdot \color{blue}{\left(\frac{-1}{8} \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}\right)} + w0 \]
                                                6. associate-*r*N/A

                                                  \[\leadsto \color{blue}{\left({D}^{2} \cdot \frac{-1}{8}\right) \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}} + w0 \]
                                                7. lower-fma.f64N/A

                                                  \[\leadsto \color{blue}{\mathsf{fma}\left({D}^{2} \cdot \frac{-1}{8}, \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}, w0\right)} \]
                                              5. Applied rewrites38.2%

                                                \[\leadsto \color{blue}{\mathsf{fma}\left(\left(D \cdot D\right) \cdot -0.125, \frac{\left(M \cdot M\right) \cdot h}{\ell} \cdot \frac{w0}{d \cdot d}, w0\right)} \]
                                              6. Step-by-step derivation
                                                1. Applied rewrites38.4%

                                                  \[\leadsto \mathsf{fma}\left(\left(D \cdot D\right) \cdot -0.125, \left(\left(M \cdot M\right) \cdot h\right) \cdot \color{blue}{\frac{w0}{\left(d \cdot d\right) \cdot \ell}}, w0\right) \]
                                                2. Step-by-step derivation
                                                  1. Applied rewrites65.1%

                                                    \[\leadsto \mathsf{fma}\left(\left(\left(\frac{w0}{\left(d \cdot d\right) \cdot \ell} \cdot M\right) \cdot \left(h \cdot M\right)\right) \cdot D, \color{blue}{-0.125 \cdot D}, w0\right) \]
                                                3. Recombined 2 regimes into one program.
                                                4. Final simplification72.9%

                                                  \[\leadsto \begin{array}{l} \mathbf{if}\;D \cdot M \leq 4 \cdot 10^{-109}:\\ \;\;\;\;1 \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(\left(\left(h \cdot M\right) \cdot \left(\frac{w0}{\left(d \cdot d\right) \cdot \ell} \cdot M\right)\right) \cdot D, -0.125 \cdot D, w0\right)\\ \end{array} \]
                                                5. Add Preprocessing

                                                Alternative 11: 77.2% accurate, 2.5× speedup?

                                                \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ \begin{array}{l} \mathbf{if}\;D\_m \cdot M\_m \leq 10^{-72}:\\ \;\;\;\;1 \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(D\_m, \frac{\left(-0.125 \cdot w0\right) \cdot \left(\left(M\_m \cdot M\_m\right) \cdot h\right)}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot D\_m, w0\right)\\ \end{array} \end{array} \]
                                                d_m = (fabs.f64 d)
                                                D_m = (fabs.f64 D)
                                                M_m = (fabs.f64 M)
                                                NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                                (FPCore (w0 M_m D_m h l d_m)
                                                 :precision binary64
                                                 (if (<= (* D_m M_m) 1e-72)
                                                   (* 1.0 w0)
                                                   (fma
                                                    D_m
                                                    (* (/ (* (* -0.125 w0) (* (* M_m M_m) h)) (* (* d_m d_m) l)) D_m)
                                                    w0)))
                                                d_m = fabs(d);
                                                D_m = fabs(D);
                                                M_m = fabs(M);
                                                assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                                                double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                                	double tmp;
                                                	if ((D_m * M_m) <= 1e-72) {
                                                		tmp = 1.0 * w0;
                                                	} else {
                                                		tmp = fma(D_m, ((((-0.125 * w0) * ((M_m * M_m) * h)) / ((d_m * d_m) * l)) * D_m), w0);
                                                	}
                                                	return tmp;
                                                }
                                                
                                                d_m = abs(d)
                                                D_m = abs(D)
                                                M_m = abs(M)
                                                w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                                                function code(w0, M_m, D_m, h, l, d_m)
                                                	tmp = 0.0
                                                	if (Float64(D_m * M_m) <= 1e-72)
                                                		tmp = Float64(1.0 * w0);
                                                	else
                                                		tmp = fma(D_m, Float64(Float64(Float64(Float64(-0.125 * w0) * Float64(Float64(M_m * M_m) * h)) / Float64(Float64(d_m * d_m) * l)) * D_m), w0);
                                                	end
                                                	return tmp
                                                end
                                                
                                                d_m = N[Abs[d], $MachinePrecision]
                                                D_m = N[Abs[D], $MachinePrecision]
                                                M_m = N[Abs[M], $MachinePrecision]
                                                NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                                code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := If[LessEqual[N[(D$95$m * M$95$m), $MachinePrecision], 1e-72], N[(1.0 * w0), $MachinePrecision], N[(D$95$m * N[(N[(N[(N[(-0.125 * w0), $MachinePrecision] * N[(N[(M$95$m * M$95$m), $MachinePrecision] * h), $MachinePrecision]), $MachinePrecision] / N[(N[(d$95$m * d$95$m), $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision] * D$95$m), $MachinePrecision] + w0), $MachinePrecision]]
                                                
                                                \begin{array}{l}
                                                d_m = \left|d\right|
                                                \\
                                                D_m = \left|D\right|
                                                \\
                                                M_m = \left|M\right|
                                                \\
                                                [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                                                \\
                                                \begin{array}{l}
                                                \mathbf{if}\;D\_m \cdot M\_m \leq 10^{-72}:\\
                                                \;\;\;\;1 \cdot w0\\
                                                
                                                \mathbf{else}:\\
                                                \;\;\;\;\mathsf{fma}\left(D\_m, \frac{\left(-0.125 \cdot w0\right) \cdot \left(\left(M\_m \cdot M\_m\right) \cdot h\right)}{\left(d\_m \cdot d\_m\right) \cdot \ell} \cdot D\_m, w0\right)\\
                                                
                                                
                                                \end{array}
                                                \end{array}
                                                
                                                Derivation
                                                1. Split input into 2 regimes
                                                2. if (*.f64 M D) < 9.9999999999999997e-73

                                                  1. Initial program 80.6%

                                                    \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                                  2. Add Preprocessing
                                                  3. Taylor expanded in M around 0

                                                    \[\leadsto w0 \cdot \color{blue}{1} \]
                                                  4. Step-by-step derivation
                                                    1. Applied rewrites78.1%

                                                      \[\leadsto w0 \cdot \color{blue}{1} \]

                                                    if 9.9999999999999997e-73 < (*.f64 M D)

                                                    1. Initial program 69.5%

                                                      \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                                    2. Add Preprocessing
                                                    3. Taylor expanded in M around 0

                                                      \[\leadsto \color{blue}{w0 + \frac{-1}{8} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell}} \]
                                                    4. Step-by-step derivation
                                                      1. +-commutativeN/A

                                                        \[\leadsto \color{blue}{\frac{-1}{8} \cdot \frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell} + w0} \]
                                                      2. *-commutativeN/A

                                                        \[\leadsto \color{blue}{\frac{{D}^{2} \cdot \left({M}^{2} \cdot \left(h \cdot w0\right)\right)}{{d}^{2} \cdot \ell} \cdot \frac{-1}{8}} + w0 \]
                                                      3. associate-/l*N/A

                                                        \[\leadsto \color{blue}{\left({D}^{2} \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}\right)} \cdot \frac{-1}{8} + w0 \]
                                                      4. associate-*r*N/A

                                                        \[\leadsto \color{blue}{{D}^{2} \cdot \left(\frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell} \cdot \frac{-1}{8}\right)} + w0 \]
                                                      5. *-commutativeN/A

                                                        \[\leadsto {D}^{2} \cdot \color{blue}{\left(\frac{-1}{8} \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}\right)} + w0 \]
                                                      6. associate-*r*N/A

                                                        \[\leadsto \color{blue}{\left({D}^{2} \cdot \frac{-1}{8}\right) \cdot \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}} + w0 \]
                                                      7. lower-fma.f64N/A

                                                        \[\leadsto \color{blue}{\mathsf{fma}\left({D}^{2} \cdot \frac{-1}{8}, \frac{{M}^{2} \cdot \left(h \cdot w0\right)}{{d}^{2} \cdot \ell}, w0\right)} \]
                                                    5. Applied rewrites36.4%

                                                      \[\leadsto \color{blue}{\mathsf{fma}\left(\left(D \cdot D\right) \cdot -0.125, \frac{\left(M \cdot M\right) \cdot h}{\ell} \cdot \frac{w0}{d \cdot d}, w0\right)} \]
                                                    6. Step-by-step derivation
                                                      1. Applied rewrites50.8%

                                                        \[\leadsto \mathsf{fma}\left(D, \color{blue}{D \cdot \left(\left(-0.125 \cdot \frac{w0}{d \cdot d}\right) \cdot \frac{\left(M \cdot M\right) \cdot h}{\ell}\right)}, w0\right) \]
                                                      2. Step-by-step derivation
                                                        1. Applied rewrites50.9%

                                                          \[\leadsto \mathsf{fma}\left(D, D \cdot \frac{\left(\left(M \cdot M\right) \cdot h\right) \cdot \left(-0.125 \cdot w0\right)}{\color{blue}{\left(d \cdot d\right) \cdot \ell}}, w0\right) \]
                                                      3. Recombined 2 regimes into one program.
                                                      4. Final simplification68.1%

                                                        \[\leadsto \begin{array}{l} \mathbf{if}\;D \cdot M \leq 10^{-72}:\\ \;\;\;\;1 \cdot w0\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(D, \frac{\left(-0.125 \cdot w0\right) \cdot \left(\left(M \cdot M\right) \cdot h\right)}{\left(d \cdot d\right) \cdot \ell} \cdot D, w0\right)\\ \end{array} \]
                                                      5. Add Preprocessing

                                                      Alternative 12: 68.0% accurate, 26.2× speedup?

                                                      \[\begin{array}{l} d_m = \left|d\right| \\ D_m = \left|D\right| \\ M_m = \left|M\right| \\ [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\ \\ 1 \cdot w0 \end{array} \]
                                                      d_m = (fabs.f64 d)
                                                      D_m = (fabs.f64 D)
                                                      M_m = (fabs.f64 M)
                                                      NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                                      (FPCore (w0 M_m D_m h l d_m) :precision binary64 (* 1.0 w0))
                                                      d_m = fabs(d);
                                                      D_m = fabs(D);
                                                      M_m = fabs(M);
                                                      assert(w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m);
                                                      double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                                      	return 1.0 * w0;
                                                      }
                                                      
                                                      d_m = abs(d)
                                                      D_m = abs(d)
                                                      M_m = abs(m)
                                                      NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                                      real(8) function code(w0, m_m, d_m, h, l, d_m_1)
                                                          real(8), intent (in) :: w0
                                                          real(8), intent (in) :: m_m
                                                          real(8), intent (in) :: d_m
                                                          real(8), intent (in) :: h
                                                          real(8), intent (in) :: l
                                                          real(8), intent (in) :: d_m_1
                                                          code = 1.0d0 * w0
                                                      end function
                                                      
                                                      d_m = Math.abs(d);
                                                      D_m = Math.abs(D);
                                                      M_m = Math.abs(M);
                                                      assert w0 < M_m && M_m < D_m && D_m < h && h < l && l < d_m;
                                                      public static double code(double w0, double M_m, double D_m, double h, double l, double d_m) {
                                                      	return 1.0 * w0;
                                                      }
                                                      
                                                      d_m = math.fabs(d)
                                                      D_m = math.fabs(D)
                                                      M_m = math.fabs(M)
                                                      [w0, M_m, D_m, h, l, d_m] = sort([w0, M_m, D_m, h, l, d_m])
                                                      def code(w0, M_m, D_m, h, l, d_m):
                                                      	return 1.0 * w0
                                                      
                                                      d_m = abs(d)
                                                      D_m = abs(D)
                                                      M_m = abs(M)
                                                      w0, M_m, D_m, h, l, d_m = sort([w0, M_m, D_m, h, l, d_m])
                                                      function code(w0, M_m, D_m, h, l, d_m)
                                                      	return Float64(1.0 * w0)
                                                      end
                                                      
                                                      d_m = abs(d);
                                                      D_m = abs(D);
                                                      M_m = abs(M);
                                                      w0, M_m, D_m, h, l, d_m = num2cell(sort([w0, M_m, D_m, h, l, d_m])){:}
                                                      function tmp = code(w0, M_m, D_m, h, l, d_m)
                                                      	tmp = 1.0 * w0;
                                                      end
                                                      
                                                      d_m = N[Abs[d], $MachinePrecision]
                                                      D_m = N[Abs[D], $MachinePrecision]
                                                      M_m = N[Abs[M], $MachinePrecision]
                                                      NOTE: w0, M_m, D_m, h, l, and d_m should be sorted in increasing order before calling this function.
                                                      code[w0_, M$95$m_, D$95$m_, h_, l_, d$95$m_] := N[(1.0 * w0), $MachinePrecision]
                                                      
                                                      \begin{array}{l}
                                                      d_m = \left|d\right|
                                                      \\
                                                      D_m = \left|D\right|
                                                      \\
                                                      M_m = \left|M\right|
                                                      \\
                                                      [w0, M_m, D_m, h, l, d_m] = \mathsf{sort}([w0, M_m, D_m, h, l, d_m])\\
                                                      \\
                                                      1 \cdot w0
                                                      \end{array}
                                                      
                                                      Derivation
                                                      1. Initial program 76.5%

                                                        \[w0 \cdot \sqrt{1 - {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \]
                                                      2. Add Preprocessing
                                                      3. Taylor expanded in M around 0

                                                        \[\leadsto w0 \cdot \color{blue}{1} \]
                                                      4. Step-by-step derivation
                                                        1. Applied rewrites67.2%

                                                          \[\leadsto w0 \cdot \color{blue}{1} \]
                                                        2. Final simplification67.2%

                                                          \[\leadsto 1 \cdot w0 \]
                                                        3. Add Preprocessing

                                                        Reproduce

                                                        ?
                                                        herbie shell --seed 2024325 
                                                        (FPCore (w0 M D h l d)
                                                          :name "Henrywood and Agarwal, Equation (9a)"
                                                          :precision binary64
                                                          (* w0 (sqrt (- 1.0 (* (pow (/ (* M D) (* 2.0 d)) 2.0) (/ h l))))))