Real State::ek(const _3Vec& X) const {
#ifdef USE_LZ
	return 0.5 * (sx() * sx() + sy() * sy() / (X[0] * X[0] + X[1] * X[1]) + sz() * sz()) / rho();
#else
	return 0.5 * (sx() * sx() + sy() * sy()+ sz() * sz()) / rho();
#endif
}

      std::vector<Botan::BigInt> factorize(const Botan::BigInt& n_in,
                                           Botan::RandomNumberGenerator& rng)
         {
         Botan::BigInt n = n_in;
         std::vector<Botan::BigInt> factors = remove_small_factors(n);

         while(n != 1)
            {
            if(Botan::is_prime(n, rng))
               {
               factors.push_back(n);
               break;
               }

            Botan::BigInt a_factor = 0;
            while(a_factor == 0)
               {
               a_factor = rho(n, rng);
               }

            std::vector<Botan::BigInt> rho_factored = factorize(a_factor, rng);
            for(size_t j = 0; j != rho_factored.size(); j++)
               {
               factors.push_back(rho_factored[j]);
               }

            n /= a_factor;
            }

         return factors;
         }

 void BatesModel::generateArguments() {
     process_.reset(new BatesProcess(
          process_->riskFreeRate(), process_->dividendYield(),
          process_->s0(), v0(), 
          kappa(), theta(), sigma(), rho(),
          lambda(), nu(), delta()));
 }

  ProbabilityDistribution* 
  GaussianProcessNormal::prediction(const vectord &query)
  {
    double kq = (*mKernel)(query, query);;
    vectord kn = computeCrossCorrelation(query);
    vectord phi = mMean->getFeatures(query);
  
    vectord v(kn);
    inplace_solve(mL,v,ublas::lower_tag());

    vectord rq = phi - prod(v,mKF);

    vectord rho(rq);
    inplace_solve(mD,rho,ublas::lower_tag());
    
    double yPred = inner_prod(phi,mWMap) + inner_prod(v,mVf);
    double sPred = sqrt( mSigma * (kq - inner_prod(v,v) 
			        + inner_prod(rho,rho)));

    if ((boost::math::isnan(yPred)) || (boost::math::isnan(sPred)))
      {
	FILE_LOG(logERROR) << "Error in prediction. NaN found.";
	exit(EXIT_FAILURE);
      }
					

    d_->setMeanAndStd(yPred,sPred);
    return d_;
  }

void KeccakPermutationOnWords(UINT64 *state)
{
    unsigned int i;

    //displayStateAsWords(3, "Same, as words", state);

    for(i=0; i<nrRounds; i++) {
        //displayRoundNumber(3, i);

        theta(state);
        //displayStateAsWords(3, "After theta", state);

        rho(state);
        //displayStateAsWords(3, "After rho", state);

        pi(state);
        //displayStateAsWords(3, "After pi", state);

        chi(state);
        //displayStateAsWords(3, "After chi", state);

        iota(state, i);
        //displayStateAsWords(3, "After iota", state);
    }
}

/**
 * Calculate "surface energy" by integrating the binding energy with respect to
 * moisture content at constant temperature. Since binding energy has units of
 * J/mol, the moisture content is converted to units of moles per unit volume.
 * @param Xdb Dry basis moisture content [-]
 * @param T Temperature [K]
 * @returns [J/m^3]
 */
double Esurf(double Xdb, double T)
{
    double Mw = 0.0180153, /* kg/mol */
           rhos,
           Xmin = 1e-9, /* Since binding energy is undefined at Xdb = 0, use a
                           larger minimum moisture content instead. */
           dX, /* Size of each integration slice */
           X, /* Current moisture content for calculations */
           E = 0; /* Set the energy to zero initially. */
    int nslice = 100, /* Number of slices to integrate with */
        i; /* Loop index */
    oswin *d;
    choi_okos *co;

    /* Get the property data we need */
    d = CreateOswinData();
    co = CreateChoiOkos(PASTACOMP);
    
    /* Calculate slice width */
    dX = (Xdb-Xmin)/nslice;
    /* Integrate from Xmin to Xdb */
    for(i=0; i<nslice; i++) {
        X = Xmin+i*dX;
        E += BindingEnergyOswin(d, X, T) * dX;
    }

    /* Calculate solid density for unit conversion */
    rhos = rho(co, T);
    DestroyChoiOkos(co);
    DestroyOswinData(d);
    /* Ensure the units are consitent and return the result. */
    return E*Xdb*rhos/Mw;
}

void rho(LL n)
{
	int flag=0;
	while (true)
	{
		LL c=rand()%n+1;
		LL x=rand()%n+1;
		LL y=x;
		if (isPrime(n))
		{
			num++;
			a[num]=n;
			return;
		}
		LL k=2;
		LL i=0;
		while (true)
		{
			i++;
			LL d=gcd(n+y-x,n);
			if (d>1&&d<n)
			{
				flag=1;
				if (isPrime(d))
				{
					num++;
					a[num]=d;
				}
				else rho(d);
				d=n/d;
				if (isPrime(d))
				{
					num++;
					a[num]=d;
				}
				else rho(d);
			}
			if (i==k)
			{
				y=x;
				k*=2;
			}
			x=(pow_mod(x,x,n)+n-c)%n;
			if (y==x) break;
		}
	}
}

Real
BrineFluidProperties::e(Real pressure, Real temperature, Real xnacl) const
{
  Real enthalpy = h(pressure, temperature, xnacl);
  Real density = rho(pressure, temperature, xnacl);

  return enthalpy - pressure / density;
}

int diffrpmw(double t, const double rpmw[], double drpmw[], void *params)
{
  drpmw[0] = 1.0;
  drpmw[1] = rpmw[2] / (rpmw[0]*rpmw[0]);
  drpmw[2] = FOUR_PI * (rpmw[0]*rpmw[0]) * rho(rpmw[0], rpmw[1]);
  drpmw[3] = 0.5 * rpmw[1] * drpmw[2];
  return (GSL_SUCCESS);
}

void oneEqEddy::updateSubGridScaleFields()
{
    muSgs_ = ck_*rho()*sqrt(k_)*delta();
    muSgs_.correctBoundaryConditions();

    alphaSgs_ = muSgs_/Prt_;
    alphaSgs_.correctBoundaryConditions();
}

void dynOneEqEddy::updateSubGridScaleFields(const volSymmTensorField& D)
{
    muSgs_ = ck_(D)*rho()*sqrt(k_)*delta();
    muSgs_.correctBoundaryConditions();

    alphaSgs_ = muSgs_/Prt_;
    alphaSgs_.correctBoundaryConditions();
}

void tractionDisplacementFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const dictionary& mechanicalProperties =
        db().lookupObject<IOdictionary>("mechanicalProperties");

    const dictionary& thermalProperties =
        db().lookupObject<IOdictionary>("thermalProperties");

    dimensionedScalar rho(mechanicalProperties.lookup("rho"));
    dimensionedScalar rhoE(mechanicalProperties.lookup("E"));
    dimensionedScalar nu(mechanicalProperties.lookup("nu"));

    dimensionedScalar E = rhoE/rho;
    dimensionedScalar mu = E/(2.0*(1.0 + nu));
    dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
    dimensionedScalar threeK = E/(1.0 - 2.0*nu);

    Switch planeStress(mechanicalProperties.lookup("planeStress"));

    if (planeStress)
    {
        lambda = nu*E/((1.0 + nu)*(1.0 - nu));
        threeK = E/(1.0 - nu);
    }

    scalar twoMuLambda = (2*mu + lambda).value();

    vectorField n(patch().nf());

    const fvPatchField<symmTensor>& sigmaD =
        patch().lookupPatchField<volSymmTensorField, symmTensor>("sigmaD");

    gradient() =
    (
        (traction_ + pressure_*n)/rho.value()
      + twoMuLambda*fvPatchField<vector>::snGrad() - (n & sigmaD)
    )/twoMuLambda;

    Switch thermalStress(thermalProperties.lookup("thermalStress"));

    if (thermalStress)
    {
        dimensionedScalar alpha(thermalProperties.lookup("alpha"));
        dimensionedScalar threeKalpha = threeK*alpha;

        const fvPatchField<scalar>& T =
            patch().lookupPatchField<volScalarField, scalar>("T");

        gradient() += n*threeKalpha.value()*T/twoMuLambda;
    }

    fixedGradientFvPatchVectorField::updateCoeffs();
}

// Update the coefficients associated with the patch field
void tractionDisplacementFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const dictionary& mechanicalProperties =
        db().lookupObject<IOdictionary>("mechanicalProperties");

    const dictionary& thermalProperties =
        db().lookupObject<IOdictionary>("thermalProperties");

    dimensionedScalar rho(mechanicalProperties.lookup("rho"));
    dimensionedScalar rhoE(mechanicalProperties.lookup("E"));
    dimensionedScalar nu(mechanicalProperties.lookup("nu"));

    dimensionedScalar E = rhoE/rho;
    dimensionedScalar mu = E/(2.0*(1.0 + nu));
    dimensionedScalar lambda = nu*E/((1.0 + nu)*(1.0 - 2.0*nu));
    dimensionedScalar threeK = E/(1.0 - 2.0*nu);

    Switch planeStress(mechanicalProperties.lookup("planeStress"));

    if (planeStress)
    {
        lambda = nu*E/((1.0 + nu)*(1.0 - nu));
        threeK = E/(1.0 - nu);
    }

    vectorField n = patch().nf();

    const fvPatchField<tensor>& gradU =
        patch().lookupPatchField<volTensorField, tensor>("grad(U)");

    gradient() =
    (
        (traction_ - pressure_*n)/rho.value()
      - (n & (mu.value()*gradU.T() - (mu + lambda).value()*gradU))
      - n*tr(gradU)*lambda.value()
    )/(2.0*mu + lambda).value();


    Switch thermalStress(thermalProperties.lookup("thermalStress"));

    if (thermalStress)
    {
        dimensionedScalar alpha(thermalProperties.lookup("alpha"));
        dimensionedScalar threeKalpha = threeK*alpha;

        const fvPatchField<scalar>& T =
            patch().lookupPatchField<volScalarField, scalar>("T");

        gradient() += n*threeKalpha.value()*T/(2.0*mu + lambda).value();
    }

    fixedGradientFvPatchVectorField::updateCoeffs();
}

bool IterativeSolvers::pcg(const IRCMatrix &A,
                           Vector &x,
                           const Vector &b,
                           const Preconditioner &M) {
    /*!
      Solves Ax=b using the preconditioned conjugate gradient method.
      */
    const idx N = x.getLength();
    real resid(100.0);
    Vector p(N), z(N), q(N);
    real alpha;
    real normr(0);
    real normb = norm(b);
    real rho(0), rho_1(0), beta(0);
    Vector r = b - A * x;
    if (normb == 0.0)
        normb = 1;
    resid = norm(r) / normb;
    if (resid <= IterativeSolvers::toler) {
        IterativeSolvers::toler = resid;
        IterativeSolvers::maxIter = 0;
        return true;
    }
    // MAIN LOOP
    idx i = 1;
    for (; i <= IterativeSolvers::maxIter; i++) {
        M.solveMxb(z, r);
        rho = dot(r, z);
        if (i == 1)
            p = z;
        else {
            beta = rho / rho_1;
            aypx(beta, p, z); // p = beta*p + z;
        }
        // CALCULATES q = A*p AND dp = dot(q,p)
        real dp = multiply_dot(A, p, q);
        alpha = rho / dp;
        normr = 0;
#ifdef USES_OPENMP
        #pragma omp parallel for reduction(+:normr)
#endif
        for (idx j = 0 ; j < N ; ++j) {
            x[j] += alpha * p[j]; // x + alpha(0) * p;
            r[j] -= alpha * q[j]; // r - alpha(0) * q;
            normr += r[j] * r[j];
        }
        normr = sqrt(normr);
        resid = normr / normb;
        if (resid <= IterativeSolvers::toler) {
            IterativeSolvers::toler = resid;
            IterativeSolvers::maxIter = i;
            return true;
        }
        rho_1 = rho;
    }
    IterativeSolvers::toler = resid;
    return false;
}

	itpp::Mat<std::complex<double> > PartialTrace(const int size_H_1){
  		int size_H_2=size() /size_H_1;
		itpp::Mat<std::complex<double> > rho(size_H_1,size_H_1);
		itpp::Vec<std::complex<double> > vec_row(size_H_2),vec_col(size_H_2);
		if (size_H_2*size_H_1 != size()){
			std::cerr  << "error PartialTrace" ; 
			exit(1);
		}
		for (int j_row=0 ; j_row< size_H_1 ; j_row++){
		vec_row=coefficients.mid(size_H_2*j_row,size_H_2);
			for (int j_col=j_row; j_col< size_H_1 ; j_col++){
				vec_col=coefficients.mid(size_H_2*j_col,size_H_2);
				rho(j_row,j_col)=dot(itpp::conj(vec_row),vec_col);
    				rho(j_col,j_row)= conj(rho(j_row,j_col));  
			}
		}
		return rho ;
	}

        void MakeRhoEnergy::saverho(std::string const & filename)
        {
            fp_.reset(std::fopen(filename.c_str(), "w"));

            for (auto i = 1; i <= max_; i++) {
                auto const r = static_cast<double>(i) * dx_;
                std::fprintf(fp_.get(), "%.15f, %.15f, %.15f\n", r, rho(alpha_ * r), exactrho(r));
            }
        }

int main (void) {
   int Nr = 10001;
   double rmax = 3.615*2.5;
   double dr = rmax/(double)Nr;
   int Nrho = 10001;
   double rhomax = rho (0.);
   double drho = rhomax/(double)Nrho;

   int atomic_number = 1;
   double mass = 63.55;
   double lattice_constant = 3.615;
   char lattice_type[] = "FCC";

   int i;

   char LAMMPSFilename[] = "Cu_Zhou.eam";
   FILE *LAMMPSFile = fopen (LAMMPSFilename, "w");
   if (!LAMMPSFile) exit (ERROR);

   // Header for setfl format
   fprintf (LAMMPSFile, \
         "#-> LAMMPS Potential File in DYNAMO 86 setfl Format <-#\n"\
         "# Zhou Cu Acta mater(2001)49:4005\n"\
         "# Implemented by G. Ziegenhain (2007) [email protected]\n"\
         "%d Cu\n"\
         "%d %20.20f %d %20.20f %20.20f\n"\
         "%d %20.20f %20.20f %s\n",
         atomic_number, 
         Nrho, drho, Nr, dr, rmax,
         atomic_number, mass, lattice_constant, lattice_type);

   // Embedding function
   for (i = 0; i < Nrho; i++) 
      fprintf (LAMMPSFile, "%20.20f\n", F ((double)i*drho));
   // Density function
   for (i = 0; i < Nr; i++) 
      fprintf (LAMMPSFile, "%20.20f\n", rho ((double)i*dr));
   // Pair potential
   for (i = 0; i < Nr; i++)   
      fprintf (LAMMPSFile, "%20.20f\n", V ((double)i*dr) * (double)i*dr);
   
   fclose (LAMMPSFile);
   return (EOK);
}

  vector<double> autocovariance(const vector<double>& x,unsigned max)
  {
    const int N = x.size();

    // specify sample size if unspecified
    if (max==0) 
      max = std::max(1+N/4,N-15);

    if (max >= N) 
      max = N-1;

    // compute mean of X
    double mean = 0;
    for(int i=0;i<N;i++)
      mean += x[i];
    mean /= N;

    // allocate covariances
    vector<double> rho(max);

    // Run iteration 0 separately - to use rho[0] in the limit calculation
    {
      int k=0;
      double total = 0;
      for(int i=0;i<N-k;i++)
	total += (x[i]-mean)*(x[i+k]-mean);

      rho[k] = total/(N-k);
    }

    // compute each autocorrelation rho[k]
    double limit = rho[0]*(0.01/N);
    for(int k=1;k<max;k++) 
    {
      double total = 0;
      for(int i=0;i<N-k;i++)
	total += (x[i]-mean)*(x[i+k]-mean);

      rho[k] = total/(N-k);

      if (rho[k] < limit) {
	if (rho[k] < 0)
	{
	  rho.resize(k-1);
	  break;
	}
	else if (k > 3)
	{
	  rho.resize(k);
	  break;
	}
      }
    }

    return rho;
  }

	int CG(const MMatrix &A, MVector &x, const MVector &b, const Preconditioner &M, int &max_iter, Real &tol)
	{
		Real resid;
		MVector p, z, q;
		MVector alpha(1), beta(1), rho(1), rho_1(1);

		MVector r = b - A*x;

		Real normb = norm(b);

		if (normb == 0.0) normb = 1;

		if ((resid = norm(r) / normb) <= tol){
			tol = resid;
			max_iter = 0;
			return 0;
		}

		for (int i = 1; i <= max_iter; i++) 
		{
			// Assign Z
			z = M.solve(r);

			rho.p_[0] = dot(r, z);

			// Assign P
			if (i == 1)
				p = z;
			else {
				beta.p_[0] = rho.p_[0] / rho_1.p_[0];
				p = z + beta.p_[0] * p;
			}

			// Assign Q
			q = A*p;

			alpha.p_[0] = rho.p_[0] / dot(p, q);

			// Change X and R
			x += alpha.p_[0] * p;
			r -= alpha.p_[0] * q;

			// Check tol
			if ((resid = norm(r) / normb) <= tol)
			{
				tol = resid;
				max_iter = i;
				return 0;     
			}

			rho_1.p_[0] = rho.p_[0];
		}

		tol = resid;
		return 1;
	}

void Foam::relativeVelocityModels::general::correct()
{
    Udm_ =
        (rhoc_/rho())
       *V0_
       *(
            exp(-a_*max(alphad_ - residualAlpha_, scalar(0)))
          - exp(-a1_*max(alphad_ - residualAlpha_, scalar(0)))
        );
}