[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: gammaln(0) (octave-3.3.51+, MinGW)
|
From: |
Tatsuro MATSUOKA |
|
Subject: |
Re: gammaln(0) (octave-3.3.51+, MinGW) |
|
Date: |
Mon, 2 Aug 2010 10:25:19 +0900 (JST) |
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/