bugGNU Octave - Bugs: bug #62329, Incorrect value (regression) for...


bug #62329: Incorrect value (regression) for betainc

Submitter:  Julien Bect <jbect>
Submitted:  Mon 18 Apr 2022 05:33:52 AM UTC
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Regression
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 7.1.0 Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup

Jump to the original submission

Wed 20 Apr 2022 10:24:54 PM UTC, comment #8: 

Excellent analysis and confirmatory testing by Michael in comment #7.  I think this bug is resolved so I am marking as Fixed and closing report.

Rik <rik5>
Group administrator
Wed 20 Apr 2022 03:06:49 PM UTC, comment #7: 

I did some research, and std::pow() seems to not even come from the library, but typically be provided by the compiler itself (at sufficient optimization level). Thus, I did not find anything authoritative on how it is really implemented, and therefore I tested it:

I used


to have nice doubles that also with high powers do not become sub- or supranormal (or inf), which could necessitate special treatment. Then it is just a question of


for various b. What I found on a Intel i5-2500 with 64-bit Debian is that for any (also negative) non-integer exponent (of course real) it takes about 180 clock cycles per such power operation, which is perfectly consistent with the combined cost of exp() and log() (the latter on its own takes about 100 cycles, exp() about 40, and then some multiplying and adding). Integer exponents of large absolute value are the same. The following cases are treated differently (more efficiently): exponents of 2 and 3, which take about 15 cycles, exponents of 0 and 1, which take about 45 cycles, and exponent -1 (equivalent in terms of speed to 1./a) takes 40 cycles. This seems a bit much, as in compiled C programs on my architecture I can do inversion in about 20 cycles. And I have no idea whatsoever how the cases of 0 and 1 are implemented -- of course both are trivial and should involve no calculation whatsoever, but on the other hand there is no reason ever to use them.

To compare, doing a.*a in octave is also about 15 cycles, a.*a.*a correspondingly 30 cycles.

Bottom lines:
- the code path that is used by the x.^y operator (equivalent to the power(x,y) function) on my architecture results in about a factor of ten more efficient computation for y==2 or 3 compared to the general case, and thus also compared to the proposed fix. But the expm1 and log1p is unconditionally stable, thus I would propose to leave it.
- it seems that there is an opportunity for improvement of some of the lowest-level operations in octave, which could be interesting for Rik's long-term ambition: First, could inversion possibly be sped up? As I see it, there are not more checks to be performed than for multiplication, because 1./0 is inf automatically, isn't it? But the overhead of octave compared to C is 20 cycles per operation, more than .^2 takes in total (in octave). And further, .^ with integer exponents could definitely be made more efficient for a wide range of exponents: by repeated squaring I can take the power of, e.g., 127 within 11 cycles. Perhaps sometime I come round to write an oct-file implementation as demonstration.

Michael Leitner <mleitner>
Tue 19 Apr 2022 10:31:09 PM UTC, comment #6: 

I believe for Octave

x .^ y

is not computed in any specific way but is instead handed off to the C/C++ standard library and specifically the std::pow() routine.  See the file libinter/corefcn/xpow.cc.

This is a bit of a team effort, but I gave the checkin to Nir.  See http://hg.savannah.gnu.org/hgweb/octave/rev/fe898ae23e0e.

Marking as Ready for Test.

Rik <rik5>
Group administrator
Tue 19 Apr 2022 03:09:43 PM UTC, comment #5: 

Thanks Michael, I forgot about those Octave functions.

Here's a revised patch.

(file #53113, file #53114)

Nir Krakauer <nir_krakauer>
Tue 19 Apr 2022 10:04:54 AM UTC, comment #4: 

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).

Michael Leitner <mleitner>
Tue 19 Apr 2022 12:35:17 AM UTC, comment #3: 

Here is a patch that avoids subtractive cancellation for small x when a = 1.

(file #53107, file #53108)

Nir Krakauer <nir_krakauer>
Mon 18 Apr 2022 09:03:31 PM UTC, comment #2: 

This is certainly caused by the change introduced for bug #60682.  Adding Michele Ginesi to the CC list.

The code in question is the fast path code when a==1.

I(a_one) = 1 - (1 - x(a_one)) .^ b(a_one);

The specific issue is

(1 - x(a_one))

x is 2e-20 which is less than eps(1) so in the subtraction the result is just a perfect 1.

Rik <rik5>
Group administrator
Mon 18 Apr 2022 02:02:30 PM UTC, comment #1: 

Confirmed also on Windows with Octave 7.1.0.

Could be related to the changes for bug #60682.

Fwiw, in Matlab R2021a:

>> betainc(2e-20, 1, 0.5)

ans =


Markus M├╝tzel <mmuetzel>
Group administrator
Mon 18 Apr 2022 05:33:52 AM UTC, original submission:  

The following value is incorrect:

>> betainc(2e-20, 1, 0.5)
ans = 0

The result was correct in Octave 6.2.0:

octave:5> betainc(2e-20, 1, 0.5)
ans = 1.0000e-20

The value 1e-20 is also the one I get in Python/scipy (scipy.special.betainc) and R (pbeta).

Julien Bect <jbect>


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

Attached Files
file #53113:  betainc.patch added by nir_krakauer (2KiB - text/x-patch)
file #53114:  betainc.m added by nir_krakauer (9KiB - text/x-objcsrc)
file #53107:  betainc.patch added by nir_krakauer (1KiB - text/x-patch)
file #53108:  betainc.m added by nir_krakauer (9KiB - text/x-objcsrc)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by svillemot
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by nir_krakauer (Updated the item)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by jbect (Submitted the item)

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.


    Follow 12 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-04-20 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-04-19 rik5 StatusConfirmed Ready For Test
    2022-04-19 nir_krakauer Attached File- Added betainc.patch, #53113
        Attached File- Added betainc.m, #53114
    2022-04-19 svillemot Carbon-Copy- Added svillemot
    2022-04-19 nir_krakauer Attached File- Added betainc.patch, #53107
        Attached File- Added betainc.m, #53108
    2022-04-18 rik5 Carbon-Copy- Added m_ginesi
    2022-04-18 mmuetzel CategoryLibraries Octave Function
        StatusNone Confirmed
        Operating SystemGNU/Linux Any

    Back to the top

    Powered by Savane 3.13-ee12.
    Corresponding source code