
(FPCore (a b) :precision binary64 (sqrt (- (* a a) (* b b))))
double code(double a, double b) {
return sqrt(((a * a) - (b * b)));
}
real(8) function code(a, b)
real(8), intent (in) :: a
real(8), intent (in) :: b
code = sqrt(((a * a) - (b * b)))
end function
public static double code(double a, double b) {
return Math.sqrt(((a * a) - (b * b)));
}
def code(a, b): return math.sqrt(((a * a) - (b * b)))
function code(a, b) return sqrt(Float64(Float64(a * a) - Float64(b * b))) end
function tmp = code(a, b) tmp = sqrt(((a * a) - (b * b))); end
code[a_, b_] := N[Sqrt[N[(N[(a * a), $MachinePrecision] - N[(b * b), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}
\\
\sqrt{a \cdot a - b \cdot b}
\end{array}
Sampling outcomes in binary64 precision:
Herbie found 3 alternatives:
| Alternative | Accuracy | Speedup |
|---|
(FPCore (a b) :precision binary64 (sqrt (- (* a a) (* b b))))
double code(double a, double b) {
return sqrt(((a * a) - (b * b)));
}
real(8) function code(a, b)
real(8), intent (in) :: a
real(8), intent (in) :: b
code = sqrt(((a * a) - (b * b)))
end function
public static double code(double a, double b) {
return Math.sqrt(((a * a) - (b * b)));
}
def code(a, b): return math.sqrt(((a * a) - (b * b)))
function code(a, b) return sqrt(Float64(Float64(a * a) - Float64(b * b))) end
function tmp = code(a, b) tmp = sqrt(((a * a) - (b * b))); end
code[a_, b_] := N[Sqrt[N[(N[(a * a), $MachinePrecision] - N[(b * b), $MachinePrecision]), $MachinePrecision]], $MachinePrecision]
\begin{array}{l}
\\
\sqrt{a \cdot a - b \cdot b}
\end{array}
a_m = (fabs.f64 a) (FPCore (a_m b) :precision binary64 (/ (sqrt (+ b a_m)) (pow (- a_m b) -0.5)))
a_m = fabs(a);
double code(double a_m, double b) {
return sqrt((b + a_m)) / pow((a_m - b), -0.5);
}
a_m = abs(a)
real(8) function code(a_m, b)
real(8), intent (in) :: a_m
real(8), intent (in) :: b
code = sqrt((b + a_m)) / ((a_m - b) ** (-0.5d0))
end function
a_m = Math.abs(a);
public static double code(double a_m, double b) {
return Math.sqrt((b + a_m)) / Math.pow((a_m - b), -0.5);
}
a_m = math.fabs(a) def code(a_m, b): return math.sqrt((b + a_m)) / math.pow((a_m - b), -0.5)
a_m = abs(a) function code(a_m, b) return Float64(sqrt(Float64(b + a_m)) / (Float64(a_m - b) ^ -0.5)) end
a_m = abs(a); function tmp = code(a_m, b) tmp = sqrt((b + a_m)) / ((a_m - b) ^ -0.5); end
a_m = N[Abs[a], $MachinePrecision] code[a$95$m_, b_] := N[(N[Sqrt[N[(b + a$95$m), $MachinePrecision]], $MachinePrecision] / N[Power[N[(a$95$m - b), $MachinePrecision], -0.5], $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
a_m = \left|a\right|
\\
\frac{\sqrt{b + a\_m}}{{\left(a\_m - b\right)}^{-0.5}}
\end{array}
Initial program 51.7%
lift-sqrt.f64N/A
pow1/2N/A
lift--.f64N/A
lift-*.f64N/A
lift-*.f64N/A
difference-of-squaresN/A
*-commutativeN/A
unpow-prod-downN/A
lower-*.f64N/A
pow1/2N/A
lower-sqrt.f64N/A
lower--.f64N/A
pow1/2N/A
lower-sqrt.f64N/A
+-commutativeN/A
lower-+.f6452.3
Applied rewrites52.3%
lift-sqrt.f64N/A
lift-+.f64N/A
flip3-+N/A
clear-numN/A
sqrt-divN/A
metadata-evalN/A
clear-numN/A
flip3-+N/A
lift-+.f64N/A
lower-/.f64N/A
lower-sqrt.f64N/A
inv-powN/A
lower-pow.f6452.4
Applied rewrites52.4%
lift-*.f64N/A
lift-/.f64N/A
un-div-invN/A
clear-numN/A
lower-/.f64N/A
lower-/.f6452.4
lift-sqrt.f64N/A
lift-pow.f64N/A
sqrt-pow1N/A
lower-pow.f64N/A
metadata-eval52.4
Applied rewrites52.4%
lift-/.f64N/A
lift-/.f64N/A
div-invN/A
associate-/r*N/A
lower-/.f64N/A
lift-pow.f64N/A
pow-flipN/A
metadata-evalN/A
pow1/2N/A
lower-sqrt.f64N/A
lift-+.f64N/A
+-commutativeN/A
lower-+.f64N/A
lift-sqrt.f64N/A
pow1/2N/A
pow-flipN/A
metadata-evalN/A
lower-pow.f6452.5
Applied rewrites52.5%
Final simplification52.5%
a_m = (fabs.f64 a) (FPCore (a_m b) :precision binary64 (fma (* (/ -0.5 a_m) b) b a_m))
a_m = fabs(a);
double code(double a_m, double b) {
return fma(((-0.5 / a_m) * b), b, a_m);
}
a_m = abs(a) function code(a_m, b) return fma(Float64(Float64(-0.5 / a_m) * b), b, a_m) end
a_m = N[Abs[a], $MachinePrecision] code[a$95$m_, b_] := N[(N[(N[(-0.5 / a$95$m), $MachinePrecision] * b), $MachinePrecision] * b + a$95$m), $MachinePrecision]
\begin{array}{l}
a_m = \left|a\right|
\\
\mathsf{fma}\left(\frac{-0.5}{a\_m} \cdot b, b, a\_m\right)
\end{array}
Initial program 51.7%
Taylor expanded in b around 0
+-commutativeN/A
associate-*r/N/A
associate-*l/N/A
metadata-evalN/A
distribute-neg-fracN/A
metadata-evalN/A
associate-*r/N/A
unpow2N/A
associate-*r*N/A
lower-fma.f64N/A
lower-*.f64N/A
associate-*r/N/A
metadata-evalN/A
distribute-neg-fracN/A
metadata-evalN/A
lower-/.f6453.0
Applied rewrites53.0%
a_m = (fabs.f64 a) (FPCore (a_m b) :precision binary64 a_m)
a_m = fabs(a);
double code(double a_m, double b) {
return a_m;
}
a_m = abs(a)
real(8) function code(a_m, b)
real(8), intent (in) :: a_m
real(8), intent (in) :: b
code = a_m
end function
a_m = Math.abs(a);
public static double code(double a_m, double b) {
return a_m;
}
a_m = math.fabs(a) def code(a_m, b): return a_m
a_m = abs(a) function code(a_m, b) return a_m end
a_m = abs(a); function tmp = code(a_m, b) tmp = a_m; end
a_m = N[Abs[a], $MachinePrecision] code[a$95$m_, b_] := a$95$m
\begin{array}{l}
a_m = \left|a\right|
\\
a\_m
\end{array}
Initial program 51.7%
lift--.f64N/A
lift-*.f64N/A
lift-*.f64N/A
difference-of-squaresN/A
*-commutativeN/A
lower-*.f64N/A
lower--.f64N/A
+-commutativeN/A
lower-+.f6452.4
Applied rewrites52.4%
Taylor expanded in a around inf
+-commutativeN/A
associate-*r/N/A
distribute-rgt1-inN/A
metadata-evalN/A
mul0-lftN/A
metadata-evalN/A
mul0-lftN/A
associate-*r/N/A
mul0-lftN/A
metadata-evalN/A
*-rgt-identity52.8
Applied rewrites52.8%
(FPCore (a b) :precision binary64 (* (sqrt (+ (fabs a) (fabs b))) (sqrt (- (fabs a) (fabs b)))))
double code(double a, double b) {
return sqrt((fabs(a) + fabs(b))) * sqrt((fabs(a) - fabs(b)));
}
real(8) function code(a, b)
real(8), intent (in) :: a
real(8), intent (in) :: b
code = sqrt((abs(a) + abs(b))) * sqrt((abs(a) - abs(b)))
end function
public static double code(double a, double b) {
return Math.sqrt((Math.abs(a) + Math.abs(b))) * Math.sqrt((Math.abs(a) - Math.abs(b)));
}
def code(a, b): return math.sqrt((math.fabs(a) + math.fabs(b))) * math.sqrt((math.fabs(a) - math.fabs(b)))
function code(a, b) return Float64(sqrt(Float64(abs(a) + abs(b))) * sqrt(Float64(abs(a) - abs(b)))) end
function tmp = code(a, b) tmp = sqrt((abs(a) + abs(b))) * sqrt((abs(a) - abs(b))); end
code[a_, b_] := N[(N[Sqrt[N[(N[Abs[a], $MachinePrecision] + N[Abs[b], $MachinePrecision]), $MachinePrecision]], $MachinePrecision] * N[Sqrt[N[(N[Abs[a], $MachinePrecision] - N[Abs[b], $MachinePrecision]), $MachinePrecision]], $MachinePrecision]), $MachinePrecision]
\begin{array}{l}
\\
\sqrt{\left|a\right| + \left|b\right|} \cdot \sqrt{\left|a\right| - \left|b\right|}
\end{array}
herbie shell --seed 2024331
(FPCore (a b)
:name "bug366, discussion (missed optimization)"
:precision binary64
:alt
(! :herbie-platform default (let* ((fa (fabs a)) (fb (fabs b))) (* (sqrt (+ fa fb)) (sqrt (- fa fb)))))
(sqrt (- (* a a) (* b b))))