
(FPCore (re im) :precision binary64 (* 0.5 (sqrt (* 2.0 (+ (sqrt (+ (* re re) (* im im))) re)))))
double code(double re, double im) {
return 0.5 * sqrt((2.0 * (sqrt(((re * re) + (im * im))) + re)));
}
real(8) function code(re, im)
real(8), intent (in) :: re
real(8), intent (in) :: im
code = 0.5d0 * sqrt((2.0d0 * (sqrt(((re * re) + (im * im))) + re)))
end function
public static double code(double re, double im) {
return 0.5 * Math.sqrt((2.0 * (Math.sqrt(((re * re) + (im * im))) + re)));
}
def code(re, im): return 0.5 * math.sqrt((2.0 * (math.sqrt(((re * re) + (im * im))) + re)))
function code(re, im) return Float64(0.5 * sqrt(Float64(2.0 * Float64(sqrt(Float64(Float64(re * re) + Float64(im * im))) + re)))) end
function tmp = code(re, im) tmp = 0.5 * sqrt((2.0 * (sqrt(((re * re) + (im * im))) + re))); end
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]
\begin{array}{l}
\\
0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} + re\right)}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 4 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (re im) :precision binary64 (* 0.5 (sqrt (* 2.0 (+ (sqrt (+ (* re re) (* im im))) re)))))
double code(double re, double im) {
return 0.5 * sqrt((2.0 * (sqrt(((re * re) + (im * im))) + re)));
}
real(8) function code(re, im)
real(8), intent (in) :: re
real(8), intent (in) :: im
code = 0.5d0 * sqrt((2.0d0 * (sqrt(((re * re) + (im * im))) + re)))
end function
public static double code(double re, double im) {
return 0.5 * Math.sqrt((2.0 * (Math.sqrt(((re * re) + (im * im))) + re)));
}
def code(re, im): return 0.5 * math.sqrt((2.0 * (math.sqrt(((re * re) + (im * im))) + re)))
function code(re, im) return Float64(0.5 * sqrt(Float64(2.0 * Float64(sqrt(Float64(Float64(re * re) + Float64(im * im))) + re)))) end
function tmp = code(re, im) tmp = 0.5 * sqrt((2.0 * (sqrt(((re * re) + (im * im))) + re))); end
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]
\begin{array}{l}
\\
0.5 \cdot \sqrt{2 \cdot \left(\sqrt{re \cdot re + im \cdot im} + re\right)}
\end{array}
(FPCore (re im) :precision binary64 (if (<= re -1.3e+53) (* 0.5 (sqrt (* im (- (/ im re))))) (sqrt (* 0.5 (+ re (hypot im re))))))
double code(double re, double im) {
double tmp;
if (re <= -1.3e+53) {
tmp = 0.5 * sqrt((im * -(im / re)));
} else {
tmp = sqrt((0.5 * (re + hypot(im, re))));
}
return tmp;
}
public static double code(double re, double im) {
double tmp;
if (re <= -1.3e+53) {
tmp = 0.5 * Math.sqrt((im * -(im / re)));
} else {
tmp = Math.sqrt((0.5 * (re + Math.hypot(im, re))));
}
return tmp;
}
def code(re, im): tmp = 0 if re <= -1.3e+53: tmp = 0.5 * math.sqrt((im * -(im / re))) else: tmp = math.sqrt((0.5 * (re + math.hypot(im, re)))) return tmp
function code(re, im) tmp = 0.0 if (re <= -1.3e+53) tmp = Float64(0.5 * sqrt(Float64(im * Float64(-Float64(im / re))))); else tmp = sqrt(Float64(0.5 * Float64(re + hypot(im, re)))); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (re <= -1.3e+53) tmp = 0.5 * sqrt((im * -(im / re))); else tmp = sqrt((0.5 * (re + hypot(im, re)))); end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[re, -1.3e+53], N[(0.5 * N[Sqrt[N[(im * (-N[(im / re), $MachinePrecision])), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[N[(0.5 * N[(re + N[Sqrt[im ^ 2 + re ^ 2], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;re \leq -1.3 \cdot 10^{+53}:\\
\;\;\;\;0.5 \cdot \sqrt{im \cdot \left(-\frac{im}{re}\right)}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{0.5 \cdot \left(re + \mathsf{hypot}\left(im, re\right)\right)}\\
\end{array}
\end{array}
if re < -1.29999999999999999e53Initial program 8.8%
sqr-neg8.8%
+-commutative8.8%
sqr-neg8.8%
+-commutative8.8%
distribute-rgt-in8.8%
cancel-sign-sub8.8%
distribute-rgt-out--8.8%
sub-neg8.8%
remove-double-neg8.8%
+-commutative8.8%
hypot-define30.1%
Simplified30.1%
Taylor expanded in re around -inf 53.8%
mul-1-neg53.8%
distribute-neg-frac253.8%
+-commutative53.8%
unpow253.8%
*-commutative53.8%
metadata-eval53.8%
pow-sqr53.8%
unpow253.8%
unpow253.8%
unpow253.8%
neg-sub053.8%
Simplified53.8%
Taylor expanded in im around 0 60.7%
mul-1-neg60.7%
unpow260.7%
associate-*l/66.5%
distribute-rgt-neg-in66.5%
neg-sub066.5%
Simplified66.5%
if -1.29999999999999999e53 < re Initial program 44.3%
sqr-neg44.3%
+-commutative44.3%
sqr-neg44.3%
+-commutative44.3%
distribute-rgt-in44.3%
cancel-sign-sub44.3%
distribute-rgt-out--44.3%
sub-neg44.3%
remove-double-neg44.3%
+-commutative44.3%
hypot-define91.0%
Simplified91.0%
*-commutative91.0%
hypot-define44.3%
+-commutative44.3%
*-commutative44.3%
add-sqr-sqrt44.1%
sqrt-unprod44.3%
*-commutative44.3%
*-commutative44.3%
swap-sqr44.3%
Applied egg-rr91.0%
*-commutative91.0%
associate-*r*91.0%
metadata-eval91.0%
hypot-define44.3%
unpow244.3%
+-commutative44.3%
unpow244.3%
hypot-undefine91.0%
Simplified91.0%
Final simplification86.0%
(FPCore (re im) :precision binary64 (if (<= re -1.9e+53) (* 0.5 (sqrt (* im (- (/ im re))))) (if (<= re 3.2e-15) (* 0.5 (sqrt (* 2.0 (+ re im)))) (sqrt re))))
double code(double re, double im) {
double tmp;
if (re <= -1.9e+53) {
tmp = 0.5 * sqrt((im * -(im / re)));
} else if (re <= 3.2e-15) {
tmp = 0.5 * sqrt((2.0 * (re + im)));
} else {
tmp = sqrt(re);
}
return tmp;
}
real(8) function code(re, im)
real(8), intent (in) :: re
real(8), intent (in) :: im
real(8) :: tmp
if (re <= (-1.9d+53)) then
tmp = 0.5d0 * sqrt((im * -(im / re)))
else if (re <= 3.2d-15) then
tmp = 0.5d0 * sqrt((2.0d0 * (re + im)))
else
tmp = sqrt(re)
end if
code = tmp
end function
public static double code(double re, double im) {
double tmp;
if (re <= -1.9e+53) {
tmp = 0.5 * Math.sqrt((im * -(im / re)));
} else if (re <= 3.2e-15) {
tmp = 0.5 * Math.sqrt((2.0 * (re + im)));
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
def code(re, im): tmp = 0 if re <= -1.9e+53: tmp = 0.5 * math.sqrt((im * -(im / re))) elif re <= 3.2e-15: tmp = 0.5 * math.sqrt((2.0 * (re + im))) else: tmp = math.sqrt(re) return tmp
function code(re, im) tmp = 0.0 if (re <= -1.9e+53) tmp = Float64(0.5 * sqrt(Float64(im * Float64(-Float64(im / re))))); elseif (re <= 3.2e-15) tmp = Float64(0.5 * sqrt(Float64(2.0 * Float64(re + im)))); else tmp = sqrt(re); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (re <= -1.9e+53) tmp = 0.5 * sqrt((im * -(im / re))); elseif (re <= 3.2e-15) tmp = 0.5 * sqrt((2.0 * (re + im))); else tmp = sqrt(re); end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[re, -1.9e+53], N[(0.5 * N[Sqrt[N[(im * (-N[(im / re), $MachinePrecision])), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], If[LessEqual[re, 3.2e-15], N[(0.5 * N[Sqrt[N[(2.0 * N[(re + im), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;re \leq -1.9 \cdot 10^{+53}:\\
\;\;\;\;0.5 \cdot \sqrt{im \cdot \left(-\frac{im}{re}\right)}\\
\mathbf{elif}\;re \leq 3.2 \cdot 10^{-15}:\\
\;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re + im\right)}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < -1.89999999999999999e53Initial program 8.8%
sqr-neg8.8%
+-commutative8.8%
sqr-neg8.8%
+-commutative8.8%
distribute-rgt-in8.8%
cancel-sign-sub8.8%
distribute-rgt-out--8.8%
sub-neg8.8%
remove-double-neg8.8%
+-commutative8.8%
hypot-define30.1%
Simplified30.1%
Taylor expanded in re around -inf 53.8%
mul-1-neg53.8%
distribute-neg-frac253.8%
+-commutative53.8%
unpow253.8%
*-commutative53.8%
metadata-eval53.8%
pow-sqr53.8%
unpow253.8%
unpow253.8%
unpow253.8%
neg-sub053.8%
Simplified53.8%
Taylor expanded in im around 0 60.7%
mul-1-neg60.7%
unpow260.7%
associate-*l/66.5%
distribute-rgt-neg-in66.5%
neg-sub066.5%
Simplified66.5%
if -1.89999999999999999e53 < re < 3.1999999999999999e-15Initial program 49.6%
sqr-neg49.6%
+-commutative49.6%
sqr-neg49.6%
+-commutative49.6%
distribute-rgt-in49.6%
cancel-sign-sub49.6%
distribute-rgt-out--49.6%
sub-neg49.6%
remove-double-neg49.6%
+-commutative49.6%
hypot-define86.6%
Simplified86.6%
Taylor expanded in re around 0 37.1%
if 3.1999999999999999e-15 < re Initial program 33.8%
sqr-neg33.8%
+-commutative33.8%
sqr-neg33.8%
+-commutative33.8%
distribute-rgt-in33.8%
cancel-sign-sub33.8%
distribute-rgt-out--33.8%
sub-neg33.8%
remove-double-neg33.8%
+-commutative33.8%
hypot-define100.0%
Simplified100.0%
Taylor expanded in re around inf 75.1%
unpow275.1%
rem-square-sqrt76.6%
*-commutative76.6%
associate-*r*76.6%
metadata-eval76.6%
Simplified76.6%
Final simplification53.5%
(FPCore (re im) :precision binary64 (if (<= re 2e-15) (* 0.5 (sqrt (* im 2.0))) (sqrt re)))
double code(double re, double im) {
double tmp;
if (re <= 2e-15) {
tmp = 0.5 * sqrt((im * 2.0));
} else {
tmp = sqrt(re);
}
return tmp;
}
real(8) function code(re, im)
real(8), intent (in) :: re
real(8), intent (in) :: im
real(8) :: tmp
if (re <= 2d-15) then
tmp = 0.5d0 * sqrt((im * 2.0d0))
else
tmp = sqrt(re)
end if
code = tmp
end function
public static double code(double re, double im) {
double tmp;
if (re <= 2e-15) {
tmp = 0.5 * Math.sqrt((im * 2.0));
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
def code(re, im): tmp = 0 if re <= 2e-15: tmp = 0.5 * math.sqrt((im * 2.0)) else: tmp = math.sqrt(re) return tmp
function code(re, im) tmp = 0.0 if (re <= 2e-15) tmp = Float64(0.5 * sqrt(Float64(im * 2.0))); else tmp = sqrt(re); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (re <= 2e-15) tmp = 0.5 * sqrt((im * 2.0)); else tmp = sqrt(re); end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[re, 2e-15], N[(0.5 * N[Sqrt[N[(im * 2.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;re \leq 2 \cdot 10^{-15}:\\
\;\;\;\;0.5 \cdot \sqrt{im \cdot 2}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < 2.0000000000000002e-15Initial program 38.3%
sqr-neg38.3%
+-commutative38.3%
sqr-neg38.3%
+-commutative38.3%
distribute-rgt-in38.3%
cancel-sign-sub38.3%
distribute-rgt-out--38.3%
sub-neg38.3%
remove-double-neg38.3%
+-commutative38.3%
hypot-define70.9%
Simplified70.9%
Taylor expanded in re around 0 27.3%
*-commutative27.3%
Simplified27.3%
if 2.0000000000000002e-15 < re Initial program 33.8%
sqr-neg33.8%
+-commutative33.8%
sqr-neg33.8%
+-commutative33.8%
distribute-rgt-in33.8%
cancel-sign-sub33.8%
distribute-rgt-out--33.8%
sub-neg33.8%
remove-double-neg33.8%
+-commutative33.8%
hypot-define100.0%
Simplified100.0%
Taylor expanded in re around inf 75.1%
unpow275.1%
rem-square-sqrt76.6%
*-commutative76.6%
associate-*r*76.6%
metadata-eval76.6%
Simplified76.6%
Final simplification40.4%
(FPCore (re im) :precision binary64 (sqrt re))
double code(double re, double im) {
return sqrt(re);
}
real(8) function code(re, im)
real(8), intent (in) :: re
real(8), intent (in) :: im
code = sqrt(re)
end function
public static double code(double re, double im) {
return Math.sqrt(re);
}
def code(re, im): return math.sqrt(re)
function code(re, im) return sqrt(re) end
function tmp = code(re, im) tmp = sqrt(re); end
code[re_, im_] := N[Sqrt[re], $MachinePrecision]
\begin{array}{l}
\\
\sqrt{re}
\end{array}
Initial program 37.1%
sqr-neg37.1%
+-commutative37.1%
sqr-neg37.1%
+-commutative37.1%
distribute-rgt-in37.1%
cancel-sign-sub37.1%
distribute-rgt-out--37.1%
sub-neg37.1%
remove-double-neg37.1%
+-commutative37.1%
hypot-define78.7%
Simplified78.7%
Taylor expanded in re around inf 27.6%
unpow227.6%
rem-square-sqrt28.1%
*-commutative28.1%
associate-*r*28.1%
metadata-eval28.1%
Simplified28.1%
Final simplification28.1%
(FPCore (re im)
:precision binary64
(let* ((t_0 (sqrt (+ (* re re) (* im im)))))
(if (< re 0.0)
(* 0.5 (* (sqrt 2.0) (sqrt (/ (* im im) (- t_0 re)))))
(* 0.5 (sqrt (* 2.0 (+ t_0 re)))))))
double code(double re, double im) {
double t_0 = sqrt(((re * re) + (im * im)));
double tmp;
if (re < 0.0) {
tmp = 0.5 * (sqrt(2.0) * sqrt(((im * im) / (t_0 - re))));
} else {
tmp = 0.5 * sqrt((2.0 * (t_0 + re)));
}
return tmp;
}
real(8) function code(re, im)
real(8), intent (in) :: re
real(8), intent (in) :: im
real(8) :: t_0
real(8) :: tmp
t_0 = sqrt(((re * re) + (im * im)))
if (re < 0.0d0) then
tmp = 0.5d0 * (sqrt(2.0d0) * sqrt(((im * im) / (t_0 - re))))
else
tmp = 0.5d0 * sqrt((2.0d0 * (t_0 + re)))
end if
code = tmp
end function
public static double code(double re, double im) {
double t_0 = Math.sqrt(((re * re) + (im * im)));
double tmp;
if (re < 0.0) {
tmp = 0.5 * (Math.sqrt(2.0) * Math.sqrt(((im * im) / (t_0 - re))));
} else {
tmp = 0.5 * Math.sqrt((2.0 * (t_0 + re)));
}
return tmp;
}
def code(re, im): t_0 = math.sqrt(((re * re) + (im * im))) tmp = 0 if re < 0.0: tmp = 0.5 * (math.sqrt(2.0) * math.sqrt(((im * im) / (t_0 - re)))) else: tmp = 0.5 * math.sqrt((2.0 * (t_0 + re))) return tmp
function code(re, im) t_0 = sqrt(Float64(Float64(re * re) + Float64(im * im))) tmp = 0.0 if (re < 0.0) tmp = Float64(0.5 * Float64(sqrt(2.0) * sqrt(Float64(Float64(im * im) / Float64(t_0 - re))))); else tmp = Float64(0.5 * sqrt(Float64(2.0 * Float64(t_0 + re)))); end return tmp end
function tmp_2 = code(re, im) t_0 = sqrt(((re * re) + (im * im))); tmp = 0.0; if (re < 0.0) tmp = 0.5 * (sqrt(2.0) * sqrt(((im * im) / (t_0 - re)))); else tmp = 0.5 * sqrt((2.0 * (t_0 + re))); end tmp_2 = tmp; end
code[re_, im_] := Block[{t$95$0 = N[Sqrt[N[(N[(re * re), $MachinePrecision] + N[(im * im), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]}, If[Less[re, 0.0], N[(0.5 * N[(N[Sqrt[2.0], $MachinePrecision] * N[Sqrt[N[(N[(im * im), $MachinePrecision] / N[(t$95$0 - re), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision], N[(0.5 * N[Sqrt[N[(2.0 * N[(t$95$0 + re), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
t_0 := \sqrt{re \cdot re + im \cdot im}\\
\mathbf{if}\;re < 0:\\
\;\;\;\;0.5 \cdot \left(\sqrt{2} \cdot \sqrt{\frac{im \cdot im}{t\_0 - re}}\right)\\
\mathbf{else}:\\
\;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(t\_0 + re\right)}\\
\end{array}
\end{array}
herbie shell --seed 2024107
(FPCore (re im)
:name "math.sqrt on complex, real part"
:precision binary64
:alt
(if (< re 0.0) (* 0.5 (* (sqrt 2.0) (sqrt (/ (* im im) (- (sqrt (+ (* re re) (* im im))) re))))) (* 0.5 (sqrt (* 2.0 (+ (sqrt (+ (* re re) (* im im))) re)))))
(* 0.5 (sqrt (* 2.0 (+ (sqrt (+ (* re re) (* im im))) re)))))