// Check that the small body has the right degrees of freedom in the frame of
// the big body.
TEST_F(BodyCentredNonRotatingDynamicFrameTest, SmallBodyInBigFrame) {
  int const kSteps = 100;
  Bivector<double, ICRFJ2000Equator> const axis({0, 0, 1});

  RelativeDegreesOfFreedom<ICRFJ2000Equator> const initial_big_to_small =
      small_initial_state_ - big_initial_state_;
  ContinuousTrajectory<ICRFJ2000Equator>::Hint hint;
  for (Instant t = t0_; t < t0_ + 1 * period_; t += period_ / kSteps) {
    DegreesOfFreedom<ICRFJ2000Equator> const small_in_inertial_frame_at_t =
        solar_system_.trajectory(*ephemeris_, kSmall).
            EvaluateDegreesOfFreedom(t, &hint);

    auto const rotation_in_big_frame_at_t =
        Rotation<ICRFJ2000Equator, Big>(-2 * π * (t - t0_) * Radian / period_,
                                        axis);
    DegreesOfFreedom<Big> const small_in_big_frame_at_t(
        rotation_in_big_frame_at_t(initial_big_to_small.displacement()) +
            Big::origin,
        rotation_in_big_frame_at_t(initial_big_to_small.velocity()));

    auto const to_big_frame_at_t = big_frame_->ToThisFrameAtTime(t);
    EXPECT_THAT(AbsoluteError(
                    to_big_frame_at_t(small_in_inertial_frame_at_t).position(),
                    small_in_big_frame_at_t.position()),
                Lt(0.3 * Milli(Metre)));
    EXPECT_THAT(AbsoluteError(
                    to_big_frame_at_t(small_in_inertial_frame_at_t).velocity(),
                    small_in_big_frame_at_t.velocity()),
                Lt(4 * Milli(Metre) / Second));
  }
}

    void
    Matrix::lu(Matrix& L, Matrix& U) const
    {
      if (m_nrows != m_ncols)
        throw Error(" matrix is not square ");

      Matrix A = *this;
      Matrix Lf(m_nrows), dI(m_nrows), P(m_nrows);
      dI = dI + dI;

      for (size_t i = 0; i < m_nrows - 1; i++)
      {
        /*
        // check if pivot == 0
        // if so, try to find a valid pivot
        // if no valid pivot found, matrix isn't invertible (quit)
        // update permutation matrix
        //
        */
        Matrix Lt(m_nrows);
        for (size_t j = i + 1; j < m_nrows; j++)
          Lt(j, i) = -A(j, i) / A(i, i);

        A = Lt * A;
        Lf = Lf * (dI - Lt);
      }

      U = A;
      L = Lf;
    }

void TestSymplecticity(Integrator const& integrator,
                       Energy const& expected_energy_error) {
  Length const q_initial = 1 * Metre;
  Speed const v_initial = 0 * Metre / Second;
  Instant const t_initial;
  Instant const t_final = t_initial + 500 * Second;
  Time const step = 0.2 * Second;

  Mass const m = 1 * Kilogram;
  Stiffness const k = SIUnit<Stiffness>();
  Energy const initial_energy =
      0.5 * m * Pow<2>(v_initial) + 0.5 * k * Pow<2>(q_initial);

  std::vector<ODE::SystemState> solution;
  ODE harmonic_oscillator;
  harmonic_oscillator.compute_acceleration =
      std::bind(ComputeHarmonicOscillatorAcceleration,
                _1, _2, _3, /*evaluations=*/nullptr);
  IntegrationProblem<ODE> problem;
  problem.equation = harmonic_oscillator;
  ODE::SystemState const initial_state = {{q_initial}, {v_initial}, t_initial};
  problem.initial_state = &initial_state;
  auto append_state = [&solution](ODE::SystemState const& state) {
    solution.push_back(state);
  };

  auto const instance =
      integrator.NewInstance(problem, std::move(append_state), step);
  integrator.Solve(t_final, *instance);

  std::size_t const length = solution.size();
  std::vector<Energy> energy_error(length);
  std::vector<Time> time(length);
  Energy max_energy_error;
  for (std::size_t i = 0; i < length; ++i) {
    Length const q_i   = solution[i].positions[0].value;
    Speed const v_i = solution[i].velocities[0].value;
    time[i] = solution[i].time.value - t_initial;
    energy_error[i] =
        AbsoluteError(initial_energy,
                      0.5 * m * Pow<2>(v_i) + 0.5 * k * Pow<2>(q_i));
    max_energy_error = std::max(energy_error[i], max_energy_error);
  }
  double const correlation =
      PearsonProductMomentCorrelationCoefficient(time, energy_error);
  LOG(INFO) << "Correlation between time and energy error : " << correlation;
  EXPECT_THAT(correlation, Lt(1e-2));
  Power const slope = Slope(time, energy_error);
  LOG(INFO) << "Slope                                     : " << slope;
  EXPECT_THAT(Abs(slope), Lt(2e-6 * SIUnit<Power>()));
  LOG(INFO) << "Maximum energy error                      : " <<
      max_energy_error;
  EXPECT_EQ(expected_energy_error, max_energy_error);
}

TEST_F(RandomTests, NextNumber_TwoArgument) {
    int a = rnd.nextInt(10, 20);
    EXPECT_THAT(a, Ge(10));
    EXPECT_THAT(a, Lt(20));

    long long b = rnd.nextLongLong(1000000000000ll, 2000000000000ll);
    EXPECT_THAT(b, Ge(1000000000000ll));
    EXPECT_THAT(b, Lt(2000000000000ll));

    double c = rnd.nextDouble(100.0, 200.0);
    EXPECT_THAT(c, Ge(100.0));
    EXPECT_THAT(c, Le(200.0));
}

Spectrum ExPhotonIntegrator::LPhoton(
		KdTree<Photon, PhotonProcess> *map,
		int nPaths, int nLookup, BSDF *bsdf,
		const Intersection &isect, const Vector &wo,
		float maxDistSquared) {
	Spectrum L(0.);
	if (!map) return L;
	BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
		BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
	if (bsdf->NumComponents(nonSpecular) == 0)
		return L;
	static StatsCounter lookups("Photon Map", "Total lookups"); // NOBOOK
	// Initialize _PhotonProcess_ object, _proc_, for photon map lookups
	PhotonProcess proc(nLookup, isect.dg.p);
	proc.photons =
		(ClosePhoton *)alloca(nLookup * sizeof(ClosePhoton));
	// Do photon map lookup
	++lookups;  // NOBOOK
	map->Lookup(isect.dg.p, proc, maxDistSquared);
	// Accumulate light from nearby photons
	static StatsRatio foundRate("Photon Map", "Photons found per lookup"); // NOBOOK
	foundRate.Add(proc.foundPhotons, 1); // NOBOOK
	// Estimate reflected light from photons
	ClosePhoton *photons = proc.photons;
	int nFound = proc.foundPhotons;
	Normal Nf = Dot(wo, bsdf->dgShading.nn) < 0 ? -bsdf->dgShading.nn :
		bsdf->dgShading.nn;

	if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |
			BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {
		// Compute exitant radiance from photons for glossy surface
		for (int i = 0; i < nFound; ++i) {
			const Photon *p = photons[i].photon;
			BxDFType flag = Dot(Nf, p->wi) > 0.f ?
				BSDF_ALL_REFLECTION : BSDF_ALL_TRANSMISSION;
			float k = kernel(p, isect.dg.p, maxDistSquared);
			L += (k / nPaths) * bsdf->f(wo, p->wi, flag) * p->alpha;
		}
	}
	else {
		// Compute exitant radiance from photons for diffuse surface
		Spectrum Lr(0.), Lt(0.);
		for (int i = 0; i < nFound; ++i) {
			if (Dot(Nf, photons[i].photon->wi) > 0.f) {
				float k = kernel(photons[i].photon, isect.dg.p,
					maxDistSquared);
				Lr += (k / nPaths) * photons[i].photon->alpha;
			}
			else {
				float k = kernel(photons[i].photon, isect.dg.p,
					maxDistSquared);
				Lt += (k / nPaths) * photons[i].photon->alpha;
			}
		}
		L += Lr * bsdf->rho(wo, BSDF_ALL_REFLECTION) * INV_PI +
			Lt * bsdf->rho(wo, BSDF_ALL_TRANSMISSION) * INV_PI;
	}
	return L;
}

TEST_P(SimpleHarmonicMotionTest, Symplecticity) {
  parameters_.initial.positions.emplace_back(SIUnit<Length>());
  parameters_.initial.momenta.emplace_back(Speed());
  parameters_.initial.time = Time();
  Stiffness const k = SIUnit<Stiffness>();
  Mass const m      = SIUnit<Mass>();
  Length const q0   = parameters_.initial.positions[0].value;
  Speed const v0 = parameters_.initial.momenta[0].value;
  Energy const initial_energy = 0.5 * m * Pow<2>(v0) + 0.5 * k * Pow<2>(q0);
  parameters_.tmax = 500.0 * SIUnit<Time>();
  parameters_.Δt = 0.2 * Second;
  parameters_.sampling_period = 1;
  integrator_->SolveTrivialKineticEnergyIncrement<Length>(
      &ComputeHarmonicOscillatorAcceleration,
      parameters_,
      &solution_);
  std::size_t const length = solution_.size();
  std::vector<Energy> energy_error(length);
  std::vector<Time> time_steps(length);
  Energy max_energy_error = 0 * SIUnit<Energy>();
  for (std::size_t i = 0; i < length; ++i) {
    Length const q_i   = solution_[i].positions[0].value;
    Speed const v_i = solution_[i].momenta[0].value;
    time_steps[i] = solution_[i].time.value;
    energy_error[i] = Abs(0.5 * m * Pow<2>(v_i) + 0.5 * k * Pow<2>(q_i) -
                          initial_energy);
    max_energy_error = std::max(energy_error[i], max_energy_error);
  }
#if 0
  LOG(INFO) << "Energy error as a function of time:\n" <<
      BidimensionalDatasetMathematicaInput(time_steps, energy_error);
#endif
  double const correlation =
      PearsonProductMomentCorrelationCoefficient(time_steps, energy_error);
  LOG(INFO) << GetParam();
  LOG(INFO) << "Correlation between time and energy error : " << correlation;
  EXPECT_THAT(correlation, Lt(2E-3));
  Power const slope = Slope(time_steps, energy_error);
  LOG(INFO) << "Slope                                     : " << slope;
  EXPECT_THAT(Abs(slope), Lt(2E-6 * SIUnit<Power>()));
  LOG(INFO) << "Maximum energy error                      : " <<
      max_energy_error;
  EXPECT_EQ(GetParam().expected_energy_error, max_energy_error);
}

Spectrum LPhoton(KdTree<Photon> *map, int nPaths, int nLookup,
      MemoryArena &arena, BSDF *bsdf, RNG &rng, const Intersection &isect,
      const Vector &wo, float maxDistSquared) {
    Spectrum L(0.);
    BxDFType nonSpecular = BxDFType(BSDF_REFLECTION |
        BSDF_TRANSMISSION | BSDF_DIFFUSE | BSDF_GLOSSY);
    if (map && bsdf->NumComponents(nonSpecular) > 0) {
        PBRT_PHOTON_MAP_STARTED_LOOKUP(const_cast<DifferentialGeometry *>(&isect.dg));
        // Do photon map lookup at intersection point
        PhotonProcess proc(nLookup, arena.Alloc<ClosePhoton>(nLookup));
        map->Lookup(isect.dg.p, proc, maxDistSquared);

        // Estimate reflected radiance due to incident photons
        ClosePhoton *photons = proc.photons;
        int nFound = proc.nFound;
        Normal Nf = Faceforward(bsdf->dgShading.nn, wo);
        if (bsdf->NumComponents(BxDFType(BSDF_REFLECTION |
                BSDF_TRANSMISSION | BSDF_GLOSSY)) > 0) {
            // Compute exitant radiance from photons for glossy surface
            for (int i = 0; i < nFound; ++i) {
                const Photon *p = photons[i].photon;
                float k = kernel(p, isect.dg.p, maxDistSquared);
                L += (k / (nPaths * maxDistSquared)) * bsdf->f(wo, p->wi) *
                     p->alpha;
            }
        }
        else {
            // Compute exitant radiance from photons for diffuse surface
            Spectrum Lr(0.), Lt(0.);
            for (int i = 0; i < nFound; ++i) {
                if (Dot(Nf, photons[i].photon->wi) > 0.f) {
                    float k = kernel(photons[i].photon, isect.dg.p,
                        maxDistSquared);
                    Lr += (k / (nPaths * maxDistSquared)) * photons[i].photon->alpha;
                }
                else {
                    float k = kernel(photons[i].photon, isect.dg.p,
                        maxDistSquared);
                    Lt += (k / (nPaths * maxDistSquared)) * photons[i].photon->alpha;
                }
            }
            const int sqrtRhoSamples = 4;
            float rhoRSamples[2*sqrtRhoSamples*sqrtRhoSamples];
            StratifiedSample2D(rhoRSamples, sqrtRhoSamples, sqrtRhoSamples, rng);
            float rhoTSamples[2*sqrtRhoSamples*sqrtRhoSamples];
            StratifiedSample2D(rhoTSamples, sqrtRhoSamples, sqrtRhoSamples, rng);
            L += Lr * bsdf->rho(wo, sqrtRhoSamples*sqrtRhoSamples, rhoRSamples,
                                BSDF_ALL_REFLECTION) * INV_PI +
                 Lt * bsdf->rho(wo, sqrtRhoSamples*sqrtRhoSamples, rhoTSamples,
                                BSDF_ALL_TRANSMISSION) * INV_PI;
        }
        PBRT_PHOTON_MAP_FINISHED_LOOKUP(const_cast<DifferentialGeometry *>(&isect.dg),
            proc.nFound, proc.nLookup, &L);
    }
    return L;
}

tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::FDES
(
    const volScalarField::Internal& F1,
    const volScalarField::Internal& F2
) const
{
    switch (FSST_)
    {
        case 0:
            return max(Lt()/(CDES_*this->delta()()), scalar(1));
        case 1:
            return max(Lt()*(1 - F1)/(CDES_*this->delta()()), scalar(1));
        case 2:
            return max(Lt()*(1 - F2)/(CDES_*this->delta()()), scalar(1));
        default:
            FatalErrorInFunction
                << "Incorrect FSST = " << FSST_ << ", should be 0, 1 or 2"
                << exit(FatalError);
            return F1;
    }
}

bool cellLess(Cell cell, const ArrayData* val) {
    return cellRelOp(Lt(), cell, val);
}

 static bool Eq(Value lhs, Value rhs) { return !Lt(lhs, rhs) && !Gt(lhs, rhs); }

bool cellLess(const Cell* cell, double val) {
  return cellRelOp(Lt(), cell, val);
}

TEST_F(ManœuvreTest, Apollo8SIVB) {
    // Data from NASA's Saturn V Launch Vehicle, Flight Evaluation Report AS-503,
    // Apollo 8 Mission (1969),
    // http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19690015314.pdf.
    // We use the reconstructed or actual values.

    // Table 2-2. Significant Event Times Summary.
    Instant const range_zero;
    Instant const s_ivb_1st_90_percent_thrust = range_zero +    530.53 * Second;
    Instant const s_ivb_1st_eco               = range_zero +    684.98 * Second;
    // Initiate S-IVB Restart Sequence and Start of Time Base 6 (T6).
    Instant const t6                          = range_zero +   9659.54 * Second;
    Instant const s_ivb_2nd_90_percent_thrust = range_zero + 10'240.02 * Second;
            Instant const s_ivb_2nd_eco               = range_zero + 10'555.51 * Second;

    // From Table 7-2. S-IVB Steady State Performance - First Burn.
    Force thrust_1st                  = 901'557 * Newton;
                                        Speed specific_impulse_1st        = 4'204.1 * Newton * Second / Kilogram;
    Variation<Mass> lox_flowrate_1st  = 178.16 * Kilogram / Second;
    Variation<Mass> fuel_flowrate_1st = 36.30 * Kilogram / Second;

    // From Table 7-7. S-IVB Steady State Performance - Second Burn.
    Force thrust_2nd                  = 897'548 * Newton;
                                        Speed specific_impulse_2nd        = 4199.2 * Newton * Second / Kilogram;
                                        Variation<Mass> lox_flowrate_2nd  = 177.70 * Kilogram / Second;
                                        Variation<Mass> fuel_flowrate_2nd = 36.01 * Kilogram / Second;

                                        // Table 21-5. Total Vehicle Mass, S-IVB First Burn Phase, Kilograms.
                                        Mass total_vehicle_at_s_ivb_1st_90_percent_thrust = 161143 * Kilogram;
                                        Mass total_vehicle_at_s_ivb_1st_eco               = 128095 * Kilogram;

                                        // Table 21-7. Total Vehicle Mass, S-IVB Second Burn Phase, Kilograms.
                                        Mass total_vehicle_at_s_ivb_2nd_90_percent_thrust = 126780 * Kilogram;
                                        Mass total_vehicle_at_s_ivb_2nd_eco               =  59285 * Kilogram;

                                        // An arbitrary direction, we're not testing this.
                                        Vector<double, World> e_y({0, 1, 0});

    Manœuvre<World> first_burn(thrust_1st,
                                total_vehicle_at_s_ivb_1st_90_percent_thrust,
                                specific_impulse_1st, e_y);
    EXPECT_THAT(RelativeError(lox_flowrate_1st + fuel_flowrate_1st,
                              first_burn.mass_flow()),
                Lt(1E-4));

    first_burn.set_duration(s_ivb_1st_eco - s_ivb_1st_90_percent_thrust);
    EXPECT_THAT(
        RelativeError(total_vehicle_at_s_ivb_1st_eco, first_burn.final_mass()),
        Lt(1E-3));

    first_burn.set_initial_time(s_ivb_1st_90_percent_thrust);
    EXPECT_EQ(s_ivb_1st_eco, first_burn.final_time());

    // Accelerations from Figure 4-4. Ascent Trajectory Acceleration Comparison.
    // Final acceleration from Table 4-2. Comparison of Significant Trajectory
    // Events.
    EXPECT_THAT(
        first_burn.acceleration()(first_burn.initial_time()).Norm(),
        AllOf(Gt(5 * Metre / Pow<2>(Second)), Lt(6.25 * Metre / Pow<2>(Second))));
    EXPECT_THAT(first_burn.acceleration()(range_zero + 600 * Second).Norm(),
                AllOf(Gt(6.15 * Metre / Pow<2>(Second)),
                      Lt(6.35 * Metre / Pow<2>(Second))));
    EXPECT_THAT(first_burn.acceleration()(first_burn.final_time()).Norm(),
                AllOf(Gt(7.03 * Metre / Pow<2>(Second)),
                      Lt(7.05 * Metre / Pow<2>(Second))));

    Manœuvre<World> second_burn(thrust_2nd,
                                 total_vehicle_at_s_ivb_2nd_90_percent_thrust,
                                 specific_impulse_2nd, e_y);
    EXPECT_THAT(RelativeError(lox_flowrate_2nd + fuel_flowrate_2nd,
                              second_burn.mass_flow()),
                Lt(2E-4));

    second_burn.set_duration(s_ivb_2nd_eco - s_ivb_2nd_90_percent_thrust);
    EXPECT_THAT(
        RelativeError(total_vehicle_at_s_ivb_2nd_eco, second_burn.final_mass()),
        Lt(2E-3));

    second_burn.set_initial_time(s_ivb_2nd_90_percent_thrust);
    EXPECT_EQ(s_ivb_2nd_eco, second_burn.final_time());

    // Accelerations from Figure 4-9. Injection Phase Acceleration Comparison.
    // Final acceleration from Table 4-2. Comparison of Significant Trajectory
    // Events.
    EXPECT_THAT(second_burn.acceleration()(second_burn.initial_time()).Norm(),
                AllOf(Gt(7 * Metre / Pow<2>(Second)),
                      Lt(7.5 * Metre / Pow<2>(Second))));
    EXPECT_THAT(second_burn.acceleration()(t6 + 650 * Second).Norm(),
                AllOf(Gt(8 * Metre / Pow<2>(Second)),
                      Lt(8.02 * Metre / Pow<2>(Second))));
    EXPECT_THAT(second_burn.acceleration()(t6 + 700 * Second).Norm(),
                AllOf(Gt(8.8 * Metre / Pow<2>(Second)),
                      Lt(9 * Metre / Pow<2>(Second))));
    EXPECT_THAT(second_burn.acceleration()(t6 + 750 * Second).Norm(),
                AllOf(Gt(9.9 * Metre / Pow<2>(Second)),
                      Lt(10 * Metre / Pow<2>(Second))));
    EXPECT_THAT(second_burn.acceleration()(t6 + 850 * Second).Norm(),
                AllOf(Gt(12.97 * Metre / Pow<2>(Second)),
                      Lt(13 * Metre / Pow<2>(Second))));
    EXPECT_THAT(second_burn.acceleration()(second_burn.final_time()).Norm(),
//.........这里部分代码省略.........

void TestConvergence(Integrator const& integrator,
                     Time const& beginning_of_convergence) {
  Length const q_initial = 1 * Metre;
  Speed const v_initial = 0 * Metre / Second;
  Speed const v_amplitude = 1 * Metre / Second;
  AngularFrequency const ω = 1 * Radian / Second;
  Instant const t_initial;
  Instant const t_final = t_initial + 100 * Second;

  Time step = beginning_of_convergence;
  int const step_sizes = 50;
  double const step_reduction = 1.1;
  std::vector<double> log_step_sizes;
  log_step_sizes.reserve(step_sizes);
  std::vector<double> log_q_errors;
  log_step_sizes.reserve(step_sizes);
  std::vector<double> log_p_errors;
  log_step_sizes.reserve(step_sizes);

  std::vector<ODE::SystemState> solution;
  ODE harmonic_oscillator;
  harmonic_oscillator.compute_acceleration =
      std::bind(ComputeHarmonicOscillatorAcceleration,
                _1, _2, _3, /*evaluations=*/nullptr);
  IntegrationProblem<ODE> problem;
  problem.equation = harmonic_oscillator;
  ODE::SystemState const initial_state = {{q_initial}, {v_initial}, t_initial};
  problem.initial_state = &initial_state;
  ODE::SystemState final_state;
  auto const append_state = [&final_state](ODE::SystemState const& state) {
    final_state = state;
  };

  for (int i = 0; i < step_sizes; ++i, step /= step_reduction) {
    auto const instance = integrator.NewInstance(problem, append_state, step);
    integrator.Solve(t_final, *instance);
    Time const t = final_state.time.value - t_initial;
    Length const& q = final_state.positions[0].value;
    Speed const& v = final_state.velocities[0].value;
    double const log_q_error = std::log10(
        AbsoluteError(q / q_initial, Cos(ω * t)));
    double const log_p_error = std::log10(
        AbsoluteError(v / v_amplitude, -Sin(ω * t)));
    if (log_q_error <= -13 || log_p_error <= -13) {
      // If we keep going the effects of finite precision will drown out
      // convergence.
      break;
    }
    log_step_sizes.push_back(std::log10(step / Second));
    log_q_errors.push_back(log_q_error);
    log_p_errors.push_back(log_p_error);
  }
  double const q_convergence_order = Slope(log_step_sizes, log_q_errors);
  double const q_correlation =
      PearsonProductMomentCorrelationCoefficient(log_step_sizes, log_q_errors);
  LOG(INFO) << "Convergence order in q : " << q_convergence_order;
  LOG(INFO) << "Correlation            : " << q_correlation;

#if !defined(_DEBUG)
  EXPECT_THAT(RelativeError(integrator.order, q_convergence_order),
              Lt(0.05));
  EXPECT_THAT(q_correlation, AllOf(Gt(0.99), Lt(1.01)));
#endif
  double const v_convergence_order = Slope(log_step_sizes, log_p_errors);
  double const v_correlation =
      PearsonProductMomentCorrelationCoefficient(log_step_sizes, log_p_errors);
  LOG(INFO) << "Convergence order in p : " << v_convergence_order;
  LOG(INFO) << "Correlation            : " << v_correlation;
#if !defined(_DEBUG)
  // SPRKs with odd convergence order have a higher convergence order in p.
  EXPECT_THAT(
      RelativeError(integrator.order + (integrator.order % 2),
                    v_convergence_order),
      Lt(0.03));
  EXPECT_THAT(v_correlation, AllOf(Gt(0.99), Lt(1.01)));
#endif
}

void Parser::executeAction(int production) try {
    if (d_token__ != _UNDETERMINED_)
        pushToken__(d_token__); // save an already available token

    // $insert defaultactionreturn
    // save default non-nested block $$
    if (int size = s_productionInfo[production].d_size)
        d_val__ = d_vsp__[1 - size];

    switch (production) {
        // $insert actioncases

        case 1:
#line 43 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = d_vsp__[0].data<Tag__::basic>();
            res = d_val__.get<Tag__::basic>();
        } break;

        case 2:
#line 51 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = add(d_vsp__[-2].data<Tag__::basic>(),
                                              d_vsp__[0].data<Tag__::basic>());
        } break;

        case 3:
#line 54 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = sub(d_vsp__[-2].data<Tag__::basic>(),
                                              d_vsp__[0].data<Tag__::basic>());
        } break;

        case 4:
#line 57 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = mul(d_vsp__[-2].data<Tag__::basic>(),
                                              d_vsp__[0].data<Tag__::basic>());
        } break;

        case 5:
#line 60 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = div(d_vsp__[-2].data<Tag__::basic>(),
                                              d_vsp__[0].data<Tag__::basic>());
        } break;

        case 6:
#line 63 "parser.yy"
        {
            auto tup = parse_implicit_mul(d_vsp__[-2].data<Tag__::string>());
            if (neq(*std::get<1>(tup), *one)) {
                d_val__.get<Tag__::basic>() = mul(
                    std::get<0>(tup),
                    pow(std::get<1>(tup), d_vsp__[0].data<Tag__::basic>()));
            } else {
                d_val__.get<Tag__::basic>()
                    = pow(std::get<0>(tup), d_vsp__[0].data<Tag__::basic>());
            }
        } break;

        case 7:
#line 73 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = pow(d_vsp__[-2].data<Tag__::basic>(),
                                              d_vsp__[0].data<Tag__::basic>());
        } break;

        case 8:
#line 76 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(
                Lt(d_vsp__[-2].data<Tag__::basic>(),
                   d_vsp__[0].data<Tag__::basic>()));
        } break;

        case 9:
#line 79 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(
                Gt(d_vsp__[-2].data<Tag__::basic>(),
                   d_vsp__[0].data<Tag__::basic>()));
        } break;

        case 10:
#line 82 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(
                Le(d_vsp__[-2].data<Tag__::basic>(),
                   d_vsp__[0].data<Tag__::basic>()));
        } break;

        case 11:
#line 85 "parser.yy"
        {
            d_val__.get<Tag__::basic>() = rcp_static_cast<const Basic>(
                Ge(d_vsp__[-2].data<Tag__::basic>(),
                   d_vsp__[0].data<Tag__::basic>()));
        } break;

//.........这里部分代码省略.........

HOT_FUNC
bool tvLess(const TypedValue* tv1, const TypedValue* tv2) {
  return tvRelOp(Lt(), tv1, tv2);
}

TEST_F(EmbeddedExplicitRungeKuttaNyströmIntegratorTest,
       MaxSteps) {
  AdaptiveStepSizeIntegrator<ODE> const& integrator =
      DormandElMikkawyPrince1986RKN434FM<Length>();
  Length const x_initial = 1 * Metre;
  Speed const v_initial = 0 * Metre / Second;
  Speed const v_amplitude = 1 * Metre / Second;
  Time const period = 2 * π * Second;
  AngularFrequency const ω = 1 * Radian / Second;
  Instant const t_initial;
  Instant const t_final = t_initial + 10 * period;
  Length const length_tolerance = 1 * Milli(Metre);
  Speed const speed_tolerance = 1 * Milli(Metre) / Second;
  // The number of steps if no step limit is set.
  std::int64_t const steps_forward = 132;

  int evaluations = 0;
  auto const step_size_callback = [](bool tolerable) {};

  std::vector<ODE::SystemState> solution;
  ODE harmonic_oscillator;
  harmonic_oscillator.compute_acceleration =
      std::bind(ComputeHarmonicOscillatorAcceleration,
                _1, _2, _3, &evaluations);
  IntegrationProblem<ODE> problem;
  problem.equation = harmonic_oscillator;
  ODE::SystemState const initial_state = {{x_initial}, {v_initial}, t_initial};
  problem.initial_state = &initial_state;
  problem.t_final = t_final;
  problem.append_state = [&solution](ODE::SystemState const& state) {
    solution.push_back(state);
  };
  AdaptiveStepSize<ODE> adaptive_step_size;
  adaptive_step_size.first_time_step = t_final - t_initial;
  adaptive_step_size.safety_factor = 0.9;
  adaptive_step_size.tolerance_to_error_ratio =
      std::bind(HarmonicOscillatorToleranceRatio,
                _1, _2, length_tolerance, speed_tolerance, step_size_callback);
  adaptive_step_size.max_steps = 100;

  auto const outcome = integrator.Solve(problem, adaptive_step_size);
  EXPECT_EQ(termination_condition::ReachedMaximalStepCount, outcome.error());
  EXPECT_THAT(AbsoluteError(
                  x_initial * Cos(ω * (solution.back().time.value - t_initial)),
                      solution.back().positions[0].value),
              AllOf(Ge(8e-4 * Metre), Le(9e-4 * Metre)));
  EXPECT_THAT(AbsoluteError(
                  -v_amplitude *
                      Sin(ω * (solution.back().time.value - t_initial)),
                  solution.back().velocities[0].value),
              AllOf(Ge(1e-3 * Metre / Second), Le(2e-3 * Metre / Second)));
  EXPECT_THAT(solution.back().time.value, Lt(t_final));
  EXPECT_EQ(100, solution.size());

  // Check that a |max_steps| greater than or equal to the unconstrained number
  // of steps has no effect.
  for (std::int64_t const max_steps :
       {steps_forward, steps_forward + 1234}) {
    solution.clear();
    adaptive_step_size.max_steps = steps_forward;
    auto const outcome = integrator.Solve(problem, adaptive_step_size);
    EXPECT_EQ(termination_condition::Done, outcome.error());
    EXPECT_THAT(AbsoluteError(x_initial, solution.back().positions[0].value),
                AllOf(Ge(3e-4 * Metre), Le(4e-4 * Metre)));
    EXPECT_THAT(AbsoluteError(v_initial, solution.back().velocities[0].value),
                AllOf(Ge(2e-3 * Metre / Second), Le(3e-3 * Metre / Second)));
    EXPECT_EQ(t_final, solution.back().time.value);
    EXPECT_EQ(steps_forward, solution.size());
  }
}

bool tvLess(TypedValue tv1, TypedValue tv2) {
    return tvRelOp(Lt(), tv1, tv2);
}