add gamma function an other version

git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1485 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 229 insertions, 0 deletions

diff --git a/gammaPiAGMformula.c b/gammaPiAGMformula.c
new file mode 100644
index 000000000..5b580f023
--- /dev/null
+++ b/gammaPiAGMformula.c
@@ -0,0 +1,229 @@
+/* mpfr_gamma -- gamma function
+
+Copyright (C) 2001 Free Software Foundation.
+
+This file is part of the MPFR Library, and was contributed by Mathieu Dutour.
+
+The MPFR Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Library General Public License as published by
+the Free Software Foundation; either version 2 of the License, or (at your
+option) any later version.
+
+The MPFR Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
+License for more details.
+
+You should have received a copy of the GNU Library General Public License
+along with the MPFR Library; see the file COPYING.LIB.  If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <math.h>
+#include "gmp.h"
+#include "gmp-impl.h"
+#include "mpfr.h"
+#include "mpfr-impl.h"
+
+int mpfr_gamma _PROTO ((mpfr_ptr, mpfr_srcptr, mp_rnd_t));
+
+/* We use the reflection formula 
+  Gamma(1+x)Gamma(1-x)=\pi x/(sin(\pi x))
+  in order to treat the case x<=1  */
+
+#define CST   0.38  /* CST=ln(2)/(ln(2*pi)) */
+#define zCST  0.26  /* zCST=1/(2*ln(2*pi)) */
+
+
+int
+#if __STDC__
+mpfr_gamma (mpfr_ptr gamma, mpfr_srcptr x, mp_rnd_t rnd_mode)
+#else
+mpfr_gamma (gamma, x, rnd_mode)
+     mpfr_ptr gamma;
+     mpfr_srcptr x; 
+     mp_rnd_t rnd_mode;
+#endif
+{
+  mpfr_t xp;
+  mpfr_t product;
+  mpfr_t driv;
+  mpfr_t xminusone;
+  mpfr_t the_pi_constant;
+  mpfr_t Csin, Ccos;
+  mpfr_t GammaTrial;
+  int reflex;
+
+  mpfr_t tmp, tmp2;
+  int Prec;
+  int Prec_nec;
+  int prec_gamma;
+  int prec_nec;
+  int good = 0;
+  double C;
+  long A, N;
+  int realprec;
+  int estimated_delta;
+  int compared; 
+  int compar; 
+  int factorial_case;
+  int k;
+  
+  /* Trivial cases */
+  if (MPFR_IS_NAN(x))
+    {
+      MPFR_SET_NAN(gamma);
+      return 1;
+    }
+  if (!MPFR_NOTZERO(x))
+    {
+      MPFR_SET_INF(gamma);
+      return 1;
+    }
+  if (MPFR_IS_INF(x))
+    {
+      MPFR_SET_INF(gamma);
+      return 1;
+    }
+  /* Set x_p=x if x> 1 else set x_p=2-x */
+  prec_gamma = MPFR_PREC(gamma);
+  compared = mpfr_cmp_ui(x, 1);
+  if (compared == 0)
+    {
+      mpfr_set_ui(gamma, 1, rnd_mode);
+      return 1;
+    }
+  if (compared == -1)
+    {
+      prec_nec = 2+prec_gamma;
+      reflex = 1;
+      return 1;
+    }
+  else
+    {
+      prec_nec = prec_gamma;
+      reflex = 0;
+      return 1;
+    }
+
+  realprec = prec_nec+10;
+
+  while (!good){
+    C = ((double) realprec)*CST-0.5;
+    A = (long) ceil(C-zCST*log(C));
+    N = A-1;
+    /* Compute the correct estimated_delta as a function of realprec */
+    Prec = realprec+estimated_delta;
+    mpfr_init2(xp, Prec);
+    if (compared == -1)
+      {
+	mpfr_ui_sub(xp, 2, x, GMP_RNDN);
+      }
+    else
+      {
+	mpfr_set(xp, x, GMP_RNDN);
+      }
+    /* Initialisation    */
+    mpfr_init2(product, Prec);
+    mpfr_init2(driv, Prec);
+    mpfr_init2(tmp, Prec);
+    mpfr_set(driv, xp, GMP_RNDN);
+    mpfr_set_ui(product, 1, GMP_RNDN);
+
+    /* We use a naugthy algorithm to reduce to the domain 1<=x <2*/
+    while (1)
+      {
+	compar = mpfr_cmp_ui(driv, 2);
+	if (compar == 0) /* the factorial in fact */
+	  {
+	    mpfr_mul_ui(GammaTrial, product, 2, GMP_RNDN);
+	    factorial_case = 1;
+	    break;
+	  }
+	if (compar == -1)
+	  {
+	    factorial_case = 0;
+	    break;
+	  }
+	mpfr_sub_ui(driv, driv, 1, GMP_RNDN);
+	mpfr_mul(product, product, driv, GMP_RNDN);
+      }
+    /* now we run into the trivial case */
+    if (factorial_case == 0 && reflex == 1)
+      {
+	mpfr_clear(product);
+	mpfr_clear(driv);
+	mpfr_clear(tmp);
+	MPFR_SET_INF(gamma);
+	return 1;
+      }
+    mpfr_init2(GammaTrial, Prec);
+    if (factorial_case == 0)
+      {
+	/* compute the Gamma for 1< driv < 2 */
+	mpfr_init2(tmp2, Prec);
+
+	mpfr_add_ui(tmp2, driv,6*N, GMP_RNDN);
+	mpfr_ui_div(tmp, 1, tmp2, GMP_RNDN);
+	for(k=6*N; k>=1; k--)
+	  {
+	    mpfr_mul_ui(tmp, tmp, N, GMP_RNDN);
+	    mpfr_neg(tmp, tmp, GMP_RNDN);
+	    mpfr_div_ui(tmp, tmp, k, GMP_RNDN);
+	    
+	    mpfr_add_ui(tmp2, driv, k-1, GMP_RNDN);
+	    mpfr_ui_div(tmp2, 1, tmp2, GMP_RNDN);
+	    mpfr_add(tmp, tmp, tmp2, GMP_RNDN);
+	  }
+	mpfr_set_ui(tmp2, N, GMP_RNDN);
+	mpfr_log(tmp2, tmp2, GMP_RNDN);
+	mpfr_mul(tmp2, tmp2, driv, GMP_RNDN);
+	mpfr_exp(tmp2, tmp2, GMP_RNDN);
+	mpfr_mul(tmp, tmp2, tmp, GMP_RNDN);
+	mpfr_mul(GammaTrial, tmp, product, GMP_RNDN);
+      }
+    if (reflex == 1)
+      {
+	mpfr_init2(xminusone, Prec);
+	mpfr_init2(Csin, Prec);
+	mpfr_init2(Ccos, Prec);
+	mpfr_sub_ui(xminusone, x, 1, GMP_RNDN);
+	mpfr_const_pi(the_pi_constant, Prec_nec);
+	mpfr_mul(tmp, the_pi_constant, xminusone, GMP_RNDN);
+	mpfr_sin(Csin, tmp, GMP_RNDN);
+	mpfr_cos(Ccos, tmp, GMP_RNDN);
+	mpfr_div(tmp, tmp, Csin, GMP_RNDN);
+	mpfr_neg(tmp, tmp, GMP_RNDN);
+	mpfr_div(GammaTrial, tmp, GammaTrial, GMP_RNDN);
+      }
+    if (mpfr_can_round(GammaTrial, realprec, GMP_RNDD, rnd_mode, MPFR_PREC(gamma)))
+      {
+	mpfr_set(gamma, GammaTrial, rnd_mode);
+	good = 1;
+      }
+    else
+      {
+	realprec += _mpfr_ceil_log2 ((double) realprec);
+#ifdef DEBUG
+	printf("RETRY\n");
+#endif
+      }
+    mpfr_clear(tmp);
+    mpfr_clear(driv);
+    mpfr_clear(product);
+    mpfr_clear(GammaTrial);
+    if (reflex == 1)
+      {
+	mpfr_clear(xminusone);
+	mpfr_clear(Csin);
+	mpfr_clear(Ccos);
+      }
+  }
+
+  mpfr_clear(xp);
+
+  return 1; /* inexact result */
+}
+