This is the mail archive of the
gcc-patches@gcc.gnu.org
mailing list for the GCC project.
[Patch, Fortran] PR32980 - Add GAMMA and LGAMMA intrinsics
- From: Tobias Burnus <burnus at net-b dot de>
- To: gcc-patches <gcc-patches at gcc dot gnu dot org>, "'fortran at gcc dot gnu dot org'" <fortran at gcc dot gnu dot org>
- Date: Tue, 21 Aug 2007 21:00:18 +0200
- Subject: [Patch, Fortran] PR32980 - Add GAMMA and LGAMMA intrinsics
:ADDPATCH fortran:
This patch adds:
a) The gamma function as vendor intrinsic GAMMA with the specific names
GAMMA and DGAMMA. (Only) "GAMMA" is also part of the current Fortran
2008 draft.
b) The logarithm of the absolute value of the Gamma function as vendor
intrinsic LGAMMA with the specific names AGAMA and DGAMA. The Fortran
2008 draft has this function under the name LOG_GAMMA.
While LGAMMA is compatible to LOG_GAMMA, I did not add the F2008/draft
version. Anyone thinks that I should?
Note: The functions tgamma (t = true) and lgamma are part of C99.
Regression tested on x86-64/Linux. Ok for the trunk?
Tobias
PS: GAMMA seems to be the last missing popular vendor intrinsic -
ignoring things which should be rather implemented in modules such as
functions to interface with POSIX or similar.
PPS: I think BUILTIN_GAMMA and BUILTIN_TGAMMA are identically, if
someone has reason to believe that this is not the case, please tell me.
2007-08-21 Tobias Burnus <burnus@net-b.de>
PR fortran/32980
* intrinsic.h (gfc_simplify_gamma,gfc_simplify_lgamma,
gfc_resolve_gamma,gfc_resolve_lgamma): New function declations.
* mathbuiltins.def: Define GAMMA and LGAMMA.
* intrinsic.c (add_functions): Add GAMMA, DGAMMA, LGAMMA, ALGAMA
and DLGAMA.
* simplify.c (gfc_simplify_gamma,gfc_simplify_lgamma): New functions.
* iresolve.c (gfc_resolve_gamma,gfc_resolve_lgamma): New functions.
* intrinsic.texi: Add documentation for GAMMA and LGAMMA.
2007-08-21 Tobias Burnus <burnus@net-b.de>
PR fortran/32980
* gfortran.dg/gamma_1.f90: New.
* gfortran.dg/gamma_2.f90: New.
* gfortran.dg/gamma_3.f90: New.
Index: gcc/fortran/intrinsic.h
===================================================================
--- gcc/fortran/intrinsic.h (revision 127666)
+++ gcc/fortran/intrinsic.h (working copy)
@@ -221,6 +221,7 @@ gfc_expr *gfc_simplify_exponent (gfc_exp
gfc_expr *gfc_simplify_float (gfc_expr *);
gfc_expr *gfc_simplify_floor (gfc_expr *, gfc_expr *);
gfc_expr *gfc_simplify_fraction (gfc_expr *);
+gfc_expr *gfc_simplify_gamma (gfc_expr *);
gfc_expr *gfc_simplify_huge (gfc_expr *);
gfc_expr *gfc_simplify_iachar (gfc_expr *, gfc_expr *);
gfc_expr *gfc_simplify_iand (gfc_expr *, gfc_expr *);
@@ -243,6 +244,7 @@ gfc_expr *gfc_simplify_kind (gfc_expr *)
gfc_expr *gfc_simplify_lbound (gfc_expr *, gfc_expr *, gfc_expr *);
gfc_expr *gfc_simplify_len (gfc_expr *, gfc_expr *);
gfc_expr *gfc_simplify_len_trim (gfc_expr *, gfc_expr *);
+gfc_expr *gfc_simplify_lgamma (gfc_expr *);
gfc_expr *gfc_simplify_lge (gfc_expr *, gfc_expr *);
gfc_expr *gfc_simplify_lgt (gfc_expr *, gfc_expr *);
gfc_expr *gfc_simplify_lle (gfc_expr *, gfc_expr *);
@@ -354,6 +356,7 @@ void gfc_resolve_fget (gfc_expr *, gfc_e
void gfc_resolve_fputc (gfc_expr *, gfc_expr *, gfc_expr *);
void gfc_resolve_fput (gfc_expr *, gfc_expr *);
void gfc_resolve_g77_math1 (gfc_expr *, gfc_expr *);
+void gfc_resolve_gamma (gfc_expr *, gfc_expr *);
void gfc_resolve_getcwd (gfc_expr *, gfc_expr *);
void gfc_resolve_getgid (gfc_expr *);
void gfc_resolve_getpid (gfc_expr *);
@@ -384,6 +387,7 @@ void gfc_resolve_kill (gfc_expr *, gfc_e
void gfc_resolve_lbound (gfc_expr *, gfc_expr *, gfc_expr *, gfc_expr *);
void gfc_resolve_len (gfc_expr *, gfc_expr *, gfc_expr *);
void gfc_resolve_len_trim (gfc_expr *, gfc_expr *, gfc_expr *);
+void gfc_resolve_lgamma (gfc_expr *, gfc_expr *);
void gfc_resolve_link (gfc_expr *, gfc_expr *, gfc_expr *);
void gfc_resolve_loc (gfc_expr *, gfc_expr *);
void gfc_resolve_log (gfc_expr *, gfc_expr *);
Index: gcc/fortran/mathbuiltins.def
===================================================================
--- gcc/fortran/mathbuiltins.def (revision 127666)
+++ gcc/fortran/mathbuiltins.def (working copy)
@@ -30,3 +30,5 @@ DEFINE_MATH_BUILTIN (Y1, "y1",
DEFINE_MATH_BUILTIN (YN, "yn", 2)
DEFINE_MATH_BUILTIN (ERF, "erf", 0)
DEFINE_MATH_BUILTIN (ERFC, "erfc", 0)
+DEFINE_MATH_BUILTIN (GAMMA, "gamma", 0)
+DEFINE_MATH_BUILTIN (LGAMMA,"lgamma", 0)
Index: gcc/fortran/intrinsic.c
===================================================================
--- gcc/fortran/intrinsic.c (revision 127666)
+++ gcc/fortran/intrinsic.c (working copy)
@@ -1453,6 +1453,16 @@ add_functions (void)
make_generic ("fput", GFC_ISYM_FPUT, GFC_STD_GNU);
+ add_sym_1 ("gamma", GFC_ISYM_GAMMA, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
+ gfc_check_fn_r, gfc_simplify_gamma, gfc_resolve_gamma,
+ x, BT_REAL, dr, REQUIRED);
+
+ add_sym_1 ("dgamma", GFC_ISYM_GAMMA, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
+ gfc_check_fn_r, gfc_simplify_gamma, gfc_resolve_gamma,
+ x, BT_REAL, dr, REQUIRED);
+
+ make_generic ("gamma", GFC_ISYM_GAMMA, GFC_STD_GNU);
+
/* Unix IDs (g77 compatibility) */
add_sym_1 ("getcwd", GFC_ISYM_GETCWD, NO_CLASS, ACTUAL_NO, BT_INTEGER, di, GFC_STD_GNU,
NULL, NULL, gfc_resolve_getcwd,
@@ -1690,6 +1700,21 @@ add_functions (void)
make_generic ("len_trim", GFC_ISYM_LEN_TRIM, GFC_STD_F95);
+ add_sym_1 ("lgamma", GFC_ISYM_LGAMMA, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
+ gfc_check_fn_r, gfc_simplify_lgamma, gfc_resolve_lgamma,
+ x, BT_REAL, dr, REQUIRED);
+
+ add_sym_1 ("algama", GFC_ISYM_LGAMMA, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
+ gfc_check_fn_r, gfc_simplify_lgamma, gfc_resolve_lgamma,
+ x, BT_REAL, dr, REQUIRED);
+
+ add_sym_1 ("dlgama", GFC_ISYM_LGAMMA, CLASS_ELEMENTAL, ACTUAL_NO, BT_REAL, dr, GFC_STD_GNU,
+ gfc_check_fn_r, gfc_simplify_lgamma, gfc_resolve_lgamma,
+ x, BT_REAL, dr, REQUIRED);
+
+ make_generic ("lgamma", GFC_ISYM_LGAMMA, GFC_STD_GNU);
+
+
add_sym_2 ("lge", GFC_ISYM_LGE, CLASS_ELEMENTAL, ACTUAL_NO, BT_LOGICAL, dl, GFC_STD_F77,
NULL, gfc_simplify_lge, NULL,
sta, BT_CHARACTER, dc, REQUIRED, stb, BT_CHARACTER, dc, REQUIRED);
Index: gcc/fortran/simplify.c
===================================================================
--- gcc/fortran/simplify.c (revision 127666)
+++ gcc/fortran/simplify.c (working copy)
@@ -1183,6 +1183,24 @@ gfc_simplify_fraction (gfc_expr *x)
gfc_expr *
+gfc_simplify_gamma (gfc_expr *x)
+{
+ gfc_expr *result;
+
+ if (x->expr_type != EXPR_CONSTANT)
+ return NULL;
+
+ result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+
+ gfc_set_model_kind (x->ts.kind);
+
+ mpfr_gamma (result->value.real, x->value.real, GFC_RND_MODE);
+
+ return range_check (result, "GAMMA");
+}
+
+
+gfc_expr *
gfc_simplify_huge (gfc_expr *e)
{
gfc_expr *result;
@@ -2212,6 +2230,27 @@ gfc_simplify_len_trim (gfc_expr *e, gfc_
return range_check (result, "LEN_TRIM");
}
+gfc_expr *
+gfc_simplify_lgamma (gfc_expr *x __attribute__((unused)))
+{
+#if MPFR_VERSION >= MPFR_VERSION_NUM(2,3,0)
+ gfc_expr *result;
+
+ if (x->expr_type != EXPR_CONSTANT)
+ return NULL;
+
+ result = gfc_constant_result (x->ts.type, x->ts.kind, &x->where);
+
+ gfc_set_model_kind (x->ts.kind);
+
+ mpfr_lgamma (result->value.real, x->value.real, GFC_RND_MODE);
+
+ return range_check (result, "LGAMMA");
+#else
+ return NULL;
+#endif
+}
+
gfc_expr *
gfc_simplify_lge (gfc_expr *a, gfc_expr *b)
Index: gcc/fortran/iresolve.c
===================================================================
--- gcc/fortran/iresolve.c (revision 127666)
+++ gcc/fortran/iresolve.c (working copy)
@@ -771,6 +771,15 @@ gfc_resolve_g77_math1 (gfc_expr *f, gfc_
void
+gfc_resolve_gamma (gfc_expr *f, gfc_expr *x)
+{
+ f->ts = x->ts;
+ f->value.function.name
+ = gfc_get_string ("__gamma_%d", x->ts.kind);
+}
+
+
+void
gfc_resolve_getcwd (gfc_expr *f, gfc_expr *n ATTRIBUTE_UNUSED)
{
f->ts.type = BT_INTEGER;
@@ -1128,6 +1137,15 @@ gfc_resolve_len_trim (gfc_expr *f, gfc_e
void
+gfc_resolve_lgamma (gfc_expr *f, gfc_expr *x)
+{
+ f->ts = x->ts;
+ f->value.function.name
+ = gfc_get_string ("__lgamma_%d", x->ts.kind);
+}
+
+
+void
gfc_resolve_link (gfc_expr *f, gfc_expr *p1 ATTRIBUTE_UNUSED,
gfc_expr *p2 ATTRIBUTE_UNUSED)
{
Index: gcc/fortran/intrinsic.texi
===================================================================
--- gcc/fortran/intrinsic.texi (revision 127666)
+++ gcc/fortran/intrinsic.texi (working copy)
@@ -121,6 +121,7 @@ Some basic guidelines for editing this d
* @code{FSEEK}: FSEEK, Low level file positioning subroutine
* @code{FSTAT}: FSTAT, Get file status
* @code{FTELL}: FTELL, Current stream position
+* @code{GAMMA}: GAMMA, Gamma function
* @code{GERROR}: GERROR, Get last system error message
* @code{GETARG}: GETARG, Get command line arguments
* @code{GET_COMMAND}: GET_COMMAND, Get the entire command line
@@ -161,6 +164,7 @@ Some basic guidelines for editing this d
* @code{LBOUND}: LBOUND, Lower dimension bounds of an array
* @code{LEN}: LEN, Length of a character entity
* @code{LEN_TRIM}: LEN_TRIM, Length of a character entity without trailing blank characters
+* @code{LGAMMA}: LGAMMA, Logarithm of the Gamma function
* @code{LGE}: LGE, Lexical greater than or equal
* @code{LGT}: LGT, Lexical greater than
* @code{LINK}: LINK, Create a hard link
@@ -4484,6 +4488,65 @@ END PROGRAM
+@node GAMMA
+@section @code{GAMMA} --- Gamma function
+@fnindex GAMMA
+@fnindex DGAMMA
+@cindex Gamma function
+@cindex Factorial function
+
+@table @asis
+@item @emph{Description}:
+@code{GAMMA(X)} computes Gamma (@math{\Gamma}) of @var{X}. For positive,
+integer values of @var{X} the Gamma function simplifies to the factorial
+function @math{\Gamma(x)=(x-1)!}.
+
+@tex
+$$
+\Gamma(x) = \int_0^\infty t^{x-1}{\rm e}^{-t}\,{\rm d}t
+$$
+@end tex
+
+@item @emph{Standard}:
+GNU Extension
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{X = GAMMA(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab Shall be of type @code{REAL} and neither zero
+nor a negative integer.
+@end multitable
+
+@item @emph{Return value}:
+The return value is of type @code{REAL} of the same kind as @var{X}.
+
+@item @emph{Example}:
+@smallexample
+program test_gamma
+ real :: x = 1.0
+ x = gamma(x) ! returns 1.0
+end program test_gamma
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name @tab Argument @tab Return type @tab Standard
+@item @code{GAMMA(X)} @tab @code{REAL(4) X} @tab @code{REAL(4)} @tab GNU Extension
+@item @code{DGAMMA(X)} @tab @code{REAL(8) X} @tab @code{REAL(8)} @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Logarithm of the Gamma function: @ref{LGAMMA}
+
+@end table
+
+
+
@node GERROR
@section @code{GERROR} --- Get last system error message
@fnindex GERROR
@@ -6230,6 +6381,60 @@ The return value is of type @code{INTEGE
+@node LGAMMA
+@section @code{LGAMMA} --- Logarithm of the Gamma function
+@fnindex GAMMA
+@fnindex ALGAMA
+@fnindex DLGAMA
+@cindex Gamma function, logarithm of
+@cindex
+
+@table @asis
+@item @emph{Description}:
+@code{GAMMA(X)} computes the natural logrithm of the absolute value of the
+Gamma (@math{\Gamma}) function.
+
+@item @emph{Standard}:
+GNU Extension
+
+@item @emph{Class}:
+Elemental function
+
+@item @emph{Syntax}:
+@code{X = LGAMMA(X)}
+
+@item @emph{Arguments}:
+@multitable @columnfractions .15 .70
+@item @var{X} @tab Shall be of type @code{REAL} and neither zero
+nor a negative integer.
+@end multitable
+
+@item @emph{Return value}:
+The return value is of type @code{REAL} of the same kind as @var{X}.
+
+@item @emph{Example}:
+@smallexample
+program test_log_gamma
+ real :: x = 1.0
+ x = lgamma(x) ! returns 0.0
+end program test_log_gamma
+@end smallexample
+
+@item @emph{Specific names}:
+@multitable @columnfractions .20 .20 .20 .25
+@item Name @tab Argument @tab Return type @tab Standard
+@item @code{LGAMMA(X)} @tab @code{REAL(4) X} @tab @code{REAL(4)} @tab GNU Extension
+@item @code{ALGAMA(X)} @tab @code{REAL(4) X} @tab @code{REAL(4)} @tab GNU Extension
+@item @code{DLGAMA(X)} @tab @code{REAL(8) X} @tab @code{REAL(8)} @tab GNU Extension
+@end multitable
+
+@item @emph{See also}:
+Gamma function: @ref{GAMMA}
+
+@end table
+
+
+
@node LGE
@section @code{LGE} --- Lexical greater than or equal
@fnindex LGE
Index: gcc/testsuite/gfortran.dg/gamma_1.f90
===================================================================
--- gcc/testsuite/gfortran.dg/gamma_1.f90 (revision 0)
+++ gcc/testsuite/gfortran.dg/gamma_1.f90 (revision 0)
@@ -0,0 +1,31 @@
+! { dg-do run }
+!
+! Test the vendor intrinsic (d)gamma, lgamma and algama/dlgama
+! gamma is also part of the Fortran 2008 draft; lgamma is called
+! log_gamma in the Fortran 2008 draft.
+!
+! PR fortran/32980
+!
+program gamma_test
+implicit none
+intrinsic :: gamma, lgamma
+integer, parameter :: sp = kind(1.0)
+integer, parameter :: dp = kind(1.0d0)
+integer, parameter :: qp = selected_real_kind(p=15,r=900)
+
+real(sp) :: rsp
+real(dp) :: rdp
+real(dp) :: rqp
+
+if (abs(gamma(1.0_sp) - 1.0_sp) > tiny(1.0_sp)) call abort()
+if (abs(gamma(1.0_dp) - 1.0_dp) > tiny(1.0_dp)) call abort()
+if (abs(gamma(1.0_qp) - 1.0_qp) > tiny(1.0_qp)) call abort()
+if (abs(dgamma(1.0_dp) - 1.0_dp) > tiny(1.0_dp)) call abort()
+
+if (abs(lgamma(1.0_sp)) > tiny(1.0_sp)) call abort()
+if (abs(lgamma(1.0_dp)) > tiny(1.0_dp)) call abort()
+if (abs(lgamma(1.0_qp)) > tiny(1.0_qp)) call abort()
+if (abs(algama(1.0_sp)) > tiny(1.0_sp)) call abort()
+if (abs(dlgama(1.0_dp)) > tiny(1.0_dp)) call abort()
+end program gamma_test
+
Index: gcc/testsuite/gfortran.dg/gamma_2.f90
===================================================================
--- gcc/testsuite/gfortran.dg/gamma_2.f90 (revision 0)
+++ gcc/testsuite/gfortran.dg/gamma_2.f90 (revision 0)
@@ -0,0 +1,36 @@
+! { dg-do compile }
+! { dg-options "-std=f2003 -Wall" }
+!
+! Test the vendor intrinsic (d)gamma, lgamma and algama/dlgama
+! gamma is also part of the Fortran 2008 draft; lgamma is called
+! log_gamma in the Fortran 2008 draft.
+!
+! PR fortran/32980
+!
+subroutine foo()
+intrinsic :: gamma
+intrinsic :: dgamma
+intrinsic :: lgamma
+intrinsic :: algama
+intrinsic :: dlgama
+
+integer, parameter :: sp = kind(1.0)
+integer, parameter :: dp = kind(1.0d0)
+integer, parameter :: qp = selected_real_kind(p=15,r=900)
+
+real(sp) :: rsp = 1.0_sp
+real(dp) :: rdp = 1.0_dp
+real(dp) :: rqp = 1.0_qp
+
+rsp = gamma(rsp) ! FIXME: "is not included in the selected standard"
+rdp = gamma(rdp) ! FIXME: "is not included in the selected standard"
+rqp = gamma(rqp) ! FIXME: "is not included in the selected standard"
+rdp = dgamma(rdp) ! { dg-error "is not included in the selected standard" }
+
+rsp = lgamma(rsp) ! FIXME: "is not included in the selected standard"
+rdp = lgamma(rdp) ! FIXME: "is not included in the selected standard"
+rqp = lgamma(rqp) ! FIXME: "is not included in the selected standard"
+rsp = algama(rsp) ! { dg-error "is not included in the selected standard" }
+rdp = dlgama(rdp) ! { dg-error "is not included in the selected standard" }
+end subroutine foo
+end
Index: gcc/testsuite/gfortran.dg/gamma_3.f90
===================================================================
--- gcc/testsuite/gfortran.dg/gamma_3.f90 (revision 0)
+++ gcc/testsuite/gfortran.dg/gamma_3.f90 (revision 0)
@@ -0,0 +1,27 @@
+! { dg-do compile }
+!
+! Test the vendor intrinsic (d)gamma, lgamma and algama/dlgama
+! gamma is also part of the Fortran 2008 draft; lgamma is called
+! log_gamma in the Fortran 2008 draft.
+!
+! PR fortran/32980
+!
+program gamma_test
+implicit none
+intrinsic :: gamma, lgamma
+real :: x
+
+x = gamma(cmplx(1.0,0.0)) ! { dg-error "is not consistent with a specific intrinsic interface" }
+x = dgamma(cmplx(1.0,0.0,kind(0d0))) ! { dg-error "must be REAL" }
+x = gamma(int(1)) ! { dg-error "is not consistent with a specific intrinsic interface" }
+x = dgamma(int(1)) ! { dg-error "must be REAL" }
+
+x = lgamma(cmplx(1.0,0.0)) ! { dg-error "is not consistent with a specific intrinsic interface" }
+x = algama(cmplx(1.0,0.0)) ! { dg-error "must be REAL" }
+x = dlgama(cmplx(1.0,0.0,kind(0d0))) ! { dg-error "must be REAL" }
+
+x = lgamma(int(1)) ! { dg-error "is not consistent with a specific intrinsic interface" }
+x = algama(int(1)) ! { dg-error "must be REAL" }
+x = dlgama(int(1)) ! { dg-error "must be REAL" }
+end program gamma_test
+