Falkner and Boettcher, Appendix A

Average Error: 2.2 → 0.3

Time: 7.0s

Precision: binary64

Math

\[\frac{a \cdot {k}^{m}}{\left(1 + 10 \cdot k\right) + k \cdot k} \]

↓

\[\begin{array}{l} t_0 := \frac{a \cdot {k}^{m}}{\left(1 + k \cdot 10\right) + k \cdot k}\\ \mathbf{if}\;t_0 \leq -9.662942266681646 \cdot 10^{-284}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;t_0 \leq 0:\\ \;\;\;\;\frac{a \cdot \frac{{k}^{m}}{k}}{k}\\ \mathbf{else}:\\ \;\;\;\;t_0\\ \end{array} \]

FPCore

(FPCore (a k m)
 :precision binary64
 (/ (* a (pow k m)) (+ (+ 1.0 (* 10.0 k)) (* k k))))

↓

(FPCore (a k m)
 :precision binary64
 (let* ((t_0 (/ (* a (pow k m)) (+ (+ 1.0 (* k 10.0)) (* k k)))))
   (if (<= t_0 -9.662942266681646e-284)
     t_0
     (if (<= t_0 0.0) (/ (* a (/ (pow k m) k)) k) t_0))))

C

double code(double a, double k, double m) {
	return (a * pow(k, m)) / ((1.0 + (10.0 * k)) + (k * k));
}

↓

double code(double a, double k, double m) {
	double t_0 = (a * pow(k, m)) / ((1.0 + (k * 10.0)) + (k * k));
	double tmp;
	if (t_0 <= -9.662942266681646e-284) {
		tmp = t_0;
	} else if (t_0 <= 0.0) {
		tmp = (a * (pow(k, m) / k)) / k;
	} else {
		tmp = t_0;
	}
	return tmp;
}

Fortran

real(8) function code(a, k, m)
    real(8), intent (in) :: a
    real(8), intent (in) :: k
    real(8), intent (in) :: m
    code = (a * (k ** m)) / ((1.0d0 + (10.0d0 * k)) + (k * k))
end function

↓

real(8) function code(a, k, m)
    real(8), intent (in) :: a
    real(8), intent (in) :: k
    real(8), intent (in) :: m
    real(8) :: t_0
    real(8) :: tmp
    t_0 = (a * (k ** m)) / ((1.0d0 + (k * 10.0d0)) + (k * k))
    if (t_0 <= (-9.662942266681646d-284)) then
        tmp = t_0
    else if (t_0 <= 0.0d0) then
        tmp = (a * ((k ** m) / k)) / k
    else
        tmp = t_0
    end if
    code = tmp
end function

Java

public static double code(double a, double k, double m) {
	return (a * Math.pow(k, m)) / ((1.0 + (10.0 * k)) + (k * k));
}

↓

public static double code(double a, double k, double m) {
	double t_0 = (a * Math.pow(k, m)) / ((1.0 + (k * 10.0)) + (k * k));
	double tmp;
	if (t_0 <= -9.662942266681646e-284) {
		tmp = t_0;
	} else if (t_0 <= 0.0) {
		tmp = (a * (Math.pow(k, m) / k)) / k;
	} else {
		tmp = t_0;
	}
	return tmp;
}

Python

def code(a, k, m):
	return (a * math.pow(k, m)) / ((1.0 + (10.0 * k)) + (k * k))

↓

def code(a, k, m):
	t_0 = (a * math.pow(k, m)) / ((1.0 + (k * 10.0)) + (k * k))
	tmp = 0
	if t_0 <= -9.662942266681646e-284:
		tmp = t_0
	elif t_0 <= 0.0:
		tmp = (a * (math.pow(k, m) / k)) / k
	else:
		tmp = t_0
	return tmp

Julia

function code(a, k, m)
	return Float64(Float64(a * (k ^ m)) / Float64(Float64(1.0 + Float64(10.0 * k)) + Float64(k * k)))
end

↓

function code(a, k, m)
	t_0 = Float64(Float64(a * (k ^ m)) / Float64(Float64(1.0 + Float64(k * 10.0)) + Float64(k * k)))
	tmp = 0.0
	if (t_0 <= -9.662942266681646e-284)
		tmp = t_0;
	elseif (t_0 <= 0.0)
		tmp = Float64(Float64(a * Float64((k ^ m) / k)) / k);
	else
		tmp = t_0;
	end
	return tmp
end

MATLAB

function tmp = code(a, k, m)
	tmp = (a * (k ^ m)) / ((1.0 + (10.0 * k)) + (k * k));
end

↓

function tmp_2 = code(a, k, m)
	t_0 = (a * (k ^ m)) / ((1.0 + (k * 10.0)) + (k * k));
	tmp = 0.0;
	if (t_0 <= -9.662942266681646e-284)
		tmp = t_0;
	elseif (t_0 <= 0.0)
		tmp = (a * ((k ^ m) / k)) / k;
	else
		tmp = t_0;
	end
	tmp_2 = tmp;
end

Wolfram

code[a_, k_, m_] := N[(N[(a * N[Power[k, m], $MachinePrecision]), $MachinePrecision] / N[(N[(1.0 + N[(10.0 * k), $MachinePrecision]), $MachinePrecision] + N[(k * k), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]

↓

code[a_, k_, m_] := Block[{t$95$0 = N[(N[(a * N[Power[k, m], $MachinePrecision]), $MachinePrecision] / N[(N[(1.0 + N[(k * 10.0), $MachinePrecision]), $MachinePrecision] + N[(k * k), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]}, If[LessEqual[t$95$0, -9.662942266681646e-284], t$95$0, If[LessEqual[t$95$0, 0.0], N[(N[(a * N[(N[Power[k, m], $MachinePrecision] / k), $MachinePrecision]), $MachinePrecision] / k), $MachinePrecision], t$95$0]]]

Error

Bits error versus a

Bits error versus k

Bits error versus m

Derivation

1. Split input into 2 regimes

2. if (/.f64 (*.f64 a (pow.f64 k m)) (+.f64 (+.f64 1 (*.f64 10 k)) (*.f64 k k))) < -9.6629422666816458e-284 or 0.0 < (/.f64 (*.f64 a (pow.f64 k m)) (+.f64 (+.f64 1 (*.f64 10 k)) (*.f64 k k)))

   1. Initial program 0.2

      \[\frac{a \cdot {k}^{m}}{\left(1 + 10 \cdot k\right) + k \cdot k} \]

   if -9.6629422666816458e-284 < (/.f64 (*.f64 a (pow.f64 k m)) (+.f64 (+.f64 1 (*.f64 10 k)) (*.f64 k k))) < 0.0

   1. Initial program 3.2

      \[\frac{a \cdot {k}^{m}}{\left(1 + 10 \cdot k\right) + k \cdot k} \]

   2. Simplified 3.2

      \[\leadsto \color{blue}{\frac{a \cdot {k}^{m}}{\mathsf{fma}\left(k, k + 10, 1\right)}} \]

   3. Taylor expanded in k around inf 33.2

      \[\leadsto \color{blue}{\frac{a \cdot e^{-1 \cdot \left(\log \left(\frac{1}{k}\right) \cdot m\right)}}{{k}^{2}}} \]

   4. Simplified 15.5

      \[\leadsto \color{blue}{\frac{a \cdot {k}^{m}}{k \cdot k}} \]

   5. Applied associate-/r*_binary64 0.4

      \[\leadsto \color{blue}{\frac{\frac{a \cdot {k}^{m}}{k}}{k}} \]

   6. Simplified 0.4

      \[\leadsto \frac{\color{blue}{a \cdot \frac{{k}^{m}}{k}}}{k} \]
3. Recombined 2 regimes into one program.

4. Final simplification 0.3

   \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{a \cdot {k}^{m}}{\left(1 + k \cdot 10\right) + k \cdot k} \leq -9.662942266681646 \cdot 10^{-284}:\\ \;\;\;\;\frac{a \cdot {k}^{m}}{\left(1 + k \cdot 10\right) + k \cdot k}\\ \mathbf{elif}\;\frac{a \cdot {k}^{m}}{\left(1 + k \cdot 10\right) + k \cdot k} \leq 0:\\ \;\;\;\;\frac{a \cdot \frac{{k}^{m}}{k}}{k}\\ \mathbf{else}:\\ \;\;\;\;\frac{a \cdot {k}^{m}}{\left(1 + k \cdot 10\right) + k \cdot k}\\ \end{array} \]

Reproduce

herbie shell --seed 2022129 
(FPCore (a k m)
  :name "Falkner and Boettcher, Appendix A"
  :precision binary64
  (/ (* a (pow k m)) (+ (+ 1.0 (* 10.0 k)) (* k k))))