math.sqrt on complex, imaginary part, im greater than 0 branch

Average Error: 38.8 → 19.3

Time: 3.6s

Precision: binary64

\[im > 0\]

Math

\[0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} - re\right)}\]

↓

\[\begin{array}{l} \mathbf{if}\;\sqrt{re \cdot re + im \cdot im} - re \leq 0:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \frac{im \cdot im}{re + \sqrt{re \cdot re + im \cdot im}}}\\ \mathbf{elif}\;\sqrt{re \cdot re + im \cdot im} - re \leq 3.7554936952845314 \cdot 10^{-197}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(im - re\right)}\\ \mathbf{elif}\;\sqrt{re \cdot re + im \cdot im} - re \leq 1.4005610381443633 \cdot 10^{-169}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re \cdot -2\right)}\\ \mathbf{elif}\;\sqrt{re \cdot re + im \cdot im} - re \leq 5.0967300915531066 \cdot 10^{+129}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(\sqrt{\sqrt{re \cdot re + im \cdot im}} \cdot \sqrt{\sqrt{re \cdot re + im \cdot im}} - re\right)}\\ \mathbf{else}:\\ \;\;\;\;0.5 \cdot \left(\sqrt{2} \cdot \sqrt{im - re}\right)\\ \end{array}\]

FPCore

(FPCore (re im)
 :precision binary64
 (* 0.5 (sqrt (* 2.0 (- (sqrt (+ (* re re) (* im im))) re)))))

↓

(FPCore (re im)
 :precision binary64
 (if (<= (- (sqrt (+ (* re re) (* im im))) re) 0.0)
   (* 0.5 (sqrt (* 2.0 (/ (* im im) (+ re (sqrt (+ (* re re) (* im im))))))))
   (if (<= (- (sqrt (+ (* re re) (* im im))) re) 3.7554936952845314e-197)
     (* 0.5 (sqrt (* 2.0 (- im re))))
     (if (<= (- (sqrt (+ (* re re) (* im im))) re) 1.4005610381443633e-169)
       (* 0.5 (sqrt (* 2.0 (* re -2.0))))
       (if (<= (- (sqrt (+ (* re re) (* im im))) re) 5.0967300915531066e+129)
         (*
          0.5
          (sqrt
           (*
            2.0
            (-
             (*
              (sqrt (sqrt (+ (* re re) (* im im))))
              (sqrt (sqrt (+ (* re re) (* im im)))))
             re))))
         (* 0.5 (* (sqrt 2.0) (sqrt (- im re)))))))))

C

double code(double re, double im) {
	return ((double) (0.5 * ((double) sqrt(((double) (2.0 * ((double) (((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))) - re))))))));
}

↓

double code(double re, double im) {
	double tmp;
	if ((((double) (((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))) - re)) <= 0.0)) {
		tmp = ((double) (0.5 * ((double) sqrt(((double) (2.0 * (((double) (im * im)) / ((double) (re + ((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))))))))))));
	} else {
		double tmp_1;
		if ((((double) (((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))) - re)) <= 3.7554936952845314e-197)) {
			tmp_1 = ((double) (0.5 * ((double) sqrt(((double) (2.0 * ((double) (im - re))))))));
		} else {
			double tmp_2;
			if ((((double) (((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))) - re)) <= 1.4005610381443633e-169)) {
				tmp_2 = ((double) (0.5 * ((double) sqrt(((double) (2.0 * ((double) (re * -2.0))))))));
			} else {
				double tmp_3;
				if ((((double) (((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))) - re)) <= 5.0967300915531066e+129)) {
					tmp_3 = ((double) (0.5 * ((double) sqrt(((double) (2.0 * ((double) (((double) (((double) sqrt(((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))))) * ((double) sqrt(((double) sqrt(((double) (((double) (re * re)) + ((double) (im * im)))))))))) - re))))))));
				} else {
					tmp_3 = ((double) (0.5 * ((double) (((double) sqrt(2.0)) * ((double) sqrt(((double) (im - re))))))));
				}
				tmp_2 = tmp_3;
			}
			tmp_1 = tmp_2;
		}
		tmp = tmp_1;
	}
	return tmp;
}

Error

Bits error versus re

Bits error versus im

Derivation

1. Split input into 5 regimes

2. if (-.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re) < 0.0

   1. Initial program 58.7

      \[0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} - re\right)}\]
   2. Using strategy rm

   3. Applied flip--_binary64 58.2

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \color{blue}{\frac{\sqrt{re \cdot re + im \cdot im} \cdot \sqrt{re \cdot re + im \cdot im} - re \cdot re}{\sqrt{re \cdot re + im \cdot im} + re}}}\]

   4. Simplified 36.0

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \frac{\color{blue}{im \cdot im}}{\sqrt{re \cdot re + im \cdot im} + re}}\]

   5. Simplified 36.0

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \frac{im \cdot im}{\color{blue}{re + \sqrt{re \cdot re + im \cdot im}}}}\]

   if 0.0 < (-.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re) < 3.7554936952845314e-197

   1. Initial program 55.3

      \[0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} - re\right)}\]

   2. Taylor expanded around 0 20.1

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \left(\color{blue}{im} - re\right)}\]

   if 3.7554936952845314e-197 < (-.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re) < 1.4005610381443633e-169

   1. Initial program 51.2

      \[0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} - re\right)}\]

   2. Taylor expanded around -inf 4.0

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \color{blue}{\left(-2 \cdot re\right)}}\]

   3. Simplified 4.0

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \color{blue}{\left(re \cdot -2\right)}}\]

   if 1.4005610381443633e-169 < (-.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re) < 5.0967300915531066e129

   1. Initial program 0.7

      \[0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} - re\right)}\]
   2. Using strategy rm

   3. Applied add-sqr-sqrt_binary64 0.7

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \left(\sqrt{\color{blue}{\sqrt{re \cdot re + im \cdot im} \cdot \sqrt{re \cdot re + im \cdot im}}} - re\right)}\]

   4. Applied sqrt-prod_binary64 0.9

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \left(\color{blue}{\sqrt{\sqrt{re \cdot re + im \cdot im}} \cdot \sqrt{\sqrt{re \cdot re + im \cdot im}}} - re\right)}\]

   if 5.0967300915531066e129 < (-.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re)

   1. Initial program 58.5

      \[0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} - re\right)}\]

   2. Taylor expanded around 0 27.2

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \left(\color{blue}{im} - re\right)}\]
   3. Using strategy rm

   4. Applied sqrt-prod_binary64 27.3

      \[\leadsto 0.5 \cdot \color{blue}{\left(\sqrt{2} \cdot \sqrt{im - re}\right)}\]
3. Recombined 5 regimes into one program.

4. Final simplification 19.3

   \[\leadsto \begin{array}{l} \mathbf{if}\;\sqrt{re \cdot re + im \cdot im} - re \leq 0:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \frac{im \cdot im}{re + \sqrt{re \cdot re + im \cdot im}}}\\ \mathbf{elif}\;\sqrt{re \cdot re + im \cdot im} - re \leq 3.7554936952845314 \cdot 10^{-197}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(im - re\right)}\\ \mathbf{elif}\;\sqrt{re \cdot re + im \cdot im} - re \leq 1.4005610381443633 \cdot 10^{-169}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re \cdot -2\right)}\\ \mathbf{elif}\;\sqrt{re \cdot re + im \cdot im} - re \leq 5.0967300915531066 \cdot 10^{+129}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(\sqrt{\sqrt{re \cdot re + im \cdot im}} \cdot \sqrt{\sqrt{re \cdot re + im \cdot im}} - re\right)}\\ \mathbf{else}:\\ \;\;\;\;0.5 \cdot \left(\sqrt{2} \cdot \sqrt{im - re}\right)\\ \end{array}\]

Reproduce

herbie shell --seed 2020210 
(FPCore (re im)
  :name "math.sqrt on complex, imaginary part, im greater than 0 branch"
  :precision binary64
  :pre (> im 0.0)
  (* 0.5 (sqrt (* 2.0 (- (sqrt (+ (* re re) (* im im))) re)))))