Henrywood and Agarwal, Equation (12)

Percentage Accurate: 66.6% → 71.0%

Time: 22.5s

Alternatives: 17

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.6% 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: 71.0% accurate, 0.8× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} t_0 := \sqrt{\frac{d}{\ell}}\\ t_1 := \left(t_0 \cdot \frac{\sqrt{d}}{\sqrt{h}}\right) \cdot \left(1 - 0.5 \cdot \left(\frac{h}{\ell} \cdot {\left(\frac{D}{d} \cdot \frac{M}{2}\right)}^{2}\right)\right)\\ t_2 := \sqrt{\frac{d}{h}} \cdot t_0\\ \mathbf{if}\;\ell \leq 5 \cdot 10^{-308}:\\ \;\;\;\;t_2 \cdot \left(1 - 0.5 \cdot \frac{h \cdot {\left(\frac{D}{d} \cdot \left(0.5 \cdot M\right)\right)}^{2}}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 1.1 \cdot 10^{-201}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;\ell \leq 5.5 \cdot 10^{-29}:\\ \;\;\;\;t_2 \cdot \left(1 - 0.5 \cdot \frac{0.25 \cdot \left(\frac{\left(D \cdot M\right) \cdot \left(D \cdot M\right)}{d} \cdot \frac{h}{d}\right)}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 3 \cdot 10^{+122}:\\ \;\;\;\;t_1\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (let* ((t_0 (sqrt (/ d l)))
        (t_1
         (*
          (* t_0 (/ (sqrt d) (sqrt h)))
          (- 1.0 (* 0.5 (* (/ h l) (pow (* (/ D d) (/ M 2.0)) 2.0))))))
        (t_2 (* (sqrt (/ d h)) t_0)))
   (if (<= l 5e-308)
     (* t_2 (- 1.0 (* 0.5 (/ (* h (pow (* (/ D d) (* 0.5 M)) 2.0)) l))))
     (if (<= l 1.1e-201)
       t_1
       (if (<= l 5.5e-29)
         (*
          t_2
          (- 1.0 (* 0.5 (/ (* 0.25 (* (/ (* (* D M) (* D M)) d) (/ h d))) l))))
         (if (<= l 3e+122) t_1 (* d (* (pow h -0.5) (pow l -0.5)))))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double t_0 = sqrt((d / l));
	double t_1 = (t_0 * (sqrt(d) / sqrt(h))) * (1.0 - (0.5 * ((h / l) * pow(((D / d) * (M / 2.0)), 2.0))));
	double t_2 = sqrt((d / h)) * t_0;
	double tmp;
	if (l <= 5e-308) {
		tmp = t_2 * (1.0 - (0.5 * ((h * pow(((D / d) * (0.5 * M)), 2.0)) / l)));
	} else if (l <= 1.1e-201) {
		tmp = t_1;
	} else if (l <= 5.5e-29) {
		tmp = t_2 * (1.0 - (0.5 * ((0.25 * ((((D * M) * (D * M)) / d) * (h / d))) / l)));
	} else if (l <= 3e+122) {
		tmp = t_1;
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: t_0
    real(8) :: t_1
    real(8) :: t_2
    real(8) :: tmp
    t_0 = sqrt((d / l))
    t_1 = (t_0 * (sqrt(d) / sqrt(h))) * (1.0d0 - (0.5d0 * ((h / l) * (((d_1 / d) * (m / 2.0d0)) ** 2.0d0))))
    t_2 = sqrt((d / h)) * t_0
    if (l <= 5d-308) then
        tmp = t_2 * (1.0d0 - (0.5d0 * ((h * (((d_1 / d) * (0.5d0 * m)) ** 2.0d0)) / l)))
    else if (l <= 1.1d-201) then
        tmp = t_1
    else if (l <= 5.5d-29) then
        tmp = t_2 * (1.0d0 - (0.5d0 * ((0.25d0 * ((((d_1 * m) * (d_1 * m)) / d) * (h / d))) / l)))
    else if (l <= 3d+122) then
        tmp = t_1
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double t_0 = Math.sqrt((d / l));
	double t_1 = (t_0 * (Math.sqrt(d) / Math.sqrt(h))) * (1.0 - (0.5 * ((h / l) * Math.pow(((D / d) * (M / 2.0)), 2.0))));
	double t_2 = Math.sqrt((d / h)) * t_0;
	double tmp;
	if (l <= 5e-308) {
		tmp = t_2 * (1.0 - (0.5 * ((h * Math.pow(((D / d) * (0.5 * M)), 2.0)) / l)));
	} else if (l <= 1.1e-201) {
		tmp = t_1;
	} else if (l <= 5.5e-29) {
		tmp = t_2 * (1.0 - (0.5 * ((0.25 * ((((D * M) * (D * M)) / d) * (h / d))) / l)));
	} else if (l <= 3e+122) {
		tmp = t_1;
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	t_0 = math.sqrt((d / l))
	t_1 = (t_0 * (math.sqrt(d) / math.sqrt(h))) * (1.0 - (0.5 * ((h / l) * math.pow(((D / d) * (M / 2.0)), 2.0))))
	t_2 = math.sqrt((d / h)) * t_0
	tmp = 0
	if l <= 5e-308:
		tmp = t_2 * (1.0 - (0.5 * ((h * math.pow(((D / d) * (0.5 * M)), 2.0)) / l)))
	elif l <= 1.1e-201:
		tmp = t_1
	elif l <= 5.5e-29:
		tmp = t_2 * (1.0 - (0.5 * ((0.25 * ((((D * M) * (D * M)) / d) * (h / d))) / l)))
	elif l <= 3e+122:
		tmp = t_1
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	t_0 = sqrt(Float64(d / l))
	t_1 = Float64(Float64(t_0 * Float64(sqrt(d) / sqrt(h))) * Float64(1.0 - Float64(0.5 * Float64(Float64(h / l) * (Float64(Float64(D / d) * Float64(M / 2.0)) ^ 2.0)))))
	t_2 = Float64(sqrt(Float64(d / h)) * t_0)
	tmp = 0.0
	if (l <= 5e-308)
		tmp = Float64(t_2 * Float64(1.0 - Float64(0.5 * Float64(Float64(h * (Float64(Float64(D / d) * Float64(0.5 * M)) ^ 2.0)) / l))));
	elseif (l <= 1.1e-201)
		tmp = t_1;
	elseif (l <= 5.5e-29)
		tmp = Float64(t_2 * Float64(1.0 - Float64(0.5 * Float64(Float64(0.25 * Float64(Float64(Float64(Float64(D * M) * Float64(D * M)) / d) * Float64(h / d))) / l))));
	elseif (l <= 3e+122)
		tmp = t_1;
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

MATLAB

M = abs(M)
D = abs(D)
M, D = num2cell(sort([M, D])){:}
function tmp_2 = code(d, h, l, M, D)
	t_0 = sqrt((d / l));
	t_1 = (t_0 * (sqrt(d) / sqrt(h))) * (1.0 - (0.5 * ((h / l) * (((D / d) * (M / 2.0)) ^ 2.0))));
	t_2 = sqrt((d / h)) * t_0;
	tmp = 0.0;
	if (l <= 5e-308)
		tmp = t_2 * (1.0 - (0.5 * ((h * (((D / d) * (0.5 * M)) ^ 2.0)) / l)));
	elseif (l <= 1.1e-201)
		tmp = t_1;
	elseif (l <= 5.5e-29)
		tmp = t_2 * (1.0 - (0.5 * ((0.25 * ((((D * M) * (D * M)) / d) * (h / d))) / l)));
	elseif (l <= 3e+122)
		tmp = t_1;
	else
		tmp = d * ((h ^ -0.5) * (l ^ -0.5));
	end
	tmp_2 = tmp;
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := Block[{t$95$0 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[(N[(t$95$0 * N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(h / l), $MachinePrecision] * N[Power[N[(N[(D / d), $MachinePrecision] * N[(M / 2.0), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, Block[{t$95$2 = N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * t$95$0), $MachinePrecision]}, If[LessEqual[l, 5e-308], N[(t$95$2 * N[(1.0 - N[(0.5 * N[(N[(h * N[Power[N[(N[(D / d), $MachinePrecision] * N[(0.5 * M), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 1.1e-201], t$95$1, If[LessEqual[l, 5.5e-29], N[(t$95$2 * N[(1.0 - N[(0.5 * N[(N[(0.25 * N[(N[(N[(N[(D * M), $MachinePrecision] * N[(D * M), $MachinePrecision]), $MachinePrecision] / d), $MachinePrecision] * N[(h / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 3e+122], t$95$1, N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]]]

Derivation

1. Initial program 74.7%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 5 \cdot 10^{-308}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \frac{h \cdot {\left(\frac{D}{d} \cdot \left(0.5 \cdot M\right)\right)}^{2}}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 1.1 \cdot 10^{-201}:\\ \;\;\;\;\left(\sqrt{\frac{d}{\ell}} \cdot \frac{\sqrt{d}}{\sqrt{h}}\right) \cdot \left(1 - 0.5 \cdot \left(\frac{h}{\ell} \cdot {\left(\frac{D}{d} \cdot \frac{M}{2}\right)}^{2}\right)\right)\\ \mathbf{elif}\;\ell \leq 5.5 \cdot 10^{-29}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \frac{0.25 \cdot \left(\frac{\left(D \cdot M\right) \cdot \left(D \cdot M\right)}{d} \cdot \frac{h}{d}\right)}{\ell}\right)\\ \mathbf{elif}\;\ell \leq 3 \cdot 10^{+122}:\\ \;\;\;\;\left(\sqrt{\frac{d}{\ell}} \cdot \frac{\sqrt{d}}{\sqrt{h}}\right) \cdot \left(1 - 0.5 \cdot \left(\frac{h}{\ell} \cdot {\left(\frac{D}{d} \cdot \frac{M}{2}\right)}^{2}\right)\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 2: 67.7% accurate, 0.6× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} t_0 := \sqrt{\frac{d}{h}}\\ t_1 := \sqrt{\frac{d}{\ell}}\\ \mathbf{if}\;h \leq 1.85 \cdot 10^{-278}:\\ \;\;\;\;\left(t_0 \cdot t_1\right) \cdot \left(1 - 0.5 \cdot \frac{h \cdot {\left(\frac{D}{d} \cdot \left(0.5 \cdot M\right)\right)}^{2}}{\ell}\right)\\ \mathbf{elif}\;h \leq 3 \cdot 10^{-252}:\\ \;\;\;\;\frac{d}{\sqrt{h \cdot \ell}}\\ \mathbf{elif}\;h \leq 1.6 \cdot 10^{+62}:\\ \;\;\;\;\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \left(t_0 \cdot \mathsf{fma}\left({\left(0.5 \cdot \frac{M}{\frac{d}{D}}\right)}^{2}, -0.5 \cdot \frac{h}{\ell}, 1\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\left(t_1 \cdot \frac{\sqrt{d}}{\sqrt{h}}\right) \cdot \left(1 - 0.5 \cdot \frac{0.25 \cdot \left(D \cdot D\right)}{\frac{d \cdot d}{\frac{M \cdot \left(h \cdot M\right)}{\ell}}}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (let* ((t_0 (sqrt (/ d h))) (t_1 (sqrt (/ d l))))
   (if (<= h 1.85e-278)
     (*
      (* t_0 t_1)
      (- 1.0 (* 0.5 (/ (* h (pow (* (/ D d) (* 0.5 M)) 2.0)) l))))
     (if (<= h 3e-252)
       (/ d (sqrt (* h l)))
       (if (<= h 1.6e+62)
         (*
          (/ (sqrt d) (sqrt l))
          (* t_0 (fma (pow (* 0.5 (/ M (/ d D))) 2.0) (* -0.5 (/ h l)) 1.0)))
         (*
          (* t_1 (/ (sqrt d) (sqrt h)))
          (-
           1.0
           (* 0.5 (/ (* 0.25 (* D D)) (/ (* d d) (/ (* M (* h M)) l)))))))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double t_0 = sqrt((d / h));
	double t_1 = sqrt((d / l));
	double tmp;
	if (h <= 1.85e-278) {
		tmp = (t_0 * t_1) * (1.0 - (0.5 * ((h * pow(((D / d) * (0.5 * M)), 2.0)) / l)));
	} else if (h <= 3e-252) {
		tmp = d / sqrt((h * l));
	} else if (h <= 1.6e+62) {
		tmp = (sqrt(d) / sqrt(l)) * (t_0 * fma(pow((0.5 * (M / (d / D))), 2.0), (-0.5 * (h / l)), 1.0));
	} else {
		tmp = (t_1 * (sqrt(d) / sqrt(h))) * (1.0 - (0.5 * ((0.25 * (D * D)) / ((d * d) / ((M * (h * M)) / l)))));
	}
	return tmp;
}

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	t_0 = sqrt(Float64(d / h))
	t_1 = sqrt(Float64(d / l))
	tmp = 0.0
	if (h <= 1.85e-278)
		tmp = Float64(Float64(t_0 * t_1) * Float64(1.0 - Float64(0.5 * Float64(Float64(h * (Float64(Float64(D / d) * Float64(0.5 * M)) ^ 2.0)) / l))));
	elseif (h <= 3e-252)
		tmp = Float64(d / sqrt(Float64(h * l)));
	elseif (h <= 1.6e+62)
		tmp = Float64(Float64(sqrt(d) / sqrt(l)) * Float64(t_0 * fma((Float64(0.5 * Float64(M / Float64(d / D))) ^ 2.0), Float64(-0.5 * Float64(h / l)), 1.0)));
	else
		tmp = Float64(Float64(t_1 * Float64(sqrt(d) / sqrt(h))) * Float64(1.0 - Float64(0.5 * Float64(Float64(0.25 * Float64(D * D)) / Float64(Float64(d * d) / Float64(Float64(M * Float64(h * M)) / l))))));
	end
	return tmp
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := Block[{t$95$0 = N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[h, 1.85e-278], N[(N[(t$95$0 * t$95$1), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(h * N[Power[N[(N[(D / d), $MachinePrecision] * N[(0.5 * M), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[h, 3e-252], N[(d / N[Sqrt[N[(h * l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[h, 1.6e+62], N[(N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[l], $MachinePrecision]), $MachinePrecision] * N[(t$95$0 * N[(N[Power[N[(0.5 * N[(M / N[(d / D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision] * N[(-0.5 * N[(h / l), $MachinePrecision]), $MachinePrecision] + 1.0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(N[(t$95$1 * N[(N[Sqrt[d], $MachinePrecision] / N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(0.25 * N[(D * D), $MachinePrecision]), $MachinePrecision] / N[(N[(d * d), $MachinePrecision] / N[(N[(M * N[(h * M), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]]

Derivation

1. Initial program 71.6%

   \[\begin{array}{l} \mathbf{if}\;h \leq 1.85 \cdot 10^{-278}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \frac{h \cdot {\left(\frac{D}{d} \cdot \left(0.5 \cdot M\right)\right)}^{2}}{\ell}\right)\\ \mathbf{elif}\;h \leq 3 \cdot 10^{-252}:\\ \;\;\;\;\frac{d}{\sqrt{h \cdot \ell}}\\ \mathbf{elif}\;h \leq 1.6 \cdot 10^{+62}:\\ \;\;\;\;\frac{\sqrt{d}}{\sqrt{\ell}} \cdot \left(\sqrt{\frac{d}{h}} \cdot \mathsf{fma}\left({\left(0.5 \cdot \frac{M}{\frac{d}{D}}\right)}^{2}, -0.5 \cdot \frac{h}{\ell}, 1\right)\right)\\ \mathbf{else}:\\ \;\;\;\;\left(\sqrt{\frac{d}{\ell}} \cdot \frac{\sqrt{d}}{\sqrt{h}}\right) \cdot \left(1 - 0.5 \cdot \frac{0.25 \cdot \left(D \cdot D\right)}{\frac{d \cdot d}{\frac{M \cdot \left(h \cdot M\right)}{\ell}}}\right)\\ \end{array} \]

Alternative 3: 69.9% accurate, 1.0× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 1.15e+122)
   (*
    (* (sqrt (/ d h)) (sqrt (/ d l)))
    (- 1.0 (* 0.5 (/ (* h (pow (* (/ D d) (* 0.5 M)) 2.0)) l))))
   (* d (* (pow h -0.5) (pow l -0.5)))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 1.15e+122) {
		tmp = (sqrt((d / h)) * sqrt((d / l))) * (1.0 - (0.5 * ((h * pow(((D / d) * (0.5 * M)), 2.0)) / l)));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= 1.15d+122) then
        tmp = (sqrt((d / h)) * sqrt((d / l))) * (1.0d0 - (0.5d0 * ((h * (((d_1 / d) * (0.5d0 * m)) ** 2.0d0)) / l)))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

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

Python

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

Julia

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

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 1.15e+122], N[(N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(h * N[Power[N[(N[(D / d), $MachinePrecision] * N[(0.5 * M), $MachinePrecision]), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 70.2%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 1.15 \cdot 10^{+122}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \frac{h \cdot {\left(\frac{D}{d} \cdot \left(0.5 \cdot M\right)\right)}^{2}}{\ell}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 4: 61.4% accurate, 1.4× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} t_0 := \left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \left(0.25 \cdot \left(\left(\left(h \cdot M\right) \cdot \frac{M}{\ell}\right) \cdot \left(\frac{D}{d} \cdot \frac{D}{d}\right)\right)\right)\right)\\ \mathbf{if}\;\ell \leq 3.9 \cdot 10^{-305}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;\ell \leq 1.22 \cdot 10^{-237}:\\ \;\;\;\;\sqrt{\frac{h}{{\ell}^{3}}} \cdot \left(-0.125 \cdot \frac{D}{\frac{\frac{\frac{d}{M}}{M}}{D}}\right)\\ \mathbf{elif}\;\ell \leq 5.5 \cdot 10^{-26}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (let* ((t_0
         (*
          (* (sqrt (/ d h)) (sqrt (/ d l)))
          (-
           1.0
           (* 0.5 (* 0.25 (* (* (* h M) (/ M l)) (* (/ D d) (/ D d)))))))))
   (if (<= l 3.9e-305)
     t_0
     (if (<= l 1.22e-237)
       (* (sqrt (/ h (pow l 3.0))) (* -0.125 (/ D (/ (/ (/ d M) M) D))))
       (if (<= l 5.5e-26) t_0 (* d (* (pow h -0.5) (pow l -0.5))))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double t_0 = (sqrt((d / h)) * sqrt((d / l))) * (1.0 - (0.5 * (0.25 * (((h * M) * (M / l)) * ((D / d) * (D / d))))));
	double tmp;
	if (l <= 3.9e-305) {
		tmp = t_0;
	} else if (l <= 1.22e-237) {
		tmp = sqrt((h / pow(l, 3.0))) * (-0.125 * (D / (((d / M) / M) / D)));
	} else if (l <= 5.5e-26) {
		tmp = t_0;
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (sqrt((d / h)) * sqrt((d / l))) * (1.0d0 - (0.5d0 * (0.25d0 * (((h * m) * (m / l)) * ((d_1 / d) * (d_1 / d))))))
    if (l <= 3.9d-305) then
        tmp = t_0
    else if (l <= 1.22d-237) then
        tmp = sqrt((h / (l ** 3.0d0))) * ((-0.125d0) * (d_1 / (((d / m) / m) / d_1)))
    else if (l <= 5.5d-26) then
        tmp = t_0
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double t_0 = (Math.sqrt((d / h)) * Math.sqrt((d / l))) * (1.0 - (0.5 * (0.25 * (((h * M) * (M / l)) * ((D / d) * (D / d))))));
	double tmp;
	if (l <= 3.9e-305) {
		tmp = t_0;
	} else if (l <= 1.22e-237) {
		tmp = Math.sqrt((h / Math.pow(l, 3.0))) * (-0.125 * (D / (((d / M) / M) / D)));
	} else if (l <= 5.5e-26) {
		tmp = t_0;
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	t_0 = (math.sqrt((d / h)) * math.sqrt((d / l))) * (1.0 - (0.5 * (0.25 * (((h * M) * (M / l)) * ((D / d) * (D / d))))))
	tmp = 0
	if l <= 3.9e-305:
		tmp = t_0
	elif l <= 1.22e-237:
		tmp = math.sqrt((h / math.pow(l, 3.0))) * (-0.125 * (D / (((d / M) / M) / D)))
	elif l <= 5.5e-26:
		tmp = t_0
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	t_0 = Float64(Float64(sqrt(Float64(d / h)) * sqrt(Float64(d / l))) * Float64(1.0 - Float64(0.5 * Float64(0.25 * Float64(Float64(Float64(h * M) * Float64(M / l)) * Float64(Float64(D / d) * Float64(D / d)))))))
	tmp = 0.0
	if (l <= 3.9e-305)
		tmp = t_0;
	elseif (l <= 1.22e-237)
		tmp = Float64(sqrt(Float64(h / (l ^ 3.0))) * Float64(-0.125 * Float64(D / Float64(Float64(Float64(d / M) / M) / D))));
	elseif (l <= 5.5e-26)
		tmp = t_0;
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

MATLAB

M = abs(M)
D = abs(D)
M, D = num2cell(sort([M, D])){:}
function tmp_2 = code(d, h, l, M, D)
	t_0 = (sqrt((d / h)) * sqrt((d / l))) * (1.0 - (0.5 * (0.25 * (((h * M) * (M / l)) * ((D / d) * (D / d))))));
	tmp = 0.0;
	if (l <= 3.9e-305)
		tmp = t_0;
	elseif (l <= 1.22e-237)
		tmp = sqrt((h / (l ^ 3.0))) * (-0.125 * (D / (((d / M) / M) / D)));
	elseif (l <= 5.5e-26)
		tmp = t_0;
	else
		tmp = d * ((h ^ -0.5) * (l ^ -0.5));
	end
	tmp_2 = tmp;
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := Block[{t$95$0 = N[(N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(0.25 * N[(N[(N[(h * M), $MachinePrecision] * N[(M / l), $MachinePrecision]), $MachinePrecision] * N[(N[(D / d), $MachinePrecision] * N[(D / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[l, 3.9e-305], t$95$0, If[LessEqual[l, 1.22e-237], N[(N[Sqrt[N[(h / N[Power[l, 3.0], $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(-0.125 * N[(D / N[(N[(N[(d / M), $MachinePrecision] / M), $MachinePrecision] / D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 5.5e-26], t$95$0, N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

1. Initial program 62.8%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 3.9 \cdot 10^{-305}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \left(0.25 \cdot \left(\left(\left(h \cdot M\right) \cdot \frac{M}{\ell}\right) \cdot \left(\frac{D}{d} \cdot \frac{D}{d}\right)\right)\right)\right)\\ \mathbf{elif}\;\ell \leq 1.22 \cdot 10^{-237}:\\ \;\;\;\;\sqrt{\frac{h}{{\ell}^{3}}} \cdot \left(-0.125 \cdot \frac{D}{\frac{\frac{\frac{d}{M}}{M}}{D}}\right)\\ \mathbf{elif}\;\ell \leq 5.5 \cdot 10^{-26}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \left(0.25 \cdot \left(\left(\left(h \cdot M\right) \cdot \frac{M}{\ell}\right) \cdot \left(\frac{D}{d} \cdot \frac{D}{d}\right)\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 5: 68.4% accurate, 1.4× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 1.55e+122)
   (*
    (* (sqrt (/ d h)) (sqrt (/ d l)))
    (- 1.0 (* 0.5 (/ (* 0.25 (* D (* (/ h d) (* M (/ (* D M) d))))) l))))
   (* d (* (pow h -0.5) (pow l -0.5)))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 1.55e+122) {
		tmp = (sqrt((d / h)) * sqrt((d / l))) * (1.0 - (0.5 * ((0.25 * (D * ((h / d) * (M * ((D * M) / d))))) / l)));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= 1.55d+122) then
        tmp = (sqrt((d / h)) * sqrt((d / l))) * (1.0d0 - (0.5d0 * ((0.25d0 * (d_1 * ((h / d) * (m * ((d_1 * m) / d))))) / l)))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 1.55e+122) {
		tmp = (Math.sqrt((d / h)) * Math.sqrt((d / l))) * (1.0 - (0.5 * ((0.25 * (D * ((h / d) * (M * ((D * M) / d))))) / l)));
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= 1.55e+122:
		tmp = (math.sqrt((d / h)) * math.sqrt((d / l))) * (1.0 - (0.5 * ((0.25 * (D * ((h / d) * (M * ((D * M) / d))))) / l)))
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

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

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 1.55e+122], N[(N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 - N[(0.5 * N[(N[(0.25 * N[(D * N[(N[(h / d), $MachinePrecision] * N[(M * N[(N[(D * M), $MachinePrecision] / d), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] / l), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 68.9%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 1.55 \cdot 10^{+122}:\\ \;\;\;\;\left(\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\right) \cdot \left(1 - 0.5 \cdot \frac{0.25 \cdot \left(D \cdot \left(\frac{h}{d} \cdot \left(M \cdot \frac{D \cdot M}{d}\right)\right)\right)}{\ell}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 6: 59.4% accurate, 1.5× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l -2.4e+113)
   (* (sqrt (/ d h)) (sqrt (/ d l)))
   (if (<= l 3e-7)
     (*
      (sqrt (* (/ d h) (/ d l)))
      (+ 1.0 (* -0.5 (* (/ h l) (pow (/ (* D (* 0.5 M)) d) 2.0)))))
     (* d (* (pow h -0.5) (pow l -0.5))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= -2.4e+113) {
		tmp = sqrt((d / h)) * sqrt((d / l));
	} else if (l <= 3e-7) {
		tmp = sqrt(((d / h) * (d / l))) * (1.0 + (-0.5 * ((h / l) * pow(((D * (0.5 * M)) / d), 2.0))));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= (-2.4d+113)) then
        tmp = sqrt((d / h)) * sqrt((d / l))
    else if (l <= 3d-7) then
        tmp = sqrt(((d / h) * (d / l))) * (1.0d0 + ((-0.5d0) * ((h / l) * (((d_1 * (0.5d0 * m)) / d) ** 2.0d0))))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

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

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= -2.4e+113:
		tmp = math.sqrt((d / h)) * math.sqrt((d / l))
	elif l <= 3e-7:
		tmp = math.sqrt(((d / h) * (d / l))) * (1.0 + (-0.5 * ((h / l) * math.pow(((D * (0.5 * M)) / d), 2.0))))
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

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

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, -2.4e+113], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 3e-7], N[(N[Sqrt[N[(N[(d / h), $MachinePrecision] * N[(d / l), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(1.0 + N[(-0.5 * N[(N[(h / l), $MachinePrecision] * N[Power[N[(N[(D * N[(0.5 * M), $MachinePrecision]), $MachinePrecision] / d), $MachinePrecision], 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Initial program 59.3%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq -2.4 \cdot 10^{+113}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 3 \cdot 10^{-7}:\\ \;\;\;\;\sqrt{\frac{d}{h} \cdot \frac{d}{\ell}} \cdot \left(1 + -0.5 \cdot \left(\frac{h}{\ell} \cdot {\left(\frac{D \cdot \left(0.5 \cdot M\right)}{d}\right)}^{2}\right)\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 7: 45.8% accurate, 1.5× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq 9.2 \cdot 10^{-293}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 2.1 \cdot 10^{-195}:\\ \;\;\;\;-0.125 \cdot \left(D \cdot \left(\sqrt{\frac{h}{{\ell}^{3}}} \cdot \left(\frac{D}{d} \cdot \left(M \cdot M\right)\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 9.2e-293)
   (* (sqrt (/ d h)) (sqrt (/ d l)))
   (if (<= l 2.1e-195)
     (* -0.125 (* D (* (sqrt (/ h (pow l 3.0))) (* (/ D d) (* M M)))))
     (* d (* (pow h -0.5) (pow l -0.5))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 9.2e-293) {
		tmp = sqrt((d / h)) * sqrt((d / l));
	} else if (l <= 2.1e-195) {
		tmp = -0.125 * (D * (sqrt((h / pow(l, 3.0))) * ((D / d) * (M * M))));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= 9.2d-293) then
        tmp = sqrt((d / h)) * sqrt((d / l))
    else if (l <= 2.1d-195) then
        tmp = (-0.125d0) * (d_1 * (sqrt((h / (l ** 3.0d0))) * ((d_1 / d) * (m * m))))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 9.2e-293) {
		tmp = Math.sqrt((d / h)) * Math.sqrt((d / l));
	} else if (l <= 2.1e-195) {
		tmp = -0.125 * (D * (Math.sqrt((h / Math.pow(l, 3.0))) * ((D / d) * (M * M))));
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= 9.2e-293:
		tmp = math.sqrt((d / h)) * math.sqrt((d / l))
	elif l <= 2.1e-195:
		tmp = -0.125 * (D * (math.sqrt((h / math.pow(l, 3.0))) * ((D / d) * (M * M))))
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	tmp = 0.0
	if (l <= 9.2e-293)
		tmp = Float64(sqrt(Float64(d / h)) * sqrt(Float64(d / l)));
	elseif (l <= 2.1e-195)
		tmp = Float64(-0.125 * Float64(D * Float64(sqrt(Float64(h / (l ^ 3.0))) * Float64(Float64(D / d) * Float64(M * M)))));
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

MATLAB

M = abs(M)
D = abs(D)
M, D = num2cell(sort([M, D])){:}
function tmp_2 = code(d, h, l, M, D)
	tmp = 0.0;
	if (l <= 9.2e-293)
		tmp = sqrt((d / h)) * sqrt((d / l));
	elseif (l <= 2.1e-195)
		tmp = -0.125 * (D * (sqrt((h / (l ^ 3.0))) * ((D / d) * (M * M))));
	else
		tmp = d * ((h ^ -0.5) * (l ^ -0.5));
	end
	tmp_2 = tmp;
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 9.2e-293], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 2.1e-195], N[(-0.125 * N[(D * N[(N[Sqrt[N[(h / N[Power[l, 3.0], $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(N[(D / d), $MachinePrecision] * N[(M * M), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Initial program 50.6%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 9.2 \cdot 10^{-293}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 2.1 \cdot 10^{-195}:\\ \;\;\;\;-0.125 \cdot \left(D \cdot \left(\sqrt{\frac{h}{{\ell}^{3}}} \cdot \left(\frac{D}{d} \cdot \left(M \cdot M\right)\right)\right)\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 8: 45.7% accurate, 1.5× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq 7.4 \cdot 10^{-293}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 1.45 \cdot 10^{-142}:\\ \;\;\;\;-0.125 \cdot \left(\frac{D \cdot D}{\frac{d}{M \cdot M}} \cdot \frac{\sqrt{h}}{{\ell}^{1.5}}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 7.4e-293)
   (* (sqrt (/ d h)) (sqrt (/ d l)))
   (if (<= l 1.45e-142)
     (* -0.125 (* (/ (* D D) (/ d (* M M))) (/ (sqrt h) (pow l 1.5))))
     (* d (* (pow h -0.5) (pow l -0.5))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 7.4e-293) {
		tmp = sqrt((d / h)) * sqrt((d / l));
	} else if (l <= 1.45e-142) {
		tmp = -0.125 * (((D * D) / (d / (M * M))) * (sqrt(h) / pow(l, 1.5)));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= 7.4d-293) then
        tmp = sqrt((d / h)) * sqrt((d / l))
    else if (l <= 1.45d-142) then
        tmp = (-0.125d0) * (((d_1 * d_1) / (d / (m * m))) * (sqrt(h) / (l ** 1.5d0)))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 7.4e-293) {
		tmp = Math.sqrt((d / h)) * Math.sqrt((d / l));
	} else if (l <= 1.45e-142) {
		tmp = -0.125 * (((D * D) / (d / (M * M))) * (Math.sqrt(h) / Math.pow(l, 1.5)));
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= 7.4e-293:
		tmp = math.sqrt((d / h)) * math.sqrt((d / l))
	elif l <= 1.45e-142:
		tmp = -0.125 * (((D * D) / (d / (M * M))) * (math.sqrt(h) / math.pow(l, 1.5)))
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	tmp = 0.0
	if (l <= 7.4e-293)
		tmp = Float64(sqrt(Float64(d / h)) * sqrt(Float64(d / l)));
	elseif (l <= 1.45e-142)
		tmp = Float64(-0.125 * Float64(Float64(Float64(D * D) / Float64(d / Float64(M * M))) * Float64(sqrt(h) / (l ^ 1.5))));
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

MATLAB

M = abs(M)
D = abs(D)
M, D = num2cell(sort([M, D])){:}
function tmp_2 = code(d, h, l, M, D)
	tmp = 0.0;
	if (l <= 7.4e-293)
		tmp = sqrt((d / h)) * sqrt((d / l));
	elseif (l <= 1.45e-142)
		tmp = -0.125 * (((D * D) / (d / (M * M))) * (sqrt(h) / (l ^ 1.5)));
	else
		tmp = d * ((h ^ -0.5) * (l ^ -0.5));
	end
	tmp_2 = tmp;
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 7.4e-293], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 1.45e-142], N[(-0.125 * N[(N[(N[(D * D), $MachinePrecision] / N[(d / N[(M * M), $MachinePrecision]), $MachinePrecision]), $MachinePrecision] * N[(N[Sqrt[h], $MachinePrecision] / N[Power[l, 1.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Initial program 51.0%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 7.4 \cdot 10^{-293}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 1.45 \cdot 10^{-142}:\\ \;\;\;\;-0.125 \cdot \left(\frac{D \cdot D}{\frac{d}{M \cdot M}} \cdot \frac{\sqrt{h}}{{\ell}^{1.5}}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 9: 46.2% accurate, 1.5× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq -2 \cdot 10^{-310}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 1.25 \cdot 10^{-176}:\\ \;\;\;\;\sqrt{\frac{h}{{\ell}^{3}}} \cdot \left(-0.125 \cdot \frac{D}{\frac{\frac{\frac{d}{M}}{M}}{D}}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l -2e-310)
   (* (sqrt (/ d h)) (sqrt (/ d l)))
   (if (<= l 1.25e-176)
     (* (sqrt (/ h (pow l 3.0))) (* -0.125 (/ D (/ (/ (/ d M) M) D))))
     (* d (* (pow h -0.5) (pow l -0.5))))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= -2e-310) {
		tmp = sqrt((d / h)) * sqrt((d / l));
	} else if (l <= 1.25e-176) {
		tmp = sqrt((h / pow(l, 3.0))) * (-0.125 * (D / (((d / M) / M) / D)));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= (-2d-310)) then
        tmp = sqrt((d / h)) * sqrt((d / l))
    else if (l <= 1.25d-176) then
        tmp = sqrt((h / (l ** 3.0d0))) * ((-0.125d0) * (d_1 / (((d / m) / m) / d_1)))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= -2e-310) {
		tmp = Math.sqrt((d / h)) * Math.sqrt((d / l));
	} else if (l <= 1.25e-176) {
		tmp = Math.sqrt((h / Math.pow(l, 3.0))) * (-0.125 * (D / (((d / M) / M) / D)));
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= -2e-310:
		tmp = math.sqrt((d / h)) * math.sqrt((d / l))
	elif l <= 1.25e-176:
		tmp = math.sqrt((h / math.pow(l, 3.0))) * (-0.125 * (D / (((d / M) / M) / D)))
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	tmp = 0.0
	if (l <= -2e-310)
		tmp = Float64(sqrt(Float64(d / h)) * sqrt(Float64(d / l)));
	elseif (l <= 1.25e-176)
		tmp = Float64(sqrt(Float64(h / (l ^ 3.0))) * Float64(-0.125 * Float64(D / Float64(Float64(Float64(d / M) / M) / D))));
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

MATLAB

M = abs(M)
D = abs(D)
M, D = num2cell(sort([M, D])){:}
function tmp_2 = code(d, h, l, M, D)
	tmp = 0.0;
	if (l <= -2e-310)
		tmp = sqrt((d / h)) * sqrt((d / l));
	elseif (l <= 1.25e-176)
		tmp = sqrt((h / (l ^ 3.0))) * (-0.125 * (D / (((d / M) / M) / D)));
	else
		tmp = d * ((h ^ -0.5) * (l ^ -0.5));
	end
	tmp_2 = tmp;
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, -2e-310], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[l, 1.25e-176], N[(N[Sqrt[N[(h / N[Power[l, 3.0], $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[(-0.125 * N[(D / N[(N[(N[(d / M), $MachinePrecision] / M), $MachinePrecision] / D), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]

Derivation

1. Initial program 51.8%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq -2 \cdot 10^{-310}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{elif}\;\ell \leq 1.25 \cdot 10^{-176}:\\ \;\;\;\;\sqrt{\frac{h}{{\ell}^{3}}} \cdot \left(-0.125 \cdot \frac{D}{\frac{\frac{\frac{d}{M}}{M}}{D}}\right)\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 10: 32.1% accurate, 1.6× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ \begin{array}{l} \mathbf{if}\;\ell \leq 7.4 \cdot 10^{-293}:\\ \;\;\;\;\frac{d}{\sqrt[3]{{\left(h \cdot \ell\right)}^{1.5}}}\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 7.4e-293)
   (/ d (cbrt (pow (* h l) 1.5)))
   (* d (* (pow h -0.5) (pow l -0.5)))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 7.4e-293) {
		tmp = d / cbrt(pow((h * l), 1.5));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 7.4e-293) {
		tmp = d / Math.cbrt(Math.pow((h * l), 1.5));
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	tmp = 0.0
	if (l <= 7.4e-293)
		tmp = Float64(d / cbrt((Float64(h * l) ^ 1.5)));
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 7.4e-293], N[(d / N[Power[N[Power[N[(h * l), $MachinePrecision], 1.5], $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 30.8%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 7.4 \cdot 10^{-293}:\\ \;\;\;\;\frac{d}{\sqrt[3]{{\left(h \cdot \ell\right)}^{1.5}}}\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 11: 45.4% accurate, 1.6× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 1.55e-296)
   (* (sqrt (/ d h)) (sqrt (/ d l)))
   (* d (* (pow h -0.5) (pow l -0.5)))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 1.55e-296) {
		tmp = sqrt((d / h)) * sqrt((d / l));
	} else {
		tmp = d * (pow(h, -0.5) * pow(l, -0.5));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= 1.55d-296) then
        tmp = sqrt((d / h)) * sqrt((d / l))
    else
        tmp = d * ((h ** (-0.5d0)) * (l ** (-0.5d0)))
    end if
    code = tmp
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 1.55e-296) {
		tmp = Math.sqrt((d / h)) * Math.sqrt((d / l));
	} else {
		tmp = d * (Math.pow(h, -0.5) * Math.pow(l, -0.5));
	}
	return tmp;
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= 1.55e-296:
		tmp = math.sqrt((d / h)) * math.sqrt((d / l))
	else:
		tmp = d * (math.pow(h, -0.5) * math.pow(l, -0.5))
	return tmp

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	tmp = 0.0
	if (l <= 1.55e-296)
		tmp = Float64(sqrt(Float64(d / h)) * sqrt(Float64(d / l)));
	else
		tmp = Float64(d * Float64((h ^ -0.5) * (l ^ -0.5)));
	end
	return tmp
end

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 1.55e-296], N[(N[Sqrt[N[(d / h), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(d / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(d * N[(N[Power[h, -0.5], $MachinePrecision] * N[Power[l, -0.5], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 45.5%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 1.55 \cdot 10^{-296}:\\ \;\;\;\;\sqrt{\frac{d}{h}} \cdot \sqrt{\frac{d}{\ell}}\\ \mathbf{else}:\\ \;\;\;\;d \cdot \left({h}^{-0.5} \cdot {\ell}^{-0.5}\right)\\ \end{array} \]

Alternative 12: 32.1% accurate, 1.6× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 7.4e-293)
   (/ d (cbrt (pow (* h l) 1.5)))
   (/ d (* (sqrt l) (sqrt h)))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 7.4e-293) {
		tmp = d / cbrt(pow((h * l), 1.5));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 7.4e-293) {
		tmp = d / Math.cbrt(Math.pow((h * l), 1.5));
	} else {
		tmp = d / (Math.sqrt(l) * Math.sqrt(h));
	}
	return tmp;
}

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	tmp = 0.0
	if (l <= 7.4e-293)
		tmp = Float64(d / cbrt((Float64(h * l) ^ 1.5)));
	else
		tmp = Float64(d / Float64(sqrt(l) * sqrt(h)));
	end
	return tmp
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 7.4e-293], N[(d / N[Power[N[Power[N[(h * l), $MachinePrecision], 1.5], $MachinePrecision], 1/3], $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 30.7%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 7.4 \cdot 10^{-293}:\\ \;\;\;\;\frac{d}{\sqrt[3]{{\left(h \cdot \ell\right)}^{1.5}}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \]

Alternative 13: 30.8% accurate, 1.6× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D)
 :precision binary64
 (if (<= l 8e-293) (/ d (sqrt (* h l))) (/ d (* (sqrt l) (sqrt h)))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	double tmp;
	if (l <= 8e-293) {
		tmp = d / sqrt((h * l));
	} else {
		tmp = d / (sqrt(l) * sqrt(h));
	}
	return tmp;
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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
    real(8) :: tmp
    if (l <= 8d-293) then
        tmp = d / sqrt((h * l))
    else
        tmp = d / (sqrt(l) * sqrt(h))
    end if
    code = tmp
end function

Java

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

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	tmp = 0
	if l <= 8e-293:
		tmp = d / math.sqrt((h * l))
	else:
		tmp = d / (math.sqrt(l) * math.sqrt(h))
	return tmp

Julia

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

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := If[LessEqual[l, 8e-293], N[(d / N[Sqrt[N[(h * l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(d / N[(N[Sqrt[l], $MachinePrecision] * N[Sqrt[h], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]

Derivation

1. Initial program 30.7%

   \[\begin{array}{l} \mathbf{if}\;\ell \leq 8 \cdot 10^{-293}:\\ \;\;\;\;\frac{d}{\sqrt{h \cdot \ell}}\\ \mathbf{else}:\\ \;\;\;\;\frac{d}{\sqrt{\ell} \cdot \sqrt{h}}\\ \end{array} \]

Alternative 14: 27.2% accurate, 3.1× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D) :precision binary64 (* d (sqrt (/ (/ 1.0 h) l))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	return d * sqrt(((1.0 / h) / l));
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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 * sqrt(((1.0d0 / h) / l))
end function

Java

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

Python

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

Julia

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

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := N[(d * N[Sqrt[N[(N[(1.0 / h), $MachinePrecision] / l), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

Derivation

1. Initial program 27.0%

   \[d \cdot \sqrt{\frac{\frac{1}{h}}{\ell}} \]

Alternative 15: 27.2% accurate, 3.1× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D) :precision binary64 (* d (sqrt (/ (/ 1.0 l) h))))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	return d * sqrt(((1.0 / l) / h));
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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 * sqrt(((1.0d0 / l) / h))
end function

Java

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

Python

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

Julia

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

MATLAB

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

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := N[(d * N[Sqrt[N[(N[(1.0 / l), $MachinePrecision] / h), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

Derivation

1. Initial program 27.0%

   \[d \cdot \sqrt{\frac{\frac{1}{\ell}}{h}} \]

Alternative 16: 27.0% accurate, 3.1× speedup

Math

\[\begin{array}{l} M = |M|\\ D = |D|\\ [M, D] = \mathsf{sort}([M, D])\\ \\ d \cdot {\left(h \cdot \ell\right)}^{-0.5} \end{array} \]

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D) :precision binary64 (* d (pow (* h l) -0.5)))

C

M = abs(M);
D = abs(D);
assert(M < D);
double code(double d, double h, double l, double M, double D) {
	return d * pow((h * l), -0.5);
}

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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 * l) ** (-0.5d0))
end function

Java

M = Math.abs(M);
D = Math.abs(D);
assert M < D;
public static double code(double d, double h, double l, double M, double D) {
	return d * Math.pow((h * l), -0.5);
}

Python

M = abs(M)
D = abs(D)
[M, D] = sort([M, D])
def code(d, h, l, M, D):
	return d * math.pow((h * l), -0.5)

Julia

M = abs(M)
D = abs(D)
M, D = sort([M, D])
function code(d, h, l, M, D)
	return Float64(d * (Float64(h * l) ^ -0.5))
end

MATLAB

M = abs(M)
D = abs(D)
M, D = num2cell(sort([M, D])){:}
function tmp = code(d, h, l, M, D)
	tmp = d * ((h * l) ^ -0.5);
end

Wolfram

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
code[d_, h_, l_, M_, D_] := N[(d * N[Power[N[(h * l), $MachinePrecision], -0.5], $MachinePrecision]), $MachinePrecision]

Derivation

1. Initial program 26.6%

   \[d \cdot {\left(h \cdot \ell\right)}^{-0.5} \]

Alternative 17: 27.0% accurate, 3.2× speedup

Math

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

FPCore

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
(FPCore (d h l M D) :precision binary64 (/ d (sqrt (* h l))))

C

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

Fortran

NOTE: M should be positive before calling this function
NOTE: D should be positive before calling this function
NOTE: M and D should be sorted in increasing order before calling this function.
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 / sqrt((h * l))
end function

Java

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

Python

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

Julia

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

MATLAB

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

Wolfram

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

Derivation

1. Initial program 26.6%

   \[\frac{d}{\sqrt{h \cdot \ell}} \]

Reproduce

herbie shell --seed 2023297 
(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)))))