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.