Henrywood and Agarwal, Equation (12)

Percentage Accurate: 66.4% → 77.8%

Time: 26.5s

Alternatives: 24

Speedup: 1.5×

Specification

Math

\[\begin{array}{l} \\ \left({\left(\frac{d}{h}\right)}^{\left(\frac{1}{2}\right)} \cdot {\left(\frac{d}{\ell}\right)}^{\left(\frac{1}{2}\right)}\right) \cdot \left(1 - \left(\frac{1}{2} \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right) \cdot \frac{h}{\ell}\right) \end{array} \]

FPCore

(FPCore (d h l M D)
 :precision binary64
 (*
  (* (pow (/ d h) (/ 1.0 2.0)) (pow (/ d l) (/ 1.0 2.0)))
  (- 1.0 (* (* (/ 1.0 2.0) (pow (/ (* M D) (* 2.0 d)) 2.0)) (/ h l)))))

C

double code(double d, double h, double l, double M, double D) {
	return (pow((d / h), (1.0 / 2.0)) * pow((d / l), (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * pow(((M * D) / (2.0 * d)), 2.0)) * (h / l)));
}

Fortran

real(8) function code(d, h, l, m, d_1)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m
    real(8), intent (in) :: d_1
    code = (((d / h) ** (1.0d0 / 2.0d0)) * ((d / l) ** (1.0d0 / 2.0d0))) * (1.0d0 - (((1.0d0 / 2.0d0) * (((m * d_1) / (2.0d0 * d)) ** 2.0d0)) * (h / l)))
end function

Java

public static double code(double d, double h, double l, double M, double D) {
	return (Math.pow((d / h), (1.0 / 2.0)) * Math.pow((d / l), (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * Math.pow(((M * D) / (2.0 * d)), 2.0)) * (h / l)));
}

Python

def code(d, h, l, M, D):
	return (math.pow((d / h), (1.0 / 2.0)) * math.pow((d / l), (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * math.pow(((M * D) / (2.0 * d)), 2.0)) * (h / l)))

Julia

function code(d, h, l, M, D)
	return Float64(Float64((Float64(d / h) ^ Float64(1.0 / 2.0)) * (Float64(d / l) ^ Float64(1.0 / 2.0))) * Float64(1.0 - Float64(Float64(Float64(1.0 / 2.0) * (Float64(Float64(M * D) / Float64(2.0 * d)) ^ 2.0)) * Float64(h / l))))
end

MATLAB

function tmp = code(d, h, l, M, D)
	tmp = (((d / h) ^ (1.0 / 2.0)) * ((d / l) ^ (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * (((M * D) / (2.0 * d)) ^ 2.0)) * (h / l)));
end

Wolfram

code[d_, h_, l_, M_, D_] := N[(N[(N[Power[N[(d / h), $MachinePrecision], N[(1.0 / 2.0), $MachinePrecision]], $MachinePrecision] * N[Power[N[(d / l), $MachinePrecision], N[(1.0 / 2.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(N[(N[(1.0 / 2.0), $MachinePrecision] * N[Power[N[(N[(M * D), $MachinePrecision] / N[(2.0 * d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

Initial Program: 66.4% accurate, 1.0× speedup

Math

\[\begin{array}{l} \\ \left({\left(\frac{d}{h}\right)}^{\left(\frac{1}{2}\right)} \cdot {\left(\frac{d}{\ell}\right)}^{\left(\frac{1}{2}\right)}\right) \cdot \left(1 - \left(\frac{1}{2} \cdot {\left(\frac{M \cdot D}{2 \cdot d}\right)}^{2}\right) \cdot \frac{h}{\ell}\right) \end{array} \]

FPCore

(FPCore (d h l M D)
 :precision binary64
 (*
  (* (pow (/ d h) (/ 1.0 2.0)) (pow (/ d l) (/ 1.0 2.0)))
  (- 1.0 (* (* (/ 1.0 2.0) (pow (/ (* M D) (* 2.0 d)) 2.0)) (/ h l)))))

C

double code(double d, double h, double l, double M, double D) {
	return (pow((d / h), (1.0 / 2.0)) * pow((d / l), (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * pow(((M * D) / (2.0 * d)), 2.0)) * (h / l)));
}

Fortran

real(8) function code(d, h, l, m, d_1)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m
    real(8), intent (in) :: d_1
    code = (((d / h) ** (1.0d0 / 2.0d0)) * ((d / l) ** (1.0d0 / 2.0d0))) * (1.0d0 - (((1.0d0 / 2.0d0) * (((m * d_1) / (2.0d0 * d)) ** 2.0d0)) * (h / l)))
end function

Java

public static double code(double d, double h, double l, double M, double D) {
	return (Math.pow((d / h), (1.0 / 2.0)) * Math.pow((d / l), (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * Math.pow(((M * D) / (2.0 * d)), 2.0)) * (h / l)));
}

Python

def code(d, h, l, M, D):
	return (math.pow((d / h), (1.0 / 2.0)) * math.pow((d / l), (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * math.pow(((M * D) / (2.0 * d)), 2.0)) * (h / l)))

Julia

function code(d, h, l, M, D)
	return Float64(Float64((Float64(d / h) ^ Float64(1.0 / 2.0)) * (Float64(d / l) ^ Float64(1.0 / 2.0))) * Float64(1.0 - Float64(Float64(Float64(1.0 / 2.0) * (Float64(Float64(M * D) / Float64(2.0 * d)) ^ 2.0)) * Float64(h / l))))
end

MATLAB

function tmp = code(d, h, l, M, D)
	tmp = (((d / h) ^ (1.0 / 2.0)) * ((d / l) ^ (1.0 / 2.0))) * (1.0 - (((1.0 / 2.0) * (((M * D) / (2.0 * d)) ^ 2.0)) * (h / l)));
end

Wolfram

code[d_, h_, l_, M_, D_] := N[(N[(N[Power[N[(d / h), $MachinePrecision], N[(1.0 / 2.0), $MachinePrecision]], $MachinePrecision] * N[Power[N[(d / l), $MachinePrecision], N[(1.0 / 2.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(N[(N[(1.0 / 2.0), $MachinePrecision] * N[Power[N[(N[(M * D), $MachinePrecision] / N[(2.0 * d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

Alternative 1: 77.8% accurate, 0.3× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := \sqrt{\frac{d}{h}}\\ t_1 := \left({\left(\frac{d}{h}\right)}^{0.5} \cdot {\left(\frac{d}{\ell}\right)}^{0.5}\right) \cdot \left(1 - \frac{h}{\ell} \cdot \left(0.5 \cdot {\left(\frac{D_m \cdot M_m}{d \cdot 2}\right)}^{2}\right)\right)\\ \mathbf{if}\;t_1 \leq -1 \cdot 10^{-250}:\\ \;\;\;\;\left(t_0 \cdot \frac{1}{\sqrt{\frac{\ell}{d}}}\right) \cdot \left(1 - 0.5 \cdot \left(h \cdot \frac{{\left(D_m \cdot \left(M_m \cdot \frac{0.5}{d}\right)\right)}^{2}}{\ell}\right)\right)\\ \mathbf{elif}\;t_1 \leq 0 \lor \neg \left(t_1 \leq 10^{+241}\right):\\ \;\;\;\;\left|\frac{d}{\sqrt{\ell \cdot h}}\right|\\ \mathbf{else}:\\ \;\;\;\;t_0 \cdot \left(\sqrt{\frac{d}{\ell}} \cdot \left(1 + h \cdot \left({\left(D_m \cdot \frac{M_m \cdot 0.5}{d}\right)}^{2} \cdot \frac{-0.5}{\ell}\right)\right)\right)\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (sqrt (/ d h)))
        (t_1
         (*
          (* (pow (/ d h) 0.5) (pow (/ d l) 0.5))
          (- 1.0 (* (/ h l) (* 0.5 (pow (/ (* D_m M_m) (* d 2.0)) 2.0)))))))
   (if (<= t_1 -1e-250)
     (*
      (* t_0 (/ 1.0 (sqrt (/ l d))))
      (- 1.0 (* 0.5 (* h (/ (pow (* D_m (* M_m (/ 0.5 d))) 2.0) l)))))
     (if (or (<= t_1 0.0) (not (<= t_1 1e+241)))
       (fabs (/ d (sqrt (* l h))))
       (*
        t_0
        (*
         (sqrt (/ d l))
         (+ 1.0 (* h (* (pow (* D_m (/ (* M_m 0.5) d)) 2.0) (/ -0.5 l))))))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = sqrt((d / h));
	double t_1 = (pow((d / h), 0.5) * pow((d / l), 0.5)) * (1.0 - ((h / l) * (0.5 * pow(((D_m * M_m) / (d * 2.0)), 2.0))));
	double tmp;
	if (t_1 <= -1e-250) {
		tmp = (t_0 * (1.0 / sqrt((l / d)))) * (1.0 - (0.5 * (h * (pow((D_m * (M_m * (0.5 / d))), 2.0) / l))));
	} else if ((t_1 <= 0.0) || !(t_1 <= 1e+241)) {
		tmp = fabs((d / sqrt((l * h))));
	} else {
		tmp = t_0 * (sqrt((d / l)) * (1.0 + (h * (pow((D_m * ((M_m * 0.5) / d)), 2.0) * (-0.5 / l)))));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: t_1
    real(8) :: tmp
    t_0 = sqrt((d / h))
    t_1 = (((d / h) ** 0.5d0) * ((d / l) ** 0.5d0)) * (1.0d0 - ((h / l) * (0.5d0 * (((d_m * m_m) / (d * 2.0d0)) ** 2.0d0))))
    if (t_1 <= (-1d-250)) then
        tmp = (t_0 * (1.0d0 / sqrt((l / d)))) * (1.0d0 - (0.5d0 * (h * (((d_m * (m_m * (0.5d0 / d))) ** 2.0d0) / l))))
    else if ((t_1 <= 0.0d0) .or. (.not. (t_1 <= 1d+241))) then
        tmp = abs((d / sqrt((l * h))))
    else
        tmp = t_0 * (sqrt((d / l)) * (1.0d0 + (h * (((d_m * ((m_m * 0.5d0) / d)) ** 2.0d0) * ((-0.5d0) / l)))))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = Math.sqrt((d / h));
	double t_1 = (Math.pow((d / h), 0.5) * Math.pow((d / l), 0.5)) * (1.0 - ((h / l) * (0.5 * Math.pow(((D_m * M_m) / (d * 2.0)), 2.0))));
	double tmp;
	if (t_1 <= -1e-250) {
		tmp = (t_0 * (1.0 / Math.sqrt((l / d)))) * (1.0 - (0.5 * (h * (Math.pow((D_m * (M_m * (0.5 / d))), 2.0) / l))));
	} else if ((t_1 <= 0.0) || !(t_1 <= 1e+241)) {
		tmp = Math.abs((d / Math.sqrt((l * h))));
	} else {
		tmp = t_0 * (Math.sqrt((d / l)) * (1.0 + (h * (Math.pow((D_m * ((M_m * 0.5) / d)), 2.0) * (-0.5 / l)))));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = math.sqrt((d / h))
	t_1 = (math.pow((d / h), 0.5) * math.pow((d / l), 0.5)) * (1.0 - ((h / l) * (0.5 * math.pow(((D_m * M_m) / (d * 2.0)), 2.0))))
	tmp = 0
	if t_1 <= -1e-250:
		tmp = (t_0 * (1.0 / math.sqrt((l / d)))) * (1.0 - (0.5 * (h * (math.pow((D_m * (M_m * (0.5 / d))), 2.0) / l))))
	elif (t_1 <= 0.0) or not (t_1 <= 1e+241):
		tmp = math.fabs((d / math.sqrt((l * h))))
	else:
		tmp = t_0 * (math.sqrt((d / l)) * (1.0 + (h * (math.pow((D_m * ((M_m * 0.5) / d)), 2.0) * (-0.5 / l)))))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = sqrt(Float64(d / h))
	t_1 = Float64(Float64((Float64(d / h) ^ 0.5) * (Float64(d / l) ^ 0.5)) * Float64(1.0 - Float64(Float64(h / l) * Float64(0.5 * (Float64(Float64(D_m * M_m) / Float64(d * 2.0)) ^ 2.0)))))
	tmp = 0.0
	if (t_1 <= -1e-250)
		tmp = Float64(Float64(t_0 * Float64(1.0 / sqrt(Float64(l / d)))) * Float64(1.0 - Float64(0.5 * Float64(h * Float64((Float64(D_m * Float64(M_m * Float64(0.5 / d))) ^ 2.0) / l)))));
	elseif ((t_1 <= 0.0) || !(t_1 <= 1e+241))
		tmp = abs(Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(t_0 * Float64(sqrt(Float64(d / l)) * Float64(1.0 + Float64(h * Float64((Float64(D_m * Float64(Float64(M_m * 0.5) / d)) ^ 2.0) * Float64(-0.5 / l))))));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = sqrt((d / h));
	t_1 = (((d / h) ^ 0.5) * ((d / l) ^ 0.5)) * (1.0 - ((h / l) * (0.5 * (((D_m * M_m) / (d * 2.0)) ^ 2.0))));
	tmp = 0.0;
	if (t_1 <= -1e-250)
		tmp = (t_0 * (1.0 / sqrt((l / d)))) * (1.0 - (0.5 * (h * (((D_m * (M_m * (0.5 / d))) ^ 2.0) / l))));
	elseif ((t_1 <= 0.0) || ~((t_1 <= 1e+241)))
		tmp = abs((d / sqrt((l * h))));
	else
		tmp = t_0 * (sqrt((d / l)) * (1.0 + (h * (((D_m * ((M_m * 0.5) / d)) ^ 2.0) * (-0.5 / l)))));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[(N[(N[Power[N[(d / h), $MachinePrecision], 0.5], $MachinePrecision] * N[Power[N[(d / l), $MachinePrecision], 0.5], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(N[(h / l), $MachinePrecision] * N[(0.5 * N[Power[N[(N[(D$95$m * M$95$m), $MachinePrecision] / N[(d * 2.0), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$1, -1e-250], N[(N[(t$95$0 * N[(1.0 / N[Sqrt[N[(l / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(h * N[(N[Power[N[(D$95$m * N[(M$95$m * N[(0.5 / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[Or[LessEqual[t$95$1, 0.0], N[Not[LessEqual[t$95$1, 1e+241]], $MachinePrecision]], N[Abs[N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision], N[(t$95$0 * N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(h * N[(N[Power[N[(D$95$m * N[(N[(M$95$m * 0.5), $MachinePrecision] / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(-0.5 / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 2: 79.0% accurate, 0.8× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq -1 \cdot 10^{-312}:\\ \;\;\;\;\frac{\sqrt{-d}}{\sqrt{-h}} \cdot \left(\sqrt{\frac{d}{\ell}} \cdot \left(1 + h \cdot \frac{{\left(D_m \cdot \left(M_m \cdot \frac{0.5}{d}\right)\right)}^{2}}{\frac{\ell}{-0.5}}\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{d}{\sqrt{\ell}}}{\sqrt{h}} \cdot \left(1 + -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\right)\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= l -1e-312)
   (*
    (/ (sqrt (- d)) (sqrt (- h)))
    (*
     (sqrt (/ d l))
     (+ 1.0 (* h (/ (pow (* D_m (* M_m (/ 0.5 d))) 2.0) (/ l -0.5))))))
   (*
    (/ (/ d (sqrt l)) (sqrt h))
    (+ 1.0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (l <= -1e-312) {
		tmp = (sqrt(-d) / sqrt(-h)) * (sqrt((d / l)) * (1.0 + (h * (pow((D_m * (M_m * (0.5 / d))), 2.0) / (l / -0.5)))));
	} else {
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + (-0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (l <= (-1d-312)) then
        tmp = (sqrt(-d) / sqrt(-h)) * (sqrt((d / l)) * (1.0d0 + (h * (((d_m * (m_m * (0.5d0 / d))) ** 2.0d0) / (l / (-0.5d0))))))
    else
        tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0d0 + ((-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (l <= -1e-312) {
		tmp = (Math.sqrt(-d) / Math.sqrt(-h)) * (Math.sqrt((d / l)) * (1.0 + (h * (Math.pow((D_m * (M_m * (0.5 / d))), 2.0) / (l / -0.5)))));
	} else {
		tmp = ((d / Math.sqrt(l)) / Math.sqrt(h)) * (1.0 + (-0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if l <= -1e-312:
		tmp = (math.sqrt(-d) / math.sqrt(-h)) * (math.sqrt((d / l)) * (1.0 + (h * (math.pow((D_m * (M_m * (0.5 / d))), 2.0) / (l / -0.5)))))
	else:
		tmp = ((d / math.sqrt(l)) / math.sqrt(h)) * (1.0 + (-0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (l <= -1e-312)
		tmp = Float64(Float64(sqrt(Float64(-d)) / sqrt(Float64(-h))) * Float64(sqrt(Float64(d / l)) * Float64(1.0 + Float64(h * Float64((Float64(D_m * Float64(M_m * Float64(0.5 / d))) ^ 2.0) / Float64(l / -0.5))))));
	else
		tmp = Float64(Float64(Float64(d / sqrt(l)) / sqrt(h)) * Float64(1.0 + Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (l <= -1e-312)
		tmp = (sqrt(-d) / sqrt(-h)) * (sqrt((d / l)) * (1.0 + (h * (((D_m * (M_m * (0.5 / d))) ^ 2.0) / (l / -0.5)))));
	else
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + (-0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l))));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[l, -1e-312], N[(N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-h)], $MachinePrecision]), $MachinePrecision] * N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(h * N[(N[Power[N[(D$95$m * N[(M$95$m * N[(0.5 / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / N[(l / -0.5), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(d / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 3: 77.5% accurate, 0.8× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;h \leq -2.6 \cdot 10^{-261}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \left(\left(1 + h \cdot \frac{{\left(D_m \cdot \left(M_m \cdot \frac{0.5}{d}\right)\right)}^{2}}{\frac{\ell}{-0.5}}\right) \cdot \frac{\sqrt{-d}}{\sqrt{-\ell}}\right)\\ \mathbf{elif}\;h \leq -5 \cdot 10^{-310}:\\ \;\;\;\;\left(d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\right) \cdot \left(-1 - t_0\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{d}{\sqrt{\ell}}}{\sqrt{h}} \cdot \left(1 + t_0\right)\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))
   (if (<= h -2.6e-261)
     (*
      (sqrt (/ d h))
      (*
       (+ 1.0 (* h (/ (pow (* D_m (* M_m (/ 0.5 d))) 2.0) (/ l -0.5))))
       (/ (sqrt (- d)) (sqrt (- l)))))
     (if (<= h -5e-310)
       (* (* d (sqrt (/ (/ 1.0 l) h))) (- -1.0 t_0))
       (* (/ (/ d (sqrt l)) (sqrt h)) (+ 1.0 t_0))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -2.6e-261) {
		tmp = sqrt((d / h)) * ((1.0 + (h * (pow((D_m * (M_m * (0.5 / d))), 2.0) / (l / -0.5)))) * (sqrt(-d) / sqrt(-l)));
	} else if (h <= -5e-310) {
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))
    if (h <= (-2.6d-261)) then
        tmp = sqrt((d / h)) * ((1.0d0 + (h * (((d_m * (m_m * (0.5d0 / d))) ** 2.0d0) / (l / (-0.5d0))))) * (sqrt(-d) / sqrt(-l)))
    else if (h <= (-5d-310)) then
        tmp = (d * sqrt(((1.0d0 / l) / h))) * ((-1.0d0) - t_0)
    else
        tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0d0 + t_0)
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -2.6e-261) {
		tmp = Math.sqrt((d / h)) * ((1.0 + (h * (Math.pow((D_m * (M_m * (0.5 / d))), 2.0) / (l / -0.5)))) * (Math.sqrt(-d) / Math.sqrt(-l)));
	} else if (h <= -5e-310) {
		tmp = (d * Math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / Math.sqrt(l)) / Math.sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = -0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))
	tmp = 0
	if h <= -2.6e-261:
		tmp = math.sqrt((d / h)) * ((1.0 + (h * (math.pow((D_m * (M_m * (0.5 / d))), 2.0) / (l / -0.5)))) * (math.sqrt(-d) / math.sqrt(-l)))
	elif h <= -5e-310:
		tmp = (d * math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0)
	else:
		tmp = ((d / math.sqrt(l)) / math.sqrt(h)) * (1.0 + t_0)
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))
	tmp = 0.0
	if (h <= -2.6e-261)
		tmp = Float64(sqrt(Float64(d / h)) * Float64(Float64(1.0 + Float64(h * Float64((Float64(D_m * Float64(M_m * Float64(0.5 / d))) ^ 2.0) / Float64(l / -0.5)))) * Float64(sqrt(Float64(-d)) / sqrt(Float64(-l)))));
	elseif (h <= -5e-310)
		tmp = Float64(Float64(d * sqrt(Float64(Float64(1.0 / l) / h))) * Float64(-1.0 - t_0));
	else
		tmp = Float64(Float64(Float64(d / sqrt(l)) / sqrt(h)) * Float64(1.0 + t_0));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = -0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l));
	tmp = 0.0;
	if (h <= -2.6e-261)
		tmp = sqrt((d / h)) * ((1.0 + (h * (((D_m * (M_m * (0.5 / d))) ^ 2.0) / (l / -0.5)))) * (sqrt(-d) / sqrt(-l)));
	elseif (h <= -5e-310)
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	else
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[h, -2.6e-261], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[(N[(1.0 + N[(h * N[(N[Power[N[(D$95$m * N[(M$95$m * N[(0.5 / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / N[(l / -0.5), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-l)], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[h, -5e-310], N[(N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(-1.0 - t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[(N[(d / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(1.0 + t$95$0), $MachinePrecision]), $MachinePrecision]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 4: 73.3% accurate, 1.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;h \leq -1.56 \cdot 10^{+116}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \frac{\sqrt{-d}}{\sqrt{-\ell}}\\ \mathbf{elif}\;h \leq -5 \cdot 10^{-310}:\\ \;\;\;\;\left(d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\right) \cdot \left(-1 - t_0\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{d}{\sqrt{\ell}}}{\sqrt{h}} \cdot \left(1 + t_0\right)\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))
   (if (<= h -1.56e+116)
     (* (sqrt (/ d h)) (/ (sqrt (- d)) (sqrt (- l))))
     (if (<= h -5e-310)
       (* (* d (sqrt (/ (/ 1.0 l) h))) (- -1.0 t_0))
       (* (/ (/ d (sqrt l)) (sqrt h)) (+ 1.0 t_0))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -1.56e+116) {
		tmp = sqrt((d / h)) * (sqrt(-d) / sqrt(-l));
	} else if (h <= -5e-310) {
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))
    if (h <= (-1.56d+116)) then
        tmp = sqrt((d / h)) * (sqrt(-d) / sqrt(-l))
    else if (h <= (-5d-310)) then
        tmp = (d * sqrt(((1.0d0 / l) / h))) * ((-1.0d0) - t_0)
    else
        tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0d0 + t_0)
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -1.56e+116) {
		tmp = Math.sqrt((d / h)) * (Math.sqrt(-d) / Math.sqrt(-l));
	} else if (h <= -5e-310) {
		tmp = (d * Math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / Math.sqrt(l)) / Math.sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = -0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))
	tmp = 0
	if h <= -1.56e+116:
		tmp = math.sqrt((d / h)) * (math.sqrt(-d) / math.sqrt(-l))
	elif h <= -5e-310:
		tmp = (d * math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0)
	else:
		tmp = ((d / math.sqrt(l)) / math.sqrt(h)) * (1.0 + t_0)
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))
	tmp = 0.0
	if (h <= -1.56e+116)
		tmp = Float64(sqrt(Float64(d / h)) * Float64(sqrt(Float64(-d)) / sqrt(Float64(-l))));
	elseif (h <= -5e-310)
		tmp = Float64(Float64(d * sqrt(Float64(Float64(1.0 / l) / h))) * Float64(-1.0 - t_0));
	else
		tmp = Float64(Float64(Float64(d / sqrt(l)) / sqrt(h)) * Float64(1.0 + t_0));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = -0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l));
	tmp = 0.0;
	if (h <= -1.56e+116)
		tmp = sqrt((d / h)) * (sqrt(-d) / sqrt(-l));
	elseif (h <= -5e-310)
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	else
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[h, -1.56e+116], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-l)], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[h, -5e-310], N[(N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(-1.0 - t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[(N[(d / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(1.0 + t$95$0), $MachinePrecision]), $MachinePrecision]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 5: 76.5% accurate, 1.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;h \leq -9.4 \cdot 10^{+114}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \left(\sqrt{\frac{d}{\ell}} \cdot \left(1 + h \cdot \left({\left(D_m \cdot \frac{M_m \cdot 0.5}{d}\right)}^{2} \cdot \frac{-0.5}{\ell}\right)\right)\right)\\ \mathbf{elif}\;h \leq -5 \cdot 10^{-310}:\\ \;\;\;\;\left(d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\right) \cdot \left(-1 - t_0\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{d}{\sqrt{\ell}}}{\sqrt{h}} \cdot \left(1 + t_0\right)\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))
   (if (<= h -9.4e+114)
     (*
      (sqrt (/ d h))
      (*
       (sqrt (/ d l))
       (+ 1.0 (* h (* (pow (* D_m (/ (* M_m 0.5) d)) 2.0) (/ -0.5 l))))))
     (if (<= h -5e-310)
       (* (* d (sqrt (/ (/ 1.0 l) h))) (- -1.0 t_0))
       (* (/ (/ d (sqrt l)) (sqrt h)) (+ 1.0 t_0))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -9.4e+114) {
		tmp = sqrt((d / h)) * (sqrt((d / l)) * (1.0 + (h * (pow((D_m * ((M_m * 0.5) / d)), 2.0) * (-0.5 / l)))));
	} else if (h <= -5e-310) {
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))
    if (h <= (-9.4d+114)) then
        tmp = sqrt((d / h)) * (sqrt((d / l)) * (1.0d0 + (h * (((d_m * ((m_m * 0.5d0) / d)) ** 2.0d0) * ((-0.5d0) / l)))))
    else if (h <= (-5d-310)) then
        tmp = (d * sqrt(((1.0d0 / l) / h))) * ((-1.0d0) - t_0)
    else
        tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0d0 + t_0)
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -9.4e+114) {
		tmp = Math.sqrt((d / h)) * (Math.sqrt((d / l)) * (1.0 + (h * (Math.pow((D_m * ((M_m * 0.5) / d)), 2.0) * (-0.5 / l)))));
	} else if (h <= -5e-310) {
		tmp = (d * Math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / Math.sqrt(l)) / Math.sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = -0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))
	tmp = 0
	if h <= -9.4e+114:
		tmp = math.sqrt((d / h)) * (math.sqrt((d / l)) * (1.0 + (h * (math.pow((D_m * ((M_m * 0.5) / d)), 2.0) * (-0.5 / l)))))
	elif h <= -5e-310:
		tmp = (d * math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0)
	else:
		tmp = ((d / math.sqrt(l)) / math.sqrt(h)) * (1.0 + t_0)
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))
	tmp = 0.0
	if (h <= -9.4e+114)
		tmp = Float64(sqrt(Float64(d / h)) * Float64(sqrt(Float64(d / l)) * Float64(1.0 + Float64(h * Float64((Float64(D_m * Float64(Float64(M_m * 0.5) / d)) ^ 2.0) * Float64(-0.5 / l))))));
	elseif (h <= -5e-310)
		tmp = Float64(Float64(d * sqrt(Float64(Float64(1.0 / l) / h))) * Float64(-1.0 - t_0));
	else
		tmp = Float64(Float64(Float64(d / sqrt(l)) / sqrt(h)) * Float64(1.0 + t_0));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = -0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l));
	tmp = 0.0;
	if (h <= -9.4e+114)
		tmp = sqrt((d / h)) * (sqrt((d / l)) * (1.0 + (h * (((D_m * ((M_m * 0.5) / d)) ^ 2.0) * (-0.5 / l)))));
	elseif (h <= -5e-310)
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	else
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[h, -9.4e+114], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(h * N[(N[Power[N[(D$95$m * N[(N[(M$95$m * 0.5), $MachinePrecision] / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(-0.5 / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[h, -5e-310], N[(N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(-1.0 - t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[(N[(d / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(1.0 + t$95$0), $MachinePrecision]), $MachinePrecision]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 6: 76.1% accurate, 1.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;h \leq -3100000000000:\\ \;\;\;\;\left(\sqrt{\frac{d}{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot \left(1 - 0.5 \cdot \frac{h \cdot {\left(\left(M_m \cdot 0.5\right) \cdot \frac{D_m}{d}\right)}^{2}}{\ell}\right)\\ \mathbf{elif}\;h \leq -5 \cdot 10^{-310}:\\ \;\;\;\;\left(d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\right) \cdot \left(-1 - t_0\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{\frac{d}{\sqrt{\ell}}}{\sqrt{h}} \cdot \left(1 + t_0\right)\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))
   (if (<= h -3100000000000.0)
     (*
      (* (sqrt (/ d l)) (sqrt (/ d h)))
      (- 1.0 (* 0.5 (/ (* h (pow (* (* M_m 0.5) (/ D_m d)) 2.0)) l))))
     (if (<= h -5e-310)
       (* (* d (sqrt (/ (/ 1.0 l) h))) (- -1.0 t_0))
       (* (/ (/ d (sqrt l)) (sqrt h)) (+ 1.0 t_0))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -3100000000000.0) {
		tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0 - (0.5 * ((h * pow(((M_m * 0.5) * (D_m / d)), 2.0)) / l)));
	} else if (h <= -5e-310) {
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))
    if (h <= (-3100000000000.0d0)) then
        tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0d0 - (0.5d0 * ((h * (((m_m * 0.5d0) * (d_m / d)) ** 2.0d0)) / l)))
    else if (h <= (-5d-310)) then
        tmp = (d * sqrt(((1.0d0 / l) / h))) * ((-1.0d0) - t_0)
    else
        tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0d0 + t_0)
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (h <= -3100000000000.0) {
		tmp = (Math.sqrt((d / l)) * Math.sqrt((d / h))) * (1.0 - (0.5 * ((h * Math.pow(((M_m * 0.5) * (D_m / d)), 2.0)) / l)));
	} else if (h <= -5e-310) {
		tmp = (d * Math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else {
		tmp = ((d / Math.sqrt(l)) / Math.sqrt(h)) * (1.0 + t_0);
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = -0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))
	tmp = 0
	if h <= -3100000000000.0:
		tmp = (math.sqrt((d / l)) * math.sqrt((d / h))) * (1.0 - (0.5 * ((h * math.pow(((M_m * 0.5) * (D_m / d)), 2.0)) / l)))
	elif h <= -5e-310:
		tmp = (d * math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0)
	else:
		tmp = ((d / math.sqrt(l)) / math.sqrt(h)) * (1.0 + t_0)
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))
	tmp = 0.0
	if (h <= -3100000000000.0)
		tmp = Float64(Float64(sqrt(Float64(d / l)) * sqrt(Float64(d / h))) * Float64(1.0 - Float64(0.5 * Float64(Float64(h * (Float64(Float64(M_m * 0.5) * Float64(D_m / d)) ^ 2.0)) / l))));
	elseif (h <= -5e-310)
		tmp = Float64(Float64(d * sqrt(Float64(Float64(1.0 / l) / h))) * Float64(-1.0 - t_0));
	else
		tmp = Float64(Float64(Float64(d / sqrt(l)) / sqrt(h)) * Float64(1.0 + t_0));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = -0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l));
	tmp = 0.0;
	if (h <= -3100000000000.0)
		tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0 - (0.5 * ((h * (((M_m * 0.5) * (D_m / d)) ^ 2.0)) / l)));
	elseif (h <= -5e-310)
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	else
		tmp = ((d / sqrt(l)) / sqrt(h)) * (1.0 + t_0);
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[h, -3100000000000.0], N[(N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(h * N[Power[N[(N[(M$95$m * 0.5), $MachinePrecision] * N[(D$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[h, -5e-310], N[(N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(-1.0 - t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[(N[(d / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(1.0 + t$95$0), $MachinePrecision]), $MachinePrecision]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 7: 70.9% accurate, 1.1× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;\ell \leq -3.3 \cdot 10^{+106}:\\ \;\;\;\;\frac{\sqrt{-d}}{\sqrt{-h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq -2.95 \cdot 10^{-167}:\\ \;\;\;\;\left(d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\right) \cdot \left(-1 - t_0\right)\\ \mathbf{elif}\;\ell \leq 1.2 \cdot 10^{-153}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}} \cdot \left(1 + -0.5 \cdot \frac{h \cdot \left({\left(M_m \cdot \frac{D_m}{d}\right)}^{2} \cdot 0.25\right)}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 9.5 \cdot 10^{+234}:\\ \;\;\;\;\left(1 + t_0\right) \cdot \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))
   (if (<= l -3.3e+106)
     (* (/ (sqrt (- d)) (sqrt (- h))) (sqrt (/ d l)))
     (if (<= l -2.95e-167)
       (* (* d (sqrt (/ (/ 1.0 l) h))) (- -1.0 t_0))
       (if (<= l 1.2e-153)
         (*
          (sqrt (* (/ d l) (/ d h)))
          (+ 1.0 (* -0.5 (/ (* h (* (pow (* M_m (/ D_m d)) 2.0) 0.25)) l))))
         (if (<= l 9.5e+234)
           (* (+ 1.0 t_0) (/ d (sqrt (* l h))))
           (/ d (* (sqrt l) (sqrt h)))))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (l <= -3.3e+106) {
		tmp = (sqrt(-d) / sqrt(-h)) * sqrt((d / l));
	} else if (l <= -2.95e-167) {
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else if (l <= 1.2e-153) {
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)));
	} else if (l <= 9.5e+234) {
		tmp = (1.0 + t_0) * (d / sqrt((l * h)));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))
    if (l <= (-3.3d+106)) then
        tmp = (sqrt(-d) / sqrt(-h)) * sqrt((d / l))
    else if (l <= (-2.95d-167)) then
        tmp = (d * sqrt(((1.0d0 / l) / h))) * ((-1.0d0) - t_0)
    else if (l <= 1.2d-153) then
        tmp = sqrt(((d / l) * (d / h))) * (1.0d0 + ((-0.5d0) * ((h * (((m_m * (d_m / d)) ** 2.0d0) * 0.25d0)) / l)))
    else if (l <= 9.5d+234) then
        tmp = (1.0d0 + t_0) * (d / sqrt((l * h)))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (l <= -3.3e+106) {
		tmp = (Math.sqrt(-d) / Math.sqrt(-h)) * Math.sqrt((d / l));
	} else if (l <= -2.95e-167) {
		tmp = (d * Math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else if (l <= 1.2e-153) {
		tmp = Math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (Math.pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)));
	} else if (l <= 9.5e+234) {
		tmp = (1.0 + t_0) * (d / Math.sqrt((l * h)));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = -0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))
	tmp = 0
	if l <= -3.3e+106:
		tmp = (math.sqrt(-d) / math.sqrt(-h)) * math.sqrt((d / l))
	elif l <= -2.95e-167:
		tmp = (d * math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0)
	elif l <= 1.2e-153:
		tmp = math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (math.pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)))
	elif l <= 9.5e+234:
		tmp = (1.0 + t_0) * (d / math.sqrt((l * h)))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))
	tmp = 0.0
	if (l <= -3.3e+106)
		tmp = Float64(Float64(sqrt(Float64(-d)) / sqrt(Float64(-h))) * sqrt(Float64(d / l)));
	elseif (l <= -2.95e-167)
		tmp = Float64(Float64(d * sqrt(Float64(Float64(1.0 / l) / h))) * Float64(-1.0 - t_0));
	elseif (l <= 1.2e-153)
		tmp = Float64(sqrt(Float64(Float64(d / l) * Float64(d / h))) * Float64(1.0 + Float64(-0.5 * Float64(Float64(h * Float64((Float64(M_m * Float64(D_m / d)) ^ 2.0) * 0.25)) / l))));
	elseif (l <= 9.5e+234)
		tmp = Float64(Float64(1.0 + t_0) * Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = -0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l));
	tmp = 0.0;
	if (l <= -3.3e+106)
		tmp = (sqrt(-d) / sqrt(-h)) * sqrt((d / l));
	elseif (l <= -2.95e-167)
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	elseif (l <= 1.2e-153)
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (((M_m * (D_m / d)) ^ 2.0) * 0.25)) / l)));
	elseif (l <= 9.5e+234)
		tmp = (1.0 + t_0) * (d / sqrt((l * h)));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[l, -3.3e+106], N[(N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-h)], $MachinePrecision]), $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[l, -2.95e-167], N[(N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(-1.0 - t$95$0), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 1.2e-153], N[(N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[(h * N[(N[Power[N[(M$95$m * N[(D$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * 0.25), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 9.5e+234], N[(N[(1.0 + t$95$0), $MachinePrecision] * N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 8: 70.7% accurate, 1.5× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;\ell \leq -8.8 \cdot 10^{-165}:\\ \;\;\;\;\left(d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\right) \cdot \left(-1 - t_0\right)\\ \mathbf{elif}\;\ell \leq 2.4 \cdot 10^{-154}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}} \cdot \left(1 + -0.5 \cdot \frac{h \cdot \left({\left(M_m \cdot \frac{D_m}{d}\right)}^{2} \cdot 0.25\right)}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 9 \cdot 10^{+234}:\\ \;\;\;\;\left(1 + t_0\right) \cdot \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l)))))
   (if (<= l -8.8e-165)
     (* (* d (sqrt (/ (/ 1.0 l) h))) (- -1.0 t_0))
     (if (<= l 2.4e-154)
       (*
        (sqrt (* (/ d l) (/ d h)))
        (+ 1.0 (* -0.5 (/ (* h (* (pow (* M_m (/ D_m d)) 2.0) 0.25)) l))))
       (if (<= l 9e+234)
         (* (+ 1.0 t_0) (/ d (sqrt (* l h))))
         (/ d (* (sqrt l) (sqrt h))))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (l <= -8.8e-165) {
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else if (l <= 2.4e-154) {
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)));
	} else if (l <= 9e+234) {
		tmp = (1.0 + t_0) * (d / sqrt((l * h)));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l))
    if (l <= (-8.8d-165)) then
        tmp = (d * sqrt(((1.0d0 / l) / h))) * ((-1.0d0) - t_0)
    else if (l <= 2.4d-154) then
        tmp = sqrt(((d / l) * (d / h))) * (1.0d0 + ((-0.5d0) * ((h * (((m_m * (d_m / d)) ** 2.0d0) * 0.25d0)) / l)))
    else if (l <= 9d+234) then
        tmp = (1.0d0 + t_0) * (d / sqrt((l * h)))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = -0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l));
	double tmp;
	if (l <= -8.8e-165) {
		tmp = (d * Math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	} else if (l <= 2.4e-154) {
		tmp = Math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (Math.pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)));
	} else if (l <= 9e+234) {
		tmp = (1.0 + t_0) * (d / Math.sqrt((l * h)));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = -0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l))
	tmp = 0
	if l <= -8.8e-165:
		tmp = (d * math.sqrt(((1.0 / l) / h))) * (-1.0 - t_0)
	elif l <= 2.4e-154:
		tmp = math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (math.pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)))
	elif l <= 9e+234:
		tmp = (1.0 + t_0) * (d / math.sqrt((l * h)))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))
	tmp = 0.0
	if (l <= -8.8e-165)
		tmp = Float64(Float64(d * sqrt(Float64(Float64(1.0 / l) / h))) * Float64(-1.0 - t_0));
	elseif (l <= 2.4e-154)
		tmp = Float64(sqrt(Float64(Float64(d / l) * Float64(d / h))) * Float64(1.0 + Float64(-0.5 * Float64(Float64(h * Float64((Float64(M_m * Float64(D_m / d)) ^ 2.0) * 0.25)) / l))));
	elseif (l <= 9e+234)
		tmp = Float64(Float64(1.0 + t_0) * Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = -0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l));
	tmp = 0.0;
	if (l <= -8.8e-165)
		tmp = (d * sqrt(((1.0 / l) / h))) * (-1.0 - t_0);
	elseif (l <= 2.4e-154)
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (((M_m * (D_m / d)) ^ 2.0) * 0.25)) / l)));
	elseif (l <= 9e+234)
		tmp = (1.0 + t_0) * (d / sqrt((l * h)));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[l, -8.8e-165], N[(N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(-1.0 - t$95$0), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 2.4e-154], N[(N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[(h * N[(N[Power[N[(M$95$m * N[(D$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * 0.25), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 9e+234], N[(N[(1.0 + t$95$0), $MachinePrecision] * N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 9: 56.5% accurate, 1.5× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq 1.9 \cdot 10^{-300}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{elif}\;\ell \leq 10^{+235}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(1 + -0.5 \cdot \left({\left(M_m \cdot \frac{D_m}{d}\right)}^{2} \cdot \left(\frac{h}{\ell} \cdot 0.25\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= l 1.9e-300)
   (* d (- (sqrt (/ (/ 1.0 l) h))))
   (if (<= l 1e+235)
     (*
      (/ d (sqrt (* l h)))
      (+ 1.0 (* -0.5 (* (pow (* M_m (/ D_m d)) 2.0) (* (/ h l) 0.25)))))
     (/ d (* (sqrt l) (sqrt h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (l <= 1.9e-300) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else if (l <= 1e+235) {
		tmp = (d / sqrt((l * h))) * (1.0 + (-0.5 * (pow((M_m * (D_m / d)), 2.0) * ((h / l) * 0.25))));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (l <= 1.9d-300) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else if (l <= 1d+235) then
        tmp = (d / sqrt((l * h))) * (1.0d0 + ((-0.5d0) * (((m_m * (d_m / d)) ** 2.0d0) * ((h / l) * 0.25d0))))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (l <= 1.9e-300) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else if (l <= 1e+235) {
		tmp = (d / Math.sqrt((l * h))) * (1.0 + (-0.5 * (Math.pow((M_m * (D_m / d)), 2.0) * ((h / l) * 0.25))));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if l <= 1.9e-300:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	elif l <= 1e+235:
		tmp = (d / math.sqrt((l * h))) * (1.0 + (-0.5 * (math.pow((M_m * (D_m / d)), 2.0) * ((h / l) * 0.25))))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (l <= 1.9e-300)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	elseif (l <= 1e+235)
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(1.0 + Float64(-0.5 * Float64((Float64(M_m * Float64(D_m / d)) ^ 2.0) * Float64(Float64(h / l) * 0.25)))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (l <= 1.9e-300)
		tmp = d * -sqrt(((1.0 / l) / h));
	elseif (l <= 1e+235)
		tmp = (d / sqrt((l * h))) * (1.0 + (-0.5 * (((M_m * (D_m / d)) ^ 2.0) * ((h / l) * 0.25))));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[l, 1.9e-300], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[l, 1e+235], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[Power[N[(M$95$m * N[(D$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(N[(h / l), $MachinePrecision] * 0.25), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 10: 55.8% accurate, 1.5× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;d \leq 1.02 \cdot 10^{-305}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{elif}\;d \leq 7 \cdot 10^{+145}:\\ \;\;\;\;\left(1 + -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) \cdot \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= d 1.02e-305)
   (* d (- (sqrt (/ (/ 1.0 l) h))))
   (if (<= d 7e+145)
     (*
      (+ 1.0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l))))
      (/ d (sqrt (* l h))))
     (/ d (* (sqrt l) (sqrt h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= 1.02e-305) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else if (d <= 7e+145) {
		tmp = (1.0 + (-0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / sqrt((l * h)));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (d <= 1.02d-305) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else if (d <= 7d+145) then
        tmp = (1.0d0 + ((-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l)))) * (d / sqrt((l * h)))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= 1.02e-305) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else if (d <= 7e+145) {
		tmp = (1.0 + (-0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / Math.sqrt((l * h)));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if d <= 1.02e-305:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	elif d <= 7e+145:
		tmp = (1.0 + (-0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / math.sqrt((l * h)))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (d <= 1.02e-305)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	elseif (d <= 7e+145)
		tmp = Float64(Float64(1.0 + Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))) * Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (d <= 1.02e-305)
		tmp = d * -sqrt(((1.0 / l) / h));
	elseif (d <= 7e+145)
		tmp = (1.0 + (-0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l)))) * (d / sqrt((l * h)));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[d, 1.02e-305], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[d, 7e+145], N[(N[(1.0 + N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 11: 63.1% accurate, 1.5× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} t_0 := 1 + -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\\ \mathbf{if}\;d \leq -6 \cdot 10^{-308}:\\ \;\;\;\;t_0 \cdot \sqrt{\frac{d}{\ell} \cdot \frac{d}{h}}\\ \mathbf{elif}\;d \leq 1.7 \cdot 10^{+152}:\\ \;\;\;\;t_0 \cdot \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (let* ((t_0 (+ 1.0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l))))))
   (if (<= d -6e-308)
     (* t_0 (sqrt (* (/ d l) (/ d h))))
     (if (<= d 1.7e+152)
       (* t_0 (/ d (sqrt (* l h))))
       (/ d (* (sqrt l) (sqrt h)))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = 1.0 + (-0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)));
	double tmp;
	if (d <= -6e-308) {
		tmp = t_0 * sqrt(((d / l) * (d / h)));
	} else if (d <= 1.7e+152) {
		tmp = t_0 * (d / sqrt((l * h)));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = 1.0d0 + ((-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l)))
    if (d <= (-6d-308)) then
        tmp = t_0 * sqrt(((d / l) * (d / h)))
    else if (d <= 1.7d+152) then
        tmp = t_0 * (d / sqrt((l * h)))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double t_0 = 1.0 + (-0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)));
	double tmp;
	if (d <= -6e-308) {
		tmp = t_0 * Math.sqrt(((d / l) * (d / h)));
	} else if (d <= 1.7e+152) {
		tmp = t_0 * (d / Math.sqrt((l * h)));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	t_0 = 1.0 + (-0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))
	tmp = 0
	if d <= -6e-308:
		tmp = t_0 * math.sqrt(((d / l) * (d / h)))
	elif d <= 1.7e+152:
		tmp = t_0 * (d / math.sqrt((l * h)))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	t_0 = Float64(1.0 + Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l))))
	tmp = 0.0
	if (d <= -6e-308)
		tmp = Float64(t_0 * sqrt(Float64(Float64(d / l) * Float64(d / h))));
	elseif (d <= 1.7e+152)
		tmp = Float64(t_0 * Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	t_0 = 1.0 + (-0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l)));
	tmp = 0.0;
	if (d <= -6e-308)
		tmp = t_0 * sqrt(((d / l) * (d / h)));
	elseif (d <= 1.7e+152)
		tmp = t_0 * (d / sqrt((l * h)));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := Block[{t$95$0 = N[(1.0 + N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[d, -6e-308], N[(t$95$0 * N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 1.7e+152], N[(t$95$0 * N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 12: 63.2% accurate, 1.5× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;d \leq -3.55 \cdot 10^{-307}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}} \cdot \left(1 + -0.5 \cdot \left(\frac{h}{\ell} \cdot {\left(\frac{\frac{D_m \cdot M_m}{d}}{2}\right)}^{2}\right)\right)\\ \mathbf{elif}\;d \leq 3.2 \cdot 10^{+144}:\\ \;\;\;\;\left(1 + -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) \cdot \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= d -3.55e-307)
   (*
    (sqrt (* (/ d l) (/ d h)))
    (+ 1.0 (* -0.5 (* (/ h l) (pow (/ (/ (* D_m M_m) d) 2.0) 2.0)))))
   (if (<= d 3.2e+144)
     (*
      (+ 1.0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l))))
      (/ d (sqrt (* l h))))
     (/ d (* (sqrt l) (sqrt h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -3.55e-307) {
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h / l) * pow((((D_m * M_m) / d) / 2.0), 2.0))));
	} else if (d <= 3.2e+144) {
		tmp = (1.0 + (-0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / sqrt((l * h)));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (d <= (-3.55d-307)) then
        tmp = sqrt(((d / l) * (d / h))) * (1.0d0 + ((-0.5d0) * ((h / l) * ((((d_m * m_m) / d) / 2.0d0) ** 2.0d0))))
    else if (d <= 3.2d+144) then
        tmp = (1.0d0 + ((-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l)))) * (d / sqrt((l * h)))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -3.55e-307) {
		tmp = Math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h / l) * Math.pow((((D_m * M_m) / d) / 2.0), 2.0))));
	} else if (d <= 3.2e+144) {
		tmp = (1.0 + (-0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / Math.sqrt((l * h)));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if d <= -3.55e-307:
		tmp = math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h / l) * math.pow((((D_m * M_m) / d) / 2.0), 2.0))))
	elif d <= 3.2e+144:
		tmp = (1.0 + (-0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / math.sqrt((l * h)))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (d <= -3.55e-307)
		tmp = Float64(sqrt(Float64(Float64(d / l) * Float64(d / h))) * Float64(1.0 + Float64(-0.5 * Float64(Float64(h / l) * (Float64(Float64(Float64(D_m * M_m) / d) / 2.0) ^ 2.0)))));
	elseif (d <= 3.2e+144)
		tmp = Float64(Float64(1.0 + Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))) * Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (d <= -3.55e-307)
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h / l) * ((((D_m * M_m) / d) / 2.0) ^ 2.0))));
	elseif (d <= 3.2e+144)
		tmp = (1.0 + (-0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l)))) * (d / sqrt((l * h)));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[d, -3.55e-307], N[(N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[(h / l), $MachinePrecision] * N[Power[N[(N[(N[(D$95$m * M$95$m), $MachinePrecision] / d), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 3.2e+144], N[(N[(1.0 + N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 13: 64.8% accurate, 1.5× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq 9.6 \cdot 10^{-150}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}} \cdot \left(1 + -0.5 \cdot \frac{h \cdot \left({\left(M_m \cdot \frac{D_m}{d}\right)}^{2} \cdot 0.25\right)}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 9 \cdot 10^{+234}:\\ \;\;\;\;\left(1 + -0.5 \cdot \left({\left(\frac{\frac{D_m}{\frac{d}{M_m}}}{2}\right)}^{2} \cdot \frac{h}{\ell}\right)\right) \cdot \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= l 9.6e-150)
   (*
    (sqrt (* (/ d l) (/ d h)))
    (+ 1.0 (* -0.5 (/ (* h (* (pow (* M_m (/ D_m d)) 2.0) 0.25)) l))))
   (if (<= l 9e+234)
     (*
      (+ 1.0 (* -0.5 (* (pow (/ (/ D_m (/ d M_m)) 2.0) 2.0) (/ h l))))
      (/ d (sqrt (* l h))))
     (/ d (* (sqrt l) (sqrt h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (l <= 9.6e-150) {
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)));
	} else if (l <= 9e+234) {
		tmp = (1.0 + (-0.5 * (pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / sqrt((l * h)));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (l <= 9.6d-150) then
        tmp = sqrt(((d / l) * (d / h))) * (1.0d0 + ((-0.5d0) * ((h * (((m_m * (d_m / d)) ** 2.0d0) * 0.25d0)) / l)))
    else if (l <= 9d+234) then
        tmp = (1.0d0 + ((-0.5d0) * ((((d_m / (d / m_m)) / 2.0d0) ** 2.0d0) * (h / l)))) * (d / sqrt((l * h)))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (l <= 9.6e-150) {
		tmp = Math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (Math.pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)));
	} else if (l <= 9e+234) {
		tmp = (1.0 + (-0.5 * (Math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / Math.sqrt((l * h)));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if l <= 9.6e-150:
		tmp = math.sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (math.pow((M_m * (D_m / d)), 2.0) * 0.25)) / l)))
	elif l <= 9e+234:
		tmp = (1.0 + (-0.5 * (math.pow(((D_m / (d / M_m)) / 2.0), 2.0) * (h / l)))) * (d / math.sqrt((l * h)))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (l <= 9.6e-150)
		tmp = Float64(sqrt(Float64(Float64(d / l) * Float64(d / h))) * Float64(1.0 + Float64(-0.5 * Float64(Float64(h * Float64((Float64(M_m * Float64(D_m / d)) ^ 2.0) * 0.25)) / l))));
	elseif (l <= 9e+234)
		tmp = Float64(Float64(1.0 + Float64(-0.5 * Float64((Float64(Float64(D_m / Float64(d / M_m)) / 2.0) ^ 2.0) * Float64(h / l)))) * Float64(d / sqrt(Float64(l * h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (l <= 9.6e-150)
		tmp = sqrt(((d / l) * (d / h))) * (1.0 + (-0.5 * ((h * (((M_m * (D_m / d)) ^ 2.0) * 0.25)) / l)));
	elseif (l <= 9e+234)
		tmp = (1.0 + (-0.5 * ((((D_m / (d / M_m)) / 2.0) ^ 2.0) * (h / l)))) * (d / sqrt((l * h)));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[l, 9.6e-150], N[(N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[(h * N[(N[Power[N[(M$95$m * N[(D$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * 0.25), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 9e+234], N[(N[(1.0 + N[(-0.5 * N[(N[Power[N[(N[(D$95$m / N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / 2.0), $MachinePrecision], 2.0], $MachinePrecision] * N[(h / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 14: 47.0% accurate, 1.6× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;h \leq -4.8 \cdot 10^{+151}:\\ \;\;\;\;\sqrt{\frac{d}{\ell}} \cdot \frac{1}{\sqrt{\frac{h}{d}}}\\ \mathbf{elif}\;h \leq 1.6 \cdot 10^{-299}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= h -4.8e+151)
   (* (sqrt (/ d l)) (/ 1.0 (sqrt (/ h d))))
   (if (<= h 1.6e-299)
     (* d (- (sqrt (/ (/ 1.0 l) h))))
     (/ d (* (sqrt l) (sqrt h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= -4.8e+151) {
		tmp = sqrt((d / l)) * (1.0 / sqrt((h / d)));
	} else if (h <= 1.6e-299) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (h <= (-4.8d+151)) then
        tmp = sqrt((d / l)) * (1.0d0 / sqrt((h / d)))
    else if (h <= 1.6d-299) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= -4.8e+151) {
		tmp = Math.sqrt((d / l)) * (1.0 / Math.sqrt((h / d)));
	} else if (h <= 1.6e-299) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if h <= -4.8e+151:
		tmp = math.sqrt((d / l)) * (1.0 / math.sqrt((h / d)))
	elif h <= 1.6e-299:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (h <= -4.8e+151)
		tmp = Float64(sqrt(Float64(d / l)) * Float64(1.0 / sqrt(Float64(h / d))));
	elseif (h <= 1.6e-299)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (h <= -4.8e+151)
		tmp = sqrt((d / l)) * (1.0 / sqrt((h / d)));
	elseif (h <= 1.6e-299)
		tmp = d * -sqrt(((1.0 / l) / h));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[h, -4.8e+151], N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(1.0 / N[Sqrt[N[(h / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[h, 1.6e-299], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 15: 41.9% accurate, 1.6× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;h \leq 1.25 \cdot 10^{-287}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{elif}\;h \leq 2.9 \cdot 10^{+84}:\\ \;\;\;\;\left|\frac{d}{\sqrt{\ell \cdot h}}\right|\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= h 1.25e-287)
   (* d (- (sqrt (/ (/ 1.0 l) h))))
   (if (<= h 2.9e+84) (fabs (/ d (sqrt (* l h)))) (sqrt (* (/ d l) (/ d h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 1.25e-287) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else if (h <= 2.9e+84) {
		tmp = fabs((d / sqrt((l * h))));
	} else {
		tmp = sqrt(((d / l) * (d / h)));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (h <= 1.25d-287) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else if (h <= 2.9d+84) then
        tmp = abs((d / sqrt((l * h))))
    else
        tmp = sqrt(((d / l) * (d / h)))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 1.25e-287) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else if (h <= 2.9e+84) {
		tmp = Math.abs((d / Math.sqrt((l * h))));
	} else {
		tmp = Math.sqrt(((d / l) * (d / h)));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if h <= 1.25e-287:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	elif h <= 2.9e+84:
		tmp = math.fabs((d / math.sqrt((l * h))))
	else:
		tmp = math.sqrt(((d / l) * (d / h)))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (h <= 1.25e-287)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	elseif (h <= 2.9e+84)
		tmp = abs(Float64(d / sqrt(Float64(l * h))));
	else
		tmp = sqrt(Float64(Float64(d / l) * Float64(d / h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (h <= 1.25e-287)
		tmp = d * -sqrt(((1.0 / l) / h));
	elseif (h <= 2.9e+84)
		tmp = abs((d / sqrt((l * h))));
	else
		tmp = sqrt(((d / l) * (d / h)));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[h, 1.25e-287], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[h, 2.9e+84], N[Abs[N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]], $MachinePrecision], N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 16: 47.1% accurate, 1.6× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;h \leq -4.7 \cdot 10^{+149}:\\ \;\;\;\;\sqrt{\frac{d}{\ell}} \cdot \sqrt{\frac{d}{h}}\\ \mathbf{elif}\;h \leq 8.5 \cdot 10^{-302}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= h -4.7e+149)
   (* (sqrt (/ d l)) (sqrt (/ d h)))
   (if (<= h 8.5e-302)
     (* d (- (sqrt (/ (/ 1.0 l) h))))
     (/ d (* (sqrt l) (sqrt h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= -4.7e+149) {
		tmp = sqrt((d / l)) * sqrt((d / h));
	} else if (h <= 8.5e-302) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (h <= (-4.7d+149)) then
        tmp = sqrt((d / l)) * sqrt((d / h))
    else if (h <= 8.5d-302) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= -4.7e+149) {
		tmp = Math.sqrt((d / l)) * Math.sqrt((d / h));
	} else if (h <= 8.5e-302) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if h <= -4.7e+149:
		tmp = math.sqrt((d / l)) * math.sqrt((d / h))
	elif h <= 8.5e-302:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (h <= -4.7e+149)
		tmp = Float64(sqrt(Float64(d / l)) * sqrt(Float64(d / h)));
	elseif (h <= 8.5e-302)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (h <= -4.7e+149)
		tmp = sqrt((d / l)) * sqrt((d / h));
	elseif (h <= 8.5e-302)
		tmp = d * -sqrt(((1.0 / l) / h));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[h, -4.7e+149], N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[h, 8.5e-302], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 17: 46.2% accurate, 1.6× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;h \leq 8.5 \cdot 10^{-302}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= h 8.5e-302)
   (* d (- (sqrt (/ (/ 1.0 l) h))))
   (/ d (* (sqrt l) (sqrt h)))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 8.5e-302) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (h <= 8.5d-302) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 8.5e-302) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if h <= 8.5e-302:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (h <= 8.5e-302)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (h <= 8.5e-302)
		tmp = d * -sqrt(((1.0 / l) / h));
	else
		tmp = d / (sqrt(l) * sqrt(h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[h, 8.5e-302], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 18: 41.7% accurate, 3.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;h \leq 1.25 \cdot 10^{-287}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{elif}\;h \leq 1.8 \cdot 10^{+91}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= h 1.25e-287)
   (* d (- (sqrt (/ 1.0 (* l h)))))
   (if (<= h 1.8e+91) (/ d (sqrt (* l h))) (sqrt (* (/ d l) (/ d h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 1.25e-287) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else if (h <= 1.8e+91) {
		tmp = d / sqrt((l * h));
	} else {
		tmp = sqrt(((d / l) * (d / h)));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (h <= 1.25d-287) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else if (h <= 1.8d+91) then
        tmp = d / sqrt((l * h))
    else
        tmp = sqrt(((d / l) * (d / h)))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 1.25e-287) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else if (h <= 1.8e+91) {
		tmp = d / Math.sqrt((l * h));
	} else {
		tmp = Math.sqrt(((d / l) * (d / h)));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if h <= 1.25e-287:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	elif h <= 1.8e+91:
		tmp = d / math.sqrt((l * h))
	else:
		tmp = math.sqrt(((d / l) * (d / h)))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (h <= 1.25e-287)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	elseif (h <= 1.8e+91)
		tmp = Float64(d / sqrt(Float64(l * h)));
	else
		tmp = sqrt(Float64(Float64(d / l) * Float64(d / h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (h <= 1.25e-287)
		tmp = d * -sqrt((1.0 / (l * h)));
	elseif (h <= 1.8e+91)
		tmp = d / sqrt((l * h));
	else
		tmp = sqrt(((d / l) * (d / h)));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[h, 1.25e-287], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[h, 1.8e+91], N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 19: 41.9% accurate, 3.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;h \leq 1.25 \cdot 10^{-287}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{\frac{1}{\ell}}{h}}\right)\\ \mathbf{elif}\;h \leq 3.6 \cdot 10^{+86}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{else}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= h 1.25e-287)
   (* d (- (sqrt (/ (/ 1.0 l) h))))
   (if (<= h 3.6e+86) (/ d (sqrt (* l h))) (sqrt (* (/ d l) (/ d h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 1.25e-287) {
		tmp = d * -sqrt(((1.0 / l) / h));
	} else if (h <= 3.6e+86) {
		tmp = d / sqrt((l * h));
	} else {
		tmp = sqrt(((d / l) * (d / h)));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (h <= 1.25d-287) then
        tmp = d * -sqrt(((1.0d0 / l) / h))
    else if (h <= 3.6d+86) then
        tmp = d / sqrt((l * h))
    else
        tmp = sqrt(((d / l) * (d / h)))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (h <= 1.25e-287) {
		tmp = d * -Math.sqrt(((1.0 / l) / h));
	} else if (h <= 3.6e+86) {
		tmp = d / Math.sqrt((l * h));
	} else {
		tmp = Math.sqrt(((d / l) * (d / h)));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if h <= 1.25e-287:
		tmp = d * -math.sqrt(((1.0 / l) / h))
	elif h <= 3.6e+86:
		tmp = d / math.sqrt((l * h))
	else:
		tmp = math.sqrt(((d / l) * (d / h)))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (h <= 1.25e-287)
		tmp = Float64(d * Float64(-sqrt(Float64(Float64(1.0 / l) / h))));
	elseif (h <= 3.6e+86)
		tmp = Float64(d / sqrt(Float64(l * h)));
	else
		tmp = sqrt(Float64(Float64(d / l) * Float64(d / h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (h <= 1.25e-287)
		tmp = d * -sqrt(((1.0 / l) / h));
	elseif (h <= 3.6e+86)
		tmp = d / sqrt((l * h));
	else
		tmp = sqrt(((d / l) * (d / h)));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[h, 1.25e-287], N[(d * (-N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[h, 3.6e+86], N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 20: 38.1% accurate, 3.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;d \leq -7.6 \cdot 10^{-164}:\\ \;\;\;\;\sqrt{\frac{d}{\ell} \cdot \frac{d}{h}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= d -7.6e-164) (sqrt (* (/ d l) (/ d h))) (/ d (sqrt (* l h)))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -7.6e-164) {
		tmp = sqrt(((d / l) * (d / h)));
	} else {
		tmp = d / sqrt((l * h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (d <= (-7.6d-164)) then
        tmp = sqrt(((d / l) * (d / h)))
    else
        tmp = d / sqrt((l * h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -7.6e-164) {
		tmp = Math.sqrt(((d / l) * (d / h)));
	} else {
		tmp = d / Math.sqrt((l * h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if d <= -7.6e-164:
		tmp = math.sqrt(((d / l) * (d / h)))
	else:
		tmp = d / math.sqrt((l * h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (d <= -7.6e-164)
		tmp = sqrt(Float64(Float64(d / l) * Float64(d / h)));
	else
		tmp = Float64(d / sqrt(Float64(l * h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (d <= -7.6e-164)
		tmp = sqrt(((d / l) * (d / h)));
	else
		tmp = d / sqrt((l * h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[d, -7.6e-164], N[Sqrt[N[(N[(d / l), $MachinePrecision] * N[(d / h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 21: 36.5% accurate, 3.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;d \leq -2.9 \cdot 10^{-161}:\\ \;\;\;\;\sqrt{\frac{d \cdot \frac{d}{h}}{\ell}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= d -2.9e-161) (sqrt (/ (* d (/ d h)) l)) (/ d (sqrt (* l h)))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -2.9e-161) {
		tmp = sqrt(((d * (d / h)) / l));
	} else {
		tmp = d / sqrt((l * h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (d <= (-2.9d-161)) then
        tmp = sqrt(((d * (d / h)) / l))
    else
        tmp = d / sqrt((l * h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -2.9e-161) {
		tmp = Math.sqrt(((d * (d / h)) / l));
	} else {
		tmp = d / Math.sqrt((l * h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if d <= -2.9e-161:
		tmp = math.sqrt(((d * (d / h)) / l))
	else:
		tmp = d / math.sqrt((l * h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (d <= -2.9e-161)
		tmp = sqrt(Float64(Float64(d * Float64(d / h)) / l));
	else
		tmp = Float64(d / sqrt(Float64(l * h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (d <= -2.9e-161)
		tmp = sqrt(((d * (d / h)) / l));
	else
		tmp = d / sqrt((l * h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[d, -2.9e-161], N[Sqrt[N[(N[(d * N[(d / h), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision], N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 22: 36.5% accurate, 3.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;d \leq -1.3 \cdot 10^{-165}:\\ \;\;\;\;\sqrt{\frac{d \cdot \frac{d}{h}}{\ell}}\\ \mathbf{else}:\\ \;\;\;\;d \cdot \sqrt{\frac{1}{\ell \cdot h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= d -1.3e-165) (sqrt (/ (* d (/ d h)) l)) (* d (sqrt (/ 1.0 (* l h))))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -1.3e-165) {
		tmp = sqrt(((d * (d / h)) / l));
	} else {
		tmp = d * sqrt((1.0 / (l * h)));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (d <= (-1.3d-165)) then
        tmp = sqrt(((d * (d / h)) / l))
    else
        tmp = d * sqrt((1.0d0 / (l * h)))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -1.3e-165) {
		tmp = Math.sqrt(((d * (d / h)) / l));
	} else {
		tmp = d * Math.sqrt((1.0 / (l * h)));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if d <= -1.3e-165:
		tmp = math.sqrt(((d * (d / h)) / l))
	else:
		tmp = d * math.sqrt((1.0 / (l * h)))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (d <= -1.3e-165)
		tmp = sqrt(Float64(Float64(d * Float64(d / h)) / l));
	else
		tmp = Float64(d * sqrt(Float64(1.0 / Float64(l * h))));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (d <= -1.3e-165)
		tmp = sqrt(((d * (d / h)) / l));
	else
		tmp = d * sqrt((1.0 / (l * h)));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[d, -1.3e-165], N[Sqrt[N[(N[(d * N[(d / h), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision], N[(d * N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 23: 36.6% accurate, 3.0× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \begin{array}{l} \mathbf{if}\;d \leq -2.5 \cdot 10^{-163}:\\ \;\;\;\;\sqrt{\frac{d \cdot \frac{d}{h}}{\ell}}\\ \mathbf{else}:\\ \;\;\;\;d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}}\\ \end{array} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m)
 :precision binary64
 (if (<= d -2.5e-163) (sqrt (/ (* d (/ d h)) l)) (* d (sqrt (/ (/ 1.0 l) h)))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -2.5e-163) {
		tmp = sqrt(((d * (d / h)) / l));
	} else {
		tmp = d * sqrt(((1.0 / l) / h));
	}
	return tmp;
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    real(8) :: tmp
    if (d <= (-2.5d-163)) then
        tmp = sqrt(((d * (d / h)) / l))
    else
        tmp = d * sqrt(((1.0d0 / l) / h))
    end if
    code = tmp
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	double tmp;
	if (d <= -2.5e-163) {
		tmp = Math.sqrt(((d * (d / h)) / l));
	} else {
		tmp = d * Math.sqrt(((1.0 / l) / h));
	}
	return tmp;
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	tmp = 0
	if d <= -2.5e-163:
		tmp = math.sqrt(((d * (d / h)) / l))
	else:
		tmp = d * math.sqrt(((1.0 / l) / h))
	return tmp

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	tmp = 0.0
	if (d <= -2.5e-163)
		tmp = sqrt(Float64(Float64(d * Float64(d / h)) / l));
	else
		tmp = Float64(d * sqrt(Float64(Float64(1.0 / l) / h)));
	end
	return tmp
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp_2 = code(d, h, l, M_m, D_m)
	tmp = 0.0;
	if (d <= -2.5e-163)
		tmp = sqrt(((d * (d / h)) / l));
	else
		tmp = d * sqrt(((1.0 / l) / h));
	end
	tmp_2 = tmp;
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := If[LessEqual[d, -2.5e-163], N[Sqrt[N[(N[(d * N[(d / h), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision], N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Alternative 24: 27.0% accurate, 3.2× speedup

Math

\[\begin{array}{l} M_m = \left|M\right| \\ D_m = \left|D\right| \\ [d, h, l, M_m, D_m] = \mathsf{sort}([d, h, l, M_m, D_m])\\ \\ \frac{d}{\sqrt{\ell \cdot h}} \end{array} \]

FPCore

M_m = (fabs.f64 M)
D_m = (fabs.f64 D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
(FPCore (d h l M_m D_m) :precision binary64 (/ d (sqrt (* l h))))

C

M_m = fabs(M);
D_m = fabs(D);
assert(d < h && h < l && l < M_m && M_m < D_m);
double code(double d, double h, double l, double M_m, double D_m) {
	return d / sqrt((l * h));
}

Fortran

M_m = abs(M)
D_m = abs(D)
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
real(8) function code(d, h, l, m_m, d_m)
    real(8), intent (in) :: d
    real(8), intent (in) :: h
    real(8), intent (in) :: l
    real(8), intent (in) :: m_m
    real(8), intent (in) :: d_m
    code = d / sqrt((l * h))
end function

Java

M_m = Math.abs(M);
D_m = Math.abs(D);
assert d < h && h < l && l < M_m && M_m < D_m;
public static double code(double d, double h, double l, double M_m, double D_m) {
	return d / Math.sqrt((l * h));
}

Python

M_m = math.fabs(M)
D_m = math.fabs(D)
[d, h, l, M_m, D_m] = sort([d, h, l, M_m, D_m])
def code(d, h, l, M_m, D_m):
	return d / math.sqrt((l * h))

Julia

M_m = abs(M)
D_m = abs(D)
d, h, l, M_m, D_m = sort([d, h, l, M_m, D_m])
function code(d, h, l, M_m, D_m)
	return Float64(d / sqrt(Float64(l * h)))
end

MATLAB

M_m = abs(M);
D_m = abs(D);
d, h, l, M_m, D_m = num2cell(sort([d, h, l, M_m, D_m])){:}
function tmp = code(d, h, l, M_m, D_m)
	tmp = d / sqrt((l * h));
end

Wolfram

M_m = N[Abs[M], $MachinePrecision]
D_m = N[Abs[D], $MachinePrecision]
NOTE: d, h, l, M_m, and D_m should be sorted in increasing order before calling this function.
code[d_, h_, l_, M$95$m_, D$95$m_] := N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

Derivation

&prev;&pcontext;&pcontext2;&ctx;
1. Add Preprocessing

Reproduce

herbie shell --seed 2023343 
(FPCore (d h l M D)
  :name "Henrywood and Agarwal, Equation (12)"
  :precision binary64
  (* (* (pow (/ d h) (/ 1.0 2.0)) (pow (/ d l) (/ 1.0 2.0))) (- 1.0 (* (* (/ 1.0 2.0) (pow (/ (* M D) (* 2.0 d)) 2.0)) (/ h l)))))