
(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 13 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}
im_m = (fabs.f64 im) (FPCore (re im_m) :precision binary64 (if (<= (sqrt (* 2.0 (+ re (sqrt (+ (* re re) (* im_m im_m)))))) 0.0) (* 0.5 (exp (* (+ (log (/ im_m (- re))) (log im_m)) 0.5))) (* 0.5 (sqrt (* 2.0 (+ re (hypot re im_m)))))))
im_m = fabs(im);
double code(double re, double im_m) {
double tmp;
if (sqrt((2.0 * (re + sqrt(((re * re) + (im_m * im_m)))))) <= 0.0) {
tmp = 0.5 * exp(((log((im_m / -re)) + log(im_m)) * 0.5));
} else {
tmp = 0.5 * sqrt((2.0 * (re + hypot(re, im_m))));
}
return tmp;
}
im_m = Math.abs(im);
public static double code(double re, double im_m) {
double tmp;
if (Math.sqrt((2.0 * (re + Math.sqrt(((re * re) + (im_m * im_m)))))) <= 0.0) {
tmp = 0.5 * Math.exp(((Math.log((im_m / -re)) + Math.log(im_m)) * 0.5));
} else {
tmp = 0.5 * Math.sqrt((2.0 * (re + Math.hypot(re, im_m))));
}
return tmp;
}
im_m = math.fabs(im) def code(re, im_m): tmp = 0 if math.sqrt((2.0 * (re + math.sqrt(((re * re) + (im_m * im_m)))))) <= 0.0: tmp = 0.5 * math.exp(((math.log((im_m / -re)) + math.log(im_m)) * 0.5)) else: tmp = 0.5 * math.sqrt((2.0 * (re + math.hypot(re, im_m)))) return tmp
im_m = abs(im) function code(re, im_m) tmp = 0.0 if (sqrt(Float64(2.0 * Float64(re + sqrt(Float64(Float64(re * re) + Float64(im_m * im_m)))))) <= 0.0) tmp = Float64(0.5 * exp(Float64(Float64(log(Float64(im_m / Float64(-re))) + log(im_m)) * 0.5))); else tmp = Float64(0.5 * sqrt(Float64(2.0 * Float64(re + hypot(re, im_m))))); end return tmp end
im_m = abs(im); function tmp_2 = code(re, im_m) tmp = 0.0; if (sqrt((2.0 * (re + sqrt(((re * re) + (im_m * im_m)))))) <= 0.0) tmp = 0.5 * exp(((log((im_m / -re)) + log(im_m)) * 0.5)); else tmp = 0.5 * sqrt((2.0 * (re + hypot(re, im_m)))); end tmp_2 = tmp; end
im_m = N[Abs[im], $MachinePrecision] code[re_, im$95$m_] := If[LessEqual[N[Sqrt[N[(2.0 * N[(re + N[Sqrt[N[(N[(re * re), $MachinePrecision] + N[(im$95$m * im$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], 0.0], N[(0.5 * N[Exp[N[(N[(N[Log[N[(im$95$m / (-re)), $MachinePrecision]], $MachinePrecision] + N[Log[im$95$m], $MachinePrecision]), $MachinePrecision] * 0.5), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(0.5 * N[Sqrt[N[(2.0 * N[(re + N[Sqrt[re ^ 2 + im$95$m ^ 2], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
im_m = \left|im\right|
\\
\begin{array}{l}
\mathbf{if}\;\sqrt{2 \cdot \left(re + \sqrt{re \cdot re + im\_m \cdot im\_m}\right)} \leq 0:\\
\;\;\;\;0.5 \cdot e^{\left(\log \left(\frac{im\_m}{-re}\right) + \log im\_m\right) \cdot 0.5}\\
\mathbf{else}:\\
\;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re + \mathsf{hypot}\left(re, im\_m\right)\right)}\\
\end{array}
\end{array}
if (sqrt.f64 (*.f64 #s(literal 2 binary64) (+.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re))) < 0.0Initial program 8.3%
lift-*.f64N/A
*-commutativeN/A
lower-*.f648.3
lift-+.f64N/A
+-commutativeN/A
lower-+.f648.3
lift-+.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f648.3
Applied rewrites8.3%
lift-sqrt.f64N/A
pow1/2N/A
lift-+.f64N/A
+-commutativeN/A
lift-sqrt.f64N/A
lift-fma.f64N/A
+-commutativeN/A
lift-*.f64N/A
lift-hypot.f64N/A
lift-+.f64N/A
pow-to-expN/A
lower-exp.f64N/A
lower-*.f64N/A
Applied rewrites8.3%
Taylor expanded in re around -inf
lower-+.f64N/A
lower-log.f64N/A
lower-/.f64N/A
lower-log.f64N/A
unpow2N/A
lower-*.f6443.6
Applied rewrites43.6%
Applied rewrites90.7%
if 0.0 < (sqrt.f64 (*.f64 #s(literal 2 binary64) (+.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re))) Initial program 43.9%
lift-sqrt.f64N/A
lift-+.f64N/A
lift-*.f64N/A
lift-*.f64N/A
lower-hypot.f6488.7
Applied rewrites88.7%
Final simplification88.9%
im_m = (fabs.f64 im) (FPCore (re im_m) :precision binary64 (if (<= (sqrt (* 2.0 (+ re (sqrt (+ (* re re) (* im_m im_m)))))) 0.0) (* 0.5 (exp (* 0.5 (fma (log im_m) 2.0 (- (log (- re))))))) (* 0.5 (sqrt (* 2.0 (+ re (hypot re im_m)))))))
im_m = fabs(im);
double code(double re, double im_m) {
double tmp;
if (sqrt((2.0 * (re + sqrt(((re * re) + (im_m * im_m)))))) <= 0.0) {
tmp = 0.5 * exp((0.5 * fma(log(im_m), 2.0, -log(-re))));
} else {
tmp = 0.5 * sqrt((2.0 * (re + hypot(re, im_m))));
}
return tmp;
}
im_m = abs(im) function code(re, im_m) tmp = 0.0 if (sqrt(Float64(2.0 * Float64(re + sqrt(Float64(Float64(re * re) + Float64(im_m * im_m)))))) <= 0.0) tmp = Float64(0.5 * exp(Float64(0.5 * fma(log(im_m), 2.0, Float64(-log(Float64(-re))))))); else tmp = Float64(0.5 * sqrt(Float64(2.0 * Float64(re + hypot(re, im_m))))); end return tmp end
im_m = N[Abs[im], $MachinePrecision] code[re_, im$95$m_] := If[LessEqual[N[Sqrt[N[(2.0 * N[(re + N[Sqrt[N[(N[(re * re), $MachinePrecision] + N[(im$95$m * im$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision], 0.0], N[(0.5 * N[Exp[N[(0.5 * N[(N[Log[im$95$m], $MachinePrecision] * 2.0 + (-N[Log[(-re)], $MachinePrecision])), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(0.5 * N[Sqrt[N[(2.0 * N[(re + N[Sqrt[re ^ 2 + im$95$m ^ 2], $MachinePrecision]), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]]
\begin{array}{l}
im_m = \left|im\right|
\\
\begin{array}{l}
\mathbf{if}\;\sqrt{2 \cdot \left(re + \sqrt{re \cdot re + im\_m \cdot im\_m}\right)} \leq 0:\\
\;\;\;\;0.5 \cdot e^{0.5 \cdot \mathsf{fma}\left(\log im\_m, 2, -\log \left(-re\right)\right)}\\
\mathbf{else}:\\
\;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re + \mathsf{hypot}\left(re, im\_m\right)\right)}\\
\end{array}
\end{array}
if (sqrt.f64 (*.f64 #s(literal 2 binary64) (+.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re))) < 0.0Initial program 9.6%
lift-*.f64N/A
*-commutativeN/A
lower-*.f649.6
lift-+.f64N/A
+-commutativeN/A
lower-+.f649.6
lift-+.f64N/A
+-commutativeN/A
lift-*.f64N/A
lower-fma.f649.6
Applied rewrites9.6%
lift-sqrt.f64N/A
pow1/2N/A
lift-+.f64N/A
+-commutativeN/A
lift-sqrt.f64N/A
lift-fma.f64N/A
+-commutativeN/A
lift-*.f64N/A
lift-hypot.f64N/A
lift-+.f64N/A
pow-to-expN/A
lower-exp.f64N/A
lower-*.f64N/A
Applied rewrites9.6%
Taylor expanded in re around -inf
lower-+.f64N/A
lower-log.f64N/A
lower-/.f64N/A
lower-log.f64N/A
unpow2N/A
lower-*.f6452.3
Applied rewrites52.3%
Applied rewrites90.7%
if 0.0 < (sqrt.f64 (*.f64 #s(literal 2 binary64) (+.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re))) Initial program 45.7%
lift-sqrt.f64N/A
lift-+.f64N/A
lift-*.f64N/A
lift-*.f64N/A
lower-hypot.f6489.0
Applied rewrites89.0%
Final simplification89.2%
(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 2024228
(FPCore (re im)
:name "math.sqrt on complex, real part"
:precision binary64
:alt
(! :herbie-platform default (if (< re 0) (* 1/2 (* (sqrt 2) (sqrt (/ (* im im) (- (modulus re im) re))))) (* 1/2 (sqrt (* 2 (+ (modulus re im) re))))))
(* 0.5 (sqrt (* 2.0 (+ (sqrt (+ (* re re) (* im im))) re)))))