add derivation feature

src/analysis/PolynomialP.hpp

@@ -217,7 +217,6 @@ class PolynomialP
     if constexpr (Dimension == 1) {
       abs_pos = relative_pos[0];
     } else if constexpr (Dimension == 2) {
-      size_t total_degree = relative_pos[0] + relative_pos[1];
       abs_pos = total_degree * (total_degree + 1) / 2 + relative_pos[1];
     } else {
       throw NotImplementedError("Not yet Available in 3D");
@@ -316,9 +315,97 @@ class PolynomialP
         }
         return Q;
       } else {
-        throw UnexpectedError("Not yet Available in 3D");
+        throw NotImplementedError("Not yet Available in 3D");
+      }
+    }
+  }
+
+  PUGS_INLINE
+  constexpr friend std::ostream&
+  operator<<(std::ostream& os, const PolynomialP<N, Dimension>& P)
+  {
+    //    os << "P(x) = ";
+    bool all_coef_zero = true;
+    if (N == 0) {
+      os << P.coefficients()[0];
+      return os;
+    }
+    if constexpr (Dimension == 1) {
+      if (N != 1) {
+        if (P.coefficients()[N] != 0.) {
+          if (P.coefficients()[N] < 0.) {
+            os << "- ";
+          }
+          if (P.coefficients()[N] != 1 && P.coefficients()[N] != -1) {
+            os << std::abs(P.coefficients()[N]);
+          }
+          os << "x^" << N;
+          all_coef_zero = false;
+        }
+      }
+      for (size_t i = N - 1; i > 1; --i) {
+        if (P.coefficients()[i] != 0.) {
+          if (P.coefficients()[i] > 0.) {
+            os << " + ";
+          } else if (P.coefficients()[i] < 0.) {
+            os << " - ";
+          }
+          if (P.coefficients()[i] != 1 && P.coefficients()[i] != -1) {
+            os << std::abs(P.coefficients()[i]);
+          }
+          os << "x^" << i;
+          all_coef_zero = false;
+        }
+      }
+      if (P.coefficients()[1] != 0.) {
+        if (P.coefficients()[1] > 0. && N != 1) {
+          os << " + ";
+        } else if (P.coefficients()[1] < 0.) {
+          os << " - ";
+        }
+        if (P.coefficients()[1] != 1 && P.coefficients()[1] != -1) {
+          os << std::abs(P.coefficients()[1]);
+        }
+        os << "x";
+        all_coef_zero = false;
+      }
+      if (P.coefficients()[0] != 0. || all_coef_zero) {
+        if (P.coefficients()[0] > 0.) {
+          os << " + ";
+        } else if (P.coefficients()[0] < 0.) {
+          os << " - ";
+        }
+        os << std::abs(P.coefficients()[0]);
+      }
+      return os;
+    } else if constexpr (Dimension == 2) {
+      for (size_t j = N; j > 0; --j) {
+        for (size_t i = N - j; i > 0; --i) {
+          double coef = P[i, j];
+          if (coef != 0.) {
+            if (coef > 0.) {
+              os << " + ";
+            } else if (coef < 0.) {
+              os << " - ";
+            }
+            if (coef != 1 && coef != -1) {
+              os << std::abs(coef);
+            }
+            if (i == 0) {
+              os << "y" << j;
+            } else if (j == 0) {
+              os << "x^" << i;
+            } else {
+              os << "x^" << i << "y^" << j;
+            }
+            all_coef_zero = false;
+          }
         }
       }
+      return os;
+    } else {
+      throw NotImplementedError("Not yet Available in 3D");
+    }
   }
 
   PUGS_INLINE constexpr PolynomialP() noexcept = default;