4

I need to to evaluate the function

$f(x) = {1 - (1-A)^x \over A}$,

where $0 < A \leq 1$ and $0 \leq x \leq 1$.

A straightforward C implementation of $f(x)$ with floating-point arithmetic works fine as long as A is reasonably large, but it breaks down when A is tiny. $(1-A)^x$ is very close to 1, and subtracting it from 1 yields a small number with hardly any significant digits, which is then divided by the tiny A.

$f(x)$ appears to approach an identity function as A approaches zero, but I cannot prove that.

Is there some way to rearrange the definition of $f$ that avoids the loss of precision when A is tiny? Even better, is there a way to make it so that $f(x) = x$ for $A = 0$, without introducing a special case?

abnry
  • 14,664
Stumped
  • 41
  • 1
  • You should do a case structure and taylor expand your function using the binomial taylor expansion. Choose as many terms as you can until you get below machine error. – abnry Aug 30 '13 at 22:22
  • How do I do that? Could you point me to an example? – Stumped Aug 30 '13 at 22:24
  • Can you give an example of $x,A$ for which your computation breaks down? What precision are you working with? – copper.hat Aug 30 '13 at 23:20
  • Here is an idea: Use Ross' suggestion below, and arrange the computation slightly differently. (Just a suggestion, I haven't spent any time analyzing this.) Instead of $\frac{1-(1-A)^x}{A}$, use $x+\frac{1-((1-A)^x+Ax)}{A}$. – copper.hat Aug 30 '13 at 23:36
  • Copper.hat, your suggestion seems to work for a number of $x,A$, but it fails for some. For example, with double precision and $x=0.1,A=10^-16$ I get 0.1, as expected, but for $x=0.5,A=10^-16$ I get 1.6102. (The original version of $f(x)$ produces 0 and 1.1102 respectively.) – Stumped Aug 31 '13 at 00:47
  • That's strange, using Octave I have octave:7> A=1e-16 ; x=.5; x+(1-((1-A)^x+A*x))/A ans = 0.50000 – copper.hat Aug 31 '13 at 01:08

2 Answers2

2

The simple thing to do is to say $(1-A)^x \approx 1-xA$ for $A \ll 1$. This leads to $f(x)=x$ as you say. You can test $A$ and change to $f(x)$ when it is small. The problem this leads to is you have a jump in $f(x)$ as you cross the transition. Maybe this is a problem, maybe not. You can use more terms of the Taylor series: $(1-A)^x \approx 1+x\log(1-A)+\frac 12 (x \log (1-A))^2 + \frac 16 (x \log (1-A))^3 +\ldots$ Then $\log (1-A)=-A-\frac 12A^2-\frac 13A^3-\ldots$. Keeping more terms will reduce the jump. If you need continuity, you can keep a bunch of terms, then near the transition do an interpolation.

Ross Millikan
  • 374,822
1

Construct a polynomial (or rational) approximation of your function that gives good accuracy over the entire interval $[0,1]$, and evaluate this approximation instead of evaluating your original function.

Constructing a good approximation is a fair bit of work, so it may not be worthwhile for you. The best approach is Chebyshev approximation. You should be able to construct a polynomial of reasonably small degree that approximates your function adequately throughout $[0,1]$. With a bit more work, you can arrange for the approximation to be exactly correct at $x=0$ and $x=1$. No special-case branching is needed in your code, so this avoids the discontinuity problem mentioned in Ross Millikan's answer. In general, Chebshev approximations are nice because they are accurate throughout an entire interval, whereas Taylor series approximations are very good at a single point (and nearby) but poor elsewhere.

Constructing the approximation is easy if you have access to Matlab -- you can use the Chebfun package.

There is the added benefit that evaluating a polynomial will probably be faster than evaluating your original function.

bubba
  • 43,483
  • 3
  • 61
  • 122