
(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 5 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 (<= (sqrt (* (+ (sqrt (+ (* im im) (* re re))) re) 2.0)) 0.0) (* 0.5 (sqrt (/ (- im) (/ re im)))) (* (sqrt (* (+ (hypot im re) re) 2.0)) 0.5)))
double code(double re, double im) {
double tmp;
if (sqrt(((sqrt(((im * im) + (re * re))) + re) * 2.0)) <= 0.0) {
tmp = 0.5 * sqrt((-im / (re / im)));
} else {
tmp = sqrt(((hypot(im, re) + re) * 2.0)) * 0.5;
}
return tmp;
}
public static double code(double re, double im) {
double tmp;
if (Math.sqrt(((Math.sqrt(((im * im) + (re * re))) + re) * 2.0)) <= 0.0) {
tmp = 0.5 * Math.sqrt((-im / (re / im)));
} else {
tmp = Math.sqrt(((Math.hypot(im, re) + re) * 2.0)) * 0.5;
}
return tmp;
}
def code(re, im): tmp = 0 if math.sqrt(((math.sqrt(((im * im) + (re * re))) + re) * 2.0)) <= 0.0: tmp = 0.5 * math.sqrt((-im / (re / im))) else: tmp = math.sqrt(((math.hypot(im, re) + re) * 2.0)) * 0.5 return tmp
function code(re, im) tmp = 0.0 if (sqrt(Float64(Float64(sqrt(Float64(Float64(im * im) + Float64(re * re))) + re) * 2.0)) <= 0.0) tmp = Float64(0.5 * sqrt(Float64(Float64(-im) / Float64(re / im)))); else tmp = Float64(sqrt(Float64(Float64(hypot(im, re) + re) * 2.0)) * 0.5); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (sqrt(((sqrt(((im * im) + (re * re))) + re) * 2.0)) <= 0.0) tmp = 0.5 * sqrt((-im / (re / im))); else tmp = sqrt(((hypot(im, re) + re) * 2.0)) * 0.5; end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[N[Sqrt[N[(N[(N[Sqrt[N[(N[(im * im), $MachinePrecision] + N[(re * re), $MachinePrecision]), $MachinePrecision]], $MachinePrecision] + re), $MachinePrecision] * 2.0), $MachinePrecision]], $MachinePrecision], 0.0], N[(0.5 * N[Sqrt[N[((-im) / N[(re / im), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[(N[Sqrt[N[(N[(N[Sqrt[im ^ 2 + re ^ 2], $MachinePrecision] + re), $MachinePrecision] * 2.0), $MachinePrecision]], $MachinePrecision] * 0.5), $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;\sqrt{\left(\sqrt{im \cdot im + re \cdot re} + re\right) \cdot 2} \leq 0:\\
\;\;\;\;0.5 \cdot \sqrt{\frac{-im}{\frac{re}{im}}}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{\left(\mathsf{hypot}\left(im, re\right) + re\right) \cdot 2} \cdot 0.5\\
\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 10.5%
lift-*.f64N/A
*-commutativeN/A
lower-*.f6410.5
lift-*.f64N/A
*-commutativeN/A
lower-*.f6410.5
lift-sqrt.f64N/A
lift-+.f64N/A
+-commutativeN/A
lift-*.f64N/A
lift-*.f64N/A
lower-hypot.f6410.5
Applied rewrites10.5%
Taylor expanded in re around -inf
mul-1-negN/A
unpow2N/A
associate-/l*N/A
distribute-lft-neg-inN/A
mul-1-negN/A
lower-*.f64N/A
mul-1-negN/A
lower-neg.f64N/A
lower-/.f6462.9
Applied rewrites62.9%
Applied rewrites63.0%
if 0.0 < (sqrt.f64 (*.f64 #s(literal 2 binary64) (+.f64 (sqrt.f64 (+.f64 (*.f64 re re) (*.f64 im im))) re))) Initial program 46.4%
lift-*.f64N/A
*-commutativeN/A
lower-*.f6446.4
lift-*.f64N/A
*-commutativeN/A
lower-*.f6446.4
lift-sqrt.f64N/A
lift-+.f64N/A
+-commutativeN/A
lift-*.f64N/A
lift-*.f64N/A
lower-hypot.f6490.4
Applied rewrites90.4%
Final simplification86.5%
(FPCore (re im) :precision binary64 (if (<= re -4.5e-28) (* (sqrt (* (/ (- im) re) im)) 0.5) (if (<= re 3.1e+86) (* (sqrt (* (+ im re) 2.0)) 0.5) (sqrt re))))
double code(double re, double im) {
double tmp;
if (re <= -4.5e-28) {
tmp = sqrt(((-im / re) * im)) * 0.5;
} else if (re <= 3.1e+86) {
tmp = sqrt(((im + re) * 2.0)) * 0.5;
} 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 <= (-4.5d-28)) then
tmp = sqrt(((-im / re) * im)) * 0.5d0
else if (re <= 3.1d+86) then
tmp = sqrt(((im + re) * 2.0d0)) * 0.5d0
else
tmp = sqrt(re)
end if
code = tmp
end function
public static double code(double re, double im) {
double tmp;
if (re <= -4.5e-28) {
tmp = Math.sqrt(((-im / re) * im)) * 0.5;
} else if (re <= 3.1e+86) {
tmp = Math.sqrt(((im + re) * 2.0)) * 0.5;
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
def code(re, im): tmp = 0 if re <= -4.5e-28: tmp = math.sqrt(((-im / re) * im)) * 0.5 elif re <= 3.1e+86: tmp = math.sqrt(((im + re) * 2.0)) * 0.5 else: tmp = math.sqrt(re) return tmp
function code(re, im) tmp = 0.0 if (re <= -4.5e-28) tmp = Float64(sqrt(Float64(Float64(Float64(-im) / re) * im)) * 0.5); elseif (re <= 3.1e+86) tmp = Float64(sqrt(Float64(Float64(im + re) * 2.0)) * 0.5); else tmp = sqrt(re); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (re <= -4.5e-28) tmp = sqrt(((-im / re) * im)) * 0.5; elseif (re <= 3.1e+86) tmp = sqrt(((im + re) * 2.0)) * 0.5; else tmp = sqrt(re); end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[re, -4.5e-28], N[(N[Sqrt[N[(N[((-im) / re), $MachinePrecision] * im), $MachinePrecision]], $MachinePrecision] * 0.5), $MachinePrecision], If[LessEqual[re, 3.1e+86], N[(N[Sqrt[N[(N[(im + re), $MachinePrecision] * 2.0), $MachinePrecision]], $MachinePrecision] * 0.5), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;re \leq -4.5 \cdot 10^{-28}:\\
\;\;\;\;\sqrt{\frac{-im}{re} \cdot im} \cdot 0.5\\
\mathbf{elif}\;re \leq 3.1 \cdot 10^{+86}:\\
\;\;\;\;\sqrt{\left(im + re\right) \cdot 2} \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < -4.4999999999999998e-28Initial program 17.8%
Taylor expanded in re around -inf
mul-1-negN/A
unpow2N/A
associate-/l*N/A
distribute-lft-neg-inN/A
mul-1-negN/A
lower-*.f64N/A
mul-1-negN/A
lower-neg.f64N/A
lower-/.f6456.7
Applied rewrites56.7%
if -4.4999999999999998e-28 < re < 3.1000000000000002e86Initial program 59.5%
Taylor expanded in re around 0
+-commutativeN/A
lower-+.f6441.1
Applied rewrites41.1%
if 3.1000000000000002e86 < re Initial program 24.9%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lft-identityN/A
lower-sqrt.f6487.3
Applied rewrites87.3%
Final simplification54.5%
(FPCore (re im) :precision binary64 (if (<= re -3.2e+132) (* (sqrt (* (+ (- re) re) 2.0)) 0.5) (if (<= re 3.1e+86) (* (sqrt (* (+ im re) 2.0)) 0.5) (sqrt re))))
double code(double re, double im) {
double tmp;
if (re <= -3.2e+132) {
tmp = sqrt(((-re + re) * 2.0)) * 0.5;
} else if (re <= 3.1e+86) {
tmp = sqrt(((im + re) * 2.0)) * 0.5;
} 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 <= (-3.2d+132)) then
tmp = sqrt(((-re + re) * 2.0d0)) * 0.5d0
else if (re <= 3.1d+86) then
tmp = sqrt(((im + re) * 2.0d0)) * 0.5d0
else
tmp = sqrt(re)
end if
code = tmp
end function
public static double code(double re, double im) {
double tmp;
if (re <= -3.2e+132) {
tmp = Math.sqrt(((-re + re) * 2.0)) * 0.5;
} else if (re <= 3.1e+86) {
tmp = Math.sqrt(((im + re) * 2.0)) * 0.5;
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
def code(re, im): tmp = 0 if re <= -3.2e+132: tmp = math.sqrt(((-re + re) * 2.0)) * 0.5 elif re <= 3.1e+86: tmp = math.sqrt(((im + re) * 2.0)) * 0.5 else: tmp = math.sqrt(re) return tmp
function code(re, im) tmp = 0.0 if (re <= -3.2e+132) tmp = Float64(sqrt(Float64(Float64(Float64(-re) + re) * 2.0)) * 0.5); elseif (re <= 3.1e+86) tmp = Float64(sqrt(Float64(Float64(im + re) * 2.0)) * 0.5); else tmp = sqrt(re); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (re <= -3.2e+132) tmp = sqrt(((-re + re) * 2.0)) * 0.5; elseif (re <= 3.1e+86) tmp = sqrt(((im + re) * 2.0)) * 0.5; else tmp = sqrt(re); end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[re, -3.2e+132], N[(N[Sqrt[N[(N[((-re) + re), $MachinePrecision] * 2.0), $MachinePrecision]], $MachinePrecision] * 0.5), $MachinePrecision], If[LessEqual[re, 3.1e+86], N[(N[Sqrt[N[(N[(im + re), $MachinePrecision] * 2.0), $MachinePrecision]], $MachinePrecision] * 0.5), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;re \leq -3.2 \cdot 10^{+132}:\\
\;\;\;\;\sqrt{\left(\left(-re\right) + re\right) \cdot 2} \cdot 0.5\\
\mathbf{elif}\;re \leq 3.1 \cdot 10^{+86}:\\
\;\;\;\;\sqrt{\left(im + re\right) \cdot 2} \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < -3.1999999999999997e132Initial program 5.8%
Taylor expanded in re around -inf
mul-1-negN/A
lower-neg.f6427.3
Applied rewrites27.3%
if -3.1999999999999997e132 < re < 3.1000000000000002e86Initial program 53.1%
Taylor expanded in re around 0
+-commutativeN/A
lower-+.f6437.8
Applied rewrites37.8%
if 3.1000000000000002e86 < re Initial program 24.9%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lft-identityN/A
lower-sqrt.f6487.3
Applied rewrites87.3%
Final simplification46.3%
(FPCore (re im) :precision binary64 (if (<= re 1.75e-61) (* (sqrt (* im 2.0)) 0.5) (sqrt re)))
double code(double re, double im) {
double tmp;
if (re <= 1.75e-61) {
tmp = sqrt((im * 2.0)) * 0.5;
} 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.75d-61) then
tmp = sqrt((im * 2.0d0)) * 0.5d0
else
tmp = sqrt(re)
end if
code = tmp
end function
public static double code(double re, double im) {
double tmp;
if (re <= 1.75e-61) {
tmp = Math.sqrt((im * 2.0)) * 0.5;
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
def code(re, im): tmp = 0 if re <= 1.75e-61: tmp = math.sqrt((im * 2.0)) * 0.5 else: tmp = math.sqrt(re) return tmp
function code(re, im) tmp = 0.0 if (re <= 1.75e-61) tmp = Float64(sqrt(Float64(im * 2.0)) * 0.5); else tmp = sqrt(re); end return tmp end
function tmp_2 = code(re, im) tmp = 0.0; if (re <= 1.75e-61) tmp = sqrt((im * 2.0)) * 0.5; else tmp = sqrt(re); end tmp_2 = tmp; end
code[re_, im_] := If[LessEqual[re, 1.75e-61], N[(N[Sqrt[N[(im * 2.0), $MachinePrecision]], $MachinePrecision] * 0.5), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]
\begin{array}{l}
\\
\begin{array}{l}
\mathbf{if}\;re \leq 1.75 \cdot 10^{-61}:\\
\;\;\;\;\sqrt{im \cdot 2} \cdot 0.5\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < 1.7500000000000002e-61Initial program 42.4%
Taylor expanded in re around 0
lower-*.f6431.1
Applied rewrites31.1%
if 1.7500000000000002e-61 < re Initial program 39.0%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lft-identityN/A
lower-sqrt.f6475.3
Applied rewrites75.3%
Final simplification44.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 41.4%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lft-identityN/A
lower-sqrt.f6425.8
Applied rewrites25.8%
(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 2024255
(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)))))