Henrywood and Agarwal, Equation (9a)

Percentage Accurate: 81.3% → 85.8%
Time: 10.4s
Alternatives: 11
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 11 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.3% 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: 85.8% accurate, 1.6× speedup?

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

\mathbf{elif}\;M\_m \cdot D\_m \leq 5 \cdot 10^{+106}:\\
\;\;\;\;w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(\left(D\_m \cdot M\_m\right) \cdot M\_m\right) \cdot D\_m}{\left(\ell \cdot d\right) \cdot d}, 1\right)}\\

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


\end{array}
\end{array}
Derivation
  1. Split input into 3 regimes
  2. if (*.f64 M D) < 5.0000000000000002e-136

    1. Initial program 82.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 h around inf

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

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

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

        \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
      4. 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}}} \]
      5. 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}} \]
      6. 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}} \]
      7. 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 rewrites62.6%

      \[\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 rewrites74.6%

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

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

        if 5.0000000000000002e-136 < (*.f64 M D) < 4.9999999999999998e106

        1. Initial program 78.7%

          \[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 h around inf

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

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

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

            \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
          4. 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}}} \]
          5. 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}} \]
          6. 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}} \]
          7. 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 rewrites72.5%

          \[\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 rewrites76.7%

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

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

            if 4.9999999999999998e106 < (*.f64 M D)

            1. Initial program 68.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 h around inf

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

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

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

                \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
              4. 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}}} \]
              5. 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}} \]
              6. 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}} \]
              7. 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 rewrites56.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 rewrites57.0%

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

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

              Alternative 2: 82.0% accurate, 0.7× speedup?

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

                1. Initial program 98.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 rewrites98.5%

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

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

                  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 h around inf

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

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

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

                      \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                    4. 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}}} \]
                    5. 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}} \]
                    6. 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}} \]
                    7. 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 rewrites54.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 rewrites56.1%

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

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

                    Alternative 3: 82.8% accurate, 0.7× speedup?

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

                      1. Initial program 98.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 rewrites98.5%

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

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

                        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 h around inf

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

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

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

                            \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                          4. 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}}} \]
                          5. 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}} \]
                          6. 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}} \]
                          7. 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 rewrites54.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 rewrites56.1%

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

                              \[\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(\ell \cdot d\right) \cdot d}\right)}, 1\right)} \]
                          3. Recombined 2 regimes into one program.
                          4. Add Preprocessing

                          Alternative 4: 80.4% accurate, 0.8× speedup?

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

                            1. Initial program 63.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 rewrites44.3%

                              \[\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 rewrites43.3%

                                \[\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(\ell \cdot d\right) \cdot d}}, w0\right) \]
                              2. Step-by-step derivation
                                1. Applied rewrites46.9%

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

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

                                1. Initial program 86.7%

                                  \[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 rewrites93.7%

                                    \[\leadsto w0 \cdot \color{blue}{1} \]
                                5. Recombined 2 regimes into one program.
                                6. Add Preprocessing

                                Alternative 5: 83.8% accurate, 1.7× speedup?

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

                                  1. Initial program 82.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 rewrites75.7%

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

                                    if 5.0000000000000002e-136 < (*.f64 M D) < 4.99999999999999967e145

                                    1. Initial program 79.2%

                                      \[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 h around inf

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

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

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

                                        \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                      4. 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}}} \]
                                      5. 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}} \]
                                      6. 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}} \]
                                      7. 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 rewrites69.3%

                                      \[\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 rewrites72.8%

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

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

                                        if 4.99999999999999967e145 < (*.f64 M D)

                                        1. Initial program 62.7%

                                          \[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 h around inf

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

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

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

                                            \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                          4. 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}}} \]
                                          5. 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}} \]
                                          6. 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}} \]
                                          7. 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 rewrites57.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 rewrites58.3%

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

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

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

                                          Alternative 6: 83.2% accurate, 1.7× speedup?

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

                                            1. Initial program 82.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 rewrites75.7%

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

                                              if 5.0000000000000002e-136 < (*.f64 M D) < 1.9999999999999999e99

                                              1. Initial program 77.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 h around inf

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

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

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

                                                  \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                                4. 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}}} \]
                                                5. 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}} \]
                                                6. 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}} \]
                                                7. 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 rewrites73.4%

                                                \[\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 rewrites77.8%

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

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

                                                  if 1.9999999999999999e99 < (*.f64 M D)

                                                  1. Initial program 70.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 h around inf

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

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

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

                                                      \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                                    4. 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}}} \]
                                                    5. 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}} \]
                                                    6. 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}} \]
                                                    7. 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 rewrites56.5%

                                                    \[\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 rewrites59.9%

                                                      \[\leadsto w0 \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{D}{d} \cdot \color{blue}{\frac{\left(M \cdot M\right) \cdot D}{\ell \cdot d}}, 1\right)} \]
                                                  7. Recombined 3 regimes into one program.
                                                  8. Add Preprocessing

                                                  Alternative 7: 88.4% accurate, 2.1× speedup?

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

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

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

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

                                                      \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                                    4. 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}}} \]
                                                    5. 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}} \]
                                                    6. 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}} \]
                                                    7. 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 rewrites63.6%

                                                    \[\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 rewrites72.7%

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

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

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

                                                        Alternative 8: 86.3% accurate, 2.1× speedup?

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

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

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

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

                                                            \[\leadsto w0 \cdot \sqrt{h \cdot \color{blue}{\left(\frac{-1}{4} \cdot \frac{{D}^{2} \cdot {M}^{2}}{{d}^{2} \cdot \ell} + \frac{1}{h}\right)}} \]
                                                          4. 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}}} \]
                                                          5. 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}} \]
                                                          6. 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}} \]
                                                          7. 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 rewrites63.6%

                                                          \[\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 rewrites72.7%

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

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

                                                            Alternative 9: 76.8% accurate, 2.2× speedup?

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

                                                              1. Initial program 81.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 rewrites76.5%

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

                                                                if 1.9999999999999999e99 < (*.f64 M D)

                                                                1. Initial program 70.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 \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 rewrites41.3%

                                                                  \[\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 rewrites56.7%

                                                                    \[\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) \]
                                                                7. Recombined 2 regimes into one program.
                                                                8. Add Preprocessing

                                                                Alternative 10: 76.2% accurate, 2.4× speedup?

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

                                                                  1. Initial program 83.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 rewrites78.7%

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

                                                                    if 1.16000000000000001e-94 < M

                                                                    1. Initial program 72.2%

                                                                      \[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 rewrites46.1%

                                                                      \[\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 rewrites49.1%

                                                                        \[\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(\ell \cdot d\right) \cdot d}}, w0\right) \]
                                                                      2. Applied rewrites44.3%

                                                                        \[\leadsto \mathsf{fma}\left(D, \color{blue}{D \cdot \left(\frac{-0.125}{d \cdot d} \cdot \frac{\left(\left(M \cdot M\right) \cdot h\right) \cdot w0}{\ell}\right)}, w0\right) \]
                                                                    7. Recombined 2 regimes into one program.
                                                                    8. Add Preprocessing

                                                                    Alternative 11: 68.2% accurate, 26.2× speedup?

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

                                                                        \[\leadsto w0 \cdot \color{blue}{1} \]
                                                                      2. Add Preprocessing

                                                                      Reproduce

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