[Feature] Add integration by parts

src/Exponent.cpp

@@ -51,6 +51,18 @@ auto Exponent<Expression>::Integrate(const Expression& integrationVariable) cons
                 return adder.Accept(simplifyVisitor).value();
             }
         }
+
+        // EulerNumber raised to a (degree-one) Variable
+        if (auto eulerBase = RecursiveCast<Exponent<EulerNumber, Variable>>(*simplifiedExponent); eulerBase != nullptr) {
+            const Variable& expPow = eulerBase->GetLeastSigOp();
+
+            Add adder {
+                Exponent { EulerNumber {}, expPow },
+                Variable { "C" }
+            };
+
+            return adder.Accept(simplifyVisitor).value();
+        }
     }
 
     Integral<Expression, Expression> integral { *(this->Copy()), *(integrationVariable.Copy()) };

src/Multiply.cpp

@@ -5,11 +5,14 @@
 #include "Oasis/Multiply.hpp"
 #include "Oasis/Add.hpp"
 #include "Oasis/Divide.hpp"
+#include "Oasis/EulerNumber.hpp"
 #include "Oasis/Exponent.hpp"
 #include "Oasis/Imaginary.hpp"
 #include "Oasis/Integral.hpp"
+#include "Oasis/Log.hpp"
 #include "Oasis/Matrix.hpp"
 #include "Oasis/RecursiveCast.hpp"
+#include "Oasis/Subtract.hpp"
 
 #define EPSILON 10E-6
 
@@ -42,7 +45,101 @@ auto Multiply<Expression>::Integrate(const Expression& integrationVariable) cons
             }
         }
 
-        // TODO: Implement integration by parts
+        // Apply Integration By Parts
+        else {
+            // Start by assuming u is MostSigOp and dv is LeastSigOp
+            auto multiplyCopy = std::make_unique<Multiply<Expression, Expression>>(Multiply { this->GetMostSigOp(), this->GetLeastSigOp() });
+
+            // Check by LIPET if it's necessary to swap operands so u is LeastSigOp and dv is MostSigOp
+            // LIPET: Logarithm case
+            if (!(multiplyCopy->GetMostSigOp().Is<Log<Expression, Expression>>())
+                && multiplyCopy->GetLeastSigOp().Is<Log<Expression, Expression>>()) {
+                multiplyCopy = std::make_unique<Multiply<Expression, Expression>>(multiplyCopy->SwapOperands());
+            }
+
+            // LIPET: Inverse trigonometry is not implemented yet in Oasis
+
+            // LIPET: Swap polynomial (linear variable) and EulerNumber
+            if ((multiplyCopy->GetLeastSigOp().Is<Variable>()
+                    && multiplyCopy->GetMostSigOp().Is<Exponent<EulerNumber, Expression>>())) {
+                multiplyCopy = std::make_unique<Multiply<Expression, Expression>>(multiplyCopy->SwapOperands());
+            }
+
+            // LIPET: Swap polynomial (exponential) and EulerNumber
+            if (auto polynomial = RecursiveCast<Exponent<Variable, Real>>(multiplyCopy->GetLeastSigOp()); polynomial != nullptr) {
+                if (auto euler = RecursiveCast<Exponent<EulerNumber, Variable>>(multiplyCopy->GetMostSigOp()); euler != nullptr) {
+                    multiplyCopy = std::make_unique<Multiply<Expression, Expression>>(multiplyCopy->SwapOperands());
+                }
+            }
+
+            // LIPET: Trigonometry is not implemented yet in Oasis
+
+            // Integrate dv and differentiate u to attain v and du
+            auto v = multiplyCopy->GetLeastSigOp().Integrate(integrationVariable);
+            auto du = multiplyCopy->GetMostSigOp().Differentiate(integrationVariable);
+
+            // Remove the +C from v
+            auto vPlusC = RecursiveCast<Add<Expression, Variable>>(*v);
+            v = vPlusC->GetMostSigOp().Copy();
+
+            // Initialize the coefficients of v and du to 1
+            auto vCoefficient = std::make_unique<Real>(Real { 1 });
+            auto duCoefficient = std::make_unique<Real>(Real { 1 });
+
+            // Make a copy of v so that the coefficient can be factored out when computing vdu,
+            // while preserving the original v for the multiplication of u*v
+            auto v_VDU = v->Copy();
+
+            // Factor out the coefficient for v; v_VDU now has a coefficient of one
+            if (auto coefficientMultiply = RecursiveCast<Multiply<Real, Expression>>(*v); coefficientMultiply != nullptr) {
+                vCoefficient = RecursiveCast<Real>(coefficientMultiply->GetMostSigOp());
+                v_VDU = coefficientMultiply->GetLeastSigOp().Accept(simplifyVisitor).value();
+            }
+
+            // Factor out the coefficient for du; du now has a coefficient of one
+            if (auto coefficientMultiply = RecursiveCast<Multiply<Real, Expression>>(*du); coefficientMultiply != nullptr) {
+                duCoefficient = RecursiveCast<Real>(coefficientMultiply->GetMostSigOp());
+                du = coefficientMultiply->GetLeastSigOp().Accept(simplifyVisitor).value();
+            }
+
+            // Attain a constant for the coefficient of vdu
+            auto vduCoefficient = std::make_unique<Multiply<Expression, Expression>>(Multiply { *vCoefficient, *duCoefficient })->Accept(simplifyVisitor).value();
+
+            // Multiply v and du to attain partVDU, which is vdu with a coefficient of one
+            auto partVDU = std::make_unique<Multiply<Expression, Expression>>(Multiply { *v_VDU, *du })->Accept(simplifyVisitor).value();
+
+            // Prevent infinite recursion
+            // If fullVDU is equal to the original Multiply instance,
+            // then IBP cannot reduce the integrand, so return an integral of the Multiply instance
+            auto fullVDU = std::make_unique<Multiply<Expression, Expression>>(Multiply { *vduCoefficient, *partVDU });
+            if (multiplyCopy->Equals(*fullVDU)) {
+                Integral<Expression, Expression> integral { *(multiplyCopy->Copy()), *(integrationVariable.Copy()) };
+                return integral.Copy();
+            }
+
+            // If partVDU is a Multiply instance it will use IBP again,
+            // otherwise it will find another integration method to use.
+            // Implicit base case: IBP will eventually reduce the integrand such that it is not a Multiply instance.
+            auto integratedPartVDU = partVDU->Integrate(integrationVariable);
+
+            // Remove the +C from the integration result
+            if (auto removeC = RecursiveCast<Add<Expression, Expression>>(*integratedPartVDU); removeC != nullptr) {
+                integratedPartVDU = removeC->GetMostSigOp().Copy();
+            }
+
+            // Apply coefficient correction for integratedPartVDU, thus building integratedFullVDU
+            auto integratedFullVDU = std::make_unique<Multiply<Expression, Expression>>(Multiply { *vduCoefficient, *integratedPartVDU });
+
+            // Apply the Integration By Parts formula, and add the +C at the end
+            Add<Subtract<Multiply<Expression, Expression>, Expression>, Variable> adder {
+                Subtract<Multiply<Expression, Expression>, Expression> {
+                    Multiply<Expression, Expression> { multiplyCopy->GetMostSigOp(), *v },
+                    *integratedFullVDU },
+                Variable { "C" }
+            };
+
+            return adder.Accept(simplifyVisitor).value();
+        }
     }
     Integral<Expression, Expression> integral { *(this->Copy()), *(integrationVariable.Copy()) };
 

tests/IntegrateTests.cpp

@@ -2,7 +2,9 @@
 
 #include "Oasis/Add.hpp"
 #include "Oasis/Divide.hpp"
+#include "Oasis/EulerNumber.hpp"
 #include "Oasis/Exponent.hpp"
+#include "Oasis/Log.hpp"
 #include "Oasis/Multiply.hpp"
 #include "Oasis/Real.hpp"
 #include "Oasis/Subtract.hpp"
@@ -181,4 +183,160 @@ TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]")
     auto simplified = integrated->Accept(simplifyVisitor).value();
 
     REQUIRE(simplified->Equals(*(integral.Accept(simplifyVisitor).value())));
+}
+
+TEST_CASE("Integrate EulerNumber Raised to Power x", "[Integrate][Euler]")
+{
+    Oasis::Variable var { "x" };
+
+    Oasis::Exponent integrand { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } };
+
+    Oasis::Add<Oasis::Exponent<Oasis::EulerNumber, Oasis::Variable>, Oasis::Variable> integral {
+        Oasis::Add {
+            Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } },
+            Oasis::Variable { "C" }
+        }
+    };
+
+    auto ptr = integral.Accept(simplifyVisitor).value();
+    auto integrated = integrand.Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+}
+
+TEST_CASE("Integration By Parts: Variable and Euler's Number", "[Integrate][Variable][Euler]")
+{
+    Oasis::Variable var { "x" };
+
+    Oasis::Multiply<Oasis::Variable, Oasis::Exponent<Oasis::EulerNumber, Oasis::Variable>> integrand {
+        Oasis::Variable { var.GetName() },
+        Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }
+    };
+
+    Oasis::Add<Oasis::Subtract<Oasis::Multiply<Oasis::Variable, Oasis::Exponent<Oasis::EulerNumber, Oasis::Variable>>, Oasis::Exponent<Oasis::EulerNumber, Oasis::Variable>>, Oasis::Variable> integral {
+        Oasis::Add {
+            Oasis::Subtract {
+                Oasis::Multiply {
+                    Oasis::Variable { var.GetName() },
+                    Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } },
+                Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } },
+            Oasis::Variable { "C" }
+        }
+    };
+
+    auto ptr = integral.Accept(simplifyVisitor).value();
+    auto integrated = integrand.Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+
+    integrated = integrand.SwapOperands().Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+}
+
+TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Exponent][Euler]")
+{
+    Oasis::Variable var { "x" };
+
+    Oasis::Multiply<Oasis::Exponent<Oasis::Variable, Oasis::Real>, Oasis::Exponent<Oasis::EulerNumber, Oasis::Variable>> integrand {
+        Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } },
+        Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }
+    };
+
+    Oasis::Add<Oasis::Multiply<Oasis::Exponent<Oasis::EulerNumber, Oasis::Variable>,
+    Oasis::Add<Oasis::Exponent<Oasis::Variable, Oasis::Real>, Oasis::Multiply<Oasis::Real,
+    Oasis::Add<Oasis::Variable, Oasis::Real>>>>, Oasis::Variable> integral {
+        Oasis::Add {
+            Oasis::Multiply {
+                Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } },
+                Oasis::Add {
+                    Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } },
+                    Oasis::Multiply {
+                        Oasis::Real {-2},
+                        Oasis::Add {
+                            Oasis::Variable { var.GetName() },
+                            Oasis::Real { -1 }
+                        }
+                    },
+                }
+            },
+            Oasis::Variable { "C" }
+        }
+    };
+
+
+    auto ptr = integral.Accept(simplifyVisitor).value();
+    auto integrated = integrand.Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+
+    integrated = integrand.SwapOperands().Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+
+}
+
+TEST_CASE("Integration By Parts: Variable and Logarithm", "[Integrate][Variable][Logarithm]")
+{
+    Oasis::Variable var { "x" };
+
+    Oasis::Multiply<Oasis::Variable, Oasis::Log<Oasis::EulerNumber, Oasis::Variable>> integrand {
+        Oasis::Variable { var.GetName() },
+        Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }
+    };
+
+    Oasis::Add<Oasis::Multiply<Oasis::Exponent<Oasis::Variable, Oasis::Real>,
+    Oasis::Add<Oasis::Multiply< Oasis::Divide<Oasis::Real, Oasis::Real>,
+    Oasis::Log<Oasis::EulerNumber, Oasis::Variable>>, Oasis::Divide<Oasis::Real, Oasis::Real>>>, Oasis::Variable> integral {
+        Oasis::Add {
+            Oasis::Multiply {
+                Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 }  },
+                Oasis::Add {
+                    Oasis::Multiply {
+                        Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 2 } },
+                        Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }
+                    },
+                    Oasis::Divide { Oasis::Real { -1 }, Oasis::Real { 4 } },
+                 }
+            },
+            Oasis::Variable { "C" }
+        }
+    };
+
+    auto ptr = integral.Accept(simplifyVisitor).value();
+    auto integrated = integrand.Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+
+    integrated = integrand.SwapOperands().Integrate(var);
+    REQUIRE((integrated->Equals(*ptr)));
+}
+
+TEST_CASE("Integration By Parts: Exponent and Logarithm", "[Integrate][Exponent][Logarithm]")
+{
+    Oasis::Variable var { "x" };
+
+    Oasis::Multiply<Oasis::Exponent<Oasis::Variable, Oasis::Real>, Oasis::Log<Oasis::EulerNumber, Oasis::Variable>> integrand {
+        Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } },
+        Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }
+    };
+
+    Oasis::Add<Oasis::Multiply<Oasis::Exponent<Oasis::Variable, Oasis::Real>,
+    Oasis::Add<Oasis::Multiply< Oasis::Divide<Oasis::Real, Oasis::Real>,
+    Oasis::Log<Oasis::EulerNumber, Oasis::Variable>>, Oasis::Divide<Oasis::Real, Oasis::Real>>>, Oasis::Variable> integral {
+        Oasis::Add {
+            Oasis::Multiply {
+                Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 }  },
+                Oasis::Add {
+                   Oasis::Multiply {
+                       Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 3 } },
+                       Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }
+                   },
+                   Oasis::Divide { Oasis::Real { -1 }, Oasis::Real { 9 } },
+                }
+            },
+            Oasis::Variable { "C" }
+        }
+    };
+
+     auto ptr = integral.Accept(simplifyVisitor).value();
+     auto integrated = integrand.Integrate(var);
+     REQUIRE((integrated->Equals(*ptr)));
+
+     integrated = integrand.SwapOperands().Integrate(var);
+     REQUIRE((integrated->Equals(*ptr)));
 }