
(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}
im_m = (fabs.f64 im) (FPCore (re im_m) :precision binary64 (if (<= re -2.5e-65) (* (pow (- 0.0 re) -0.5) (* 0.5 im_m)) (* 0.5 (sqrt (* 2.0 (+ re (hypot re im_m)))))))
im_m = fabs(im);
double code(double re, double im_m) {
double tmp;
if (re <= -2.5e-65) {
tmp = pow((0.0 - re), -0.5) * (0.5 * im_m);
} 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 (re <= -2.5e-65) {
tmp = Math.pow((0.0 - re), -0.5) * (0.5 * im_m);
} 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 re <= -2.5e-65: tmp = math.pow((0.0 - re), -0.5) * (0.5 * im_m) 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 (re <= -2.5e-65) tmp = Float64((Float64(0.0 - re) ^ -0.5) * Float64(0.5 * im_m)); 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 (re <= -2.5e-65) tmp = ((0.0 - re) ^ -0.5) * (0.5 * im_m); 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[re, -2.5e-65], N[(N[Power[N[(0.0 - re), $MachinePrecision], -0.5], $MachinePrecision] * N[(0.5 * im$95$m), $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}\;re \leq -2.5 \cdot 10^{-65}:\\
\;\;\;\;{\left(0 - re\right)}^{-0.5} \cdot \left(0.5 \cdot im\_m\right)\\
\mathbf{else}:\\
\;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re + \mathsf{hypot}\left(re, im\_m\right)\right)}\\
\end{array}
\end{array}
if re < -2.49999999999999991e-65Initial program 13.3%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6434.9%
Simplified34.9%
Taylor expanded in re around -inf
mul-1-negN/A
neg-sub0N/A
--lowering--.f64N/A
/-lowering-/.f64N/A
unpow2N/A
*-lowering-*.f6441.3%
Simplified41.3%
sub0-negN/A
clear-numN/A
distribute-neg-fracN/A
metadata-evalN/A
/-lowering-/.f64N/A
/-lowering-/.f64N/A
*-lowering-*.f6441.4%
Applied egg-rr41.4%
*-commutativeN/A
associate-/r/N/A
sqrt-prodN/A
pow1/2N/A
sqrt-prodN/A
rem-square-sqrtN/A
associate-*l*N/A
*-commutativeN/A
rem-square-sqrtN/A
sqrt-prodN/A
pow1/2N/A
*-lowering-*.f64N/A
frac-2negN/A
metadata-evalN/A
inv-powN/A
pow-powN/A
pow-lowering-pow.f64N/A
neg-sub0N/A
--lowering--.f64N/A
metadata-evalN/A
pow1/2N/A
sqrt-prodN/A
rem-square-sqrtN/A
*-lowering-*.f6445.6%
Applied egg-rr45.6%
sub0-negN/A
neg-lowering-neg.f6445.6%
Applied egg-rr45.6%
if -2.49999999999999991e-65 < re Initial program 53.4%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6495.3%
Simplified95.3%
Final simplification81.3%
im_m = (fabs.f64 im) (FPCore (re im_m) :precision binary64 (if (<= re -4.5e-63) (* (pow (- 0.0 re) -0.5) (* 0.5 im_m)) (if (<= re 2.2e-25) (* 0.5 (sqrt (* 2.0 (+ re im_m)))) (sqrt re))))
im_m = fabs(im);
double code(double re, double im_m) {
double tmp;
if (re <= -4.5e-63) {
tmp = pow((0.0 - re), -0.5) * (0.5 * im_m);
} else if (re <= 2.2e-25) {
tmp = 0.5 * sqrt((2.0 * (re + im_m)));
} else {
tmp = sqrt(re);
}
return tmp;
}
im_m = abs(im)
real(8) function code(re, im_m)
real(8), intent (in) :: re
real(8), intent (in) :: im_m
real(8) :: tmp
if (re <= (-4.5d-63)) then
tmp = ((0.0d0 - re) ** (-0.5d0)) * (0.5d0 * im_m)
else if (re <= 2.2d-25) then
tmp = 0.5d0 * sqrt((2.0d0 * (re + im_m)))
else
tmp = sqrt(re)
end if
code = tmp
end function
im_m = Math.abs(im);
public static double code(double re, double im_m) {
double tmp;
if (re <= -4.5e-63) {
tmp = Math.pow((0.0 - re), -0.5) * (0.5 * im_m);
} else if (re <= 2.2e-25) {
tmp = 0.5 * Math.sqrt((2.0 * (re + im_m)));
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
im_m = math.fabs(im) def code(re, im_m): tmp = 0 if re <= -4.5e-63: tmp = math.pow((0.0 - re), -0.5) * (0.5 * im_m) elif re <= 2.2e-25: tmp = 0.5 * math.sqrt((2.0 * (re + im_m))) else: tmp = math.sqrt(re) return tmp
im_m = abs(im) function code(re, im_m) tmp = 0.0 if (re <= -4.5e-63) tmp = Float64((Float64(0.0 - re) ^ -0.5) * Float64(0.5 * im_m)); elseif (re <= 2.2e-25) tmp = Float64(0.5 * sqrt(Float64(2.0 * Float64(re + im_m)))); else tmp = sqrt(re); end return tmp end
im_m = abs(im); function tmp_2 = code(re, im_m) tmp = 0.0; if (re <= -4.5e-63) tmp = ((0.0 - re) ^ -0.5) * (0.5 * im_m); elseif (re <= 2.2e-25) tmp = 0.5 * sqrt((2.0 * (re + im_m))); else tmp = sqrt(re); end tmp_2 = tmp; end
im_m = N[Abs[im], $MachinePrecision] code[re_, im$95$m_] := If[LessEqual[re, -4.5e-63], N[(N[Power[N[(0.0 - re), $MachinePrecision], -0.5], $MachinePrecision] * N[(0.5 * im$95$m), $MachinePrecision]), $MachinePrecision], If[LessEqual[re, 2.2e-25], N[(0.5 * N[Sqrt[N[(2.0 * N[(re + im$95$m), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]]
\begin{array}{l}
im_m = \left|im\right|
\\
\begin{array}{l}
\mathbf{if}\;re \leq -4.5 \cdot 10^{-63}:\\
\;\;\;\;{\left(0 - re\right)}^{-0.5} \cdot \left(0.5 \cdot im\_m\right)\\
\mathbf{elif}\;re \leq 2.2 \cdot 10^{-25}:\\
\;\;\;\;0.5 \cdot \sqrt{2 \cdot \left(re + im\_m\right)}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < -4.5e-63Initial program 13.3%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6434.9%
Simplified34.9%
Taylor expanded in re around -inf
mul-1-negN/A
neg-sub0N/A
--lowering--.f64N/A
/-lowering-/.f64N/A
unpow2N/A
*-lowering-*.f6441.3%
Simplified41.3%
sub0-negN/A
clear-numN/A
distribute-neg-fracN/A
metadata-evalN/A
/-lowering-/.f64N/A
/-lowering-/.f64N/A
*-lowering-*.f6441.4%
Applied egg-rr41.4%
*-commutativeN/A
associate-/r/N/A
sqrt-prodN/A
pow1/2N/A
sqrt-prodN/A
rem-square-sqrtN/A
associate-*l*N/A
*-commutativeN/A
rem-square-sqrtN/A
sqrt-prodN/A
pow1/2N/A
*-lowering-*.f64N/A
frac-2negN/A
metadata-evalN/A
inv-powN/A
pow-powN/A
pow-lowering-pow.f64N/A
neg-sub0N/A
--lowering--.f64N/A
metadata-evalN/A
pow1/2N/A
sqrt-prodN/A
rem-square-sqrtN/A
*-lowering-*.f6445.6%
Applied egg-rr45.6%
sub0-negN/A
neg-lowering-neg.f6445.6%
Applied egg-rr45.6%
if -4.5e-63 < re < 2.2000000000000002e-25Initial program 59.7%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6492.4%
Simplified92.4%
Taylor expanded in re around 0
distribute-lft-outN/A
*-lowering-*.f64N/A
+-lowering-+.f6446.6%
Simplified46.6%
if 2.2000000000000002e-25 < re Initial program 43.4%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6499.9%
Simplified99.9%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f6477.1%
Simplified77.1%
*-lft-identityN/A
sqrt-lowering-sqrt.f6477.1%
Applied egg-rr77.1%
Final simplification54.8%
im_m = (fabs.f64 im) (FPCore (re im_m) :precision binary64 (if (<= re 5.5e-20) (* 0.5 (sqrt (* im_m 2.0))) (sqrt re)))
im_m = fabs(im);
double code(double re, double im_m) {
double tmp;
if (re <= 5.5e-20) {
tmp = 0.5 * sqrt((im_m * 2.0));
} else {
tmp = sqrt(re);
}
return tmp;
}
im_m = abs(im)
real(8) function code(re, im_m)
real(8), intent (in) :: re
real(8), intent (in) :: im_m
real(8) :: tmp
if (re <= 5.5d-20) then
tmp = 0.5d0 * sqrt((im_m * 2.0d0))
else
tmp = sqrt(re)
end if
code = tmp
end function
im_m = Math.abs(im);
public static double code(double re, double im_m) {
double tmp;
if (re <= 5.5e-20) {
tmp = 0.5 * Math.sqrt((im_m * 2.0));
} else {
tmp = Math.sqrt(re);
}
return tmp;
}
im_m = math.fabs(im) def code(re, im_m): tmp = 0 if re <= 5.5e-20: tmp = 0.5 * math.sqrt((im_m * 2.0)) else: tmp = math.sqrt(re) return tmp
im_m = abs(im) function code(re, im_m) tmp = 0.0 if (re <= 5.5e-20) tmp = Float64(0.5 * sqrt(Float64(im_m * 2.0))); else tmp = sqrt(re); end return tmp end
im_m = abs(im); function tmp_2 = code(re, im_m) tmp = 0.0; if (re <= 5.5e-20) tmp = 0.5 * sqrt((im_m * 2.0)); else tmp = sqrt(re); end tmp_2 = tmp; end
im_m = N[Abs[im], $MachinePrecision] code[re_, im$95$m_] := If[LessEqual[re, 5.5e-20], N[(0.5 * N[Sqrt[N[(im$95$m * 2.0), $MachinePrecision]], $MachinePrecision]), $MachinePrecision], N[Sqrt[re], $MachinePrecision]]
\begin{array}{l}
im_m = \left|im\right|
\\
\begin{array}{l}
\mathbf{if}\;re \leq 5.5 \cdot 10^{-20}:\\
\;\;\;\;0.5 \cdot \sqrt{im\_m \cdot 2}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{re}\\
\end{array}
\end{array}
if re < 5.4999999999999996e-20Initial program 41.7%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6470.0%
Simplified70.0%
Taylor expanded in re around 0
*-commutativeN/A
*-lowering-*.f6431.3%
Simplified31.3%
if 5.4999999999999996e-20 < re Initial program 43.4%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6499.9%
Simplified99.9%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f6477.1%
Simplified77.1%
*-lft-identityN/A
sqrt-lowering-sqrt.f6477.1%
Applied egg-rr77.1%
im_m = (fabs.f64 im) (FPCore (re im_m) :precision binary64 (sqrt re))
im_m = fabs(im);
double code(double re, double im_m) {
return sqrt(re);
}
im_m = abs(im)
real(8) function code(re, im_m)
real(8), intent (in) :: re
real(8), intent (in) :: im_m
code = sqrt(re)
end function
im_m = Math.abs(im);
public static double code(double re, double im_m) {
return Math.sqrt(re);
}
im_m = math.fabs(im) def code(re, im_m): return math.sqrt(re)
im_m = abs(im) function code(re, im_m) return sqrt(re) end
im_m = abs(im); function tmp = code(re, im_m) tmp = sqrt(re); end
im_m = N[Abs[im], $MachinePrecision] code[re_, im$95$m_] := N[Sqrt[re], $MachinePrecision]
\begin{array}{l}
im_m = \left|im\right|
\\
\sqrt{re}
\end{array}
Initial program 42.1%
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f64N/A
*-lowering-*.f64N/A
+-commutativeN/A
+-lowering-+.f64N/A
hypot-defineN/A
hypot-lowering-hypot.f6478.3%
Simplified78.3%
Taylor expanded in re around inf
*-commutativeN/A
unpow2N/A
rem-square-sqrtN/A
associate-*r*N/A
metadata-evalN/A
*-lowering-*.f64N/A
sqrt-lowering-sqrt.f6428.0%
Simplified28.0%
*-lft-identityN/A
sqrt-lowering-sqrt.f6428.0%
Applied egg-rr28.0%
(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 2024164
(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)))))