public inbox for gcc-bugs@sourceware.org
help / color / mirror / Atom feed
* [Bug fortran/33698]  New: FAIL: gfortran.dg/gamma_5.f90
@ 2007-10-08 16:51 danglin at gcc dot gnu dot org
  2007-10-09  1:32 ` [Bug fortran/33698] " jvdelisle at gcc dot gnu dot org
                   ` (21 more replies)
  0 siblings, 22 replies; 23+ messages in thread
From: danglin at gcc dot gnu dot org @ 2007-10-08 16:51 UTC (permalink / raw)
  To: gcc-bugs

Executing on host: /test/gnu/gcc/objdir/gcc/testsuite/gfortran/../../gfortran
-B
/test/gnu/gcc/objdir/gcc/testsuite/gfortran/../../
/test/gnu/gcc/gcc/gcc/testsui
te/gfortran.dg/gamma_5.f90   -O0   -pedantic-errors 
-L/test/gnu/gcc/objdir/hppa
64-hp-hpux11.11/./libgfortran/.libs
-L/test/gnu/gcc/objdir/hppa64-hp-hpux11.11/.
/libgfortran/.libs -L/test/gnu/gcc/objdir/hppa64-hp-hpux11.11/./libiberty  -lm
 -o ./gamma_5.exe    (timeout = 300)
ld: Unsatisfied symbol "tgammaf" in file /var/tmp//ccUIJdfk.o
ld: Unsatisfied symbol "tgamma" in file /var/tmp//ccUIJdfk.o
2 errors.
collect2: ld returned 1 exit status
compiler exited with status 1
output is:
ld: Unsatisfied symbol "tgammaf" in file /var/tmp//ccUIJdfk.o
ld: Unsatisfied symbol "tgamma" in file /var/tmp//ccUIJdfk.o
2 errors.
collect2: ld returned 1 exit status

FAIL: gfortran.dg/gamma_5.f90  -O0  (test for excess errors)
Excess errors:
ld: Unsatisfied symbol "tgammaf" in file /var/tmp//ccUIJdfk.o
ld: Unsatisfied symbol "tgamma" in file /var/tmp//ccUIJdfk.o


-- 
           Summary: FAIL: gfortran.dg/gamma_5.f90
           Product: gcc
           Version: 4.3.0
            Status: UNCONFIRMED
          Severity: normal
          Priority: P3
         Component: fortran
        AssignedTo: unassigned at gcc dot gnu dot org
        ReportedBy: danglin at gcc dot gnu dot org
 GCC build triplet: hppa64-hp-hpux11.11
  GCC host triplet: hppa64-hp-hpux11.11
GCC target triplet: hppa64-hp-hpux11.11


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
@ 2007-10-09  1:32 ` jvdelisle at gcc dot gnu dot org
  2007-10-09  1:45 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (20 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: jvdelisle at gcc dot gnu dot org @ 2007-10-09  1:32 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #1 from jvdelisle at gcc dot gnu dot org  2007-10-09 01:32 -------
Just before the #include in trans-intrinsic.c we could do:

#ifndef tgamma
#define tgamma gamma
#endif

Could you try this and see it works?


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
  2007-10-09  1:32 ` [Bug fortran/33698] " jvdelisle at gcc dot gnu dot org
@ 2007-10-09  1:45 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-10-09  1:52 ` jvdelisle at gcc dot gnu dot org
                   ` (19 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-10-09  1:45 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #2 from dave at hiauly1 dot hia dot nrc dot ca  2007-10-09 01:45 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> ------- Comment #1 from jvdelisle at gcc dot gnu dot org  2007-10-09 01:32 -------
> Just before the #include in trans-intrinsic.c we could do:

Which include?

> #ifndef tgamma
> #define tgamma gamma
> #endif

Don't understand.

> Could you try this and see it works?

Sure.

Dave


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
  2007-10-09  1:32 ` [Bug fortran/33698] " jvdelisle at gcc dot gnu dot org
  2007-10-09  1:45 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-10-09  1:52 ` jvdelisle at gcc dot gnu dot org
  2007-10-09  2:12 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (18 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: jvdelisle at gcc dot gnu dot org @ 2007-10-09  1:52 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #3 from jvdelisle at gcc dot gnu dot org  2007-10-09 01:52 -------
Index: trans-intrinsic.c
===================================================================
--- trans-intrinsic.c   (revision 129058)
+++ trans-intrinsic.c   (working copy)
@@ -119,6 +119,10 @@ gfc_intrinsic_map_t;
 static GTY(()) gfc_intrinsic_map_t gfc_intrinsic_map[] =
 {
   /* Functions built into gcc itself.  */
+#ifndef tgamma
+#define tgamma gamma
+#endif
+
 #include "mathbuiltins.def"

   /* Functions in libm.  */


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (2 preceding siblings ...)
  2007-10-09  1:52 ` jvdelisle at gcc dot gnu dot org
@ 2007-10-09  2:12 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-10-09  2:19 ` jvdelisle at gcc dot gnu dot org
                   ` (17 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-10-09  2:12 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #4 from dave at hiauly1 dot hia dot nrc dot ca  2007-10-09 02:11 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

>    /* Functions built into gcc itself.  */
> +#ifndef tgamma
> +#define tgamma gamma
> +#endif
> +
>  #include "mathbuiltins.def"

The HP-UX 11 manpage says:

      lgamma() and gamma() return ln(|Gamma(x)|), where Gamma(x) is defined
      as the integral, as t goes from zero to infinity, of exp(-t) times t
      to the power (x-1).

      The sign of Gamma(x) is returned in the external integer signgam.

      The following C program fragment can be used to calculate
      Gamma(x):

          if ((y = lgamma(x)) > LN_MAXDOUBLE)
            error();
          y = signgam * exp(y);

       where if y is greater than LN_MAXDOUBLE, as defined in the <values.h>
       header file, exp() returns a range error (see exp(3M)).

There's also a reentrant version lgamma_r:

       double lgamma_r(double x, int *sign);

The main point is the HP-UX gamma function isn't the true gamma function.

Dave


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (3 preceding siblings ...)
  2007-10-09  2:12 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-10-09  2:19 ` jvdelisle at gcc dot gnu dot org
  2007-10-09  2:39 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (16 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: jvdelisle at gcc dot gnu dot org @ 2007-10-09  2:19 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #5 from jvdelisle at gcc dot gnu dot org  2007-10-09 02:18 -------
I was just looking at a DOC search and for HP-UX 11i Version 3: February 2007

tgamma is suppose to be provided.

So I now wonder what is really going on.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (4 preceding siblings ...)
  2007-10-09  2:19 ` jvdelisle at gcc dot gnu dot org
@ 2007-10-09  2:39 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-10-09 20:22 ` tkoenig at gcc dot gnu dot org
                   ` (15 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-10-09  2:39 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #6 from dave at hiauly1 dot hia dot nrc dot ca  2007-10-09 02:39 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> ------- Comment #5 from jvdelisle at gcc dot gnu dot org  2007-10-09 02:18 -------
> I was just looking at a DOC search and for HP-UX 11i Version 3: February 2007
> 
> tgamma is suppose to be provided.
> 
> So I now wonder what is really going on.

For parisc, 11i Version 3 is only supported on recent servers (mainly the
big ones).  No workstations are supported.  11i Version 3 is supposed to
be certified UNIX 2003 compliant on ia64.  However, this isn't the case
on parisc.  I think that is mainly because HP doesn't provide a compliant
compiler.  It is interesting that 11i V3 appears to have complete C99
function support even on parisc, although I wouldn't be surprised if
long double support was missing.  Unfortunately, older HP-UX versions
don't have the tgamma function. 

Dave


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (5 preceding siblings ...)
  2007-10-09  2:39 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-10-09 20:22 ` tkoenig at gcc dot gnu dot org
  2007-10-09 22:06 ` burnus at gcc dot gnu dot org
                   ` (14 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at gcc dot gnu dot org @ 2007-10-09 20:22 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #7 from tkoenig at gcc dot gnu dot org  2007-10-09 20:22 -------
We'll probably need to roll our own tgamma function:  To
cover cases like this, where the system doesn't provide
one, and to get numerically better answers.


-- 

tkoenig at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |tkoenig at gcc dot gnu dot
                   |                            |org
             Status|UNCONFIRMED                 |NEW
     Ever Confirmed|0                           |1
   Last reconfirmed|0000-00-00 00:00:00         |2007-10-09 20:22:06
               date|                            |


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (6 preceding siblings ...)
  2007-10-09 20:22 ` tkoenig at gcc dot gnu dot org
@ 2007-10-09 22:06 ` burnus at gcc dot gnu dot org
  2007-10-20 10:06 ` tkoenig at gcc dot gnu dot org
                   ` (13 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: burnus at gcc dot gnu dot org @ 2007-10-09 22:06 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #8 from burnus at gcc dot gnu dot org  2007-10-09 22:06 -------
> We'll probably need to roll our own tgamma function:  To
> cover cases like this, where the system doesn't provide
> one, and to get numerically better answers.

g95 uses a C version of W. J. Cody and L. Stoltz' Fortran procedure available
from http://www.netlib.org/specfun/gamma

If I read GLIBC correctly, it implements gamma as exp(lgamma(x)), which
explains why the results are not optimal, cf.

http://sources.redhat.com/cgi-bin/cvsweb.cgi/libc/math/w_tgamma.c?cvsroot=glibc
http://sources.redhat.com/cgi-bin/cvsweb.cgi/libc/sysdeps/ieee754/dbl-64/e_gamma_r.c?cvsroot=glibc

Initially, I had claimed that it is better to use the C99/POSIX-2001 tgamma
function of the system C library (if available), but now I am not so sure
anymore.


-- 

burnus at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
                 CC|                            |burnus at gcc dot gnu dot
                   |                            |org


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (7 preceding siblings ...)
  2007-10-09 22:06 ` burnus at gcc dot gnu dot org
@ 2007-10-20 10:06 ` tkoenig at gcc dot gnu dot org
  2007-10-20 10:17 ` tkoenig at gcc dot gnu dot org
                   ` (12 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at gcc dot gnu dot org @ 2007-10-20 10:06 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #9 from tkoenig at gcc dot gnu dot org  2007-10-20 10:06 -------
Mine.


-- 

tkoenig at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
         AssignedTo|unassigned at gcc dot gnu   |tkoenig at gcc dot gnu dot
                   |dot org                     |org
             Status|NEW                         |ASSIGNED
   Last reconfirmed|2007-10-09 20:22:06         |2007-10-20 10:06:20
               date|                            |


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (8 preceding siblings ...)
  2007-10-20 10:06 ` tkoenig at gcc dot gnu dot org
@ 2007-10-20 10:17 ` tkoenig at gcc dot gnu dot org
  2007-11-02 23:07 ` tkoenig at gcc dot gnu dot org
                   ` (11 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at gcc dot gnu dot org @ 2007-10-20 10:17 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #10 from tkoenig at gcc dot gnu dot org  2007-10-20 10:17 -------
I am tempted not to use gamma*() at all, but rather
check for tgamma*() and lgamma*() and use that
if available.  If none of them are present, use
a fallback implementation.

This avoids potential problems with cross-compilation
when gamma() has a different meaning for host and
target system.

What a mess...


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (9 preceding siblings ...)
  2007-10-20 10:17 ` tkoenig at gcc dot gnu dot org
@ 2007-11-02 23:07 ` tkoenig at gcc dot gnu dot org
  2007-11-02 23:45 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (10 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at gcc dot gnu dot org @ 2007-11-02 23:07 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #11 from tkoenig at gcc dot gnu dot org  2007-11-02 23:07 -------
Created an attachment (id=14475)
 --> (http://gcc.gnu.org/bugzilla/attachment.cgi?id=14475&action=view)
proposed patch

This implements the fallback functions, but naturally
doesn't do anything on my linux system (which has all tgamma*
and lgamma* functions).

Can anybody test this?  John?


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (10 preceding siblings ...)
  2007-11-02 23:07 ` tkoenig at gcc dot gnu dot org
@ 2007-11-02 23:45 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-11-03 19:37 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (9 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-11-02 23:45 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #12 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-02 23:45 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> Can anybody test this?  John?

I'm on it.

Dave


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (11 preceding siblings ...)
  2007-11-02 23:45 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-11-03 19:37 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-11-03 20:19 ` tkoenig at netcologne dot de
                   ` (8 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-11-03 19:37 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #13 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-03 19:37 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> This implements the fallback functions, but naturally
> doesn't do anything on my linux system (which has all tgamma*
> and lgamma* functions).

You might hack config.h to test.

> Can anybody test this?  John?

c99_functions.c didn't compile due to some typos, etc.  I've attached
ammended version of the patch which builds on hppa2.0w-hp-hpux11.11.

Now, the hard part.  The gfortran.dg/gamma_5.f90 test fails at n = 16.

Breakpoint 2, tgamma (z=16.5)
    at ../../../gcc/libgfortran/intrinsics/c99_functions.c:1730
1730      if (z < 0.0)
(gdb)
Continuing.

Breakpoint 1, _gfortran_abort ()
    at ../../../gcc/libgfortran/intrinsics/abort.c:38
38        close_units ();

It fails here:
     if (abs(gamma(xd)-td)/td > 5e-14) call abort

I hacked the test a bit to make it easier to see the return value
from gamma(xd)-td):

(gdb) p gd
$1 = 5189998453040.5889
(gdb) p td
$2 = 5189998453040.125
(gdb) p (gd-td)/td
$4 = 8.9377134058350696e-14

I guess the first thing to check is the coefficients.

Dave


------- Comment #14 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-03 19:37 -------
Created an attachment (id=14479)
 --> (http://gcc.gnu.org/bugzilla/attachment.cgi?id=14479&action=view)


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (12 preceding siblings ...)
  2007-11-03 19:37 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-11-03 20:19 ` tkoenig at netcologne dot de
  2007-11-03 20:47 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (7 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at netcologne dot de @ 2007-11-03 20:19 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #15 from tkoenig at netcologne dot de  2007-11-03 20:19 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90


> ------- Comment #13 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-03 19:37 -------

> Now, the hard part.  The gfortran.dg/gamma_5.f90 test fails at n = 16.

This works on my i686-pc-linux-gnu system, and also fails when I use
-ffloat-store.  Seems like we have a roundoff problem with normal ieee
double precision.

I'll continue to work on it...


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (13 preceding siblings ...)
  2007-11-03 20:19 ` tkoenig at netcologne dot de
@ 2007-11-03 20:47 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-11-03 22:49 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (6 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-11-03 20:47 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #16 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-03 20:47 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> This works on my i686-pc-linux-gnu system, and also fails when I use
> -ffloat-store.  Seems like we have a roundoff problem with normal ieee
> double precision.

Test passes if I increase 5e-14 to 17e-14.

Dave


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (14 preceding siblings ...)
  2007-11-03 20:47 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-11-03 22:49 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-11-04 19:52 ` tkoenig at gcc dot gnu dot org
                   ` (5 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-11-03 22:49 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #17 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-03 22:49 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> > Now, the hard part.  The gfortran.dg/gamma_5.f90 test fails at n = 16.
> 
> This works on my i686-pc-linux-gnu system, and also fails when I use
> -ffloat-store.  Seems like we have a roundoff problem with normal ieee
> double precision.
> 
> I'll continue to work on it...

Looking at Pugh's paper, he shows coefficients for n = 10 and
r = 10.900511.  This is the same as you are using for the double
case.  However, you seem to be missing coefficient d10.  I believe
that Pugh's coefficients and your coefficients should differ by
the factor

     2 * sqrt (e / pi) * exp (r)

due to Toth's rearrangement of the approximation.  This should
provide a check for your coefficient calculation.

Dave


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (15 preceding siblings ...)
  2007-11-03 22:49 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-11-04 19:52 ` tkoenig at gcc dot gnu dot org
  2007-11-04 21:30 ` dave at hiauly1 dot hia dot nrc dot ca
                   ` (4 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at gcc dot gnu dot org @ 2007-11-04 19:52 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #18 from tkoenig at gcc dot gnu dot org  2007-11-04 19:51 -------

> Looking at Pugh's paper, he shows coefficients for n = 10 and
> r = 10.900511.  This is the same as you are using for the double
> case.  However, you seem to be missing coefficient d10.

Good catch, thanks!

The main problem with the Lanczos approximation are alternating-sign
terms with nearly cancel each other, which leads to massive precision
loss.

For z=16.5 and r=10.900511, the terms in the sum are

1 6425.81075191694890236249
2 -19919.53511610527857556008
3 24595.63902224190678680316
4 -15425.21437829293790855445
5 5196.45802113903846475296
6 -916.60901318718765651283
7   76.62541745991659070114
8   -2.45417886377221794447
9    0.01907340042601936639
10   -0.00001080118945830201

and the sum is

 33.24621718250205049117

so I'd expect about log2(24000/33) ~ 9.5 bits of precision loss.
Not good.

Some rearrangement can help here, possibly.  OTOH, maybe using
straight Netlib code would be better.

Ouch.

Maybe it's better to fall back on the Netlib code.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (16 preceding siblings ...)
  2007-11-04 19:52 ` tkoenig at gcc dot gnu dot org
@ 2007-11-04 21:30 ` dave at hiauly1 dot hia dot nrc dot ca
  2007-11-09 18:41 ` fxcoudert at gcc dot gnu dot org
                   ` (3 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: dave at hiauly1 dot hia dot nrc dot ca @ 2007-11-04 21:30 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #19 from dave at hiauly1 dot hia dot nrc dot ca  2007-11-04 21:30 -------
Subject: Re:  FAIL: gfortran.dg/gamma_5.f90

> The main problem with the Lanczos approximation are alternating-sign
> terms with nearly cancel each other, which leads to massive precision
> loss.
> 
> For z=16.5 and r=10.900511, the terms in the sum are
> 
> 1 6425.81075191694890236249
> 2 -19919.53511610527857556008
> 3 24595.63902224190678680316
> 4 -15425.21437829293790855445
> 5 5196.45802113903846475296
> 6 -916.60901318718765651283
> 7   76.62541745991659070114
> 8   -2.45417886377221794447
> 9    0.01907340042601936639
> 10   -0.00001080118945830201
> 
> and the sum is
> 
>  33.24621718250205049117
> 
> so I'd expect about log2(24000/33) ~ 9.5 bits of precision loss.
> Not good.

I don't think this can be avoided without a huge performance hit
(i.e., use long double arithmetic for the sum).  Possibly, the
float versions would benefit from doing the sum with doubles.
That shouldn't hurt.

> Some rearrangement can help here, possibly.  OTOH, maybe using
> straight Netlib code would be better.

On systems with lgamma, it's possible to use exp (lgamma (x)) for
tgamma (x).  I checked and gamma_5 doesn't fail on HP-UX with this
implementation.  I retained the transformation for negative arguments.
Float variants could easily be implemented using lgamma.  I think
using the system lgamma implementation might be best when it's available.

Found an old fortran implementation on one of my machines!  I believe
that I wrote an implementation using a rational approximation (Pade?)
about 35 years ago, but I couldn't find it.

Dave
-- 
J. David Anglin                                  dave.anglin@nrc-cnrc.gc.ca
National Research Council of Canada              (613) 990-0752 (FAX: 952-6602)

      double precision function dgamma(z)
c
c     Source for routine:  J.F. Hart et al, Computer Approximations,
c     Wiley,1968.
c     This is a double precision routine for the gamma function of
c     real positive arguments.  Values for negative arguments (non-
c     integral or zero) may be obtained from this program and
c     standard functional relations for the gamma function.
c     argument range           method of computation
c     z>14                  gamma(z)=dexp(ln(gamma(z))
c     3<= z <=14        argument reduction z*gam(z)=gam(z+1)
c     2<= z <3          rational approx. w/o argument reduction
c     0<  z <2          argument expansion z*gam(z)=gam(z+1)
c
      implicit double precision (a-h,o-z)
c
      double precision nlsqp
c
      dimension pnum(11),qden(11),phnum(4),qhden(4)
c
      data pnum,qden,phnum,qhden/
     1-.2983543278574342138830437659d+6,
     2-.2384953970018198872468734423d+6,
     3-.1170494760121780688403854445d+6,
     4-.3949445048301571936421824091d+5,
     5-.1046699423827521405330650531d+5,
     6-.2188218110071816359394795998d+4,
     7-.3805112208641734657584922631d+3,
     8-.5283123755635845383718978382d+2,
     9-.6128571763704498306889428212d+1,
     9-.502801805441681246736419875d00,
     9-.3343060322330595274515660112d-1,
     1-.2983543278574342138830438524d+6,
     2-.1123558608748644911342306408d+6,
     3+.5332716689118142157485686311d+5,
     4+.8571160498907043851961147763d+4,
     5-.473486597702821170655681977d+4,
     6+.1960497612885585838997039621d+3,
     7+.1257733367869888645966647426d+3,
     8-.2053126153100672764513929067d+2,
     9+.1d1,
     90.d0,
     90.d0,
     1+.12398282342474941538685913d0,
     2+.670827838343321349614617d0,
     3+.6450730291289920251389d0,
     4+.666629070402007526d-1,
     1+.148779388109699298468156d+1,
     2+.809952718948975574728214d+1,
     3+.7996691123663644194772d+1,
     4+.1d1/
c
      az=z
      pi=4.d0*datan(1.d0)
      accum=1.d0
      if (az.gt.14.d0) go to 44
      if (az.lt.2.d0) go to 333
      if (az.lt.3.d0) go to 33
      do 22 i=1,12
      az=az-1.d0
      accum=accum*az
      if (az.lt.3.d0) go to 33
   22 continue
  333 do 222 i=1,12
      accum=accum*(1.d0/az)
      az=az+1.d0
      if (az.ge.2.d0) go to 33
  222 continue
   33 az=az-2.d0
      anum=0.d0
      aden=0.d0
      do 11 i=0,10
      ii=i+1
      if (i.eq.0) power=1.d0
      if (i.eq.0) go to 55
      power=az**i
   55 anum=anum+pnum(ii)*power
      aden=aden+qden(ii)*power
   11 continue
      dgamma=accum*(anum/aden)
      go to 999
   44 anum=0.d0
      aden=0.d0
      do 1 i=0,3
      ii=i+1
      if (i.eq.0) power=1.d0
      if (i.eq.0) go to 5
      power=(1.d0/(az*az))**i
    5 anum=anum+phnum(ii)*power
      aden=aden+qhden(ii)*power
    1 continue
      riemn=anum/(aden*az)
      nlsqp=dlog(dsqrt(2.d0*pi))
      daleg=(az-0.5d0)*dlog(az)-az+nlsqp+riemn
      dgamma=dexp(daleg)
  999 return
      end


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (17 preceding siblings ...)
  2007-11-04 21:30 ` dave at hiauly1 dot hia dot nrc dot ca
@ 2007-11-09 18:41 ` fxcoudert at gcc dot gnu dot org
  2007-11-09 21:51 ` tkoenig at gcc dot gnu dot org
                   ` (2 subsequent siblings)
  21 siblings, 0 replies; 23+ messages in thread
From: fxcoudert at gcc dot gnu dot org @ 2007-11-09 18:41 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #20 from fxcoudert at gcc dot gnu dot org  2007-11-09 18:40 -------
Created an attachment (id=14519)
 --> (http://gcc.gnu.org/bugzilla/attachment.cgi?id=14519&action=view)
gamma from netlib, converted to C

(In reply to comment #18)
> OTOH, maybe using straight Netlib code would be better.

Thomas, here is my C version of the netlib gamma, including some testing code.
I also think using the netlib version is rather safe and we shouldn't take any
risk for a fallback implementation of a math function.

I can work on the C version of lgamma at some point in the week-end, if you
want me to.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug fortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (18 preceding siblings ...)
  2007-11-09 18:41 ` fxcoudert at gcc dot gnu dot org
@ 2007-11-09 21:51 ` tkoenig at gcc dot gnu dot org
  2007-11-16 22:32 ` [Bug libfortran/33698] " fxcoudert at gcc dot gnu dot org
  2007-11-16 22:46 ` fxcoudert at gcc dot gnu dot org
  21 siblings, 0 replies; 23+ messages in thread
From: tkoenig at gcc dot gnu dot org @ 2007-11-09 21:51 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #21 from tkoenig at gcc dot gnu dot org  2007-11-09 21:50 -------
(In reply to comment #20)

Hi FX,

> I can work on the C version of lgamma at some point in the week-end, if you
> want me to.

That would be great.  Real Life has caught up with me in a big
way the last few weeks, so I couldn't do anything in this direction.


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug libfortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (19 preceding siblings ...)
  2007-11-09 21:51 ` tkoenig at gcc dot gnu dot org
@ 2007-11-16 22:32 ` fxcoudert at gcc dot gnu dot org
  2007-11-16 22:46 ` fxcoudert at gcc dot gnu dot org
  21 siblings, 0 replies; 23+ messages in thread
From: fxcoudert at gcc dot gnu dot org @ 2007-11-16 22:32 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #22 from fxcoudert at gcc dot gnu dot org  2007-11-16 22:31 -------
Subject: Bug 33698

Author: fxcoudert
Date: Fri Nov 16 22:31:28 2007
New Revision: 130245

URL: http://gcc.gnu.org/viewcvs?root=gcc&view=rev&rev=130245
Log:
        PR libfortran/33583
        PR libfortran/33698

        * intrinsics/c99_functions.c (tgamma, tgammaf, lgamma, lgammaf):
        New fallback functions.
        * c99_protos.h (tgamma, tgammaf, lgamma, lgammaf): New prototypes.
        * configure.ac: Add checks for tgamma, tgammaf, tgammal, lgamma,
        lgammaf and lgammal.
        * config.h.in: Regenerate.
        * configure: Regenerate.

Modified:
    trunk/libgfortran/ChangeLog
    trunk/libgfortran/c99_protos.h
    trunk/libgfortran/config.h.in
    trunk/libgfortran/configure
    trunk/libgfortran/configure.ac
    trunk/libgfortran/intrinsics/c99_functions.c


-- 


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

* [Bug libfortran/33698] FAIL: gfortran.dg/gamma_5.f90
  2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
                   ` (20 preceding siblings ...)
  2007-11-16 22:32 ` [Bug libfortran/33698] " fxcoudert at gcc dot gnu dot org
@ 2007-11-16 22:46 ` fxcoudert at gcc dot gnu dot org
  21 siblings, 0 replies; 23+ messages in thread
From: fxcoudert at gcc dot gnu dot org @ 2007-11-16 22:46 UTC (permalink / raw)
  To: gcc-bugs



------- Comment #23 from fxcoudert at gcc dot gnu dot org  2007-11-16 22:46 -------
I hope it's now fixed. If there is something more to it, please reopen the PR
or file a new one.


-- 

fxcoudert at gcc dot gnu dot org changed:

           What    |Removed                     |Added
----------------------------------------------------------------------------
             Status|ASSIGNED                    |RESOLVED
         Resolution|                            |FIXED
   Target Milestone|---                         |4.3.0


http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33698


^ permalink raw reply	[flat|nested] 23+ messages in thread

end of thread, other threads:[~2007-11-16 22:46 UTC | newest]

Thread overview: 23+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2007-10-08 16:51 [Bug fortran/33698] New: FAIL: gfortran.dg/gamma_5.f90 danglin at gcc dot gnu dot org
2007-10-09  1:32 ` [Bug fortran/33698] " jvdelisle at gcc dot gnu dot org
2007-10-09  1:45 ` dave at hiauly1 dot hia dot nrc dot ca
2007-10-09  1:52 ` jvdelisle at gcc dot gnu dot org
2007-10-09  2:12 ` dave at hiauly1 dot hia dot nrc dot ca
2007-10-09  2:19 ` jvdelisle at gcc dot gnu dot org
2007-10-09  2:39 ` dave at hiauly1 dot hia dot nrc dot ca
2007-10-09 20:22 ` tkoenig at gcc dot gnu dot org
2007-10-09 22:06 ` burnus at gcc dot gnu dot org
2007-10-20 10:06 ` tkoenig at gcc dot gnu dot org
2007-10-20 10:17 ` tkoenig at gcc dot gnu dot org
2007-11-02 23:07 ` tkoenig at gcc dot gnu dot org
2007-11-02 23:45 ` dave at hiauly1 dot hia dot nrc dot ca
2007-11-03 19:37 ` dave at hiauly1 dot hia dot nrc dot ca
2007-11-03 20:19 ` tkoenig at netcologne dot de
2007-11-03 20:47 ` dave at hiauly1 dot hia dot nrc dot ca
2007-11-03 22:49 ` dave at hiauly1 dot hia dot nrc dot ca
2007-11-04 19:52 ` tkoenig at gcc dot gnu dot org
2007-11-04 21:30 ` dave at hiauly1 dot hia dot nrc dot ca
2007-11-09 18:41 ` fxcoudert at gcc dot gnu dot org
2007-11-09 21:51 ` tkoenig at gcc dot gnu dot org
2007-11-16 22:32 ` [Bug libfortran/33698] " fxcoudert at gcc dot gnu dot org
2007-11-16 22:46 ` fxcoudert at gcc dot gnu dot org

This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).