
(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))));
}
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(w0, m, d, h, l, d_1)
use fmin_fmax_functions
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:
Herbie found 14 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(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))));
}
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(w0, m, d, h, l, d_1)
use fmin_fmax_functions
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}
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<=
(* w0_m (sqrt (- 1.0 (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)))))
1e+187)
(* w0_m (sqrt (- 1.0 (* (pow (/ (* (/ M_m d) D_m) 2.0) 2.0) (/ h l)))))
(*
w0_m
(sqrt
(-
1.0
(*
(/ D_m 2.0)
(* (/ M_m d) (/ (* (* h D_m) M_m) (* l (* d 2.0)))))))))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
double tmp;
if ((w0_m * sqrt((1.0 - (pow(((M_m * D_m) / (2.0 * d)), 2.0) * (h / l))))) <= 1e+187) {
tmp = w0_m * sqrt((1.0 - (pow((((M_m / d) * D_m) / 2.0), 2.0) * (h / l))));
} else {
tmp = w0_m * sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0)))))));
}
return w0_s * tmp;
}
D_m = private
M_m = private
w0\_m = private
w0\_s = private
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(w0_s, w0_m, m_m, d_m, h, l, d)
use fmin_fmax_functions
real(8), intent (in) :: w0_s
real(8), intent (in) :: w0_m
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
real(8) :: tmp
if ((w0_m * sqrt((1.0d0 - ((((m_m * d_m) / (2.0d0 * d)) ** 2.0d0) * (h / l))))) <= 1d+187) then
tmp = w0_m * sqrt((1.0d0 - (((((m_m / d) * d_m) / 2.0d0) ** 2.0d0) * (h / l))))
else
tmp = w0_m * sqrt((1.0d0 - ((d_m / 2.0d0) * ((m_m / d) * (((h * d_m) * m_m) / (l * (d * 2.0d0)))))))
end if
code = w0_s * tmp
end function
D_m = Math.abs(D);
M_m = Math.abs(M);
w0\_m = Math.abs(w0);
w0\_s = Math.copySign(1.0, w0);
assert w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d;
public static double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
double tmp;
if ((w0_m * Math.sqrt((1.0 - (Math.pow(((M_m * D_m) / (2.0 * d)), 2.0) * (h / l))))) <= 1e+187) {
tmp = w0_m * Math.sqrt((1.0 - (Math.pow((((M_m / d) * D_m) / 2.0), 2.0) * (h / l))));
} else {
tmp = w0_m * Math.sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0)))))));
}
return w0_s * tmp;
}
D_m = math.fabs(D) M_m = math.fabs(M) w0\_m = math.fabs(w0) w0\_s = math.copysign(1.0, w0) [w0_m, M_m, D_m, h, l, d] = sort([w0_m, M_m, D_m, h, l, d]) def code(w0_s, w0_m, M_m, D_m, h, l, d): tmp = 0 if (w0_m * math.sqrt((1.0 - (math.pow(((M_m * D_m) / (2.0 * d)), 2.0) * (h / l))))) <= 1e+187: tmp = w0_m * math.sqrt((1.0 - (math.pow((((M_m / d) * D_m) / 2.0), 2.0) * (h / l)))) else: tmp = w0_m * math.sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0))))))) return w0_s * tmp
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, M_m, D_m, h, l, d) tmp = 0.0 if (Float64(w0_m * sqrt(Float64(1.0 - Float64((Float64(Float64(M_m * D_m) / Float64(2.0 * d)) ^ 2.0) * Float64(h / l))))) <= 1e+187) tmp = Float64(w0_m * sqrt(Float64(1.0 - Float64((Float64(Float64(Float64(M_m / d) * D_m) / 2.0) ^ 2.0) * Float64(h / l))))); else tmp = Float64(w0_m * sqrt(Float64(1.0 - Float64(Float64(D_m / 2.0) * Float64(Float64(M_m / d) * Float64(Float64(Float64(h * D_m) * M_m) / Float64(l * Float64(d * 2.0)))))))); end return Float64(w0_s * tmp) end
D_m = abs(D);
M_m = abs(M);
w0\_m = abs(w0);
w0\_s = sign(w0) * abs(1.0);
w0_m, M_m, D_m, h, l, d = num2cell(sort([w0_m, M_m, D_m, h, l, d])){:}
function tmp_2 = code(w0_s, w0_m, M_m, D_m, h, l, d)
tmp = 0.0;
if ((w0_m * sqrt((1.0 - ((((M_m * D_m) / (2.0 * d)) ^ 2.0) * (h / l))))) <= 1e+187)
tmp = w0_m * sqrt((1.0 - (((((M_m / d) * D_m) / 2.0) ^ 2.0) * (h / l))));
else
tmp = w0_m * sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0)))))));
end
tmp_2 = w0_s * tmp;
end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * If[LessEqual[N[(w0$95$m * 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]), $MachinePrecision], 1e+187], N[(w0$95$m * N[Sqrt[N[(1.0 - N[(N[Power[N[(N[(N[(M$95$m / d), $MachinePrecision] * D$95$m), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * N[Sqrt[N[(1.0 - N[(N[(D$95$m / 2.0), $MachinePrecision] * N[(N[(M$95$m / d), $MachinePrecision] * N[(N[(N[(h * D$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] / N[(l * N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;w0\_m \cdot \sqrt{1 - {\left(\frac{M\_m \cdot D\_m}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}} \leq 10^{+187}:\\
\;\;\;\;w0\_m \cdot \sqrt{1 - {\left(\frac{\frac{M\_m}{d} \cdot D\_m}{2}\right)}^{2} \cdot \frac{h}{\ell}}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot \sqrt{1 - \frac{D\_m}{2} \cdot \left(\frac{M\_m}{d} \cdot \frac{\left(h \cdot D\_m\right) \cdot M\_m}{\ell \cdot \left(d \cdot 2\right)}\right)}\\
\end{array}
\end{array}
if (*.f64 w0 (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))))) < 9.99999999999999907e186Initial program 93.0%
lift-/.f64N/A
lift-*.f64N/A
lift-*.f64N/A
*-commutativeN/A
times-fracN/A
associate-*r/N/A
lower-/.f64N/A
lower-*.f64N/A
lower-/.f6492.4
Applied rewrites92.4%
if 9.99999999999999907e186 < (*.f64 w0 (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))))) Initial program 61.7%
lift-*.f64N/A
lift-pow.f64N/A
unpow2N/A
associate-*l*N/A
lift-/.f64N/A
lift-*.f64N/A
*-commutativeN/A
lift-*.f64N/A
times-fracN/A
associate-*l*N/A
lower-*.f64N/A
lower-/.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower-/.f64N/A
lift-/.f64N/A
lift-*.f64N/A
*-commutativeN/A
associate-/l*N/A
associate-*r*N/A
lower-*.f64N/A
Applied rewrites74.4%
lift-*.f64N/A
lift-*.f64N/A
lift-/.f64N/A
associate-*l/N/A
lift-/.f64N/A
lift-/.f64N/A
associate-/l/N/A
*-commutativeN/A
frac-timesN/A
lower-/.f64N/A
lower-*.f64N/A
lower-*.f64N/A
lower-*.f64N/A
*-commutativeN/A
lower-*.f6483.6
Applied rewrites83.6%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(let* ((t_0
(*
w0_m
(sqrt (- 1.0 (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)))))))
(*
w0_s
(if (<= t_0 1e+187)
t_0
(*
w0_m
(sqrt
(-
1.0
(*
(/ D_m 2.0)
(* (/ M_m d) (/ (* (* h D_m) M_m) (* l (* d 2.0))))))))))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
double t_0 = w0_m * sqrt((1.0 - (pow(((M_m * D_m) / (2.0 * d)), 2.0) * (h / l))));
double tmp;
if (t_0 <= 1e+187) {
tmp = t_0;
} else {
tmp = w0_m * sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0)))))));
}
return w0_s * tmp;
}
D_m = private
M_m = private
w0\_m = private
w0\_s = private
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(w0_s, w0_m, m_m, d_m, h, l, d)
use fmin_fmax_functions
real(8), intent (in) :: w0_s
real(8), intent (in) :: w0_m
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
real(8) :: t_0
real(8) :: tmp
t_0 = w0_m * sqrt((1.0d0 - ((((m_m * d_m) / (2.0d0 * d)) ** 2.0d0) * (h / l))))
if (t_0 <= 1d+187) then
tmp = t_0
else
tmp = w0_m * sqrt((1.0d0 - ((d_m / 2.0d0) * ((m_m / d) * (((h * d_m) * m_m) / (l * (d * 2.0d0)))))))
end if
code = w0_s * tmp
end function
D_m = Math.abs(D);
M_m = Math.abs(M);
w0\_m = Math.abs(w0);
w0\_s = Math.copySign(1.0, w0);
assert w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d;
public static double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
double t_0 = w0_m * Math.sqrt((1.0 - (Math.pow(((M_m * D_m) / (2.0 * d)), 2.0) * (h / l))));
double tmp;
if (t_0 <= 1e+187) {
tmp = t_0;
} else {
tmp = w0_m * Math.sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0)))))));
}
return w0_s * tmp;
}
D_m = math.fabs(D) M_m = math.fabs(M) w0\_m = math.fabs(w0) w0\_s = math.copysign(1.0, w0) [w0_m, M_m, D_m, h, l, d] = sort([w0_m, M_m, D_m, h, l, d]) def code(w0_s, w0_m, M_m, D_m, h, l, d): t_0 = w0_m * math.sqrt((1.0 - (math.pow(((M_m * D_m) / (2.0 * d)), 2.0) * (h / l)))) tmp = 0 if t_0 <= 1e+187: tmp = t_0 else: tmp = w0_m * math.sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0))))))) return w0_s * tmp
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, M_m, D_m, h, l, d) t_0 = Float64(w0_m * sqrt(Float64(1.0 - Float64((Float64(Float64(M_m * D_m) / Float64(2.0 * d)) ^ 2.0) * Float64(h / l))))) tmp = 0.0 if (t_0 <= 1e+187) tmp = t_0; else tmp = Float64(w0_m * sqrt(Float64(1.0 - Float64(Float64(D_m / 2.0) * Float64(Float64(M_m / d) * Float64(Float64(Float64(h * D_m) * M_m) / Float64(l * Float64(d * 2.0)))))))); end return Float64(w0_s * tmp) end
D_m = abs(D);
M_m = abs(M);
w0\_m = abs(w0);
w0\_s = sign(w0) * abs(1.0);
w0_m, M_m, D_m, h, l, d = num2cell(sort([w0_m, M_m, D_m, h, l, d])){:}
function tmp_2 = code(w0_s, w0_m, M_m, D_m, h, l, d)
t_0 = w0_m * sqrt((1.0 - ((((M_m * D_m) / (2.0 * d)) ^ 2.0) * (h / l))));
tmp = 0.0;
if (t_0 <= 1e+187)
tmp = t_0;
else
tmp = w0_m * sqrt((1.0 - ((D_m / 2.0) * ((M_m / d) * (((h * D_m) * M_m) / (l * (d * 2.0)))))));
end
tmp_2 = w0_s * tmp;
end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := Block[{t$95$0 = N[(w0$95$m * 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]), $MachinePrecision]}, N[(w0$95$s * If[LessEqual[t$95$0, 1e+187], t$95$0, N[(w0$95$m * N[Sqrt[N[(1.0 - N[(N[(D$95$m / 2.0), $MachinePrecision] * N[(N[(M$95$m / d), $MachinePrecision] * N[(N[(N[(h * D$95$m), $MachinePrecision] * M$95$m), $MachinePrecision] / N[(l * N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
\begin{array}{l}
t_0 := w0\_m \cdot \sqrt{1 - {\left(\frac{M\_m \cdot D\_m}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell}}\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;t\_0 \leq 10^{+187}:\\
\;\;\;\;t\_0\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot \sqrt{1 - \frac{D\_m}{2} \cdot \left(\frac{M\_m}{d} \cdot \frac{\left(h \cdot D\_m\right) \cdot M\_m}{\ell \cdot \left(d \cdot 2\right)}\right)}\\
\end{array}
\end{array}
\end{array}
if (*.f64 w0 (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))))) < 9.99999999999999907e186Initial program 93.0%
if 9.99999999999999907e186 < (*.f64 w0 (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))))) Initial program 61.7%
lift-*.f64N/A
lift-pow.f64N/A
unpow2N/A
associate-*l*N/A
lift-/.f64N/A
lift-*.f64N/A
*-commutativeN/A
lift-*.f64N/A
times-fracN/A
associate-*l*N/A
lower-*.f64N/A
lower-/.f64N/A
*-commutativeN/A
lower-*.f64N/A
lower-/.f64N/A
lift-/.f64N/A
lift-*.f64N/A
*-commutativeN/A
associate-/l*N/A
associate-*r*N/A
lower-*.f64N/A
Applied rewrites74.4%
lift-*.f64N/A
lift-*.f64N/A
lift-/.f64N/A
associate-*l/N/A
lift-/.f64N/A
lift-/.f64N/A
associate-/l/N/A
*-commutativeN/A
frac-timesN/A
lower-/.f64N/A
lower-*.f64N/A
lower-*.f64N/A
lower-*.f64N/A
*-commutativeN/A
lower-*.f6483.6
Applied rewrites83.6%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (sqrt (- 1.0 (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)))) 1.0)
(* w0_m 1.0)
(*
w0_m
(sqrt
(fma (* (* -0.25 h) (* (* (/ M_m d) D_m) (/ M_m d))) (/ D_m l) 1.0))))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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_m * 1.0;
} else {
tmp = w0_m * sqrt(fma(((-0.25 * h) * (((M_m / d) * D_m) * (M_m / d))), (D_m / l), 1.0));
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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_m * 1.0); else tmp = Float64(w0_m * sqrt(fma(Float64(Float64(-0.25 * h) * Float64(Float64(Float64(M_m / d) * D_m) * Float64(M_m / d))), Float64(D_m / l), 1.0))); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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$95$m * 1.0), $MachinePrecision], N[(w0$95$m * N[Sqrt[N[(N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(M$95$m / d), $MachinePrecision] * D$95$m), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(D$95$m / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \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\_m \cdot 1\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(\left(-0.25 \cdot h\right) \cdot \left(\left(\frac{M\_m}{d} \cdot D\_m\right) \cdot \frac{M\_m}{d}\right), \frac{D\_m}{\ell}, 1\right)}\\
\end{array}
\end{array}
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)))) < 1Initial program 99.0%
Taylor expanded in M around 0
Applied rewrites98.6%
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)))) Initial program 59.8%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites45.3%
Applied rewrites49.9%
Applied rewrites69.1%
Applied rewrites69.0%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -1e-9)
(*
w0_m
(sqrt
(fma (* h -0.25) (* (/ (* D_m M_m) (* l d)) (* (/ M_m d) D_m)) 1.0)))
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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)) <= -1e-9) {
tmp = w0_m * sqrt(fma((h * -0.25), (((D_m * M_m) / (l * d)) * ((M_m / d) * D_m)), 1.0));
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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)) <= -1e-9) tmp = Float64(w0_m * sqrt(fma(Float64(h * -0.25), Float64(Float64(Float64(D_m * M_m) / Float64(l * d)) * Float64(Float64(M_m / d) * D_m)), 1.0))); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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], -1e-9], N[(w0$95$m * N[Sqrt[N[(N[(h * -0.25), $MachinePrecision] * N[(N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(l * d), $MachinePrecision]), $MachinePrecision] * N[(N[(M$95$m / d), $MachinePrecision] * D$95$m), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;{\left(\frac{M\_m \cdot D\_m}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell} \leq -1 \cdot 10^{-9}:\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{D\_m \cdot M\_m}{\ell \cdot d} \cdot \left(\frac{M\_m}{d} \cdot D\_m\right), 1\right)}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -1.00000000000000006e-9Initial program 70.8%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites44.4%
Applied rewrites48.5%
Applied rewrites67.3%
Applied rewrites66.2%
if -1.00000000000000006e-9 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.3%
Taylor expanded in M around 0
Applied rewrites97.1%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -1e-9)
(*
w0_m
(sqrt
(fma (* h -0.25) (* (* (* (/ M_m d) M_m) D_m) (/ D_m (* l d))) 1.0)))
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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)) <= -1e-9) {
tmp = w0_m * sqrt(fma((h * -0.25), ((((M_m / d) * M_m) * D_m) * (D_m / (l * d))), 1.0));
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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)) <= -1e-9) tmp = Float64(w0_m * sqrt(fma(Float64(h * -0.25), Float64(Float64(Float64(Float64(M_m / d) * M_m) * D_m) * Float64(D_m / Float64(l * d))), 1.0))); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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], -1e-9], N[(w0$95$m * N[Sqrt[N[(N[(h * -0.25), $MachinePrecision] * N[(N[(N[(N[(M$95$m / d), $MachinePrecision] * M$95$m), $MachinePrecision] * D$95$m), $MachinePrecision] * N[(D$95$m / N[(l * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;{\left(\frac{M\_m \cdot D\_m}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell} \leq -1 \cdot 10^{-9}:\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \left(\left(\frac{M\_m}{d} \cdot M\_m\right) \cdot D\_m\right) \cdot \frac{D\_m}{\ell \cdot d}, 1\right)}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -1.00000000000000006e-9Initial program 70.8%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites44.4%
Applied rewrites54.8%
if -1.00000000000000006e-9 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.3%
Taylor expanded in M around 0
Applied rewrites97.1%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -500.0)
(*
w0_m
(sqrt
(fma (* h -0.25) (* M_m (* (* D_m M_m) (/ D_m (* (* d l) d)))) 1.0)))
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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)) <= -500.0) {
tmp = w0_m * sqrt(fma((h * -0.25), (M_m * ((D_m * M_m) * (D_m / ((d * l) * d)))), 1.0));
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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)) <= -500.0) tmp = Float64(w0_m * sqrt(fma(Float64(h * -0.25), Float64(M_m * Float64(Float64(D_m * M_m) * Float64(D_m / Float64(Float64(d * l) * d)))), 1.0))); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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], -500.0], N[(w0$95$m * 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[(d * l), $MachinePrecision] * d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;{\left(\frac{M\_m \cdot D\_m}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell} \leq -500:\\
\;\;\;\;w0\_m \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(d \cdot \ell\right) \cdot d}\right), 1\right)}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -500Initial program 70.1%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites45.5%
Applied rewrites48.5%
Applied rewrites54.9%
if -500 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.3%
Taylor expanded in M around 0
Applied rewrites96.7%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -500.0)
(*
w0_m
(sqrt
(fma (* h -0.25) (/ (* (* M_m D_m) (* M_m D_m)) (* (* d d) l)) 1.0)))
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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)) <= -500.0) {
tmp = w0_m * sqrt(fma((h * -0.25), (((M_m * D_m) * (M_m * D_m)) / ((d * d) * l)), 1.0));
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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)) <= -500.0) tmp = Float64(w0_m * sqrt(fma(Float64(h * -0.25), Float64(Float64(Float64(M_m * D_m) * Float64(M_m * D_m)) / Float64(Float64(d * d) * l)), 1.0))); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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], -500.0], N[(w0$95$m * N[Sqrt[N[(N[(h * -0.25), $MachinePrecision] * N[(N[(N[(M$95$m * D$95$m), $MachinePrecision] * N[(M$95$m * D$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(d * d), $MachinePrecision] * l), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;{\left(\frac{M\_m \cdot D\_m}{2 \cdot d}\right)}^{2} \cdot \frac{h}{\ell} \leq -500:\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{\left(M\_m \cdot D\_m\right) \cdot \left(M\_m \cdot D\_m\right)}{\left(d \cdot d\right) \cdot \ell}, 1\right)}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -500Initial program 70.1%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites45.5%
Taylor expanded in M around 0
Applied rewrites49.4%
if -500 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.3%
Taylor expanded in M around 0
Applied rewrites96.7%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -4e+125)
(fma
(* (* D_m D_m) -0.125)
(* (/ M_m (* l d)) (/ (* (* h M_m) w0_m) d))
w0_m)
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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+125) {
tmp = fma(((D_m * D_m) * -0.125), ((M_m / (l * d)) * (((h * M_m) * w0_m) / d)), w0_m);
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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+125) tmp = fma(Float64(Float64(D_m * D_m) * -0.125), Float64(Float64(M_m / Float64(l * d)) * Float64(Float64(Float64(h * M_m) * w0_m) / d)), w0_m); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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+125], N[(N[(N[(D$95$m * D$95$m), $MachinePrecision] * -0.125), $MachinePrecision] * N[(N[(M$95$m / N[(l * d), $MachinePrecision]), $MachinePrecision] * N[(N[(N[(h * M$95$m), $MachinePrecision] * w0$95$m), $MachinePrecision] / d), $MachinePrecision]), $MachinePrecision] + w0$95$m), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \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^{+125}:\\
\;\;\;\;\mathsf{fma}\left(\left(D\_m \cdot D\_m\right) \cdot -0.125, \frac{M\_m}{\ell \cdot d} \cdot \frac{\left(h \cdot M\_m\right) \cdot w0\_m}{d}, w0\_m\right)\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -3.9999999999999997e125Initial program 66.0%
Taylor expanded in M around 0
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
associate-*r*N/A
*-commutativeN/A
associate-*r*N/A
lower-fma.f64N/A
Applied rewrites38.4%
Applied rewrites41.4%
Applied rewrites45.0%
if -3.9999999999999997e125 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.7%
Taylor expanded in M around 0
Applied rewrites92.6%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -4e+153)
(fma
(* (* D_m D_m) -0.125)
(* (* (* M_m (/ h d)) M_m) (/ w0_m (* d l)))
w0_m)
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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+153) {
tmp = fma(((D_m * D_m) * -0.125), (((M_m * (h / d)) * M_m) * (w0_m / (d * l))), w0_m);
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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+153) tmp = fma(Float64(Float64(D_m * D_m) * -0.125), Float64(Float64(Float64(M_m * Float64(h / d)) * M_m) * Float64(w0_m / Float64(d * l))), w0_m); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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+153], N[(N[(N[(D$95$m * D$95$m), $MachinePrecision] * -0.125), $MachinePrecision] * N[(N[(N[(M$95$m * N[(h / d), $MachinePrecision]), $MachinePrecision] * M$95$m), $MachinePrecision] * N[(w0$95$m / N[(d * l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + w0$95$m), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \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^{+153}:\\
\;\;\;\;\mathsf{fma}\left(\left(D\_m \cdot D\_m\right) \cdot -0.125, \left(\left(M\_m \cdot \frac{h}{d}\right) \cdot M\_m\right) \cdot \frac{w0\_m}{d \cdot \ell}, w0\_m\right)\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -4e153Initial program 65.5%
Taylor expanded in M around 0
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
associate-*r*N/A
*-commutativeN/A
associate-*r*N/A
lower-fma.f64N/A
Applied rewrites39.0%
Applied rewrites45.5%
if -4e153 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.7%
Taylor expanded in M around 0
Applied rewrites92.2%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -4e+153)
(fma
(* (* D_m D_m) -0.125)
(* (* (/ w0_m (* (* l d) d)) M_m) (* h M_m))
w0_m)
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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+153) {
tmp = fma(((D_m * D_m) * -0.125), (((w0_m / ((l * d) * d)) * M_m) * (h * M_m)), w0_m);
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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+153) tmp = fma(Float64(Float64(D_m * D_m) * -0.125), Float64(Float64(Float64(w0_m / Float64(Float64(l * d) * d)) * M_m) * Float64(h * M_m)), w0_m); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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+153], N[(N[(N[(D$95$m * D$95$m), $MachinePrecision] * -0.125), $MachinePrecision] * N[(N[(N[(w0$95$m / N[(N[(l * d), $MachinePrecision] * d), $MachinePrecision]), $MachinePrecision] * M$95$m), $MachinePrecision] * N[(h * M$95$m), $MachinePrecision]), $MachinePrecision] + w0$95$m), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \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^{+153}:\\
\;\;\;\;\mathsf{fma}\left(\left(D\_m \cdot D\_m\right) \cdot -0.125, \left(\frac{w0\_m}{\left(\ell \cdot d\right) \cdot d} \cdot M\_m\right) \cdot \left(h \cdot M\_m\right), w0\_m\right)\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -4e153Initial program 65.5%
Taylor expanded in M around 0
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
associate-*r*N/A
*-commutativeN/A
associate-*r*N/A
lower-fma.f64N/A
Applied rewrites39.0%
Applied rewrites42.0%
Applied rewrites43.8%
if -4e153 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.7%
Taylor expanded in M around 0
Applied rewrites92.2%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= (* (pow (/ (* M_m D_m) (* 2.0 d)) 2.0) (/ h l)) -4e-8)
(fma
(* (* D_m D_m) -0.125)
(* M_m (* M_m (* (/ w0_m (* (* l d) d)) h)))
w0_m)
(* w0_m 1.0))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, 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-8) {
tmp = fma(((D_m * D_m) * -0.125), (M_m * (M_m * ((w0_m / ((l * d) * d)) * h))), w0_m);
} else {
tmp = w0_m * 1.0;
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, 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-8) tmp = fma(Float64(Float64(D_m * D_m) * -0.125), Float64(M_m * Float64(M_m * Float64(Float64(w0_m / Float64(Float64(l * d) * d)) * h))), w0_m); else tmp = Float64(w0_m * 1.0); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * 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-8], N[(N[(N[(D$95$m * D$95$m), $MachinePrecision] * -0.125), $MachinePrecision] * N[(M$95$m * N[(M$95$m * N[(N[(w0$95$m / N[(N[(l * d), $MachinePrecision] * d), $MachinePrecision]), $MachinePrecision] * h), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + w0$95$m), $MachinePrecision], N[(w0$95$m * 1.0), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \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^{-8}:\\
\;\;\;\;\mathsf{fma}\left(\left(D\_m \cdot D\_m\right) \cdot -0.125, M\_m \cdot \left(M\_m \cdot \left(\frac{w0\_m}{\left(\ell \cdot d\right) \cdot d} \cdot h\right)\right), w0\_m\right)\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot 1\\
\end{array}
\end{array}
if (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) < -4.0000000000000001e-8Initial program 70.4%
Taylor expanded in M around 0
+-commutativeN/A
*-commutativeN/A
associate-/l*N/A
associate-*r*N/A
*-commutativeN/A
associate-*r*N/A
lower-fma.f64N/A
Applied rewrites33.7%
Applied rewrites37.6%
Applied rewrites39.7%
if -4.0000000000000001e-8 < (*.f64 (pow.f64 (/.f64 (*.f64 M D) (*.f64 #s(literal 2 binary64) d)) #s(literal 2 binary64)) (/.f64 h l)) Initial program 91.3%
Taylor expanded in M around 0
Applied rewrites96.9%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (or (<= (* M_m D_m) 5e-204) (not (<= (* M_m D_m) 5e+225)))
(*
w0_m
(sqrt
(fma (* h -0.25) (* (* (/ D_m d) (/ M_m d)) (* M_m (/ D_m l))) 1.0)))
(*
w0_m
(sqrt
(fma (* h -0.25) (* (/ (* D_m M_m) (* l d)) (* (/ M_m d) D_m)) 1.0))))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
double tmp;
if (((M_m * D_m) <= 5e-204) || !((M_m * D_m) <= 5e+225)) {
tmp = w0_m * sqrt(fma((h * -0.25), (((D_m / d) * (M_m / d)) * (M_m * (D_m / l))), 1.0));
} else {
tmp = w0_m * sqrt(fma((h * -0.25), (((D_m * M_m) / (l * d)) * ((M_m / d) * D_m)), 1.0));
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, M_m, D_m, h, l, d) tmp = 0.0 if ((Float64(M_m * D_m) <= 5e-204) || !(Float64(M_m * D_m) <= 5e+225)) tmp = Float64(w0_m * sqrt(fma(Float64(h * -0.25), Float64(Float64(Float64(D_m / d) * Float64(M_m / d)) * Float64(M_m * Float64(D_m / l))), 1.0))); else tmp = Float64(w0_m * sqrt(fma(Float64(h * -0.25), Float64(Float64(Float64(D_m * M_m) / Float64(l * d)) * Float64(Float64(M_m / d) * D_m)), 1.0))); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * If[Or[LessEqual[N[(M$95$m * D$95$m), $MachinePrecision], 5e-204], N[Not[LessEqual[N[(M$95$m * D$95$m), $MachinePrecision], 5e+225]], $MachinePrecision]], N[(w0$95$m * N[Sqrt[N[(N[(h * -0.25), $MachinePrecision] * N[(N[(N[(D$95$m / d), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision] * N[(M$95$m * N[(D$95$m / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * N[Sqrt[N[(N[(h * -0.25), $MachinePrecision] * N[(N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(l * d), $MachinePrecision]), $MachinePrecision] * N[(N[(M$95$m / d), $MachinePrecision] * D$95$m), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;M\_m \cdot D\_m \leq 5 \cdot 10^{-204} \lor \neg \left(M\_m \cdot D\_m \leq 5 \cdot 10^{+225}\right):\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \left(\frac{D\_m}{d} \cdot \frac{M\_m}{d}\right) \cdot \left(M\_m \cdot \frac{D\_m}{\ell}\right), 1\right)}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(h \cdot -0.25, \frac{D\_m \cdot M\_m}{\ell \cdot d} \cdot \left(\frac{M\_m}{d} \cdot D\_m\right), 1\right)}\\
\end{array}
\end{array}
if (*.f64 M D) < 5.0000000000000002e-204 or 4.99999999999999981e225 < (*.f64 M D) Initial program 82.8%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites66.6%
Applied rewrites68.9%
Applied rewrites85.7%
if 5.0000000000000002e-204 < (*.f64 M D) < 4.99999999999999981e225Initial program 90.8%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites63.1%
Applied rewrites65.9%
Applied rewrites84.4%
Applied rewrites95.1%
Final simplification88.5%
D_m = (fabs.f64 D)
M_m = (fabs.f64 M)
w0\_m = (fabs.f64 w0)
w0\_s = (copysign.f64 #s(literal 1 binary64) w0)
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
(FPCore (w0_s w0_m M_m D_m h l d)
:precision binary64
(*
w0_s
(if (<= d 3.6e+92)
(*
w0_m
(sqrt
(fma (* (* -0.25 h) (* (* (/ M_m d) D_m) (/ M_m d))) (/ D_m l) 1.0)))
(*
w0_m
(sqrt
(/ (- l (* (* (* (* (/ h d) M_m) (/ M_m d)) (* 0.25 D_m)) D_m)) l))))))D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
double tmp;
if (d <= 3.6e+92) {
tmp = w0_m * sqrt(fma(((-0.25 * h) * (((M_m / d) * D_m) * (M_m / d))), (D_m / l), 1.0));
} else {
tmp = w0_m * sqrt(((l - (((((h / d) * M_m) * (M_m / d)) * (0.25 * D_m)) * D_m)) / l));
}
return w0_s * tmp;
}
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, M_m, D_m, h, l, d) tmp = 0.0 if (d <= 3.6e+92) tmp = Float64(w0_m * sqrt(fma(Float64(Float64(-0.25 * h) * Float64(Float64(Float64(M_m / d) * D_m) * Float64(M_m / d))), Float64(D_m / l), 1.0))); else tmp = Float64(w0_m * sqrt(Float64(Float64(l - Float64(Float64(Float64(Float64(Float64(h / d) * M_m) * Float64(M_m / d)) * Float64(0.25 * D_m)) * D_m)) / l))); end return Float64(w0_s * tmp) end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * If[LessEqual[d, 3.6e+92], N[(w0$95$m * N[Sqrt[N[(N[(N[(-0.25 * h), $MachinePrecision] * N[(N[(N[(M$95$m / d), $MachinePrecision] * D$95$m), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(D$95$m / l), $MachinePrecision] + 1.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(w0$95$m * N[Sqrt[N[(N[(l - N[(N[(N[(N[(N[(h / d), $MachinePrecision] * M$95$m), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision] * N[(0.25 * D$95$m), $MachinePrecision]), $MachinePrecision] * D$95$m), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \begin{array}{l}
\mathbf{if}\;d \leq 3.6 \cdot 10^{+92}:\\
\;\;\;\;w0\_m \cdot \sqrt{\mathsf{fma}\left(\left(-0.25 \cdot h\right) \cdot \left(\left(\frac{M\_m}{d} \cdot D\_m\right) \cdot \frac{M\_m}{d}\right), \frac{D\_m}{\ell}, 1\right)}\\
\mathbf{else}:\\
\;\;\;\;w0\_m \cdot \sqrt{\frac{\ell - \left(\left(\left(\frac{h}{d} \cdot M\_m\right) \cdot \frac{M\_m}{d}\right) \cdot \left(0.25 \cdot D\_m\right)\right) \cdot D\_m}{\ell}}\\
\end{array}
\end{array}
if d < 3.6e92Initial program 84.2%
Taylor expanded in h around inf
fp-cancel-sub-sign-invN/A
metadata-evalN/A
+-commutativeN/A
distribute-lft-inN/A
associate-*r*N/A
rgt-mult-inverseN/A
lower-fma.f64N/A
Applied rewrites63.5%
Applied rewrites66.2%
Applied rewrites82.2%
Applied rewrites82.7%
if 3.6e92 < d Initial program 88.1%
Taylor expanded in l around 0
lower-/.f64N/A
lower--.f64N/A
associate-/l*N/A
associate-*r*N/A
lower-*.f64N/A
lower-*.f64N/A
unpow2N/A
lower-*.f64N/A
unpow2N/A
times-fracN/A
lower-*.f64N/A
lower-/.f64N/A
unpow2N/A
lower-*.f64N/A
lower-/.f6471.4
Applied rewrites71.4%
Applied rewrites92.5%
D_m = (fabs.f64 D) M_m = (fabs.f64 M) w0\_m = (fabs.f64 w0) w0\_s = (copysign.f64 #s(literal 1 binary64) w0) NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function. (FPCore (w0_s w0_m M_m D_m h l d) :precision binary64 (* w0_s (* w0_m 1.0)))
D_m = fabs(D);
M_m = fabs(M);
w0\_m = fabs(w0);
w0\_s = copysign(1.0, w0);
assert(w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d);
double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
return w0_s * (w0_m * 1.0);
}
D_m = private
M_m = private
w0\_m = private
w0\_s = private
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
module fmin_fmax_functions
implicit none
private
public fmax
public fmin
interface fmax
module procedure fmax88
module procedure fmax44
module procedure fmax84
module procedure fmax48
end interface
interface fmin
module procedure fmin88
module procedure fmin44
module procedure fmin84
module procedure fmin48
end interface
contains
real(8) function fmax88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(4) function fmax44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, max(x, y), y /= y), x /= x)
end function
real(8) function fmax84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, max(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmax48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), max(dble(x), y), y /= y), x /= x)
end function
real(8) function fmin88(x, y) result (res)
real(8), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(4) function fmin44(x, y) result (res)
real(4), intent (in) :: x
real(4), intent (in) :: y
res = merge(y, merge(x, min(x, y), y /= y), x /= x)
end function
real(8) function fmin84(x, y) result(res)
real(8), intent (in) :: x
real(4), intent (in) :: y
res = merge(dble(y), merge(x, min(x, dble(y)), y /= y), x /= x)
end function
real(8) function fmin48(x, y) result(res)
real(4), intent (in) :: x
real(8), intent (in) :: y
res = merge(y, merge(dble(x), min(dble(x), y), y /= y), x /= x)
end function
end module
real(8) function code(w0_s, w0_m, m_m, d_m, h, l, d)
use fmin_fmax_functions
real(8), intent (in) :: w0_s
real(8), intent (in) :: w0_m
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_s * (w0_m * 1.0d0)
end function
D_m = Math.abs(D);
M_m = Math.abs(M);
w0\_m = Math.abs(w0);
w0\_s = Math.copySign(1.0, w0);
assert w0_m < M_m && M_m < D_m && D_m < h && h < l && l < d;
public static double code(double w0_s, double w0_m, double M_m, double D_m, double h, double l, double d) {
return w0_s * (w0_m * 1.0);
}
D_m = math.fabs(D) M_m = math.fabs(M) w0\_m = math.fabs(w0) w0\_s = math.copysign(1.0, w0) [w0_m, M_m, D_m, h, l, d] = sort([w0_m, M_m, D_m, h, l, d]) def code(w0_s, w0_m, M_m, D_m, h, l, d): return w0_s * (w0_m * 1.0)
D_m = abs(D) M_m = abs(M) w0\_m = abs(w0) w0\_s = copysign(1.0, w0) w0_m, M_m, D_m, h, l, d = sort([w0_m, M_m, D_m, h, l, d]) function code(w0_s, w0_m, M_m, D_m, h, l, d) return Float64(w0_s * Float64(w0_m * 1.0)) end
D_m = abs(D);
M_m = abs(M);
w0\_m = abs(w0);
w0\_s = sign(w0) * abs(1.0);
w0_m, M_m, D_m, h, l, d = num2cell(sort([w0_m, M_m, D_m, h, l, d])){:}
function tmp = code(w0_s, w0_m, M_m, D_m, h, l, d)
tmp = w0_s * (w0_m * 1.0);
end
D_m = N[Abs[D], $MachinePrecision]
M_m = N[Abs[M], $MachinePrecision]
w0\_m = N[Abs[w0], $MachinePrecision]
w0\_s = N[With[{TMP1 = Abs[1.0], TMP2 = Sign[w0]}, TMP1 * If[TMP2 == 0, 1, TMP2]], $MachinePrecision]
NOTE: w0_m, M_m, D_m, h, l, and d should be sorted in increasing order before calling this function.
code[w0$95$s_, w0$95$m_, M$95$m_, D$95$m_, h_, l_, d_] := N[(w0$95$s * N[(w0$95$m * 1.0), $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
D_m = \left|D\right|
\\
M_m = \left|M\right|
\\
w0\_m = \left|w0\right|
\\
w0\_s = \mathsf{copysign}\left(1, w0\right)
\\
[w0_m, M_m, D_m, h, l, d] = \mathsf{sort}([w0_m, M_m, D_m, h, l, d])\\
\\
w0\_s \cdot \left(w0\_m \cdot 1\right)
\end{array}
Initial program 85.2%
Taylor expanded in M around 0
Applied rewrites70.5%
herbie shell --seed 2024356
(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))))))