C++ rho函数代码示例

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 ;
	}