2isqrt (example 3.6)

Average Error: 19.4 → 5.3

Time: 8.5s

Precision: binary64

Math

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

↓

\[\begin{array}{l} t_0 := {\left(1 + x\right)}^{-0.5}\\ t_1 := {\left(1 + x\right)}^{-0.25}\\ t_2 := t_1 \cdot \sqrt[3]{t_1}\\ t_3 := t_2 \cdot \sqrt[3]{t_0}\\ \mathbf{if}\;\frac{1}{\sqrt{x}} - \frac{1}{\sqrt{1 + x}} \leq 4.621890052781269 \cdot 10^{-7}:\\ \;\;\;\;\frac{\left(\frac{1.5}{x \cdot x} + \frac{1.484375}{{x}^{4}}\right) - \left(\frac{1.5}{{x}^{3}} + \frac{1.46875}{{x}^{5}}\right)}{{x}^{-0.5} + \left(t_0 + {\left(x \cdot \left(1 + x\right)\right)}^{-0.25}\right)}\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(1, {x}^{-0.5}, -t_3\right) + \mathsf{fma}\left(-t_2, e^{-0.16666666666666666 \cdot \mathsf{log1p}\left(x\right)}, t_3\right)\\ \end{array} \]

FPCore

(FPCore (x) :precision binary64 (- (/ 1.0 (sqrt x)) (/ 1.0 (sqrt (+ x 1.0)))))

↓

(FPCore (x)
 :precision binary64
 (let* ((t_0 (pow (+ 1.0 x) -0.5))
        (t_1 (pow (+ 1.0 x) -0.25))
        (t_2 (* t_1 (cbrt t_1)))
        (t_3 (* t_2 (cbrt t_0))))
   (if (<= (- (/ 1.0 (sqrt x)) (/ 1.0 (sqrt (+ 1.0 x)))) 4.621890052781269e-7)
     (/
      (-
       (+ (/ 1.5 (* x x)) (/ 1.484375 (pow x 4.0)))
       (+ (/ 1.5 (pow x 3.0)) (/ 1.46875 (pow x 5.0))))
      (+ (pow x -0.5) (+ t_0 (pow (* x (+ 1.0 x)) -0.25))))
     (+
      (fma 1.0 (pow x -0.5) (- t_3))
      (fma (- t_2) (exp (* -0.16666666666666666 (log1p x))) t_3)))))

C

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

↓

double code(double x) {
	double t_0 = pow((1.0 + x), -0.5);
	double t_1 = pow((1.0 + x), -0.25);
	double t_2 = t_1 * cbrt(t_1);
	double t_3 = t_2 * cbrt(t_0);
	double tmp;
	if (((1.0 / sqrt(x)) - (1.0 / sqrt((1.0 + x)))) <= 4.621890052781269e-7) {
		tmp = (((1.5 / (x * x)) + (1.484375 / pow(x, 4.0))) - ((1.5 / pow(x, 3.0)) + (1.46875 / pow(x, 5.0)))) / (pow(x, -0.5) + (t_0 + pow((x * (1.0 + x)), -0.25)));
	} else {
		tmp = fma(1.0, pow(x, -0.5), -t_3) + fma(-t_2, exp((-0.16666666666666666 * log1p(x))), t_3);
	}
	return tmp;
}

Error

Bits error versus x

Target

Original	19.4
Target	0.7
Herbie	5.3

\[\frac{1}{\left(x + 1\right) \cdot \sqrt{x} + x \cdot \sqrt{x + 1}} \]

Derivation

1. Split input into 2 regimes

2. if (-.f64 (/.f64 1 (sqrt.f64 x)) (/.f64 1 (sqrt.f64 (+.f64 x 1)))) < 4.6218900527813e-7

   1. Initial program 38.8

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

   2. Applied egg-rr 38.8

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

   3. Applied egg-rr 38.7

      \[\leadsto \color{blue}{\frac{\left({x}^{-0.75} - {\left(x + 1\right)}^{-0.75}\right) \cdot \left({x}^{-0.25} + {\left(x + 1\right)}^{-0.25}\right)}{{x}^{-0.5} + \left({\left(x + 1\right)}^{-0.5} + {\left(x \cdot \left(x + 1\right)\right)}^{-0.25}\right)}} \]

   4. Taylor expanded in x around inf 10.6

      \[\leadsto \frac{\color{blue}{\left(1.484375 \cdot \frac{1}{{x}^{4}} + 1.5 \cdot \frac{1}{{x}^{2}}\right) - \left(1.46875 \cdot \frac{1}{{x}^{5}} + 1.5 \cdot \frac{1}{{x}^{3}}\right)}}{{x}^{-0.5} + \left({\left(x + 1\right)}^{-0.5} + {\left(x \cdot \left(x + 1\right)\right)}^{-0.25}\right)} \]

   5. Simplified 10.5

      \[\leadsto \frac{\color{blue}{\left(\frac{1.5}{x \cdot x} + \frac{1.484375}{{x}^{4}}\right) - \left(\frac{1.5}{{x}^{3}} + \frac{1.46875}{{x}^{5}}\right)}}{{x}^{-0.5} + \left({\left(x + 1\right)}^{-0.5} + {\left(x \cdot \left(x + 1\right)\right)}^{-0.25}\right)} \]

   if 4.6218900527813e-7 < (-.f64 (/.f64 1 (sqrt.f64 x)) (/.f64 1 (sqrt.f64 (+.f64 x 1))))

   1. Initial program 0.3

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

   2. Applied egg-rr 0.6

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

   3. Applied egg-rr 0.1

      \[\leadsto \color{blue}{\mathsf{fma}\left(1, {x}^{-0.5}, -\left(\sqrt[3]{{\left(x + 1\right)}^{-0.25}} \cdot {\left(x + 1\right)}^{-0.25}\right) \cdot \sqrt[3]{{\left(x + 1\right)}^{-0.5}}\right) + \mathsf{fma}\left(-\sqrt[3]{{\left(x + 1\right)}^{-0.25}} \cdot {\left(x + 1\right)}^{-0.25}, \sqrt[3]{{\left(x + 1\right)}^{-0.5}}, \left(\sqrt[3]{{\left(x + 1\right)}^{-0.25}} \cdot {\left(x + 1\right)}^{-0.25}\right) \cdot \sqrt[3]{{\left(x + 1\right)}^{-0.5}}\right)} \]

   4. Applied egg-rr 0.1

      \[\leadsto \mathsf{fma}\left(1, {x}^{-0.5}, -\left(\sqrt[3]{{\left(x + 1\right)}^{-0.25}} \cdot {\left(x + 1\right)}^{-0.25}\right) \cdot \sqrt[3]{{\left(x + 1\right)}^{-0.5}}\right) + \mathsf{fma}\left(-\sqrt[3]{{\left(x + 1\right)}^{-0.25}} \cdot {\left(x + 1\right)}^{-0.25}, \color{blue}{e^{-0.16666666666666666 \cdot \mathsf{log1p}\left(x\right)}}, \left(\sqrt[3]{{\left(x + 1\right)}^{-0.25}} \cdot {\left(x + 1\right)}^{-0.25}\right) \cdot \sqrt[3]{{\left(x + 1\right)}^{-0.5}}\right) \]
3. Recombined 2 regimes into one program.

4. Final simplification 5.3

   \[\leadsto \begin{array}{l} \mathbf{if}\;\frac{1}{\sqrt{x}} - \frac{1}{\sqrt{1 + x}} \leq 4.621890052781269 \cdot 10^{-7}:\\ \;\;\;\;\frac{\left(\frac{1.5}{x \cdot x} + \frac{1.484375}{{x}^{4}}\right) - \left(\frac{1.5}{{x}^{3}} + \frac{1.46875}{{x}^{5}}\right)}{{x}^{-0.5} + \left({\left(1 + x\right)}^{-0.5} + {\left(x \cdot \left(1 + x\right)\right)}^{-0.25}\right)}\\ \mathbf{else}:\\ \;\;\;\;\mathsf{fma}\left(1, {x}^{-0.5}, -\left({\left(1 + x\right)}^{-0.25} \cdot \sqrt[3]{{\left(1 + x\right)}^{-0.25}}\right) \cdot \sqrt[3]{{\left(1 + x\right)}^{-0.5}}\right) + \mathsf{fma}\left(-{\left(1 + x\right)}^{-0.25} \cdot \sqrt[3]{{\left(1 + x\right)}^{-0.25}}, e^{-0.16666666666666666 \cdot \mathsf{log1p}\left(x\right)}, \left({\left(1 + x\right)}^{-0.25} \cdot \sqrt[3]{{\left(1 + x\right)}^{-0.25}}\right) \cdot \sqrt[3]{{\left(1 + x\right)}^{-0.5}}\right)\\ \end{array} \]

Reproduce

herbie shell --seed 2022127 
(FPCore (x)
  :name "2isqrt (example 3.6)"
  :precision binary64

  :herbie-target
  (/ 1.0 (+ (* (+ x 1.0) (sqrt x)) (* x (sqrt (+ x 1.0)))))

  (- (/ 1.0 (sqrt x)) (/ 1.0 (sqrt (+ x 1.0)))))