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 file
void 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 */