bugGNU Octave - Bugs: bug #43551, Matlab incompatibility for gamma ...

 
 

bug #43551: Matlab incompatibility for gamma function for negative integer arguments

Submitter:  John W. Eaton <jwe>
Submitted:  Thu 06 Nov 2014 08:07:42 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Matlab Compatibility
Status:  Fixed Assigned to:  None
Originator Name:  jwe Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 12 Nov 2014 03:54:29 PM UTC, comment #9: 

@jwe: I think you can go ahead and apply the patch and close down this report.

Rik <rik5>
Group administrator
Sun 09 Nov 2014 01:26:16 AM UTC, comment #8: 

Or, to quote Abramowitz and Stegun:

Gamma(z) is single valued and analytic over the entire
complex plane, save for the points z = -n (n+0, 1, 2,...)
where it possesses simple poles with residue (-1)^n/n!.

As said before, Inf will do.

Michael Godfrey <godfrey>
Group Member
Fri 07 Nov 2014 09:07:06 PM UTC, comment #7: 

The value is undefined at negative integer values in the same way that 1/0 is undefined.

lim    1/x == -Inf
x->0-

lim
x->0+  1/x == +Inf

The magnitude is infinite, but the phase can't be determined.  We could return NaN in this case to be mathematically precise, but nobody bothers to.  This suggests that we drop mathematical correctness and just go with the rest of the crowd.

Rik <rik5>
Group administrator
Fri 07 Nov 2014 08:53:26 PM UTC, comment #6: 

John:

The patch looks right to me, assuming that it
uses Inf or the negative poles. This is the usual
assumption in the math/computation (ACM, etc.)
literature and coded algorithms. Not just Matlab.

NaN is not obviously a good substitute for "not defined"
in this case.

Michael Godfrey <godfrey>
Group Member
Fri 07 Nov 2014 08:36:13 PM UTC, comment #5: 

I propose the attached patch.  Comments?

(file #32391)

John W. Eaton <jwe>
Group administrator
Fri 07 Nov 2014 05:11:04 PM UTC, comment #4: 

Isn't the value really undefined at exact negative integer values?  So for the best math, we should have NaN.  If Inf is more convenient (but mathematically incorrect regardless of sign) then we might as well just be compatible with Matlab and make it +Inf always.

John W. Eaton <jwe>
Group administrator
Fri 07 Nov 2014 05:06:29 PM UTC, comment #3: 

That's kind of what I expected.  It does seem like it would be more useful to return Inf and have terms go to zero.

I was just picking up the patch that someone had supplied so I don't know much about the code.

The mod 2 thing seems to be trying to get the phase to be consistent for negative integers.  The result of the gamma function is only infinite when x is a negative integer (or when X itself is +Inf).  So the code, by checking xisinf (result), is effectively deciding that x itself was a negative integer.  The next test will determine whether x was an even or odd integer and will flip the phase (+1,-1,+1,-1) to go with the absolute magnitude of infinity.  This corresponds to consistently taking the limit from the right.

Rik <rik5>
Group administrator
Fri 07 Nov 2014 04:37:26 PM UTC, comment #2: 

The example that resulted in someone pointing out the change to me was one in which gamma(x) was in the denominator of an expression inside a summation.  Previously, these terms were then zero and not causing trouble.  With the new behavior, they are NaN and cause trouble because then the summation is NaN as well.

If we choose Inf as compatible with Matlab and perhaps more useful, then I think we can still treat -0 specially.

One more question:  what is the purpose of the lines


      if (xisinf (result) && (static_cast<int> (gnulib::floor (x)) % 2))
        result = -octave_Inf;


now?  Given the current checks at the top of the function, is this condition ever met?  And why would we make the result negative only at odd values?

John W. Eaton <jwe>
Group administrator
Fri 07 Nov 2014 04:25:51 PM UTC, comment #1: 

I guess we should think about which form is more useful.

The Wikipedia entry for the gamma function has a nice plot of the function (http://en.wikipedia.org/wiki/Gamma_function).  The magnitude of the function at a negative integer is infinite.  But, if you choose a consistent direction, such as always approaching the pole from the left, then the sign of the result flips back and forth between +1 and -1.  Wolfram Alpha, a.k.a Mathematica, returns complex infinity for gamma (negative integer).  A complex infinity has infinite magnitude, but no determined phase.  This is consistent with the plot which shows that you are either getting +Inf or -Inf in the limit.

So, is it more useful to simply pick a phase (+1) and return Inf as Matlab does?  Or is it more useful to throw a NaN which will clog up downstream calculations and indicate that we have gone beyond our ability to truly calculate the result?

Rik <rik5>
Group administrator
Thu 06 Nov 2014 08:07:42 PM UTC, original submission:  

The change made because of the following patch submission results in Matlab-incompatible behavior.

  https://savannah.gnu.org/patch/?func=detailitem&item_id=8101

Apparently Matlab's gamma function returns Inf for negative integer values and for 0 (Matlab has no -0).

Is that important for us?

If we decide to keep the current behavior, then I think we should at least mention that the difference with Matlab is intentional.

John W. Eaton <jwe>
Group administrator

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #32391:  diffs.txt added by jwe (2KiB - text/plain)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by godfrey (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by jwe (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 3 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-12-01 rik5 StatusNone Fixed
        Open/ClosedOpen Closed
    2014-11-07 jwe Attached File- Added diffs.txt, #32391

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code