Toniolo and Linder, Equation (3a)

Percentage Accurate: 98.6% → 99.5%
Time: 8.6s
Alternatives: 6
Speedup: 2.1×

Specification

?
\[\begin{array}{l} \\ \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \end{array} \]
(FPCore (l Om kx ky)
 :precision binary64
 (sqrt
  (*
   (/ 1.0 2.0)
   (+
    1.0
    (/
     1.0
     (sqrt
      (+
       1.0
       (*
        (pow (/ (* 2.0 l) Om) 2.0)
        (+ (pow (sin kx) 2.0) (pow (sin ky) 2.0))))))))))
double code(double l, double Om, double kx, double ky) {
	return sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt((1.0 + (pow(((2.0 * l) / Om), 2.0) * (pow(sin(kx), 2.0) + pow(sin(ky), 2.0)))))))));
}
real(8) function code(l, om, kx, ky)
    real(8), intent (in) :: l
    real(8), intent (in) :: om
    real(8), intent (in) :: kx
    real(8), intent (in) :: ky
    code = sqrt(((1.0d0 / 2.0d0) * (1.0d0 + (1.0d0 / sqrt((1.0d0 + ((((2.0d0 * l) / om) ** 2.0d0) * ((sin(kx) ** 2.0d0) + (sin(ky) ** 2.0d0)))))))))
end function
public static double code(double l, double Om, double kx, double ky) {
	return Math.sqrt(((1.0 / 2.0) * (1.0 + (1.0 / Math.sqrt((1.0 + (Math.pow(((2.0 * l) / Om), 2.0) * (Math.pow(Math.sin(kx), 2.0) + Math.pow(Math.sin(ky), 2.0)))))))));
}
def code(l, Om, kx, ky):
	return math.sqrt(((1.0 / 2.0) * (1.0 + (1.0 / math.sqrt((1.0 + (math.pow(((2.0 * l) / Om), 2.0) * (math.pow(math.sin(kx), 2.0) + math.pow(math.sin(ky), 2.0)))))))))
function code(l, Om, kx, ky)
	return sqrt(Float64(Float64(1.0 / 2.0) * Float64(1.0 + Float64(1.0 / sqrt(Float64(1.0 + Float64((Float64(Float64(2.0 * l) / Om) ^ 2.0) * Float64((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0)))))))))
end
function tmp = code(l, Om, kx, ky)
	tmp = sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt((1.0 + ((((2.0 * l) / Om) ^ 2.0) * ((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0)))))))));
end
code[l_, Om_, kx_, ky_] := N[Sqrt[N[(N[(1.0 / 2.0), $MachinePrecision] * N[(1.0 + N[(1.0 / N[Sqrt[N[(1.0 + N[(N[Power[N[(N[(2.0 * l), $MachinePrecision] / Om), $MachinePrecision], 2.0], $MachinePrecision] * N[(N[Power[N[Sin[kx], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[ky], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}

\\
\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)}
\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 6 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: 98.6% accurate, 1.0× speedup?

\[\begin{array}{l} \\ \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \end{array} \]
(FPCore (l Om kx ky)
 :precision binary64
 (sqrt
  (*
   (/ 1.0 2.0)
   (+
    1.0
    (/
     1.0
     (sqrt
      (+
       1.0
       (*
        (pow (/ (* 2.0 l) Om) 2.0)
        (+ (pow (sin kx) 2.0) (pow (sin ky) 2.0))))))))))
double code(double l, double Om, double kx, double ky) {
	return sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt((1.0 + (pow(((2.0 * l) / Om), 2.0) * (pow(sin(kx), 2.0) + pow(sin(ky), 2.0)))))))));
}
real(8) function code(l, om, kx, ky)
    real(8), intent (in) :: l
    real(8), intent (in) :: om
    real(8), intent (in) :: kx
    real(8), intent (in) :: ky
    code = sqrt(((1.0d0 / 2.0d0) * (1.0d0 + (1.0d0 / sqrt((1.0d0 + ((((2.0d0 * l) / om) ** 2.0d0) * ((sin(kx) ** 2.0d0) + (sin(ky) ** 2.0d0)))))))))
end function
public static double code(double l, double Om, double kx, double ky) {
	return Math.sqrt(((1.0 / 2.0) * (1.0 + (1.0 / Math.sqrt((1.0 + (Math.pow(((2.0 * l) / Om), 2.0) * (Math.pow(Math.sin(kx), 2.0) + Math.pow(Math.sin(ky), 2.0)))))))));
}
def code(l, Om, kx, ky):
	return math.sqrt(((1.0 / 2.0) * (1.0 + (1.0 / math.sqrt((1.0 + (math.pow(((2.0 * l) / Om), 2.0) * (math.pow(math.sin(kx), 2.0) + math.pow(math.sin(ky), 2.0)))))))))
function code(l, Om, kx, ky)
	return sqrt(Float64(Float64(1.0 / 2.0) * Float64(1.0 + Float64(1.0 / sqrt(Float64(1.0 + Float64((Float64(Float64(2.0 * l) / Om) ^ 2.0) * Float64((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0)))))))))
end
function tmp = code(l, Om, kx, ky)
	tmp = sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt((1.0 + ((((2.0 * l) / Om) ^ 2.0) * ((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0)))))))));
end
code[l_, Om_, kx_, ky_] := N[Sqrt[N[(N[(1.0 / 2.0), $MachinePrecision] * N[(1.0 + N[(1.0 / N[Sqrt[N[(1.0 + N[(N[Power[N[(N[(2.0 * l), $MachinePrecision] / Om), $MachinePrecision], 2.0], $MachinePrecision] * N[(N[Power[N[Sin[kx], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[ky], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}

\\
\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)}
\end{array}

Alternative 1: 99.5% accurate, 2.1× speedup?

\[\begin{array}{l} ky_m = \left|ky\right| \\ kx_m = \left|kx\right| \\ [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\ \\ \sqrt{\left(\sqrt{\frac{1}{{\left(\frac{\ell \cdot 2}{Om} \cdot \sin ky\_m\right)}^{2} + 1}} + 1\right) \cdot \frac{1}{2}} \end{array} \]
ky_m = (fabs.f64 ky)
kx_m = (fabs.f64 kx)
NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
(FPCore (l Om kx_m ky_m)
 :precision binary64
 (sqrt
  (*
   (+ (sqrt (/ 1.0 (+ (pow (* (/ (* l 2.0) Om) (sin ky_m)) 2.0) 1.0))) 1.0)
   (/ 1.0 2.0))))
ky_m = fabs(ky);
kx_m = fabs(kx);
assert(l < Om && Om < kx_m && kx_m < ky_m);
double code(double l, double Om, double kx_m, double ky_m) {
	return sqrt(((sqrt((1.0 / (pow((((l * 2.0) / Om) * sin(ky_m)), 2.0) + 1.0))) + 1.0) * (1.0 / 2.0)));
}
ky_m = abs(ky)
kx_m = abs(kx)
NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
real(8) function code(l, om, kx_m, ky_m)
    real(8), intent (in) :: l
    real(8), intent (in) :: om
    real(8), intent (in) :: kx_m
    real(8), intent (in) :: ky_m
    code = sqrt(((sqrt((1.0d0 / (((((l * 2.0d0) / om) * sin(ky_m)) ** 2.0d0) + 1.0d0))) + 1.0d0) * (1.0d0 / 2.0d0)))
end function
ky_m = Math.abs(ky);
kx_m = Math.abs(kx);
assert l < Om && Om < kx_m && kx_m < ky_m;
public static double code(double l, double Om, double kx_m, double ky_m) {
	return Math.sqrt(((Math.sqrt((1.0 / (Math.pow((((l * 2.0) / Om) * Math.sin(ky_m)), 2.0) + 1.0))) + 1.0) * (1.0 / 2.0)));
}
ky_m = math.fabs(ky)
kx_m = math.fabs(kx)
[l, Om, kx_m, ky_m] = sort([l, Om, kx_m, ky_m])
def code(l, Om, kx_m, ky_m):
	return math.sqrt(((math.sqrt((1.0 / (math.pow((((l * 2.0) / Om) * math.sin(ky_m)), 2.0) + 1.0))) + 1.0) * (1.0 / 2.0)))
ky_m = abs(ky)
kx_m = abs(kx)
l, Om, kx_m, ky_m = sort([l, Om, kx_m, ky_m])
function code(l, Om, kx_m, ky_m)
	return sqrt(Float64(Float64(sqrt(Float64(1.0 / Float64((Float64(Float64(Float64(l * 2.0) / Om) * sin(ky_m)) ^ 2.0) + 1.0))) + 1.0) * Float64(1.0 / 2.0)))
end
ky_m = abs(ky);
kx_m = abs(kx);
l, Om, kx_m, ky_m = num2cell(sort([l, Om, kx_m, ky_m])){:}
function tmp = code(l, Om, kx_m, ky_m)
	tmp = sqrt(((sqrt((1.0 / (((((l * 2.0) / Om) * sin(ky_m)) ^ 2.0) + 1.0))) + 1.0) * (1.0 / 2.0)));
end
ky_m = N[Abs[ky], $MachinePrecision]
kx_m = N[Abs[kx], $MachinePrecision]
NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
code[l_, Om_, kx$95$m_, ky$95$m_] := N[Sqrt[N[(N[(N[Sqrt[N[(1.0 / N[(N[Power[N[(N[(N[(l * 2.0), $MachinePrecision] / Om), $MachinePrecision] * N[Sin[ky$95$m], $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] + 1.0), $MachinePrecision] * N[(1.0 / 2.0), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}
ky_m = \left|ky\right|
\\
kx_m = \left|kx\right|
\\
[l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\
\\
\sqrt{\left(\sqrt{\frac{1}{{\left(\frac{\ell \cdot 2}{Om} \cdot \sin ky\_m\right)}^{2} + 1}} + 1\right) \cdot \frac{1}{2}}
\end{array}
Derivation
  1. Initial program 99.2%

    \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
  2. Add Preprocessing
  3. Taylor expanded in kx around 0

    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
  4. Step-by-step derivation
    1. lower-sqrt.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
    2. lower-/.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\color{blue}{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
    3. +-commutativeN/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}} + 1}}}\right)} \]
    4. associate-/l*N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{4 \cdot \color{blue}{\left({\ell}^{2} \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}\right)} + 1}}\right)} \]
    5. associate-*r*N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\left(4 \cdot {\ell}^{2}\right) \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}} + 1}}\right)} \]
    6. lower-fma.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\mathsf{fma}\left(4 \cdot {\ell}^{2}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}}\right)} \]
    7. lower-*.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(\color{blue}{4 \cdot {\ell}^{2}}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
    8. unpow2N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
    9. lower-*.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
    10. lower-/.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \color{blue}{\frac{{\sin ky}^{2}}{{Om}^{2}}}, 1\right)}}\right)} \]
    11. lower-pow.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{\color{blue}{{\sin ky}^{2}}}{{Om}^{2}}, 1\right)}}\right)} \]
    12. lower-sin.f64N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\color{blue}{\sin ky}}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
    13. unpow2N/A

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
    14. lower-*.f6481.3

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
  5. Applied rewrites81.3%

    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{Om \cdot Om}, 1\right)}}}\right)} \]
  6. Step-by-step derivation
    1. Applied rewrites93.0%

      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{{\left(\sin ky \cdot \frac{\ell \cdot 2}{Om}\right)}^{2} + 1}}\right)} \]
    2. Final simplification93.0%

      \[\leadsto \sqrt{\left(\sqrt{\frac{1}{{\left(\frac{\ell \cdot 2}{Om} \cdot \sin ky\right)}^{2} + 1}} + 1\right) \cdot \frac{1}{2}} \]
    3. Add Preprocessing

    Alternative 2: 98.4% accurate, 0.8× speedup?

    \[\begin{array}{l} ky_m = \left|ky\right| \\ kx_m = \left|kx\right| \\ [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\ \\ \begin{array}{l} \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 0.4:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\left(\frac{Om}{\sin ky\_m \cdot \ell} \cdot 0.5 + 1\right) \cdot \frac{1}{2}}\\ \end{array} \end{array} \]
    ky_m = (fabs.f64 ky)
    kx_m = (fabs.f64 kx)
    NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
    (FPCore (l Om kx_m ky_m)
     :precision binary64
     (if (<=
          (*
           (+ (pow (sin ky_m) 2.0) (pow (sin kx_m) 2.0))
           (pow (/ (* l 2.0) Om) 2.0))
          0.4)
       (sqrt (+ 0.5 0.5))
       (sqrt (* (+ (* (/ Om (* (sin ky_m) l)) 0.5) 1.0) (/ 1.0 2.0)))))
    ky_m = fabs(ky);
    kx_m = fabs(kx);
    assert(l < Om && Om < kx_m && kx_m < ky_m);
    double code(double l, double Om, double kx_m, double ky_m) {
    	double tmp;
    	if (((pow(sin(ky_m), 2.0) + pow(sin(kx_m), 2.0)) * pow(((l * 2.0) / Om), 2.0)) <= 0.4) {
    		tmp = sqrt((0.5 + 0.5));
    	} else {
    		tmp = sqrt(((((Om / (sin(ky_m) * l)) * 0.5) + 1.0) * (1.0 / 2.0)));
    	}
    	return tmp;
    }
    
    ky_m = abs(ky)
    kx_m = abs(kx)
    NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
    real(8) function code(l, om, kx_m, ky_m)
        real(8), intent (in) :: l
        real(8), intent (in) :: om
        real(8), intent (in) :: kx_m
        real(8), intent (in) :: ky_m
        real(8) :: tmp
        if ((((sin(ky_m) ** 2.0d0) + (sin(kx_m) ** 2.0d0)) * (((l * 2.0d0) / om) ** 2.0d0)) <= 0.4d0) then
            tmp = sqrt((0.5d0 + 0.5d0))
        else
            tmp = sqrt(((((om / (sin(ky_m) * l)) * 0.5d0) + 1.0d0) * (1.0d0 / 2.0d0)))
        end if
        code = tmp
    end function
    
    ky_m = Math.abs(ky);
    kx_m = Math.abs(kx);
    assert l < Om && Om < kx_m && kx_m < ky_m;
    public static double code(double l, double Om, double kx_m, double ky_m) {
    	double tmp;
    	if (((Math.pow(Math.sin(ky_m), 2.0) + Math.pow(Math.sin(kx_m), 2.0)) * Math.pow(((l * 2.0) / Om), 2.0)) <= 0.4) {
    		tmp = Math.sqrt((0.5 + 0.5));
    	} else {
    		tmp = Math.sqrt(((((Om / (Math.sin(ky_m) * l)) * 0.5) + 1.0) * (1.0 / 2.0)));
    	}
    	return tmp;
    }
    
    ky_m = math.fabs(ky)
    kx_m = math.fabs(kx)
    [l, Om, kx_m, ky_m] = sort([l, Om, kx_m, ky_m])
    def code(l, Om, kx_m, ky_m):
    	tmp = 0
    	if ((math.pow(math.sin(ky_m), 2.0) + math.pow(math.sin(kx_m), 2.0)) * math.pow(((l * 2.0) / Om), 2.0)) <= 0.4:
    		tmp = math.sqrt((0.5 + 0.5))
    	else:
    		tmp = math.sqrt(((((Om / (math.sin(ky_m) * l)) * 0.5) + 1.0) * (1.0 / 2.0)))
    	return tmp
    
    ky_m = abs(ky)
    kx_m = abs(kx)
    l, Om, kx_m, ky_m = sort([l, Om, kx_m, ky_m])
    function code(l, Om, kx_m, ky_m)
    	tmp = 0.0
    	if (Float64(Float64((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (Float64(Float64(l * 2.0) / Om) ^ 2.0)) <= 0.4)
    		tmp = sqrt(Float64(0.5 + 0.5));
    	else
    		tmp = sqrt(Float64(Float64(Float64(Float64(Om / Float64(sin(ky_m) * l)) * 0.5) + 1.0) * Float64(1.0 / 2.0)));
    	end
    	return tmp
    end
    
    ky_m = abs(ky);
    kx_m = abs(kx);
    l, Om, kx_m, ky_m = num2cell(sort([l, Om, kx_m, ky_m])){:}
    function tmp_2 = code(l, Om, kx_m, ky_m)
    	tmp = 0.0;
    	if ((((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (((l * 2.0) / Om) ^ 2.0)) <= 0.4)
    		tmp = sqrt((0.5 + 0.5));
    	else
    		tmp = sqrt(((((Om / (sin(ky_m) * l)) * 0.5) + 1.0) * (1.0 / 2.0)));
    	end
    	tmp_2 = tmp;
    end
    
    ky_m = N[Abs[ky], $MachinePrecision]
    kx_m = N[Abs[kx], $MachinePrecision]
    NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
    code[l_, Om_, kx$95$m_, ky$95$m_] := If[LessEqual[N[(N[(N[Power[N[Sin[ky$95$m], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[kx$95$m], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] * N[Power[N[(N[(l * 2.0), $MachinePrecision] / Om), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision], 0.4], N[Sqrt[N[(0.5 + 0.5), $MachinePrecision]], $MachinePrecision], N[Sqrt[N[(N[(N[(N[(Om / N[(N[Sin[ky$95$m], $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision] + 1.0), $MachinePrecision] * N[(1.0 / 2.0), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]
    
    \begin{array}{l}
    ky_m = \left|ky\right|
    \\
    kx_m = \left|kx\right|
    \\
    [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\
    \\
    \begin{array}{l}
    \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 0.4:\\
    \;\;\;\;\sqrt{0.5 + 0.5}\\
    
    \mathbf{else}:\\
    \;\;\;\;\sqrt{\left(\frac{Om}{\sin ky\_m \cdot \ell} \cdot 0.5 + 1\right) \cdot \frac{1}{2}}\\
    
    
    \end{array}
    \end{array}
    
    Derivation
    1. Split input into 2 regimes
    2. if (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64)))) < 0.40000000000000002

      1. Initial program 100.0%

        \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
      2. Add Preprocessing
      3. Taylor expanded in l around 0

        \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
      4. Step-by-step derivation
        1. Applied rewrites98.0%

          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
        2. Step-by-step derivation
          1. lift-*.f64N/A

            \[\leadsto \sqrt{\color{blue}{\frac{1}{2} \cdot \left(1 + \frac{1}{1}\right)}} \]
          2. lift-/.f64N/A

            \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
          3. metadata-evalN/A

            \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
          4. lift-+.f64N/A

            \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(1 + \frac{1}{1}\right)}} \]
          5. +-commutativeN/A

            \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(\frac{1}{1} + 1\right)}} \]
          6. distribute-rgt-inN/A

            \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + 1 \cdot \frac{1}{2}}} \]
          7. metadata-evalN/A

            \[\leadsto \sqrt{\frac{1}{1} \cdot \frac{1}{2} + \color{blue}{\frac{1}{2}}} \]
          8. lower-+.f64N/A

            \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + \frac{1}{2}}} \]
        3. Applied rewrites98.0%

          \[\leadsto \color{blue}{\sqrt{\frac{0.5}{1} + 0.5}} \]
        4. Taylor expanded in l around 0

          \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} + \frac{1}{2}} \]
        5. Step-by-step derivation
          1. Applied rewrites98.0%

            \[\leadsto \sqrt{\color{blue}{0.5} + 0.5} \]

          if 0.40000000000000002 < (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64))))

          1. Initial program 98.1%

            \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
          2. Add Preprocessing
          3. Taylor expanded in kx around 0

            \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
          4. Step-by-step derivation
            1. lower-sqrt.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
            2. lower-/.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\color{blue}{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
            3. +-commutativeN/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}} + 1}}}\right)} \]
            4. associate-/l*N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{4 \cdot \color{blue}{\left({\ell}^{2} \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}\right)} + 1}}\right)} \]
            5. associate-*r*N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\left(4 \cdot {\ell}^{2}\right) \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}} + 1}}\right)} \]
            6. lower-fma.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\mathsf{fma}\left(4 \cdot {\ell}^{2}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}}\right)} \]
            7. lower-*.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(\color{blue}{4 \cdot {\ell}^{2}}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
            8. unpow2N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
            9. lower-*.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
            10. lower-/.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \color{blue}{\frac{{\sin ky}^{2}}{{Om}^{2}}}, 1\right)}}\right)} \]
            11. lower-pow.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{\color{blue}{{\sin ky}^{2}}}{{Om}^{2}}, 1\right)}}\right)} \]
            12. lower-sin.f64N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\color{blue}{\sin ky}}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
            13. unpow2N/A

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
            14. lower-*.f6469.7

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
          5. Applied rewrites69.7%

            \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{Om \cdot Om}, 1\right)}}}\right)} \]
          6. Taylor expanded in l around inf

            \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{2} \cdot \color{blue}{\frac{Om}{\ell \cdot \sin ky}}\right)} \]
          7. Step-by-step derivation
            1. Applied rewrites80.6%

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{Om}{\sin ky \cdot \ell} \cdot \color{blue}{0.5}\right)} \]
          8. Recombined 2 regimes into one program.
          9. Final simplification90.7%

            \[\leadsto \begin{array}{l} \mathbf{if}\;\left({\sin ky}^{2} + {\sin kx}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 0.4:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\left(\frac{Om}{\sin ky \cdot \ell} \cdot 0.5 + 1\right) \cdot \frac{1}{2}}\\ \end{array} \]
          10. Add Preprocessing

          Alternative 3: 98.4% accurate, 0.8× speedup?

          \[\begin{array}{l} ky_m = \left|ky\right| \\ kx_m = \left|kx\right| \\ [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\ \\ \begin{array}{l} \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 2:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\left(-0.5 \cdot \frac{Om}{\sin ky\_m \cdot \ell} + 1\right) \cdot \frac{1}{2}}\\ \end{array} \end{array} \]
          ky_m = (fabs.f64 ky)
          kx_m = (fabs.f64 kx)
          NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
          (FPCore (l Om kx_m ky_m)
           :precision binary64
           (if (<=
                (*
                 (+ (pow (sin ky_m) 2.0) (pow (sin kx_m) 2.0))
                 (pow (/ (* l 2.0) Om) 2.0))
                2.0)
             (sqrt (+ 0.5 0.5))
             (sqrt (* (+ (* -0.5 (/ Om (* (sin ky_m) l))) 1.0) (/ 1.0 2.0)))))
          ky_m = fabs(ky);
          kx_m = fabs(kx);
          assert(l < Om && Om < kx_m && kx_m < ky_m);
          double code(double l, double Om, double kx_m, double ky_m) {
          	double tmp;
          	if (((pow(sin(ky_m), 2.0) + pow(sin(kx_m), 2.0)) * pow(((l * 2.0) / Om), 2.0)) <= 2.0) {
          		tmp = sqrt((0.5 + 0.5));
          	} else {
          		tmp = sqrt((((-0.5 * (Om / (sin(ky_m) * l))) + 1.0) * (1.0 / 2.0)));
          	}
          	return tmp;
          }
          
          ky_m = abs(ky)
          kx_m = abs(kx)
          NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
          real(8) function code(l, om, kx_m, ky_m)
              real(8), intent (in) :: l
              real(8), intent (in) :: om
              real(8), intent (in) :: kx_m
              real(8), intent (in) :: ky_m
              real(8) :: tmp
              if ((((sin(ky_m) ** 2.0d0) + (sin(kx_m) ** 2.0d0)) * (((l * 2.0d0) / om) ** 2.0d0)) <= 2.0d0) then
                  tmp = sqrt((0.5d0 + 0.5d0))
              else
                  tmp = sqrt(((((-0.5d0) * (om / (sin(ky_m) * l))) + 1.0d0) * (1.0d0 / 2.0d0)))
              end if
              code = tmp
          end function
          
          ky_m = Math.abs(ky);
          kx_m = Math.abs(kx);
          assert l < Om && Om < kx_m && kx_m < ky_m;
          public static double code(double l, double Om, double kx_m, double ky_m) {
          	double tmp;
          	if (((Math.pow(Math.sin(ky_m), 2.0) + Math.pow(Math.sin(kx_m), 2.0)) * Math.pow(((l * 2.0) / Om), 2.0)) <= 2.0) {
          		tmp = Math.sqrt((0.5 + 0.5));
          	} else {
          		tmp = Math.sqrt((((-0.5 * (Om / (Math.sin(ky_m) * l))) + 1.0) * (1.0 / 2.0)));
          	}
          	return tmp;
          }
          
          ky_m = math.fabs(ky)
          kx_m = math.fabs(kx)
          [l, Om, kx_m, ky_m] = sort([l, Om, kx_m, ky_m])
          def code(l, Om, kx_m, ky_m):
          	tmp = 0
          	if ((math.pow(math.sin(ky_m), 2.0) + math.pow(math.sin(kx_m), 2.0)) * math.pow(((l * 2.0) / Om), 2.0)) <= 2.0:
          		tmp = math.sqrt((0.5 + 0.5))
          	else:
          		tmp = math.sqrt((((-0.5 * (Om / (math.sin(ky_m) * l))) + 1.0) * (1.0 / 2.0)))
          	return tmp
          
          ky_m = abs(ky)
          kx_m = abs(kx)
          l, Om, kx_m, ky_m = sort([l, Om, kx_m, ky_m])
          function code(l, Om, kx_m, ky_m)
          	tmp = 0.0
          	if (Float64(Float64((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (Float64(Float64(l * 2.0) / Om) ^ 2.0)) <= 2.0)
          		tmp = sqrt(Float64(0.5 + 0.5));
          	else
          		tmp = sqrt(Float64(Float64(Float64(-0.5 * Float64(Om / Float64(sin(ky_m) * l))) + 1.0) * Float64(1.0 / 2.0)));
          	end
          	return tmp
          end
          
          ky_m = abs(ky);
          kx_m = abs(kx);
          l, Om, kx_m, ky_m = num2cell(sort([l, Om, kx_m, ky_m])){:}
          function tmp_2 = code(l, Om, kx_m, ky_m)
          	tmp = 0.0;
          	if ((((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (((l * 2.0) / Om) ^ 2.0)) <= 2.0)
          		tmp = sqrt((0.5 + 0.5));
          	else
          		tmp = sqrt((((-0.5 * (Om / (sin(ky_m) * l))) + 1.0) * (1.0 / 2.0)));
          	end
          	tmp_2 = tmp;
          end
          
          ky_m = N[Abs[ky], $MachinePrecision]
          kx_m = N[Abs[kx], $MachinePrecision]
          NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
          code[l_, Om_, kx$95$m_, ky$95$m_] := If[LessEqual[N[(N[(N[Power[N[Sin[ky$95$m], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[kx$95$m], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] * N[Power[N[(N[(l * 2.0), $MachinePrecision] / Om), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision], 2.0], N[Sqrt[N[(0.5 + 0.5), $MachinePrecision]], $MachinePrecision], N[Sqrt[N[(N[(N[(-0.5 * N[(Om / N[(N[Sin[ky$95$m], $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] * N[(1.0 / 2.0), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]
          
          \begin{array}{l}
          ky_m = \left|ky\right|
          \\
          kx_m = \left|kx\right|
          \\
          [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\
          \\
          \begin{array}{l}
          \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 2:\\
          \;\;\;\;\sqrt{0.5 + 0.5}\\
          
          \mathbf{else}:\\
          \;\;\;\;\sqrt{\left(-0.5 \cdot \frac{Om}{\sin ky\_m \cdot \ell} + 1\right) \cdot \frac{1}{2}}\\
          
          
          \end{array}
          \end{array}
          
          Derivation
          1. Split input into 2 regimes
          2. if (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64)))) < 2

            1. Initial program 100.0%

              \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
            2. Add Preprocessing
            3. Taylor expanded in l around 0

              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
            4. Step-by-step derivation
              1. Applied rewrites97.5%

                \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
              2. Step-by-step derivation
                1. lift-*.f64N/A

                  \[\leadsto \sqrt{\color{blue}{\frac{1}{2} \cdot \left(1 + \frac{1}{1}\right)}} \]
                2. lift-/.f64N/A

                  \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
                3. metadata-evalN/A

                  \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
                4. lift-+.f64N/A

                  \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(1 + \frac{1}{1}\right)}} \]
                5. +-commutativeN/A

                  \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(\frac{1}{1} + 1\right)}} \]
                6. distribute-rgt-inN/A

                  \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + 1 \cdot \frac{1}{2}}} \]
                7. metadata-evalN/A

                  \[\leadsto \sqrt{\frac{1}{1} \cdot \frac{1}{2} + \color{blue}{\frac{1}{2}}} \]
                8. lower-+.f64N/A

                  \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + \frac{1}{2}}} \]
              3. Applied rewrites97.5%

                \[\leadsto \color{blue}{\sqrt{\frac{0.5}{1} + 0.5}} \]
              4. Taylor expanded in l around 0

                \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} + \frac{1}{2}} \]
              5. Step-by-step derivation
                1. Applied rewrites97.5%

                  \[\leadsto \sqrt{\color{blue}{0.5} + 0.5} \]

                if 2 < (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64))))

                1. Initial program 98.1%

                  \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
                2. Add Preprocessing
                3. Taylor expanded in kx around 0

                  \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
                4. Step-by-step derivation
                  1. lower-sqrt.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
                  2. lower-/.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\color{blue}{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
                  3. +-commutativeN/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}} + 1}}}\right)} \]
                  4. associate-/l*N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{4 \cdot \color{blue}{\left({\ell}^{2} \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}\right)} + 1}}\right)} \]
                  5. associate-*r*N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\left(4 \cdot {\ell}^{2}\right) \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}} + 1}}\right)} \]
                  6. lower-fma.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\mathsf{fma}\left(4 \cdot {\ell}^{2}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}}\right)} \]
                  7. lower-*.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(\color{blue}{4 \cdot {\ell}^{2}}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                  8. unpow2N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                  9. lower-*.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                  10. lower-/.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \color{blue}{\frac{{\sin ky}^{2}}{{Om}^{2}}}, 1\right)}}\right)} \]
                  11. lower-pow.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{\color{blue}{{\sin ky}^{2}}}{{Om}^{2}}, 1\right)}}\right)} \]
                  12. lower-sin.f64N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\color{blue}{\sin ky}}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                  13. unpow2N/A

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
                  14. lower-*.f6470.1

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
                5. Applied rewrites70.1%

                  \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{Om \cdot Om}, 1\right)}}}\right)} \]
                6. Taylor expanded in l around -inf

                  \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{-1}{2} \cdot \color{blue}{\frac{Om}{\ell \cdot \sin ky}}\right)} \]
                7. Step-by-step derivation
                  1. Applied rewrites81.4%

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{Om}{\sin ky \cdot \ell} \cdot \color{blue}{-0.5}\right)} \]
                8. Recombined 2 regimes into one program.
                9. Final simplification90.8%

                  \[\leadsto \begin{array}{l} \mathbf{if}\;\left({\sin ky}^{2} + {\sin kx}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 2:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\left(-0.5 \cdot \frac{Om}{\sin ky \cdot \ell} + 1\right) \cdot \frac{1}{2}}\\ \end{array} \]
                10. Add Preprocessing

                Alternative 4: 98.3% accurate, 1.0× speedup?

                \[\begin{array}{l} ky_m = \left|ky\right| \\ kx_m = \left|kx\right| \\ [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\ \\ \begin{array}{l} \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 0.4:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\left(\frac{Om}{ky\_m \cdot \ell} \cdot 0.5 + 1\right) \cdot \frac{1}{2}}\\ \end{array} \end{array} \]
                ky_m = (fabs.f64 ky)
                kx_m = (fabs.f64 kx)
                NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                (FPCore (l Om kx_m ky_m)
                 :precision binary64
                 (if (<=
                      (*
                       (+ (pow (sin ky_m) 2.0) (pow (sin kx_m) 2.0))
                       (pow (/ (* l 2.0) Om) 2.0))
                      0.4)
                   (sqrt (+ 0.5 0.5))
                   (sqrt (* (+ (* (/ Om (* ky_m l)) 0.5) 1.0) (/ 1.0 2.0)))))
                ky_m = fabs(ky);
                kx_m = fabs(kx);
                assert(l < Om && Om < kx_m && kx_m < ky_m);
                double code(double l, double Om, double kx_m, double ky_m) {
                	double tmp;
                	if (((pow(sin(ky_m), 2.0) + pow(sin(kx_m), 2.0)) * pow(((l * 2.0) / Om), 2.0)) <= 0.4) {
                		tmp = sqrt((0.5 + 0.5));
                	} else {
                		tmp = sqrt(((((Om / (ky_m * l)) * 0.5) + 1.0) * (1.0 / 2.0)));
                	}
                	return tmp;
                }
                
                ky_m = abs(ky)
                kx_m = abs(kx)
                NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                real(8) function code(l, om, kx_m, ky_m)
                    real(8), intent (in) :: l
                    real(8), intent (in) :: om
                    real(8), intent (in) :: kx_m
                    real(8), intent (in) :: ky_m
                    real(8) :: tmp
                    if ((((sin(ky_m) ** 2.0d0) + (sin(kx_m) ** 2.0d0)) * (((l * 2.0d0) / om) ** 2.0d0)) <= 0.4d0) then
                        tmp = sqrt((0.5d0 + 0.5d0))
                    else
                        tmp = sqrt(((((om / (ky_m * l)) * 0.5d0) + 1.0d0) * (1.0d0 / 2.0d0)))
                    end if
                    code = tmp
                end function
                
                ky_m = Math.abs(ky);
                kx_m = Math.abs(kx);
                assert l < Om && Om < kx_m && kx_m < ky_m;
                public static double code(double l, double Om, double kx_m, double ky_m) {
                	double tmp;
                	if (((Math.pow(Math.sin(ky_m), 2.0) + Math.pow(Math.sin(kx_m), 2.0)) * Math.pow(((l * 2.0) / Om), 2.0)) <= 0.4) {
                		tmp = Math.sqrt((0.5 + 0.5));
                	} else {
                		tmp = Math.sqrt(((((Om / (ky_m * l)) * 0.5) + 1.0) * (1.0 / 2.0)));
                	}
                	return tmp;
                }
                
                ky_m = math.fabs(ky)
                kx_m = math.fabs(kx)
                [l, Om, kx_m, ky_m] = sort([l, Om, kx_m, ky_m])
                def code(l, Om, kx_m, ky_m):
                	tmp = 0
                	if ((math.pow(math.sin(ky_m), 2.0) + math.pow(math.sin(kx_m), 2.0)) * math.pow(((l * 2.0) / Om), 2.0)) <= 0.4:
                		tmp = math.sqrt((0.5 + 0.5))
                	else:
                		tmp = math.sqrt(((((Om / (ky_m * l)) * 0.5) + 1.0) * (1.0 / 2.0)))
                	return tmp
                
                ky_m = abs(ky)
                kx_m = abs(kx)
                l, Om, kx_m, ky_m = sort([l, Om, kx_m, ky_m])
                function code(l, Om, kx_m, ky_m)
                	tmp = 0.0
                	if (Float64(Float64((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (Float64(Float64(l * 2.0) / Om) ^ 2.0)) <= 0.4)
                		tmp = sqrt(Float64(0.5 + 0.5));
                	else
                		tmp = sqrt(Float64(Float64(Float64(Float64(Om / Float64(ky_m * l)) * 0.5) + 1.0) * Float64(1.0 / 2.0)));
                	end
                	return tmp
                end
                
                ky_m = abs(ky);
                kx_m = abs(kx);
                l, Om, kx_m, ky_m = num2cell(sort([l, Om, kx_m, ky_m])){:}
                function tmp_2 = code(l, Om, kx_m, ky_m)
                	tmp = 0.0;
                	if ((((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (((l * 2.0) / Om) ^ 2.0)) <= 0.4)
                		tmp = sqrt((0.5 + 0.5));
                	else
                		tmp = sqrt(((((Om / (ky_m * l)) * 0.5) + 1.0) * (1.0 / 2.0)));
                	end
                	tmp_2 = tmp;
                end
                
                ky_m = N[Abs[ky], $MachinePrecision]
                kx_m = N[Abs[kx], $MachinePrecision]
                NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                code[l_, Om_, kx$95$m_, ky$95$m_] := If[LessEqual[N[(N[(N[Power[N[Sin[ky$95$m], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[kx$95$m], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] * N[Power[N[(N[(l * 2.0), $MachinePrecision] / Om), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision], 0.4], N[Sqrt[N[(0.5 + 0.5), $MachinePrecision]], $MachinePrecision], N[Sqrt[N[(N[(N[(N[(Om / N[(ky$95$m * l), $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision] + 1.0), $MachinePrecision] * N[(1.0 / 2.0), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]
                
                \begin{array}{l}
                ky_m = \left|ky\right|
                \\
                kx_m = \left|kx\right|
                \\
                [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\
                \\
                \begin{array}{l}
                \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 0.4:\\
                \;\;\;\;\sqrt{0.5 + 0.5}\\
                
                \mathbf{else}:\\
                \;\;\;\;\sqrt{\left(\frac{Om}{ky\_m \cdot \ell} \cdot 0.5 + 1\right) \cdot \frac{1}{2}}\\
                
                
                \end{array}
                \end{array}
                
                Derivation
                1. Split input into 2 regimes
                2. if (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64)))) < 0.40000000000000002

                  1. Initial program 100.0%

                    \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
                  2. Add Preprocessing
                  3. Taylor expanded in l around 0

                    \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                  4. Step-by-step derivation
                    1. Applied rewrites98.0%

                      \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                    2. Step-by-step derivation
                      1. lift-*.f64N/A

                        \[\leadsto \sqrt{\color{blue}{\frac{1}{2} \cdot \left(1 + \frac{1}{1}\right)}} \]
                      2. lift-/.f64N/A

                        \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
                      3. metadata-evalN/A

                        \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
                      4. lift-+.f64N/A

                        \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(1 + \frac{1}{1}\right)}} \]
                      5. +-commutativeN/A

                        \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(\frac{1}{1} + 1\right)}} \]
                      6. distribute-rgt-inN/A

                        \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + 1 \cdot \frac{1}{2}}} \]
                      7. metadata-evalN/A

                        \[\leadsto \sqrt{\frac{1}{1} \cdot \frac{1}{2} + \color{blue}{\frac{1}{2}}} \]
                      8. lower-+.f64N/A

                        \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + \frac{1}{2}}} \]
                    3. Applied rewrites98.0%

                      \[\leadsto \color{blue}{\sqrt{\frac{0.5}{1} + 0.5}} \]
                    4. Taylor expanded in l around 0

                      \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} + \frac{1}{2}} \]
                    5. Step-by-step derivation
                      1. Applied rewrites98.0%

                        \[\leadsto \sqrt{\color{blue}{0.5} + 0.5} \]

                      if 0.40000000000000002 < (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64))))

                      1. Initial program 98.1%

                        \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
                      2. Add Preprocessing
                      3. Taylor expanded in kx around 0

                        \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
                      4. Step-by-step derivation
                        1. lower-sqrt.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
                        2. lower-/.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\color{blue}{\frac{1}{1 + 4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}}}}}\right)} \]
                        3. +-commutativeN/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{4 \cdot \frac{{\ell}^{2} \cdot {\sin ky}^{2}}{{Om}^{2}} + 1}}}\right)} \]
                        4. associate-/l*N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{4 \cdot \color{blue}{\left({\ell}^{2} \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}\right)} + 1}}\right)} \]
                        5. associate-*r*N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\left(4 \cdot {\ell}^{2}\right) \cdot \frac{{\sin ky}^{2}}{{Om}^{2}}} + 1}}\right)} \]
                        6. lower-fma.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\color{blue}{\mathsf{fma}\left(4 \cdot {\ell}^{2}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}}\right)} \]
                        7. lower-*.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(\color{blue}{4 \cdot {\ell}^{2}}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                        8. unpow2N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                        9. lower-*.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \color{blue}{\left(\ell \cdot \ell\right)}, \frac{{\sin ky}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                        10. lower-/.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \color{blue}{\frac{{\sin ky}^{2}}{{Om}^{2}}}, 1\right)}}\right)} \]
                        11. lower-pow.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{\color{blue}{{\sin ky}^{2}}}{{Om}^{2}}, 1\right)}}\right)} \]
                        12. lower-sin.f64N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\color{blue}{\sin ky}}^{2}}{{Om}^{2}}, 1\right)}}\right)} \]
                        13. unpow2N/A

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
                        14. lower-*.f6469.7

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{\color{blue}{Om \cdot Om}}, 1\right)}}\right)} \]
                      5. Applied rewrites69.7%

                        \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \color{blue}{\sqrt{\frac{1}{\mathsf{fma}\left(4 \cdot \left(\ell \cdot \ell\right), \frac{{\sin ky}^{2}}{Om \cdot Om}, 1\right)}}}\right)} \]
                      6. Taylor expanded in l around inf

                        \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{2} \cdot \color{blue}{\frac{Om}{\ell \cdot \sin ky}}\right)} \]
                      7. Step-by-step derivation
                        1. Applied rewrites80.6%

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{Om}{\sin ky \cdot \ell} \cdot \color{blue}{0.5}\right)} \]
                        2. Taylor expanded in ky around 0

                          \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{Om}{ky \cdot \ell} \cdot \frac{1}{2}\right)} \]
                        3. Step-by-step derivation
                          1. Applied rewrites80.1%

                            \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{Om}{\ell \cdot ky} \cdot 0.5\right)} \]
                        4. Recombined 2 regimes into one program.
                        5. Final simplification90.5%

                          \[\leadsto \begin{array}{l} \mathbf{if}\;\left({\sin ky}^{2} + {\sin kx}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 0.4:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\left(\frac{Om}{ky \cdot \ell} \cdot 0.5 + 1\right) \cdot \frac{1}{2}}\\ \end{array} \]
                        6. Add Preprocessing

                        Alternative 5: 98.2% accurate, 1.1× speedup?

                        \[\begin{array}{l} ky_m = \left|ky\right| \\ kx_m = \left|kx\right| \\ [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\ \\ \begin{array}{l} \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 3.8:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{0.5}\\ \end{array} \end{array} \]
                        ky_m = (fabs.f64 ky)
                        kx_m = (fabs.f64 kx)
                        NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                        (FPCore (l Om kx_m ky_m)
                         :precision binary64
                         (if (<=
                              (*
                               (+ (pow (sin ky_m) 2.0) (pow (sin kx_m) 2.0))
                               (pow (/ (* l 2.0) Om) 2.0))
                              3.8)
                           (sqrt (+ 0.5 0.5))
                           (sqrt 0.5)))
                        ky_m = fabs(ky);
                        kx_m = fabs(kx);
                        assert(l < Om && Om < kx_m && kx_m < ky_m);
                        double code(double l, double Om, double kx_m, double ky_m) {
                        	double tmp;
                        	if (((pow(sin(ky_m), 2.0) + pow(sin(kx_m), 2.0)) * pow(((l * 2.0) / Om), 2.0)) <= 3.8) {
                        		tmp = sqrt((0.5 + 0.5));
                        	} else {
                        		tmp = sqrt(0.5);
                        	}
                        	return tmp;
                        }
                        
                        ky_m = abs(ky)
                        kx_m = abs(kx)
                        NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                        real(8) function code(l, om, kx_m, ky_m)
                            real(8), intent (in) :: l
                            real(8), intent (in) :: om
                            real(8), intent (in) :: kx_m
                            real(8), intent (in) :: ky_m
                            real(8) :: tmp
                            if ((((sin(ky_m) ** 2.0d0) + (sin(kx_m) ** 2.0d0)) * (((l * 2.0d0) / om) ** 2.0d0)) <= 3.8d0) then
                                tmp = sqrt((0.5d0 + 0.5d0))
                            else
                                tmp = sqrt(0.5d0)
                            end if
                            code = tmp
                        end function
                        
                        ky_m = Math.abs(ky);
                        kx_m = Math.abs(kx);
                        assert l < Om && Om < kx_m && kx_m < ky_m;
                        public static double code(double l, double Om, double kx_m, double ky_m) {
                        	double tmp;
                        	if (((Math.pow(Math.sin(ky_m), 2.0) + Math.pow(Math.sin(kx_m), 2.0)) * Math.pow(((l * 2.0) / Om), 2.0)) <= 3.8) {
                        		tmp = Math.sqrt((0.5 + 0.5));
                        	} else {
                        		tmp = Math.sqrt(0.5);
                        	}
                        	return tmp;
                        }
                        
                        ky_m = math.fabs(ky)
                        kx_m = math.fabs(kx)
                        [l, Om, kx_m, ky_m] = sort([l, Om, kx_m, ky_m])
                        def code(l, Om, kx_m, ky_m):
                        	tmp = 0
                        	if ((math.pow(math.sin(ky_m), 2.0) + math.pow(math.sin(kx_m), 2.0)) * math.pow(((l * 2.0) / Om), 2.0)) <= 3.8:
                        		tmp = math.sqrt((0.5 + 0.5))
                        	else:
                        		tmp = math.sqrt(0.5)
                        	return tmp
                        
                        ky_m = abs(ky)
                        kx_m = abs(kx)
                        l, Om, kx_m, ky_m = sort([l, Om, kx_m, ky_m])
                        function code(l, Om, kx_m, ky_m)
                        	tmp = 0.0
                        	if (Float64(Float64((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (Float64(Float64(l * 2.0) / Om) ^ 2.0)) <= 3.8)
                        		tmp = sqrt(Float64(0.5 + 0.5));
                        	else
                        		tmp = sqrt(0.5);
                        	end
                        	return tmp
                        end
                        
                        ky_m = abs(ky);
                        kx_m = abs(kx);
                        l, Om, kx_m, ky_m = num2cell(sort([l, Om, kx_m, ky_m])){:}
                        function tmp_2 = code(l, Om, kx_m, ky_m)
                        	tmp = 0.0;
                        	if ((((sin(ky_m) ^ 2.0) + (sin(kx_m) ^ 2.0)) * (((l * 2.0) / Om) ^ 2.0)) <= 3.8)
                        		tmp = sqrt((0.5 + 0.5));
                        	else
                        		tmp = sqrt(0.5);
                        	end
                        	tmp_2 = tmp;
                        end
                        
                        ky_m = N[Abs[ky], $MachinePrecision]
                        kx_m = N[Abs[kx], $MachinePrecision]
                        NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                        code[l_, Om_, kx$95$m_, ky$95$m_] := If[LessEqual[N[(N[(N[Power[N[Sin[ky$95$m], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[kx$95$m], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] * N[Power[N[(N[(l * 2.0), $MachinePrecision] / Om), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision], 3.8], N[Sqrt[N[(0.5 + 0.5), $MachinePrecision]], $MachinePrecision], N[Sqrt[0.5], $MachinePrecision]]
                        
                        \begin{array}{l}
                        ky_m = \left|ky\right|
                        \\
                        kx_m = \left|kx\right|
                        \\
                        [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\
                        \\
                        \begin{array}{l}
                        \mathbf{if}\;\left({\sin ky\_m}^{2} + {\sin kx\_m}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 3.8:\\
                        \;\;\;\;\sqrt{0.5 + 0.5}\\
                        
                        \mathbf{else}:\\
                        \;\;\;\;\sqrt{0.5}\\
                        
                        
                        \end{array}
                        \end{array}
                        
                        Derivation
                        1. Split input into 2 regimes
                        2. if (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64)))) < 3.7999999999999998

                          1. Initial program 100.0%

                            \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
                          2. Add Preprocessing
                          3. Taylor expanded in l around 0

                            \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                          4. Step-by-step derivation
                            1. Applied rewrites97.5%

                              \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                            2. Step-by-step derivation
                              1. lift-*.f64N/A

                                \[\leadsto \sqrt{\color{blue}{\frac{1}{2} \cdot \left(1 + \frac{1}{1}\right)}} \]
                              2. lift-/.f64N/A

                                \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
                              3. metadata-evalN/A

                                \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} \cdot \left(1 + \frac{1}{1}\right)} \]
                              4. lift-+.f64N/A

                                \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(1 + \frac{1}{1}\right)}} \]
                              5. +-commutativeN/A

                                \[\leadsto \sqrt{\frac{1}{2} \cdot \color{blue}{\left(\frac{1}{1} + 1\right)}} \]
                              6. distribute-rgt-inN/A

                                \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + 1 \cdot \frac{1}{2}}} \]
                              7. metadata-evalN/A

                                \[\leadsto \sqrt{\frac{1}{1} \cdot \frac{1}{2} + \color{blue}{\frac{1}{2}}} \]
                              8. lower-+.f64N/A

                                \[\leadsto \sqrt{\color{blue}{\frac{1}{1} \cdot \frac{1}{2} + \frac{1}{2}}} \]
                            3. Applied rewrites97.5%

                              \[\leadsto \color{blue}{\sqrt{\frac{0.5}{1} + 0.5}} \]
                            4. Taylor expanded in l around 0

                              \[\leadsto \sqrt{\color{blue}{\frac{1}{2}} + \frac{1}{2}} \]
                            5. Step-by-step derivation
                              1. Applied rewrites97.5%

                                \[\leadsto \sqrt{\color{blue}{0.5} + 0.5} \]

                              if 3.7999999999999998 < (*.f64 (pow.f64 (/.f64 (*.f64 #s(literal 2 binary64) l) Om) #s(literal 2 binary64)) (+.f64 (pow.f64 (sin.f64 kx) #s(literal 2 binary64)) (pow.f64 (sin.f64 ky) #s(literal 2 binary64))))

                              1. Initial program 98.1%

                                \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
                              2. Add Preprocessing
                              3. Taylor expanded in l around 0

                                \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                              4. Step-by-step derivation
                                1. Applied rewrites20.0%

                                  \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                                2. Taylor expanded in l around inf

                                  \[\leadsto \color{blue}{\sqrt{\frac{1}{2}}} \]
                                3. Step-by-step derivation
                                  1. lower-sqrt.f6497.1

                                    \[\leadsto \color{blue}{\sqrt{0.5}} \]
                                4. Applied rewrites97.1%

                                  \[\leadsto \color{blue}{\sqrt{0.5}} \]
                              5. Recombined 2 regimes into one program.
                              6. Final simplification97.4%

                                \[\leadsto \begin{array}{l} \mathbf{if}\;\left({\sin ky}^{2} + {\sin kx}^{2}\right) \cdot {\left(\frac{\ell \cdot 2}{Om}\right)}^{2} \leq 3.8:\\ \;\;\;\;\sqrt{0.5 + 0.5}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{0.5}\\ \end{array} \]
                              7. Add Preprocessing

                              Alternative 6: 55.9% accurate, 52.8× speedup?

                              \[\begin{array}{l} ky_m = \left|ky\right| \\ kx_m = \left|kx\right| \\ [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\ \\ \sqrt{0.5} \end{array} \]
                              ky_m = (fabs.f64 ky)
                              kx_m = (fabs.f64 kx)
                              NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                              (FPCore (l Om kx_m ky_m) :precision binary64 (sqrt 0.5))
                              ky_m = fabs(ky);
                              kx_m = fabs(kx);
                              assert(l < Om && Om < kx_m && kx_m < ky_m);
                              double code(double l, double Om, double kx_m, double ky_m) {
                              	return sqrt(0.5);
                              }
                              
                              ky_m = abs(ky)
                              kx_m = abs(kx)
                              NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                              real(8) function code(l, om, kx_m, ky_m)
                                  real(8), intent (in) :: l
                                  real(8), intent (in) :: om
                                  real(8), intent (in) :: kx_m
                                  real(8), intent (in) :: ky_m
                                  code = sqrt(0.5d0)
                              end function
                              
                              ky_m = Math.abs(ky);
                              kx_m = Math.abs(kx);
                              assert l < Om && Om < kx_m && kx_m < ky_m;
                              public static double code(double l, double Om, double kx_m, double ky_m) {
                              	return Math.sqrt(0.5);
                              }
                              
                              ky_m = math.fabs(ky)
                              kx_m = math.fabs(kx)
                              [l, Om, kx_m, ky_m] = sort([l, Om, kx_m, ky_m])
                              def code(l, Om, kx_m, ky_m):
                              	return math.sqrt(0.5)
                              
                              ky_m = abs(ky)
                              kx_m = abs(kx)
                              l, Om, kx_m, ky_m = sort([l, Om, kx_m, ky_m])
                              function code(l, Om, kx_m, ky_m)
                              	return sqrt(0.5)
                              end
                              
                              ky_m = abs(ky);
                              kx_m = abs(kx);
                              l, Om, kx_m, ky_m = num2cell(sort([l, Om, kx_m, ky_m])){:}
                              function tmp = code(l, Om, kx_m, ky_m)
                              	tmp = sqrt(0.5);
                              end
                              
                              ky_m = N[Abs[ky], $MachinePrecision]
                              kx_m = N[Abs[kx], $MachinePrecision]
                              NOTE: l, Om, kx_m, and ky_m should be sorted in increasing order before calling this function.
                              code[l_, Om_, kx$95$m_, ky$95$m_] := N[Sqrt[0.5], $MachinePrecision]
                              
                              \begin{array}{l}
                              ky_m = \left|ky\right|
                              \\
                              kx_m = \left|kx\right|
                              \\
                              [l, Om, kx_m, ky_m] = \mathsf{sort}([l, Om, kx_m, ky_m])\\
                              \\
                              \sqrt{0.5}
                              \end{array}
                              
                              Derivation
                              1. Initial program 99.2%

                                \[\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {\left(\frac{2 \cdot \ell}{Om}\right)}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)} \]
                              2. Add Preprocessing
                              3. Taylor expanded in l around 0

                                \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                              4. Step-by-step derivation
                                1. Applied rewrites65.1%

                                  \[\leadsto \sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\color{blue}{1}}\right)} \]
                                2. Taylor expanded in l around inf

                                  \[\leadsto \color{blue}{\sqrt{\frac{1}{2}}} \]
                                3. Step-by-step derivation
                                  1. lower-sqrt.f6452.2

                                    \[\leadsto \color{blue}{\sqrt{0.5}} \]
                                4. Applied rewrites52.2%

                                  \[\leadsto \color{blue}{\sqrt{0.5}} \]
                                5. Add Preprocessing

                                Reproduce

                                ?
                                herbie shell --seed 2024312 
                                (FPCore (l Om kx ky)
                                  :name "Toniolo and Linder, Equation (3a)"
                                  :precision binary64
                                  (sqrt (* (/ 1.0 2.0) (+ 1.0 (/ 1.0 (sqrt (+ 1.0 (* (pow (/ (* 2.0 l) Om) 2.0) (+ (pow (sin kx) 2.0) (pow (sin ky) 2.0))))))))))