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

Average Error: 38.9 → 9.2

Time: 10.8s

Precision: binary64

Cost: 39692

\[im > 0\]

Math

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

↓

\[\begin{array}{l} t_0 := \mathsf{hypot}\left(re, im\right) - re\\ \mathbf{if}\;re \leq 3.322933833770276 \cdot 10^{+20}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot t_0}\\ \mathbf{elif}\;re \leq 4.926861461571793 \cdot 10^{+88}:\\ \;\;\;\;0.5 \cdot \frac{1}{\frac{\sqrt{re}}{im}}\\ \mathbf{elif}\;re \leq 3.680589009025248 \cdot 10^{+151}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(t_0 + \mathsf{fma}\left(-\sqrt[3]{re}, {\left(\sqrt[3]{re}\right)}^{2}, re\right)\right)}\\ \mathbf{else}:\\ \;\;\;\;0.5 \cdot \left(im \cdot \sqrt{\frac{1}{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
 (let* ((t_0 (- (hypot re im) re)))
   (if (<= re 3.322933833770276e+20)
     (* 0.5 (sqrt (* 2.0 t_0)))
     (if (<= re 4.926861461571793e+88)
       (* 0.5 (/ 1.0 (/ (sqrt re) im)))
       (if (<= re 3.680589009025248e+151)
         (*
          0.5
          (sqrt (* 2.0 (+ t_0 (fma (- (cbrt re)) (pow (cbrt re) 2.0) re)))))
         (* 0.5 (* im (sqrt (/ 1.0 re)))))))))

C

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

↓

double code(double re, double im) {
	double t_0 = hypot(re, im) - re;
	double tmp;
	if (re <= 3.322933833770276e+20) {
		tmp = 0.5 * sqrt((2.0 * t_0));
	} else if (re <= 4.926861461571793e+88) {
		tmp = 0.5 * (1.0 / (sqrt(re) / im));
	} else if (re <= 3.680589009025248e+151) {
		tmp = 0.5 * sqrt((2.0 * (t_0 + fma(-cbrt(re), pow(cbrt(re), 2.0), re))));
	} else {
		tmp = 0.5 * (im * sqrt((1.0 / re)));
	}
	return tmp;
}

Julia

function code(re, im)
	return Float64(0.5 * sqrt(Float64(2.0 * Float64(sqrt(Float64(Float64(re * re) + Float64(im * im))) - re))))
end

↓

function code(re, im)
	t_0 = Float64(hypot(re, im) - re)
	tmp = 0.0
	if (re <= 3.322933833770276e+20)
		tmp = Float64(0.5 * sqrt(Float64(2.0 * t_0)));
	elseif (re <= 4.926861461571793e+88)
		tmp = Float64(0.5 * Float64(1.0 / Float64(sqrt(re) / im)));
	elseif (re <= 3.680589009025248e+151)
		tmp = Float64(0.5 * sqrt(Float64(2.0 * Float64(t_0 + fma(Float64(-cbrt(re)), (cbrt(re) ^ 2.0), re)))));
	else
		tmp = Float64(0.5 * Float64(im * sqrt(Float64(1.0 / re))));
	end
	return tmp
end

Wolfram

code[re_, im_] := N[(0.5 * N[Sqrt[N[(2.0 * N[(N[Sqrt[N[(N[(re * re), $MachinePrecision] + N[(im * im), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] - re), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]

↓

code[re_, im_] := Block[{t$95$0 = N[(N[Sqrt[re ^ 2 + im ^ 2], $MachinePrecision] - re), $MachinePrecision]}, If[LessEqual[re, 3.322933833770276e+20], N[(0.5 * N[Sqrt[N[(2.0 * t$95$0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[re, 4.926861461571793e+88], N[(0.5 * N[(1.0 / N[(N[Sqrt[re], $MachinePrecision] / im), $MachinePrecision]), $MachinePrecision]), $MachinePrecision], If[LessEqual[re, 3.680589009025248e+151], N[(0.5 * N[Sqrt[N[(2.0 * N[(t$95$0 + N[((-N[Power[re, 1/3], $MachinePrecision]) * N[Power[N[Power[re, 1/3], $MachinePrecision], 2.0], $MachinePrecision] + re), $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(0.5 * N[(im * N[Sqrt[N[(1.0 / re), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]]]]]

Derivation

1. Split input into 4 regimes

2. if re < 3.322933833770276e20

   1. Initial program 33.0

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

   2. Simplified 5.9

      \[\leadsto \color{blue}{0.5 \cdot \sqrt{2 \cdot \left(\mathsf{hypot}\left(re, im\right) - re\right)}} \]
      Proof
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (hypot.f64 re im) re)))): 0 points increase in error, 0 points decrease in error
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (Rewrite<= hypot-def_binary64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im)))) re)))): 141 points increase in error, 0 points decrease in error

   if 3.322933833770276e20 < re < 4.92686146157179325e88

   1. Initial program 47.0

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

   2. Simplified 31.0

      \[\leadsto \color{blue}{0.5 \cdot \sqrt{2 \cdot \left(\mathsf{hypot}\left(re, im\right) - re\right)}} \]
      Proof
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (hypot.f64 re im) re)))): 0 points increase in error, 0 points decrease in error
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (Rewrite<= hypot-def_binary64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im)))) re)))): 141 points increase in error, 0 points decrease in error

   3. Taylor expanded in im around 0 26.2

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

   4. Simplified 26.2

      \[\leadsto 0.5 \cdot \color{blue}{\left(im \cdot \left(\left(\sqrt{0.5} \cdot \sqrt{\frac{1}{re}}\right) \cdot \sqrt{2}\right)\right)} \]
      Proof
      (*.f64 im (*.f64 (*.f64 (sqrt.f64 1/2) (sqrt.f64 (/.f64 1 re))) (sqrt.f64 2))): 0 points increase in error, 0 points decrease in error
      (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 im (*.f64 (sqrt.f64 1/2) (sqrt.f64 (/.f64 1 re)))) (sqrt.f64 2))): 13 points increase in error, 21 points decrease in error
      (*.f64 (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 im (sqrt.f64 1/2)) (sqrt.f64 (/.f64 1 re)))) (sqrt.f64 2)): 11 points increase in error, 22 points decrease in error
      (*.f64 (*.f64 (Rewrite<= *-commutative_binary64 (*.f64 (sqrt.f64 1/2) im)) (sqrt.f64 (/.f64 1 re))) (sqrt.f64 2)): 0 points increase in error, 0 points decrease in error
      (Rewrite<= *-commutative_binary64 (*.f64 (sqrt.f64 2) (*.f64 (*.f64 (sqrt.f64 1/2) im) (sqrt.f64 (/.f64 1 re))))): 0 points increase in error, 0 points decrease in error
      (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 (sqrt.f64 2) (*.f64 (sqrt.f64 1/2) im)) (sqrt.f64 (/.f64 1 re)))): 22 points increase in error, 7 points decrease in error

   5. Applied egg-rr 26.0

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

   6. Applied egg-rr 26.7

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

   if 4.92686146157179325e88 < re < 3.6805890090252481e151

   1. Initial program 53.5

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

   2. Simplified 38.5

      \[\leadsto \color{blue}{0.5 \cdot \sqrt{2 \cdot \left(\mathsf{hypot}\left(re, im\right) - re\right)}} \]
      Proof
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (hypot.f64 re im) re)))): 0 points increase in error, 0 points decrease in error
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (Rewrite<= hypot-def_binary64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im)))) re)))): 141 points increase in error, 0 points decrease in error

   3. Applied egg-rr 44.1

      \[\leadsto 0.5 \cdot \sqrt{2 \cdot \color{blue}{\left(\left(\mathsf{hypot}\left(re, im\right) - re\right) + \mathsf{fma}\left(-\sqrt[3]{re}, {\left(\sqrt[3]{re}\right)}^{2}, re\right)\right)}} \]

   if 3.6805890090252481e151 < re

   1. Initial program 63.8

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

   2. Simplified 41.9

      \[\leadsto \color{blue}{0.5 \cdot \sqrt{2 \cdot \left(\mathsf{hypot}\left(re, im\right) - re\right)}} \]
      Proof
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (hypot.f64 re im) re)))): 0 points increase in error, 0 points decrease in error
      (*.f64 1/2 (sqrt.f64 (*.f64 2 (-.f64 (Rewrite<= hypot-def_binary64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im)))) re)))): 141 points increase in error, 0 points decrease in error

   3. Taylor expanded in im around 0 7.6

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

   4. Simplified 7.7

      \[\leadsto 0.5 \cdot \color{blue}{\left(im \cdot \left(\left(\sqrt{0.5} \cdot \sqrt{\frac{1}{re}}\right) \cdot \sqrt{2}\right)\right)} \]
      Proof
      (*.f64 im (*.f64 (*.f64 (sqrt.f64 1/2) (sqrt.f64 (/.f64 1 re))) (sqrt.f64 2))): 0 points increase in error, 0 points decrease in error
      (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 im (*.f64 (sqrt.f64 1/2) (sqrt.f64 (/.f64 1 re)))) (sqrt.f64 2))): 13 points increase in error, 21 points decrease in error
      (*.f64 (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 im (sqrt.f64 1/2)) (sqrt.f64 (/.f64 1 re)))) (sqrt.f64 2)): 11 points increase in error, 22 points decrease in error
      (*.f64 (*.f64 (Rewrite<= *-commutative_binary64 (*.f64 (sqrt.f64 1/2) im)) (sqrt.f64 (/.f64 1 re))) (sqrt.f64 2)): 0 points increase in error, 0 points decrease in error
      (Rewrite<= *-commutative_binary64 (*.f64 (sqrt.f64 2) (*.f64 (*.f64 (sqrt.f64 1/2) im) (sqrt.f64 (/.f64 1 re))))): 0 points increase in error, 0 points decrease in error
      (Rewrite<= associate-*l*_binary64 (*.f64 (*.f64 (sqrt.f64 2) (*.f64 (sqrt.f64 1/2) im)) (sqrt.f64 (/.f64 1 re)))): 22 points increase in error, 7 points decrease in error

   5. Applied egg-rr 7.2

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

4. Final simplification 9.2

   \[\leadsto \begin{array}{l} \mathbf{if}\;re \leq 3.322933833770276 \cdot 10^{+20}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(\mathsf{hypot}\left(re, im\right) - re\right)}\\ \mathbf{elif}\;re \leq 4.926861461571793 \cdot 10^{+88}:\\ \;\;\;\;0.5 \cdot \frac{1}{\frac{\sqrt{re}}{im}}\\ \mathbf{elif}\;re \leq 3.680589009025248 \cdot 10^{+151}:\\ \;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(\left(\mathsf{hypot}\left(re, im\right) - re\right) + \mathsf{fma}\left(-\sqrt[3]{re}, {\left(\sqrt[3]{re}\right)}^{2}, re\right)\right)}\\ \mathbf{else}:\\ \;\;\;\;0.5 \cdot \left(im \cdot \sqrt{\frac{1}{re}}\right)\\ \end{array} \]

Alternatives

Alternative 1
Error	24.4
Cost	7180

\[\begin{array}{l} t_0 := 0.5 \cdot \sqrt{2 \cdot im}\\ t_1 := 0.5 \cdot \left(im \cdot {re}^{-0.5}\right)\\ \mathbf{if}\;re \leq 7.568300857283601 \cdot 10^{+35}:\\ \;\;\;\;t_0\\ \mathbf{elif}\;re \leq 4.926861461571793 \cdot 10^{+88}:\\ \;\;\;\;t_1\\ \mathbf{elif}\;re \leq 3.680589009025248 \cdot 10^{+151}:\\ \;\;\;\;t_0\\ \mathbf{else}:\\ \;\;\;\;t_1\\ \end{array} \]

Reproduce

herbie shell --seed 2022291 
(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)))))