Henrywood and Agarwal, Equation (12)

Percentage Accurate: 66.1% → 78.7%

Time: 38.1s

Alternatives: 27

Speedup: 1.4×

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.1% 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: 78.7% 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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := 0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\\ t_2 := \sqrt{\frac{d}{\ell}}\\ \mathbf{if}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_1 \cdot \left(t_2 \cdot \frac{\sqrt{-d}}{\sqrt{-h}}\right)\\ \mathbf{elif}\;d \leq 7.4 \cdot 10^{-121}:\\ \;\;\;\;-0.125 \cdot \frac{1}{\frac{\frac{1}{M_m} \cdot \frac{d}{M_m}}{\frac{\sqrt{h}}{{\ell}^{1.5}} \cdot {D_m}^{2}}}\\ \mathbf{elif}\;d \leq 2.7 \cdot 10^{+153}:\\ \;\;\;\;t_1 \cdot \left(t_2 \cdot \frac{\sqrt{d}}{\sqrt{h}}\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot t_1\\ \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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0))
        (t_2 (sqrt (/ d l))))
   (if (<= d -4e-310)
     (* t_1 (* t_2 (/ (sqrt (- d)) (sqrt (- h)))))
     (if (<= d 7.4e-121)
       (*
        -0.125
        (/
         1.0
         (/
          (* (/ 1.0 M_m) (/ d M_m))
          (* (/ (sqrt h) (pow l 1.5)) (pow D_m 2.0)))))
       (if (<= d 2.7e+153)
         (* t_1 (* t_2 (/ (sqrt d) (sqrt h))))
         (* (* (/ (sqrt d) (sqrt l)) (sqrt (/ d h))) t_1))))))

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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = sqrt((d / l));
	double tmp;
	if (d <= -4e-310) {
		tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-h)));
	} else if (d <= 7.4e-121) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / pow(l, 1.5)) * pow(D_m, 2.0))));
	} else if (d <= 2.7e+153) {
		tmp = t_1 * (t_2 * (sqrt(d) / sqrt(h)));
	} else {
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	}
	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) :: t_2
    real(8) :: tmp
    t_0 = m_m * (d_m / (d * 2.0d0))
    t_1 = (0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0
    t_2 = sqrt((d / l))
    if (d <= (-4d-310)) then
        tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-h)))
    else if (d <= 7.4d-121) then
        tmp = (-0.125d0) * (1.0d0 / (((1.0d0 / m_m) * (d / m_m)) / ((sqrt(h) / (l ** 1.5d0)) * (d_m ** 2.0d0))))
    else if (d <= 2.7d+153) then
        tmp = t_1 * (t_2 * (sqrt(d) / sqrt(h)))
    else
        tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = Math.sqrt((d / l));
	double tmp;
	if (d <= -4e-310) {
		tmp = t_1 * (t_2 * (Math.sqrt(-d) / Math.sqrt(-h)));
	} else if (d <= 7.4e-121) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((Math.sqrt(h) / Math.pow(l, 1.5)) * Math.pow(D_m, 2.0))));
	} else if (d <= 2.7e+153) {
		tmp = t_1 * (t_2 * (Math.sqrt(d) / Math.sqrt(h)));
	} else {
		tmp = ((Math.sqrt(d) / Math.sqrt(l)) * Math.sqrt((d / h))) * t_1;
	}
	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 = M_m * (D_m / (d * 2.0))
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0
	t_2 = math.sqrt((d / l))
	tmp = 0
	if d <= -4e-310:
		tmp = t_1 * (t_2 * (math.sqrt(-d) / math.sqrt(-h)))
	elif d <= 7.4e-121:
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((math.sqrt(h) / math.pow(l, 1.5)) * math.pow(D_m, 2.0))))
	elif d <= 2.7e+153:
		tmp = t_1 * (t_2 * (math.sqrt(d) / math.sqrt(h)))
	else:
		tmp = ((math.sqrt(d) / math.sqrt(l)) * math.sqrt((d / h))) * t_1
	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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0)
	t_2 = sqrt(Float64(d / l))
	tmp = 0.0
	if (d <= -4e-310)
		tmp = Float64(t_1 * Float64(t_2 * Float64(sqrt(Float64(-d)) / sqrt(Float64(-h)))));
	elseif (d <= 7.4e-121)
		tmp = Float64(-0.125 * Float64(1.0 / Float64(Float64(Float64(1.0 / M_m) * Float64(d / M_m)) / Float64(Float64(sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0)))));
	elseif (d <= 2.7e+153)
		tmp = Float64(t_1 * Float64(t_2 * Float64(sqrt(d) / sqrt(h))));
	else
		tmp = Float64(Float64(Float64(sqrt(d) / sqrt(l)) * sqrt(Float64(d / h))) * t_1);
	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 = M_m * (D_m / (d * 2.0));
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	t_2 = sqrt((d / l));
	tmp = 0.0;
	if (d <= -4e-310)
		tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-h)));
	elseif (d <= 7.4e-121)
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0))));
	elseif (d <= 2.7e+153)
		tmp = t_1 * (t_2 * (sqrt(d) / sqrt(h)));
	else
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]}, Block[{t$95$2 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[d, -4e-310], N[(t$95$1 * N[(t$95$2 * N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-h)], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 7.4e-121], N[(-0.125 * N[(1.0 / N[(N[(N[(1.0 / M$95$m), $MachinePrecision] * N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 2.7e+153], N[(t$95$1 * N[(t$95$2 * N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * t$95$1), $MachinePrecision]]]]]]]

Derivation

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

Alternative 2: 83.6% accurate, 0.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} t_0 := h \cdot \frac{{\left(D_m \cdot \frac{M_m}{d \cdot 2}\right)}^{2}}{\frac{\ell}{-0.5}} + 1\\ \mathbf{if}\;\ell \leq -5 \cdot 10^{-310}:\\ \;\;\;\;\sqrt{\frac{d}{\ell}} \cdot \left(\frac{\sqrt{-d}}{\sqrt{-h}} \cdot t_0\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \left(t_0 \cdot \frac{\sqrt{d}}{\sqrt{h}}\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 (+ (* h (/ (pow (* D_m (/ M_m (* d 2.0))) 2.0) (/ l -0.5))) 1.0)))
   (if (<= l -5e-310)
     (* (sqrt (/ d l)) (* (/ (sqrt (- d)) (sqrt (- h))) t_0))
     (* (/ (sqrt d) (sqrt l)) (* t_0 (/ (sqrt d) (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 = (h * (pow((D_m * (M_m / (d * 2.0))), 2.0) / (l / -0.5))) + 1.0;
	double tmp;
	if (l <= -5e-310) {
		tmp = sqrt((d / l)) * ((sqrt(-d) / sqrt(-h)) * t_0);
	} else {
		tmp = (sqrt(d) / sqrt(l)) * (t_0 * (sqrt(d) / 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 = (h * (((d_m * (m_m / (d * 2.0d0))) ** 2.0d0) / (l / (-0.5d0)))) + 1.0d0
    if (l <= (-5d-310)) then
        tmp = sqrt((d / l)) * ((sqrt(-d) / sqrt(-h)) * t_0)
    else
        tmp = (sqrt(d) / sqrt(l)) * (t_0 * (sqrt(d) / 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 = (h * (Math.pow((D_m * (M_m / (d * 2.0))), 2.0) / (l / -0.5))) + 1.0;
	double tmp;
	if (l <= -5e-310) {
		tmp = Math.sqrt((d / l)) * ((Math.sqrt(-d) / Math.sqrt(-h)) * t_0);
	} else {
		tmp = (Math.sqrt(d) / Math.sqrt(l)) * (t_0 * (Math.sqrt(d) / 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 = (h * (math.pow((D_m * (M_m / (d * 2.0))), 2.0) / (l / -0.5))) + 1.0
	tmp = 0
	if l <= -5e-310:
		tmp = math.sqrt((d / l)) * ((math.sqrt(-d) / math.sqrt(-h)) * t_0)
	else:
		tmp = (math.sqrt(d) / math.sqrt(l)) * (t_0 * (math.sqrt(d) / 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(Float64(h * Float64((Float64(D_m * Float64(M_m / Float64(d * 2.0))) ^ 2.0) / Float64(l / -0.5))) + 1.0)
	tmp = 0.0
	if (l <= -5e-310)
		tmp = Float64(sqrt(Float64(d / l)) * Float64(Float64(sqrt(Float64(-d)) / sqrt(Float64(-h))) * t_0));
	else
		tmp = Float64(Float64(sqrt(d) / sqrt(l)) * Float64(t_0 * Float64(sqrt(d) / 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 = (h * (((D_m * (M_m / (d * 2.0))) ^ 2.0) / (l / -0.5))) + 1.0;
	tmp = 0.0;
	if (l <= -5e-310)
		tmp = sqrt((d / l)) * ((sqrt(-d) / sqrt(-h)) * t_0);
	else
		tmp = (sqrt(d) / sqrt(l)) * (t_0 * (sqrt(d) / 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[(N[(h * N[(N[Power[N[(D$95$m * N[(M$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / N[(l / -0.5), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]}, If[LessEqual[l, -5e-310], N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-h)], $MachinePrecision]), $MachinePrecision] * t$95$0), $MachinePrecision]), $MachinePrecision], N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[(t$95$0 * N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

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

Alternative 3: 78.9% 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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := 0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\\ t_2 := \sqrt{\frac{d}{\ell}}\\ \mathbf{if}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_1 \cdot \left(t_2 \cdot \frac{\sqrt{-d}}{\sqrt{-h}}\right)\\ \mathbf{elif}\;d \leq 6.8 \cdot 10^{-69}:\\ \;\;\;\;t_2 \cdot \left(\frac{\sqrt{d}}{\sqrt{h}} \cdot \left(\frac{h}{\ell} \cdot \left(-0.5 \cdot {\left(\frac{D_m}{2} \cdot \frac{M_m}{d}\right)}^{2}\right) + 1\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot t_1\\ \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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0))
        (t_2 (sqrt (/ d l))))
   (if (<= d -4e-310)
     (* t_1 (* t_2 (/ (sqrt (- d)) (sqrt (- h)))))
     (if (<= d 6.8e-69)
       (*
        t_2
        (*
         (/ (sqrt d) (sqrt h))
         (+ (* (/ h l) (* -0.5 (pow (* (/ D_m 2.0) (/ M_m d)) 2.0))) 1.0)))
       (* (* (/ (sqrt d) (sqrt l)) (sqrt (/ d h))) t_1)))))

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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = sqrt((d / l));
	double tmp;
	if (d <= -4e-310) {
		tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-h)));
	} else if (d <= 6.8e-69) {
		tmp = t_2 * ((sqrt(d) / sqrt(h)) * (((h / l) * (-0.5 * pow(((D_m / 2.0) * (M_m / d)), 2.0))) + 1.0));
	} else {
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	}
	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) :: t_2
    real(8) :: tmp
    t_0 = m_m * (d_m / (d * 2.0d0))
    t_1 = (0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0
    t_2 = sqrt((d / l))
    if (d <= (-4d-310)) then
        tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-h)))
    else if (d <= 6.8d-69) then
        tmp = t_2 * ((sqrt(d) / sqrt(h)) * (((h / l) * ((-0.5d0) * (((d_m / 2.0d0) * (m_m / d)) ** 2.0d0))) + 1.0d0))
    else
        tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = Math.sqrt((d / l));
	double tmp;
	if (d <= -4e-310) {
		tmp = t_1 * (t_2 * (Math.sqrt(-d) / Math.sqrt(-h)));
	} else if (d <= 6.8e-69) {
		tmp = t_2 * ((Math.sqrt(d) / Math.sqrt(h)) * (((h / l) * (-0.5 * Math.pow(((D_m / 2.0) * (M_m / d)), 2.0))) + 1.0));
	} else {
		tmp = ((Math.sqrt(d) / Math.sqrt(l)) * Math.sqrt((d / h))) * t_1;
	}
	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 = M_m * (D_m / (d * 2.0))
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0
	t_2 = math.sqrt((d / l))
	tmp = 0
	if d <= -4e-310:
		tmp = t_1 * (t_2 * (math.sqrt(-d) / math.sqrt(-h)))
	elif d <= 6.8e-69:
		tmp = t_2 * ((math.sqrt(d) / math.sqrt(h)) * (((h / l) * (-0.5 * math.pow(((D_m / 2.0) * (M_m / d)), 2.0))) + 1.0))
	else:
		tmp = ((math.sqrt(d) / math.sqrt(l)) * math.sqrt((d / h))) * t_1
	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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0)
	t_2 = sqrt(Float64(d / l))
	tmp = 0.0
	if (d <= -4e-310)
		tmp = Float64(t_1 * Float64(t_2 * Float64(sqrt(Float64(-d)) / sqrt(Float64(-h)))));
	elseif (d <= 6.8e-69)
		tmp = Float64(t_2 * Float64(Float64(sqrt(d) / sqrt(h)) * Float64(Float64(Float64(h / l) * Float64(-0.5 * (Float64(Float64(D_m / 2.0) * Float64(M_m / d)) ^ 2.0))) + 1.0)));
	else
		tmp = Float64(Float64(Float64(sqrt(d) / sqrt(l)) * sqrt(Float64(d / h))) * t_1);
	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 = M_m * (D_m / (d * 2.0));
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	t_2 = sqrt((d / l));
	tmp = 0.0;
	if (d <= -4e-310)
		tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-h)));
	elseif (d <= 6.8e-69)
		tmp = t_2 * ((sqrt(d) / sqrt(h)) * (((h / l) * (-0.5 * (((D_m / 2.0) * (M_m / d)) ^ 2.0))) + 1.0));
	else
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]}, Block[{t$95$2 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[d, -4e-310], N[(t$95$1 * N[(t$95$2 * N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-h)], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 6.8e-69], N[(t$95$2 * N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(N[(N[(h / l), $MachinePrecision] * N[(-0.5 * N[Power[N[(N[(D$95$m / 2.0), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * t$95$1), $MachinePrecision]]]]]]

Derivation

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

Alternative 4: 78.8% 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 := \sqrt{\frac{d}{\ell}}\\ t_1 := M_m \cdot \frac{D_m}{d \cdot 2}\\ \mathbf{if}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_0 \cdot \left(\frac{\sqrt{-d}}{\sqrt{-h}} \cdot \left(h \cdot \frac{{\left(D_m \cdot \frac{M_m}{d \cdot 2}\right)}^{2}}{\frac{\ell}{-0.5}} + 1\right)\right)\\ \mathbf{elif}\;d \leq 3.9 \cdot 10^{-69}:\\ \;\;\;\;t_0 \cdot \left(\frac{\sqrt{d}}{\sqrt{h}} \cdot \left(\frac{h}{\ell} \cdot \left(-0.5 \cdot {\left(\frac{D_m}{2} \cdot \frac{M_m}{d}\right)}^{2}\right) + 1\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot \left(0.5 \cdot \frac{h}{\frac{-1}{t_1} \cdot \frac{\ell}{t_1}} + 1\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 l))) (t_1 (* M_m (/ D_m (* d 2.0)))))
   (if (<= d -4e-310)
     (*
      t_0
      (*
       (/ (sqrt (- d)) (sqrt (- h)))
       (+ (* h (/ (pow (* D_m (/ M_m (* d 2.0))) 2.0) (/ l -0.5))) 1.0)))
     (if (<= d 3.9e-69)
       (*
        t_0
        (*
         (/ (sqrt d) (sqrt h))
         (+ (* (/ h l) (* -0.5 (pow (* (/ D_m 2.0) (/ M_m d)) 2.0))) 1.0)))
       (*
        (* (/ (sqrt d) (sqrt l)) (sqrt (/ d h)))
        (+ (* 0.5 (/ h (* (/ -1.0 t_1) (/ l t_1)))) 1.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 = sqrt((d / l));
	double t_1 = M_m * (D_m / (d * 2.0));
	double tmp;
	if (d <= -4e-310) {
		tmp = t_0 * ((sqrt(-d) / sqrt(-h)) * ((h * (pow((D_m * (M_m / (d * 2.0))), 2.0) / (l / -0.5))) + 1.0));
	} else if (d <= 3.9e-69) {
		tmp = t_0 * ((sqrt(d) / sqrt(h)) * (((h / l) * (-0.5 * pow(((D_m / 2.0) * (M_m / d)), 2.0))) + 1.0));
	} else {
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * ((0.5 * (h / ((-1.0 / t_1) * (l / t_1)))) + 1.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) :: t_1
    real(8) :: tmp
    t_0 = sqrt((d / l))
    t_1 = m_m * (d_m / (d * 2.0d0))
    if (d <= (-4d-310)) then
        tmp = t_0 * ((sqrt(-d) / sqrt(-h)) * ((h * (((d_m * (m_m / (d * 2.0d0))) ** 2.0d0) / (l / (-0.5d0)))) + 1.0d0))
    else if (d <= 3.9d-69) then
        tmp = t_0 * ((sqrt(d) / sqrt(h)) * (((h / l) * ((-0.5d0) * (((d_m / 2.0d0) * (m_m / d)) ** 2.0d0))) + 1.0d0))
    else
        tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * ((0.5d0 * (h / (((-1.0d0) / t_1) * (l / t_1)))) + 1.0d0)
    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 / l));
	double t_1 = M_m * (D_m / (d * 2.0));
	double tmp;
	if (d <= -4e-310) {
		tmp = t_0 * ((Math.sqrt(-d) / Math.sqrt(-h)) * ((h * (Math.pow((D_m * (M_m / (d * 2.0))), 2.0) / (l / -0.5))) + 1.0));
	} else if (d <= 3.9e-69) {
		tmp = t_0 * ((Math.sqrt(d) / Math.sqrt(h)) * (((h / l) * (-0.5 * Math.pow(((D_m / 2.0) * (M_m / d)), 2.0))) + 1.0));
	} else {
		tmp = ((Math.sqrt(d) / Math.sqrt(l)) * Math.sqrt((d / h))) * ((0.5 * (h / ((-1.0 / t_1) * (l / t_1)))) + 1.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 = math.sqrt((d / l))
	t_1 = M_m * (D_m / (d * 2.0))
	tmp = 0
	if d <= -4e-310:
		tmp = t_0 * ((math.sqrt(-d) / math.sqrt(-h)) * ((h * (math.pow((D_m * (M_m / (d * 2.0))), 2.0) / (l / -0.5))) + 1.0))
	elif d <= 3.9e-69:
		tmp = t_0 * ((math.sqrt(d) / math.sqrt(h)) * (((h / l) * (-0.5 * math.pow(((D_m / 2.0) * (M_m / d)), 2.0))) + 1.0))
	else:
		tmp = ((math.sqrt(d) / math.sqrt(l)) * math.sqrt((d / h))) * ((0.5 * (h / ((-1.0 / t_1) * (l / t_1)))) + 1.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 = sqrt(Float64(d / l))
	t_1 = Float64(M_m * Float64(D_m / Float64(d * 2.0)))
	tmp = 0.0
	if (d <= -4e-310)
		tmp = Float64(t_0 * Float64(Float64(sqrt(Float64(-d)) / sqrt(Float64(-h))) * Float64(Float64(h * Float64((Float64(D_m * Float64(M_m / Float64(d * 2.0))) ^ 2.0) / Float64(l / -0.5))) + 1.0)));
	elseif (d <= 3.9e-69)
		tmp = Float64(t_0 * Float64(Float64(sqrt(d) / sqrt(h)) * Float64(Float64(Float64(h / l) * Float64(-0.5 * (Float64(Float64(D_m / 2.0) * Float64(M_m / d)) ^ 2.0))) + 1.0)));
	else
		tmp = Float64(Float64(Float64(sqrt(d) / sqrt(l)) * sqrt(Float64(d / h))) * Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_1) * Float64(l / t_1)))) + 1.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 = sqrt((d / l));
	t_1 = M_m * (D_m / (d * 2.0));
	tmp = 0.0;
	if (d <= -4e-310)
		tmp = t_0 * ((sqrt(-d) / sqrt(-h)) * ((h * (((D_m * (M_m / (d * 2.0))) ^ 2.0) / (l / -0.5))) + 1.0));
	elseif (d <= 3.9e-69)
		tmp = t_0 * ((sqrt(d) / sqrt(h)) * (((h / l) * (-0.5 * (((D_m / 2.0) * (M_m / d)) ^ 2.0))) + 1.0));
	else
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * ((0.5 * (h / ((-1.0 / t_1) * (l / t_1)))) + 1.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[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[d, -4e-310], N[(t$95$0 * N[(N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-h)], $MachinePrecision]), $MachinePrecision] * N[(N[(h * N[(N[Power[N[(D$95$m * N[(M$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / N[(l / -0.5), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 3.9e-69], N[(t$95$0 * N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision] * N[(N[(N[(h / l), $MachinePrecision] * N[(-0.5 * N[Power[N[(N[(D$95$m / 2.0), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$1), $MachinePrecision] * N[(l / t$95$1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

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

Alternative 5: 75.2% 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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := 0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\\ t_2 := \sqrt{\frac{d}{\ell}}\\ \mathbf{if}\;d \leq -4.8 \cdot 10^{-141}:\\ \;\;\;\;t_1 \cdot \left(t_2 \cdot \frac{1}{\sqrt{\frac{h}{d}}}\right)\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;\frac{d \cdot \left(\frac{h}{\ell} \cdot \left(0.5 \cdot {\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}\right) + 1\right)}{\sqrt{\ell \cdot h}}\\ \mathbf{elif}\;d \leq 7.4 \cdot 10^{-121}:\\ \;\;\;\;-0.125 \cdot \frac{1}{\frac{\frac{1}{M_m} \cdot \frac{d}{M_m}}{\frac{\sqrt{h}}{{\ell}^{1.5}} \cdot {D_m}^{2}}}\\ \mathbf{elif}\;d \leq 2.5 \cdot 10^{+153}:\\ \;\;\;\;t_1 \cdot \left(t_2 \cdot \frac{\sqrt{d}}{\sqrt{h}}\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot t_1\\ \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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0))
        (t_2 (sqrt (/ d l))))
   (if (<= d -4.8e-141)
     (* t_1 (* t_2 (/ 1.0 (sqrt (/ h d)))))
     (if (<= d -4e-310)
       (/
        (* d (+ (* (/ h l) (* 0.5 (pow (* 0.5 (* D_m (/ M_m d))) 2.0))) 1.0))
        (sqrt (* l h)))
       (if (<= d 7.4e-121)
         (*
          -0.125
          (/
           1.0
           (/
            (* (/ 1.0 M_m) (/ d M_m))
            (* (/ (sqrt h) (pow l 1.5)) (pow D_m 2.0)))))
         (if (<= d 2.5e+153)
           (* t_1 (* t_2 (/ (sqrt d) (sqrt h))))
           (* (* (/ (sqrt d) (sqrt l)) (sqrt (/ d h))) t_1)))))))

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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = sqrt((d / l));
	double tmp;
	if (d <= -4.8e-141) {
		tmp = t_1 * (t_2 * (1.0 / sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / sqrt((l * h));
	} else if (d <= 7.4e-121) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / pow(l, 1.5)) * pow(D_m, 2.0))));
	} else if (d <= 2.5e+153) {
		tmp = t_1 * (t_2 * (sqrt(d) / sqrt(h)));
	} else {
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	}
	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) :: t_2
    real(8) :: tmp
    t_0 = m_m * (d_m / (d * 2.0d0))
    t_1 = (0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0
    t_2 = sqrt((d / l))
    if (d <= (-4.8d-141)) then
        tmp = t_1 * (t_2 * (1.0d0 / sqrt((h / d))))
    else if (d <= (-4d-310)) then
        tmp = (d * (((h / l) * (0.5d0 * ((0.5d0 * (d_m * (m_m / d))) ** 2.0d0))) + 1.0d0)) / sqrt((l * h))
    else if (d <= 7.4d-121) then
        tmp = (-0.125d0) * (1.0d0 / (((1.0d0 / m_m) * (d / m_m)) / ((sqrt(h) / (l ** 1.5d0)) * (d_m ** 2.0d0))))
    else if (d <= 2.5d+153) then
        tmp = t_1 * (t_2 * (sqrt(d) / sqrt(h)))
    else
        tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = Math.sqrt((d / l));
	double tmp;
	if (d <= -4.8e-141) {
		tmp = t_1 * (t_2 * (1.0 / Math.sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * Math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / Math.sqrt((l * h));
	} else if (d <= 7.4e-121) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((Math.sqrt(h) / Math.pow(l, 1.5)) * Math.pow(D_m, 2.0))));
	} else if (d <= 2.5e+153) {
		tmp = t_1 * (t_2 * (Math.sqrt(d) / Math.sqrt(h)));
	} else {
		tmp = ((Math.sqrt(d) / Math.sqrt(l)) * Math.sqrt((d / h))) * t_1;
	}
	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 = M_m * (D_m / (d * 2.0))
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0
	t_2 = math.sqrt((d / l))
	tmp = 0
	if d <= -4.8e-141:
		tmp = t_1 * (t_2 * (1.0 / math.sqrt((h / d))))
	elif d <= -4e-310:
		tmp = (d * (((h / l) * (0.5 * math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / math.sqrt((l * h))
	elif d <= 7.4e-121:
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((math.sqrt(h) / math.pow(l, 1.5)) * math.pow(D_m, 2.0))))
	elif d <= 2.5e+153:
		tmp = t_1 * (t_2 * (math.sqrt(d) / math.sqrt(h)))
	else:
		tmp = ((math.sqrt(d) / math.sqrt(l)) * math.sqrt((d / h))) * t_1
	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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0)
	t_2 = sqrt(Float64(d / l))
	tmp = 0.0
	if (d <= -4.8e-141)
		tmp = Float64(t_1 * Float64(t_2 * Float64(1.0 / sqrt(Float64(h / d)))));
	elseif (d <= -4e-310)
		tmp = Float64(Float64(d * Float64(Float64(Float64(h / l) * Float64(0.5 * (Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0))) + 1.0)) / sqrt(Float64(l * h)));
	elseif (d <= 7.4e-121)
		tmp = Float64(-0.125 * Float64(1.0 / Float64(Float64(Float64(1.0 / M_m) * Float64(d / M_m)) / Float64(Float64(sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0)))));
	elseif (d <= 2.5e+153)
		tmp = Float64(t_1 * Float64(t_2 * Float64(sqrt(d) / sqrt(h))));
	else
		tmp = Float64(Float64(Float64(sqrt(d) / sqrt(l)) * sqrt(Float64(d / h))) * t_1);
	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 = M_m * (D_m / (d * 2.0));
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	t_2 = sqrt((d / l));
	tmp = 0.0;
	if (d <= -4.8e-141)
		tmp = t_1 * (t_2 * (1.0 / sqrt((h / d))));
	elseif (d <= -4e-310)
		tmp = (d * (((h / l) * (0.5 * ((0.5 * (D_m * (M_m / d))) ^ 2.0))) + 1.0)) / sqrt((l * h));
	elseif (d <= 7.4e-121)
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0))));
	elseif (d <= 2.5e+153)
		tmp = t_1 * (t_2 * (sqrt(d) / sqrt(h)));
	else
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]}, Block[{t$95$2 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[d, -4.8e-141], N[(t$95$1 * N[(t$95$2 * N[(1.0 / N[Sqrt[N[(h / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -4e-310], N[(N[(d * N[(N[(N[(h / l), $MachinePrecision] * N[(0.5 * N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 7.4e-121], N[(-0.125 * N[(1.0 / N[(N[(N[(1.0 / M$95$m), $MachinePrecision] * N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 2.5e+153], N[(t$95$1 * N[(t$95$2 * N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * t$95$1), $MachinePrecision]]]]]]]]

Derivation

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

Alternative 6: 77.8% 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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := 0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\\ t_2 := \sqrt{\frac{d}{h}}\\ \mathbf{if}\;d \leq -1.95 \cdot 10^{-236}:\\ \;\;\;\;t_1 \cdot \left(t_2 \cdot \frac{\sqrt{-d}}{\sqrt{-\ell}}\right)\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;\frac{d \cdot \left(\frac{h}{\ell} \cdot \left(0.5 \cdot {\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}\right) + 1\right)}{\sqrt{\ell \cdot h}}\\ \mathbf{elif}\;d \leq 7.4 \cdot 10^{-121}:\\ \;\;\;\;-0.125 \cdot \frac{1}{\frac{\frac{1}{M_m} \cdot \frac{d}{M_m}}{\frac{\sqrt{h}}{{\ell}^{1.5}} \cdot {D_m}^{2}}}\\ \mathbf{elif}\;d \leq 2.15 \cdot 10^{+153}:\\ \;\;\;\;t_1 \cdot \left(\sqrt{\frac{d}{\ell}} \cdot \frac{\sqrt{d}}{\sqrt{h}}\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\sqrt{d}}{\sqrt{\ell}} \cdot t_2\right) \cdot t_1\\ \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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0))
        (t_2 (sqrt (/ d h))))
   (if (<= d -1.95e-236)
     (* t_1 (* t_2 (/ (sqrt (- d)) (sqrt (- l)))))
     (if (<= d -4e-310)
       (/
        (* d (+ (* (/ h l) (* 0.5 (pow (* 0.5 (* D_m (/ M_m d))) 2.0))) 1.0))
        (sqrt (* l h)))
       (if (<= d 7.4e-121)
         (*
          -0.125
          (/
           1.0
           (/
            (* (/ 1.0 M_m) (/ d M_m))
            (* (/ (sqrt h) (pow l 1.5)) (pow D_m 2.0)))))
         (if (<= d 2.15e+153)
           (* t_1 (* (sqrt (/ d l)) (/ (sqrt d) (sqrt h))))
           (* (* (/ (sqrt d) (sqrt l)) t_2) t_1)))))))

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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = sqrt((d / h));
	double tmp;
	if (d <= -1.95e-236) {
		tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-l)));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / sqrt((l * h));
	} else if (d <= 7.4e-121) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / pow(l, 1.5)) * pow(D_m, 2.0))));
	} else if (d <= 2.15e+153) {
		tmp = t_1 * (sqrt((d / l)) * (sqrt(d) / sqrt(h)));
	} else {
		tmp = ((sqrt(d) / sqrt(l)) * t_2) * t_1;
	}
	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) :: t_2
    real(8) :: tmp
    t_0 = m_m * (d_m / (d * 2.0d0))
    t_1 = (0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0
    t_2 = sqrt((d / h))
    if (d <= (-1.95d-236)) then
        tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-l)))
    else if (d <= (-4d-310)) then
        tmp = (d * (((h / l) * (0.5d0 * ((0.5d0 * (d_m * (m_m / d))) ** 2.0d0))) + 1.0d0)) / sqrt((l * h))
    else if (d <= 7.4d-121) then
        tmp = (-0.125d0) * (1.0d0 / (((1.0d0 / m_m) * (d / m_m)) / ((sqrt(h) / (l ** 1.5d0)) * (d_m ** 2.0d0))))
    else if (d <= 2.15d+153) then
        tmp = t_1 * (sqrt((d / l)) * (sqrt(d) / sqrt(h)))
    else
        tmp = ((sqrt(d) / sqrt(l)) * t_2) * t_1
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double t_2 = Math.sqrt((d / h));
	double tmp;
	if (d <= -1.95e-236) {
		tmp = t_1 * (t_2 * (Math.sqrt(-d) / Math.sqrt(-l)));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * Math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / Math.sqrt((l * h));
	} else if (d <= 7.4e-121) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((Math.sqrt(h) / Math.pow(l, 1.5)) * Math.pow(D_m, 2.0))));
	} else if (d <= 2.15e+153) {
		tmp = t_1 * (Math.sqrt((d / l)) * (Math.sqrt(d) / Math.sqrt(h)));
	} else {
		tmp = ((Math.sqrt(d) / Math.sqrt(l)) * t_2) * t_1;
	}
	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 = M_m * (D_m / (d * 2.0))
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0
	t_2 = math.sqrt((d / h))
	tmp = 0
	if d <= -1.95e-236:
		tmp = t_1 * (t_2 * (math.sqrt(-d) / math.sqrt(-l)))
	elif d <= -4e-310:
		tmp = (d * (((h / l) * (0.5 * math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / math.sqrt((l * h))
	elif d <= 7.4e-121:
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((math.sqrt(h) / math.pow(l, 1.5)) * math.pow(D_m, 2.0))))
	elif d <= 2.15e+153:
		tmp = t_1 * (math.sqrt((d / l)) * (math.sqrt(d) / math.sqrt(h)))
	else:
		tmp = ((math.sqrt(d) / math.sqrt(l)) * t_2) * t_1
	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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0)
	t_2 = sqrt(Float64(d / h))
	tmp = 0.0
	if (d <= -1.95e-236)
		tmp = Float64(t_1 * Float64(t_2 * Float64(sqrt(Float64(-d)) / sqrt(Float64(-l)))));
	elseif (d <= -4e-310)
		tmp = Float64(Float64(d * Float64(Float64(Float64(h / l) * Float64(0.5 * (Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0))) + 1.0)) / sqrt(Float64(l * h)));
	elseif (d <= 7.4e-121)
		tmp = Float64(-0.125 * Float64(1.0 / Float64(Float64(Float64(1.0 / M_m) * Float64(d / M_m)) / Float64(Float64(sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0)))));
	elseif (d <= 2.15e+153)
		tmp = Float64(t_1 * Float64(sqrt(Float64(d / l)) * Float64(sqrt(d) / sqrt(h))));
	else
		tmp = Float64(Float64(Float64(sqrt(d) / sqrt(l)) * t_2) * t_1);
	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 = M_m * (D_m / (d * 2.0));
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	t_2 = sqrt((d / h));
	tmp = 0.0;
	if (d <= -1.95e-236)
		tmp = t_1 * (t_2 * (sqrt(-d) / sqrt(-l)));
	elseif (d <= -4e-310)
		tmp = (d * (((h / l) * (0.5 * ((0.5 * (D_m * (M_m / d))) ^ 2.0))) + 1.0)) / sqrt((l * h));
	elseif (d <= 7.4e-121)
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0))));
	elseif (d <= 2.15e+153)
		tmp = t_1 * (sqrt((d / l)) * (sqrt(d) / sqrt(h)));
	else
		tmp = ((sqrt(d) / sqrt(l)) * t_2) * t_1;
	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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]}, Block[{t$95$2 = N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[d, -1.95e-236], N[(t$95$1 * N[(t$95$2 * N[(N[Sqrt[(-d)], $MachinePrecision] / N[Sqrt[(-l)], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -4e-310], N[(N[(d * N[(N[(N[(h / l), $MachinePrecision] * N[(0.5 * N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 7.4e-121], N[(-0.125 * N[(1.0 / N[(N[(N[(1.0 / M$95$m), $MachinePrecision] * N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 2.15e+153], N[(t$95$1 * N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * t$95$2), $MachinePrecision] * t$95$1), $MachinePrecision]]]]]]]]

Derivation

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

Alternative 7: 74.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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := 0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\\ \mathbf{if}\;d \leq -1.2 \cdot 10^{-142}:\\ \;\;\;\;t_1 \cdot \left(\sqrt{\frac{d}{\ell}} \cdot \frac{1}{\sqrt{\frac{h}{d}}}\right)\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;\frac{d \cdot \left(\frac{h}{\ell} \cdot \left(0.5 \cdot {\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}\right) + 1\right)}{\sqrt{\ell \cdot h}}\\ \mathbf{elif}\;d \leq 2 \cdot 10^{-109}:\\ \;\;\;\;-0.125 \cdot \frac{1}{\frac{\frac{1}{M_m} \cdot \frac{d}{M_m}}{\frac{\sqrt{h}}{{\ell}^{1.5}} \cdot {D_m}^{2}}}\\ \mathbf{else}:\\ \;\;\;\;\left(\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot t_1\\ \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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0)))
   (if (<= d -1.2e-142)
     (* t_1 (* (sqrt (/ d l)) (/ 1.0 (sqrt (/ h d)))))
     (if (<= d -4e-310)
       (/
        (* d (+ (* (/ h l) (* 0.5 (pow (* 0.5 (* D_m (/ M_m d))) 2.0))) 1.0))
        (sqrt (* l h)))
       (if (<= d 2e-109)
         (*
          -0.125
          (/
           1.0
           (/
            (* (/ 1.0 M_m) (/ d M_m))
            (* (/ (sqrt h) (pow l 1.5)) (pow D_m 2.0)))))
         (* (* (/ (sqrt d) (sqrt l)) (sqrt (/ d h))) t_1))))))

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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double tmp;
	if (d <= -1.2e-142) {
		tmp = t_1 * (sqrt((d / l)) * (1.0 / sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / sqrt((l * h));
	} else if (d <= 2e-109) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / pow(l, 1.5)) * pow(D_m, 2.0))));
	} else {
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	}
	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 = m_m * (d_m / (d * 2.0d0))
    t_1 = (0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0
    if (d <= (-1.2d-142)) then
        tmp = t_1 * (sqrt((d / l)) * (1.0d0 / sqrt((h / d))))
    else if (d <= (-4d-310)) then
        tmp = (d * (((h / l) * (0.5d0 * ((0.5d0 * (d_m * (m_m / d))) ** 2.0d0))) + 1.0d0)) / sqrt((l * h))
    else if (d <= 2d-109) then
        tmp = (-0.125d0) * (1.0d0 / (((1.0d0 / m_m) * (d / m_m)) / ((sqrt(h) / (l ** 1.5d0)) * (d_m ** 2.0d0))))
    else
        tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	double tmp;
	if (d <= -1.2e-142) {
		tmp = t_1 * (Math.sqrt((d / l)) * (1.0 / Math.sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * Math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / Math.sqrt((l * h));
	} else if (d <= 2e-109) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((Math.sqrt(h) / Math.pow(l, 1.5)) * Math.pow(D_m, 2.0))));
	} else {
		tmp = ((Math.sqrt(d) / Math.sqrt(l)) * Math.sqrt((d / h))) * t_1;
	}
	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 = M_m * (D_m / (d * 2.0))
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0
	tmp = 0
	if d <= -1.2e-142:
		tmp = t_1 * (math.sqrt((d / l)) * (1.0 / math.sqrt((h / d))))
	elif d <= -4e-310:
		tmp = (d * (((h / l) * (0.5 * math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / math.sqrt((l * h))
	elif d <= 2e-109:
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((math.sqrt(h) / math.pow(l, 1.5)) * math.pow(D_m, 2.0))))
	else:
		tmp = ((math.sqrt(d) / math.sqrt(l)) * math.sqrt((d / h))) * t_1
	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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0)
	tmp = 0.0
	if (d <= -1.2e-142)
		tmp = Float64(t_1 * Float64(sqrt(Float64(d / l)) * Float64(1.0 / sqrt(Float64(h / d)))));
	elseif (d <= -4e-310)
		tmp = Float64(Float64(d * Float64(Float64(Float64(h / l) * Float64(0.5 * (Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0))) + 1.0)) / sqrt(Float64(l * h)));
	elseif (d <= 2e-109)
		tmp = Float64(-0.125 * Float64(1.0 / Float64(Float64(Float64(1.0 / M_m) * Float64(d / M_m)) / Float64(Float64(sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0)))));
	else
		tmp = Float64(Float64(Float64(sqrt(d) / sqrt(l)) * sqrt(Float64(d / h))) * t_1);
	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 = M_m * (D_m / (d * 2.0));
	t_1 = (0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0;
	tmp = 0.0;
	if (d <= -1.2e-142)
		tmp = t_1 * (sqrt((d / l)) * (1.0 / sqrt((h / d))));
	elseif (d <= -4e-310)
		tmp = (d * (((h / l) * (0.5 * ((0.5 * (D_m * (M_m / d))) ^ 2.0))) + 1.0)) / sqrt((l * h));
	elseif (d <= 2e-109)
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0))));
	else
		tmp = ((sqrt(d) / sqrt(l)) * sqrt((d / h))) * t_1;
	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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]}, If[LessEqual[d, -1.2e-142], N[(t$95$1 * N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(1.0 / N[Sqrt[N[(h / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -4e-310], N[(N[(d * N[(N[(N[(h / l), $MachinePrecision] * N[(0.5 * N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 2e-109], N[(-0.125 * N[(1.0 / N[(N[(N[(1.0 / M$95$m), $MachinePrecision] * N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * t$95$1), $MachinePrecision]]]]]]

Derivation

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

Alternative 8: 74.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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ \mathbf{if}\;d \leq -3.5 \cdot 10^{-142}:\\ \;\;\;\;\left(0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\right) \cdot \left(\sqrt{\frac{d}{\ell}} \cdot \frac{1}{\sqrt{\frac{h}{d}}}\right)\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;\frac{d \cdot \left(\frac{h}{\ell} \cdot \left(0.5 \cdot {\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}\right) + 1\right)}{\sqrt{\ell \cdot h}}\\ \mathbf{elif}\;d \leq 2.1 \cdot 10^{-109}:\\ \;\;\;\;-0.125 \cdot \frac{1}{\frac{\frac{1}{M_m} \cdot \frac{d}{M_m}}{\frac{\sqrt{h}}{{\ell}^{1.5}} \cdot {D_m}^{2}}}\\ \mathbf{else}:\\ \;\;\;\;\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \left(\sqrt{\frac{d}{h}} \cdot \left(h \cdot \left(\frac{t_0}{\ell} \cdot \frac{t_0}{-2}\right) + 1\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 (* M_m (/ D_m (* d 2.0)))))
   (if (<= d -3.5e-142)
     (*
      (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0)
      (* (sqrt (/ d l)) (/ 1.0 (sqrt (/ h d)))))
     (if (<= d -4e-310)
       (/
        (* d (+ (* (/ h l) (* 0.5 (pow (* 0.5 (* D_m (/ M_m d))) 2.0))) 1.0))
        (sqrt (* l h)))
       (if (<= d 2.1e-109)
         (*
          -0.125
          (/
           1.0
           (/
            (* (/ 1.0 M_m) (/ d M_m))
            (* (/ (sqrt h) (pow l 1.5)) (pow D_m 2.0)))))
         (*
          (/ (sqrt d) (sqrt l))
          (* (sqrt (/ d h)) (+ (* h (* (/ t_0 l) (/ t_0 -2.0))) 1.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 = M_m * (D_m / (d * 2.0));
	double tmp;
	if (d <= -3.5e-142) {
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (sqrt((d / l)) * (1.0 / sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / sqrt((l * h));
	} else if (d <= 2.1e-109) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / pow(l, 1.5)) * pow(D_m, 2.0))));
	} else {
		tmp = (sqrt(d) / sqrt(l)) * (sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.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 = m_m * (d_m / (d * 2.0d0))
    if (d <= (-3.5d-142)) then
        tmp = ((0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0) * (sqrt((d / l)) * (1.0d0 / sqrt((h / d))))
    else if (d <= (-4d-310)) then
        tmp = (d * (((h / l) * (0.5d0 * ((0.5d0 * (d_m * (m_m / d))) ** 2.0d0))) + 1.0d0)) / sqrt((l * h))
    else if (d <= 2.1d-109) then
        tmp = (-0.125d0) * (1.0d0 / (((1.0d0 / m_m) * (d / m_m)) / ((sqrt(h) / (l ** 1.5d0)) * (d_m ** 2.0d0))))
    else
        tmp = (sqrt(d) / sqrt(l)) * (sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / (-2.0d0)))) + 1.0d0))
    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 = M_m * (D_m / (d * 2.0));
	double tmp;
	if (d <= -3.5e-142) {
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (Math.sqrt((d / l)) * (1.0 / Math.sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * Math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / Math.sqrt((l * h));
	} else if (d <= 2.1e-109) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((Math.sqrt(h) / Math.pow(l, 1.5)) * Math.pow(D_m, 2.0))));
	} else {
		tmp = (Math.sqrt(d) / Math.sqrt(l)) * (Math.sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.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 = M_m * (D_m / (d * 2.0))
	tmp = 0
	if d <= -3.5e-142:
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (math.sqrt((d / l)) * (1.0 / math.sqrt((h / d))))
	elif d <= -4e-310:
		tmp = (d * (((h / l) * (0.5 * math.pow((0.5 * (D_m * (M_m / d))), 2.0))) + 1.0)) / math.sqrt((l * h))
	elif d <= 2.1e-109:
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((math.sqrt(h) / math.pow(l, 1.5)) * math.pow(D_m, 2.0))))
	else:
		tmp = (math.sqrt(d) / math.sqrt(l)) * (math.sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.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(M_m * Float64(D_m / Float64(d * 2.0)))
	tmp = 0.0
	if (d <= -3.5e-142)
		tmp = Float64(Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0) * Float64(sqrt(Float64(d / l)) * Float64(1.0 / sqrt(Float64(h / d)))));
	elseif (d <= -4e-310)
		tmp = Float64(Float64(d * Float64(Float64(Float64(h / l) * Float64(0.5 * (Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0))) + 1.0)) / sqrt(Float64(l * h)));
	elseif (d <= 2.1e-109)
		tmp = Float64(-0.125 * Float64(1.0 / Float64(Float64(Float64(1.0 / M_m) * Float64(d / M_m)) / Float64(Float64(sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0)))));
	else
		tmp = Float64(Float64(sqrt(d) / sqrt(l)) * Float64(sqrt(Float64(d / h)) * Float64(Float64(h * Float64(Float64(t_0 / l) * Float64(t_0 / -2.0))) + 1.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 = M_m * (D_m / (d * 2.0));
	tmp = 0.0;
	if (d <= -3.5e-142)
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (sqrt((d / l)) * (1.0 / sqrt((h / d))));
	elseif (d <= -4e-310)
		tmp = (d * (((h / l) * (0.5 * ((0.5 * (D_m * (M_m / d))) ^ 2.0))) + 1.0)) / sqrt((l * h));
	elseif (d <= 2.1e-109)
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0))));
	else
		tmp = (sqrt(d) / sqrt(l)) * (sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[d, -3.5e-142], N[(N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] * N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(1.0 / N[Sqrt[N[(h / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -4e-310], N[(N[(d * N[(N[(N[(h / l), $MachinePrecision] * N[(0.5 * N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 2.1e-109], N[(-0.125 * N[(1.0 / N[(N[(N[(1.0 / M$95$m), $MachinePrecision] * N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[(N[(h * N[(N[(t$95$0 / l), $MachinePrecision] * N[(t$95$0 / -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

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

Alternative 9: 73.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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := \sqrt{\ell \cdot h}\\ t_2 := \sqrt{\frac{d}{\ell}}\\ t_3 := D_m \cdot \frac{M_m}{d}\\ t_4 := 0.5 \cdot t_3\\ t_5 := {t_4}^{2}\\ \mathbf{if}\;d \leq -8 \cdot 10^{-143}:\\ \;\;\;\;\left(0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\right) \cdot \left(t_2 \cdot \frac{1}{\sqrt{\frac{h}{d}}}\right)\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;\frac{d \cdot \left(\frac{h}{\ell} \cdot \left(0.5 \cdot t_5\right) + 1\right)}{t_1}\\ \mathbf{elif}\;d \leq 2 \cdot 10^{-109}:\\ \;\;\;\;-0.125 \cdot \frac{1}{\frac{\frac{1}{M_m} \cdot \frac{d}{M_m}}{\frac{\sqrt{h}}{{\ell}^{1.5}} \cdot {D_m}^{2}}}\\ \mathbf{elif}\;d \leq 3.7 \cdot 10^{+121}:\\ \;\;\;\;\left(t_2 \cdot \sqrt{\frac{d}{h}}\right) \cdot \left(1 - 0.5 \cdot \left(\frac{0.5}{\frac{\ell}{t_3}} \cdot \left(h \cdot t_4\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{t_1} \cdot \left(h \cdot \left(-0.5 \cdot \frac{t_5}{\ell}\right) + 1\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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (sqrt (* l h)))
        (t_2 (sqrt (/ d l)))
        (t_3 (* D_m (/ M_m d)))
        (t_4 (* 0.5 t_3))
        (t_5 (pow t_4 2.0)))
   (if (<= d -8e-143)
     (*
      (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0)
      (* t_2 (/ 1.0 (sqrt (/ h d)))))
     (if (<= d -4e-310)
       (/ (* d (+ (* (/ h l) (* 0.5 t_5)) 1.0)) t_1)
       (if (<= d 2e-109)
         (*
          -0.125
          (/
           1.0
           (/
            (* (/ 1.0 M_m) (/ d M_m))
            (* (/ (sqrt h) (pow l 1.5)) (pow D_m 2.0)))))
         (if (<= d 3.7e+121)
           (*
            (* t_2 (sqrt (/ d h)))
            (- 1.0 (* 0.5 (* (/ 0.5 (/ l t_3)) (* h t_4)))))
           (* (/ d t_1) (+ (* h (* -0.5 (/ t_5 l))) 1.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 = M_m * (D_m / (d * 2.0));
	double t_1 = sqrt((l * h));
	double t_2 = sqrt((d / l));
	double t_3 = D_m * (M_m / d);
	double t_4 = 0.5 * t_3;
	double t_5 = pow(t_4, 2.0);
	double tmp;
	if (d <= -8e-143) {
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_2 * (1.0 / sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * t_5)) + 1.0)) / t_1;
	} else if (d <= 2e-109) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / pow(l, 1.5)) * pow(D_m, 2.0))));
	} else if (d <= 3.7e+121) {
		tmp = (t_2 * sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_3)) * (h * t_4))));
	} else {
		tmp = (d / t_1) * ((h * (-0.5 * (t_5 / l))) + 1.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) :: t_1
    real(8) :: t_2
    real(8) :: t_3
    real(8) :: t_4
    real(8) :: t_5
    real(8) :: tmp
    t_0 = m_m * (d_m / (d * 2.0d0))
    t_1 = sqrt((l * h))
    t_2 = sqrt((d / l))
    t_3 = d_m * (m_m / d)
    t_4 = 0.5d0 * t_3
    t_5 = t_4 ** 2.0d0
    if (d <= (-8d-143)) then
        tmp = ((0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0) * (t_2 * (1.0d0 / sqrt((h / d))))
    else if (d <= (-4d-310)) then
        tmp = (d * (((h / l) * (0.5d0 * t_5)) + 1.0d0)) / t_1
    else if (d <= 2d-109) then
        tmp = (-0.125d0) * (1.0d0 / (((1.0d0 / m_m) * (d / m_m)) / ((sqrt(h) / (l ** 1.5d0)) * (d_m ** 2.0d0))))
    else if (d <= 3.7d+121) then
        tmp = (t_2 * sqrt((d / h))) * (1.0d0 - (0.5d0 * ((0.5d0 / (l / t_3)) * (h * t_4))))
    else
        tmp = (d / t_1) * ((h * ((-0.5d0) * (t_5 / l))) + 1.0d0)
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = Math.sqrt((l * h));
	double t_2 = Math.sqrt((d / l));
	double t_3 = D_m * (M_m / d);
	double t_4 = 0.5 * t_3;
	double t_5 = Math.pow(t_4, 2.0);
	double tmp;
	if (d <= -8e-143) {
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_2 * (1.0 / Math.sqrt((h / d))));
	} else if (d <= -4e-310) {
		tmp = (d * (((h / l) * (0.5 * t_5)) + 1.0)) / t_1;
	} else if (d <= 2e-109) {
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((Math.sqrt(h) / Math.pow(l, 1.5)) * Math.pow(D_m, 2.0))));
	} else if (d <= 3.7e+121) {
		tmp = (t_2 * Math.sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_3)) * (h * t_4))));
	} else {
		tmp = (d / t_1) * ((h * (-0.5 * (t_5 / l))) + 1.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 = M_m * (D_m / (d * 2.0))
	t_1 = math.sqrt((l * h))
	t_2 = math.sqrt((d / l))
	t_3 = D_m * (M_m / d)
	t_4 = 0.5 * t_3
	t_5 = math.pow(t_4, 2.0)
	tmp = 0
	if d <= -8e-143:
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_2 * (1.0 / math.sqrt((h / d))))
	elif d <= -4e-310:
		tmp = (d * (((h / l) * (0.5 * t_5)) + 1.0)) / t_1
	elif d <= 2e-109:
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((math.sqrt(h) / math.pow(l, 1.5)) * math.pow(D_m, 2.0))))
	elif d <= 3.7e+121:
		tmp = (t_2 * math.sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_3)) * (h * t_4))))
	else:
		tmp = (d / t_1) * ((h * (-0.5 * (t_5 / l))) + 1.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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = sqrt(Float64(l * h))
	t_2 = sqrt(Float64(d / l))
	t_3 = Float64(D_m * Float64(M_m / d))
	t_4 = Float64(0.5 * t_3)
	t_5 = t_4 ^ 2.0
	tmp = 0.0
	if (d <= -8e-143)
		tmp = Float64(Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0) * Float64(t_2 * Float64(1.0 / sqrt(Float64(h / d)))));
	elseif (d <= -4e-310)
		tmp = Float64(Float64(d * Float64(Float64(Float64(h / l) * Float64(0.5 * t_5)) + 1.0)) / t_1);
	elseif (d <= 2e-109)
		tmp = Float64(-0.125 * Float64(1.0 / Float64(Float64(Float64(1.0 / M_m) * Float64(d / M_m)) / Float64(Float64(sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0)))));
	elseif (d <= 3.7e+121)
		tmp = Float64(Float64(t_2 * sqrt(Float64(d / h))) * Float64(1.0 - Float64(0.5 * Float64(Float64(0.5 / Float64(l / t_3)) * Float64(h * t_4)))));
	else
		tmp = Float64(Float64(d / t_1) * Float64(Float64(h * Float64(-0.5 * Float64(t_5 / l))) + 1.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 = M_m * (D_m / (d * 2.0));
	t_1 = sqrt((l * h));
	t_2 = sqrt((d / l));
	t_3 = D_m * (M_m / d);
	t_4 = 0.5 * t_3;
	t_5 = t_4 ^ 2.0;
	tmp = 0.0;
	if (d <= -8e-143)
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_2 * (1.0 / sqrt((h / d))));
	elseif (d <= -4e-310)
		tmp = (d * (((h / l) * (0.5 * t_5)) + 1.0)) / t_1;
	elseif (d <= 2e-109)
		tmp = -0.125 * (1.0 / (((1.0 / M_m) * (d / M_m)) / ((sqrt(h) / (l ^ 1.5)) * (D_m ^ 2.0))));
	elseif (d <= 3.7e+121)
		tmp = (t_2 * sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_3)) * (h * t_4))));
	else
		tmp = (d / t_1) * ((h * (-0.5 * (t_5 / l))) + 1.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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$2 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$3 = N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$4 = N[(0.5 * t$95$3), $MachinePrecision]}, Block[{t$95$5 = N[Power[t$95$4, 2.0], $MachinePrecision]}, If[LessEqual[d, -8e-143], N[(N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] * N[(t$95$2 * N[(1.0 / N[Sqrt[N[(h / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -4e-310], N[(N[(d * N[(N[(N[(h / l), $MachinePrecision] * N[(0.5 * t$95$5), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision] / t$95$1), $MachinePrecision], If[LessEqual[d, 2e-109], N[(-0.125 * N[(1.0 / N[(N[(N[(1.0 / M$95$m), $MachinePrecision] * N[(d / M$95$m), $MachinePrecision]), $MachinePrecision] / N[(N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision] * N[Power[D$95$m, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, 3.7e+121], N[(N[(t$95$2 * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(0.5 / N[(l / t$95$3), $MachinePrecision]), $MachinePrecision] * N[(h * t$95$4), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(d / t$95$1), $MachinePrecision] * N[(N[(h * N[(-0.5 * N[(t$95$5 / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]]]]]]]]]]]

Derivation

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

Alternative 10: 71.9% accurate, 1.4× 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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ t_1 := \sqrt{\frac{d}{\ell}}\\ t_2 := D_m \cdot \frac{M_m}{d}\\ t_3 := 0.5 \cdot t_2\\ \mathbf{if}\;\ell \leq 10^{-273}:\\ \;\;\;\;\left(0.5 \cdot \frac{h}{\frac{-1}{t_0} \cdot \frac{\ell}{t_0}} + 1\right) \cdot \left(t_1 \cdot \frac{1}{\sqrt{\frac{h}{d}}}\right)\\ \mathbf{elif}\;\ell \leq 1.4 \cdot 10^{+133}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(h \cdot \left(-0.5 \cdot \frac{{t_3}^{2}}{\ell}\right) + 1\right)\\ \mathbf{else}:\\ \;\;\;\;\left(t_1 \cdot \sqrt{\frac{d}{h}}\right) \cdot \left(1 - 0.5 \cdot \left(\frac{0.5}{\frac{\ell}{t_2}} \cdot \left(h \cdot t_3\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 (* M_m (/ D_m (* d 2.0))))
        (t_1 (sqrt (/ d l)))
        (t_2 (* D_m (/ M_m d)))
        (t_3 (* 0.5 t_2)))
   (if (<= l 1e-273)
     (*
      (+ (* 0.5 (/ h (* (/ -1.0 t_0) (/ l t_0)))) 1.0)
      (* t_1 (/ 1.0 (sqrt (/ h d)))))
     (if (<= l 1.4e+133)
       (* (/ d (sqrt (* l h))) (+ (* h (* -0.5 (/ (pow t_3 2.0) l))) 1.0))
       (*
        (* t_1 (sqrt (/ d h)))
        (- 1.0 (* 0.5 (* (/ 0.5 (/ l t_2)) (* h t_3)))))))))

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 = M_m * (D_m / (d * 2.0));
	double t_1 = sqrt((d / l));
	double t_2 = D_m * (M_m / d);
	double t_3 = 0.5 * t_2;
	double tmp;
	if (l <= 1e-273) {
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_1 * (1.0 / sqrt((h / d))));
	} else if (l <= 1.4e+133) {
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * (pow(t_3, 2.0) / l))) + 1.0);
	} else {
		tmp = (t_1 * sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_2)) * (h * t_3))));
	}
	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) :: t_2
    real(8) :: t_3
    real(8) :: tmp
    t_0 = m_m * (d_m / (d * 2.0d0))
    t_1 = sqrt((d / l))
    t_2 = d_m * (m_m / d)
    t_3 = 0.5d0 * t_2
    if (l <= 1d-273) then
        tmp = ((0.5d0 * (h / (((-1.0d0) / t_0) * (l / t_0)))) + 1.0d0) * (t_1 * (1.0d0 / sqrt((h / d))))
    else if (l <= 1.4d+133) then
        tmp = (d / sqrt((l * h))) * ((h * ((-0.5d0) * ((t_3 ** 2.0d0) / l))) + 1.0d0)
    else
        tmp = (t_1 * sqrt((d / h))) * (1.0d0 - (0.5d0 * ((0.5d0 / (l / t_2)) * (h * t_3))))
    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 = M_m * (D_m / (d * 2.0));
	double t_1 = Math.sqrt((d / l));
	double t_2 = D_m * (M_m / d);
	double t_3 = 0.5 * t_2;
	double tmp;
	if (l <= 1e-273) {
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_1 * (1.0 / Math.sqrt((h / d))));
	} else if (l <= 1.4e+133) {
		tmp = (d / Math.sqrt((l * h))) * ((h * (-0.5 * (Math.pow(t_3, 2.0) / l))) + 1.0);
	} else {
		tmp = (t_1 * Math.sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_2)) * (h * t_3))));
	}
	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 = M_m * (D_m / (d * 2.0))
	t_1 = math.sqrt((d / l))
	t_2 = D_m * (M_m / d)
	t_3 = 0.5 * t_2
	tmp = 0
	if l <= 1e-273:
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_1 * (1.0 / math.sqrt((h / d))))
	elif l <= 1.4e+133:
		tmp = (d / math.sqrt((l * h))) * ((h * (-0.5 * (math.pow(t_3, 2.0) / l))) + 1.0)
	else:
		tmp = (t_1 * math.sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_2)) * (h * t_3))))
	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(M_m * Float64(D_m / Float64(d * 2.0)))
	t_1 = sqrt(Float64(d / l))
	t_2 = Float64(D_m * Float64(M_m / d))
	t_3 = Float64(0.5 * t_2)
	tmp = 0.0
	if (l <= 1e-273)
		tmp = Float64(Float64(Float64(0.5 * Float64(h / Float64(Float64(-1.0 / t_0) * Float64(l / t_0)))) + 1.0) * Float64(t_1 * Float64(1.0 / sqrt(Float64(h / d)))));
	elseif (l <= 1.4e+133)
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(Float64(h * Float64(-0.5 * Float64((t_3 ^ 2.0) / l))) + 1.0));
	else
		tmp = Float64(Float64(t_1 * sqrt(Float64(d / h))) * Float64(1.0 - Float64(0.5 * Float64(Float64(0.5 / Float64(l / t_2)) * Float64(h * t_3)))));
	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 = M_m * (D_m / (d * 2.0));
	t_1 = sqrt((d / l));
	t_2 = D_m * (M_m / d);
	t_3 = 0.5 * t_2;
	tmp = 0.0;
	if (l <= 1e-273)
		tmp = ((0.5 * (h / ((-1.0 / t_0) * (l / t_0)))) + 1.0) * (t_1 * (1.0 / sqrt((h / d))));
	elseif (l <= 1.4e+133)
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * ((t_3 ^ 2.0) / l))) + 1.0);
	else
		tmp = (t_1 * sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_2)) * (h * t_3))));
	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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$2 = N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$3 = N[(0.5 * t$95$2), $MachinePrecision]}, If[LessEqual[l, 1e-273], N[(N[(N[(0.5 * N[(h / N[(N[(-1.0 / t$95$0), $MachinePrecision] * N[(l / t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision] * N[(t$95$1 * N[(1.0 / N[Sqrt[N[(h / d), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 1.4e+133], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(h * N[(-0.5 * N[(N[Power[t$95$3, 2.0], $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(N[(t$95$1 * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(0.5 / N[(l / t$95$2), $MachinePrecision]), $MachinePrecision] * N[(h * t$95$3), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]]

Derivation

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

Alternative 11: 71.8% accurate, 1.4× 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 := M_m \cdot \frac{D_m}{d \cdot 2}\\ \mathbf{if}\;\ell \leq 5.8 \cdot 10^{-277} \lor \neg \left(\ell \leq 1.75 \cdot 10^{+137}\right):\\ \;\;\;\;\sqrt{\frac{d}{\ell}} \cdot \left(\sqrt{\frac{d}{h}} \cdot \left(h \cdot \left(\frac{t_0}{\ell} \cdot \frac{t_0}{-2}\right) + 1\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(h \cdot \left(-0.5 \cdot \frac{{\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}}{\ell}\right) + 1\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 (* M_m (/ D_m (* d 2.0)))))
   (if (or (<= l 5.8e-277) (not (<= l 1.75e+137)))
     (*
      (sqrt (/ d l))
      (* (sqrt (/ d h)) (+ (* h (* (/ t_0 l) (/ t_0 -2.0))) 1.0)))
     (*
      (/ d (sqrt (* l h)))
      (+ (* h (* -0.5 (/ (pow (* 0.5 (* D_m (/ M_m d))) 2.0) l))) 1.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 = M_m * (D_m / (d * 2.0));
	double tmp;
	if ((l <= 5.8e-277) || !(l <= 1.75e+137)) {
		tmp = sqrt((d / l)) * (sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.0));
	} else {
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * (pow((0.5 * (D_m * (M_m / d))), 2.0) / l))) + 1.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 = m_m * (d_m / (d * 2.0d0))
    if ((l <= 5.8d-277) .or. (.not. (l <= 1.75d+137))) then
        tmp = sqrt((d / l)) * (sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / (-2.0d0)))) + 1.0d0))
    else
        tmp = (d / sqrt((l * h))) * ((h * ((-0.5d0) * (((0.5d0 * (d_m * (m_m / d))) ** 2.0d0) / l))) + 1.0d0)
    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 = M_m * (D_m / (d * 2.0));
	double tmp;
	if ((l <= 5.8e-277) || !(l <= 1.75e+137)) {
		tmp = Math.sqrt((d / l)) * (Math.sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.0));
	} else {
		tmp = (d / Math.sqrt((l * h))) * ((h * (-0.5 * (Math.pow((0.5 * (D_m * (M_m / d))), 2.0) / l))) + 1.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 = M_m * (D_m / (d * 2.0))
	tmp = 0
	if (l <= 5.8e-277) or not (l <= 1.75e+137):
		tmp = math.sqrt((d / l)) * (math.sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.0))
	else:
		tmp = (d / math.sqrt((l * h))) * ((h * (-0.5 * (math.pow((0.5 * (D_m * (M_m / d))), 2.0) / l))) + 1.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(M_m * Float64(D_m / Float64(d * 2.0)))
	tmp = 0.0
	if ((l <= 5.8e-277) || !(l <= 1.75e+137))
		tmp = Float64(sqrt(Float64(d / l)) * Float64(sqrt(Float64(d / h)) * Float64(Float64(h * Float64(Float64(t_0 / l) * Float64(t_0 / -2.0))) + 1.0)));
	else
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(Float64(h * Float64(-0.5 * Float64((Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0) / l))) + 1.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 = M_m * (D_m / (d * 2.0));
	tmp = 0.0;
	if ((l <= 5.8e-277) || ~((l <= 1.75e+137)))
		tmp = sqrt((d / l)) * (sqrt((d / h)) * ((h * ((t_0 / l) * (t_0 / -2.0))) + 1.0));
	else
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * (((0.5 * (D_m * (M_m / d))) ^ 2.0) / l))) + 1.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[(M$95$m * N[(D$95$m / N[(d * 2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[Or[LessEqual[l, 5.8e-277], N[Not[LessEqual[l, 1.75e+137]], $MachinePrecision]], N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[(N[(h * N[(N[(t$95$0 / l), $MachinePrecision] * N[(t$95$0 / -2.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(h * N[(-0.5 * N[(N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]]]

Derivation

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

Alternative 12: 71.0% accurate, 1.4× 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(D_m \cdot \frac{M_m}{d}\right)\\ \mathbf{if}\;\ell \leq 5 \cdot 10^{-277} \lor \neg \left(\ell \leq 5.2 \cdot 10^{+134}\right):\\ \;\;\;\;\left(\sqrt{\frac{d}{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot \left(1 - 0.5 \cdot \left(t_0 \cdot \left(\frac{h}{\ell} \cdot t_0\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(h \cdot \left(-0.5 \cdot \frac{{t_0}^{2}}{\ell}\right) + 1\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 (* D_m (/ M_m d)))))
   (if (or (<= l 5e-277) (not (<= l 5.2e+134)))
     (*
      (* (sqrt (/ d l)) (sqrt (/ d h)))
      (- 1.0 (* 0.5 (* t_0 (* (/ h l) t_0)))))
     (* (/ d (sqrt (* l h))) (+ (* h (* -0.5 (/ (pow t_0 2.0) l))) 1.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 * (D_m * (M_m / d));
	double tmp;
	if ((l <= 5e-277) || !(l <= 5.2e+134)) {
		tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0 - (0.5 * (t_0 * ((h / l) * t_0))));
	} else {
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * (pow(t_0, 2.0) / l))) + 1.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 * (m_m / d))
    if ((l <= 5d-277) .or. (.not. (l <= 5.2d+134))) then
        tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0d0 - (0.5d0 * (t_0 * ((h / l) * t_0))))
    else
        tmp = (d / sqrt((l * h))) * ((h * ((-0.5d0) * ((t_0 ** 2.0d0) / l))) + 1.0d0)
    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 * (D_m * (M_m / d));
	double tmp;
	if ((l <= 5e-277) || !(l <= 5.2e+134)) {
		tmp = (Math.sqrt((d / l)) * Math.sqrt((d / h))) * (1.0 - (0.5 * (t_0 * ((h / l) * t_0))));
	} else {
		tmp = (d / Math.sqrt((l * h))) * ((h * (-0.5 * (Math.pow(t_0, 2.0) / l))) + 1.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 * (D_m * (M_m / d))
	tmp = 0
	if (l <= 5e-277) or not (l <= 5.2e+134):
		tmp = (math.sqrt((d / l)) * math.sqrt((d / h))) * (1.0 - (0.5 * (t_0 * ((h / l) * t_0))))
	else:
		tmp = (d / math.sqrt((l * h))) * ((h * (-0.5 * (math.pow(t_0, 2.0) / l))) + 1.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(D_m * Float64(M_m / d)))
	tmp = 0.0
	if ((l <= 5e-277) || !(l <= 5.2e+134))
		tmp = Float64(Float64(sqrt(Float64(d / l)) * sqrt(Float64(d / h))) * Float64(1.0 - Float64(0.5 * Float64(t_0 * Float64(Float64(h / l) * t_0)))));
	else
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(Float64(h * Float64(-0.5 * Float64((t_0 ^ 2.0) / l))) + 1.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 * (M_m / d));
	tmp = 0.0;
	if ((l <= 5e-277) || ~((l <= 5.2e+134)))
		tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0 - (0.5 * (t_0 * ((h / l) * t_0))));
	else
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * ((t_0 ^ 2.0) / l))) + 1.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[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[Or[LessEqual[l, 5e-277], N[Not[LessEqual[l, 5.2e+134]], $MachinePrecision]], N[(N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(t$95$0 * N[(N[(h / l), $MachinePrecision] * t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(h * N[(-0.5 * N[(N[Power[t$95$0, 2.0], $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]]]

Derivation

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

Alternative 13: 72.7% accurate, 1.4× 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 := D_m \cdot \frac{M_m}{d}\\ t_1 := 0.5 \cdot t_0\\ \mathbf{if}\;\ell \leq 5.2 \cdot 10^{-277} \lor \neg \left(\ell \leq 4 \cdot 10^{+132}\right):\\ \;\;\;\;\left(\sqrt{\frac{d}{\ell}} \cdot \sqrt{\frac{d}{h}}\right) \cdot \left(1 - 0.5 \cdot \left(\frac{0.5}{\frac{\ell}{t_0}} \cdot \left(h \cdot t_1\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(h \cdot \left(-0.5 \cdot \frac{{t_1}^{2}}{\ell}\right) + 1\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 (* D_m (/ M_m d))) (t_1 (* 0.5 t_0)))
   (if (or (<= l 5.2e-277) (not (<= l 4e+132)))
     (*
      (* (sqrt (/ d l)) (sqrt (/ d h)))
      (- 1.0 (* 0.5 (* (/ 0.5 (/ l t_0)) (* h t_1)))))
     (* (/ d (sqrt (* l h))) (+ (* h (* -0.5 (/ (pow t_1 2.0) l))) 1.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 = D_m * (M_m / d);
	double t_1 = 0.5 * t_0;
	double tmp;
	if ((l <= 5.2e-277) || !(l <= 4e+132)) {
		tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_0)) * (h * t_1))));
	} else {
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * (pow(t_1, 2.0) / l))) + 1.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) :: t_1
    real(8) :: tmp
    t_0 = d_m * (m_m / d)
    t_1 = 0.5d0 * t_0
    if ((l <= 5.2d-277) .or. (.not. (l <= 4d+132))) then
        tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0d0 - (0.5d0 * ((0.5d0 / (l / t_0)) * (h * t_1))))
    else
        tmp = (d / sqrt((l * h))) * ((h * ((-0.5d0) * ((t_1 ** 2.0d0) / l))) + 1.0d0)
    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 = D_m * (M_m / d);
	double t_1 = 0.5 * t_0;
	double tmp;
	if ((l <= 5.2e-277) || !(l <= 4e+132)) {
		tmp = (Math.sqrt((d / l)) * Math.sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_0)) * (h * t_1))));
	} else {
		tmp = (d / Math.sqrt((l * h))) * ((h * (-0.5 * (Math.pow(t_1, 2.0) / l))) + 1.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 = D_m * (M_m / d)
	t_1 = 0.5 * t_0
	tmp = 0
	if (l <= 5.2e-277) or not (l <= 4e+132):
		tmp = (math.sqrt((d / l)) * math.sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_0)) * (h * t_1))))
	else:
		tmp = (d / math.sqrt((l * h))) * ((h * (-0.5 * (math.pow(t_1, 2.0) / l))) + 1.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(D_m * Float64(M_m / d))
	t_1 = Float64(0.5 * t_0)
	tmp = 0.0
	if ((l <= 5.2e-277) || !(l <= 4e+132))
		tmp = Float64(Float64(sqrt(Float64(d / l)) * sqrt(Float64(d / h))) * Float64(1.0 - Float64(0.5 * Float64(Float64(0.5 / Float64(l / t_0)) * Float64(h * t_1)))));
	else
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(Float64(h * Float64(-0.5 * Float64((t_1 ^ 2.0) / l))) + 1.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 = D_m * (M_m / d);
	t_1 = 0.5 * t_0;
	tmp = 0.0;
	if ((l <= 5.2e-277) || ~((l <= 4e+132)))
		tmp = (sqrt((d / l)) * sqrt((d / h))) * (1.0 - (0.5 * ((0.5 / (l / t_0)) * (h * t_1))));
	else
		tmp = (d / sqrt((l * h))) * ((h * (-0.5 * ((t_1 ^ 2.0) / l))) + 1.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[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(0.5 * t$95$0), $MachinePrecision]}, If[Or[LessEqual[l, 5.2e-277], N[Not[LessEqual[l, 4e+132]], $MachinePrecision]], 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[(0.5 / N[(l / t$95$0), $MachinePrecision]), $MachinePrecision] * N[(h * t$95$1), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(h * N[(-0.5 * N[(N[Power[t$95$1, 2.0], $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]]]]

Derivation

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

Alternative 14: 63.0% 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 := \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{if}\;\ell \leq -1.75 \cdot 10^{-15}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{elif}\;\ell \leq -6.4 \cdot 10^{-301}:\\ \;\;\;\;t_0 \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;\ell \leq 4.4 \cdot 10^{+160}:\\ \;\;\;\;t_0 \cdot \left(h \cdot \left(-0.5 \cdot \frac{{\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}}{\ell}\right) + 1\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \frac{\sqrt{\frac{1}{\ell}}}{\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 (/ d (sqrt (* l h)))))
   (if (<= l -1.75e-15)
     (* d (- (sqrt (/ 1.0 (* l h)))))
     (if (<= l -6.4e-301)
       (*
        t_0
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= l 4.4e+160)
         (*
          t_0
          (+ (* h (* -0.5 (/ (pow (* 0.5 (* D_m (/ M_m d))) 2.0) l))) 1.0))
         (* d (/ (sqrt (/ 1.0 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 = d / sqrt((l * h));
	double tmp;
	if (l <= -1.75e-15) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else if (l <= -6.4e-301) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (l <= 4.4e+160) {
		tmp = t_0 * ((h * (-0.5 * (pow((0.5 * (D_m * (M_m / d))), 2.0) / l))) + 1.0);
	} else {
		tmp = d * (sqrt((1.0 / 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 = d / sqrt((l * h))
    if (l <= (-1.75d-15)) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else if (l <= (-6.4d-301)) then
        tmp = t_0 * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (l <= 4.4d+160) then
        tmp = t_0 * ((h * ((-0.5d0) * (((0.5d0 * (d_m * (m_m / d))) ** 2.0d0) / l))) + 1.0d0)
    else
        tmp = d * (sqrt((1.0d0 / 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 = d / Math.sqrt((l * h));
	double tmp;
	if (l <= -1.75e-15) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else if (l <= -6.4e-301) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (l <= 4.4e+160) {
		tmp = t_0 * ((h * (-0.5 * (Math.pow((0.5 * (D_m * (M_m / d))), 2.0) / l))) + 1.0);
	} else {
		tmp = d * (Math.sqrt((1.0 / 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 = d / math.sqrt((l * h))
	tmp = 0
	if l <= -1.75e-15:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	elif l <= -6.4e-301:
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif l <= 4.4e+160:
		tmp = t_0 * ((h * (-0.5 * (math.pow((0.5 * (D_m * (M_m / d))), 2.0) / l))) + 1.0)
	else:
		tmp = d * (math.sqrt((1.0 / 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(d / sqrt(Float64(l * h)))
	tmp = 0.0
	if (l <= -1.75e-15)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	elseif (l <= -6.4e-301)
		tmp = Float64(t_0 * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (l <= 4.4e+160)
		tmp = Float64(t_0 * Float64(Float64(h * Float64(-0.5 * Float64((Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0) / l))) + 1.0));
	else
		tmp = Float64(d * Float64(sqrt(Float64(1.0 / 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 = d / sqrt((l * h));
	tmp = 0.0;
	if (l <= -1.75e-15)
		tmp = d * -sqrt((1.0 / (l * h)));
	elseif (l <= -6.4e-301)
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (l <= 4.4e+160)
		tmp = t_0 * ((h * (-0.5 * (((0.5 * (D_m * (M_m / d))) ^ 2.0) / l))) + 1.0);
	else
		tmp = d * (sqrt((1.0 / 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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[l, -1.75e-15], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[l, -6.4e-301], N[(t$95$0 * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 4.4e+160], N[(t$95$0 * N[(N[(h * N[(-0.5 * N[(N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Sqrt[N[(1.0 / l), $MachinePrecision]], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

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

Alternative 15: 63.6% 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 := \frac{d}{\sqrt{\ell \cdot h}}\\ t_1 := {\left(0.5 \cdot \left(D_m \cdot \frac{M_m}{d}\right)\right)}^{2}\\ \mathbf{if}\;\ell \leq -7.5 \cdot 10^{-16}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{elif}\;\ell \leq -6.4 \cdot 10^{-301}:\\ \;\;\;\;t_0 \cdot \left(\frac{t_1 \cdot \left(h \cdot 0.5\right)}{\ell} + 1\right)\\ \mathbf{elif}\;\ell \leq 4.2 \cdot 10^{+160}:\\ \;\;\;\;t_0 \cdot \left(h \cdot \left(-0.5 \cdot \frac{t_1}{\ell}\right) + 1\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \frac{\sqrt{\frac{1}{\ell}}}{\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 (/ d (sqrt (* l h)))) (t_1 (pow (* 0.5 (* D_m (/ M_m d))) 2.0)))
   (if (<= l -7.5e-16)
     (* d (- (sqrt (/ 1.0 (* l h)))))
     (if (<= l -6.4e-301)
       (* t_0 (+ (/ (* t_1 (* h 0.5)) l) 1.0))
       (if (<= l 4.2e+160)
         (* t_0 (+ (* h (* -0.5 (/ t_1 l))) 1.0))
         (* d (/ (sqrt (/ 1.0 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 = d / sqrt((l * h));
	double t_1 = pow((0.5 * (D_m * (M_m / d))), 2.0);
	double tmp;
	if (l <= -7.5e-16) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else if (l <= -6.4e-301) {
		tmp = t_0 * (((t_1 * (h * 0.5)) / l) + 1.0);
	} else if (l <= 4.2e+160) {
		tmp = t_0 * ((h * (-0.5 * (t_1 / l))) + 1.0);
	} else {
		tmp = d * (sqrt((1.0 / 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) :: t_1
    real(8) :: tmp
    t_0 = d / sqrt((l * h))
    t_1 = (0.5d0 * (d_m * (m_m / d))) ** 2.0d0
    if (l <= (-7.5d-16)) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else if (l <= (-6.4d-301)) then
        tmp = t_0 * (((t_1 * (h * 0.5d0)) / l) + 1.0d0)
    else if (l <= 4.2d+160) then
        tmp = t_0 * ((h * ((-0.5d0) * (t_1 / l))) + 1.0d0)
    else
        tmp = d * (sqrt((1.0d0 / 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 = d / Math.sqrt((l * h));
	double t_1 = Math.pow((0.5 * (D_m * (M_m / d))), 2.0);
	double tmp;
	if (l <= -7.5e-16) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else if (l <= -6.4e-301) {
		tmp = t_0 * (((t_1 * (h * 0.5)) / l) + 1.0);
	} else if (l <= 4.2e+160) {
		tmp = t_0 * ((h * (-0.5 * (t_1 / l))) + 1.0);
	} else {
		tmp = d * (Math.sqrt((1.0 / 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 = d / math.sqrt((l * h))
	t_1 = math.pow((0.5 * (D_m * (M_m / d))), 2.0)
	tmp = 0
	if l <= -7.5e-16:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	elif l <= -6.4e-301:
		tmp = t_0 * (((t_1 * (h * 0.5)) / l) + 1.0)
	elif l <= 4.2e+160:
		tmp = t_0 * ((h * (-0.5 * (t_1 / l))) + 1.0)
	else:
		tmp = d * (math.sqrt((1.0 / 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(d / sqrt(Float64(l * h)))
	t_1 = Float64(0.5 * Float64(D_m * Float64(M_m / d))) ^ 2.0
	tmp = 0.0
	if (l <= -7.5e-16)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	elseif (l <= -6.4e-301)
		tmp = Float64(t_0 * Float64(Float64(Float64(t_1 * Float64(h * 0.5)) / l) + 1.0));
	elseif (l <= 4.2e+160)
		tmp = Float64(t_0 * Float64(Float64(h * Float64(-0.5 * Float64(t_1 / l))) + 1.0));
	else
		tmp = Float64(d * Float64(sqrt(Float64(1.0 / 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 = d / sqrt((l * h));
	t_1 = (0.5 * (D_m * (M_m / d))) ^ 2.0;
	tmp = 0.0;
	if (l <= -7.5e-16)
		tmp = d * -sqrt((1.0 / (l * h)));
	elseif (l <= -6.4e-301)
		tmp = t_0 * (((t_1 * (h * 0.5)) / l) + 1.0);
	elseif (l <= 4.2e+160)
		tmp = t_0 * ((h * (-0.5 * (t_1 / l))) + 1.0);
	else
		tmp = d * (sqrt((1.0 / 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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[Power[N[(0.5 * N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]}, If[LessEqual[l, -7.5e-16], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[l, -6.4e-301], N[(t$95$0 * N[(N[(N[(t$95$1 * N[(h * 0.5), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 4.2e+160], N[(t$95$0 * N[(N[(h * N[(-0.5 * N[(t$95$1 / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Sqrt[N[(1.0 / l), $MachinePrecision]], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

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

Alternative 16: 51.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} t_0 := \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{if}\;d \leq -5.5 \cdot 10^{+90}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{elif}\;d \leq -4.5 \cdot 10^{+42}:\\ \;\;\;\;t_0 \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;d \leq -6.8 \cdot 10^{-13}:\\ \;\;\;\;\sqrt{\frac{{d}^{2}}{\ell \cdot h}}\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_0 \cdot \left(\frac{\frac{h \cdot -0.5}{\frac{-2}{D_m \cdot \frac{M_m}{d}}}}{\left(d \cdot 2\right) \cdot \frac{\ell}{D_m \cdot M_m}} + 1\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \frac{\sqrt{\frac{1}{\ell}}}{\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 (/ d (sqrt (* l h)))))
   (if (<= d -5.5e+90)
     (* d (- (sqrt (/ 1.0 (* l h)))))
     (if (<= d -4.5e+42)
       (*
        t_0
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= d -6.8e-13)
         (sqrt (/ (pow d 2.0) (* l h)))
         (if (<= d -4e-310)
           (*
            t_0
            (+
             (/
              (/ (* h -0.5) (/ (- 2.0) (* D_m (/ M_m d))))
              (* (* d 2.0) (/ l (* D_m M_m))))
             1.0))
           (* d (/ (sqrt (/ 1.0 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 = d / sqrt((l * h));
	double tmp;
	if (d <= -5.5e+90) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else if (d <= -4.5e+42) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -6.8e-13) {
		tmp = sqrt((pow(d, 2.0) / (l * h)));
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} else {
		tmp = d * (sqrt((1.0 / 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 = d / sqrt((l * h))
    if (d <= (-5.5d+90)) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else if (d <= (-4.5d+42)) then
        tmp = t_0 * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (d <= (-6.8d-13)) then
        tmp = sqrt(((d ** 2.0d0) / (l * h)))
    else if (d <= (-4d-310)) then
        tmp = t_0 * ((((h * (-0.5d0)) / (-2.0d0 / (d_m * (m_m / d)))) / ((d * 2.0d0) * (l / (d_m * m_m)))) + 1.0d0)
    else
        tmp = d * (sqrt((1.0d0 / 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 = d / Math.sqrt((l * h));
	double tmp;
	if (d <= -5.5e+90) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else if (d <= -4.5e+42) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -6.8e-13) {
		tmp = Math.sqrt((Math.pow(d, 2.0) / (l * h)));
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} else {
		tmp = d * (Math.sqrt((1.0 / 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 = d / math.sqrt((l * h))
	tmp = 0
	if d <= -5.5e+90:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	elif d <= -4.5e+42:
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif d <= -6.8e-13:
		tmp = math.sqrt((math.pow(d, 2.0) / (l * h)))
	elif d <= -4e-310:
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0)
	else:
		tmp = d * (math.sqrt((1.0 / 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(d / sqrt(Float64(l * h)))
	tmp = 0.0
	if (d <= -5.5e+90)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	elseif (d <= -4.5e+42)
		tmp = Float64(t_0 * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (d <= -6.8e-13)
		tmp = sqrt(Float64((d ^ 2.0) / Float64(l * h)));
	elseif (d <= -4e-310)
		tmp = Float64(t_0 * Float64(Float64(Float64(Float64(h * -0.5) / Float64(Float64(-2.0) / Float64(D_m * Float64(M_m / d)))) / Float64(Float64(d * 2.0) * Float64(l / Float64(D_m * M_m)))) + 1.0));
	else
		tmp = Float64(d * Float64(sqrt(Float64(1.0 / 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 = d / sqrt((l * h));
	tmp = 0.0;
	if (d <= -5.5e+90)
		tmp = d * -sqrt((1.0 / (l * h)));
	elseif (d <= -4.5e+42)
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (d <= -6.8e-13)
		tmp = sqrt(((d ^ 2.0) / (l * h)));
	elseif (d <= -4e-310)
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	else
		tmp = d * (sqrt((1.0 / 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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[d, -5.5e+90], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[d, -4.5e+42], N[(t$95$0 * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -6.8e-13], N[Sqrt[N[(N[Power[d, 2.0], $MachinePrecision] / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], If[LessEqual[d, -4e-310], N[(t$95$0 * N[(N[(N[(N[(h * -0.5), $MachinePrecision] / N[((-2.0) / N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(d * 2.0), $MachinePrecision] * N[(l / N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Sqrt[N[(1.0 / l), $MachinePrecision]], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

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

Alternative 17: 51.6% 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 := \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{if}\;d \leq -5.4 \cdot 10^{+81}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{elif}\;d \leq -2.45 \cdot 10^{+43}:\\ \;\;\;\;t_0 \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;d \leq -2.5 \cdot 10^{-93}:\\ \;\;\;\;\sqrt{\frac{d}{\ell}} \cdot \sqrt{\frac{d}{h}}\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_0 \cdot \left(\frac{\frac{h \cdot -0.5}{\frac{-2}{D_m \cdot \frac{M_m}{d}}}}{\left(d \cdot 2\right) \cdot \frac{\ell}{D_m \cdot M_m}} + 1\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \frac{\sqrt{\frac{1}{\ell}}}{\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 (/ d (sqrt (* l h)))))
   (if (<= d -5.4e+81)
     (* d (- (sqrt (/ 1.0 (* l h)))))
     (if (<= d -2.45e+43)
       (*
        t_0
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= d -2.5e-93)
         (* (sqrt (/ d l)) (sqrt (/ d h)))
         (if (<= d -4e-310)
           (*
            t_0
            (+
             (/
              (/ (* h -0.5) (/ (- 2.0) (* D_m (/ M_m d))))
              (* (* d 2.0) (/ l (* D_m M_m))))
             1.0))
           (* d (/ (sqrt (/ 1.0 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 = d / sqrt((l * h));
	double tmp;
	if (d <= -5.4e+81) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else if (d <= -2.45e+43) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -2.5e-93) {
		tmp = sqrt((d / l)) * sqrt((d / h));
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} else {
		tmp = d * (sqrt((1.0 / 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 = d / sqrt((l * h))
    if (d <= (-5.4d+81)) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else if (d <= (-2.45d+43)) then
        tmp = t_0 * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (d <= (-2.5d-93)) then
        tmp = sqrt((d / l)) * sqrt((d / h))
    else if (d <= (-4d-310)) then
        tmp = t_0 * ((((h * (-0.5d0)) / (-2.0d0 / (d_m * (m_m / d)))) / ((d * 2.0d0) * (l / (d_m * m_m)))) + 1.0d0)
    else
        tmp = d * (sqrt((1.0d0 / 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 = d / Math.sqrt((l * h));
	double tmp;
	if (d <= -5.4e+81) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else if (d <= -2.45e+43) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -2.5e-93) {
		tmp = Math.sqrt((d / l)) * Math.sqrt((d / h));
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} else {
		tmp = d * (Math.sqrt((1.0 / 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 = d / math.sqrt((l * h))
	tmp = 0
	if d <= -5.4e+81:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	elif d <= -2.45e+43:
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif d <= -2.5e-93:
		tmp = math.sqrt((d / l)) * math.sqrt((d / h))
	elif d <= -4e-310:
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0)
	else:
		tmp = d * (math.sqrt((1.0 / 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(d / sqrt(Float64(l * h)))
	tmp = 0.0
	if (d <= -5.4e+81)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	elseif (d <= -2.45e+43)
		tmp = Float64(t_0 * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (d <= -2.5e-93)
		tmp = Float64(sqrt(Float64(d / l)) * sqrt(Float64(d / h)));
	elseif (d <= -4e-310)
		tmp = Float64(t_0 * Float64(Float64(Float64(Float64(h * -0.5) / Float64(Float64(-2.0) / Float64(D_m * Float64(M_m / d)))) / Float64(Float64(d * 2.0) * Float64(l / Float64(D_m * M_m)))) + 1.0));
	else
		tmp = Float64(d * Float64(sqrt(Float64(1.0 / 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 = d / sqrt((l * h));
	tmp = 0.0;
	if (d <= -5.4e+81)
		tmp = d * -sqrt((1.0 / (l * h)));
	elseif (d <= -2.45e+43)
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (d <= -2.5e-93)
		tmp = sqrt((d / l)) * sqrt((d / h));
	elseif (d <= -4e-310)
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	else
		tmp = d * (sqrt((1.0 / 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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[d, -5.4e+81], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[d, -2.45e+43], N[(t$95$0 * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -2.5e-93], N[(N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -4e-310], N[(t$95$0 * N[(N[(N[(N[(h * -0.5), $MachinePrecision] / N[((-2.0) / N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(d * 2.0), $MachinePrecision] * N[(l / N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Sqrt[N[(1.0 / l), $MachinePrecision]], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

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

Alternative 18: 51.8% 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} t_0 := \frac{d}{\sqrt{\ell \cdot h}}\\ t_1 := d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{if}\;d \leq -1.3 \cdot 10^{+93}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq -1.55 \cdot 10^{+44}:\\ \;\;\;\;t_0 \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;d \leq -6.5 \cdot 10^{-9}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_0 \cdot \left(\frac{\frac{h \cdot -0.5}{\frac{-2}{D_m \cdot \frac{M_m}{d}}}}{\left(d \cdot 2\right) \cdot \frac{\ell}{D_m \cdot M_m}} + 1\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
 (let* ((t_0 (/ d (sqrt (* l h)))) (t_1 (* d (- (sqrt (/ 1.0 (* l h)))))))
   (if (<= d -1.3e+93)
     t_1
     (if (<= d -1.55e+44)
       (*
        t_0
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= d -6.5e-9)
         t_1
         (if (<= d -4e-310)
           (*
            t_0
            (+
             (/
              (/ (* h -0.5) (/ (- 2.0) (* D_m (/ M_m d))))
              (* (* d 2.0) (/ l (* D_m M_m))))
             1.0))
           (/ 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 = d / sqrt((l * h));
	double t_1 = d * -sqrt((1.0 / (l * h)));
	double tmp;
	if (d <= -1.3e+93) {
		tmp = t_1;
	} else if (d <= -1.55e+44) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -6.5e-9) {
		tmp = t_1;
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} 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) :: t_1
    real(8) :: tmp
    t_0 = d / sqrt((l * h))
    t_1 = d * -sqrt((1.0d0 / (l * h)))
    if (d <= (-1.3d+93)) then
        tmp = t_1
    else if (d <= (-1.55d+44)) then
        tmp = t_0 * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (d <= (-6.5d-9)) then
        tmp = t_1
    else if (d <= (-4d-310)) then
        tmp = t_0 * ((((h * (-0.5d0)) / (-2.0d0 / (d_m * (m_m / d)))) / ((d * 2.0d0) * (l / (d_m * m_m)))) + 1.0d0)
    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 = d / Math.sqrt((l * h));
	double t_1 = d * -Math.sqrt((1.0 / (l * h)));
	double tmp;
	if (d <= -1.3e+93) {
		tmp = t_1;
	} else if (d <= -1.55e+44) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -6.5e-9) {
		tmp = t_1;
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} 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 = d / math.sqrt((l * h))
	t_1 = d * -math.sqrt((1.0 / (l * h)))
	tmp = 0
	if d <= -1.3e+93:
		tmp = t_1
	elif d <= -1.55e+44:
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif d <= -6.5e-9:
		tmp = t_1
	elif d <= -4e-310:
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0)
	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(d / sqrt(Float64(l * h)))
	t_1 = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))))
	tmp = 0.0
	if (d <= -1.3e+93)
		tmp = t_1;
	elseif (d <= -1.55e+44)
		tmp = Float64(t_0 * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (d <= -6.5e-9)
		tmp = t_1;
	elseif (d <= -4e-310)
		tmp = Float64(t_0 * Float64(Float64(Float64(Float64(h * -0.5) / Float64(Float64(-2.0) / Float64(D_m * Float64(M_m / d)))) / Float64(Float64(d * 2.0) * Float64(l / Float64(D_m * M_m)))) + 1.0));
	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 = d / sqrt((l * h));
	t_1 = d * -sqrt((1.0 / (l * h)));
	tmp = 0.0;
	if (d <= -1.3e+93)
		tmp = t_1;
	elseif (d <= -1.55e+44)
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (d <= -6.5e-9)
		tmp = t_1;
	elseif (d <= -4e-310)
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision]}, If[LessEqual[d, -1.3e+93], t$95$1, If[LessEqual[d, -1.55e+44], N[(t$95$0 * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -6.5e-9], t$95$1, If[LessEqual[d, -4e-310], N[(t$95$0 * N[(N[(N[(N[(h * -0.5), $MachinePrecision] / N[((-2.0) / N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(d * 2.0), $MachinePrecision] * N[(l / N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]]

Derivation

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

Alternative 19: 52.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} t_0 := \frac{d}{\sqrt{\ell \cdot h}}\\ \mathbf{if}\;d \leq -4.5 \cdot 10^{+83}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{elif}\;d \leq -1.95 \cdot 10^{+42}:\\ \;\;\;\;t_0 \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;d \leq -0.0021:\\ \;\;\;\;\sqrt{\frac{{d}^{2}}{\ell \cdot h}}\\ \mathbf{elif}\;d \leq -4 \cdot 10^{-310}:\\ \;\;\;\;t_0 \cdot \left(\frac{\frac{h \cdot -0.5}{\frac{-2}{D_m \cdot \frac{M_m}{d}}}}{\left(d \cdot 2\right) \cdot \frac{\ell}{D_m \cdot M_m}} + 1\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
 (let* ((t_0 (/ d (sqrt (* l h)))))
   (if (<= d -4.5e+83)
     (* d (- (sqrt (/ 1.0 (* l h)))))
     (if (<= d -1.95e+42)
       (*
        t_0
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= d -0.0021)
         (sqrt (/ (pow d 2.0) (* l h)))
         (if (<= d -4e-310)
           (*
            t_0
            (+
             (/
              (/ (* h -0.5) (/ (- 2.0) (* D_m (/ M_m d))))
              (* (* d 2.0) (/ l (* D_m M_m))))
             1.0))
           (/ 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 = d / sqrt((l * h));
	double tmp;
	if (d <= -4.5e+83) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else if (d <= -1.95e+42) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -0.0021) {
		tmp = sqrt((pow(d, 2.0) / (l * h)));
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} 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 = d / sqrt((l * h))
    if (d <= (-4.5d+83)) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else if (d <= (-1.95d+42)) then
        tmp = t_0 * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (d <= (-0.0021d0)) then
        tmp = sqrt(((d ** 2.0d0) / (l * h)))
    else if (d <= (-4d-310)) then
        tmp = t_0 * ((((h * (-0.5d0)) / (-2.0d0 / (d_m * (m_m / d)))) / ((d * 2.0d0) * (l / (d_m * m_m)))) + 1.0d0)
    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 = d / Math.sqrt((l * h));
	double tmp;
	if (d <= -4.5e+83) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else if (d <= -1.95e+42) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -0.0021) {
		tmp = Math.sqrt((Math.pow(d, 2.0) / (l * h)));
	} else if (d <= -4e-310) {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	} 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 = d / math.sqrt((l * h))
	tmp = 0
	if d <= -4.5e+83:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	elif d <= -1.95e+42:
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif d <= -0.0021:
		tmp = math.sqrt((math.pow(d, 2.0) / (l * h)))
	elif d <= -4e-310:
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0)
	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(d / sqrt(Float64(l * h)))
	tmp = 0.0
	if (d <= -4.5e+83)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	elseif (d <= -1.95e+42)
		tmp = Float64(t_0 * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (d <= -0.0021)
		tmp = sqrt(Float64((d ^ 2.0) / Float64(l * h)));
	elseif (d <= -4e-310)
		tmp = Float64(t_0 * Float64(Float64(Float64(Float64(h * -0.5) / Float64(Float64(-2.0) / Float64(D_m * Float64(M_m / d)))) / Float64(Float64(d * 2.0) * Float64(l / Float64(D_m * M_m)))) + 1.0));
	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 = d / sqrt((l * h));
	tmp = 0.0;
	if (d <= -4.5e+83)
		tmp = d * -sqrt((1.0 / (l * h)));
	elseif (d <= -1.95e+42)
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (d <= -0.0021)
		tmp = sqrt(((d ^ 2.0) / (l * h)));
	elseif (d <= -4e-310)
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.0);
	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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, If[LessEqual[d, -4.5e+83], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], If[LessEqual[d, -1.95e+42], N[(t$95$0 * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -0.0021], N[Sqrt[N[(N[Power[d, 2.0], $MachinePrecision] / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], If[LessEqual[d, -4e-310], N[(t$95$0 * N[(N[(N[(N[(h * -0.5), $MachinePrecision] / N[((-2.0) / N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(d * 2.0), $MachinePrecision] * N[(l / N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

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

Alternative 20: 47.3% accurate, 2.4× 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 := \frac{d}{\sqrt{\ell \cdot h}}\\ t_1 := d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{if}\;d \leq -3.6 \cdot 10^{+95}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;d \leq -2.6 \cdot 10^{+42}:\\ \;\;\;\;t_0 \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;d \leq -7 \cdot 10^{-6}:\\ \;\;\;\;t_1\\ \mathbf{else}:\\ \;\;\;\;t_0 \cdot \left(\frac{\frac{h \cdot -0.5}{\frac{-2}{D_m \cdot \frac{M_m}{d}}}}{\left(d \cdot 2\right) \cdot \frac{\ell}{D_m \cdot M_m}} + 1\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 (/ d (sqrt (* l h)))) (t_1 (* d (- (sqrt (/ 1.0 (* l h)))))))
   (if (<= d -3.6e+95)
     t_1
     (if (<= d -2.6e+42)
       (*
        t_0
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= d -7e-6)
         t_1
         (*
          t_0
          (+
           (/
            (/ (* h -0.5) (/ (- 2.0) (* D_m (/ M_m d))))
            (* (* d 2.0) (/ l (* D_m M_m))))
           1.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 = d / sqrt((l * h));
	double t_1 = d * -sqrt((1.0 / (l * h)));
	double tmp;
	if (d <= -3.6e+95) {
		tmp = t_1;
	} else if (d <= -2.6e+42) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -7e-6) {
		tmp = t_1;
	} else {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.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) :: t_1
    real(8) :: tmp
    t_0 = d / sqrt((l * h))
    t_1 = d * -sqrt((1.0d0 / (l * h)))
    if (d <= (-3.6d+95)) then
        tmp = t_1
    else if (d <= (-2.6d+42)) then
        tmp = t_0 * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (d <= (-7d-6)) then
        tmp = t_1
    else
        tmp = t_0 * ((((h * (-0.5d0)) / (-2.0d0 / (d_m * (m_m / d)))) / ((d * 2.0d0) * (l / (d_m * m_m)))) + 1.0d0)
    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 = d / Math.sqrt((l * h));
	double t_1 = d * -Math.sqrt((1.0 / (l * h)));
	double tmp;
	if (d <= -3.6e+95) {
		tmp = t_1;
	} else if (d <= -2.6e+42) {
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (d <= -7e-6) {
		tmp = t_1;
	} else {
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.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 = d / math.sqrt((l * h))
	t_1 = d * -math.sqrt((1.0 / (l * h)))
	tmp = 0
	if d <= -3.6e+95:
		tmp = t_1
	elif d <= -2.6e+42:
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif d <= -7e-6:
		tmp = t_1
	else:
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.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(d / sqrt(Float64(l * h)))
	t_1 = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))))
	tmp = 0.0
	if (d <= -3.6e+95)
		tmp = t_1;
	elseif (d <= -2.6e+42)
		tmp = Float64(t_0 * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (d <= -7e-6)
		tmp = t_1;
	else
		tmp = Float64(t_0 * Float64(Float64(Float64(Float64(h * -0.5) / Float64(Float64(-2.0) / Float64(D_m * Float64(M_m / d)))) / Float64(Float64(d * 2.0) * Float64(l / Float64(D_m * M_m)))) + 1.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 = d / sqrt((l * h));
	t_1 = d * -sqrt((1.0 / (l * h)));
	tmp = 0.0;
	if (d <= -3.6e+95)
		tmp = t_1;
	elseif (d <= -2.6e+42)
		tmp = t_0 * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (d <= -7e-6)
		tmp = t_1;
	else
		tmp = t_0 * ((((h * -0.5) / (-2.0 / (D_m * (M_m / d)))) / ((d * 2.0) * (l / (D_m * M_m)))) + 1.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[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]}, Block[{t$95$1 = N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision]}, If[LessEqual[d, -3.6e+95], t$95$1, If[LessEqual[d, -2.6e+42], N[(t$95$0 * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[d, -7e-6], t$95$1, N[(t$95$0 * N[(N[(N[(N[(h * -0.5), $MachinePrecision] / N[((-2.0) / N[(D$95$m * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / N[(N[(d * 2.0), $MachinePrecision] * N[(l / N[(D$95$m * M$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

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

Alternative 21: 45.1% accurate, 2.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 := d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{if}\;\ell \leq -2.35 \cdot 10^{-15}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;\ell \leq -4 \cdot 10^{-287}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(\left(\left(D_m \cdot M_m\right) \cdot \left(\frac{h}{\ell} \cdot \left(\frac{D_m}{2} \cdot \frac{M_m}{d}\right)\right)\right) \cdot \frac{0.25}{d} + 1\right)\\ \mathbf{elif}\;\ell \leq 8.2 \cdot 10^{-281}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;d \cdot \sqrt{\frac{\frac{1}{h}}{\ell}}\\ \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 (* d (- (sqrt (/ 1.0 (* l h)))))))
   (if (<= l -2.35e-15)
     t_0
     (if (<= l -4e-287)
       (*
        (/ d (sqrt (* l h)))
        (+
         (* (* (* D_m M_m) (* (/ h l) (* (/ D_m 2.0) (/ M_m d)))) (/ 0.25 d))
         1.0))
       (if (<= l 8.2e-281) t_0 (* d (sqrt (/ (/ 1.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 t_0 = d * -sqrt((1.0 / (l * h)));
	double tmp;
	if (l <= -2.35e-15) {
		tmp = t_0;
	} else if (l <= -4e-287) {
		tmp = (d / sqrt((l * h))) * ((((D_m * M_m) * ((h / l) * ((D_m / 2.0) * (M_m / d)))) * (0.25 / d)) + 1.0);
	} else if (l <= 8.2e-281) {
		tmp = t_0;
	} else {
		tmp = d * sqrt(((1.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) :: t_0
    real(8) :: tmp
    t_0 = d * -sqrt((1.0d0 / (l * h)))
    if (l <= (-2.35d-15)) then
        tmp = t_0
    else if (l <= (-4d-287)) then
        tmp = (d / sqrt((l * h))) * ((((d_m * m_m) * ((h / l) * ((d_m / 2.0d0) * (m_m / d)))) * (0.25d0 / d)) + 1.0d0)
    else if (l <= 8.2d-281) then
        tmp = t_0
    else
        tmp = d * sqrt(((1.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 t_0 = d * -Math.sqrt((1.0 / (l * h)));
	double tmp;
	if (l <= -2.35e-15) {
		tmp = t_0;
	} else if (l <= -4e-287) {
		tmp = (d / Math.sqrt((l * h))) * ((((D_m * M_m) * ((h / l) * ((D_m / 2.0) * (M_m / d)))) * (0.25 / d)) + 1.0);
	} else if (l <= 8.2e-281) {
		tmp = t_0;
	} else {
		tmp = d * Math.sqrt(((1.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):
	t_0 = d * -math.sqrt((1.0 / (l * h)))
	tmp = 0
	if l <= -2.35e-15:
		tmp = t_0
	elif l <= -4e-287:
		tmp = (d / math.sqrt((l * h))) * ((((D_m * M_m) * ((h / l) * ((D_m / 2.0) * (M_m / d)))) * (0.25 / d)) + 1.0)
	elif l <= 8.2e-281:
		tmp = t_0
	else:
		tmp = d * math.sqrt(((1.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)
	t_0 = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))))
	tmp = 0.0
	if (l <= -2.35e-15)
		tmp = t_0;
	elseif (l <= -4e-287)
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(Float64(Float64(Float64(D_m * M_m) * Float64(Float64(h / l) * Float64(Float64(D_m / 2.0) * Float64(M_m / d)))) * Float64(0.25 / d)) + 1.0));
	elseif (l <= 8.2e-281)
		tmp = t_0;
	else
		tmp = Float64(d * sqrt(Float64(Float64(1.0 / 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)
	t_0 = d * -sqrt((1.0 / (l * h)));
	tmp = 0.0;
	if (l <= -2.35e-15)
		tmp = t_0;
	elseif (l <= -4e-287)
		tmp = (d / sqrt((l * h))) * ((((D_m * M_m) * ((h / l) * ((D_m / 2.0) * (M_m / d)))) * (0.25 / d)) + 1.0);
	elseif (l <= 8.2e-281)
		tmp = t_0;
	else
		tmp = d * sqrt(((1.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_] := Block[{t$95$0 = N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision]}, If[LessEqual[l, -2.35e-15], t$95$0, If[LessEqual[l, -4e-287], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(N[(N[(D$95$m * M$95$m), $MachinePrecision] * N[(N[(h / l), $MachinePrecision] * N[(N[(D$95$m / 2.0), $MachinePrecision] * N[(M$95$m / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(0.25 / d), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 8.2e-281], t$95$0, N[(d * N[Sqrt[N[(N[(1.0 / h), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]]]]

Derivation

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

Alternative 22: 44.8% accurate, 2.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 := d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{if}\;\ell \leq -4.2 \cdot 10^{-15}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;\ell \leq -1.4 \cdot 10^{-288}:\\ \;\;\;\;\frac{d}{\sqrt{\ell \cdot h}} \cdot \left(M_m \cdot \frac{\frac{h}{\frac{\frac{\ell}{\frac{D_m}{d}}}{M_m}} \cdot -0.25}{-2 \cdot \frac{d}{D_m}} + 1\right)\\ \mathbf{elif}\;\ell \leq 6.2 \cdot 10^{-282}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;d \cdot \sqrt{\frac{\frac{1}{h}}{\ell}}\\ \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 (* d (- (sqrt (/ 1.0 (* l h)))))))
   (if (<= l -4.2e-15)
     t_0
     (if (<= l -1.4e-288)
       (*
        (/ d (sqrt (* l h)))
        (+
         (* M_m (/ (* (/ h (/ (/ l (/ D_m d)) M_m)) -0.25) (* -2.0 (/ d D_m))))
         1.0))
       (if (<= l 6.2e-282) t_0 (* d (sqrt (/ (/ 1.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 t_0 = d * -sqrt((1.0 / (l * h)));
	double tmp;
	if (l <= -4.2e-15) {
		tmp = t_0;
	} else if (l <= -1.4e-288) {
		tmp = (d / sqrt((l * h))) * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (l <= 6.2e-282) {
		tmp = t_0;
	} else {
		tmp = d * sqrt(((1.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) :: t_0
    real(8) :: tmp
    t_0 = d * -sqrt((1.0d0 / (l * h)))
    if (l <= (-4.2d-15)) then
        tmp = t_0
    else if (l <= (-1.4d-288)) then
        tmp = (d / sqrt((l * h))) * ((m_m * (((h / ((l / (d_m / d)) / m_m)) * (-0.25d0)) / ((-2.0d0) * (d / d_m)))) + 1.0d0)
    else if (l <= 6.2d-282) then
        tmp = t_0
    else
        tmp = d * sqrt(((1.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 t_0 = d * -Math.sqrt((1.0 / (l * h)));
	double tmp;
	if (l <= -4.2e-15) {
		tmp = t_0;
	} else if (l <= -1.4e-288) {
		tmp = (d / Math.sqrt((l * h))) * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	} else if (l <= 6.2e-282) {
		tmp = t_0;
	} else {
		tmp = d * Math.sqrt(((1.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):
	t_0 = d * -math.sqrt((1.0 / (l * h)))
	tmp = 0
	if l <= -4.2e-15:
		tmp = t_0
	elif l <= -1.4e-288:
		tmp = (d / math.sqrt((l * h))) * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0)
	elif l <= 6.2e-282:
		tmp = t_0
	else:
		tmp = d * math.sqrt(((1.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)
	t_0 = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))))
	tmp = 0.0
	if (l <= -4.2e-15)
		tmp = t_0;
	elseif (l <= -1.4e-288)
		tmp = Float64(Float64(d / sqrt(Float64(l * h))) * Float64(Float64(M_m * Float64(Float64(Float64(h / Float64(Float64(l / Float64(D_m / d)) / M_m)) * -0.25) / Float64(-2.0 * Float64(d / D_m)))) + 1.0));
	elseif (l <= 6.2e-282)
		tmp = t_0;
	else
		tmp = Float64(d * sqrt(Float64(Float64(1.0 / 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)
	t_0 = d * -sqrt((1.0 / (l * h)));
	tmp = 0.0;
	if (l <= -4.2e-15)
		tmp = t_0;
	elseif (l <= -1.4e-288)
		tmp = (d / sqrt((l * h))) * ((M_m * (((h / ((l / (D_m / d)) / M_m)) * -0.25) / (-2.0 * (d / D_m)))) + 1.0);
	elseif (l <= 6.2e-282)
		tmp = t_0;
	else
		tmp = d * sqrt(((1.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_] := Block[{t$95$0 = N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision]}, If[LessEqual[l, -4.2e-15], t$95$0, If[LessEqual[l, -1.4e-288], N[(N[(d / N[Sqrt[N[(l * h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(N[(M$95$m * N[(N[(N[(h / N[(N[(l / N[(D$95$m / d), $MachinePrecision]), $MachinePrecision] / M$95$m), $MachinePrecision]), $MachinePrecision] * -0.25), $MachinePrecision] / N[(-2.0 * N[(d / D$95$m), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 6.2e-282], t$95$0, N[(d * N[Sqrt[N[(N[(1.0 / h), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]]]]

Derivation

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

Alternative 23: 41.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 \cdot 10^{-142}:\\ \;\;\;\;d \cdot \left(-\sqrt{\frac{1}{\ell \cdot h}}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \sqrt{\frac{\frac{1}{h}}{\ell}}\\ \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 -1e-142)
   (* d (- (sqrt (/ 1.0 (* l h)))))
   (* d (sqrt (/ (/ 1.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 (d <= -1e-142) {
		tmp = d * -sqrt((1.0 / (l * h)));
	} else {
		tmp = d * sqrt(((1.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 (d <= (-1d-142)) then
        tmp = d * -sqrt((1.0d0 / (l * h)))
    else
        tmp = d * sqrt(((1.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 (d <= -1e-142) {
		tmp = d * -Math.sqrt((1.0 / (l * h)));
	} else {
		tmp = d * Math.sqrt(((1.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 d <= -1e-142:
		tmp = d * -math.sqrt((1.0 / (l * h)))
	else:
		tmp = d * math.sqrt(((1.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 (d <= -1e-142)
		tmp = Float64(d * Float64(-sqrt(Float64(1.0 / Float64(l * h)))));
	else
		tmp = Float64(d * sqrt(Float64(Float64(1.0 / 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 (d <= -1e-142)
		tmp = d * -sqrt((1.0 / (l * h)));
	else
		tmp = d * sqrt(((1.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[d, -1e-142], N[(d * (-N[Sqrt[N[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision])), $MachinePrecision], N[(d * N[Sqrt[N[(N[(1.0 / h), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]

Derivation

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

Alternative 24: 25.9% accurate, 3.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])\\ \\ d \cdot \sqrt{\frac{1}{\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 (/ 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) {
	return d * sqrt((1.0 / (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((1.0d0 / (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((1.0 / (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((1.0 / (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(1.0 / 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((1.0 / (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[(1.0 / N[(l * h), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

Derivation

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

Alternative 25: 26.1% accurate, 3.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])\\ \\ d \cdot \sqrt{\frac{\frac{1}{h}}{\ell}} \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 (/ (/ 1.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) {
	return d * sqrt(((1.0 / h) / l));
}

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(((1.0d0 / h) / l))
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(((1.0 / h) / l));
}

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(((1.0 / h) / l))

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(Float64(1.0 / h) / l)))
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(((1.0 / h) / l));
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[(N[(1.0 / h), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

Derivation

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

Alternative 26: 25.7% accurate, 3.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])\\ \\ d \cdot {\left(\ell \cdot h\right)}^{-0.5} \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 (pow (* l h) -0.5)))

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 * pow((l * h), -0.5);
}

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 * ((l * h) ** (-0.5d0))
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.pow((l * h), -0.5);
}

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.pow((l * h), -0.5)

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 * (Float64(l * h) ^ -0.5))
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 * ((l * h) ^ -0.5);
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[Power[N[(l * h), $MachinePrecision], -0.5], $MachinePrecision]), $MachinePrecision]

Derivation

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

Alternative 27: 25.7% 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 2024010 
(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)))))