## Givens Rotation in old SLATEC Routine

Open discussion regarding features, bugs, issues, vendors, etc.

### Givens Rotation in old SLATEC Routine

I am translating an old Fortran sub-routine (R1UPDT) from the SLATEC library into C++ and came across a section of code that puzzles me.

A Givens Rotation is pretty straightforward; it comes down to computing r = 1.0/sqrt(a^2 + b^2)
Usually, this equation is put in a slightly different form, depending upon whether abs(a) < abs( b) or abs(a) > abs(b). A temporary variable, c, is used such that c is smaller than 1.

So, if abs(a) < abs(b), c = a/b and the equation is written r = 1.0 / sqrt(c^2 + 1)

On the other hand, if abs(a) > abs(b), c = b/a and the equation is written r = 1.0 / 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.

However, the Fortan code in R1UPDT is written such that r = 0.5 / sqrt(0.25 + 0.25 * c^2)

i.e.,
Code: Select all
`SIN = 0.5 / SQRT(0.25 + 0.25 * COTAN**2)`

I can't help but ask, why leave the code like this? The equation could be factored by hand simply enough: take the 0.5 factor out of the radicand in the denominator, cancel the common factor from numerator and denominator, and avoid a multiplication process:

Code: Select all
`SIN = 1.0 / SQRT(1.0 + COTAN**2)`

Since this code was originally written in the 80's, I am wondering if keeping the 0.25 factor in the radicand in the denominator is one more technique for avoiding overflow (or something else.)

Any ideas or suggestions? This routine is actually a helper routine for SNSQE, which is a multi-dimensional root-finding routine. Anybody here use SNSQE?
DavidB

Posts: 38
Joined: Thu Dec 15, 2005 1:49 pm