[Octave-bug-tracker] [bug #62329] Incorrect value (regression) for betainc

[Michael Leitner]

Follow-up Comment #4, bug #62329 (project octave):

I think we can have that easier (less code, better performance, and  more
accurate): note that a general exponentiation

++
x .^ y
--

is computed numerically as

++
exp (y .* log (x))
--

(but see the remark below). Thus, the line 

++
y(a_one) = 1 - (1 - x(a_one)) .^ b(a_one);
--

is just syntactic sugar for 

++
y(a_one) = 1 - exp (log (1 - x(a_one)) .* b(a_one));
--

And now, as we are all perfectly familiar with What Every Computer Scientist
Should Know About Floating-Point Arithmetics, whenever we see such an
expression we should immediately replace each instance of exp(x)-1 with
expm1(x) and log(1+x) with log1p(x), that is, above line is equivalent to 

++
y(a_one) = - expm1 (log1p (- x(a_one)) .* b(a_one));
--

which should always be at least as accurate, but for small arguments much more
accurate. In particular, it solves the present bug.

While we are at it, one should also do the equivalent replacement for the
a_one and b_one cases in the "upper" tail branch. 

Remark: I am not perfectly sure that x.^y is always computed as exp(y.*log(x))
internally. First, for small positive integer exponents y explicit sequential
multiplication and squaring would be faster and more accurate, you can
generalize that to negative integers by a single inversion, and even to
rational exponents with power-of-two denominator by using sqrt (which is very
fast on present architectures). I do not know how far octave internally goes
along this path when it encounters a .^ operator. Perhaps somebody can comment
on that, but if it uses always the exp(log(..)) path, as I see it my proposal
should fix accuracy issues as good as it is possible at all (without doing the
repeated squaring and multiplication here in betainc, which is definitely not
the place to do it). 


    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?62329>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/