2007-08-26 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-26  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.



git-svn-id: svn+ssh://gcc.gnu.org/svn/gcc/trunk@127809 138bc75d-0d04-0410-961f-82ee72b054a4

11 files changed, 316 insertions, 0 deletions

diff --git a/gcc/fortran/ChangeLog b/gcc/fortran/ChangeLog
index 9adaf0791c7..259716415b2 100644
--- a/gcc/fortran/ChangeLog
+++ b/gcc/fortran/ChangeLog
@@ -1,5 +1,17 @@
 2007-08-26  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-26  Tobias Burnus  <burnus@net-b.de>
+
 	PR fortran/33188
 	* parse.c (parse_derived): Support empty derived type
 	definitions for Fortran 2003.
diff --git a/gcc/fortran/intrinsic.c b/gcc/fortran/intrinsic.c
index 4128e441d63..d273f806111 100644
--- a/gcc/fortran/intrinsic.c
+++ b/gcc/fortran/intrinsic.c
@@ -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);
diff --git a/gcc/fortran/intrinsic.h b/gcc/fortran/intrinsic.h
index cf242b8995d..e284a6c7177 100644
--- a/gcc/fortran/intrinsic.h
+++ b/gcc/fortran/intrinsic.h
@@ -221,6 +221,7 @@ gfc_expr *gfc_simplify_exponent (gfc_expr *);
 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_expr *);
 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_expr *, gfc_expr *);
 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 *);
diff --git a/gcc/fortran/intrinsic.texi b/gcc/fortran/intrinsic.texi
index f09246cc95e..e94a7e3023b 100644
--- a/gcc/fortran/intrinsic.texi
+++ b/gcc/fortran/intrinsic.texi
@@ -121,6 +121,7 @@ Some basic guidelines for editing this document:
 * @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 +162,7 @@ Some basic guidelines for editing this document:
 * @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 +4486,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 +6291,60 @@ The return value is of type @code{INTEGER} and of kind @var{KIND}. If
 
 
 
+@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
diff --git a/gcc/fortran/iresolve.c b/gcc/fortran/iresolve.c
index fc837e1fcd0..7948b14eeb9 100644
--- a/gcc/fortran/iresolve.c
+++ b/gcc/fortran/iresolve.c
@@ -757,6 +757,15 @@ gfc_resolve_g77_math1 (gfc_expr *f, gfc_expr *x)
 
 
 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;
@@ -1114,6 +1123,15 @@ gfc_resolve_len_trim (gfc_expr *f, gfc_expr *string, gfc_expr *kind)
 
 
 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)
 {
diff --git a/gcc/fortran/mathbuiltins.def b/gcc/fortran/mathbuiltins.def
index 0fc73688279..33e87d14290 100644
--- a/gcc/fortran/mathbuiltins.def
+++ b/gcc/fortran/mathbuiltins.def
@@ -30,3 +30,5 @@ DEFINE_MATH_BUILTIN   (Y1,    "y1",     0)
 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)
diff --git a/gcc/fortran/simplify.c b/gcc/fortran/simplify.c
index a395b04a599..5e129ae18de 100644
--- a/gcc/fortran/simplify.c
+++ b/gcc/fortran/simplify.c
@@ -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_expr *kind)
   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)
diff --git a/gcc/testsuite/ChangeLog b/gcc/testsuite/ChangeLog
index 37d280eee28..8dbd4c2c9e7 100644
--- a/gcc/testsuite/ChangeLog
+++ b/gcc/testsuite/ChangeLog
@@ -1,5 +1,12 @@
 2007-08-26  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.
+
+2007-08-26  Tobias Burnus  <burnus@net-b.de>
+
 	PR fortran/33188
 	* gfortran.dg/used_types_18.f90: New.
 
diff --git a/gcc/testsuite/gfortran.dg/gamma_1.f90 b/gcc/testsuite/gfortran.dg/gamma_1.f90
new file mode 100644
index 00000000000..819f4672176
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/gamma_1.f90
@@ -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
+
diff --git a/gcc/testsuite/gfortran.dg/gamma_2.f90 b/gcc/testsuite/gfortran.dg/gamma_2.f90
new file mode 100644
index 00000000000..b15696fe350
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/gamma_2.f90
@@ -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
diff --git a/gcc/testsuite/gfortran.dg/gamma_3.f90 b/gcc/testsuite/gfortran.dg/gamma_3.f90
new file mode 100644
index 00000000000..b35596fa5a9
--- /dev/null
+++ b/gcc/testsuite/gfortran.dg/gamma_3.f90
@@ -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
+