# Numerically stable way to compute ugly double fraction

I am looking for a numerically stable version of this (ugly) equation $$s^2=\frac{1}{\frac{1}{\beta_1}+\frac{1}{\beta_2}W}$$ where $$\beta_1 = c_1-c_2m+(m-c_2)b\\ \beta_2 = \frac{1}{2}\left((a-m)^2-(b-m)^2\right)$$ All symbols are real numbers that are stored with floating point accuracy. In my setting, $$\beta_1$$ and $$\beta_2$$ are sometimes close to zero, but never simultaneously.

Clearly, directly computing $$\frac{1}{\beta_1}$$ and $$\frac{1}{\beta_2}$$ might be a bad idea so I am wondering if there is an indirect way, given that $$\beta_1$$ and $$\beta_2$$ share $$b$$ and $$m$$.

I can see that $$\beta_2$$ is close to zero if $$m\approx\frac{1}{2}(a+b)$$ so I was thinking that maybe multiplying $$\frac{1}{\beta_2}$$ with $$\frac{\beta_1}{\beta_1}$$ might help, but I don't see how.

• What about $s^2 = \frac{\beta_1 \beta_2}{\beta_2 + W \beta_1}$? Does that make things better? Jan 24 at 15:56

Addressing the numerical robustness of this expression in full generality, that is, without assumptions about the value range of the variables and their distribution inside those ranges, appears to be a hard problem. The basic problem is that various algebraic re-arrangements leave at least on instance of subtractive cancellation between quantities associated with numerical error.

One standard numerical trick applicable here is that $$a^{2} - b^{2}$$ is best computed as $$(a - b) (a +b)$$. Another common technique is to make use of the fused multiply-add operation (FMA) where $$\mathrm{fma} \ (a, b, c)$$ computes $$a*b+c$$ with a single rounding, and can therefore be used to mitigate instances of subtractive cancellation involving products. First introduced by IBM in the 1990s, FMA is now supported by many processors and is conveniently accessible via a function fma() in most programming environments. Example usage:

$$\beta_{1} = \mathrm{fma} \ (m, b, \mathrm{fma} \ (-c_2, b, \mathrm{fma} \ (-c_2, m, c_1)))$$
$$\beta_{2} = \mathrm{fma} \ (a - b, (a + b) \ / \ 2, (b - a) m)$$
$$d = \mathrm{fma} \ (w, \beta_{1}, \beta_{2})$$
$$s^{2} = \beta_{1} \beta_{2} \ / \ d$$

The above represents my initial attempt. I then spend five hours generating and testing numerous other re-arrangements. I eventually found this slightly more complicated alternative arrangement that appeared somewhat more robust in certain scenarios based on testing with many random operand values:

$$\beta_{1} = \mathrm{fma} \ (m, b, \mathrm{fma } \ (-c_2, b, \mathrm{fma} \ (-c_2, m, c_1)))$$
$$t = \mathrm{fma} \ (-2, m, a + b)$$
$$s^{2} = \beta_{1} (((a - b) t) \ / \ \mathrm{fma} \ (a - b, t, 2 \beta_{1} w))$$

Both of the re-arrangements presented here mitigate but do no fully solve issues of subtractive cancellation. It is possible that there is an FMA-enhanced algebraic re-arrangement that escaped my efforts and provides a universally applicable solution.

• I've tested your initial attempt for over one hour and had only a handful of numeric instabilities. This is good enough for my application, thank you. Jan 25 at 11:34