void test_erfc()
{
    static_assert((std::is_same<decltype(erfc((double)0)), double>::value), "");
    static_assert((std::is_same<decltype(erfcf(0)), float>::value), "");
    static_assert((std::is_same<decltype(erfcl(0)), long double>::value), "");
    assert(erfc(0) == 1);
}

  double JointProb::probparam(double parameter, int thecoord, double maxd0sig   ) const
  {
    
    double prob = 0;
    double resolutionparameters[3][5];
    if (_ResolutionParameterRphi.size() != 5 || _ResolutionParameterZ.size() != 5 || _ResolutionParameter3D.size() != 5) 
    std::cerr << "Warning jointprob.cpp:229 Parameters of wrong length" << std::endl;
    for( int iii=0; iii<5; iii++)
    {
	resolutionparameters[0][iii] = _ResolutionParameterRphi[iii];
	resolutionparameters[1][iii] = _ResolutionParameterZ[iii];
	resolutionparameters[2][iii] = _ResolutionParameter3D[iii];
    }
    
    // The if statement takes into account a different parametrization for different coordinates. 
    if ( thecoord < 2 )
      {
	// part one is the gaussian part
	// to understand this part better one should look at the meaning of the complementary error function
	prob = erfc( parameter / (sqrt( double(2) ) * resolutionparameters[thecoord][0] ))
	  -erfc( maxd0sig / ( sqrt( double(2) ) *  resolutionparameters[thecoord][0] ) );
	
	// part 2 is the added exponential tails
      
      prob += resolutionparameters[thecoord][1]* ( exp(- resolutionparameters[thecoord][2] * parameter )
						   -exp(- resolutionparameters[thecoord][2] * maxd0sig ) )
	+resolutionparameters[thecoord][3]* ( exp(- resolutionparameters[thecoord][4] * parameter )
					      -exp(- resolutionparameters[thecoord][4] * maxd0sig ) );
      
      
      }
    else
      {
	
	//in this view the gaussian part is just squared. 
	
	prob = exp(- ( parameter * parameter ) /  ( resolutionparameters[thecoord][0] * resolutionparameters[thecoord][0] * double ( 2 ) ) ) 
	  -exp(- ( maxd0sig * maxd0sig )  /  ( resolutionparameters[thecoord][0] * resolutionparameters[thecoord][0] * double ( 2 ) ) );

	prob +=  resolutionparameters[thecoord][1]* ( ( 1 + resolutionparameters[thecoord][2] * parameter )
						      *exp ( - resolutionparameters[thecoord][2] * parameter )
						      - ( 1 + resolutionparameters[thecoord][2] * maxd0sig )
						      *exp ( - resolutionparameters[thecoord][2] * maxd0sig ))
	  +resolutionparameters[thecoord][3]* ( ( 1 + resolutionparameters[thecoord][4] * parameter )
						*exp ( - resolutionparameters[thecoord][4] * parameter )
						- ( 1 + resolutionparameters[thecoord][4] * maxd0sig )
						*exp ( - resolutionparameters[thecoord][4] * maxd0sig ));
	
      }

    return prob;
    
  }

	double operator()(const double time, const Vect3d& r) {
		const double x = r[0];
		const double alpha = 1.0;
		const double beta = sqrt(alpha/D);
		const double exact = (mu/(D*beta)) * (
				exp(-beta*(1.0-x))
				- 0.5*exp(-beta*(1.0-x))*erfc((2.0*beta*D*time-(1.0-x))/sqrt(4.0*D*time))
				- 0.5*exp(beta*(1.0-x))*erfc((2.0*beta*D*time+(1.0-x))/sqrt(4.0*D*time))
		);

		return exact;
	}

	double energyAndGrad(std::vector<Atom>& atoms) const
	{	double eta = sqrt(0.5)/sigma, etaSq=eta*eta;
		double sigmaSq = sigma * sigma;
		double detR = fabs(det(R)); //cell volume
		//Position independent terms:
		double Ztot = 0., ZsqTot = 0.;
		for(const Atom& a: atoms)
		{	Ztot += a.Z;
			ZsqTot += a.Z * a.Z;
		}
		double E
			= 0.5 * 4*M_PI * Ztot*Ztot * (-0.5*sigmaSq) / detR //G=0 correction
			- 0.5 * ZsqTot * eta * (2./sqrt(M_PI)); //Self-energy correction
		//Reduce positions to first centered unit cell:
		for(Atom& a: atoms)
			for(int k=0; k<3; k++)
				a.pos[k] -= floor(0.5 + a.pos[k]);
		//Real space sum:
		vector3<int> iR; //integer cell number
		for(const Atom& a2: atoms)
			for(Atom& a1: atoms)
				for(iR[0]=-Nreal[0]; iR[0]<=Nreal[0]; iR[0]++)
					for(iR[1]=-Nreal[1]; iR[1]<=Nreal[1]; iR[1]++)
						for(iR[2]=-Nreal[2]; iR[2]<=Nreal[2]; iR[2]++)
						{	vector3<> x = iR + (a1.pos - a2.pos);
							double rSq = RTR.metric_length_squared(x);
							if(!rSq) continue; //exclude self-interaction
							double r = sqrt(rSq);
							E += 0.5 * a1.Z * a2.Z * erfc(eta*r)/r;
							a1.force += (RTR * x) *
								(a1.Z * a2.Z * (erfc(eta*r)/r + (2./sqrt(M_PI))*eta*exp(-etaSq*rSq))/rSq);
						}
		//Reciprocal space sum:
		vector3<int> iG; //integer reciprocal cell number
		for(iG[0]=-Nrecip[0]; iG[0]<=Nrecip[0]; iG[0]++)
			for(iG[1]=-Nrecip[1]; iG[1]<=Nrecip[1]; iG[1]++)
				for(iG[2]=-Nrecip[2]; iG[2]<=Nrecip[2]; iG[2]++)
				{	double Gsq = GGT.metric_length_squared(iG);
					if(!Gsq) continue; //skip G=0
					//Compute structure factor:
					complex SG = 0.;
					for(const Atom& a: atoms)
						SG += a.Z * cis(-2*M_PI*dot(iG,a.pos));
					//Accumulate energy:
					double eG = 4*M_PI * exp(-0.5*sigmaSq*Gsq)/(Gsq * detR);
					E += 0.5 * eG * SG.norm();
					//Accumulate forces:
					for(Atom& a: atoms)
						a.force -= (eG * a.Z * 2*M_PI * (SG.conj() * cis(-2*M_PI*dot(iG,a.pos))).imag()) * iG;
				}
		return E;
	}

double Hypermet::eval_step_tail(double x) {
  if (width_.val == 0)
    return 0;

  double xc = x - center_.val;

  double step = step_amplitude.val * erfc( xc/width_.val );

  double tail = 0;
  double lexp = exp(pow(0.5*width_.val/tail_slope.val, 2) + xc/tail_slope.val);
  if ((tail_slope.val != 0) && !isinf(lexp))
    tail = tail_amplitude.val * lexp * erfc( 0.5*width_.val/tail_slope.val + xc/width_.val);

  return height_.val * 0.5 * (step + tail);
}

int main() {
    double x = 1.0;
    double y = 1.0;
    int i = 1;
    acosh(x);
    asinh(x);
    atanh(x);
    cbrt(x);
    expm1(x);
    erf(x);
    erfc(x);
    isnan(x);
    j0(x);
    j1(x);
    jn(i,x);
    ilogb(x);
    logb(x);
    log1p(x);
    rint(x);
    y0(x);
    y1(x);
    yn(i,x);
#   ifdef _THREAD_SAFE
    gamma_r(x,&i);
    lgamma_r(x,&i);
#   else
    gamma(x);
    lgamma(x);
#   endif
    hypot(x,y);
    nextafter(x,y);
    remainder(x,y);
    scalb(x,y);
    return 0;
}

/** Calculate real space contribution of coulomb pair forces.
    If NPT is compiled in, it returns the energy, which is needed for NPT. */
inline double p3m_add_pair_force(double chgfac, double *d,double dist2,double dist,double force[3])
{
  int j;
  double fac1,fac2, adist, erfc_part_ri;
  if(dist < p3m.params.r_cut) {
    if (dist > 0.0){		//Vincent
      adist = p3m.params.alpha * dist;
#if USE_ERFC_APPROXIMATION
      erfc_part_ri = AS_erfc_part(adist) / dist;
      fac1 = coulomb.prefactor * chgfac  * exp(-adist*adist);
      fac2 = fac1 * (erfc_part_ri + 2.0*p3m.params.alpha*wupii) / dist2;
#else
      erfc_part_ri = erfc(adist) / dist;
      fac1 = coulomb.prefactor * chgfac;
      fac2 = fac1 * (erfc_part_ri + 2.0*p3m.params.alpha*wupii*exp(-adist*adist)) / dist2;
#endif
      for(j=0;j<3;j++)
	force[j] += fac2 * d[j];
      ESR_TRACE(fprintf(stderr,"%d: RSE: Pair dist=%.3f: force (%.3e,%.3e,%.3e)\n",this_node,
			dist,fac2*d[0],fac2*d[1],fac2*d[2]));
#ifdef NPT
      return fac1 * erfc_part_ri;
#endif
    }
  }
  return 0.0;
}

double Frequency(int ** bits, int tailleCollection, int tailleMot)
{
    int Sn = 0;
    double Sobs = 0;
    for(int j = 0; j < tailleCollection; j++)
    {
        for(int i = 0; i < tailleMot; i ++)
        {
            if(bits[j][i] == 0)
            {
                Sn--;
            }
            else
            {
                Sn++;
            }
        }
    }
    Sobs = fabs(Sn)/ (double) sqrt(tailleCollection*tailleMot);

    double result = erfc((double)(Sobs/sqrt(2.0)));
    printf("Sn : %d\n", Sn);
    printf("Sobs : %lf\n", Sobs);
    printf("Pval : %lf\n", result);

    return result;
}

static double
l_erfc(char *nm)
{
    extern double  erfc();

    return(erfc(argument(1)));
}

static void
yyerf(void)
{
	double d;

	p1 = pop();

	if (isdouble(p1)) {
		d = 1.0 - erfc(p1->u.d);
		push_double(d);
		return;
	}

	if (isnegativeterm(p1)) {
		push_symbol(ERF);
		push(p1);
		negate();
		list(2);
		negate();
		return;
	}
	
	push_symbol(ERF);
	push(p1);
	list(2);
	return;
}

void Foam::equationReader::evalDimsErfcDimCheck
(
    const equationReader * eqnReader,
    const label index,
    const label i,
    const label storageOffset,
    label& storeIndex,
    dimensionSet& xDims,
    dimensionSet sourceDims
) const
{
    if
    (
        !xDims.dimensionless() && dimensionSet::debug
    )
    {
        WarningIn("equationReader::evalDimsErfcDimCheck")
            << "Dimension error thrown for operation ["
            << equationOperation::opName
            (
                operator[](index)[i].operation()
            )
            << "] in equation " << operator[](index).name()
            << ", given by:" << token::NL << token::TAB
            << operator[](index).rawText();
    }
    dimensionedScalar ds("temp", xDims, 0.0);
    xDims.reset(erfc(ds).dimensions());
    operator[](index)[i].assignOpDimsFunction
    (
        &Foam::equationReader::evalDimsErfc
    );
}

double lrchisq(const double *par)
{
  // chisq from likelihood ratio
  // for more information see Baker and Cousins pg. 439 and Appendix
  double yi=0.;        // model
  double ni=0.;        // experiment
  double lrchisq = 0.; // likelihood ratio chisq
  int i=0;

  for (i=startCh[spCurrent]; i<=endCh[spCurrent]; i++)
    {
      // events in ith bin
      ni = expCurrent[i]; 
      
      // calculate model in the ith bin
      yi = par[0]*simCurrent[i] + par[1] + par[2]*erfc(((double)i-sfc[spCurrent])/sfw[spCurrent]);
      
      // evaluate chisq given input parameters
      if(ni > 0.)
	lrchisq += ( yi - ni + ni*log(ni/yi) );
      else
	lrchisq += yi; // the log(0) case
    }
  lrchisq *= 2.;
  return lrchisq;
}

double ntail(double zval) 
/** normal distribution tail area 
 uses erfc 
*/

{
  double pi,  t ;
  double p, q, d ;

  if (zval == 0.0) return 0.5 ;
  if (zval<0.0) return (1.0 - ntail(-zval)) ;
  if (zval<ZLIM) {
   t = zval/sqrt(2.0) ;
   q = erfc(t)/2.0 ;
   return q ;
  }

  pi = 2.0*acos(0.0)  ;

  t = exp(-0.5*zval*zval) ;
  t /= (sqrt(2.0*pi) * zval) ;

  return t ;

}

double erf(double x)
{
    //
    // Computation of the error function erf(x).
    //
    return (1-erfc(x));
}

/*============================================================================*/
double ighmm_rand_get_1overa (double x, double mean, double u)
{
  /* Calulates 1/a(x, mean, u), with a = the integral from x til \infty over
     the Gauss density function */
# define CUR_PROC "ighmm_rand_get_1overa"

  double erfc_value;

  if (u <= 0.0) {
    GHMM_LOG(LCONVERTED, "u <= 0.0 not allowed\n");
    goto STOP;
  }

#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
  erfc_value = erfc ((x - mean) / sqrt (u * 2));
#else
  erfc_value = ighmm_erfc ((x - mean) / sqrt (u * 2));
#endif

  if (erfc_value <= DBL_MIN) {
    ighmm_mes (MES_WIN, "a ~= 0.0 critical! (mue = %.2f, u =%.2f)\n", mean, u);
    return (erfc_value);
  }
  else
    return (2.0 / erfc_value);

STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* ighmm_rand_get_1overa */

/* cumalative distribution function of a-truncated N(mean, u) */
double ighmm_rand_normal_right_cdf (double x, double mean, double u, double a)
{
# define CUR_PROC "ighmm_rand_normal_right_cdf"

  if (x <= a)
    return (0.0);
  if (u <= a) {
    GHMM_LOG(LCONVERTED, "u <= a not allowed\n");
    goto STOP;
  }
#if defined(__STDC_VERSION__) && (__STDC_VERSION__ >= 199901L)
  /*
     Function: int erfc (double x, gsl_sf_result * result) 
     These routines compute the complementary error function
     erfc(x) = 1 - erf(x) = 2/\sqrt(\pi) \int_x^\infty \exp(-t^2). 
   */
  return 1.0 + (erf ((x - mean) / sqrt (u * 2)) -
                1.0) / erfc ((a - mean) / sqrt (u * 2));
#else
  return 1.0 + (ighmm_erf ((x - mean) / sqrt (u * 2)) -
                1.0) / ighmm_erfc ((a - mean) / sqrt (u * 2));
#endif /* Check for ISO C99 */
STOP:
  return (-1.0);
# undef CUR_PROC
}                               /* double ighmm_rand_normal_cdf */

double Runs(int n, unsigned char *epsilon)
{
	int		S, k;
	double	pi, V, erfc_arg, p_value;

	S = 0;
	for ( k=0; k<n; k++ )
		if ( epsilon[k] )
			S++;
	pi = (double)S / (double)n;

	if ( fabs(pi - 0.5) > (2.0 / sqrt(n)) ) {
		//PI ESTIMATOR CRITERIA NOT MET!
		p_value = -1.0;
	}
	else {

		V = 1;
		for ( k=1; k<n; k++ )
			if ( epsilon[k] != epsilon[k-1] )
				V++;
	
		erfc_arg = fabs(V - 2.0 * n * pi * (1-pi)) / (2.0 * pi * (1-pi) * sqrt(2*n));
		p_value = erfc(erfc_arg);
	}

  return p_value ;
}

double ndtr(double a)
{
double x, y, z;

if (isnan(a)) {
  mtherr("ndtr", DOMAIN);
  return (NAN);
}

x = a * SQRTH;
z = fabs(x);

if( z < SQRTH )
	y = 0.5 + 0.5 * erf(x);

else
	{
	y = 0.5 * erfc(z);

	if( x > 0 )
		y = 1.0 - y;
	}

return(y);
}

static void add_top_tbl(struct All_variables *E, double age_in_myrs, double mantle_temp)
{
    int m, i, j, k, node;
    int nox, noy, noz;
    double r1, dT, tmp;

    nox = E->lmesh.nox;
    noy = E->lmesh.noy;
    noz = E->lmesh.noz;

    dT = (mantle_temp - E->control.TBCtopval);
    tmp = 0.5 / sqrt(age_in_myrs / E->data.scalet);

    fprintf(stderr, "%e %e\n", dT, tmp);
    for(m=1; m<=E->sphere.caps_per_proc; m++)
        for(i=1; i<=noy; i++)
            for(j=1; j<=nox;j ++)
                for(k=1; k<=noz; k++) {
                    node = k + (j-1)*noz + (i-1)*nox*noz;
                    r1 = E->sx[m][3][node];
                    E->T[m][node] -= dT * erfc(tmp * (E->sphere.ro - r1));
                }

    return;
}

ATF_TC_BODY(erfc_inf_neg, tc)
{
	const double x = -1.0L / 0.0L;

	if (erfc(x) != 2.0)
		atf_tc_fail_nonfatal("erfc(-Inf) != 2.0");
}