* [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