symengine · rikardn · May 28, 2021 · May 30, 2021
diff --git a/symengine/CMakeLists.txt b/symengine/CMakeLists.txt
@@ -55,6 +55,7 @@ set(SRC
     functions.cpp
     infinity.cpp
     integer.cpp
+    integration.cpp
     logic.cpp
     matrix.cpp
     monomials.cpp
@@ -165,6 +166,7 @@ set(HEADERS
     functions.h
     infinity.h
     integer.h
+    integration.h
     lambda_double.h
     llvm_double.h
     logic.h

diff --git a/symengine/derivative.cpp b/symengine/derivative.cpp
@@ -192,6 +192,7 @@ void DiffVisitor::bvisit(const Basic &self)
 DIFF0(UnivariateSeries)
 DIFF0(Max)
 DIFF0(Min)
+DIFF0(Integral)
 #endif
 
 void DiffVisitor::bvisit(const Number &self)

diff --git a/symengine/derivative.h b/symengine/derivative.h
@@ -46,6 +46,7 @@ class DiffVisitor : public BaseVisitor<DiffVisitor>
     void bvisit(const UnivariateSeries &self);
     void bvisit(const Max &self);
     void bvisit(const Min &self);
+    void bvisit(const Integral &self);
 #endif
     void bvisit(const Number &self);
     void bvisit(const Constant &self);

diff --git a/symengine/functions.cpp b/symengine/functions.cpp
@@ -2025,6 +2025,48 @@ int Derivative::compare(const Basic &o) const
     return cmp;
 }
 
+Integral::Integral(const RCP<const Basic> &arg, const RCP<const Basic> &x)
+    : arg_{arg}, x_{x}
+{
+    SYMENGINE_ASSIGN_TYPEID()
+    SYMENGINE_ASSERT(is_canonical(arg, x))
+}
+
+bool Integral::is_canonical(const RCP<const Basic> &arg,
+                            const RCP<const Basic> &x) const
+{
+    if (not is_a<Symbol>(*x))
+        return false;
+    return true;
+}
+
+hash_t Integral::__hash__() const
+{
+    hash_t seed = SYMENGINE_INTEGRAL;
+    hash_combine<Basic>(seed, *arg_);
+    hash_combine<Basic>(seed, *x_);
+    return seed;
+}
+
+bool Integral::__eq__(const Basic &o) const
+{
+    if (is_a<Integral>(o) and eq(*arg_, *(down_cast<const Integral &>(o).arg_))
+        and eq(*x_, *(down_cast<const Integral &>(o).x_)))
+        return true;
+    return false;
+}
+
+int Integral::compare(const Basic &o) const
+{
+    SYMENGINE_ASSERT(is_a<Integral>(o))
+    const Integral &s = down_cast<const Integral &>(o);
+    int cmp = arg_->__cmp__(*(s.arg_));
+    if (cmp != 0)
+        return cmp;
+    cmp = x_->__cmp__(*(s.x_));
+    return cmp;
+}
+
 // Subs class
 Subs::Subs(const RCP<const Basic> &arg, const map_basic_basic &dict)
     : arg_{arg}, dict_{dict}

diff --git a/symengine/functions.h b/symengine/functions.h
@@ -731,6 +731,46 @@ class Derivative : public Basic
                       const multiset_basic &x) const;
 };
 
+/*! Integral operator
+ *  Integral(f, x) represents the integral of `f` with respect to
+ *  `x`.
+ * */
+class Integral : public Basic
+{
+private:
+    RCP<const Basic> arg_;
+    RCP<const Basic> x_;
+
+public:
+    IMPLEMENT_TYPEID(SYMENGINE_INTEGRAL)
+    Integral(const RCP<const Basic> &arg, const RCP<const Basic> &x);
+
+    static RCP<const Integral> create(const RCP<const Basic> &arg,
+                                      const RCP<const Basic> &x)
+    {
+        return make_rcp<const Integral>(arg, x);
+    }
+
+    virtual hash_t __hash__() const;
+    virtual bool __eq__(const Basic &o) const;
+    virtual int compare(const Basic &o) const;
+    inline RCP<const Basic> get_arg() const
+    {
+        return arg_;
+    }
+    inline RCP<const Basic> get_symbol() const
+    {
+        return x_;
+    }
+    virtual vec_basic get_args() const
+    {
+        vec_basic args = {arg_, x_};
+        return args;
+    }
+    bool is_canonical(const RCP<const Basic> &arg,
+                      const RCP<const Basic> &x) const;
+};
+
 /*! Subs operator
  *  Subs(f, {x1 : x2, y1: y2, ...}) represents `f` after substituting
  *  `x1` with `x2`, `y1` with `y2`, and so on.

diff --git a/symengine/integration.cpp b/symengine/integration.cpp
@@ -0,0 +1,70 @@
+#include <symengine/visitor.h>
+#include <symengine/subs.h>
+#include <symengine/symengine_casts.h>
+#include <symengine/integration.h>
+
+namespace SymEngine
+{
+
+// This is integration based on Rubi-5
+// Rubi-5 is not a finished product, but a prototype is available
+// as of March 2021 at https://github.com/RuleBasedIntegration/Rubi-5
+// Rubi-5 will rely on less than 50 pattern matching rules and below these
+// a tree of if-then-else functions.
+// The plan is to manually implement the pattern matching rules in symengine
+// and to automatically generate code for the if-then-else functions from
+// the Rubi-5 mathematica code.
+
+RCP<const Basic> integrate(const RCP<const Basic> &arg,
+                           const RCP<const Symbol> &x)
+{
+    if (not has_symbol(*arg, *x)) {
+        // Int[a_,x_Symbol] := a*x /; FreeQ[a,x]
+        return mul(arg, x);
+    }
+    if (is_a<Mul>(*arg)) {
+        // Int[a_*Fx_,x_Symbol] := a \[Star] Int[Fx,x] /; FreeQ[a,x]
+        map_basic_basic factor_out;
+        map_basic_basic keep;
+        for (const auto &p : down_cast<const Mul &>(*arg).get_dict()) {
+            if (has_symbol(*p.first, *x) or has_symbol(*p.second, *x)) {
+                keep.insert(p);
+            } else {
+                factor_out.insert(p);
+            }
+        }
+        // Could be optimized by remembering factor_out and continuing here
+        return mul(Mul::from_dict(down_cast<const Mul &>(*arg).get_coef(),
+                                  std::move(factor_out)),
+                   integrate(Mul::from_dict(one, std::move(keep)), x));
+    }
+    if (is_a<Pow>(*arg)) {
+        // Int[(a_./x_)^p_,x_Symbol] := -a*(a/x)^(p-1)/(p-1) /;
+        // FreeQ[{a,p},x] && Not[IntegerQ[p]]
+        auto base = down_cast<const Pow &>(*arg).get_base();
+        auto p = down_cast<const Pow &>(*arg).get_exp();
+        if (is_a<Mul>(*base) and not has_symbol(*p, *x)
+            and not is_a<Integer>(*p)) {
+            map_basic_basic factor_out;
+            bool found = false;
+            for (const auto &pair : down_cast<const Mul &>(*base).get_dict()) {
+                if (eq(*pair.first, *x) and eq(*pair.second, *integer(-1))) {
+                    found = true;
+                } else {
+                    factor_out.insert(pair);
+                }
+            }
+            if (found) {
+                auto a
+                    = Mul::from_dict(down_cast<const Mul &>(*base).get_coef(),
+                                     std::move(factor_out));
+                auto pm1 = sub(p, integer(1));
+                return div(mul(mul(a, integer(-1)), pow(base, pm1)), pm1);
+            }
+        }
+    }
+
+    return make_rcp<const Integral>(arg, x);
+}
+
+} // namespace SymEngine
diff --git a/symengine/integration.h b/symengine/integration.h
@@ -0,0 +1,15 @@
+#ifndef SYMENGINE_INTEGRATION_H
+#define SYMENGINE_INTEGRATION_H
+
+#include <symengine/basic.h>
+#include <symengine/visitor.h>
+
+namespace SymEngine
+{
+
+RCP<const Basic> integrate(const RCP<const Basic> &arg,
+                           const RCP<const Symbol> &x);
+
+} // namespace SymEngine
+
+#endif // SYMENGINE_INTEGRATION_H
diff --git a/symengine/tests/basic/CMakeLists.txt b/symengine/tests/basic/CMakeLists.txt
@@ -134,3 +134,7 @@ add_test(test_simplify ${PROJECT_BINARY_DIR}/test_simplify)
 add_executable(test_tuple test_tuple.cpp)
 target_link_libraries(test_tuple symengine catch)
 add_test(test_tuple ${PROJECT_BINARY_DIR}/test_tuple)
+
+add_executable(test_integration test_integration.cpp)
+target_link_libraries(test_integration symengine catch)
+add_test(test_integration ${PROJECT_BINARY_DIR}/test_integration)
diff --git a/symengine/tests/basic/test_integration.cpp b/symengine/tests/basic/test_integration.cpp
@@ -0,0 +1,50 @@
+#include "catch.hpp"
+#include <symengine/basic.h>
+#include <symengine/integration.h>
+
+using SymEngine::integer;
+using SymEngine::Integral;
+using SymEngine::integrate;
+using SymEngine::make_rcp;
+using SymEngine::RCP;
+using SymEngine::Symbol;
+using SymEngine::symbol;
+
+TEST_CASE("Test Integral", "[Integral]")
+{
+    RCP<const Symbol> x = symbol("x");
+    RCP<const Symbol> y = symbol("y");
+    auto unsolved1 = integrate(sin(x), x);
+    auto unsolved2 = integrate(sin(y), y);
+
+    REQUIRE(unsolved1->compare(*unsolved2) == -1);
+    REQUIRE(unsolved1->__eq__(*unsolved1));
+    REQUIRE(not unsolved1->__eq__(*unsolved2));
+    REQUIRE(not unsolved1->__eq__(*integer(23)));
+}
+
+TEST_CASE("Test integrate", "[integrate]")
+{
+    RCP<const Symbol> x = symbol("x");
+    RCP<const Symbol> a = symbol("a");
+    RCP<const Symbol> p = symbol("p");
+
+    // Integrals fully solved by symengine
+    REQUIRE(eq(*integrate(a, x), *mul(a, x)));
+    REQUIRE(eq(*integrate(sin(a), x), *mul(sin(a), x)));
+    auto res = integrate(pow(div(a, x), p), x);
+    auto correct
+        = div(mul(mul(a, integer(-1)), pow(div(a, x), sub(p, integer(1)))),
+              sub(p, integer(1)));
+    REQUIRE(eq(*res, *correct));
+
+    // Integrals partially solved by symengine
+    REQUIRE(eq(*integrate(mul(x, a), x), *mul(a, make_rcp<Integral>(x, x))));
+    REQUIRE(eq(*integrate(mul(mul(integer(2), x), a), x),
+               *mul(mul(integer(2), a), make_rcp<Integral>(x, x))));
+
+    // Integrals not yet solved by symengine
+    REQUIRE(eq(*integrate(sin(x), x), Integral(sin(x), x)));
+    REQUIRE(eq(*integrate(sin(x), x)->get_args()[0], *sin(x)));
+    REQUIRE(eq(*integrate(sin(x), x)->get_args()[1], *x));
+}
diff --git a/symengine/tests/printing/test_printing.cpp b/symengine/tests/printing/test_printing.cpp
@@ -729,6 +729,7 @@ TEST_CASE("test_latex_printing()", "[latex]")
     CHECK(latex(*l10) == "\\left[-3, 3\\right]");
     CHECK(latex(*l11) == "\\mathrm{True}");
     CHECK(latex(*l12) == "\\mathrm{False}");
+    // FIXME: These tests are unstable between updates of symengine
     // CHECK(latex(*l13) == "5 \\leq b \\wedge 2 \\leq a");
     //    CHECK(latex(*l14)
     //          == "b \\leq a \\wedge \\left(a \\neq c \\vee a = b\\right)");

diff --git a/symengine/type_codes.inc b/symengine/type_codes.inc
@@ -89,6 +89,7 @@ SYMENGINE_ENUM(SYMENGINE_BETA, Beta)
 SYMENGINE_ENUM(SYMENGINE_FUNCTIONSYMBOL, FunctionSymbol)
 SYMENGINE_ENUM(SYMENGINE_FUNCTIONWRAPPER, FunctionWrapper)
 SYMENGINE_ENUM(SYMENGINE_DERIVATIVE, Derivative)
+SYMENGINE_ENUM(SYMENGINE_INTEGRAL, Integral)
 SYMENGINE_ENUM(SYMENGINE_SUBS, Subs)
 SYMENGINE_ENUM(SYMENGINE_ABS, Abs)
 SYMENGINE_ENUM(SYMENGINE_MAX, Max)