bugGNU Octave - Bugs: bug #47524, lgamma (Inf) can return complex...

 
 

bug #47524: lgamma (Inf) can return complex result on Windows systems

Submitter:  Avinoam Kalma <avinoam>
Submitted:  Thu 24 Mar 2016 03:58:38 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Fixed Assigned to:  None
Originator Name:  Avinoam Open/Closed:  * Closed
Release:  * dev Operating System:  * Microsoft Windows
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 30 Mar 2016 03:37:09 PM UTC, comment #17: 

Thanks for testing.  I'll at least close this bug report.

Rik <rik5>
Group administrator
Wed 30 Mar 2016 04:56:28 AM UTC, comment #16: 

The lgamma(inf) issue was fixed, thanks :-)

Th sequence


test imread
clear all
test imread


still cause segfault



Avinoam Kalma <avinoam>
Group Member
Tue 29 Mar 2016 03:34:22 PM UTC, comment #15: 

@Avinoam: I've committed a patch to stable and merged to default here (http://hg.savannah.gnu.org/hgweb/octave/rev/6619945e4434).  It would be excellent if you could try building from a cset >= 21559:6619945e4434.  In addition to the fix for lgamma, there is also an important fix for bug #41699, bug #47372.  If you were experiencing the segfault it would be useful to know if that also disappears now.

Rik <rik5>
Group administrator
Tue 29 Mar 2016 04:08:10 AM UTC, comment #14: 

@Rik,

I can build and test a patched
version. I prefer to do it with the dev
branch.
It seems that sgngam is used
uninitialized. (Shouldn't the
compiler warn about it?),
so the first change is:
double sgngam=0;

Avinoam Kalma <avinoam>
Group Member
Mon 28 Mar 2016 10:38:00 PM UTC, comment #13: 

So it is pretty certain then that we are executing the alternate code.  I just disabled the principal code branch by commenting out HAVE_LGAMMA_R in config.h and running the build on a Linux platform.  Sure enough, the alternate code is now running and I can get it to trigger returning a complex matrix.

I think I can prepare a cset against stable.  The long term fix is probably to remove the alternate path entirely and use gnulib to provide a replacement function if the original does not exist.

Rik <rik5>
Group administrator
Mon 28 Mar 2016 08:42:47 PM UTC, comment #12: 

snippets from config.log:


configure:70179: checking for lgamma_r
configure:70179: result: no
.
.
ac_cv_func_lgamma=yes
ac_cv_func_lgamma_r=no
ac_cv_func_lgammaf=yes
ac_cv_func_lgammaf_r=no
.
.




No HAVE_LGAMMA_R defined in the log

John Donoghue <lostbard>
Group Member
Mon 28 Mar 2016 08:19:36 PM UTC, comment #11: 

@JohnD: Could you check the config log for HAVE_LGAMMA_R?  That is the conditional that I think is missing.

Rik <rik5>
Group administrator
Mon 28 Mar 2016 06:59:20 PM UTC, comment #10: 

I don't have the issue on mine.

looking in config.log in the official 4.0.1 installer,


configure:70179: checking for lgamma
configure:70179: i686-w64-mingw32-gcc -o conftest.exe -g -O2 -pthread  -I/scratch/jwe/mxe-octave/usr/i686-w64-mingw32/include -Wl,-rpath-link,/scratch/jwe/mxe-octave/usr/i686-w64-mingw32/lib -L/scratch/jwe/mxe-octave/usr/i686-w64-mingw32/lib conftest.c -lm   -lgdi32 -lws2_32 -luser32 -lkernel32 -lgdi32 -lws2_32 -luser32 -lkernel32 >&5
conftest.c:603:6: warning: conflicting types for built-in function 'lgamma'
 char lgamma ();
      ^
configure:70179: $? = 0
configure:70179: result: yes



and later on:


#define HAVE_LGAMMA 1
#define HAVE_LGAMMAF 1




John Donoghue <lostbard>
Group Member
Mon 28 Mar 2016 05:31:21 PM UTC, comment #9: 

@Avinoam: Can you build the stable branch from Mercurial sources and/or the MXE branch from a tar.gz file?

I think I know what is happening, as expressed in comment #7, but my proposed solution will need testing which requires the ability to patch and then build with MXE.

I've added John D. to the CC list for this bug because he has the most experience with the MXE system.

Finally, I've re-titled the bug report to reflect the true problem.

Rik <rik5>
Group administrator
Mon 28 Mar 2016 04:42:00 PM UTC, comment #8: 


Thanks for the detailed explanation, but I see the problem with
4.0.1 official MXE-version, which I did not compile.
Is it possible that other/wrong function enters the path?

Avinoam Kalma <avinoam>
Group Member
Mon 28 Mar 2016 04:22:36 PM UTC, comment #7: 

gammaln is an alias for lgamma.  And lgamma is a mapper on to the same function in libc.  So it it possible that the libc for cygwin, mingw, or whatever you are using is incorrect.  Eventually the function is defined in liboctave/numeric/lo-specfun.cc.  I've extracted what look like the relevant functions. 


double
xlgamma (double x)
{
#if defined (HAVE_LGAMMA)
  return lgamma (x);
#else
  double result;
  double sgngam;

  if (xisnan (x))
    result = x;
  else if ((x <= 0 && D_NINT (x) == x) || xisinf (x))
    result = octave_Inf;
  else
    F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));

  return result;
#endif
}

Complex
rc_lgamma (double x)
{
  double result;

#if defined (HAVE_LGAMMA_R)
  int sgngam;
  result = lgamma_r (x, &sgngam);
#else
  double sgngam;

  if (xisnan (x))
    result = x;
  else if ((x <= 0 && D_NINT (x) == x) || xisinf (x))
    result = octave_Inf;
  else
    F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));

#endif

  if (sgngam < 0)
    return result + Complex (0., M_PI);
  else
    return result;
}


You might check your config.h and see whether HAVE_LGAMMA and HAVE_LGAMMA_R are defined.  If they aren't defined, then you might be going through the backup definition which uses Fortran code.  As you can see, in that case there is a possibility of returning an extra complex PI value which is what you are seeing.


  if (sgngam < 0)
    return result + Complex (0., M_PI);


My suspicion is that you are taking the alternate path, and the declaration "double sgngam" doesn't initialize sgngam.  There is a shortcut path if x is Inf which initializes the return value, but doesn't initialize sgngam.  I think when you reach the if/else statement outside of the #if/#endif block you are essentially randomly sampling a memory location when you ask "if (sgngam < 0)".

So there are several things to check.  First check config.h.  Second, try altering the code to


#else
  double sgngam = 0;


and building, creating a tar.gz file, and then using that file with MXE to build a Windows version.  If that works then we can make the change, or look in to having gnulib supply these functions for cross-platform portablity.  I see that these two systems do not supply lgamma_r: mingw, MSVC 9.  That seems relevant.

Rik <rik5>
Group administrator
Mon 28 Mar 2016 03:32:57 PM UTC, comment #6: 

It works for me also on two other Win-7 machines,
but on the first one I get


>> gammaln(Inf)
ans = Inf + 3.142i


Which causes the problem. How can I debug it, and why the result is different in two other computers?

Avinoam Kalma <avinoam>
Group Member
Mon 28 Mar 2016 04:16:28 AM UTC, comment #5: 

I have checked this on Windows-7 machine.
Checking with debugger, I have found that the
complex comes from the line


   pdf(k) = exp (- a(k) .* log (b(k)) + (a(k)-1) .* log (x(k))
                  - x(k) ./ b(k) - gammaln (a(k)));


from the gammaln term

I will check later in another win7 machine

Avinoam Kalma <avinoam>
Group Member
Sun 27 Mar 2016 09:39:44 PM UTC, comment #4: 

So this may be a difference in Operating Systems, rather than a bug with Octave.  According to Hartmut, it works with Windows 7.  Avinoam, what version of Windows are you using?

Rik <rik5>
Group administrator
Sun 27 Mar 2016 06:52:12 PM UTC, comment #3: 

I have also tested "test gampdf" under Windows 7 with

  • the Octave 4.0.1 release, and
  • a mxe-build of the dev branch (4.1.0+) from yesterday


and as a result I get in both cases 19 passes out of 19.

Hartmut <hardy>
Fri 25 Mar 2016 04:27:02 AM UTC, comment #2: 



>> x = [-1, 0, 0.5, 1, Inf];
>> y = gampdf (x, [0, -Inf, NaN, Inf, 1], 1)
y =

   NaN +   0i   NaN +   0i   NaN +   0i   NaN + NaNi     0 +   0i

>> typeinfo (y)
ans = complex matrix


Avinoam Kalma <avinoam>
Group Member
Thu 24 Mar 2016 05:48:09 PM UTC, comment #1: 

Can you find out exactly what the output is?

Try


x = [-1, 0, 0.5, 1, Inf];
y = gampdf (x, [0, -Inf, NaN, Inf, 1], 1)
typeinfo (y)



Rik <rik5>
Group administrator
Thu 24 Mar 2016 03:58:38 PM UTC, original submission:  


test gampdf fails both in 4.0.1 and in the dev. branch, under windows:



>> test gampdf
***** assert (gampdf (x, [0 -Inf NaN Inf 1], 1), [NaN NaN NaN NaN y(5)])
!!!!! test failed
ASSERT errors for:  assert (gampdf (x, [0, -Inf, NaN, Inf, 1], 1),[NaN, NaN, NaN, NaN, y(5)])

  Location  |  Observed  |  Expected  |  Reason
     ()           O            E         complex != real
shared variables
  scalar structure containing the fields:

    x =

      -1.00000   0.00000   0.50000   1.00000       Inf

    y =

       0.00000   1.00000   0.60653   0.36788   0.00000


Avinoam Kalma <avinoam>
Group Member

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by lostbard (Posted a comment)
  • -email is unavailable- added by rik5
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by avinoam (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 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-03-30 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2016-03-29 rik5 StatusConfirmed Ready For Test
    2016-03-28 rik5 StatusWorks For Me Confirmed
        SummaryMXE-Octave: test gampdf fails lgamma (Inf) can return complex result on Windows systems
    2016-03-28 rik5 Carbon-Copy- Added lostbard
    2016-03-27 rik5 StatusNeed Info Works For Me
    2016-03-24 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code