2nthrt (problem 3.4.6)

Average Error: 33.3 → 7.8

Time: 15.8s

Precision: binary64

Math

\[{\left(x + 1\right)}^{\left(\frac{1}{n}\right)} - {x}^{\left(\frac{1}{n}\right)} \]

↓

\[\begin{array}{l} t_0 := e^{\frac{\mathsf{log1p}\left(x\right)}{n}}\\ t_1 := \frac{\log \left(x + 1\right) - \log x}{n}\\ t_2 := {x}^{\left(\frac{1}{n}\right)}\\ \mathbf{if}\;x \leq 2.2 \cdot 10^{-258}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;x \leq 8.5 \cdot 10^{-251}:\\ \;\;\;\;\left({t_0}^{3} - {x}^{\left(3 \cdot \frac{1}{n}\right)}\right) \cdot \frac{1}{\mathsf{fma}\left(t_2, t_0 + t_2, {t_0}^{2}\right)}\\ \mathbf{elif}\;x \leq 100000000:\\ \;\;\;\;t_1\\ \mathbf{else}:\\ \;\;\;\;\frac{e^{\frac{\log x}{n}}}{x \cdot n}\\ \end{array} \]

FPCore

(FPCore (x n)
 :precision binary64
 (- (pow (+ x 1.0) (/ 1.0 n)) (pow x (/ 1.0 n))))

↓

(FPCore (x n)
 :precision binary64
 (let* ((t_0 (exp (/ (log1p x) n)))
        (t_1 (/ (- (log (+ x 1.0)) (log x)) n))
        (t_2 (pow x (/ 1.0 n))))
   (if (<= x 2.2e-258)
     t_1
     (if (<= x 8.5e-251)
       (*
        (- (pow t_0 3.0) (pow x (* 3.0 (/ 1.0 n))))
        (/ 1.0 (fma t_2 (+ t_0 t_2) (pow t_0 2.0))))
       (if (<= x 100000000.0) t_1 (/ (exp (/ (log x) n)) (* x n)))))))

C

double code(double x, double n) {
	return pow((x + 1.0), (1.0 / n)) - pow(x, (1.0 / n));
}

↓

double code(double x, double n) {
	double t_0 = exp((log1p(x) / n));
	double t_1 = (log((x + 1.0)) - log(x)) / n;
	double t_2 = pow(x, (1.0 / n));
	double tmp;
	if (x <= 2.2e-258) {
		tmp = t_1;
	} else if (x <= 8.5e-251) {
		tmp = (pow(t_0, 3.0) - pow(x, (3.0 * (1.0 / n)))) * (1.0 / fma(t_2, (t_0 + t_2), pow(t_0, 2.0)));
	} else if (x <= 100000000.0) {
		tmp = t_1;
	} else {
		tmp = exp((log(x) / n)) / (x * n);
	}
	return tmp;
}

Julia

function code(x, n)
	return Float64((Float64(x + 1.0) ^ Float64(1.0 / n)) - (x ^ Float64(1.0 / n)))
end

↓

function code(x, n)
	t_0 = exp(Float64(log1p(x) / n))
	t_1 = Float64(Float64(log(Float64(x + 1.0)) - log(x)) / n)
	t_2 = x ^ Float64(1.0 / n)
	tmp = 0.0
	if (x <= 2.2e-258)
		tmp = t_1;
	elseif (x <= 8.5e-251)
		tmp = Float64(Float64((t_0 ^ 3.0) - (x ^ Float64(3.0 * Float64(1.0 / n)))) * Float64(1.0 / fma(t_2, Float64(t_0 + t_2), (t_0 ^ 2.0))));
	elseif (x <= 100000000.0)
		tmp = t_1;
	else
		tmp = Float64(exp(Float64(log(x) / n)) / Float64(x * n));
	end
	return tmp
end

Wolfram

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

↓

code[x_, n_] := Block[{t$95$0 = N[Exp[N[(N[Log[1 + x], $MachinePrecision] / n), $MachinePrecision]], $MachinePrecision]}, Block[{t$95$1 = N[(N[(N[Log[N[(x + 1.0), $MachinePrecision]], $MachinePrecision] - N[Log[x], $MachinePrecision]), $MachinePrecision] / n), $MachinePrecision]}, Block[{t$95$2 = N[Power[x, N[(1.0 / n), $MachinePrecision]], $MachinePrecision]}, If[LessEqual[x, 2.2e-258], t$95$1, If[LessEqual[x, 8.5e-251], N[(N[(N[Power[t$95$0, 3.0], $MachinePrecision] - N[Power[x, N[(3.0 * N[(1.0 / n), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision] * N[(1.0 / N[(t$95$2 * N[(t$95$0 + t$95$2), $MachinePrecision] + N[Power[t$95$0, 2.0], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[x, 100000000.0], t$95$1, N[(N[Exp[N[(N[Log[x], $MachinePrecision] / n), $MachinePrecision]], $MachinePrecision] / N[(x * n), $MachinePrecision]), $MachinePrecision]]]]]]]

Error

Bits error versus x

Bits error versus n

Derivation

1. Split input into 3 regimes

2. if x < 2.20000000000000015e-258 or 8.49999999999999984e-251 < x < 1e8

   1. Initial program 46.8

      \[{\left(x + 1\right)}^{\left(\frac{1}{n}\right)} - {x}^{\left(\frac{1}{n}\right)} \]

   2. Taylor expanded in n around -inf 14.4

      \[\leadsto \color{blue}{-1 \cdot \frac{-1 \cdot \log \left(1 + x\right) - -1 \cdot \log x}{n}} \]

   if 2.20000000000000015e-258 < x < 8.49999999999999984e-251

   1. Initial program 42.1

      \[{\left(x + 1\right)}^{\left(\frac{1}{n}\right)} - {x}^{\left(\frac{1}{n}\right)} \]

   2. Applied egg-rr 41.4

      \[\leadsto \color{blue}{\left({\left(e^{\frac{\mathsf{log1p}\left(x\right)}{n}}\right)}^{3} - {x}^{\left(\frac{1}{n} \cdot 3\right)}\right) \cdot \frac{1}{\mathsf{fma}\left({x}^{\left(\frac{1}{n}\right)}, e^{\frac{\mathsf{log1p}\left(x\right)}{n}} + {x}^{\left(\frac{1}{n}\right)}, {\left(e^{\frac{\mathsf{log1p}\left(x\right)}{n}}\right)}^{2}\right)}} \]

   if 1e8 < x

   1. Initial program 21.4

      \[{\left(x + 1\right)}^{\left(\frac{1}{n}\right)} - {x}^{\left(\frac{1}{n}\right)} \]

   2. Taylor expanded in x around inf 1.3

      \[\leadsto \color{blue}{\frac{e^{-1 \cdot \frac{\log \left(\frac{1}{x}\right)}{n}}}{n \cdot x}} \]

   3. Simplified 1.3

      \[\leadsto \color{blue}{\frac{e^{\frac{\log x}{n}}}{x \cdot n}} \]
3. Recombined 3 regimes into one program.

4. Final simplification 7.8

   \[\leadsto \begin{array}{l} \mathbf{if}\;x \leq 2.2 \cdot 10^{-258}:\\ \;\;\;\;\frac{\log \left(x + 1\right) - \log x}{n}\\ \mathbf{elif}\;x \leq 8.5 \cdot 10^{-251}:\\ \;\;\;\;\left({\left(e^{\frac{\mathsf{log1p}\left(x\right)}{n}}\right)}^{3} - {x}^{\left(3 \cdot \frac{1}{n}\right)}\right) \cdot \frac{1}{\mathsf{fma}\left({x}^{\left(\frac{1}{n}\right)}, e^{\frac{\mathsf{log1p}\left(x\right)}{n}} + {x}^{\left(\frac{1}{n}\right)}, {\left(e^{\frac{\mathsf{log1p}\left(x\right)}{n}}\right)}^{2}\right)}\\ \mathbf{elif}\;x \leq 100000000:\\ \;\;\;\;\frac{\log \left(x + 1\right) - \log x}{n}\\ \mathbf{else}:\\ \;\;\;\;\frac{e^{\frac{\log x}{n}}}{x \cdot n}\\ \end{array} \]

Reproduce

herbie shell --seed 2022172 
(FPCore (x n)
  :name "2nthrt (problem 3.4.6)"
  :precision binary64
  (- (pow (+ x 1.0) (/ 1.0 n)) (pow x (/ 1.0 n))))