C++ cbrt函数代码示例

int main(int argc, char *argv[]){
  // 物理定数を定義（Fe）
  material_t material = setup_material(0.4174, 1.3885, 2.845);
  // 座標長さ1における実際の長さ（格子定数）
  float length = 2.9098453;
  // float length = 2 * material.r0 / sqrt(3);
  // 座標ファイルの読み込み
  FILE *rfp = read_file();

  // ファイル一行目から総原子数を求め、atoms配列を定義
  long atoms_number;
  fscanf(rfp, "%ld", &atoms_number);
  atom_t atoms[atoms_number];

  // atomsに座標を代入（bcc_structure.h参照）
  input_structure(rfp, atoms, atoms_number);

  fclose(rfp);

  // 総原子数から構造の大きさを求める（原子間距離を原子半径の2倍と仮定）
  int structure_length = cbrt(atoms_number / 2.0);
  // float structure_size = structure_length * material.r0 * 2.0;

  // // 結合エネルギが有効数字内で変化しなくなる距離を計算し、一辺を何回ループするか決定（すべての原子が同じ物理定数を持つと仮定）
  // int loop_number = eval_loop_number(material, structure_size);

  // // 拡張したatomsを作成
  // long extended_atoms_number = atoms_number * pow(loop_number, 3);
  // atom_t extended_atoms[extended_atoms_number];

  // // 拡張したatomsに座標を定義
  // initialize_extended_structure(extended_atoms, atoms, atoms_number, structure_length, loop_number);

  // // エネルギーを計算
  // eval_energy(atoms, atoms_number, extended_atoms, extended_atoms_number, material);
  eval_energy_ex(atoms, atoms_number, material, structure_length, length);

  // 出力
  FILE *wfp = write_file();
  output_structure(wfp, atoms, atoms_number);
  fclose(wfp);

  return 0;
}

/**
 *  \brief
 *      Returns the three roots of the cubic equation with real coefficients.
 *
 *  \details
 *      Returns the three roots of the cubic equation;
 *
 *              \f[z^3 + b z^2 + c z + d = 0,\f]
 *
 *      where \f$b, c, d\f$ are real coefficients, and the coefficient on the
 *      \f$x^3\f$ term is assumed to be 1.
 *
 *      \param[in]      b   Real coefficient of the \f$z^2\f$ term.
 *      \param[in]      c   Real coefficient of the \f$z\f$ term.
 *      \param[in]      d   Real constant term.
 *      \param[out]     z1  First real root.
 *      \param[out]     z2  Second (possibly complex) root.
 *      \param[out]     z3  Third (possibly complex) root.
 *
 *      \return         The number of real roots
 *
 *      \author         Mike Henderson
 *      \date           2011
 *
 */
int  Lgm_CubicRoots( double b, double c, double d, double *z1, double complex *z2, double complex *z3 ){

    int     nReal;
    double  b2, Q, R, D2, D, RpD, RmD, S, T, Theta, SqrtNQ, f, g, h;

    b2 = b*b;

    Q  = (3.0*c - b2)/9.0;
    R  = (9.0*b*c - 27.0*d - 2.0*b2*b)/54.0;
    D2 = Q*Q*Q + R*R;

    if (D2 == 0){
        // All 3 roots are real and at least two are equal.
        if ( R > 0.0) {
            f = cbrt( R );
            *z1 = 2.0*f - b/3.0;
            *z2 = *z3 = -f - b/3.0;
        } else {
            f = -cbrt( -R );
            *z1 = 2.0*f - b/3.0;
            *z2 = *z3 = -f - b/3.0;
        }
        nReal = 3;
    } else if ( D2 > 0 ) {
        // only one real root exists. The other two are complex conjugates.
        D = sqrt(D2); RpD = R+D; RmD = R-D;
        S = (RpD > 0.0) ? cbrt( RpD ) : -cbrt( -RpD );
        T = (RmD > 0.0) ? cbrt( RmD ) : -cbrt( -RmD);
        f = S+T;
        g = sqrt(3.0)/2.0*(S-T);
        h = b/3.0;
        *z1 = f - h;
        *z2 = -0.5*f - h + g*I;
        *z3 = -0.5*f - h - g*I;
        nReal = 1;

    } else if ( D2 < 0 ) {
        // All three roots are real and unequal.
        // This uses the simplest of the alternative trig forms.
        SqrtNQ = sqrt(-Q);
        Theta = acos( R/(SqrtNQ*fabs(Q)) );
        *z1 = 2.0*SqrtNQ*cos(Theta/3.0) - b/3.0; // real (returned as real double z1)
        *z2 = 2.0*SqrtNQ*cos( (Theta + 2.0*M_PI)/3.0 ) - b/3.0; // real (but returned as double complex z2)
        *z3 = 2.0*SqrtNQ*cos( (Theta + 4.0*M_PI)/3.0 ) - b/3.0; // real (but returned as double complex z3)
        nReal = 3;

    }

    return( nReal );

}

std::pair<int, long long> recursion(long long n)
{
	if (n == 0 || cf[n].first != 0)
	{
		return cf[n];
	}
	long long a = (long long)cbrt(n);
	long long b = a-1;
	std::pair<int, long long> ap = recursion(n - a*a*a), bp = recursion(a*a*a - b*b*b - 1);
	if (ap.first >= bp.first)
	{
		cf[n] = std::make_pair(ap.first+1, ap.second + a*a*a);
	}
	else
	{
		cf[n] = std::make_pair(bp.first+1, bp.second + b*b*b);
	}
	return cf[n];
}

double dnormv (double *in, int size, int norm)
{
	double sum = 0;
	double res = 0;
	int counter = 0;

	switch (norm)
	{
		case 0:
			res = INFINITY;
			break;
		case 1:  /*Addition of all elements*/
			for (counter=0; counter < size; counter++)
			{
				sum += in[counter];
			}
			res = sum;
			break;
		case 2: /*Square root of addition of squares of all elements*/
			for (counter=0; counter < size; counter++)
			{
				sum += in[counter]*in[counter];
			}
			res = sqrt(sum);
			break;
		case 3: /*Cube root of addition of cubes of all elements*/
			for (counter=0; counter < size; counter++)
			{
				sum += pow(in[counter],3);
			}
			res = cbrt(sum);
			break;
		default : /*Nth root of addition of all elements raised to n*/
			for (counter=0; counter < size; counter++)
			{
				sum += pow(in[counter],norm);
			}
			res = pow(sum, 1./norm);
			break;
	}

	return res;
}

Foam::tmp<Foam::volScalarField> Foam::diameterModels::ADD::tauC() const
{
    //- Phase fraction
    const volScalarField& alpha = phase_;
    //- Turbulent kinetic energy
    const volScalarField& k = phase_.turbulence().k();

    return
    Cc_*
    (
        cbrt
        (
            constant::mathematical::pi/6.0
           *(alphaMax_ - alpha)/
            max(alpha,SMALL)
        )
       *d_/sqrt(2.0/3.0*k)
    );
}

void test_fp_cbrt( void )
{
#if __STDC_VERSION__ >= 199901L
    printf( "Testing C99 cube root function...\n" );

    VERIFY( CompDbl( cbrt( 0.0 ), 0.0 ) );
    VERIFY( CompDbl( cbrt( 27.0 ), 3.0 ) );
    VERIFY( CompDbl( cbrt( -27.0 ), -3.0 ) );
    
    VERIFY( isnan(cbrt( NAN )) );
    VERIFY( cbrt( INFINITY ) == INFINITY );
    VERIFY( cbrt( -INFINITY ) == -INFINITY );
#endif
}

void PelletAgglomerator::run(){
	//DemField& dem(field->cast<DemField>());
	// loop over all source particles, and loop over all contacts of each of them
	if(isnan(massIncPerRad)) throw std::runtime_error("PalletAgglomerator.massIncPerRad==NaN (must be specified)");
	if(dampHalfLife<0) dampHalfLife*=-scene->dt;
	Real sumDMass=0.;
	Real lambda=(dampHalfLife==0 || isnan(dampHalfLife))?0:(log(2)/dampHalfLife);
	for(const shared_ptr<Particle>& src: agglomSrcs){
		for(const auto& idCon: src->contacts){
			const shared_ptr<Contact>& c(idCon.second);
			if(!c->isReal()) continue;
			Particle* other(c->leakOther(src.get()));
			Real* radius=nullptr;
			// if(!other->shape->isA<Sphere>()) continue; // other particles is not a sphere
			if(other->shape->isA<Sphere>()) radius=&(other->shape->cast<Sphere>().radius);
			else if(other->shape->isA<Capsule>()) radius=&(other->shape->cast<Capsule>().radius);
			else continue;
			assert(dynamic_pointer_cast<L6Geom>(c->geom));
			// radius change
			// angVel is local already
			Real dMass=c->geom->cast<L6Geom>().angVel.tail<2>().norm()*scene->dt*massIncPerRad;
			sumDMass+=dMass;
			Real newVol=(4/3.)*M_PI*pow(*radius,3)+dMass/other->material->density;
			*radius=cbrt(3*newVol/(4*M_PI));
			other->shape->updateMassInertia(other->material->density);
			if(!other->matState) other->matState=make_shared<PelletMatState>();
			assert(dynamic_pointer_cast<PelletMatState>(other->matState));
			auto& pms=other->matState->cast<PelletMatState>();
			if(pms.stepAgglomUpdated!=scene->step){ pms.agglomRate=0.; pms.stepAgglomUpdated=scene->step; } // reset value
			pms.agglomRate+=dMass/scene->dt;
			// rotation damping
			if(lambda>0){
				other->shape->nodes[0]->getData<DemData>().angVel*=(1-lambda*scene->dt);
			}
			boost::mutex::scoped_lock l(pms.lock);
			pms.cumAgglomMass+=dMass;
			pms.cumAgglomAngle+=c->geom->cast<L6Geom>().angVel.tail<2>().norm()*scene->dt;
		}
	}
	currRate=(1-currRateSmooth)*currRate+currRateSmooth*(sumDMass/scene->dt);
	mass+=sumDMass;
};

void kdSmbhSoft(KD kd, int nSplitting)
{
  
  int i;
  double max_mass;
  
  // smbh particles are ones with the largest mass
  max_mass = kd->pInit[0].fMass;
  for (i=0;i<kd->nParticles;i++) {
    if (max_mass < kd->pInit[i].fMass) {
      max_mass = kd->pInit[i].fMass;
    }
  }
  
  fprintf(stderr, "MAX MASS = %f\n", max_mass);
  for (i=0;i<kd->nParticles;i++) 
    if(kd->pInit[i].fMass == max_mass) 
      kd->pInit[i].fSoft = kd->pInit[i].fSoft/cbrt((double)nSplitting);
    
}

/*! \brief
 * Determines a suitable grid size.
 *
 * \param[in,out] d    Grid information.
 * \param[in]     pbc  Information about the box.
 * \returns  FALSE if grid search is not suitable.
 */
static gmx_bool
grid_setup_cells(gmx_ana_nbsearch_t *d, t_pbc *pbc)
{
    real targetsize;
    int  dd;

#ifdef HAVE_CBRT
    targetsize = cbrt(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
                      * 10 / d->nref);
#else
    targetsize = pow(pbc->box[XX][XX] * pbc->box[YY][YY] * pbc->box[ZZ][ZZ]
                     * 10 / d->nref, 1./3.);
#endif

    d->ncells = 1;
    for (dd = 0; dd < DIM; ++dd)
    {
        d->ncelldim[dd] = (int)(pbc->box[dd][dd] / targetsize);
        d->ncells      *= d->ncelldim[dd];
        if (d->ncelldim[dd] < 3)
        {
            return FALSE;
        }
    }
    /* Reallocate if necessary */
    if (d->cells_nalloc < d->ncells)
    {
        int  i;

        srenew(d->ncatoms, d->ncells);
        srenew(d->catom, d->ncells);
        srenew(d->catom_nalloc, d->ncells);
        for (i = d->cells_nalloc; i < d->ncells; ++i)
        {
            d->catom[i]        = NULL;
            d->catom_nalloc[i] = 0;
        }
        d->cells_nalloc = d->ncells;
    }
    return TRUE;
}

Foam::tmp<Foam::volScalarField>
Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
(
    const phaseModel& phase1,
 const phaseModel& phase2
) const
{
    volScalarField aByaMax
    (
        cbrt
        (
            min
            (
                max(phase1, scalar(1e-3)),
                alphaMinFriction_
            )/phase1.alphaMax()
        )
    );

    return (1.0/(3*phase1.alphaMax()))/sqr(aByaMax - sqr(aByaMax));
}

// ratio = v0/v1.
inline double secant_rule(double x0,double x1,double ratio,
  int num_roots_in_cluster)
{
    switch(num_roots_in_cluster)
    {
     case 1:
      break;
     case 2:
      eprintf("disallowed");
      break;
     case 3:
      ratio = cbrt(ratio);
      break;
     default:
      assert_eq(1,num_roots_in_cluster%2);
      ratio = copysign(pow(fabs(ratio),1./num_roots_in_cluster),ratio);
    }
    double denominator = (1.-ratio);
    assert_gt(denominator,0.);
    return (x0-ratio*x1)/denominator;
}

static void Test_alpha_epsilon(void) {
  printf("\n** Test_alpha_epsilon: **\n");
  const REAL8 f = 0.01;
  const REAL8 q = 4;
  const REAL8 chil = 0.5625;
  const REAL8 chip = 0.18;

  NNLOanglecoeffs angcoeffs;
  ComputeNNLOanglecoeffs(&angcoeffs,q,chil,chip);

  const REAL8 omega = LAL_PI * f;
  const REAL8 logomega = log(omega);
  const REAL8 omega_cbrt = cbrt(omega);
  const REAL8 omega_cbrt2 = omega_cbrt*omega_cbrt;
  const REAL8 alpha = (angcoeffs.alphacoeff1/omega
                    + angcoeffs.alphacoeff2/omega_cbrt2
                    + angcoeffs.alphacoeff3/omega_cbrt
                    + angcoeffs.alphacoeff4*logomega
                    + angcoeffs.alphacoeff5*omega_cbrt);

  const REAL8 epsilon = (angcoeffs.epsiloncoeff1/omega
                      + angcoeffs.epsiloncoeff2/omega_cbrt2
                      + angcoeffs.epsiloncoeff3/omega_cbrt
                      + angcoeffs.epsiloncoeff4*logomega
                      + angcoeffs.epsiloncoeff5*omega_cbrt);

  const REAL8 alpha_expected = -11.8196;
  const REAL8 epsilon_expected = -11.936;

  print_difference("alpha", alpha, alpha_expected);
  print_difference("epsilon", epsilon, epsilon_expected);

  const REAL8 eps = 1e-5;

  assert(
       approximatelyEqual(alpha,    alpha_expected,   eps)
    && approximatelyEqual(epsilon,  epsilon_expected, eps)
    && "Test_alpha_epsilon()"
  );
}

void Foam::WallLocalSpringSliderDashpot<CloudType>::findMinMaxProperties
(
    scalar& rMin,
    scalar& rhoMax,
    scalar& UMagMax
) const
{
    rMin = VGREAT;
    rhoMax = -VGREAT;
    UMagMax = -VGREAT;

    forAllConstIter(typename CloudType, this->owner(), iter)
    {
        const typename CloudType::parcelType& p = iter();

        // Finding minimum diameter to avoid excessive arithmetic

        scalar dEff = p.d();

        if (useEquivalentSize_)
        {
            dEff *= cbrt(p.nParticle()*volumeFactor_);
        }

        rMin = min(dEff, rMin);

        rhoMax = max(p.rho(), rhoMax);

        UMagMax = max
        (
            mag(p.U()) + mag(p.omega())*dEff/2,
            UMagMax
        );
    }

    // Transform the minimum diameter into minimum radius
    //     rMin = dMin/2

    rMin /= 2.0;
}

void dump_histo_img(unsigned char* histo, unsigned int height, unsigned int width, const char *filename)
{
    RGB* pixel_map = (RGB*) malloc (height*width*sizeof(RGB));

    size_t y, x;
    for (y = 0; y < height; ++y)
    {
        for (x = 0; x < width; ++x)
        {
            unsigned char value = histo[y * width + x];

            if (value == 0){
                pixel_map[y*width+x].R = 0;
                pixel_map[y*width+x].G = 0;
                pixel_map[y*width+x].B = 0;
            } else {
                pixel_map[y*width+x] = HSVtoRGB(0.0,1.0,cbrt(1+ 63.0*((float)value)/((float)UINT8_MAX))/4);
            }
        }
    }
    create_bmp(pixel_map, height, width, filename);
    free(pixel_map);
}

static TACommandVerdict cbrt_cmd(TAThread thread,TAInputStream stream)
{
    double x, res;

    // Prepare

    x = readDouble(&stream);
    
    START_TARGET_OPERATION(thread);
    
    // Execute

    res = cbrt(x);
    
    END_TARGET_OPERATION(thread);
    
    // Response
    
    writeDouble(thread, res);

    sendResponse(thread);
    return taDefaultVerdict;
}