0

The function

$$f(a,c) = \frac{1}{a}((1-c)^a - 1), $c \in [0,1], a < 1$$

is smooth around $a = 0$ - the limit exists as $\lim_{a \rightarrow 0} f(a,c) = log(1 - c)$. However, for $|a| \ll 1$ the formula becomes numerically unstable, as for $|a| < 10^{-17}$ the numerator evaluates to exactly $0$ (and the behavior is weird from around $|a| = 10^{-15}$. Are there any tricks I could use here?

I've tried using expm1 to evaluate $f(a,c) = \frac{1}{a} ( \exp(a \log(1 - c)) - 1) = \frac{1}{a} \mathbb{expm1}(a \log(1-c))$ but this actually breaks sooner (around $|a| = 10^{-6}$)

Below is an R code showing the problematic behavior (we would expect the function to be monotonically increasing in $a$ over this range):

library(tidyverse)
tibble(a = c(10^(seq(-1, -20))), c = 0.5) %>% rowwise() %>% mutate(orig = ((1.0-c)^a - 1.0) /a, expm1 = expm1(a * log(1 - c))/a)  %>% arrange(a)

but we get:

              a     c       orig      expm1
          <dbl> <dbl>      <dbl>      <dbl>
 1 1.000000e-20   0.5  0         -0.6931472
 2 1.000000e-19   0.5  0         -0.6931472
 3 1.000000e-18   0.5  0         -0.6931472
 4 1.000000e-17   0.5  0         -0.6931472
 5 1.000000e-16   0.5 -1.110223  -0.6931472
 6 1.000000e-15   0.5 -0.6661338 -0.6931472
 7 1.000000e-14   0.5 -0.6883383 -0.6931472
 8 1.000000e-13   0.5 -0.6927792 -0.6931472
 9 1.000000e-12   0.5 -0.6931122 -0.6931472
10 1.000000e-11   0.5 -0.6931455 -0.6931472
11 1.000000e-10   0.5 -0.6931467 -0.6931472
12 1.000000e- 9   0.5 -0.6931472 -0.6931472
13 1.000000e- 8   0.5 -0.6931472 -0.6931472
14 1.000000e- 7   0.5 -0.6931472 -0.6931472
15 1.000000e- 6   0.5 -0.6931469 -0.6931469
16 1.000000e- 5   0.5 -0.6931448 -0.6931448
17 1.000000e- 4   0.5 -0.6931232 -0.6931232
18 1.000000e- 3   0.5 -0.6929070 -0.6929070
19 1.000000e- 2   0.5 -0.6907505 -0.6907505
20 1.000000e- 1   0.5 -0.6696701 -0.6696701
Martin Modrák
  • 187

1 Answers1

1

For $c\in\{0,1\}$ the calculation is trivial. Otherwise we can use $$\begin{align} \frac{(1-c)^a-1}{a} &=\sum_{k=0}^\infty\frac{a^k\ln^{k+1}{(1-c)}}{(k+1)!}\\ &=\ln{(1-c)}+\frac{a\ln^2{(1-c)}}{2}+\frac{a^2\ln^3{(1-c)}}{6}+\dots\\ \end{align}$$ which can be stably computed by adding terms until the value of the sum is unchanged in the precision used.

Peter Foreman
  • 19,947