C++ ex::degree方法代码示例

// GetCoeffs
vector<GiNaC::ex> GetCoeffs(GiNaC::ex Ex, GiNaC::symbol Sym) {
  assert(Ex.is_polynomial(Sym));
  vector<GiNaC::ex> Coeffs;
  unsigned Degree = Ex.degree(Sym);
  for (unsigned Idx = 0; Idx <= Degree; ++Idx)
    Coeffs.push_back(Ex.coeff(Sym, Idx));
  return Coeffs;
}

// Solve
vector<double> Solve(GiNaC::ex Ex, GiNaC::symbol Sym) {
  vector<double> Roots;

  unsigned Degree = Ex.degree(Sym);
  auto Coeffs = GetCoeffs(Ex, Sym);
  // Bhaskara.
  if (Degree == 2) {
    GiNaC::ex A = Coeffs[2];
    GiNaC::ex B = Coeffs[1];
    GiNaC::ex C = Coeffs[0];

    GiNaC::ex Delta = B*B - 4 * A * C;
    // Guaranteed real roots.
    if (GiNaC::is_a<GiNaC::numeric>(Delta) &&
        !GiNaC::ex_to<GiNaC::numeric>(Delta).is_negative()) {
      GiNaC::ex Delta = GiNaC::sqrt(B*B - 4 * A * C).evalf();
      GiNaC::ex One = ((-B) + Delta)/(2*A);
      GiNaC::ex Two = ((-B) - Delta)/(2*A);
      if (GiNaC::is_a<GiNaC::numeric>(One))
        Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(One.evalf()).to_double());
      if (GiNaC::is_a<GiNaC::numeric>(Two))
        Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(Two.evalf()).to_double());
    }
  }
  // Cardano.
  else if (Degree == 3) {
    GiNaC::ex A = Coeffs[3];
    GiNaC::ex B = Coeffs[2];
    GiNaC::ex C = Coeffs[1];
    GiNaC::ex D = Coeffs[1];

    GiNaC::ex Delta0 = B*B - 3 * A * C;
    GiNaC::ex Delta1 = 2 * B*B*B - 9 * A * B * C + 27 * A * A * D;
    GiNaC::ex CD = Delta1 + GiNaC::sqrt(Delta1 * Delta1 - 4 * GiNaC::pow(Delta0, 3));
    CD = CD/2;
    CD = GiNaC::pow(CD, GiNaC::numeric(1)/3);

    GiNaC::symbol U("u");
    GiNaC::ex Var = GiNaC::numeric(-1)/(3 * A) * (B + U * CD + Delta0/(U * CD));
    GiNaC::ex One   = Var.subs(U == 1);
    GiNaC::ex Two   = Var.subs(U == ((-1 + GiNaC::sqrt(GiNaC::numeric(-3)))/2));
    GiNaC::ex Three = Var.subs(U == ((-1 - GiNaC::sqrt(GiNaC::numeric(-3)))/2));
    if (GiNaC::is_a<GiNaC::numeric>(One))
      Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(One.evalf()).to_double());
    if (GiNaC::is_a<GiNaC::numeric>(Two))
      Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(Two.evalf()).to_double());
    if (GiNaC::is_a<GiNaC::numeric>(Three))
      Roots.push_back(GiNaC::ex_to<GiNaC::numeric>(Three.evalf()).to_double());
  }

  return Roots;
}

bool TestExpr::runOnModule(Module&) {
  Expr A("a");
  Expr B("b");

  assert(A + B == B + A);
  assert(A * B == B * A);
  assert(A - B != B - A);
  assert(A / B != B / A);

  //map<GiNaC::symbol, int> GetMaxExps(GiNaC::ex Ex); 
  GiNaC::symbol GA("a");
  GiNaC::symbol GB("b");
  GiNaC::ex EZ = 5 * GA * GA * GA + 4 * GB + 1;

  assert(EZ.is_polynomial(GA));
  assert(EZ.degree(GA) == 3);
  assert(EZ.coeff(GA, 3) == GiNaC::numeric(5));
  assert(EZ.coeff(GA, 2) == GiNaC::numeric(0));
  assert(EZ.coeff(GA, 1) == GiNaC::numeric(0));
  assert(EZ.coeff(GA, 0) == 4 * GB + 1);

  assert(EZ.is_polynomial(GB));
  assert(EZ.degree(GB) == 1);
  assert(EZ.coeff(GB, 1) == GiNaC::numeric(4));
  assert(EZ.coeff(GB, 0) == 5 * GA * GA * GA + 1);

  auto SEZ = GetSymbols(EZ);
  assert(SEZ.size() == 2);
  assert(SEZ.count(GA) && SEZ.count(GB));

  auto CEZGA = GetCoeffs(EZ, GA);
  assert(CEZGA.size() == 4);
  assert(CEZGA[3] == GiNaC::numeric(5));
  assert(CEZGA[2] == GiNaC::numeric(0));
  assert(CEZGA[1] == GiNaC::numeric(0));
  assert(CEZGA[0] == 4 * GB + 1);

  auto CEZGB = GetCoeffs(EZ, GB);
  assert(CEZGB.size() == 2);
  assert(CEZGB[1] == GiNaC::numeric(4));
  assert(CEZGB[0] == 5 * GA * GA * GA + 1);

  auto EF = GA * GA - 1;
  auto EFSolve = Solve(EF, GA);
  assert(EFSolve.size() == 2);
  assert(Round(EFSolve[0]) == -1 || Round(EFSolve[1]) == -1);
  assert(Round(EFSolve[0]) == 1  || Round(EFSolve[1]) == 1);
  auto EFNegs = NegativeOrdinate(EFSolve, GA, EF); 
  assert(EFNegs.empty());

  EF = GA * GA - 2 * GA + 1;
  EFSolve = Solve(EF, GA);
  assert(EFSolve.size() == 2);
  assert(Round(EFSolve[0]) == 1 && Round(EFSolve[1]) == 1);
  EFNegs = NegativeOrdinate(EFSolve, GA, EF); 
  assert(EFNegs.empty());

  EF = (-1) * GA * GA + GA - 1;
  EFSolve = Solve(EF, GA);
  assert(EFSolve.size() == 0);
  EFNegs = NegativeOrdinate(EFSolve, GA, EF); 
  assert(EFNegs.size() == 1);
  assert(EFNegs[0].first.is_equal(GiNaC::inf(-1)));
  assert(EFNegs[0].second.is_equal(GiNaC::inf(1)));

  return false;
}