bugGNU Octave - Bugs: bug #53543, regression: new gammainc...


bug #53543: regression: new gammainc inaccurate for small real inputs

Submitted by:  Colin Macdonald <cbm>
Submitted on:  Sun 01 Apr 2018 06:38:38 PM UTC  
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  Closed
Release:  dev Operating System:  Any

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


Mon 02 Apr 2018 03:39:37 AM UTC, comment #4: 

I pushed the cset here (http://hg.savannah.gnu.org/hgweb/octave/rev/f98ef2b55641).  I changed the commit message to conform more to Octave standards.  I also marked the newly added test as a regression if it shows up again.

Rik <rik5>
Project Administrator
Mon 02 Apr 2018 02:14:03 AM UTC, comment #3: 
Colin Macdonald <cbm>
Sun 01 Apr 2018 08:47:44 PM UTC, comment #2: 

Dear all,
you are right on the problem, and switching to expm1 solves it (file #43760). I've also added a test.

Michele Ginesi <m_ginesi>
Sun 01 Apr 2018 07:30:47 PM UTC, comment #1: 

I saw a case in gammainc.m:

  elseif (strcmp (tail, "scaledlower"))
    if (abs (x) < 1/2)
      y = expm1 (x) ./ x;
      y = (exp (x) - 1) ./ x;

But there is no equivalent for the case a few lines above

  if (strcmp (tail, "lower"))
    y = 1 - exp (-x);

which is the case that this bug report is about.

Mike Miller <mtmiller>
Project Administrator
Sun 01 Apr 2018 06:38:38 PM UTC, original submission:  

Downstream symbolic pkg: https://github.com/cbm755/octsympy/issues/877

@mtmiller noted:

>> format long
>> x = 1/1000;
>> a = gammainc (x, 1);
>> b = quad (@(t) exp(-t), 0, x);
>> c = double (gammainc (sym(1) / 1000, 1));
>> a, b, c
a =    9.995001666249781e-04
b =    9.995001666250085e-04
c =    9.995001666250085e-04

(in particular, note the results are worse than just using quadrature).

Previously these sorts of things passed with "assert (..., -eps)".

@mtmiller thought "this might be a case of `exp (x) - 1` vs `expm1 (x)`"

Colin Macdonald <cbm>


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

Attach Files:

Attached Files
file #43760:  gammainc.diff added by m_ginesi (1KiB - text/x-patch)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by m_ginesi (Updated the item)
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by cbm (Submitted the item)
  • -email is unavailable- added by cbm
  • -email is unavailable- added by cbm
  • -email is unavailable- added by cbm

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.




    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-04-02 rik5 StatusNone => Fixed
        Open/ClosedOpen => Closed
    2018-04-01 m_ginesi Attached File- => Added gammainc.diff, #43760
    2018-04-01 cbm Carbon-Copy- => Added m_ginesi
        Carbon-Copy- => Added mtmiller
        Carbon-Copy- => Added caliari

    Back to the top

    Powered by Savane 3.5