Symbolic polynomial fraction (RobotLocomotion#9485)

Add RationalFunction

common/BUILD.bazel

@@ -201,6 +201,8 @@ drake_cc_library(
         "symbolic_monomial_util.h",
         "symbolic_polynomial.cc",
         "symbolic_polynomial.h",
+        "symbolic_rational_function.cc",
+        "symbolic_rational_function.h",
         "symbolic_simplification.cc",
         "symbolic_simplification.h",
         "symbolic_variable.cc",
@@ -982,6 +984,15 @@ drake_cc_googletest(
     ],
 )
 
+drake_cc_googletest(
+    name = "symbolic_rational_function_test",
+    deps = [
+        ":symbolic",
+        "//common/test_utilities:expect_throws_message",
+        "//common/test_utilities:symbolic_test_util",
+    ],
+)
+
 drake_cc_googletest(
     name = "symbolic_substitution_test",
     deps = [

common/symbolic.h

@@ -37,6 +37,7 @@
 #include "drake/common/symbolic_monomial.h"
 #include "drake/common/symbolic_monomial_util.h"
 #include "drake/common/symbolic_polynomial.h"
+#include "drake/common/symbolic_rational_function.h"
 #include "drake/common/symbolic_formula.h"
 #include "drake/common/symbolic_formula_visitor.h"
 #include "drake/common/symbolic_simplification.h"

common/symbolic_rational_function.cc

@@ -0,0 +1,208 @@
+// NOLINTNEXTLINE(build/include): Its header file is included in symbolic.h.
+
+#include "drake/common/symbolic.h"
+
+namespace drake {
+namespace symbolic {
+RationalFunction::RationalFunction()
+    : numerator_{} /* zero polynomial */, denominator_{1} {}
+
+RationalFunction::RationalFunction(const Polynomial& numerator,
+                                   const Polynomial& denominator)
+    : numerator_{numerator}, denominator_{denominator} {
+  if (denominator_.EqualTo(Polynomial() /* zero polynomial */)) {
+    throw std::invalid_argument(
+        "RationalFunction: the denominator should not be 0.");
+  }
+  CheckIndeterminates();
+}
+
+bool RationalFunction::EqualTo(const RationalFunction& f) const {
+  return numerator_.EqualTo(f.numerator()) &&
+         denominator_.EqualTo(f.denominator());
+}
+
+std::ostream& operator<<(std::ostream& os, const RationalFunction& f) {
+  os << "(" << f.numerator() << ") / (" << f.denominator() << ")";
+  return os;
+}
+
+void RationalFunction::CheckIndeterminates() const {
+  const Variables vars1{intersect(numerator_.indeterminates(),
+                                  denominator_.decision_variables())};
+  const Variables vars2{intersect(numerator_.decision_variables(),
+                                  denominator_.indeterminates())};
+  if (!vars1.empty() || !vars2.empty()) {
+    std::ostringstream oss;
+    oss << "RationalFunction " << *this << " is invalid.\n";
+    if (!vars1.empty()) {
+      oss << "The following variable(s) "
+             "are used as indeterminates in the numerator and decision "
+             "variables in the denominator at the same time:\n"
+          << vars1 << ".\n";
+    }
+    if (!vars2.empty()) {
+      oss << "The following variable(s) "
+             "are used as decision variables in the numerator and "
+             "indeterminates variables in the denominator at the same time:\n"
+          << vars2 << ".\n";
+    }
+    throw std::logic_error(oss.str());
+  }
+}
+
+RationalFunction& RationalFunction::operator+=(const RationalFunction& f) {
+  numerator_ = numerator_ * f.denominator() + denominator_ * f.numerator();
+  denominator_ *= f.denominator();
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator+=(const Polynomial& p) {
+  numerator_ = p * denominator_ + numerator_;
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator+=(double c) {
+  numerator_ = c * denominator_ + numerator_;
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator-=(const RationalFunction& f) {
+  *this += -f;
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator-=(const Polynomial& p) {
+  *this += -p;
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator-=(double c) { return *this += -c; }
+
+RationalFunction& RationalFunction::operator*=(const RationalFunction& f) {
+  numerator_ *= f.numerator();
+  denominator_ *= f.denominator();
+  CheckIndeterminates();
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator*=(const Polynomial& p) {
+  numerator_ *= p;
+  CheckIndeterminates();
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator*=(double c) {
+  numerator_ *= c;
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator/=(const RationalFunction& f) {
+  if (f.numerator().EqualTo(Polynomial())) {
+    throw std::logic_error("RationalFunction: operator/=: The divider is 0.");
+  }
+  numerator_ *= f.denominator();
+  denominator_ *= f.numerator();
+  CheckIndeterminates();
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator/=(const Polynomial& p) {
+  if (p.EqualTo(Polynomial())) {
+    throw std::logic_error("RationalFunction: operator/=: The divider is 0.");
+  }
+  denominator_ *= p;
+  CheckIndeterminates();
+  return *this;
+}
+
+RationalFunction& RationalFunction::operator/=(double c) {
+  if (c == 0) {
+    throw std::logic_error("RationalFunction: operator/=: The divider is 0.");
+  }
+  denominator_ *= c;
+  return *this;
+}
+
+RationalFunction operator-(RationalFunction f) {
+  return RationalFunction(-f.numerator(), f.denominator());
+}
+
+RationalFunction operator+(RationalFunction f1, const RationalFunction& f2) {
+  return f1 += f2;
+}
+
+RationalFunction operator+(RationalFunction f, const Polynomial& p) {
+  return f += p;
+}
+
+RationalFunction operator+(const Polynomial& p, RationalFunction f) {
+  return f += p;
+}
+
+RationalFunction operator+(RationalFunction f, double c) { return f += c; }
+
+RationalFunction operator+(double c, RationalFunction f) { return f += c; }
+
+RationalFunction operator-(RationalFunction f1, const RationalFunction& f2) {
+  return f1 -= f2;
+}
+
+RationalFunction operator-(RationalFunction f, const Polynomial& p) {
+  return f -= p;
+}
+
+RationalFunction operator-(const Polynomial& p, RationalFunction f) {
+  return f = -f + p;
+}
+
+RationalFunction operator-(RationalFunction f, double c) { return f -= c; }
+
+RationalFunction operator-(double c, RationalFunction f) { return f = -f + c; }
+
+RationalFunction operator*(RationalFunction f1, const RationalFunction& f2) {
+  return f1 *= f2;
+}
+
+RationalFunction operator*(RationalFunction f, const Polynomial& p) {
+  return f *= p;
+}
+
+RationalFunction operator*(const Polynomial& p, RationalFunction f) {
+  return f *= p;
+}
+
+RationalFunction operator*(RationalFunction f, double c) { return f *= c; }
+
+RationalFunction operator*(double c, RationalFunction f) { return f *= c; }
+
+RationalFunction operator/(RationalFunction f1, const RationalFunction& f2) {
+  return f1 /= f2;
+}
+
+RationalFunction operator/(RationalFunction f, const Polynomial& p) {
+  return f /= p;
+}
+
+RationalFunction operator/(const Polynomial& p, const RationalFunction& f) {
+  return RationalFunction(p * f.denominator(), f.numerator());
+}
+
+RationalFunction operator/(RationalFunction f, double c) { return f /= c; }
+
+RationalFunction operator/(double c, const RationalFunction& f) {
+  return RationalFunction(c * f.denominator(), f.numerator());
+}
+
+RationalFunction pow(const RationalFunction& f, int n) {
+  if (n == 0) {
+    return RationalFunction(Polynomial(1), Polynomial(1));
+  } else if (n >= 1) {
+    return RationalFunction(pow(f.numerator(), n), pow(f.denominator(), n));
+  } else {
+    // n < 0
+    return RationalFunction(pow(f.denominator(), -n), pow(f.numerator(), -n));
+  }
+}
+}  // namespace symbolic
+}  // namespace drake

common/symbolic_rational_function.h

@@ -0,0 +1,126 @@
+#pragma once
+
+#ifndef DRAKE_COMMON_SYMBOLIC_HEADER
+// TODO(soonho-tri): Change to #error, when #6613 merged.
+#warning Do not directly include this file. Include "drake/common/symbolic.h".
+#endif
+
+#include <ostream>
+
+#include "drake/common/symbolic.h"
+
+namespace drake {
+namespace symbolic {
+/**
+ * Represents symbolic rational function. A function f(x) is a rational
+ * function, if f(x) = p(x) / q(x), where both p(x) and q(x) are polynomials of
+ * x. Note that rational functions are closed under (+, -, x, /). One
+ * application of rational function is in polynomial optimization, where we
+ * represent (or approximate) functions using rational functions, and then
+ * convert the constraint f(x) = h(x) (where h(x) is a polynomial) to a
+ * polynomial constraint p(x) - q(x) * h(x) = 0, or convert the inequality
+ * constraint f(x) >= h(x) as p(x) - q(x) * h(x) >= 0 if we know q(x) > 0.
+ *
+ * This class represents a special subset of the symbolic::Expression. While a
+ * symbolic::Expression can represent a rational function, extracting the
+ * numerator and denominator, generally, is quite difficult; for instance, from
+ * p1(x) / q1(x) + p2(x) / q2(x) + ... + pn(x) / qn(x). This class's explicit
+ * structure facilitates this decomposition.
+ */
+class RationalFunction {
+ public:
+  /** Constructs a zero rational function 0 / 1. */
+  RationalFunction();
+
+  DRAKE_DEFAULT_COPY_AND_MOVE_AND_ASSIGN(RationalFunction)
+
+  /**
+   * Constructs the rational function: numerator / denominator.
+   * @param numerator The numerator of the fraction.
+   * @param denominator The denominator of the fraction.
+   * @pre denominator cannot be structurally equal to 0.
+   * @pre None of the indeterminates in the numerator can be decision variables
+   * in the denominator; similarly none of the indeterminates in the denominator
+   * can be decision variables in the numerator.
+   * @throw logic_error if the precondition is not satisfied.
+   */
+  RationalFunction(const Polynomial& numerator, const Polynomial& denominator);
+
+  ~RationalFunction() = default;
+
+  /// Getter for the numerator.
+  const Polynomial& numerator() const { return numerator_; }
+
+  /// Getter for the denominator.
+  const Polynomial& denominator() const { return denominator_; }
+
+  RationalFunction& operator+=(const RationalFunction& f);
+  RationalFunction& operator+=(const Polynomial& p);
+  RationalFunction& operator+=(double c);
+
+  RationalFunction& operator-=(const RationalFunction& f);
+  RationalFunction& operator-=(const Polynomial& p);
+  RationalFunction& operator-=(double c);
+
+  RationalFunction& operator*=(const RationalFunction& f);
+  RationalFunction& operator*=(const Polynomial& p);
+  RationalFunction& operator*=(double c);
+
+  RationalFunction& operator/=(const RationalFunction& f);
+  RationalFunction& operator/=(const Polynomial& p);
+  RationalFunction& operator/=(double c);
+
+  /**
+   * Returns true if this rational function and f are structurally equal.
+   */
+  bool EqualTo(const RationalFunction& f) const;
+
+  friend std::ostream& operator<<(std::ostream&, const RationalFunction& f);
+
+ private:
+  // Throws std::logic_error if an indeterminate of the denominator (numerator,
+  // respectively) is a decision variable of the numerator (denominator).
+  void CheckIndeterminates() const;
+  Polynomial numerator_;
+  Polynomial denominator_;
+};
+
+/**
+ * Unary minus operation for rational function.
+ * if f(x) = p(x) / q(x), then -f(x) = (-p(x)) / q(x)
+ */
+RationalFunction operator-(RationalFunction f);
+
+RationalFunction operator+(RationalFunction f1, const RationalFunction& f2);
+RationalFunction operator+(RationalFunction f, const Polynomial& p);
+RationalFunction operator+(const Polynomial& p, RationalFunction f);
+RationalFunction operator+(RationalFunction f, double c);
+RationalFunction operator+(double c, RationalFunction f);
+
+RationalFunction operator-(RationalFunction f1, const RationalFunction& f2);
+RationalFunction operator-(RationalFunction f, const Polynomial& p);
+RationalFunction operator-(const Polynomial& p, RationalFunction f);
+RationalFunction operator-(RationalFunction f, double c);
+RationalFunction operator-(double c, RationalFunction f);
+
+RationalFunction operator*(RationalFunction f1, const RationalFunction& f2);
+RationalFunction operator*(RationalFunction f, const Polynomial& p);
+RationalFunction operator*(const Polynomial& p, RationalFunction f);
+RationalFunction operator*(RationalFunction f, double c);
+RationalFunction operator*(double c, RationalFunction f);
+
+RationalFunction operator/(RationalFunction f1, const RationalFunction& f2);
+RationalFunction operator/(RationalFunction f, const Polynomial& p);
+RationalFunction operator/(const Polynomial& p, const RationalFunction& f);
+RationalFunction operator/(RationalFunction f, double c);
+RationalFunction operator/(double c, const RationalFunction& f);
+
+/**
+ * Returns the rational function @p f raised to @p n.
+ * If n is positive, (f/g)ⁿ = fⁿ / gⁿ;
+ * If n is negative, (f/g)ⁿ = g⁻ⁿ / f⁻ⁿ;
+ * (f/g)⁰ = 1 / 1.
+ */
+RationalFunction pow(const RationalFunction& f, int n);
+}  // namespace symbolic
+}  // namespace drake