-1

so I have recently been working with some probability distributions of order statistics. in that area it is quite common to see lots of formulas with high powers of numbers from the real interval [0, 1].

Consider numbers a ~ b ~ 0, both positive. I want to calculate something like a^n / b^m where n, m are huge numbers.

My question is the following: If I use C code like 

double f(double a, double b, long m, long n)
{
    return( pow(a, n) / pow(b, m) );
}

will this be stable? The rule should be to first evaluate a^n, then b^m then divide, however if a^n or b^m is small enough it will just be zero or NaN. Instead I could do something like

double g(double a, double b, long m, long n)
{
    double res;
    long i;
    res = 1;
    for (i = 0; i < min(m, n); ++i)
    {
        res = res * a / b;
    }
    if ( n > m )
    {
        res = res * pow(a, n - m);
    } else {
       res = res / pow(b, m - n);
    }
    return( res );
}

Do you know whether 1) is subject to optimization in such cases? If not, how would one handle such cases for a fast and stable evaluation?

user207421
  • 305,947
  • 44
  • 307
  • 483
Jack
  • 195
  • 2
  • 10
  • 1
    To avoid such scenarios, would it be viable to compute in the log domain? i.e. `result = exp(n*log(a) - m*log(b))`. (Disclaimer: I'm definitely not a numerical expert.) – Oliver Charlesworth Oct 20 '13 at 13:11
  • You want the second form. But don't use a `for` loop. Use `pow(a/b, min(m, n))` to reduce the number of operations. (Assuming `m` and `n` must be non-negative.) – David Schwartz Oct 20 '13 at 13:35
  • Well, thanks! I guess there is no problem with combining both approaches using sum() and log(a)-log(b) :) – Jack Oct 22 '13 at 14:05
  • In your `a ~ b ~ 0`, is `~` supposed to denote approximate equality? – Keith Thompson Oct 25 '13 at 01:13

1 Answers1

3

If you're asking whether pow() implements exponentiation by repeated multiplication, the answer is no. You can assume it's optimized, so any savings you might achieve by checking for a zero numerator is probably negligible. You definitely don't want to implement this with a for() loop that iterates through a huge number of floating-point calculations.

On the other hand, it's good practice to avoid dividing by zero.

I'd implement your function like this:

double f(double a, double b, long m, long n) {
  if (b == 0) {
    // return error or infinity.
  }
  if (m >= n)
    return pow(a / b, n) * pow(a, m - n);
  else 
    return pow(a / b, m) / pow(b, n - m);
  }
}
Adam Liss
  • 47,594
  • 12
  • 108
  • 150