
(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:
Herbie found 8 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(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}
Om_m = (fabs.f64 Om)
l_m = (fabs.f64 l)
(FPCore (l_m Om_m kx ky)
:precision binary64
(let* ((t_0 (/ (* 2.0 l_m) Om_m)))
(if (<= t_0 2e+129)
(sqrt
(*
(/ 1.0 2.0)
(+
1.0
(/
1.0
(sqrt
(+
1.0
(* (pow t_0 2.0) (+ (pow (sin kx) 2.0) (pow (sin ky) 2.0)))))))))
(sqrt 0.5))))Om_m = fabs(Om);
l_m = fabs(l);
double code(double l_m, double Om_m, double kx, double ky) {
double t_0 = (2.0 * l_m) / Om_m;
double tmp;
if (t_0 <= 2e+129) {
tmp = sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt((1.0 + (pow(t_0, 2.0) * (pow(sin(kx), 2.0) + pow(sin(ky), 2.0)))))))));
} else {
tmp = sqrt(0.5);
}
return tmp;
}
Om_m = abs(om)
l_m = abs(l)
real(8) function code(l_m, om_m, kx, ky)
real(8), intent (in) :: l_m
real(8), intent (in) :: om_m
real(8), intent (in) :: kx
real(8), intent (in) :: ky
real(8) :: t_0
real(8) :: tmp
t_0 = (2.0d0 * l_m) / om_m
if (t_0 <= 2d+129) then
tmp = sqrt(((1.0d0 / 2.0d0) * (1.0d0 + (1.0d0 / sqrt((1.0d0 + ((t_0 ** 2.0d0) * ((sin(kx) ** 2.0d0) + (sin(ky) ** 2.0d0)))))))))
else
tmp = sqrt(0.5d0)
end if
code = tmp
end function
Om_m = Math.abs(Om);
l_m = Math.abs(l);
public static double code(double l_m, double Om_m, double kx, double ky) {
double t_0 = (2.0 * l_m) / Om_m;
double tmp;
if (t_0 <= 2e+129) {
tmp = Math.sqrt(((1.0 / 2.0) * (1.0 + (1.0 / Math.sqrt((1.0 + (Math.pow(t_0, 2.0) * (Math.pow(Math.sin(kx), 2.0) + Math.pow(Math.sin(ky), 2.0)))))))));
} else {
tmp = Math.sqrt(0.5);
}
return tmp;
}
Om_m = math.fabs(Om) l_m = math.fabs(l) def code(l_m, Om_m, kx, ky): t_0 = (2.0 * l_m) / Om_m tmp = 0 if t_0 <= 2e+129: tmp = math.sqrt(((1.0 / 2.0) * (1.0 + (1.0 / math.sqrt((1.0 + (math.pow(t_0, 2.0) * (math.pow(math.sin(kx), 2.0) + math.pow(math.sin(ky), 2.0))))))))) else: tmp = math.sqrt(0.5) return tmp
Om_m = abs(Om) l_m = abs(l) function code(l_m, Om_m, kx, ky) t_0 = Float64(Float64(2.0 * l_m) / Om_m) tmp = 0.0 if (t_0 <= 2e+129) tmp = sqrt(Float64(Float64(1.0 / 2.0) * Float64(1.0 + Float64(1.0 / sqrt(Float64(1.0 + Float64((t_0 ^ 2.0) * Float64((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0))))))))); else tmp = sqrt(0.5); end return tmp end
Om_m = abs(Om); l_m = abs(l); function tmp_2 = code(l_m, Om_m, kx, ky) t_0 = (2.0 * l_m) / Om_m; tmp = 0.0; if (t_0 <= 2e+129) tmp = sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt((1.0 + ((t_0 ^ 2.0) * ((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0))))))))); else tmp = sqrt(0.5); end tmp_2 = tmp; end
Om_m = N[Abs[Om], $MachinePrecision]
l_m = N[Abs[l], $MachinePrecision]
code[l$95$m_, Om$95$m_, kx_, ky_] := Block[{t$95$0 = N[(N[(2.0 * l$95$m), $MachinePrecision] / Om$95$m), $MachinePrecision]}, If[LessEqual[t$95$0, 2e+129], N[Sqrt[N[(N[(1.0 / 2.0), $MachinePrecision] * N[(1.0 + N[(1.0 / N[Sqrt[N[(1.0 + N[(N[Power[t$95$0, 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], N[Sqrt[0.5], $MachinePrecision]]]
\begin{array}{l}
Om_m = \left|Om\right|
\\
l_m = \left|\ell\right|
\\
\begin{array}{l}
t_0 := \frac{2 \cdot l\_m}{Om\_m}\\
\mathbf{if}\;t\_0 \leq 2 \cdot 10^{+129}:\\
\;\;\;\;\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{1 + {t\_0}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right)}}\right)}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{0.5}\\
\end{array}
\end{array}
if (/.f64 (*.f64 #s(literal 2 binary64) l) Om) < 2e129Initial program 100.0%
if 2e129 < (/.f64 (*.f64 #s(literal 2 binary64) l) Om) Initial program 93.9%
Taylor expanded in l around inf
Applied rewrites99.2%
Om_m = (fabs.f64 Om)
l_m = (fabs.f64 l)
(FPCore (l_m Om_m kx ky)
:precision binary64
(let* ((t_0 (/ (* 2.0 l_m) Om_m)))
(if (<=
(* (pow t_0 2.0) (+ (pow (sin kx) 2.0) (pow (sin ky) 2.0)))
2000000.0)
(sqrt
(*
(/ 1.0 2.0)
(+
1.0
(/
1.0
(sqrt
(fma
(* (/ l_m Om_m) 4.0)
(*
(/ l_m Om_m)
(+
(+ 0.5 (* -0.5 (cos (+ kx kx))))
(+ 0.5 (* -0.5 (cos (+ ky ky))))))
1.0))))))
(sqrt
(* (/ 1.0 2.0) (+ 1.0 (/ 1.0 (* t_0 (hypot (sin kx) (sin ky))))))))))Om_m = fabs(Om);
l_m = fabs(l);
double code(double l_m, double Om_m, double kx, double ky) {
double t_0 = (2.0 * l_m) / Om_m;
double tmp;
if ((pow(t_0, 2.0) * (pow(sin(kx), 2.0) + pow(sin(ky), 2.0))) <= 2000000.0) {
tmp = sqrt(((1.0 / 2.0) * (1.0 + (1.0 / sqrt(fma(((l_m / Om_m) * 4.0), ((l_m / Om_m) * ((0.5 + (-0.5 * cos((kx + kx)))) + (0.5 + (-0.5 * cos((ky + ky)))))), 1.0))))));
} else {
tmp = sqrt(((1.0 / 2.0) * (1.0 + (1.0 / (t_0 * hypot(sin(kx), sin(ky)))))));
}
return tmp;
}
Om_m = abs(Om) l_m = abs(l) function code(l_m, Om_m, kx, ky) t_0 = Float64(Float64(2.0 * l_m) / Om_m) tmp = 0.0 if (Float64((t_0 ^ 2.0) * Float64((sin(kx) ^ 2.0) + (sin(ky) ^ 2.0))) <= 2000000.0) tmp = sqrt(Float64(Float64(1.0 / 2.0) * Float64(1.0 + Float64(1.0 / sqrt(fma(Float64(Float64(l_m / Om_m) * 4.0), Float64(Float64(l_m / Om_m) * Float64(Float64(0.5 + Float64(-0.5 * cos(Float64(kx + kx)))) + Float64(0.5 + Float64(-0.5 * cos(Float64(ky + ky)))))), 1.0)))))); else tmp = sqrt(Float64(Float64(1.0 / 2.0) * Float64(1.0 + Float64(1.0 / Float64(t_0 * hypot(sin(kx), sin(ky))))))); end return tmp end
Om_m = N[Abs[Om], $MachinePrecision]
l_m = N[Abs[l], $MachinePrecision]
code[l$95$m_, Om$95$m_, kx_, ky_] := Block[{t$95$0 = N[(N[(2.0 * l$95$m), $MachinePrecision] / Om$95$m), $MachinePrecision]}, If[LessEqual[N[(N[Power[t$95$0, 2.0], $MachinePrecision] * N[(N[Power[N[Sin[kx], $MachinePrecision], 2.0], $MachinePrecision] + N[Power[N[Sin[ky], $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2000000.0], N[Sqrt[N[(N[(1.0 / 2.0), $MachinePrecision] * N[(1.0 + N[(1.0 / N[Sqrt[N[(N[(N[(l$95$m / Om$95$m), $MachinePrecision] * 4.0), $MachinePrecision] * N[(N[(l$95$m / Om$95$m), $MachinePrecision] * N[(N[(0.5 + N[(-0.5 * N[Cos[N[(kx + kx), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + N[(0.5 + N[(-0.5 * N[Cos[N[(ky + ky), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], N[Sqrt[N[(N[(1.0 / 2.0), $MachinePrecision] * N[(1.0 + N[(1.0 / N[(t$95$0 * N[Sqrt[N[Sin[kx], $MachinePrecision] ^ 2 + N[Sin[ky], $MachinePrecision] ^ 2], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]]
\begin{array}{l}
Om_m = \left|Om\right|
\\
l_m = \left|\ell\right|
\\
\begin{array}{l}
t_0 := \frac{2 \cdot l\_m}{Om\_m}\\
\mathbf{if}\;{t\_0}^{2} \cdot \left({\sin kx}^{2} + {\sin ky}^{2}\right) \leq 2000000:\\
\;\;\;\;\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{\sqrt{\mathsf{fma}\left(\frac{l\_m}{Om\_m} \cdot 4, \frac{l\_m}{Om\_m} \cdot \left(\left(0.5 + -0.5 \cdot \cos \left(kx + kx\right)\right) + \left(0.5 + -0.5 \cdot \cos \left(ky + ky\right)\right)\right), 1\right)}}\right)}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{\frac{1}{2} \cdot \left(1 + \frac{1}{t\_0 \cdot \mathsf{hypot}\left(\sin kx, \sin ky\right)}\right)}\\
\end{array}
\end{array}
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)))) < 2e6Initial program 100.0%
Applied rewrites99.7%
if 2e6 < (*.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)))) Initial program 96.5%
Taylor expanded in l around inf
associate-*r*N/A
lower-*.f64N/A
associate-*r/N/A
lower-/.f64N/A
lower-*.f64N/A
unpow2N/A
unpow2N/A
lower-hypot.f64N/A
lower-sin.f64N/A
lower-sin.f6499.3
Applied rewrites99.3%
herbie shell --seed 2024223
(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))))))))))