How can I rearrange this logarithmic formula to be computer friendly? I've had a look through the logarithmic identities on Wikipedia, but nothing fits the bill.
Basically, I have a formula which shows how much more 'risky' one number is compared to another, where 0 = safe and 1 = not safe for the variables a and b.
The formula is: $log_{1-b}(1 - a)$
The C# code looks like: Math.Log(1 - a, 1 - b) for what it's worth.
This works fine when a or b isn't an extraordinarily low number, but as soon as it becomes lower than around 1E-16, the numbers go haywire. This is of course because double types in the computer's memory are not good at holding $1.0 - 1E-17 = 0.99999...$ <17 times>.
For example, the sum:
$log_{1 - 1E-15}(1 - 0.001) = 1.001E12$ or even:
$log_{1 - 1E-16}(1 - 0.001) = 9.01E12$ is fine, but as soon as I go an order of magnitude deeper:
$log_{1 - 1E-17}(1 - 0.001)$, I get NaN. This is despite that there's plenty more room in the double datatype to allow a far greater exponent for the final answer.
That's for b dipping below 1E-16, but a similar thing happens when a dips below 1E-16. But instead of NaN, I get a (wrong) answer of zero.
So I'm trying to find a way to reformulate the expression so that it's kinder to the computer. Any ideas how?
 A: First we have the formula $\log_{1-b}(1-a)=\frac{\ln(1-a)}{\ln(1-b)}$. This reduces your question to computing ln(1+x) where $|x| \ll 1$.
The naive implementation can be most easily seen to fail by considering |x| less than machine epsilon. In this case 1+x==1 (in the sense of floating point ==). There is no way to make this not be the case, it is a property of floating point arithmetic that we are stuck with. So if |x| is less than machine epsilon then log(1+x)==0 in floating point arithmetic, which is a 100% relative error. Additionally, for non-representable numbers, there is still a roundoff problem even when |x| is just small but greater than machine epsilon. So the naive implementation is no good once |x| is fairly small, even when |x| is quite a bit larger than machine epsilon.
The way around this issue is to use a different function which computes $\ln(1+x)$ given x without computing 1+x. Such a function is called log1p in quite a few languages, including Matlab, C++ (Math.h) and Python (numpy).
A simple way to get similar (but not identical) behavior, suggested at http://www.johndcook.com/blog/2010/06/07/math-library-functions-that-seem-unnecessary/ is to use the usual log function for $|x|>10^{-4}$ and use a Taylor expansion otherwise. A four term Taylor expansion should be expected to nearly achieve floating point precision over this range, because the Lagrange remainder tells us that the relative error in the four term Taylor expansion carried out in exact arithmetic is less than $10^{-16}$. Note that this is a relative error, which is the right kind of error for your problem where you need to divide one logarithm by another.
To do better you will need to do something significantly different, because if $|x|<10^{-4}$ then $1+|x|^4$ is floating point equal to 1, so introducing a fifth term will not improve the approximation.
A: Write
$$
\log_{1-b}(1 - a) = \frac{\log(1-a)}{\log(1-b)}
$$
Many programming languages have a log1p function which computes $\log(1+x)$ accurately for $x$ small. If your language has log1p, definitely use it.
