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

 
 

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

Submitter:  Colin Macdonald <cbm>
Submitted:  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
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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>
Group 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;
    else
      y = (exp (x) - 1) ./ x;
    endif


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>
Group Member
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:
   
   
Comment:
   

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
  •  

    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 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.13-4448.
    Corresponding source code