2cbrt (problem 3.3.4)

Average Error: 30.0 → 8.8

Time: 8.9s

Precision: binary64

Math

\[\sqrt[3]{x + 1} - \sqrt[3]{x} \]

↓

\[\begin{array}{l} t_0 := \sqrt[3]{x + 1}\\ \mathbf{if}\;x \leq -0.9949826766648939:\\ \;\;\;\;\mathsf{fma}\left(0.3333333333333333, \sqrt[3]{\frac{1}{x \cdot x}}, \mathsf{fma}\left(0.06172839506172839, \sqrt[3]{\frac{1}{{x}^{8}}}, \sqrt[3]{\frac{1}{{x}^{5}}} \cdot -0.1111111111111111\right)\right)\\ \mathbf{elif}\;x \leq 1.7226177030967663 \cdot 10^{-21}:\\ \;\;\;\;{\left(x + 1\right)}^{0.3333333333333333} - \sqrt[3]{x}\\ \mathbf{else}:\\ \;\;\;\;\frac{1}{{x}^{0.6666666666666666} + t_0 \cdot \left(\sqrt[3]{x} + t_0\right)}\\ \end{array} \]

FPCore

(FPCore (x) :precision binary64 (- (cbrt (+ x 1.0)) (cbrt x)))

↓

(FPCore (x)
 :precision binary64
 (let* ((t_0 (cbrt (+ x 1.0))))
   (if (<= x -0.9949826766648939)
     (fma
      0.3333333333333333
      (cbrt (/ 1.0 (* x x)))
      (fma
       0.06172839506172839
       (cbrt (/ 1.0 (pow x 8.0)))
       (* (cbrt (/ 1.0 (pow x 5.0))) -0.1111111111111111)))
     (if (<= x 1.7226177030967663e-21)
       (- (pow (+ x 1.0) 0.3333333333333333) (cbrt x))
       (/ 1.0 (+ (pow x 0.6666666666666666) (* t_0 (+ (cbrt x) t_0))))))))

C

double code(double x) {
	return cbrt((x + 1.0)) - cbrt(x);
}

↓

double code(double x) {
	double t_0 = cbrt((x + 1.0));
	double tmp;
	if (x <= -0.9949826766648939) {
		tmp = fma(0.3333333333333333, cbrt((1.0 / (x * x))), fma(0.06172839506172839, cbrt((1.0 / pow(x, 8.0))), (cbrt((1.0 / pow(x, 5.0))) * -0.1111111111111111)));
	} else if (x <= 1.7226177030967663e-21) {
		tmp = pow((x + 1.0), 0.3333333333333333) - cbrt(x);
	} else {
		tmp = 1.0 / (pow(x, 0.6666666666666666) + (t_0 * (cbrt(x) + t_0)));
	}
	return tmp;
}

Julia

function code(x)
	return Float64(cbrt(Float64(x + 1.0)) - cbrt(x))
end

↓

function code(x)
	t_0 = cbrt(Float64(x + 1.0))
	tmp = 0.0
	if (x <= -0.9949826766648939)
		tmp = fma(0.3333333333333333, cbrt(Float64(1.0 / Float64(x * x))), fma(0.06172839506172839, cbrt(Float64(1.0 / (x ^ 8.0))), Float64(cbrt(Float64(1.0 / (x ^ 5.0))) * -0.1111111111111111)));
	elseif (x <= 1.7226177030967663e-21)
		tmp = Float64((Float64(x + 1.0) ^ 0.3333333333333333) - cbrt(x));
	else
		tmp = Float64(1.0 / Float64((x ^ 0.6666666666666666) + Float64(t_0 * Float64(cbrt(x) + t_0))));
	end
	return tmp
end

Wolfram

code[x_] := N[(N[Power[N[(x + 1.0), $MachinePrecision], 1/3], $MachinePrecision] - N[Power[x, 1/3], $MachinePrecision]), $MachinePrecision]

↓

code[x_] := Block[{t$95$0 = N[Power[N[(x + 1.0), $MachinePrecision], 1/3], $MachinePrecision]}, If[LessEqual[x, -0.9949826766648939], N[(0.3333333333333333 * N[Power[N[(1.0 / N[(x * x), $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[(0.06172839506172839 * N[Power[N[(1.0 / N[Power[x, 8.0], $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] + N[(N[Power[N[(1.0 / N[Power[x, 5.0], $MachinePrecision]), $MachinePrecision], 1/3], $MachinePrecision] * -0.1111111111111111), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[x, 1.7226177030967663e-21], N[(N[Power[N[(x + 1.0), $MachinePrecision], 0.3333333333333333], $MachinePrecision] - N[Power[x, 1/3], $MachinePrecision]), $MachinePrecision], N[(1.0 / N[(N[Power[x, 0.6666666666666666], $MachinePrecision] + N[(t$95$0 * N[(N[Power[x, 1/3], $MachinePrecision] + t$95$0), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]

Error

Bits error versus x

Derivation

1. Split input into 3 regimes

2. if x < -0.99498267666489393

   1. Initial program 59.3

      \[\sqrt[3]{x + 1} - \sqrt[3]{x} \]

   2. Taylor expanded in x around inf 44.3

      \[\leadsto \color{blue}{\left(0.3333333333333333 \cdot {\left(\frac{1}{{x}^{2}}\right)}^{0.3333333333333333} + 0.06172839506172839 \cdot {\left(\frac{1}{{x}^{8}}\right)}^{0.3333333333333333}\right) - 0.1111111111111111 \cdot {\left(\frac{1}{{x}^{5}}\right)}^{0.3333333333333333}} \]

   3. Simplified 30.8

      \[\leadsto \color{blue}{\mathsf{fma}\left(0.3333333333333333, \sqrt[3]{\frac{1}{x \cdot x}}, \mathsf{fma}\left(0.06172839506172839, \sqrt[3]{\frac{1}{{x}^{8}}}, \sqrt[3]{\frac{1}{{x}^{5}}} \cdot -0.1111111111111111\right)\right)} \]

   if -0.99498267666489393 < x < 1.7226177030967663e-21

   1. Initial program 0.0

      \[\sqrt[3]{x + 1} - \sqrt[3]{x} \]

   2. Applied pow1/3_binary64 0.0

      \[\leadsto \color{blue}{{\left(x + 1\right)}^{0.3333333333333333}} - \sqrt[3]{x} \]

   if 1.7226177030967663e-21 < x

   1. Initial program 55.9

      \[\sqrt[3]{x + 1} - \sqrt[3]{x} \]

   2. Applied flip3--_binary64 55.8

      \[\leadsto \color{blue}{\frac{{\left(\sqrt[3]{x + 1}\right)}^{3} - {\left(\sqrt[3]{x}\right)}^{3}}{\sqrt[3]{x + 1} \cdot \sqrt[3]{x + 1} + \left(\sqrt[3]{x} \cdot \sqrt[3]{x} + \sqrt[3]{x + 1} \cdot \sqrt[3]{x}\right)}} \]

   3. Simplified 1.0

      \[\leadsto \frac{\color{blue}{1}}{\sqrt[3]{x + 1} \cdot \sqrt[3]{x + 1} + \left(\sqrt[3]{x} \cdot \sqrt[3]{x} + \sqrt[3]{x + 1} \cdot \sqrt[3]{x}\right)} \]

   4. Simplified 4.2

      \[\leadsto \frac{1}{\color{blue}{{x}^{0.6666666666666666} + \sqrt[3]{1 + x} \cdot \left(\sqrt[3]{x} + \sqrt[3]{1 + x}\right)}} \]
3. Recombined 3 regimes into one program.

4. Final simplification 8.8

   \[\leadsto \begin{array}{l} \mathbf{if}\;x \leq -0.9949826766648939:\\ \;\;\;\;\mathsf{fma}\left(0.3333333333333333, \sqrt[3]{\frac{1}{x \cdot x}}, \mathsf{fma}\left(0.06172839506172839, \sqrt[3]{\frac{1}{{x}^{8}}}, \sqrt[3]{\frac{1}{{x}^{5}}} \cdot -0.1111111111111111\right)\right)\\ \mathbf{elif}\;x \leq 1.7226177030967663 \cdot 10^{-21}:\\ \;\;\;\;{\left(x + 1\right)}^{0.3333333333333333} - \sqrt[3]{x}\\ \mathbf{else}:\\ \;\;\;\;\frac{1}{{x}^{0.6666666666666666} + \sqrt[3]{x + 1} \cdot \left(\sqrt[3]{x} + \sqrt[3]{x + 1}\right)}\\ \end{array} \]

Reproduce

herbie shell --seed 2022131 
(FPCore (x)
  :name "2cbrt (problem 3.3.4)"
  :precision binary64
  (- (cbrt (+ x 1.0)) (cbrt x)))