\sqrt{0.5 \cdot \left(1 + \frac{x}{\sqrt{\left(4 \cdot p\right) \cdot p + x \cdot x}}\right)}
\begin{array}{l}
t_0 := p \cdot \left(4 \cdot p\right)\\
\mathbf{if}\;\frac{x}{\sqrt{t_0 + x \cdot x}} \leq -1:\\
\;\;\;\;\sqrt{\frac{p}{\frac{x \cdot x}{p}}}\\
\mathbf{else}:\\
\;\;\;\;\sqrt{\mathsf{fma}\left(0.5, \frac{x}{\sqrt{\mathsf{fma}\left(x, x, t_0\right)}}, 0.5\right)}\\
\end{array}
(FPCore (p x) :precision binary64 (sqrt (* 0.5 (+ 1.0 (/ x (sqrt (+ (* (* 4.0 p) p) (* x x))))))))
(FPCore (p x)
:precision binary64
(let* ((t_0 (* p (* 4.0 p))))
(if (<= (/ x (sqrt (+ t_0 (* x x)))) -1.0)
(sqrt (/ p (/ (* x x) p)))
(sqrt (fma 0.5 (/ x (sqrt (fma x x t_0))) 0.5)))))double code(double p, double x) {
return sqrt(0.5 * (1.0 + (x / sqrt(((4.0 * p) * p) + (x * x)))));
}
double code(double p, double x) {
double t_0 = p * (4.0 * p);
double tmp;
if ((x / sqrt(t_0 + (x * x))) <= -1.0) {
tmp = sqrt(p / ((x * x) / p));
} else {
tmp = sqrt(fma(0.5, (x / sqrt(fma(x, x, t_0))), 0.5));
}
return tmp;
}




Bits error versus p




Bits error versus x
| Original | 13.1 |
|---|---|
| Target | 13.1 |
| Herbie | 5.6 |
if (/.f64 x (sqrt.f64 (+.f64 (*.f64 (*.f64 4 p) p) (*.f64 x x)))) < -1Initial program 53.3
Simplified53.3
Taylor expanded in x around -inf 29.2
Simplified22.3
if -1 < (/.f64 x (sqrt.f64 (+.f64 (*.f64 (*.f64 4 p) p) (*.f64 x x)))) Initial program 0.2
Simplified0.2
Final simplification5.6
herbie shell --seed 2021215
(FPCore (p x)
:name "Given's Rotation SVD example"
:precision binary64
:pre (< 1e-150 (fabs x) 1e+150)
:herbie-target
(sqrt (+ 0.5 (/ (copysign 0.5 x) (hypot 1.0 (/ (* 2.0 p) x)))))
(sqrt (* 0.5 (+ 1.0 (/ x (sqrt (+ (* (* 4.0 p) p) (* x x))))))))