Proposal changeset of liboctave/lo-specfun.cc (was Re: gammaln(0) (octav

octave-maintainers
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Proposal changeset of liboctave/lo-specfun.cc (was Re: gammaln(0) (octav

From:	Tatsuro MATSUOKA
Subject:	Proposal changeset of liboctave/lo-specfun.cc (was Re: gammaln(0) (octave-3.3.51+, MinGW))
Date:	Tue, 3 Aug 2010 10:19:39 +0900 (JST)
Hello

To overcome error on gammaln(0) or gammaln(minus integer) using slatec gamma 
function procedures, I 
propose attached changeset for the current development source (octave-3.3.52+)

Regards

Tatsuro



--- Tatsuro MATSUOKA  wrote:

> Hello
> 
> Sorry for my frequently post,
> 
> I have misled.
> In octave-3.2.4/mingw on the octave-forge
> 
> octave:2> lgamma(-1)
>  ***MESSAGE FROM ROUTINE DGAMMA IN LIBRARY SLATEC.
>  ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
>  *  X IS A NEGATIVE INTEGER
>  *  ERROR NUMBER = 4
>  *
>  ***END OF MESSAGE
> 
>  ***JOB ABORT DUE TO FATAL ERROR.
> 0          ERROR MESSAGE SUMMARY
>  LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     COUNT
>  SLATEC     DGAMMA     X IS A NEGATIVE INTE         4         2         1
> 
> error: exception encountered in Fortran subroutine dlgams_
> 
> Perhaps, in this case slatec gamma is also used.
> 
> I cannot understand why lgamma(0) is not broken.
> 
> 
> Seeing the code, liboctave/lo-specfun.cc, I found
> 
> double
> xgamma (double x)
> {
> #if defined (HAVE_TGAMMA)
>   return tgamma (x);
> #else
>   double result;
> 
>   if (xisnan (x))
>     result = x;
>   else if ((x <= 0 && D_NINT (x) == x) || xisinf (x))
>     result = octave_Inf;
>   else
>     F77_XFCN (xdgamma, XDGAMMA, (x, result));
> 
>   return result;
> #endif
> }
> 
> For also xlgamma, is the below to be used
>   else if ((x <= 0 && D_NINT (x) == x) || xisinf (x)) 
> instead if  
>    else if (xisinf (x))
> 
> Am I right?
> 
> Regards
> 
> Tatsuro
> 
> --- Tatsuro MATSUOKA <address@hidden> wrote:
> 
> > Hello 
> > 
> > In liboctave/lo-specfun,
> > 
> > I found as:
> > 
> > double
> > xlgamma (double x)
> > {
> > #if defined (HAVE_LGAMMA)
> >   return lgamma (x);
> > #else
> >   double result;
> >   double sgngam;
> > 
> >   if (xisnan (x))
> >     result = x;
> >   else if (xisinf (x))
> >     result = octave_Inf;
> >   else
> >     F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
> > 
> >   return result;
> > #endif
> > }
> > 
> > As I posted previously,
> > 
> > I found the below definition,
> > 
> > #define HAVE_LGAMMA 1
> > #define HAVE_LGAMMAF 1
> > 
> > However it seemed to be ignored, and therefore stalec was used.
> > 
> > Hmmmm ??
> > 
> > Apart from the failure in detecting the definition, the above description 
> > should be modified,
> > because
> > stlatec gamma functions do not allow x=0.0, and beta.m assume that 
> > gammaln(0) is infinity as
> > shown in
> > their own test as
> > 
> > %!test
> > %! a = 0.25 + (0:5) * 0.5;
> > %! tol = 10 * max (a) * eps;
> > %! assert (zeros (size (a)), beta (a, -a), tol)
> > %! assert (zeros (size (a)), beta (-a, a), tol)
> > 
> > Therefore in the case of slatec gamma function, perhaps x==0 should be 
> > treated
> > 
> >   if (xisnan (x))
> >     result = x;
> >   else if (xisinf (x))
> >     result = octave_Inf;
> > +  else if (x==0.0)
> > +    result = octave_Inf;
> >   else
> >     F77_XFCN (dlgams, DLGAMS, (x, result, sgngam));
> > 
> >   return result;
> > 
> > Is the above correct ?  It should be like 'abs(x)<tol' ?
> > 
> > 
> > Regards
> > 
> > Tatsuro
> > 
> >   
> > 
> > --- Tatsuro MATSUOKA wrote:
> > 
> > > Hello
> > > 
> > > > In gfortran DLGAMA() is intrinsically defined.  
> > > > 
> > > > If it is used, perhaps the error will be disappeared.
> > > > http://gcc.gnu.org/onlinedocs/gfortran/LOG_005fGAMMA.html#LOG_005fGAMMA
> > > > 
> > > > In the case slatec gamma function, we cannot use it at x=0 so that 
> > > > special treatment is
> > > > required, for
> > > > gammaln in octave, I think.
> > > > 
> > > > Am I wrong?
> > > 
> > > In gcc, lgamma, i.e.log(abs(gamma)) is defined. (also lgammaf) 
> > > I have looked into "config.log" and found
> > > configure:47767: checking for lgamma
> > > <snip>
> > > configure:47767: result: yes
> > > configure:47767: checking for lgammaf
> > > <snip>
> > > configure:47767: result: yes
> > > 
> > > And in config.h, there exist
> > > 
> > > /* Define to 1 if you have the `lgamma' function. */
> > > #define HAVE_LGAMMA 1
> > > 
> > > /* Define to 1 if you have the `lgammaf' function. */
> > > #define HAVE_LGAMMAF 1
> > > 
> > > 
> > > Perhaps for gamma() in the octave-3.3.51+ built by me, lgamma is used.  
> > > However, for gammaln() in, subroutine dlgams in slatec is used.
> > > 
> > > Because I can see in the error message,
> > >  error: lgamma: exception encountered in Fortran subroutine dlgams_
> > > 
> > > Is there any reason that buid-in lgamma in gcc is not used but the slatec 
> > > dlgams is used,
> > > although
> > > 'configure' detects it ?
> > > 
> > > Regards
> > > 
> > > Tatsuro
> > > 
> > > --- Tatsuro MATSUOKA  wrote:
> > > 
> > > > Hello
> > > > 
> > > > > I cannot find source of xstopx in libcruft source files.
> > > > > 
> > > > > However when the code link against libcruft.dll.a, undefined symbol 
> > > > > error messages for
> > > > 'xstopx'
> > > > > were
> > > > > disappeared.
> > > > 
> > > > I found the reason. xstopx is defined in   misc/f77-fcn.h in libcruft.
> > > > 
> > > > > However
> > > > > 
> > > > >       program test
> > > > >       DOUBLE PRECISION X, DLGAM, SGNGAM, DLNGAM
> > > > >       X=0
> > > > >       write(*,*) DGAMMA(X)
> > > > >       end
> > > > > 
> > > > > The above code gives
> > > > > 
> > > > > $ ./test
> > > > >                  +Infinity
> > > > 
> > > > This is reasonable.  In this case, dgamma intrinsic in gfortran  is 
> > > > used.
> > > > If I add 'EXTERNAL DGAMMA', slatac dgamma is used and the same error as 
> > > > that by the dlgams
> > > > subroutime
> > > > is caused.
> > > > 
> > > > In gfortran DLGAMA() is intrinsically defined.  
> > > > 
> > > > If it is used, perhaps the error will be disappeared.
> > > > http://gcc.gnu.org/onlinedocs/gfortran/LOG_005fGAMMA.html#LOG_005fGAMMA
> > > > 
> > > > In the case slatec gamma function, we cannot use it at x=0 so that 
> > > > special treatment is
> > > > required, for
> > > > gammaln in octave, I think.
> > > > 
> > > > Am I wrong?
> > > > 
> > > > Regards
> > > > 
> > > > Tatsuro
> > > > 
> > > > 
> > > > 
> > > > --- Tatsuro MATSUOKA  wrote:
> > > > 
> > > > > Hello
> > > > > 
> > > > > I have reported error in test beta.m.
> > > > > Development source on July 2010 (platform MinGW, GCC-4.5.0)
> > > > > 
> > > > > 
> > > > > http://octave.1599824.n4.nabble.com/error-in-test-beta-m-tc2303186.html#a2303407
> > > > > 
> > > > > This problem reduced to 'gammaln(0)'
> > > > > 
> > > > > octave:5> gammaln(0)
> > > > >  ***MESSAGE FROM ROUTINE DGAMMA IN LIBRARY SLATEC.
> > > > >  ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
> > > > >  *  X IS 0
> > > > >  *  ERROR NUMBER = 4
> > > > >  *
> > > > >  ***END OF MESSAGE
> > > > > 
> > > > >  ***JOB ABORT DUE TO FATAL ERROR.
> > > > > 0          ERROR MESSAGE SUMMARY
> > > > >  LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     
> > > > > COUNT
> > > > >  SLATEC     DGAMMA     X IS 0                       4         2       
> > > > >   2
> > > > > 
> > > > > error: lgamma: exception encountered in Fortran subroutine dlgams_
> > > > > 
> > > > > I can fined symbol 'dlgams'
> > > > > 
> > > > > I have been trying to build simple test file using DLGAMS in 
> > > > > slatec-fn like
> > > > > 
> > > > >       program test
> > > > >       DOUBLE PRECISION X, DLGAM, SGNGAM, DLNGAM
> > > > >       X=0
> > > > >       call DLGAMS (X, DLGAM, SGNGAM)
> > > > >       write (*,*) DLGAM, SGNGAM
> > > > >       end
> > > > > 
> > > > > In the slatec-err routines (or other libcruft source files), I found 
> > > > > 'call
> > > xstopx('messeage')'
> > > > > but
> > > > > I cannot find source of xstopx in libcruft source files.
> > > > > 
> > > > > However when the code link against libcruft.dll.a, undefined symbol 
> > > > > error messages for
> > > > 'xstopx'
> > > > > were
> > > > > disappeared.
> > > > > 
> > > > > And test result is 
> > > > > 
> > > > > $ gdb test
> > > > > GNU gdb (GDB) 7.0.1
> > > > > Copyright (C) 2009 Free Software Foundation, Inc.
> > > > > License GPLv3+: GNU GPL version 3 or later 
> > > > > <http://gnu.org/licenses/gpl.html>
> > > > > This is free software: you are free to change and redistribute it.
> > > > > There is NO WARRANTY, to the extent permitted by law.  Type "show 
> > > > > copying"
> > > > > and "show warranty" for details.
> > > > > This GDB was configured as "mingw32".
> > > > > For bug reporting instructions, please see:
> > > > > <http://www.gnu.org/software/gdb/bugs/>...
> > > > > Reading symbols from 
> > > > > c:\usr\Tatsu\mingwhome\octaves\dgammatest/test.exe...done.
> > > > > (gdb) br test
> > > > > Breakpoint 1 at 0x4013c9: file test.f, line 3.
> > > > > (gdb) run
> > > > > Starting program: c:\usr\Tatsu\mingwhome\octaves\dgammatest/test.exe
> > > > > [New Thread 1048.0x8fc]
> > > > > 
> > > > > Breakpoint 1, test () at test.f:3
> > > > > 3             X=0
> > > > > Current language:  auto
> > > > > The current source language is "auto; currently fortran".
> > > > > (gdb) step
> > > > > 4             call DLGAMS (X, DLGAM, SGNGAM)
> > > > > (gdb) step
> > > > >  ***MESSAGE FROM ROUTINE DGAMMA IN LIBRARY SLATEC.
> > > > >  ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
> > > > >  *  X IS 0
> > > > >  *  ERROR NUMBER = 4
> > > > >  *
> > > > >  ***END OF MESSAGE
> > > > > 
> > > > >  ***JOB ABORT DUE TO FATAL ERROR.
> > > > > 0          ERROR MESSAGE SUMMARY
> > > > >  LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     
> > > > > COUNT
> > > > >  SLATEC     DGAMMA     X IS 0                       4         2       
> > > > >   1
> > > > > 
> > > > > gdb: unknown target exception 0xc0000027 at 0x77be5464
> > > > > 
> > > > > Program received signal ?, Unknown signal.
> > > > > 0x77be5464 in msvcrt!_global_unwind2 () from 
> > > > > C:\WINDOWS\system32\msvcrt.dll
> > > > > (gdb) bt
> > > > > #0  0x77be5464 in msvcrt!_global_unwind2 () from 
> > > > > C:\WINDOWS\system32\msvcrt.dll
> > > > > #1  0x77be6d8c in msvcrt!longjmp () from 
> > > > > C:\WINDOWS\system32\msvcrt.dll
> > > > > #2  0x00000000 in ?? ()
> > > > > 
> > > > > I have obtained the same error gcc-4.4.0 / MinGW
> > > > > 
> > > > > Unfortunately I did not install gdb on cygwin but I could the 
> > > > > following message in
> cygwin.
> > 
> > > > > 
> > > > >  ***MESSAGE FROM ROUTINE DGAMMA IN LIBRARY SLATEC.
> > > > >  ***FATAL ERROR, PROG ABORTED, TRACEBACK REQUESTED
> > > > >  *  X IS 0
> > > > >  *  ERROR NUMBER = 4
> > > > >  *
> > > > >  ***END OF MESSAGE
> > > > > 
> > > > >  ***JOB ABORT DUE TO FATAL ERROR.
> > > > > 0          ERROR MESSAGE SUMMARY
> > > > >  LIBRARY    SUBROUTINE MESSAGE START             NERR     LEVEL     
> > > > > COUNT
> > > > >  SLATEC     DGAMMA     X IS 0                       4         2       
> > > > >   1
> > > > > 
> > > > > ***********************************
> > > > > 
> > > > > The error message claim that X is 0 in DGAMMA.
> > > > > 
> > > > > However
> > > > > 
> > > > >       program test
> > > > >       DOUBLE PRECISION X, DLGAM, SGNGAM, DLNGAM
> > > > >       X=0
> > > > >       write(*,*) DGAMMA(X)
> > > > >       end
> > > > > 
> > > > > The above code gives
> > > > > 
> > > > > $ ./test
> > > > >                  +Infinity
> > > > > 
> > > > > ???????
> > > > > 
> > > > > Any suggestions?
> > > > > 
> > > > > Regards
> > > > > 
> > > > > Tatsuro
> > > > > 
> > > > > 
> > > > > 
> > > > > 
> > > > > --------------------------------------
> > > > > Get the new Internet Explorer 8 optimized for Yahoo! JAPAN
> > > > > http://pr.mail.yahoo.co.jp/ie8/
> > > > > 
> > > > 
> > > > 
> > > > --------------------------------------
> > > > Get the new Internet Explorer 8 optimized for Yahoo! JAPAN
> > > > http://pr.mail.yahoo.co.jp/ie8/
> > > > 
> > > 
> > > 
> > > --------------------------------------
> > > Get the new Internet Explorer 8 optimized for Yahoo! JAPAN
> > > http://pr.mail.yahoo.co.jp/ie8/
> > > 
> > 
> > 
> > --------------------------------------
> > Get the new Internet Explorer 8 optimized for Yahoo! JAPAN
> > http://pr.mail.yahoo.co.jp/ie8/
> > 
> 
> 
> --------------------------------------
> Get the new Internet Explorer 8 optimized for Yahoo! JAPAN
> http://pr.mail.yahoo.co.jp/ie8/
> 

--------------------------------------
Are you OK?  Online Safety Special Site - Yahoo! JAPAN
http://pr.mail.yahoo.co.jp/security/
slatec_gamma.diff
Description: 601022998-slatec_gamma.diff
[Prev in Thread]
Current Thread
[Next in Thread]
Re: gammaln(0) (octave-3.3.51+, MinGW), Tatsuro MATSUOKA, 2010/08/01
- Re: gammaln(0) (octave-3.3.51+, MinGW), Tatsuro MATSUOKA, 2010/08/01
  - Re: gammaln(0) (octave-3.3.51+, MinGW), Tatsuro MATSUOKA, 2010/08/01
    - Re: gammaln(0) (octave-3.3.51+, MinGW), Tatsuro MATSUOKA, 2010/08/01
    - Proposal changeset of liboctave/lo-specfun.cc (was Re: gammaln(0) (octave-3.3.51+, MinGW)), Tatsuro MATSUOKA <=
Prev by Date: Re: Adding functions to octave base?
Next by Date: help causes warning: delete: no such file: \octave-help-?????
Previous by thread: Re: gammaln(0) (octave-3.3.51+, MinGW)
Next by thread: Fwd: [bug #30616] isosurface fails with 2 input arguments (wrong meshgrid creation)
Index(es):
- Date
- Thread