# Numerically stable way to compute sqrt((b²*c²) / (1-c²)) for c in [-1, 1]

For some real value `b` and `c` in `[-1, 1]`, I need to compute

`sqrt( (b²*c²) / (1-c²) ) = (|b|*|c|) / sqrt((1-c)*(1+c))`

Catastrophic cancellation appears in the denominator when `c` approaches 1 or -1. The square root probably also does not help.

I was wondering if there is a clever trick I can apply here to avoid the difficult regions around c=1 and c=-1?

The most interesting part of this stability-wise is the denominator, `sqrt(1 - c*c)`. For that, all you need to do is expand it as `sqrt(1 - c) * sqrt(1 + c)`. I don’t think this really qualifies as a “clever trick”, but it’s all that’s needed.

For a typical binary floating-point format (for example IEEE 754 binary64, but other common formats should behave equally well, with the possible exception of unpleasant things like the double-double format), if `c` is close to `1` then `1 - c` will be computed exactly, by Sterbenz’ Lemma, while `1 + c` doesn’t have any stability issues. Similarly, if `c` is close to `-1` then `1 + c` will be computed exactly, and `1 - c` will be computed accurately. The square root and multiplication operations will not introduce significant new error.

Here’s a numerical demonstration, using Python on a machine with IEEE 754 binary64 floating-point and a correctly-rounded `sqrt` operation.

Let’s take a `c` close to (but smaller than) `1`:

``````>>> c = float.fromhex('0x1.ffffffff24190p-1')
>>> c
0.9999999999
``````

We have to be a little bit careful here: note that the decimal value shown, `0.999999999`, is an approximation to the exact value of `c`. The exact value is as shown in the construction from the hexadecimal string, or in fraction form, `562949953365017/562949953421312`, and it’s that exact value that we care about getting good results for.

The exact value of the expression `sqrt(1 - c*c)`, rounded to 100 decimal places after the point, is:

``````0.0000141421362084401590649378320134409069878639187055610216016949959890888003204161068184484972504813
``````

I computed this using Python’s `decimal` module, and double-checked the result using Pari/GP. Here’s the Python calculation:

``````>>> from decimal import Decimal, getcontext
>>> getcontext().prec = 1000
>>> good = (1 - Decimal(c) * Decimal(c)).sqrt().quantize(Decimal("1e-100"))
>>> print(good)
0.0000141421362084401590649378320134409069878639187055610216016949959890888003204161068184484972504813
``````

If we compute naively, we get this result:

``````>>> from math import sqrt
>>> naive = sqrt(1 - c*c)
>>> naive
1.4142136208793713e-05
``````

We can easily compute the approximate number of ulps error (with apologies for the amount of type conversion going on – `float` and `Decimal` instances can’t be mixed directly in arithmetic operations):

``````>>> from math import ulp
>>> float((Decimal(naive) - good) / Decimal(ulp(float(good))))
208701.28298527992
``````

So the naive result is out by a couple of hundred thousand ulps – roughly speaking, we’ve lost around 5 decimal places of accuracy.

Now let’s try with the expanded version:

``````>>> better = sqrt(1 - c) * sqrt(1 + c)
>>> better
1.4142136208440158e-05
>>> float((Decimal(better) - good) / Decimal(ulp(float(good))))
-0.7170147200803595
``````

So here we’re accurate to better than 1 ulp error. Not perfectly correctly rounded, but the next best thing.

With some more work, it ought to be possible to state and prove an absolute upper bound on the number of ulps error in the expression `sqrt(1 - c) * sqrt(1 + c)`, over the domain `-1 < c < 1`, assuming IEEE 754 binary floating-point, round-ties-to-even rounding mode, and correctly-rounded operations throughout. I haven’t done that, but I’d be very surprised if that upper bound turned out to be more than 10 ulps.