自学教程：C++ tgamma函数代码示例

static double		student_density(const copula		*cop,					const double		x){    double t;    t=t1_student_density;      return ( tgamma((t+1)/2)/(tgamma(t/2)*sqrt(M_PI * t)*pow(1+x*x/t,(t+1)/2)) ); //fdr de la loi student a t deg de lib}

void test_tgamma(){    static_assert((std::is_same<decltype(tgamma((double)0)), double>::value), "");    static_assert((std::is_same<decltype(tgammaf(0)), float>::value), "");    static_assert((std::is_same<decltype(tgammal(0)), long double>::value), "");    assert(tgamma(1) == 1);}

/** * @function beta */double beta( double z,double w) {  double gz, gw, gzw;    gz = tgamma(z);  gw = tgamma(w);  gzw = tgamma(z+w);  return  gz*gw/gzw;}

 Beta(Teuchos::ParameterList &parlist) {   shape1_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 1",2.);   shape2_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 2",2.);   shape1_ = (shape1_ > 0.) ? shape1_ : 2.;   shape2_ = (shape2_ > 0.) ? shape2_ : 2.;   coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_));   initializeQuadrature(); }

 Real moment(const size_t m) const {   Real val = 0., binom = 0., moment = 0.;   for (size_t i = 0; i < m; i++) {     moment = exp2_*tgamma(1.+(Real)i/exp1_)*tgamma(exp2_)/tgamma(1.+exp2_+(Real)i/exp1_);     binom  = (Real)nchoosek(m,i);     val   += binom*std::pow(a_,m-i)*std::pow(b_-a_,i+1)*moment;   }       return val;  }

double pnorm_taylor_expn(const double x,			 const int N) {    double pnorm_v = 0.0;    #pragma omp parallel for reduction(+:pnorm_v)    for (int n = 0; n < N; n++) {	pnorm_v += pow(-0.5, n) * pow(x, 2*n+1) / (tgamma(n+1) * (2*n+1));    }    pnorm_v = pnorm_v / (sqrt(2.0) * tgamma(0.5)) + 0.5;    return pnorm_v;}

void			init_cdf_doublet(const double rho,					   const double t1,					   const double t2,					   int n,					   int m){//  double 		**cdf_doublet;    double		coef1;			//coef pour la cdf    double		coef2;			//coef pour la deuxieme integrale    int 		i=0;    int			j=0;    double		a1;			//borne inf de la cdf    double		b1;			//borne sup de la cdf    double		a2;			//borne inf de la deuxieme integrale    double		b2;			//borne sup de la deuxieme integrale    double		s1;			//variable de stockage    double		s2;			//autre variable de stockage            cdf=malloc(2*sizeof(double*));		//on declare le tableau pour la cdf    cdf[0]=malloc((n+1)*sizeof(double));	//on declare la colonne des points    cdf[1]=malloc((n+1)*sizeof(double));	//on declare la colonne des valeurs aux points    coef1 = (tgamma((t1+1)/2) * tgamma((t2+1)/2));    coef1 = coef1 / (tgamma(t1/2) * tgamma(t2/2) * M_PI * rho * sqrt((1 - rho*rho) * (t1 - 2) * (t2 - 2)));     a1=dp[0];    b1=dp[1];    a2=dp[2];    b2=dp[3];    coef1 = coef1*(b1 - a1)/n;			//methode des trapezes    coef2 = (b2 - a2)/m;				//methode des trapezes    printf("test : f_cdf = %f /n",pow(a1,2));    s2=0;        printf("avant boucle/n");    for(i=0;i<n+1;i++){    cdf[0][i]= a1 + (b1 - a1)*i/( (double )n); 		//point considere = x1;    s1=0;    s1+=f_cdf(rho,t1,t2,cdf[0][i],a2)/2;    s1+=f_cdf(rho,t1,t2,cdf[0][i],b2)/2;    	for(j=1;j<m;j++)	s1+=f_cdf(rho,t1,t2,cdf[0][i],a2 + (b2 - a2)*j/((double) m));	//methode des trapezes    if (i>0) {  //  printf("s1 = %f, s2= %f ,/n",s1,s2);    cdf[1][i]=(s1+s2)/2*coef1*coef2+cdf[1][i-1];		//methode des trapezes  //  printf("X=%g, Y=%g, i=%f/n", cdf[1][i], cdf[1][i-1], cdf[0][i]);    }    else    cdf[1][i]=0;						//premiere boucle, initialisation        s2=s1;    }	  //  printf("ok %f, %f, %f/n",cdf[0][n],cdf[1][n],cdf[1][50]);    return;}

//function definitionsvoid AGGDfit(IplImage* structdis, double& lsigma_best, double& rsigma_best, double& gamma_best){	BwImage ImArr(structdis);		//int width = structdis->width;	//int height = structdis->height;	long int poscount=0, negcount=0;	double possqsum=0, negsqsum=0, abssum=0;	for(int i=0;i<structdis->height;i++)	{		for (int j =0; j<structdis->width; j++)		{			double pt = ImArr[i][j];			if(pt>0)			{				poscount++;				possqsum += pt*pt; 				abssum += pt;			}			else if(pt<0)			{				negcount++;				negsqsum += pt*pt;				abssum -= pt;			}		}	}	lsigma_best = pow(negsqsum/negcount, 0.5); //1st output parameter set	rsigma_best = pow(possqsum/poscount, 0.5); //2nd output parameter set	 	double gammahat = lsigma_best/rsigma_best;	long int totalcount = structdis->width*structdis->height;	double rhat = pow(abssum/totalcount, static_cast<double>(2))/((negsqsum + possqsum)/totalcount);	double rhatnorm = rhat*(pow(gammahat,3) +1)*(gammahat+1)/pow(pow(gammahat,2)+1,2);		double prevgamma = 0;	double prevdiff = 1e10;	        float sampling = 0.001;	for (float gam=0.2; gam<10; gam+=sampling) //possible to coarsen sampling to quicken the code, with some loss of accuracy	{		double r_gam = tgamma(2/gam)*tgamma(2/gam)/(tgamma(1/gam)*tgamma(3/gam));		double diff = abs(r_gam-rhatnorm);		if(diff> prevdiff) break;		prevdiff = diff;		prevgamma = gam;	}	gamma_best = prevgamma;}

int main(){    std::cout << "Result of erf_inv(-10) is: "              << erf_inv(-10) << std::endl;    std::cout << "Result of tgamma(-10) is: "              << tgamma(-10) << std::endl;}

static doubleneg_gam(double x){	int sgn = 1;	struct Double lg, lsine;	double y, z;	y = ceil(x);	if (y == x)		/* Negative integer. */		if (_IEEE)			return ((x - x) / zero);		else			return (infnan(ERANGE));	z = y - x;	if (z > 0.5)		z = one - z;	y = 0.5 * y;	if (y == ceil(y))		sgn = -1;	if (z < .25)		z = sin(M_PI*z);	else		z = cos(M_PI*(0.5-z));	/* Special case: G(1-x) = Inf; G(x) may be nonzero. */	if (x < -170) {		if (x < -190)			return ((double)sgn*tiny*tiny);		y = one - x;		/* exact: 128 < |x| < 255 */		lg = large_gam(y);		lsine = __log__D(M_PI/z);	/* = TRUNC(log(u)) + small */		lg.a -= lsine.a;		/* exact (opposite signs) */		lg.b -= lsine.b;		y = -(lg.a + lg.b);		z = (y + lg.a) + lg.b;		y = __exp__D(y, z);		if (sgn < 0) y = -y;		return (y);	}	y = one-x;	if (one-y == x)		y = tgamma(y);	else		/* 1-x is inexact */		y = -x*tgamma(-x);	if (sgn < 0) y = -y;	return (M_PI / (y*z));}

 Gamma(Teuchos::ParameterList &parlist) {   shape_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gamma").get("Shape",1.);   scale_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gamma").get("Scale",1.);   shape_ = (shape_ > 0.) ? shape_ : 1.;   scale_ = (scale_ > 0.) ? scale_ : 1.;   gamma_shape_ = tgamma(shape_);   coeff_       = 1./(gamma_shape_*std::pow(scale_,shape_)); }

static PetscErrorCode PetscDTGaussJacobiQuadrature1D_Internal(PetscInt npoints, PetscReal a, PetscReal b, PetscReal *x, PetscReal *w){  PetscInt       maxIter = 100;  PetscReal      eps     = 1.0e-8;  PetscReal      a1, a2, a3, a4, a5, a6;  PetscInt       k;  PetscErrorCode ierr;  PetscFunctionBegin;  a1      = pow(2, a+b+1);#if defined(PETSC_HAVE_TGAMMA)  a2      = tgamma(a + npoints + 1);  a3      = tgamma(b + npoints + 1);  a4      = tgamma(a + b + npoints + 1);#else  SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"tgamma() - math routine is unavailable.");#endif  ierr = PetscDTFactorial_Internal(npoints, &a5);CHKERRQ(ierr);  a6   = a1 * a2 * a3 / a4 / a5;  /* Computes the m roots of P_{m}^{a,b} on [-1,1] by Newton's method with Chebyshev points as initial guesses.   Algorithm implemented from the pseudocode given by Karniadakis and Sherwin and Python in FIAT */  for (k = 0; k < npoints; ++k) {    PetscReal r = -cos((2.0*k + 1.0) * PETSC_PI / (2.0 * npoints)), dP;    PetscInt  j;    if (k > 0) r = 0.5 * (r + x[k-1]);    for (j = 0; j < maxIter; ++j) {      PetscReal s = 0.0, delta, f, fp;      PetscInt  i;      for (i = 0; i < k; ++i) s = s + 1.0 / (r - x[i]);      ierr = PetscDTComputeJacobi(a, b, npoints, r, &f);CHKERRQ(ierr);      ierr = PetscDTComputeJacobiDerivative(a, b, npoints, r, &fp);CHKERRQ(ierr);      delta = f / (fp - f * s);      r     = r - delta;      if (fabs(delta) < eps) break;    }    x[k] = r;    ierr = PetscDTComputeJacobiDerivative(a, b, npoints, x[k], &dP);CHKERRQ(ierr);    w[k] = a6 / (1.0 - PetscSqr(x[k])) / PetscSqr(dP);  }  PetscFunctionReturn(0);}

 Real moment(const size_t m) const {   if ( m == 1 ) {     return shape_*scale_;   }   if ( m == 2 ) {     return shape_*scale_*scale_*(1. + shape_);   }   return std::pow(scale_,m)*tgamma(shape_+(Real)m)/gamma_shape_;  }

int main(int argc, char *argv[]){	double x = 5.5, y = 3;	double a, b;	printf("erfc(%lg): %lg, %lg/n", x, erfc(x), kf_erfc(x));	printf("upper-gamma(%lg,%lg): %lg/n", x, y, kf_gammaq(y, x)*tgamma(y));	a = 2; b = 2; x = 0.5;	printf("incomplete-beta(%lg,%lg,%lg): %lg/n", a, b, x, kf_betai(a, b, x) / exp(kf_lgamma(a+b) - kf_lgamma(a) - kf_lgamma(b)));	return 0;}

double normalize(double exponent, int ilx, int ily, int ilz) {    double coeff, ltot, power;  double glx, gly, glz;  double lx = (double) ilx;   double ly = (double) ily;   double lz = (double) ilz;   ltot = lx+ly+lz;  power = (ltot+1.5)/2.0;  glx = tgamma(lx+0.5);  gly = tgamma(ly+0.5);  glz = tgamma(lz+0.5);   coeff = pow(2.0*exponent,power)/sqrt(glx*gly*glz);  return coeff;} 

void get_cuckoos(int nnest, int lb[], int ub[], int dim, Nest *best, Nest *new_Nest, Nest *Nests, RngStream g1){	double beta = 3/2;	double sigma = pow((tgamma(1.0+beta) * sin(M_PI*beta/2.0) / ((tgamma(1.0+beta)/2.0) * beta * (pow(2.0,((beta-1.0)/2.0))))),(1.0/beta));	int i;	for (i = 0; i < nnest; i++){		int j;		for (j = 0; j < dim; j++){			double step = (RngStream_RandU01(g1)*sigma) / pow(fabs(RngStream_RandU01(g1)),(1/beta));			double stepsize = 0.01 * step * (best->pos[j]);			new_Nest[i].pos[j] = Nests[i].pos[j] + stepsize * RngStream_RandU01(g1);		}			}	simplebounds(nnest, dim, new_Nest, lb, ub);}

double pnorm_riemann_intl(const double x, const double dx) {    double pnorm_v = 0.0;    double dn = 0.0;    #pragma omp parallel for reduction(+:pnorm_v)    for (int i = 0; i < (int) (x/dx); i++) {	dn = i * dx;	pnorm_v += exp(-0.5 * pow(dn, 2.0));    }    pnorm_v = dx * pnorm_v / (sqrt(2.0) * tgamma(0.5)) + 0.5;     return pnorm_v; }

static inline doubleruby_tgamma(const double d){    const double g = tgamma(d);    if (isinf(g)) {	if (d == 0.0 && signbit(d)) return -INFINITY;    }    if (isnan(g)) {	if (!signbit(d)) return INFINITY;    }    return g;}

sl_def(t_main, void){  double x = slr_get(x)[0];  double t1 = tgamma(x);  float t2 = tgammaf(x);  output_float(x, 1, 4);  output_char(' ', 1);  output_float(t1, 1, 4);  output_char(' ', 1);  output_float(t2, 1, 4);  output_char('/n', 1);}

static VALUEmath_gamma(VALUE obj, SEL sel, VALUE x){    static const double fact_table[] = {        /* fact(0) */ 1.0,        /* fact(1) */ 1.0,        /* fact(2) */ 2.0,        /* fact(3) */ 6.0,        /* fact(4) */ 24.0,        /* fact(5) */ 120.0,        /* fact(6) */ 720.0,        /* fact(7) */ 5040.0,        /* fact(8) */ 40320.0,        /* fact(9) */ 362880.0,        /* fact(10) */ 3628800.0,        /* fact(11) */ 39916800.0,        /* fact(12) */ 479001600.0,        /* fact(13) */ 6227020800.0,        /* fact(14) */ 87178291200.0,        /* fact(15) */ 1307674368000.0,        /* fact(16) */ 20922789888000.0,        /* fact(17) */ 355687428096000.0,        /* fact(18) */ 6402373705728000.0,        /* fact(19) */ 121645100408832000.0,        /* fact(20) */ 2432902008176640000.0,        /* fact(21) */ 51090942171709440000.0,        /* fact(22) */ 1124000727777607680000.0,        /* fact(23)=25852016738884976640000 needs 56bit mantissa which is         * impossible to represent exactly in IEEE 754 double which have         * 53bit mantissa. */    };    double d0, d;    double intpart, fracpart;    Need_Float(x);    d0 = RFLOAT_VALUE(x);    /* check for domain error */    if (isinf(d0) && signbit(d0)) {	domain_error("gamma");    }    fracpart = modf(d0, &intpart);    if (fracpart == 0.0) {	if (intpart < 0) {	    domain_error("gamma");	}	if (0 < intpart &&	    intpart - 1 < (double)numberof(fact_table)) {	    return DBL2NUM(fact_table[(int)intpart - 1]);	}    }    d = tgamma(d0);    return DBL2NUM(d);}

static VALUEmath_gamma(VALUE unused_obj, VALUE x){    static const double fact_table[] = {        /* fact(0) */ 1.0,        /* fact(1) */ 1.0,        /* fact(2) */ 2.0,        /* fact(3) */ 6.0,        /* fact(4) */ 24.0,        /* fact(5) */ 120.0,        /* fact(6) */ 720.0,        /* fact(7) */ 5040.0,        /* fact(8) */ 40320.0,        /* fact(9) */ 362880.0,        /* fact(10) */ 3628800.0,        /* fact(11) */ 39916800.0,        /* fact(12) */ 479001600.0,        /* fact(13) */ 6227020800.0,        /* fact(14) */ 87178291200.0,        /* fact(15) */ 1307674368000.0,        /* fact(16) */ 20922789888000.0,        /* fact(17) */ 355687428096000.0,        /* fact(18) */ 6402373705728000.0,        /* fact(19) */ 121645100408832000.0,        /* fact(20) */ 2432902008176640000.0,        /* fact(21) */ 51090942171709440000.0,        /* fact(22) */ 1124000727777607680000.0,        /* fact(23)=25852016738884976640000 needs 56bit mantissa which is         * impossible to represent exactly in IEEE 754 double which have         * 53bit mantissa. */    };    enum {NFACT_TABLE = numberof(fact_table)};    double d;    d = Get_Double(x);    /* check for domain error */    if (isinf(d)) {	if (signbit(d)) domain_error("gamma");	return DBL2NUM(HUGE_VAL);    }    if (d == 0.0) {	return signbit(d) ? DBL2NUM(-HUGE_VAL) : DBL2NUM(HUGE_VAL);    }    if (d == floor(d)) {	if (d < 0.0) domain_error("gamma");	if (1.0 <= d && d <= (double)NFACT_TABLE) {	    return DBL2NUM(fact_table[(int)d - 1]);	}    }    return DBL2NUM(tgamma(d));}

void rec_jacobi( const double alpha,                  const double beta,                 ROL::Vector<Real> &a,                 ROL::Vector<Real> &b ) {    Teuchos::RCP<std::vector<Real> > ap =         Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(a)).getVector());     Teuchos::RCP<std::vector<Real> > bp =         Teuchos::rcp_const_cast<std::vector<Real> >((Teuchos::dyn_cast<ROL::StdVector<Real> >(b)).getVector());     int N = ap->size();    Real nu = (beta-alpha)/double(alpha+beta+2.0);    Real mu  = pow(2.0,alpha+beta+1.0)*tgamma(alpha+1.0)*tgamma(beta+1.0)/tgamma(alpha+beta+2.0);    Real nab;    Real sqdif = pow(beta,2)-pow(alpha,2);                                  (*ap)[0] = nu;    (*bp)[0] = mu;        if(N>1){             for(int n=1;n<N;++n) {            nab = 2*n+alpha+beta;            (*ap)[n] = sqdif/(nab*(nab+2));                      }        (*bp)[1] = 4.0*(alpha+1.0)*(beta+1.0)/(pow(alpha+beta+2.0,2)*(alpha+beta+3.0));        if(N>2) {            for(int n=2;n<N;++n) {                nab = 2*n+alpha+beta;                (*bp)[n] = 4.0*(n+alpha)*(n+beta)*n*(n+alpha+beta)/(nab*nab*(nab+1.0)*(nab-1.0));            }        }    }    }

static TACommandVerdict tgamma_cmd(TAThread thread,TAInputStream stream){    double x, res;    x = readDouble(&stream);    START_TARGET_OPERATION(thread);    errno = 0;    res = tgamma(x);    END_TARGET_OPERATION(thread);    writeInt(thread, errno);    writeDouble(thread, res);    sendResponse(thread);    return taDefaultVerdict;}

// get value for the leading order exponential SF double EvtPFermi::getSFBLNP(const double &what){  double SF;  double massB = 5.2792;     if ( what      > massB )   return 0;  if ( what < 0 ) return 0; #if defined(__SUNPRO_CC)  report(Severity::Error,"EvtGen") << "The tgamma function is not available on this platform/n";  report(Severity::Error,"EvtGen") <<"Presumably, you are getting the wrong answer, so I abort..";  ::abort();#else  SF = pow(_b,_b)/(tgamma(_b)*_Lambda)*pow(what/_Lambda,_b-1)*exp(-_b*what/_Lambda); #endif    return SF;}

void testSequence(){ int i; long *fib; unsigned char *primeNumbmer;  pdf *dist; dist=normalDistribution(10.0, 5.0, 4, 0.1, 400);  printf("/nSequence/n"); printf("Normal distribution/n"); for (i=0;i<400;i++){   printf("%f/t%f/n",(dist+i)->x,(dist+i)->fx); }    printf("gamma 1 =%f/n",tgamma(4)); printf("gamma 1 =%f/n",gammaFunction(4)); printf("p(a|b) =%f/n",bayes(0.8, 0.01, 0.096)); printf("factorial =%d/n",factorial(5)); printf("double factorial =%d/n",doubleFactorial(5));  fib=fibonacci(20); primeNumbmer=prime(20);  printf("fibonacci/n"); for (i=0;i<20;i++){ printf("%d ",i); }  printf("/n"); for (i=0;i<20;i++){ printf("%ld ",fib[i]); }  printf("/n"); printf("prime/n"); for (i=0;i<20;i++){  if (primeNumbmer[i] ==1)   printf("%d ",i); }   printf("/n");}

  inline var<AutodiffOrder, StrictSmoothness, ValidateIO>    tgamma(const var<AutodiffOrder, StrictSmoothness, ValidateIO>& input) {        if (ValidateIO) validate_input(input.first_val(), "tgamma");          const short partials_order = 3;    const unsigned int n_inputs = 1;        create_node<unary_var_node<AutodiffOrder, partials_order>>(n_inputs);    double val = input.first_val();    double g = tgamma(val);        try {      push_dual_numbers<AutodiffOrder, ValidateIO>(g);    } catch (nomad_error) {      throw nomad_output_value_error("tgamma");    }          push_inputs(input.dual_numbers());        try {      if (AutodiffOrder >= 1) {        double dg = digamma(val);        push_partials<ValidateIO>(g * dg);              if (AutodiffOrder >= 2) {          double tg = trigamma(val);          push_partials<ValidateIO>(g * (dg * dg + tg));                if (AutodiffOrder >= 3)            push_partials<ValidateIO>(g * (dg * dg * dg + 3.0 * dg * tg + quadrigamma(val)));        }      }    } catch (nomad_error) {      throw nomad_output_partial_error("tgamma");    }          return var<AutodiffOrder, StrictSmoothness, ValidateIO>(next_node_idx_ - 1);      }

static doubleneg_lgam(double x){	int xi;	double y, z, zero = 0.0;	/* avoid destructive cancellation as much as possible */	if (x > -170) {		xi = x;		if (xi == x)			if (_IEEE)				return(one/zero);			else				return(infnan(ERANGE));		y = tgamma(x);		if (y < 0) {			y = -y;			signgam = -1;		}		return (log(y));	}	z = floor(x + .5);	if (z == x) {		/* convention: G(-(integer)) -> +Inf */		if (_IEEE)			return (one/zero);		else			return (infnan(ERANGE));	}	y = .5*ceil(x);	if (y == ceil(y))		signgam = -1;	x = -x;	z = fabs(x + z);	/* 0 < z <= .5 */	if (z < .25)		z = sin(M_PI*z);	else		z = cos(M_PI*(0.5-z));	z = log(M_PI/(z*x));	y = large_lgam(x);	return (z - y);}

Status apply_unary_operator(const Operator *operator, Stack **operands){    double x = pop_double(operands);    switch (operator->symbol)    {        case '+':            break;        case '-':            x = -x;            break;        case '!':            x = tgamma(x + 1);            break;        default:            return ERROR_UNRECOGNIZED;    }    push_double(x, operands);    return OK;}