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

rowvec GridBdry::interp_normal(cube u) {      //TODO    rowvec temp = zeros<rowvec>(Gamma.n_rows);    for(unsigned int i = 0; i < Gamma.n_rows; i++) {        temp(i) = u(Gamma(i, 0), Gamma(i, 1), 0)*n(i, 0) + u(Gamma(i, 0), Gamma(i, 1), 1)*n(i, 1);    }    return temp;}

void ThreePhaseDecoder::makePhase() {	int n = width * height;	float i1, i2, i3;	for (int i = 0; i < n; i++) {		i1 = (float) graySequence[0][i];		i2 = (float) graySequence[1][i];		i3 = (float) graySequence[2][i];		#ifdef USE_GAMMA		i1 = Gamma(i1 / 255.0, gamma) * 255.0;		i2 = Gamma(i2 / 255.0, gamma) * 255.0;		i3 = Gamma(i3 / 255.0, gamma) * 255.0;		#endif		range[i] = findRange(i1, i2, i3);		mask[i] = range[i] <= rangeThreshold;		ready[i] = !mask[i];		reflectivity[i] = (byte) ((i1 + i2 + i3) / 3);		if(ready[i])			phase[i] = atan2f(sqrtf(3) * (i1 - i3), 2.f * i2 - i1 - i3) / (float) TWO_PI;	}	#ifdef LINEARIZE_PHASE	if(linearize) {		buildLut();		applyLut();	}	#endif}

// ===========================================================double Gamma(double u){// Only for two particular cases:// u = integer// u = integer + 0.5	int cas;	double G;	int i;	int k;		k=(int)(u+0.5);	if(fabs(k-(int)u)<0.5) cas=1; else cas=2;		switch(cas)	{		case 1: // u integer		if(k==1) return 1;			G=1;			for(i=2;i<k;i++) G=G*i;			return G;				case 2: // u = k +1/2. We use the duplication formula		k=k-1;  			G=sqrt(pi)*pow(2,1-2*k)*Gamma(2*k)/Gamma(k);		return G;		default:		ERROR("tools 79. In Gamma, u %f is neither inter nor half-integer");	}}

 /// Probabiliy density function for F distribution /// The F distribution is the ratio of two chi-square distributions with degrees /// of freedom N1 and N2, respectively, where each chi-square has first been /// divided by its degrees of freedom. /// An F-test (Snedecor and Cochran, 1983) is used to test if the standard /// deviations of two populations are equal. This test can be a two-tailed test or /// a one-tailed test. /// The F hypothesis test is defined as: ///   H0:   s1 = s2     (sN is sigma or std deviation) ///   Ha:   s1 < s2     for a lower one tailed test ///         s1 > s2     for an upper one tailed test ///         s1 != s2    for a two tailed test  /// Test Statistic: F = s1^2/s2^2 where s1^2 and s2^2 are the sample variances. /// The more this ratio deviates from 1, the stronger the evidence for unequal /// population variances. Significance Level is alpha, a probability (0<=alpha<=1). /// The hypothesis that the two standard deviations are equal is rejected if ///    F > PP(alpha,N1-1,N2-1)     for an upper one-tailed test ///    F < PP(1-alpha,N1-1,N2-1)     for a lower one-tailed test ///    F < PP(1-alpha/2,N1-1,N2-1)   for a two-tailed test ///    F > PP(alpha/2,N1-1,N2-1) /// where PP(alpha,k-1,N-1) is the percent point function of the F distribution /// [PPfunc is inverse of the CDF : PP(alpha,N1,N2) == F where alpha=CDF(F,N1,N2)] /// with N1 and N2 degrees of freedom and a significance level of alpha.  ///  /// Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.5 /// @param x probability or significance level of the test, >=0 and < 1 /// @param n1   degrees of freedom of first sample, n1 > 0 /// @param n2   degrees of freedom of second sample, n2 > 0 /// @return         the statistic (a ratio variance1/variance2) at this prob double FDistPDF(double x, int n1, int n2) throw(Exception) {    try {       double dn1(n1),dn2(n2);       double F = Gamma((dn1+dn2)/2.0) / (Gamma(dn1/2.0)*Gamma(dn2/2.0));       F *= ::pow(dn1/dn2,dn1/2.0) * ::pow(x,dn1/2.0 - 1.0);       F /= ::pow(1.0+x*dn1/dn2,(dn1+dn2)/2.0);       return F;    }    catch(Exception& e) { GPSTK_RETHROW(e); } }

// Load B-V conversion parameters from config filevoid StelSkyDrawer::initColorTableFromConfigFile(QSettings* conf){	std::map<float,Vec3f> color_map;	for (float bV=-0.5f;bV<=4.0f;bV+=0.01)	{		char entry[256];		sprintf(entry,"bv_color_%+5.2f",bV);		const QStringList s(conf->value(QString("stars/") + entry).toStringList());		if (!s.isEmpty())		{			Vec3f c;			if (s.size()==1)				c = StelUtils::strToVec3f(s[0]);			else				c =StelUtils::strToVec3f(s);			color_map[bV] = Gamma(eye->getDisplayGamma(),c);		}	}	if (color_map.size() > 1)	{		for (int i=0;i<128;i++)		{			const float bV = StelSkyDrawer::indexToBV(i);			std::map<float,Vec3f>::const_iterator greater(color_map.upper_bound(bV));			if (greater == color_map.begin())			{				colorTable[i] = greater->second;			}			else			{				std::map<float,Vec3f>::const_iterator less(greater);--less;				if (greater == color_map.end())				{					colorTable[i] = less->second;				}				else				{					colorTable[i] = Gamma(1.f/eye->getDisplayGamma(), ((bV-less->first)*greater->second + (greater->first-bV)*less->second) *(1.f/(greater->first-less->first)));				}			}		}	}// 	QString res;// 	for (int i=0;i<128;i++)// 	{// 		res += QString("Vec3f(%1,%2,%3),/n").arg(colorTable[i][0], 0, 'g', 6).arg(colorTable[i][1], 0, 'g', 6).arg(colorTable[i][2], 0, 'g', 6);// 	}// 	qDebug() << res;}

double besselpoly(double a, double lambda, double nu) {  int m, factor=0;  double Sm, relerr, Sol;  double sum=0.0;  /* Special handling for a = 0.0 */  if (a == 0.0) {    if (nu == 0.0) return 1.0/(lambda + 1);    else return 0.0;  }  /* Special handling for negative and integer nu */  if ((nu < 0) && (floor(nu)==nu)) {    nu = -nu;    factor = ((int) nu) % 2;  }      Sm = exp(nu*log(a))/(/*cephes_*/Gamma(nu+1)*(lambda+nu+1));  m = 0;  do {    sum += Sm;    Sol = Sm;    Sm *= -a*a*(lambda+nu+1+2*m)/((nu+m+1)*(m+1)*(lambda+nu+1+2*m+2));    m++;    relerr = fabs((Sm-Sol)/Sm);  } while (relerr > EPS && m < 1000);  if (!factor)    return sum;  else    return -sum;}

task autonomous(){    StartTask(descore);    StartTask(velocitycalculator);    if(SensorValue[J1]) // blue    {        if(SensorValue[J2]) //normal        {            // RAM            ArmWall();            ForwardTillStop(127);        }        else // oppo        {            Gamma();        }    }    else // red    {        if(SensorValue[J2]) // normal        {            Beta();        }        else // oppo        {            Delta();        }    }    StopTask(velocitycalculator);}

void MR::calcMagnetizations() {    for( size_t i = 0; i < G.nrNodes(); i++ ) {        if( props.updates == Properties::UpdateType::FULL ) {            // find indices in nb(i)            sub_nb _nbi( G.nb(i).size() );            _nbi.set();            // calc numerator1 and denominator1            Real sum_even, sum_odd;            sum_subs(i, _nbi, &sum_even, &sum_odd);            Mag[i] = (tanh(theta[i]) * sum_even + sum_odd) / (sum_even + tanh(theta[i]) * sum_odd);        } else if( props.updates == Properties::UpdateType::LINEAR ) {            sub_nb empty( G.nb(i).size() );            Mag[i] = T(i,empty);            for( size_t _l1 = 0; _l1 < G.nb(i).size(); _l1++ )                for( size_t _l2 = _l1 + 1; _l2 < G.nb(i).size(); _l2++ )                    Mag[i] += Gamma(i,_l1,_l2) * tJ[i][_l1] * tJ[i][_l2] * cors[i][_l1][_l2];        }        if( abs( Mag[i] ) > 1.0 )            Mag[i] = (Mag[i] > 0.0) ? 1.0 : -1.0;    }}

double rf_gauss (double d2, const double *a){   /* --- (general.) Gaussian function */    double ma;                    /* temporary buffer for m/a */    if (d2 < 0) {                 /* if to the get normalization factor */        if (a[0] <= 0) return 0;    /* check whether the integral exists */        if (a[0] == 2)              /* use simplified formula for a = 2 */            return pow(2*M_PI, 0.5*d2);        ma = -d2 /a[0];        d2 *= -0.5; /* m/a and m/2 (m = number of dims.) */        return (a[0] *Gamma(d2)) / (pow(2, ma+1) *pow(M_PI, d2) *Gamma(ma));    }                             /* return the normalization factor */    if (a[0] != 2)                /* raise distance to the given power */        d2 = pow(d2, 0.5*a[0]);     /* (note that d2 is already squared) */    return exp(-0.5 *d2);         /* compute Gaussian function */}  /* rf_gauss() */

//___________________________________________________________________________________PMathObj _Constant::IGamma (_PMathObj arg){    if (arg->ObjectClass()!=NUMBER) {        _String errMsg ("A non-numerical argument passed to IGamma(a,x)");        WarnError (errMsg);        return new _Constant (0.0);    }    _Parameter x = ((_Constant*)arg)->theValue, sum=0.0;    if (x>1e25) {        x=1e25;    } else if (x<0) {        _String errMsg ("The domain of x is {x>0} for IGamma (a,x)");        WarnError (errMsg);        return new _Constant (0.0);    } else if (x==0.0) {        return new _Constant (0.0);    }    if (x<=theValue+1) // use the series representation        // IGamma (a,x)=exp(-x) x^a /sum_{n=0}^{/infty} /frac{/Gamma((a)}{/Gamma(a+1+n)} x^n    {        _Parameter term = 1.0/theValue, den = theValue+1;        long count = 0;        while ((fabs(term)>=fabs(sum)*machineEps)&&(count<500)) {            sum+=term;            term*=x/den;            den += 1.0;            count++;        }    } else // use the continue fraction representation        // IGamma (a,x)=exp(-x) x^a 1/x+/1-a/1+/1/x+/2-a/1+/2/x+...    {        _Parameter lastTerm = 0, a0 = 1.0, a1 = x, b0 = 0.0, b1 = 1.0, factor = 1.0, an, ana, anf;        for (long count = 1; count<500; count++) {            an = count;            ana = an - theValue;            a0 = (a1+a0*ana)*factor;            b0 = (b1+b0*ana)*factor;            anf = an*factor;            a1  = x*a0+anf*a1;            b1  = x*b0+anf*b1;            if (a1!=0.0) {                factor=1.0/a1;                sum = b1*factor;                if (fabs(sum-lastTerm)/sum<machineEps) {                    break;                }                lastTerm = sum;            }        }    }    _Constant *result = (_Constant*)Gamma();    result->SetValue(sum*exp(-x+theValue*log(x))/result->theValue);    if (x>theValue+1) {        result->SetValue (1.0-result->theValue);    }    return result;}

tmp<fvMatrix<Type> >laplacian(    GeometricField<Type, fvPatchField, volMesh>& vf){    surfaceScalarField Gamma    (        IOobject        (            "1",            vf.time().constant(),            vf.mesh(),            IOobject::NO_READ        ),        vf.mesh(),        dimensionedScalar("1", dimless, 1.0)    );    return fvm::laplacian    (        Gamma,        vf,        "laplacian(" + vf.name() + ')'    );}

void FinalKick(particle_t *SPH, double dt){	#pragma omp parallel for	for(int i = 0 ; i < N_SPHP ; ++ i){		SPH[i].v = SPH[i].v_h + 0.5 * dt * SPH[i].a    ;		SPH[i].u = SPH[i].u_h + 0.5 * dt * SPH[i].u_dot;		SPH[i].Y = SPH[i].Y_h + 0.5 * dt * SPH[i].m * (Gamma(SPH[i].rho, SPH[i].u, SPH[i].p_smth) - 1.0) * SPH[i].u_dot;	}}

/*void TimeEvolve(particle_t *SPH, double dt){	#pragma omp parallel for	for(int i = 0 ; i < N_SPHP ; ++ i){		SPH[i].r   += SPH[i].v * dt + 0.5 * SPH[i].a * dt * dt;		SPH[i].v   += SPH[i].a * dt;		SPH[i].u   += SPH[i].u_dot * dt;		SPH[i].rho += - SPH[i].rho * SPH[i].div_v * dt;		SPH[i].h   += SPH[i].h * SPH[i].div_v * dt;	}}*/void InitialKick(particle_t *SPH, double dt){	#pragma omp parallel for	for(int i = 0 ; i < N_SPHP ; ++ i){		SPH[i].v_h = SPH[i].v + 0.5 * dt * SPH[i].a    ;		SPH[i].u_h = SPH[i].u + 0.5 * dt * SPH[i].u_dot;//Gammaにsmoothed pを送るのはすごく重要。		SPH[i].Y_h = SPH[i].Y + 0.5 * dt * SPH[i].m * (Gamma(SPH[i].rho, SPH[i].u, SPH[i].p_smth) - 1.0) * SPH[i].u_dot;	}}

void Predict(particle_t *SPH, double dt){	#pragma omp parallel for	for(int i = 0 ; i < N_SPHP ; ++ i){		SPH[i].v += dt * SPH[i].a    ;		SPH[i].u += dt * SPH[i].u_dot;		SPH[i].Y += dt * SPH[i].m * (Gamma(SPH[i].rho, SPH[i].u, SPH[i].p_smth) - 1.0) * SPH[i].u_dot;	}}

/* Returns a sample from Dirichlet(n,a) where n is dimensionality */int Dirichlet(RndState *S, long n, double *a, double *x){	long i;	double tot=0, z;	for (i=0; i<n; i++) { z=Gamma(S,a[i]); tot+=z; x[i]=z; }	for (i=0; i<n; i++) { x[i]/=tot; }	return 1;}

Module::ReturnType FeasibleUpwardPlanarSubgraph::call(	const Graph &G,	GraphCopy &FUPS,	adjEntry &extFaceHandle,	List<edge> &delEdges,	bool multisources){	FUPS = GraphCopy(G);	delEdges.clear();	node s_orig;	hasSingleSource(G, s_orig);	List<edge> nonTreeEdges_orig;	getSpanTree(FUPS, nonTreeEdges_orig, true, multisources);	CombinatorialEmbedding Gamma(FUPS);	nonTreeEdges_orig.permute(); // random order	//insert nonTreeEdges	while (!nonTreeEdges_orig.empty()) {		// make identical copy GC of Fups		//and insert e_orig in GC		GraphCopy GC = FUPS;		edge e_orig = nonTreeEdges_orig.popFrontRet();		//node a = GC.copy(e_orig->source());		//node b = GC.copy(e_orig->target());		GC.newEdge(e_orig);		if (UpwardPlanarity::upwardPlanarEmbed_singleSource(GC)) { //upward embedded the fups and check feasibility			CombinatorialEmbedding Beta(GC);			//choose a arbitrary feasibel ext. face			FaceSinkGraph fsg(Beta, GC.copy(s_orig));			SList<face> ext_faces;			fsg.possibleExternalFaces(ext_faces);			OGDF_ASSERT(!ext_faces.empty());			Beta.setExternalFace(ext_faces.front());			GraphCopy M = GC; // use a identical copy of GC to constrcut the merge graph of GC			adjEntry extFaceHandle_cur = getAdjEntry(Beta, GC.copy(s_orig), Beta.externalFace());			adjEntry adj_orig = GC.original(extFaceHandle_cur->theEdge())->adjSource();			if (constructMergeGraph(M, adj_orig, nonTreeEdges_orig)) {				FUPS = GC;				extFaceHandle = FUPS.copy(GC.original(extFaceHandle_cur->theEdge()))->adjSource();				continue;			}			else {				//Beta is not feasible				delEdges.pushBack(e_orig);			}		}		else {			// not ok, GC is not feasible			delEdges.pushBack(e_orig);		}	}	return Module::retFeasible;}

 void Round (u32 const * const k,u32 * const a,u8 const RC1,u8 const RC2){   a[0] ^= RC1;  Theta(k,a);   a[0] ^= RC2;  Pi1(a);   Gamma(a);   Pi2(a); }  /* Round */

double p_gamma(double a,double x){  int na = floor(a-1.0);  double dna = (a-1.0) - na;  double e = exp(-x/na)*x;  double res = 1.0;  for(int i=1;i<=na;i++){ res *= e/i; }  res *= exp(-dna)*pow(x,dna)/Gamma(dna);  return res;}

   /// Probability density function (PDF) of the Chi-square distribution.   /// The chi-square distribution results when n independent variables with   /// standard normal distributions are squared and summed; x=RSS(variables).   ///   /// A chi-square test (Snedecor and Cochran, 1983) can be used to test if the   /// standard deviation of a population is equal to a specified value. This test   /// can be either a two-sided test or a one-sided test. The two-sided version   /// tests against the alternative that the true standard deviation is either   /// less than or greater than the specified value. The one-sided version only   /// tests in one direction.   /// The chi-square hypothesis test is defined as:   ///   H0:   sigma = sigma0   ///   Ha:   sigma < sigma0    for a lower one-tailed test   ///   sigma > sigma0          for an upper one-tailed test   ///   sigma <>sigma0          for a two-tailed test   ///   Test Statistic:   T = T = (N-1)*(s/sigma0)**2   ///      where N is the sample size and s is the sample standard deviation.   /// The key element of this formula is the ratio s/sigma0 which compares the ratio   /// of the sample standard deviation to the target standard deviation. As this   /// ratio deviates from 1, the more likely is rejection of the null hypothesis.   /// Significance Level:  alpha.   /// Critical Region:  Reject the null hypothesis that the standard deviation   /// is a specified value, sigma0, if   ///   T > chisquare(alpha,N-1)     for an upper one-tailed alternative   ///   T < chisquare(1-alpha,N-1)   for a lower one-tailed alternative   ///   T < chisquare(1-alpha,N-1)   for a two-tailed test or   ///   T < chisquare(1-alpha,N-1)   /// where chi-square(p,N-1) is the critical value or inverseCDF of the chi-square   /// distribution with N-1 degrees of freedom.    ///   /// @param x input statistic, equal to an RSS(); x >= 0   /// @param n    input value for number of degrees of freedom, n > 0   /// @return         probability Chi-square probability (xsq,n)   double ChisqPDF(const double& x, const int& n) throw(Exception)   {      if(x < 0) GPSTK_THROW(Exception("Negative statistic"));      if(n <= 0)         GPSTK_THROW(Exception("Non-positive degrees of freedom"));      try {         double dn(double(n)/2.0);         return ( ::exp(-x/2.0) * ::pow(x,dn-1.0) / (::pow(2.0,dn) * Gamma(dn)) );      }      catch(Exception& e) { GPSTK_RETHROW(e); }   }

/*==================================================================================*/void Round (u32 const * const k,u32 * const a,u8 const RC1,u8 const RC2)/*----------------------------------------------------------------------------------*//* The round function, common to both encryption and decryption/* - Round constants is added to the rightmost byte of the leftmost 32-bit word (=a0)/*==================================================================================*/{   a[0] ^= RC1;  Theta(k,a);   a[0] ^= RC2;  Pi1(a);   Gamma(a);   Pi2(a); }  /* Round */

	void SemiImplicitDiffusionAlgorithm::computeDiffusivity() 	{		double c1(0.), c2(0.);		for (unsigned int i(1); i<m_Nx; ++i) 		{			for (unsigned int j(1); j<m_Ny; ++j) 			{				c1 = m_H->Staggered(i, j);				c2 = staggeredGradSurfNorm(i, j, m_H);				D(i, j) = (Gamma()*c1 + rhogn()*Sl(i, j)) * pow(c1, n() + 1)*pow(c2, n() - 1);			}		}	}

int main() {  double v;  scanf("%lf", &v);  if (v < 0 || v > 10) {    printf("%f is not valid input./n", v);    exit(1);  }  if (fmod(v, 1.0) < 1e-8) {    printf("%d/n", fact(v));  } else {    printf("%.8f/n", Gamma(v));  }}

long Poisson(RndState *S, double lambda){	long r;	if (lambda>=15) {		double m=floor(lambda*7/8);		double x=Gamma(S,m);		if (x>lambda) r=Binomial(S,lambda/x,m-1);		else          r=m+Poisson(S,lambda-x);	} else {		double p, elambda = exp(-lambda);		for (p=1, r=-1; p>elambda; p*=Uniform(S)) r++; 	}	return r;}

int main (int argc, char *argv[]){                               /* --- main function */  double x;                     /* argument */  if (argc != 2) {              /* if wrong number of arguments given */    printf("usage: %s x/n", argv[0]);    printf("compute (logarithm of) Gamma function/n");    return 0;                   /* print a usage message */  }                             /* and abort the program */  x = atof(argv[1]);            /* get argument */  if (x <= 0) { printf("%s: x must be > 0/n", argv[0]); return -1; }  printf("   Gamma(%.16g)  = % .20g/n", x, Gamma(x));  printf("ln(Gamma(%.16g)) = % .20g/n", x, logGamma(x));  return 0;                     /* compute and print Gamma function */}  /* main() */

double Gamma(double x) {    double pi=3.14159265;     if( x==0.5 )	 return sqrt(pi);  if( x== 1 )	 return 1.;  if( x==1.5 )	 return sqrt(pi)/2.;  if( x==2. )	 return 1.;  if( x==2.5 )	 return 3*sqrt(pi)/4.;  if( x==3 )	 return 2.;  if( x==3.5 )	 return 15*sqrt(pi)/8.;  if( x==4 )	 return 6.;  if( x > 3 )    return((x-1)*Gamma(x-1)); // impact::LOGGER_WRITE_ALONE(impact::msLogger::FATAL, "Pb with the Gamma function, integer or half integer positive number must be providen"   //                                 , "double Gamma(double x)");  return -1; }

double LogGamma(    double x    // x must be positive){	if (x <= 0.0)	{		std::stringstream os;        os << "Invalid input argument " << x <<  ". Argument must be positive.";        throw std::invalid_argument( os.str() ); 	}    if (x < 12.0)    {        return log(fabs(Gamma(x)));    }	// Abramowitz and Stegun 6.1.41    // Asymptotic series should be good to at least 11 or 12 figures    // For error analysis, see Whittiker and Watson    // A Course in Modern Analysis (1927), page 252    static const double c[8] =    {		 1.0/12.0,		-1.0/360.0,		1.0/1260.0,		-1.0/1680.0,		1.0/1188.0,		-691.0/360360.0,		1.0/156.0,		-3617.0/122400.0    };    double z = 1.0/(x*x);    double sum = c[7];    for (int i=6; i >= 0; i--)    {        sum *= z;        sum += c[i];    }    double series = sum/x;    static const double halfLogTwoPi = 0.91893853320467274178032973640562;    double logGamma = (x - 0.5)*log(x) - x + halfLogTwoPi + series;    	return logGamma;}

int main(int argc, char* argv[ ]) {	double v, w, x, y, z;	double timer = omp_get_wtime();	int thread_count = 1;	if (argc>1)		thread_count = strtol(argv[1], NULL, 10);#	pragma omp parallel num_threads(thread_count)  	{#		pragma omp sections		{#			pragma omp section			{				v = Alpha( );				printf("v = %lf/n", v);			}#			pragma omp section			{				w = Beta( );				printf("w = %lf/n", w);			}		}#		pragma omp sections		{#			pragma omp section			{				x = Gamma(v, w);				printf("x = %lf/n", x);			}#			pragma omp section			{				y = Delta( );				printf("y = %lf/n", y);			}		}	}	z = Epsilon(x, y);	printf("z = %lf/n", z);	timer = omp_get_wtime() - timer;	printf("timer = %lf/n", timer);	return 0;}

SimplePropertySet<string, double> propertylist() {	SimplePropertySet<string, double> result;	result.add (Property<string, double> ("Option Value", Price() ) );	result.add (Property<string, double> ("Delta",Delta() ) );	result.add (Property<string, double> ("Gamma",Gamma() ) );	result.add (Property<string, double> ("Vega",Vega() ) );	result.add (Property<string, double> ("Vega",Theta() ) );	result.add (Property<string, double> ("Rho",Rho() ) );	result.add (Property<string, double> ("Cost of Carry",Coc() ) );										// Cost of carry		cout << "counbt " << result.Count();	return result;}

double rf_cauchy (double d2, const double *a){   /* --- (generalized) Cauchy function */    double ma;                    /* temporary buffer for m/a */    assert(a);                    /* check the function arguments */    if (d2 < 0) {                 /* if to the get normalization factor */        if ((a[0] <= -d2) || (a[1] <= 0))            return 0;                 /* check whether the integral exists */        ma = -d2 /a[0];        d2 *= -0.5; /* m/a and m/2 (m = number of dims.) */        return (a[0] *Gamma(d2) *sin(ma *M_PI))               / (2 *pow(M_PI, d2+1) *pow(a[1], ma-1));    }                             /* return the normalization factor */    if (a[0] != 2)                /* raise distance to the given power */        d2 = pow(d2, 0.5*a[0]);     /* (note that d2 is already squared) */    d2 += a[1];                   /* add offset to distance */    return (d2 > MINDENOM) ? 1/d2 : 1/MINDENOM;}  /* rf_cauchy() */            /* compute Cauchy function */

int getCellScore (ProcessData * pData, ScoringData * sData, WavesData * wData, MOATypeShape * cellIndex, MOATypeElmVal * score, int * inSearchSpace, int NeighborSearch, MOATypeInd NeighbIndex) {    int ret = 0;    MOATypeDimn k;    MOATypeInd NeighbFlatIndex;    /*Check if cellIndex is found in the current scoring partition*/    if ((NeighborSearch == 1) && (IsCellInPart(cellIndex, sData->p_index, sData->seqNum, sData->seqLen, pData->partitionSize) == 0) &&         (getLocalIndex (cellIndex, sData->p_index, sData->seqNum, sData->seqLen, pData->partitionSize, &sData->neighbor) == 0)) {        NeighbFlatIndex = Gamma(sData->neighbor, sData->msaAlgn->dimn, sData->msaAlgn->shape,  sData->msaAlgn->dimn, 1);       (*score) = sData->msaAlgn->elements[NeighbFlatIndex].val;       if (sData->msaAlgn->elements[NeighbFlatIndex].prev != NULL && sData->msaAlgn->elements[NeighbFlatIndex].prev_ub > 0 &&            sData->NghbMOA != NULL && NeighbIndex >= 0 && NeighbIndex < sData->NghbMOA->elements_ub) {           sData->NghbMOA->elements[NeighbIndex].prev = mmalloc(sizeof *sData->NghbMOA->elements[NeighbIndex].prev);           sData->NghbMOA->elements[NeighbIndex].prev_ub = 1;           sData->NghbMOA->elements[NeighbIndex].prev[0] = mmalloc(sData->seqNum * sizeof *sData->NghbMOA->elements[NeighbIndex].prev[0]);           for (k=0;k<sData->seqNum;k++)                sData->NghbMOA->elements[NeighbIndex].prev[0][k] = sData->msaAlgn->elements[NeighbFlatIndex].prev[0][k];       }    }    else {        /*check if neighbor's partition is included in search space*/        MOATypeShape * partIndex = mmalloc (pData->seqNum * sizeof *partIndex);        if  (getPartitionIndex (cellIndex, pData->seqNum, pData->seqLen, wData->partitionSize, &partIndex) == 0) {            if ((*inSearchSpace) = isPartInSearchSpace(partIndex, wData) == 0) {                long waveNo, partNo;                getPartitionPosition (wData, partIndex, &waveNo, &partNo);                if (partNo >= 0) {                    if (myProcid == getProcID (wData, waveNo, partNo)) {                                                /*Check if Neighbor is found in other local partitions OCout Buffer*/                        if(checkPrevPartitions(pData, cellIndex, score) != 0) {                            /*average the neighboring (up to 2 strides) cell scores*/                            (*score) = averageNeighborsScore(pData, sData, wData, cellIndex);                        }                    }                    /*Check if Neighbor is already received from other processors in OCin Buffer*/                    else if (checkRecvOC(pData, wData, cellIndex, score, 0) != 0)                            ret = -1;                }            }        }        free (partIndex);    }    return ret;}