# SLATEC Routine Computes Givens Rotation in Unexpected Way

Some Background

I am working on a C++ translation of a SLATEC routine, R1UPDT, which performs a Givens rotation: $$r = \frac{1.0}{\sqrt{a^2 + b^2}}$$

Usually, this equation is put in a slightly different form, depending upon whether $$|a| < |b|$$. A temporary variable, $$c$$, is used such that $$c$$ is smaller than $$1$$.

$$|a| <|b| \to c =\frac{a}{b}$$ and the equation is written $$r = \frac{b}{\sqrt{c^2 + 1}}$$.

$$|a| >|b| \to c =\frac{b}{a}$$ and the equation is written $$r = \frac{a}{\sqrt{1+c^2}}$$.

Either way, the potential risk of overflow is avoided because the term squared is less than 1.

So far, so good.

On to the Question

However, the Fortran code in R1UPDT is written such that $$r = \frac{0.5}{\sqrt{0.25 + 0.25c^2}}$$

i.e.,

r = 0.5 / SQRT(0.25 + 0.25 * C**2)


I have to ask: why was the computation coded like this? The RHS of the equation could be factored by hand simply enough and coded as such.

i.e.,

r = 1.0 / SQRT(1.0 + C**2)


My first thought is that one multiplication could be eliminated and program execution made quicker. However, perhaps this is one more technique for avoiding overflow. Or perhaps ensuring that the radicand is in the range $$[0.25 - 0.5]$$ makes it quicker to calculate or yields a more accurate result. I am just speculating; perhaps somebody here can clarify.

Anybody here knows why the common factor was not taken out of this code?

• Have you looked at the blas routine drotg? See netlib.org/lapack/explore-3.1.1-html/drotg.f.html. Ive found that the lapack/blas code is more optimized and modernized than older libraries like SLATEC. Plus because givens rotations are provided in BLAS, you probably don't need to code it yourself – vibe Feb 5 '19 at 4:17
• You would have to look at the original SQRT function. It probably was not using the FPU, but some other approximation, perhaps by some approximating polynomial, that is constructed to be correct on the interval [0.25,1]. – Lutz Lehmann Feb 6 '19 at 10:41
• My best guess (thus not an answer) is that it was coded this way to mitigate the effects of wobbling precision with IBM System /360 hexadecimal floating-point format. Because this is a base-16 rather than a binary format, the mantissa can have up to three leading zero bits, and by ad-hoc rescaling once can try to avoid those cases. My second guess is that this may have something to with idiosyncrasies of pre-IEEE-754 floating-point arithmetic on systems like Cray. On modern systems, for best accuracy and intermediate overflow protection, use hypot(); or rhypot() if available. – njuffa Feb 6 '19 at 20:52

The square root is a relatively tough mathematical operation. Thus, it is reasonable to assume that in pre-IEEE-754 dominated world, the argument range could significantly influence the accuracy/performance of SQRT. Especially, as @njuffa mentioned in the comment with

...wobbling precision of IBM System/360 hexadecimal floating point format

However, I can not find a confirmation to that in this case of a relatively benign rescaling.

A lot of research has been done on accurate, fast, and efficient algorithms for square-root

In the past:

Or even in the recent future,

However, I do not see why the range $$[0.25,0.5]$$ would be preferable to $$[1,2]$$. Especially, from the point of guideline on using SQRT, even in pre-IEEE-754 era.

For example, see the manual IBM System/360 FORTRAN IV Library Subprograms, 1966. There, Table 12 "Performance Statistics" on page 39 gives accuracy/speed performance for various math functions in different IBM System/360 models.

While (in this table) certain math functions have varying accuracy/speed performance, DSQRT is listed as a single entry for the full range of argument values. Concluding, there are possibly no recommendations to limit the range of SQRT for IBM System/360 (btw, even for complex CSQRT). While there are infinitely many other platforms, I would conclude that there is not too much to read into that peculiar detail of SLATEC code.

The motivation for arranging the computation in this fashion was almost certainly the issue of wobbling precision when computing with the hexadecimal floating-point format of the IBM System/360. There, the mantissa is normalized such that $$num = \pm 16^{e} \cdot f$$, with $$\frac{1}{16} \le f \lt 1$$. This means that up to three leading bits of the mantissa can be zero.
In R1UPTD, $$c = \max (|v_{n}|, |v_{j}|) / \min (|v_{n}|, |v_{j}|) \le 1$$, which means that $$1 \le 1 + c^{2} \le 2$$, and therefore $$1 \le \sqrt{1 + c^{2}} \le \sqrt{2} \lt 2$$. Thus $$\sqrt{1 +c^{2}} = 16^{1} \cdot f$$, where $$\frac{1}{16} \le f \lt \frac{1}{8}$$. This means the leading three bits of the square root result will be zero, resulting in a computation with reduced accuracy.
Scaling with $$0.25$$ halves the result of the square root, so $$\sqrt{0.25 + 0.25 \cdot c^{2}} = 16^{0} \cdot f$$, where $$\frac{1}{2} \le f \lt 1$$, meaning there are no leading zero bits in the mantissa and full accuracy is preserved.
On modern platforms with binary floating-point arithmetic based on the IEEE-754 standard there is no need for this scaling, and most programming environments include a standard function hypot(x,y) that computes $$\sqrt{x^{2} + y^{2}}$$ with optimal accuracy and eliminating issues with overflow or underflow in intermediate computation. Some environments also offer a function rhypot(x,y) which computes $$\frac{1}{\sqrt{x^{2}+y^{2}}}$$ and which perfectly fits the needs of Givens rotation.