From 376a7534fced02722eb3d157954d8d2100d9f937 Mon Sep 17 00:00:00 2001 From: bguyott Date: Thu, 12 Feb 2026 20:24:32 -0500 Subject: [PATCH 01/21] Started work on implementing IBP Implemented a system to follow LIPET rules, differentiated u and integrated dv, started work on the tree that will perform the integration by parts according to the formula. No test cases yet. --- src/Integral.cpp | 2 ++ src/Multiply.cpp | 83 ++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+) diff --git a/src/Integral.cpp b/src/Integral.cpp index ac901ad8..92896833 100644 --- a/src/Integral.cpp +++ b/src/Integral.cpp @@ -124,6 +124,8 @@ auto Integral::Simplify(const Expression& upper, const Expression& l // Integration by Parts // Needs Trig and Euler's Number + // Work in Progress + // Partial Fraction Decomposition // Work in Progress diff --git a/src/Multiply.cpp b/src/Multiply.cpp index cbebebc3..746d1dea 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -5,12 +5,16 @@ #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/Matrix.hpp" #include "Oasis/RecursiveCast.hpp" +#include "Oasis/Log.hpp" +#include "Oasis/Subtract.hpp" + #define EPSILON 10E-6 namespace Oasis { @@ -43,6 +47,85 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons } // TODO: Implement integration by parts + // Detect whether integration by parts is appropriate + // May need to simplify before and/or after + else if (auto mult = RecursiveCast>(*simplifiedMult); mult != nullptr) { + std::unique_ptr u; + std::unique_ptr dv; + + // Check the rules of LIPET + + // May need to recursively cast u and dv before assigning them + // so that method overrides can be realized. + + // TODO: change dv to v and muliply by integration variable to correctly attain dv + + // Logarithm + if (mult->GetMostSigOp().Is>()) { + u = mult->GetMostSigOp().Copy(); + dv = mult->GetLeastSigOp().Copy(); + } else if (mult->GetLeastSigOp().Is>()) { + u = mult->GetLeastSigOp().Copy(); + dv = mult->GetMostSigOp().Copy(); + } + + // TODO: Inverse trigonometry + // Inverse trigonometry is not implemented yet in Oasis + + // Polynomial + // TODO: Could also be exponent, in the case of (x^2)*sinx + // TODO: Ensure all polynomial cases are accounted for + if (mult->GetMostSigOp().Is()) { + u = mult->GetMostSigOp().Copy(); + dv = mult->GetLeastSigOp().Copy(); + } else if (mult->GetLeastSigOp().Is()) { + u = mult->GetLeastSigOp().Copy(); + dv = mult->GetMostSigOp().Copy(); + } + + // Exponential - Euler's Number + if (mult->GetMostSigOp().Is>()) { + u = mult->GetMostSigOp().Copy(); + dv = mult->GetLeastSigOp().Copy(); + } else if (mult->GetLeastSigOp().Is>()) { + u = mult->GetLeastSigOp().Copy(); + dv = mult->GetMostSigOp().Copy(); + } + + // TODO: Trigonometry + // Trigonometry is not implemented yet in Oasis + + // Differentiate u and integrate dv + // TODO: ensure u & dv are casted as the correct types to allow method override + std::unique_ptr du = u->Differentiate(*u); + std::unique_ptr v = dv->Integrate(*dv); + + // Apply the formula: integral(udv) = uv - integral(vdu) + // Work in Progress + Subtract, Integral> subtractor { + Multiply { + *u, + *dv }, + Integral {} + }; + + + + // test code to base the structure of above tree off of + // Multiply subtract { + // Multiply { + // Real { 1.0 }, + // Real { 2.0 } }, + // Real { 3.0 } + // }; + + + + + + + } + } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; From 72052ce6d25a66976ffabb439d347240b215aedb Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 20 Feb 2026 17:57:40 -0500 Subject: [PATCH 02/21] Completed implementation of IBP formula + Corrections The implementation of the IBP formula: integral of udv = uv - integral of vdu. The function should return the correct value of the simplified subtraction done from IBP. I had also made several corrections to my prior code. I used the object's self-reference to attain its most and least significant operands rather than making a copy to ensure it remains a subclass of Expression and not an instance of Expression, where it would lose its method overrides which is critical for IBP. # Conflicts: # src/Multiply.cpp --- src/Multiply.cpp | 525 +++++++++++++++++++++++++++++++++++++++++++---- 1 file changed, 483 insertions(+), 42 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 746d1dea..d8d53673 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -18,6 +18,445 @@ #define EPSILON 10E-6 namespace Oasis { + +auto Multiply::Simplify() const -> std::unique_ptr +{ + auto simplifiedMultiplicand = mostSigOp->Simplify(); + auto simplifiedMultiplier = leastSigOp->Simplify(); + + SimplifyVisitor simplifyVisitor; + + Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier }; + if (auto onezerocase = RecursiveCast>(simplifiedMultiply); onezerocase != nullptr) { + const Real& multiplicand = onezerocase->GetMostSigOp(); + const Expression& multiplier = onezerocase->GetLeastSigOp(); + if (std::abs(multiplicand.GetValue()) <= EPSILON) { + return std::make_unique(Real { 0.0 }); + } + if (multiplicand.GetValue() == 1) { + return multiplier.Simplify(); + } + } + if (auto realCase = RecursiveCast>(simplifiedMultiply); realCase != nullptr) { + const Real& multiplicand = realCase->GetMostSigOp(); + const Real& multiplier = realCase->GetLeastSigOp(); + return std::make_unique(multiplicand.GetValue() * multiplier.GetValue()); + } + + if (auto ImgCase = RecursiveCast>(simplifiedMultiply); ImgCase != nullptr) { + return std::make_unique(-1.0); + } + + if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { + auto e = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; + auto m = e.Accept(simplifyVisitor); + if (!m) { + return Divide { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); + } + return Divide { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); + } + + if (auto exprCase = RecursiveCast>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { + return std::make_unique>(exprCase->GetMostSigOp(), Real { 2.0 }); + } + } + + if (auto rMatrixCase = RecursiveCast>(simplifiedMultiply); rMatrixCase != nullptr) { + return std::make_unique(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()); + } + + if (auto matrixCase = RecursiveCast>(simplifiedMultiply); matrixCase != nullptr) { + const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp(); + const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp(); + + if (leftTerm.GetCols() == rightTerm.GetRows()) { + return std::make_unique(leftTerm.GetMatrix() * rightTerm.GetMatrix()); + } else { + // ERROR: INVALID DIMENSION + return std::make_unique>(leftTerm, rightTerm); + } + } + + // Commented out to not cause massive problems with things that need factored expressions + // // c*(a-b) + // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { + // if (negated->GetMostSigOp().GetValue()<0){ + // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, + // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); + // } else if (negated->GetMostSigOp().GetValue()>0){ + // return Subtract{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, + // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); + // } else { + // return Real{0}.Copy(); + // } + // + // } + // + // // c*(a+b) + // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { + // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, + // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); + // } + + if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { + auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + if (!ms) { + return Divide { m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); + } + // This doesn't have to check errors because errors should have already been caught above + return Divide { *std::move(ms).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value(); + } + + if (auto multCase = RecursiveCast, Divide>>(simplifiedMultiply); multCase != nullptr) { + auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Multiply { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return this->Generalize(); + } + return Divide { *std::move(ms).value(), *std::move(ns).value() }.Accept(simplifyVisitor).value(); + } + + if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp(), e); + } + return std::make_unique>(exprCase->GetMostSigOp(), + *(std::move(s).value())); + } + } + + // x*x^n + if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp(), e); + } + return std::make_unique>(exprCase->GetMostSigOp(), *(std::move(s).value())); + } + } + + if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetLeastSigOp(), e); + } + return std::make_unique>(exprCase->GetLeastSigOp(), *(std::move(s).value())); + } + } + + // x^n*x^m + if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), e); + } + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s).value())); + } + } + + // a*x*x + if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }); + } + } + + // a*x*b*x + if (auto exprCase = RecursiveCast, Multiply>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return this->Generalize(); + } + return std::make_unique>( + *(std::move(ms).value()), + *(std::move(ns).value())); + } + } + + // a*x^n*x + if (auto exprCase = RecursiveCast>, Expression>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { + auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e }); + } + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + *(std::move(s).value()) }); + } + } + + // a*x*x^n + if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), e }); + } + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + *(std::move(s).value()) }); + } + if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), + Exponent { exprCase->GetMostSigOp().GetMostSigOp(), e }); + } + return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), + Exponent { exprCase->GetMostSigOp().GetMostSigOp(), + *(std::move(s).value()) }); + } + } + + // a*x^n*b*x + if (auto exprCase = RecursiveCast, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>( + m, Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), n }); + } + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + *(std::move(ns).value()) }); + } + } + + if (auto exprCase = RecursiveCast, Multiply, Expression>>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>(m, + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + n }); + } + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + *(std::move(ns).value()) }); + } + } + + if (auto exprCase = RecursiveCast>, Multiply>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>(m, + Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), + n }); + } + + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), + *(std::move(ns).value()) }); + } + } + + // a*x^n*x^m + if (auto exprCase = RecursiveCast>, Exponent>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>( + exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + e }); + } + return std::make_unique>( + exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + *(std::move(s).value()) }); + } + } + + // a*x^n*b*x^m + if (auto exprCase = RecursiveCast>, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>( + m, + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + n }); + } + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + *(std::move(ns).value()) }); + } + } + + // if (auto negate = RecursiveCast>>>(simplifiedMultiply); negate != nullptr){ + // return Add{Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetMostSigOp()}, + // Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetLeastSigOp()}}.Simplify(); + // } + + // multiply add like terms + std::vector> multiplies; + std::vector> vals; + simplifiedMultiply.Flatten(multiplies); + for (const auto& multiplicand : multiplies) { + size_t i = 0; + if (auto real = RecursiveCast(*multiplicand); real != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast(*vals[i]); valI != nullptr) { + vals[i] = Real { valI->GetValue() * real->GetValue() }.Generalize(); + break; + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + vals.push_back(multiplicand->Generalize()); + } + continue; + } + // single i + if (auto img = RecursiveCast(*multiplicand); img != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { Imaginary {}, e }.Generalize(); + } else { + vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); + } + + break; + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + vals.push_back(Exponent { Imaginary {}, Real { 1.0 } }.Generalize()); + } + continue; + } + // i^n + if (auto img = RecursiveCast>(*multiplicand); img != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + auto e = Add { valI->GetLeastSigOp(), img->GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { Imaginary {}, e }.Generalize(); + } else { + vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); + } + break; + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); + vals.push_back(img->Generalize()); + } + continue; + } + // expr^n + if (auto expr = RecursiveCast>(*multiplicand); expr != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) { + auto e = Add { valI->GetLeastSigOp(), expr->GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); + } else { + vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); + } + break; + } + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); + vals.push_back(expr->Generalize()); + } + continue; + } + // single expr + if (auto expr = RecursiveCast(*multiplicand); expr != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + if (valI->GetMostSigOp().Equals(*expr)) { + auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); + } else { + vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); + } + break; + } + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + vals.push_back(Exponent { *expr, Real { 1.0 } }.Generalize()); + } + continue; + } + } + + // makes all expr^1 into expr + for (auto& val : vals) { + if (auto exp = RecursiveCast>(*val); exp != nullptr) { + if (exp->GetLeastSigOp().GetValue() == 1.0) { + val = exp->GetMostSigOp().Generalize(); + } + } + if (auto mul = RecursiveCast>(*val); mul != nullptr) { + if (mul->GetMostSigOp().GetValue() == 1.0) { + val = mul->GetLeastSigOp().Generalize(); + } + } + } + + return BuildFromVector(vals); + + // return simplifiedMultiply.Copy(); +} + auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; @@ -49,24 +488,29 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons // TODO: Implement integration by parts // Detect whether integration by parts is appropriate // May need to simplify before and/or after - else if (auto mult = RecursiveCast>(*simplifiedMult); mult != nullptr) { + else { std::unique_ptr u; std::unique_ptr dv; + // Simplify the most and least significant operand + auto simplifiedMostSigOp = this->GetMostSigOp().Copy()->Accept(simplifyVisitor).value(); + auto simplifiedLeastSigOp = this->GetLeastSigOp().Copy()->Accept(simplifyVisitor).value(); + // Check the rules of LIPET // May need to recursively cast u and dv before assigning them // so that method overrides can be realized. - // TODO: change dv to v and muliply by integration variable to correctly attain dv + // TODO: change dv to v and multiply by integration variable to correctly attain dv + // ^ Not needed because the integration variable is carried over from the integration methods // Logarithm - if (mult->GetMostSigOp().Is>()) { - u = mult->GetMostSigOp().Copy(); - dv = mult->GetLeastSigOp().Copy(); - } else if (mult->GetLeastSigOp().Is>()) { - u = mult->GetLeastSigOp().Copy(); - dv = mult->GetMostSigOp().Copy(); + if (simplifiedMostSigOp->Is>()) { + u = simplifiedMostSigOp->Copy(); + dv = simplifiedLeastSigOp->Copy(); + } else if (simplifiedLeastSigOp->Is>()) { + u = simplifiedLeastSigOp->Copy(); + dv = simplifiedMostSigOp->Copy(); } // TODO: Inverse trigonometry @@ -75,53 +519,51 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons // Polynomial // TODO: Could also be exponent, in the case of (x^2)*sinx // TODO: Ensure all polynomial cases are accounted for - if (mult->GetMostSigOp().Is()) { - u = mult->GetMostSigOp().Copy(); - dv = mult->GetLeastSigOp().Copy(); - } else if (mult->GetLeastSigOp().Is()) { - u = mult->GetLeastSigOp().Copy(); - dv = mult->GetMostSigOp().Copy(); + if (simplifiedMostSigOp->Is() || simplifiedMostSigOp->Is>()) { + u = simplifiedMostSigOp->Copy(); + dv = simplifiedLeastSigOp->Copy(); + } else if (simplifiedLeastSigOp->Is() || simplifiedMostSigOp->Is>()) { + u = simplifiedLeastSigOp->Copy(); + dv = simplifiedMostSigOp->Copy(); } // Exponential - Euler's Number - if (mult->GetMostSigOp().Is>()) { - u = mult->GetMostSigOp().Copy(); - dv = mult->GetLeastSigOp().Copy(); - } else if (mult->GetLeastSigOp().Is>()) { - u = mult->GetLeastSigOp().Copy(); - dv = mult->GetMostSigOp().Copy(); + if (simplifiedMostSigOp->Is>()) { + u = simplifiedMostSigOp->Copy(); + dv = simplifiedLeastSigOp->Copy(); + } else if (simplifiedLeastSigOp->Is>()) { + u =simplifiedLeastSigOp->Copy(); + dv = simplifiedMostSigOp->Copy(); } + + // TODO: Trigonometry // Trigonometry is not implemented yet in Oasis - // Differentiate u and integrate dv - // TODO: ensure u & dv are casted as the correct types to allow method override - std::unique_ptr du = u->Differentiate(*u); - std::unique_ptr v = dv->Integrate(*dv); - - // Apply the formula: integral(udv) = uv - integral(vdu) - // Work in Progress - Subtract, Integral> subtractor { - Multiply { - *u, - *dv }, - Integral {} - }; + // Differentiate u and integrate dv + // TODO: ensure u & dv are casted as the correct types to allow method override + std::unique_ptr du = u->Differentiate(integrationVariable); + std::unique_ptr v = dv->Integrate(integrationVariable); - // test code to base the structure of above tree off of - // Multiply subtract { - // Multiply { - // Real { 1.0 }, - // Real { 2.0 } }, - // Real { 3.0 } - // }; + Integral vdu { + Multiply {*v, *du}, + integrationVariable + }; + // TODO: Make test cases in the event that IBP needs to be done again on integrated_vdu + auto integrated_vdu = vdu.Integrate(integrationVariable); + // Apply the formula: integral(udv) = uv - integral(vdu) + Subtract, Expression> subtractor { + Multiply { *u, *v }, + *integrated_vdu + }; + return subtractor.Accept(simplifyVisitor).value().Copy(); } @@ -166,8 +608,7 @@ auto Multiply::Differentiate(const Expression& differentiationVariab Multiply>>(Add, Multiply> { Multiply { *add, *right }, Multiply { *add2, *left } }) - ->Accept(simplifyVisitor) - .value(); + ->Accept(simplifyVisitor).value(); } } } From 8de01362fb2d2f7788a754c45d94d59308759a3a Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 27 Feb 2026 13:53:32 -0500 Subject: [PATCH 03/21] Started work on test cases Test cases involving Integration By Parts with Euler's Number and first/second degree variables have been completed. --- src/Multiply.cpp | 876 ++++++++++++++++++++------------------- tests/IntegrateTests.cpp | 73 +++- 2 files changed, 511 insertions(+), 438 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index d8d53673..60edc585 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -19,444 +19,448 @@ namespace Oasis { -auto Multiply::Simplify() const -> std::unique_ptr -{ - auto simplifiedMultiplicand = mostSigOp->Simplify(); - auto simplifiedMultiplier = leastSigOp->Simplify(); - - SimplifyVisitor simplifyVisitor; - - Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier }; - if (auto onezerocase = RecursiveCast>(simplifiedMultiply); onezerocase != nullptr) { - const Real& multiplicand = onezerocase->GetMostSigOp(); - const Expression& multiplier = onezerocase->GetLeastSigOp(); - if (std::abs(multiplicand.GetValue()) <= EPSILON) { - return std::make_unique(Real { 0.0 }); - } - if (multiplicand.GetValue() == 1) { - return multiplier.Simplify(); - } - } - if (auto realCase = RecursiveCast>(simplifiedMultiply); realCase != nullptr) { - const Real& multiplicand = realCase->GetMostSigOp(); - const Real& multiplier = realCase->GetLeastSigOp(); - return std::make_unique(multiplicand.GetValue() * multiplier.GetValue()); - } - - if (auto ImgCase = RecursiveCast>(simplifiedMultiply); ImgCase != nullptr) { - return std::make_unique(-1.0); - } - - if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { - auto e = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; - auto m = e.Accept(simplifyVisitor); - if (!m) { - return Divide { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); - } - return Divide { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); - } - - if (auto exprCase = RecursiveCast>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { - return std::make_unique>(exprCase->GetMostSigOp(), Real { 2.0 }); - } - } - - if (auto rMatrixCase = RecursiveCast>(simplifiedMultiply); rMatrixCase != nullptr) { - return std::make_unique(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()); - } - - if (auto matrixCase = RecursiveCast>(simplifiedMultiply); matrixCase != nullptr) { - const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp(); - const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp(); - - if (leftTerm.GetCols() == rightTerm.GetRows()) { - return std::make_unique(leftTerm.GetMatrix() * rightTerm.GetMatrix()); - } else { - // ERROR: INVALID DIMENSION - return std::make_unique>(leftTerm, rightTerm); - } - } - - // Commented out to not cause massive problems with things that need factored expressions - // // c*(a-b) - // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { - // if (negated->GetMostSigOp().GetValue()<0){ - // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, - // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); - // } else if (negated->GetMostSigOp().GetValue()>0){ - // return Subtract{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, - // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); - // } else { - // return Real{0}.Copy(); - // } - // - // } - // - // // c*(a+b) - // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { - // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, - // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); - // } - - if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { - auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - if (!ms) { - return Divide { m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); - } - // This doesn't have to check errors because errors should have already been caught above - return Divide { *std::move(ms).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value(); - } - - if (auto multCase = RecursiveCast, Divide>>(simplifiedMultiply); multCase != nullptr) { - auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Multiply { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return this->Generalize(); - } - return Divide { *std::move(ms).value(), *std::move(ns).value() }.Accept(simplifyVisitor).value(); - } - - if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp(), e); - } - return std::make_unique>(exprCase->GetMostSigOp(), - *(std::move(s).value())); - } - } - - // x*x^n - if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp(), e); - } - return std::make_unique>(exprCase->GetMostSigOp(), *(std::move(s).value())); - } - } - - if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetLeastSigOp(), e); - } - return std::make_unique>(exprCase->GetLeastSigOp(), *(std::move(s).value())); - } - } - - // x^n*x^m - if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), e); - } - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s).value())); - } - } - - // a*x*x - if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }); - } - } - - // a*x*b*x - if (auto exprCase = RecursiveCast, Multiply>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return this->Generalize(); - } - return std::make_unique>( - *(std::move(ms).value()), - *(std::move(ns).value())); - } - } - - // a*x^n*x - if (auto exprCase = RecursiveCast>, Expression>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { - auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e }); - } - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - *(std::move(s).value()) }); - } - } - - // a*x*x^n - if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), e }); - } - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - *(std::move(s).value()) }); - } - if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), - Exponent { exprCase->GetMostSigOp().GetMostSigOp(), e }); - } - return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), - Exponent { exprCase->GetMostSigOp().GetMostSigOp(), - *(std::move(s).value()) }); - } - } - - // a*x^n*b*x - if (auto exprCase = RecursiveCast, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>( - m, Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), n }); - } - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - *(std::move(ns).value()) }); - } - } - - if (auto exprCase = RecursiveCast, Multiply, Expression>>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>(m, - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - n }); - } - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - *(std::move(ns).value()) }); - } - } - - if (auto exprCase = RecursiveCast>, Multiply>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>(m, - Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), - n }); - } - - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), - *(std::move(ns).value()) }); - } - } - - // a*x^n*x^m - if (auto exprCase = RecursiveCast>, Exponent>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>( - exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - e }); - } - return std::make_unique>( - exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - *(std::move(s).value()) }); - } - } - - // a*x^n*b*x^m - if (auto exprCase = RecursiveCast>, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>( - m, - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - n }); - } - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - *(std::move(ns).value()) }); - } - } - - // if (auto negate = RecursiveCast>>>(simplifiedMultiply); negate != nullptr){ - // return Add{Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetMostSigOp()}, - // Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetLeastSigOp()}}.Simplify(); - // } - - // multiply add like terms - std::vector> multiplies; - std::vector> vals; - simplifiedMultiply.Flatten(multiplies); - for (const auto& multiplicand : multiplies) { - size_t i = 0; - if (auto real = RecursiveCast(*multiplicand); real != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast(*vals[i]); valI != nullptr) { - vals[i] = Real { valI->GetValue() * real->GetValue() }.Generalize(); - break; - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - vals.push_back(multiplicand->Generalize()); - } - continue; - } - // single i - if (auto img = RecursiveCast(*multiplicand); img != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { Imaginary {}, e }.Generalize(); - } else { - vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); - } - - break; - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - vals.push_back(Exponent { Imaginary {}, Real { 1.0 } }.Generalize()); - } - continue; - } - // i^n - if (auto img = RecursiveCast>(*multiplicand); img != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - auto e = Add { valI->GetLeastSigOp(), img->GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { Imaginary {}, e }.Generalize(); - } else { - vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); - } - break; - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); - vals.push_back(img->Generalize()); - } - continue; - } - // expr^n - if (auto expr = RecursiveCast>(*multiplicand); expr != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) { - auto e = Add { valI->GetLeastSigOp(), expr->GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); - } else { - vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); - } - break; - } - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); - vals.push_back(expr->Generalize()); - } - continue; - } - // single expr - if (auto expr = RecursiveCast(*multiplicand); expr != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - if (valI->GetMostSigOp().Equals(*expr)) { - auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); - } else { - vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); - } - break; - } - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - vals.push_back(Exponent { *expr, Real { 1.0 } }.Generalize()); - } - continue; - } - } - - // makes all expr^1 into expr - for (auto& val : vals) { - if (auto exp = RecursiveCast>(*val); exp != nullptr) { - if (exp->GetLeastSigOp().GetValue() == 1.0) { - val = exp->GetMostSigOp().Generalize(); - } - } - if (auto mul = RecursiveCast>(*val); mul != nullptr) { - if (mul->GetMostSigOp().GetValue() == 1.0) { - val = mul->GetLeastSigOp().Generalize(); - } - } - } +// Mark for deletion: +// auto Multiply::Simplify() const -> std::unique_ptr +// { +// auto simplifiedMultiplicand = mostSigOp->Simplify(); +// auto simplifiedMultiplier = leastSigOp->Simplify(); +// +// SimplifyVisitor simplifyVisitor; +// +// Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier }; +// if (auto onezerocase = RecursiveCast>(simplifiedMultiply); onezerocase != nullptr) { +// const Real& multiplicand = onezerocase->GetMostSigOp(); +// const Expression& multiplier = onezerocase->GetLeastSigOp(); +// if (std::abs(multiplicand.GetValue()) <= EPSILON) { +// return std::make_unique(Real { 0.0 }); +// } +// if (multiplicand.GetValue() == 1) { +// return multiplier.Simplify(); +// } +// } +// if (auto realCase = RecursiveCast>(simplifiedMultiply); realCase != nullptr) { +// const Real& multiplicand = realCase->GetMostSigOp(); +// const Real& multiplier = realCase->GetLeastSigOp(); +// return std::make_unique(multiplicand.GetValue() * multiplier.GetValue()); +// } +// +// if (auto ImgCase = RecursiveCast>(simplifiedMultiply); ImgCase != nullptr) { +// return std::make_unique(-1.0); +// } +// +// if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { +// auto e = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; +// auto m = e.Accept(simplifyVisitor); +// if (!m) { +// return Divide { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); +// } +// return Divide { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); +// } +// +// if (auto exprCase = RecursiveCast>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { +// return std::make_unique>(exprCase->GetMostSigOp(), Real { 2.0 }); +// } +// } +// +// if (auto rMatrixCase = RecursiveCast>(simplifiedMultiply); rMatrixCase != nullptr) { +// return std::make_unique(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()); +// } +// +// if (auto matrixCase = RecursiveCast>(simplifiedMultiply); matrixCase != nullptr) { +// const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp(); +// const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp(); +// +// if (leftTerm.GetCols() == rightTerm.GetRows()) { +// return std::make_unique(leftTerm.GetMatrix() * rightTerm.GetMatrix()); +// } else { +// // ERROR: INVALID DIMENSION +// return std::make_unique>(leftTerm, rightTerm); +// } +// } +// +// // Commented out to not cause massive problems with things that need factored expressions +// // // c*(a-b) +// // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { +// // if (negated->GetMostSigOp().GetValue()<0){ +// // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, +// // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); +// // } else if (negated->GetMostSigOp().GetValue()>0){ +// // return Subtract{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, +// // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); +// // } else { +// // return Real{0}.Copy(); +// // } +// // +// // } +// // +// // // c*(a+b) +// // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { +// // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, +// // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); +// // } +// +// if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { +// auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// if (!ms) { +// return Divide { m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); +// } +// // This doesn't have to check errors because errors should have already been caught above +// return Divide { *std::move(ms).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value(); +// } +// +// if (auto multCase = RecursiveCast, Divide>>(simplifiedMultiply); multCase != nullptr) { +// auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// auto n = Multiply { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }; +// auto ns = n.Accept(simplifyVisitor); +// if (!ms || !ns) { +// return this->Generalize(); +// } +// return Divide { *std::move(ms).value(), *std::move(ns).value() }.Accept(simplifyVisitor).value(); +// } +// +// if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { +// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetMostSigOp(), e); +// } +// return std::make_unique>(exprCase->GetMostSigOp(), +// *(std::move(s).value())); +// } +// } +// +// // x*x^n +// if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { +// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetMostSigOp(), e); +// } +// return std::make_unique>(exprCase->GetMostSigOp(), *(std::move(s).value())); +// } +// } +// +// if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) { +// auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetLeastSigOp(), e); +// } +// return std::make_unique>(exprCase->GetLeastSigOp(), *(std::move(s).value())); +// } +// } +// +// // x^n*x^m +// if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { +// auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), e); +// } +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s).value())); +// } +// } +// +// // a*x*x +// if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) { +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }); +// } +// } +// +// // a*x*b*x +// if (auto exprCase = RecursiveCast, Multiply>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { +// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// auto n = Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }; +// auto ns = n.Accept(simplifyVisitor); +// if (!ms || !ns) { +// return this->Generalize(); +// } +// return std::make_unique>( +// *(std::move(ms).value()), +// *(std::move(ns).value())); +// } +// } +// +// // a*x^n*x +// if (auto exprCase = RecursiveCast>, Expression>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { +// auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e }); +// } +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), +// *(std::move(s).value()) }); +// } +// } +// +// // a*x*x^n +// if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { +// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), e }); +// } +// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), +// *(std::move(s).value()) }); +// } +// if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { +// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), +// Exponent { exprCase->GetMostSigOp().GetMostSigOp(), e }); +// } +// return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), +// Exponent { exprCase->GetMostSigOp().GetMostSigOp(), +// *(std::move(s).value()) }); +// } +// } +// +// // a*x^n*b*x +// if (auto exprCase = RecursiveCast, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { +// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto ns = n.Accept(simplifyVisitor); +// if (!ms || !ns) { +// return std::make_unique>( +// m, Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), n }); +// } +// return std::make_unique>( +// *(std::move(ms).value()), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), +// *(std::move(ns).value()) }); +// } +// } +// +// if (auto exprCase = RecursiveCast, Multiply, Expression>>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) { +// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// auto n = Add { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto ns = n.Accept(simplifyVisitor); +// if (!ms || !ns) { +// return std::make_unique>(m, +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), +// n }); +// } +// return std::make_unique>( +// *(std::move(ms).value()), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), +// *(std::move(ns).value()) }); +// } +// } +// +// if (auto exprCase = RecursiveCast>, Multiply>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { +// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// auto n = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; +// auto ns = n.Accept(simplifyVisitor); +// if (!ms || !ns) { +// return std::make_unique>(m, +// Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), +// n }); +// } +// +// return std::make_unique>( +// *(std::move(ms).value()), +// Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), +// *(std::move(ns).value()) }); +// } +// } +// +// // a*x^n*x^m +// if (auto exprCase = RecursiveCast>, Exponent>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { +// auto e = Add { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// return std::make_unique>( +// exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), +// e }); +// } +// return std::make_unique>( +// exprCase->GetMostSigOp().GetMostSigOp(), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), +// *(std::move(s).value()) }); +// } +// } +// +// // a*x^n*b*x^m +// if (auto exprCase = RecursiveCast>, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { +// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { +// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; +// auto ms = m.Accept(simplifyVisitor); +// auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; +// auto ns = n.Accept(simplifyVisitor); +// if (!ms || !ns) { +// return std::make_unique>( +// m, +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), +// n }); +// } +// return std::make_unique>( +// *(std::move(ms).value()), +// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), +// *(std::move(ns).value()) }); +// } +// } +// +// // if (auto negate = RecursiveCast>>>(simplifiedMultiply); negate != nullptr){ +// // return Add{Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetMostSigOp()}, +// // Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetLeastSigOp()}}.Simplify(); +// // } +// +// // multiply add like terms +// std::vector> multiplies; +// std::vector> vals; +// simplifiedMultiply.Flatten(multiplies); +// for (const auto& multiplicand : multiplies) { +// size_t i = 0; +// if (auto real = RecursiveCast(*multiplicand); real != nullptr) { +// for (; i < vals.size(); i++) { +// if (auto valI = RecursiveCast(*vals[i]); valI != nullptr) { +// vals[i] = Real { valI->GetValue() * real->GetValue() }.Generalize(); +// break; +// } +// } +// if (i >= vals.size()) { +// // check to make sure it is one thing only +// vals.push_back(multiplicand->Generalize()); +// } +// continue; +// } +// // single i +// if (auto img = RecursiveCast(*multiplicand); img != nullptr) { +// for (; i < vals.size(); i++) { +// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { +// auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// vals[i] = Exponent { Imaginary {}, e }.Generalize(); +// } else { +// vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); +// } +// +// break; +// } +// } +// if (i >= vals.size()) { +// // check to make sure it is one thing only +// vals.push_back(Exponent { Imaginary {}, Real { 1.0 } }.Generalize()); +// } +// continue; +// } +// // i^n +// if (auto img = RecursiveCast>(*multiplicand); img != nullptr) { +// for (; i < vals.size(); i++) { +// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { +// auto e = Add { valI->GetLeastSigOp(), img->GetLeastSigOp() }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// vals[i] = Exponent { Imaginary {}, e }.Generalize(); +// } else { +// vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); +// } +// break; +// } +// } +// if (i >= vals.size()) { +// // check to make sure it is one thing only +// // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); +// vals.push_back(img->Generalize()); +// } +// continue; +// } +// // expr^n +// if (auto expr = RecursiveCast>(*multiplicand); expr != nullptr) { +// for (; i < vals.size(); i++) { +// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { +// if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) { +// auto e = Add { valI->GetLeastSigOp(), expr->GetLeastSigOp() }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); +// } else { +// vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); +// } +// break; +// } +// } +// } +// if (i >= vals.size()) { +// // check to make sure it is one thing only +// // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); +// vals.push_back(expr->Generalize()); +// } +// continue; +// } +// // single expr +// if (auto expr = RecursiveCast(*multiplicand); expr != nullptr) { +// for (; i < vals.size(); i++) { +// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { +// if (valI->GetMostSigOp().Equals(*expr)) { +// auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; +// auto s = e.Accept(simplifyVisitor); +// if (!s) { +// vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); +// } else { +// vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); +// } +// break; +// } +// } +// } +// if (i >= vals.size()) { +// // check to make sure it is one thing only +// vals.push_back(Exponent { *expr, Real { 1.0 } }.Generalize()); +// } +// continue; +// } +// } +// +// // makes all expr^1 into expr +// for (auto& val : vals) { +// if (auto exp = RecursiveCast>(*val); exp != nullptr) { +// if (exp->GetLeastSigOp().GetValue() == 1.0) { +// val = exp->GetMostSigOp().Generalize(); +// } +// } +// if (auto mul = RecursiveCast>(*val); mul != nullptr) { +// if (mul->GetMostSigOp().GetValue() == 1.0) { +// val = mul->GetLeastSigOp().Generalize(); +// } +// } +// } +// +// return BuildFromVector(vals); +// +// // return simplifiedMultiply.Copy(); +// } - return BuildFromVector(vals); +auto Recast_DV(Expression& expression) { - // return simplifiedMultiply.Copy(); } - auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; @@ -563,7 +567,7 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons *integrated_vdu }; - return subtractor.Accept(simplifyVisitor).value().Copy(); + return subtractor.Accept(simplifyVisitor).value(); } diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index 7e1c2ede..0a80c7b6 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -2,12 +2,13 @@ #include "Oasis/Add.hpp" #include "Oasis/Divide.hpp" +#include "Oasis/EulerNumber.hpp" #include "Oasis/Exponent.hpp" #include "Oasis/Multiply.hpp" #include "Oasis/Real.hpp" +#include "Oasis/SimplifyVisitor.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" -#include "Oasis/SimplifyVisitor.hpp" inline Oasis::SimplifyVisitor simplifyVisitor{}; @@ -181,4 +182,72 @@ TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]") auto simplified = integrated->Accept(simplifyVisitor).value(); REQUIRE(simplified->Equals(*(integral.Accept(simplifyVisitor).value()))); -} \ No newline at end of file +} + + +TEST_CASE("Integration By Parts: Variable and Euler's Number", "[Integrate][Variable][Euler]") +{ + Oasis::Variable var { "x" }; + + Oasis::Multiply> integrand { + Oasis::Variable { var.GetName() }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add>, Oasis::Exponent>, 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> integrand { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + + Oasis::Add, Oasis::Exponent>, + Oasis::Multiply>>, Oasis::Variable> integral { + Oasis::Add { + Oasis::Subtract { + Oasis::Multiply { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } }, + Oasis::Multiply { + Oasis::Real { 2 }, + 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))); +} + + + + + From cae64fb00291d8cb19cd136439671c2910c3fc74 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 27 Feb 2026 20:51:36 -0500 Subject: [PATCH 04/21] Added Exponent-Logarithm Integration Test Added a new test for integrating the multiplication (using IBP) between a variable raised to an real number and the natural logarithm. --- tests/IntegrateTests.cpp | 41 ++++++++++++++++++++++++++++++++++++++-- 1 file changed, 39 insertions(+), 2 deletions(-) diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index 0a80c7b6..f68b5728 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -9,6 +9,7 @@ #include "Oasis/SimplifyVisitor.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" +#include "Oasis/Log.hpp" inline Oasis::SimplifyVisitor simplifyVisitor{}; @@ -184,7 +185,6 @@ TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]") REQUIRE(simplified->Equals(*(integral.Accept(simplifyVisitor).value()))); } - TEST_CASE("Integration By Parts: Variable and Euler's Number", "[Integrate][Variable][Euler]") { Oasis::Variable var { "x" }; @@ -234,7 +234,7 @@ TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Expo Oasis::Multiply { Oasis::Real { 2 }, Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } }, - }, + }, Oasis::Variable { "C" } } }; @@ -247,7 +247,44 @@ TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Expo REQUIRE((integrated->Equals(*ptr))); } +TEST_CASE("Integration By Parts: Exponent and Logarithm", "[Integrate][Exponent][Logarithm]") +{ + Oasis::Variable var { "x" }; + Oasis::Multiply, Oasis::Log> integrand { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add, + Oasis::Exponent>, Oasis::Log>, + Oasis::Multiply, Oasis::Exponent>>, + Oasis::Variable> integral { + Oasis::Add { + Oasis::Subtract { + Oasis::Multiply { + Oasis::Multiply { + Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 3 } }, + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } + }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }, + Oasis::Multiply { + Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 9 } }, + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } + } + }, + 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))); +} From ffb6cd366a8ac25b7b481f4c545c54dd79db768f Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 13 Mar 2026 18:46:00 -0400 Subject: [PATCH 05/21] Added integration of Exponent with EulerNumber base and Variable power Integration of an instance of the Exponent class containing an EulerNumber base raised to a degree-one Variable has been implemented. --- src/Exponent.cpp | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/src/Exponent.cpp b/src/Exponent.cpp index 08e3bd54..65f040c9 100644 --- a/src/Exponent.cpp +++ b/src/Exponent.cpp @@ -51,6 +51,18 @@ auto Exponent::Integrate(const Expression& integrationVariable) cons return adder.Accept(simplifyVisitor).value(); } } + + // EulerNumber raised to a (degree-one) Variable + if (auto eulerBase = RecursiveCast>(*simplifiedExponent); eulerBase != nullptr) { + const Variable& expPow = eulerBase->GetLeastSigOp(); + + Add adder { + Exponent { EulerNumber{}, expPow }, + Variable { "C" } + }; + + return adder.Accept(simplifyVisitor).value(); + } } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; From a0fa0f91ee42b7002ca840cb61808ea64adc08b2 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 13 Mar 2026 18:58:57 -0400 Subject: [PATCH 06/21] Revised implementation of IBP Formula After running test cases, the code was rewritten to rely on the Multiply instance's most and least significant operands, thereby preserving method overrides needed for integration. --- src/Multiply.cpp | 985 +++++++++++++++++++++++------------------------ 1 file changed, 482 insertions(+), 503 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 60edc585..5c5e7847 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -20,447 +20,445 @@ namespace Oasis { // Mark for deletion: -// auto Multiply::Simplify() const -> std::unique_ptr -// { -// auto simplifiedMultiplicand = mostSigOp->Simplify(); -// auto simplifiedMultiplier = leastSigOp->Simplify(); -// -// SimplifyVisitor simplifyVisitor; -// -// Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier }; -// if (auto onezerocase = RecursiveCast>(simplifiedMultiply); onezerocase != nullptr) { -// const Real& multiplicand = onezerocase->GetMostSigOp(); -// const Expression& multiplier = onezerocase->GetLeastSigOp(); -// if (std::abs(multiplicand.GetValue()) <= EPSILON) { -// return std::make_unique(Real { 0.0 }); -// } -// if (multiplicand.GetValue() == 1) { -// return multiplier.Simplify(); -// } -// } -// if (auto realCase = RecursiveCast>(simplifiedMultiply); realCase != nullptr) { -// const Real& multiplicand = realCase->GetMostSigOp(); -// const Real& multiplier = realCase->GetLeastSigOp(); -// return std::make_unique(multiplicand.GetValue() * multiplier.GetValue()); -// } -// -// if (auto ImgCase = RecursiveCast>(simplifiedMultiply); ImgCase != nullptr) { -// return std::make_unique(-1.0); -// } -// -// if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { -// auto e = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; -// auto m = e.Accept(simplifyVisitor); -// if (!m) { -// return Divide { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); -// } -// return Divide { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); -// } -// -// if (auto exprCase = RecursiveCast>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { -// return std::make_unique>(exprCase->GetMostSigOp(), Real { 2.0 }); -// } -// } -// -// if (auto rMatrixCase = RecursiveCast>(simplifiedMultiply); rMatrixCase != nullptr) { -// return std::make_unique(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()); -// } -// -// if (auto matrixCase = RecursiveCast>(simplifiedMultiply); matrixCase != nullptr) { -// const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp(); -// const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp(); -// -// if (leftTerm.GetCols() == rightTerm.GetRows()) { -// return std::make_unique(leftTerm.GetMatrix() * rightTerm.GetMatrix()); -// } else { -// // ERROR: INVALID DIMENSION -// return std::make_unique>(leftTerm, rightTerm); -// } -// } -// -// // Commented out to not cause massive problems with things that need factored expressions -// // // c*(a-b) -// // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { -// // if (negated->GetMostSigOp().GetValue()<0){ -// // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, -// // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); -// // } else if (negated->GetMostSigOp().GetValue()>0){ -// // return Subtract{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, -// // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); -// // } else { -// // return Real{0}.Copy(); -// // } -// // -// // } -// // -// // // c*(a+b) -// // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { -// // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, -// // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); -// // } -// -// if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { -// auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// if (!ms) { -// return Divide { m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); -// } -// // This doesn't have to check errors because errors should have already been caught above -// return Divide { *std::move(ms).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value(); -// } -// -// if (auto multCase = RecursiveCast, Divide>>(simplifiedMultiply); multCase != nullptr) { -// auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// auto n = Multiply { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }; -// auto ns = n.Accept(simplifyVisitor); -// if (!ms || !ns) { -// return this->Generalize(); -// } -// return Divide { *std::move(ms).value(), *std::move(ns).value() }.Accept(simplifyVisitor).value(); -// } -// -// if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { -// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetMostSigOp(), e); -// } -// return std::make_unique>(exprCase->GetMostSigOp(), -// *(std::move(s).value())); -// } -// } -// -// // x*x^n -// if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { -// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetMostSigOp(), e); -// } -// return std::make_unique>(exprCase->GetMostSigOp(), *(std::move(s).value())); -// } -// } -// -// if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) { -// auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetLeastSigOp(), e); -// } -// return std::make_unique>(exprCase->GetLeastSigOp(), *(std::move(s).value())); -// } -// } -// -// // x^n*x^m -// if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { -// auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), e); -// } -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s).value())); -// } -// } -// -// // a*x*x -// if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) { -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }); -// } -// } -// -// // a*x*b*x -// if (auto exprCase = RecursiveCast, Multiply>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { -// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// auto n = Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }; -// auto ns = n.Accept(simplifyVisitor); -// if (!ms || !ns) { -// return this->Generalize(); -// } -// return std::make_unique>( -// *(std::move(ms).value()), -// *(std::move(ns).value())); -// } -// } -// -// // a*x^n*x -// if (auto exprCase = RecursiveCast>, Expression>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { -// auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e }); -// } -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), -// *(std::move(s).value()) }); -// } -// } -// -// // a*x*x^n -// if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { -// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), e }); -// } -// return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), -// *(std::move(s).value()) }); -// } -// if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { -// auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), -// Exponent { exprCase->GetMostSigOp().GetMostSigOp(), e }); -// } -// return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), -// Exponent { exprCase->GetMostSigOp().GetMostSigOp(), -// *(std::move(s).value()) }); -// } -// } -// -// // a*x^n*b*x -// if (auto exprCase = RecursiveCast, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { -// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto ns = n.Accept(simplifyVisitor); -// if (!ms || !ns) { -// return std::make_unique>( -// m, Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), n }); -// } -// return std::make_unique>( -// *(std::move(ms).value()), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), -// *(std::move(ns).value()) }); -// } -// } -// -// if (auto exprCase = RecursiveCast, Multiply, Expression>>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) { -// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// auto n = Add { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto ns = n.Accept(simplifyVisitor); -// if (!ms || !ns) { -// return std::make_unique>(m, -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), -// n }); -// } -// return std::make_unique>( -// *(std::move(ms).value()), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), -// *(std::move(ns).value()) }); -// } -// } -// -// if (auto exprCase = RecursiveCast>, Multiply>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { -// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// auto n = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; -// auto ns = n.Accept(simplifyVisitor); -// if (!ms || !ns) { -// return std::make_unique>(m, -// Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), -// n }); -// } -// -// return std::make_unique>( -// *(std::move(ms).value()), -// Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), -// *(std::move(ns).value()) }); -// } -// } -// -// // a*x^n*x^m -// if (auto exprCase = RecursiveCast>, Exponent>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { -// auto e = Add { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// return std::make_unique>( -// exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), -// e }); -// } -// return std::make_unique>( -// exprCase->GetMostSigOp().GetMostSigOp(), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), -// *(std::move(s).value()) }); -// } -// } -// -// // a*x^n*b*x^m -// if (auto exprCase = RecursiveCast>, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { -// if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { -// auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; -// auto ms = m.Accept(simplifyVisitor); -// auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; -// auto ns = n.Accept(simplifyVisitor); -// if (!ms || !ns) { -// return std::make_unique>( -// m, -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), -// n }); -// } -// return std::make_unique>( -// *(std::move(ms).value()), -// Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), -// *(std::move(ns).value()) }); -// } -// } -// -// // if (auto negate = RecursiveCast>>>(simplifiedMultiply); negate != nullptr){ -// // return Add{Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetMostSigOp()}, -// // Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetLeastSigOp()}}.Simplify(); -// // } -// -// // multiply add like terms -// std::vector> multiplies; -// std::vector> vals; -// simplifiedMultiply.Flatten(multiplies); -// for (const auto& multiplicand : multiplies) { -// size_t i = 0; -// if (auto real = RecursiveCast(*multiplicand); real != nullptr) { -// for (; i < vals.size(); i++) { -// if (auto valI = RecursiveCast(*vals[i]); valI != nullptr) { -// vals[i] = Real { valI->GetValue() * real->GetValue() }.Generalize(); -// break; -// } -// } -// if (i >= vals.size()) { -// // check to make sure it is one thing only -// vals.push_back(multiplicand->Generalize()); -// } -// continue; -// } -// // single i -// if (auto img = RecursiveCast(*multiplicand); img != nullptr) { -// for (; i < vals.size(); i++) { -// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { -// auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// vals[i] = Exponent { Imaginary {}, e }.Generalize(); -// } else { -// vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); -// } -// -// break; -// } -// } -// if (i >= vals.size()) { -// // check to make sure it is one thing only -// vals.push_back(Exponent { Imaginary {}, Real { 1.0 } }.Generalize()); -// } -// continue; -// } -// // i^n -// if (auto img = RecursiveCast>(*multiplicand); img != nullptr) { -// for (; i < vals.size(); i++) { -// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { -// auto e = Add { valI->GetLeastSigOp(), img->GetLeastSigOp() }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// vals[i] = Exponent { Imaginary {}, e }.Generalize(); -// } else { -// vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); -// } -// break; -// } -// } -// if (i >= vals.size()) { -// // check to make sure it is one thing only -// // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); -// vals.push_back(img->Generalize()); -// } -// continue; -// } -// // expr^n -// if (auto expr = RecursiveCast>(*multiplicand); expr != nullptr) { -// for (; i < vals.size(); i++) { -// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { -// if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) { -// auto e = Add { valI->GetLeastSigOp(), expr->GetLeastSigOp() }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); -// } else { -// vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); -// } -// break; -// } -// } -// } -// if (i >= vals.size()) { -// // check to make sure it is one thing only -// // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); -// vals.push_back(expr->Generalize()); -// } -// continue; -// } -// // single expr -// if (auto expr = RecursiveCast(*multiplicand); expr != nullptr) { -// for (; i < vals.size(); i++) { -// if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { -// if (valI->GetMostSigOp().Equals(*expr)) { -// auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; -// auto s = e.Accept(simplifyVisitor); -// if (!s) { -// vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); -// } else { -// vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); -// } -// break; -// } -// } -// } -// if (i >= vals.size()) { -// // check to make sure it is one thing only -// vals.push_back(Exponent { *expr, Real { 1.0 } }.Generalize()); -// } -// continue; -// } -// } -// -// // makes all expr^1 into expr -// for (auto& val : vals) { -// if (auto exp = RecursiveCast>(*val); exp != nullptr) { -// if (exp->GetLeastSigOp().GetValue() == 1.0) { -// val = exp->GetMostSigOp().Generalize(); -// } -// } -// if (auto mul = RecursiveCast>(*val); mul != nullptr) { -// if (mul->GetMostSigOp().GetValue() == 1.0) { -// val = mul->GetLeastSigOp().Generalize(); -// } -// } -// } -// -// return BuildFromVector(vals); -// -// // return simplifiedMultiply.Copy(); -// } +// Keep in until merge to resolve other methods dependent on this one +auto Multiply::Simplify() const -> std::unique_ptr +{ + auto simplifiedMultiplicand = mostSigOp->Simplify(); + auto simplifiedMultiplier = leastSigOp->Simplify(); + + SimplifyVisitor simplifyVisitor; + + Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier }; + if (auto onezerocase = RecursiveCast>(simplifiedMultiply); onezerocase != nullptr) { + const Real& multiplicand = onezerocase->GetMostSigOp(); + const Expression& multiplier = onezerocase->GetLeastSigOp(); + if (std::abs(multiplicand.GetValue()) <= EPSILON) { + return std::make_unique(Real { 0.0 }); + } + if (multiplicand.GetValue() == 1) { + return multiplier.Simplify(); + } + } + if (auto realCase = RecursiveCast>(simplifiedMultiply); realCase != nullptr) { + const Real& multiplicand = realCase->GetMostSigOp(); + const Real& multiplier = realCase->GetLeastSigOp(); + return std::make_unique(multiplicand.GetValue() * multiplier.GetValue()); + } -auto Recast_DV(Expression& expression) { + if (auto ImgCase = RecursiveCast>(simplifiedMultiply); ImgCase != nullptr) { + return std::make_unique(-1.0); + } + if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { + auto e = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; + auto m = e.Accept(simplifyVisitor); + if (!m) { + return Divide { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); + } + return Divide { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); + } + + if (auto exprCase = RecursiveCast>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { + return std::make_unique>(exprCase->GetMostSigOp(), Real { 2.0 }); + } + } + + if (auto rMatrixCase = RecursiveCast>(simplifiedMultiply); rMatrixCase != nullptr) { + return std::make_unique(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()); + } + + if (auto matrixCase = RecursiveCast>(simplifiedMultiply); matrixCase != nullptr) { + const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp(); + const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp(); + + if (leftTerm.GetCols() == rightTerm.GetRows()) { + return std::make_unique(leftTerm.GetMatrix() * rightTerm.GetMatrix()); + } else { + // ERROR: INVALID DIMENSION + return std::make_unique>(leftTerm, rightTerm); + } + } + + // Commented out to not cause massive problems with things that need factored expressions + // // c*(a-b) + // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { + // if (negated->GetMostSigOp().GetValue()<0){ + // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, + // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); + // } else if (negated->GetMostSigOp().GetValue()>0){ + // return Subtract{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, + // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); + // } else { + // return Real{0}.Copy(); + // } + // + // } + // + // // c*(a+b) + // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { + // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, + // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); + // } + + if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { + auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + if (!ms) { + return Divide { m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); + } + // This doesn't have to check errors because errors should have already been caught above + return Divide { *std::move(ms).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value(); + } + + if (auto multCase = RecursiveCast, Divide>>(simplifiedMultiply); multCase != nullptr) { + auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Multiply { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return this->Generalize(); + } + return Divide { *std::move(ms).value(), *std::move(ns).value() }.Accept(simplifyVisitor).value(); + } + + if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp(), e); + } + return std::make_unique>(exprCase->GetMostSigOp(), + *(std::move(s).value())); + } + } + + // x*x^n + if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp(), e); + } + return std::make_unique>(exprCase->GetMostSigOp(), *(std::move(s).value())); + } + } + + if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetLeastSigOp(), e); + } + return std::make_unique>(exprCase->GetLeastSigOp(), *(std::move(s).value())); + } + } + + // x^n*x^m + if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), e); + } + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s).value())); + } + } + + // a*x*x + if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }); + } + } + + // a*x*b*x + if (auto exprCase = RecursiveCast, Multiply>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return this->Generalize(); + } + return std::make_unique>( + *(std::move(ms).value()), + *(std::move(ns).value())); + } + } + + // a*x^n*x + if (auto exprCase = RecursiveCast>, Expression>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { + auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e }); + } + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + *(std::move(s).value()) }); + } + } + + // a*x*x^n + if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), e }); + } + return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + *(std::move(s).value()) }); + } + if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), + Exponent { exprCase->GetMostSigOp().GetMostSigOp(), e }); + } + return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), + Exponent { exprCase->GetMostSigOp().GetMostSigOp(), + *(std::move(s).value()) }); + } + } + + // a*x^n*b*x + if (auto exprCase = RecursiveCast, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>( + m, Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), n }); + } + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + *(std::move(ns).value()) }); + } + } + + if (auto exprCase = RecursiveCast, Multiply, Expression>>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>(m, + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + n }); + } + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), + *(std::move(ns).value()) }); + } + } + + if (auto exprCase = RecursiveCast>, Multiply>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>(m, + Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), + n }); + } + + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), + *(std::move(ns).value()) }); + } + } + + // a*x^n*x^m + if (auto exprCase = RecursiveCast>, Exponent>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { + auto e = Add { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + return std::make_unique>( + exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + e }); + } + return std::make_unique>( + exprCase->GetMostSigOp().GetMostSigOp(), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + *(std::move(s).value()) }); + } + } + + // a*x^n*b*x^m + if (auto exprCase = RecursiveCast>, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { + if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { + auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; + auto ms = m.Accept(simplifyVisitor); + auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; + auto ns = n.Accept(simplifyVisitor); + if (!ms || !ns) { + return std::make_unique>( + m, + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + n }); + } + return std::make_unique>( + *(std::move(ms).value()), + Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), + *(std::move(ns).value()) }); + } + } + + // if (auto negate = RecursiveCast>>>(simplifiedMultiply); negate != nullptr){ + // return Add{Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetMostSigOp()}, + // Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetLeastSigOp()}}.Simplify(); + // } + + // multiply add like terms + std::vector> multiplies; + std::vector> vals; + simplifiedMultiply.Flatten(multiplies); + for (const auto& multiplicand : multiplies) { + size_t i = 0; + if (auto real = RecursiveCast(*multiplicand); real != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast(*vals[i]); valI != nullptr) { + vals[i] = Real { valI->GetValue() * real->GetValue() }.Generalize(); + break; + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + vals.push_back(multiplicand->Generalize()); + } + continue; + } + // single i + if (auto img = RecursiveCast(*multiplicand); img != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { Imaginary {}, e }.Generalize(); + } else { + vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); + } + + break; + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + vals.push_back(Exponent { Imaginary {}, Real { 1.0 } }.Generalize()); + } + continue; + } + // i^n + if (auto img = RecursiveCast>(*multiplicand); img != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + auto e = Add { valI->GetLeastSigOp(), img->GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { Imaginary {}, e }.Generalize(); + } else { + vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); + } + break; + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); + vals.push_back(img->Generalize()); + } + continue; + } + // expr^n + if (auto expr = RecursiveCast>(*multiplicand); expr != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) { + auto e = Add { valI->GetLeastSigOp(), expr->GetLeastSigOp() }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); + } else { + vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); + } + break; + } + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); + vals.push_back(expr->Generalize()); + } + continue; + } + // single expr + if (auto expr = RecursiveCast(*multiplicand); expr != nullptr) { + for (; i < vals.size(); i++) { + if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { + if (valI->GetMostSigOp().Equals(*expr)) { + auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; + auto s = e.Accept(simplifyVisitor); + if (!s) { + vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); + } else { + vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); + } + break; + } + } + } + if (i >= vals.size()) { + // check to make sure it is one thing only + vals.push_back(Exponent { *expr, Real { 1.0 } }.Generalize()); + } + continue; + } + } + + // makes all expr^1 into expr + for (auto& val : vals) { + if (auto exp = RecursiveCast>(*val); exp != nullptr) { + if (exp->GetLeastSigOp().GetValue() == 1.0) { + val = exp->GetMostSigOp().Generalize(); + } + } + if (auto mul = RecursiveCast>(*val); mul != nullptr) { + if (mul->GetMostSigOp().GetValue() == 1.0) { + val = mul->GetLeastSigOp().Generalize(); + } + } + } + + return BuildFromVector(vals); + + return simplifiedMultiply.Copy(); } + auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; @@ -489,89 +487,70 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons } } - // TODO: Implement integration by parts + // TODO: Finish implementation of integration by parts // Detect whether integration by parts is appropriate // May need to simplify before and/or after else { - std::unique_ptr u; - std::unique_ptr dv; - - // Simplify the most and least significant operand - auto simplifiedMostSigOp = this->GetMostSigOp().Copy()->Accept(simplifyVisitor).value(); - auto simplifiedLeastSigOp = this->GetLeastSigOp().Copy()->Accept(simplifyVisitor).value(); - - // Check the rules of LIPET - - // May need to recursively cast u and dv before assigning them - // so that method overrides can be realized. - - // TODO: change dv to v and multiply by integration variable to correctly attain dv - // ^ Not needed because the integration variable is carried over from the integration methods - - // Logarithm - if (simplifiedMostSigOp->Is>()) { - u = simplifiedMostSigOp->Copy(); - dv = simplifiedLeastSigOp->Copy(); - } else if (simplifiedLeastSigOp->Is>()) { - u = simplifiedLeastSigOp->Copy(); - dv = simplifiedMostSigOp->Copy(); - } - - // TODO: Inverse trigonometry - // Inverse trigonometry is not implemented yet in Oasis - - // Polynomial - // TODO: Could also be exponent, in the case of (x^2)*sinx - // TODO: Ensure all polynomial cases are accounted for - if (simplifiedMostSigOp->Is() || simplifiedMostSigOp->Is>()) { - u = simplifiedMostSigOp->Copy(); - dv = simplifiedLeastSigOp->Copy(); - } else if (simplifiedLeastSigOp->Is() || simplifiedMostSigOp->Is>()) { - u = simplifiedLeastSigOp->Copy(); - dv = simplifiedMostSigOp->Copy(); - } - - // Exponential - Euler's Number - if (simplifiedMostSigOp->Is>()) { - u = simplifiedMostSigOp->Copy(); - dv = simplifiedLeastSigOp->Copy(); - } else if (simplifiedLeastSigOp->Is>()) { - u =simplifiedLeastSigOp->Copy(); - dv = simplifiedMostSigOp->Copy(); - } - + // In the form of mult{u, dv} + Multiply mult { + this->GetMostSigOp(), + this->GetLeastSigOp() + }; + // // Check the rules of LIPET + // Current code assumes MostSigOp is u and LeastSigOp is dv + // TODO: Need to rewrite if statements in the case it's necessary to swap operands + // // Logarithm + // if (simplifiedMostSigOp->Is>()) { + // u = simplifiedMostSigOp->Copy(); + // dv = simplifiedLeastSigOp->Copy(); + // } else if (simplifiedLeastSigOp->Is>()) { + // u = simplifiedLeastSigOp->Copy(); + // dv = simplifiedMostSigOp->Copy(); + // } + // + // // TODO: Inverse trigonometry + // // Inverse trigonometry is not implemented yet in Oasis + // + // // Polynomial + // // TODO: Could also be exponent, in the case of (x^2)*sinx + // // TODO: Ensure all polynomial cases are accounted for + // if (simplifiedMostSigOp->Is() || simplifiedMostSigOp->Is>()) { + // u = simplifiedMostSigOp->Copy(); + // dv = simplifiedLeastSigOp->Copy(); + // } else if (simplifiedLeastSigOp->Is() || simplifiedMostSigOp->Is>()) { + // u = simplifiedLeastSigOp->Copy(); + // dv = simplifiedMostSigOp->Copy(); + // } + // + // // Exponential - Euler's Number + // if (simplifiedMostSigOp->Is>()) { + // u = simplifiedMostSigOp->Copy(); + // dv = simplifiedLeastSigOp->Copy(); + // } else if (simplifiedLeastSigOp->Is>()) { + // u = simplifiedLeastSigOp->Copy(); + // dv = simplifiedMostSigOp->Copy(); + // } // TODO: Trigonometry // Trigonometry is not implemented yet in Oasis - - - // Differentiate u and integrate dv - // TODO: ensure u & dv are casted as the correct types to allow method override - std::unique_ptr du = u->Differentiate(integrationVariable); - std::unique_ptr v = dv->Integrate(integrationVariable); - - - Integral vdu { - Multiply {*v, *du}, - integrationVariable - }; + Multiply unsimplified_vdu { + *mult.GetMostSigOp().Differentiate(integrationVariable), + *mult.GetLeastSigOp().Differentiate(integrationVariable)}; // TODO: Make test cases in the event that IBP needs to be done again on integrated_vdu - auto integrated_vdu = vdu.Integrate(integrationVariable); + auto integrated_vdu = (unsimplified_vdu.Accept(simplifyVisitor).value())->Integrate(integrationVariable); // Apply the formula: integral(udv) = uv - integral(vdu) Subtract, Expression> subtractor { - Multiply { *u, *v }, + Multiply { mult.GetMostSigOp(), mult.GetLeastSigOp() }, *integrated_vdu }; + // TODO: Need to address the -C at the end of the integral, should be +C return subtractor.Accept(simplifyVisitor).value(); - - } - } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; From e2c032baffd7c0e09c107036f494a20e07f50a9d Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 20 Mar 2026 18:51:00 -0400 Subject: [PATCH 07/21] Added temporary fix for constant adjustment The first test case on integrating xe^x passed. However, if solving an integral by parts requires applying IBP again, then the constant will become negative. This was solved by multiplying the coefficient by -1. This does not work for higher levels of recursion. A more thorough fix that relies on recursion level tracking will be included in the next commit. --- include/Oasis/Multiply.hpp | 2 ++ src/Multiply.cpp | 16 +++++++++++++--- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/include/Oasis/Multiply.hpp b/include/Oasis/Multiply.hpp index fc0b1b9d..f4314fbf 100644 --- a/include/Oasis/Multiply.hpp +++ b/include/Oasis/Multiply.hpp @@ -22,6 +22,8 @@ class Multiply : public BinaryExpression { EXPRESSION_TYPE(Multiply) EXPRESSION_CATEGORY(Associative | Commutative | BinExp) +private: + [[nodiscard]] auto IntegrateByParts(const Expression& integrationVariable, double level) const -> std::unique_ptr; }; /// @endcond diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 5c5e7847..40908b8f 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -491,7 +491,6 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons // Detect whether integration by parts is appropriate // May need to simplify before and/or after else { - // In the form of mult{u, dv} Multiply mult { this->GetMostSigOp(), @@ -500,7 +499,7 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons // // Check the rules of LIPET // Current code assumes MostSigOp is u and LeastSigOp is dv - // TODO: Need to rewrite if statements in the case it's necessary to swap operands + // TODO: Need to rewrite if statements in the case it's necessary to swap operands of mult // // Logarithm // if (simplifiedMostSigOp->Is>()) { // u = simplifiedMostSigOp->Copy(); @@ -542,12 +541,23 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons // TODO: Make test cases in the event that IBP needs to be done again on integrated_vdu auto integrated_vdu = (unsimplified_vdu.Accept(simplifyVisitor).value())->Integrate(integrationVariable); + // Adjust the coefficient of the Constant variable "C" to be -1, since the quantity will be subtracted afterward + // TODO: change the exponent on (-1) to the level of recursion + if (auto adjusted_vdu = RecursiveCast>(*integrated_vdu); adjusted_vdu != nullptr) { + Add corrected_coefficient { + adjusted_vdu->GetMostSigOp(), + Multiply { Exponent { Real {-1}, Real { 1 }}, Variable { "C" } } + + }; + integrated_vdu = corrected_coefficient.Accept(simplifyVisitor).value(); + } + // Apply the formula: integral(udv) = uv - integral(vdu) Subtract, Expression> subtractor { Multiply { mult.GetMostSigOp(), mult.GetLeastSigOp() }, *integrated_vdu }; - // TODO: Need to address the -C at the end of the integral, should be +C + return subtractor.Accept(simplifyVisitor).value(); } From 176e08f993511fef76f4d90346342212e6ae2331 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 27 Mar 2026 16:30:14 -0400 Subject: [PATCH 08/21] Added a fix for constant "C" coefficient adjustment Implemented a solution for the coefficient of the variable C being negative. This previously occurred when subtracting a positive C variable. --- include/Oasis/Multiply.hpp | 2 +- src/Multiply.cpp | 97 ++++++++++++++++++++------------------ 2 files changed, 51 insertions(+), 48 deletions(-) diff --git a/include/Oasis/Multiply.hpp b/include/Oasis/Multiply.hpp index f4314fbf..84313d8d 100644 --- a/include/Oasis/Multiply.hpp +++ b/include/Oasis/Multiply.hpp @@ -23,7 +23,7 @@ class Multiply : public BinaryExpression { EXPRESSION_TYPE(Multiply) EXPRESSION_CATEGORY(Associative | Commutative | BinExp) private: - [[nodiscard]] auto IntegrateByParts(const Expression& integrationVariable, double level) const -> std::unique_ptr; + [[nodiscard]] auto Integrate(const Expression& integrationVariable, int recurseLevel) const -> std::unique_ptr; }; /// @endcond diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 40908b8f..1845dd76 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -459,7 +459,11 @@ auto Multiply::Simplify() const -> std::unique_ptr return simplifiedMultiply.Copy(); } -auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr +auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { + return Integrate(integrationVariable, 1); +} + +auto Multiply::Integrate(const Expression& integrationVariable, int recurseLevel) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; // Single integration variable @@ -488,77 +492,76 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons } // TODO: Finish implementation of integration by parts - // Detect whether integration by parts is appropriate - // May need to simplify before and/or after else { // In the form of mult{u, dv} - Multiply mult { + Multiply copy { this->GetMostSigOp(), this->GetLeastSigOp() }; - // // Check the rules of LIPET + auto mult = std::make_unique>(copy); // Current code assumes MostSigOp is u and LeastSigOp is dv - // TODO: Need to rewrite if statements in the case it's necessary to swap operands of mult - // // Logarithm - // if (simplifiedMostSigOp->Is>()) { - // u = simplifiedMostSigOp->Copy(); - // dv = simplifiedLeastSigOp->Copy(); - // } else if (simplifiedLeastSigOp->Is>()) { - // u = simplifiedLeastSigOp->Copy(); - // dv = simplifiedMostSigOp->Copy(); + // Check if by LIPET it's necessary to swap operands of mult + // Logarithm + // if ( !(mult->GetMostSigOp().Is>()) && mult->GetLeastSigOp().Is>()) { + // mult = std::make_unique>(mult->SwapOperands()); // } - // - // // TODO: Inverse trigonometry - // // Inverse trigonometry is not implemented yet in Oasis - // - // // Polynomial - // // TODO: Could also be exponent, in the case of (x^2)*sinx - // // TODO: Ensure all polynomial cases are accounted for - // if (simplifiedMostSigOp->Is() || simplifiedMostSigOp->Is>()) { - // u = simplifiedMostSigOp->Copy(); - // dv = simplifiedLeastSigOp->Copy(); - // } else if (simplifiedLeastSigOp->Is() || simplifiedMostSigOp->Is>()) { - // u = simplifiedLeastSigOp->Copy(); - // dv = simplifiedMostSigOp->Copy(); + + // TODO: Inverse trigonometry + // Inverse trigonometry is not implemented yet in Oasis + + // Polynomial + // TODO: Could also be exponent, in the case of (x^2)*sinx + // TODO: Ensure all polynomial cases are accounted for + // else if (mult->GetLeastSigOp().Is() || mult->GetLeastSigOp().Is>()) { + // mult = std::make_unique>(mult->SwapOperands()); // } - // - // // Exponential - Euler's Number - // if (simplifiedMostSigOp->Is>()) { - // u = simplifiedMostSigOp->Copy(); - // dv = simplifiedLeastSigOp->Copy(); - // } else if (simplifiedLeastSigOp->Is>()) { - // u = simplifiedLeastSigOp->Copy(); - // dv = simplifiedMostSigOp->Copy(); + + // Exponential - Euler's Number + // else if ( !(mult->GetMostSigOp().Is>()) && mult->GetLeastSigOp().Is>()) { + // mult = std::make_unique>(mult->SwapOperands()); // } // TODO: Trigonometry // Trigonometry is not implemented yet in Oasis Multiply unsimplified_vdu { - *mult.GetMostSigOp().Differentiate(integrationVariable), - *mult.GetLeastSigOp().Differentiate(integrationVariable)}; + *(mult->GetMostSigOp().Differentiate(integrationVariable)), + *(mult->GetLeastSigOp().Differentiate(integrationVariable))}; // TODO: Make test cases in the event that IBP needs to be done again on integrated_vdu - auto integrated_vdu = (unsimplified_vdu.Accept(simplifyVisitor).value())->Integrate(integrationVariable); - // Adjust the coefficient of the Constant variable "C" to be -1, since the quantity will be subtracted afterward - // TODO: change the exponent on (-1) to the level of recursion - if (auto adjusted_vdu = RecursiveCast>(*integrated_vdu); adjusted_vdu != nullptr) { - Add corrected_coefficient { - adjusted_vdu->GetMostSigOp(), - Multiply { Exponent { Real {-1}, Real { 1 }}, Variable { "C" } } + auto simplified_vdu = unsimplified_vdu.Accept(simplifyVisitor).value();; - }; - integrated_vdu = corrected_coefficient.Accept(simplifyVisitor).value(); + std::unique_ptr integrated_vdu; + + // Check if simplified_vdu would be integrated by parts again. + // If so, include the recursive level with the integration call + // Otherwise, call the normal integration function. + if (auto multiply = RecursiveCast>(*simplified_vdu); multiply != nullptr ) { + integrated_vdu = multiply->Integrate(integrationVariable, recurseLevel + 1); + } + else { + integrated_vdu = simplified_vdu->Integrate(integrationVariable); + } + + // Correct the coefficient of the constant variable "C", if necessary + // Subtracting a negative "C" would produce a positive "C" in the final answer with an odd final recurse level + if (recurseLevel % 2 == 1) { + if (auto adjusted_vdu = RecursiveCast>(*integrated_vdu); adjusted_vdu != nullptr) { + Subtract corrected_coefficient { + adjusted_vdu->GetMostSigOp(), + adjusted_vdu->GetLeastSigOp() + }; + integrated_vdu = corrected_coefficient.Accept(simplifyVisitor).value(); + } } // Apply the formula: integral(udv) = uv - integral(vdu) Subtract, Expression> subtractor { - Multiply { mult.GetMostSigOp(), mult.GetLeastSigOp() }, + Multiply { mult->GetMostSigOp(), mult->GetLeastSigOp() }, *integrated_vdu }; - return subtractor.Accept(simplifyVisitor).value(); } } From 4cdc578e376bf5de270b769d1a6da3d127c12cad Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 27 Mar 2026 17:34:45 -0400 Subject: [PATCH 09/21] Fixed IBP test case 2 Fixed the Integration By Parts test case with a polynomial and euler's number. The tree was not previously in accordance with the way expressions are simplified. It is now fixed. --- tests/IntegrateTests.cpp | 109 +++++++++++++++++++++------------------ 1 file changed, 58 insertions(+), 51 deletions(-) diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index f68b5728..e7d9562c 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -210,8 +210,8 @@ TEST_CASE("Integration By Parts: Variable and Euler's Number", "[Integrate][Vari auto integrated = integrand.Integrate(var); REQUIRE((integrated->Equals(*ptr))); - integrated = integrand.SwapOperands().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]") @@ -223,68 +223,75 @@ TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Expo Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } }; - - Oasis::Add, Oasis::Exponent>, - Oasis::Multiply>>, Oasis::Variable> integral { + Oasis::Add, + Oasis::Add, Oasis::Multiply, Oasis::Real>>>>, Oasis::Variable> integral { Oasis::Add { - Oasis::Subtract { - Oasis::Multiply { + Oasis::Multiply { + Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } }, + Oasis::Add { Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, - Oasis::Exponent { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } }, - Oasis::Multiply { - Oasis::Real { 2 }, - 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 Logarithm", "[Integrate][Exponent][Logarithm]") -{ - Oasis::Variable var { "x" }; - - Oasis::Multiply, Oasis::Log> integrand { - Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, - Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } - }; - - Oasis::Add, - Oasis::Exponent>, Oasis::Log>, - Oasis::Multiply, Oasis::Exponent>>, - Oasis::Variable> integral { - Oasis::Add { - Oasis::Subtract { - Oasis::Multiply { Oasis::Multiply { - Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 3 } }, - Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } + Oasis::Real {-1}, + Oasis::Add { + Oasis::Multiply { Oasis::Real { 2 }, Oasis::Variable { var.GetName() } }, + Oasis::Real { -2 } + } }, - Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } - }, - Oasis::Multiply { - Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 9 } }, - Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } } }, 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))); - 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::Log> integrand { +// Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, +// Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } +// }; +// +// Oasis::Add, +// Oasis::Exponent>, Oasis::Log>, +// Oasis::Multiply, Oasis::Exponent>>, +// Oasis::Variable> integral { +// Oasis::Add { +// Oasis::Subtract { +// Oasis::Multiply { +// Oasis::Multiply { +// Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 3 } }, +// Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } +// }, +// Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } +// }, +// Oasis::Multiply { +// Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 9 } }, +// Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } +// } +// }, +// 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))); +// } +// +// +// From f73fd86f93cc966763430de5553e7374ae895292 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 3 Apr 2026 16:20:38 -0400 Subject: [PATCH 10/21] Added detection to swap operands for LIPET Added functionality to detect whether operands should be swapped according to the rules of LIPET in determining which expression is u and which expression is dv. Also added a check to prevent infinite recursion. --- src/Multiply.cpp | 57 +++++++++++++++++++++++++++--------------------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 1845dd76..4a6cb90a 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -500,36 +500,37 @@ auto Multiply::Integrate(const Expression& integrationVariable, int }; auto mult = std::make_unique>(copy); - // Current code assumes MostSigOp is u and LeastSigOp is dv - // Check if by LIPET it's necessary to swap operands of mult - // Logarithm - // if ( !(mult->GetMostSigOp().Is>()) && mult->GetLeastSigOp().Is>()) { - // mult = std::make_unique>(mult->SwapOperands()); - // } - - // TODO: Inverse trigonometry - // Inverse trigonometry is not implemented yet in Oasis - - // Polynomial - // TODO: Could also be exponent, in the case of (x^2)*sinx - // TODO: Ensure all polynomial cases are accounted for - // else if (mult->GetLeastSigOp().Is() || mult->GetLeastSigOp().Is>()) { - // mult = std::make_unique>(mult->SwapOperands()); - // } - - // Exponential - Euler's Number - // else if ( !(mult->GetMostSigOp().Is>()) && mult->GetLeastSigOp().Is>()) { - // mult = std::make_unique>(mult->SwapOperands()); - // } - // TODO: Trigonometry - // Trigonometry is not implemented yet in Oasis + + // Assume MostSigOp is u and LeastSigOp is dv + // Check if by LIPET it's necessary to swap operands + + // LIPET: Logarithm case + if (!(mult->GetMostSigOp().Is>()) + && mult->GetLeastSigOp().Is>()) { + mult = std::make_unique>(mult->SwapOperands()); + } + + // LIPET: Inverse trigonometry is not implemented yet in Oasis + + // LIPET: Swap polynomial (linear variable) and EulerNumber + if ((mult->GetLeastSigOp().Is() + && mult->GetMostSigOp().Is>())) { + mult = std::make_unique>(mult->SwapOperands()); + } + + // LIPET: Swap polynomial (exponential) and EulerNumber + if (auto test1 = RecursiveCast>(mult->GetLeastSigOp()); test1 != nullptr) { + if (auto test2 = RecursiveCast>(mult->GetMostSigOp()); test2 != nullptr) { + mult = std::make_unique>(mult->SwapOperands()); + } + } + + // LIPET: Trigonometry is not implemented yet in Oasis Multiply unsimplified_vdu { *(mult->GetMostSigOp().Differentiate(integrationVariable)), *(mult->GetLeastSigOp().Differentiate(integrationVariable))}; - // TODO: Make test cases in the event that IBP needs to be done again on integrated_vdu - auto simplified_vdu = unsimplified_vdu.Accept(simplifyVisitor).value();; std::unique_ptr integrated_vdu; @@ -544,6 +545,12 @@ auto Multiply::Integrate(const Expression& integrationVariable, int integrated_vdu = simplified_vdu->Integrate(integrationVariable); } + // Prevent infinite recursion + // if an expression cannot be integrated, return the integral returned by the integration attempt + if (simplified_vdu->Equals(*integrated_vdu)) { + return integrated_vdu->Copy(); + } + // Correct the coefficient of the constant variable "C", if necessary // Subtracting a negative "C" would produce a positive "C" in the final answer with an odd final recurse level if (recurseLevel % 2 == 1) { From b72bd45e3be3fde2fd63a1de8db643b1ff302a49 Mon Sep 17 00:00:00 2001 From: bguyott Date: Sat, 11 Apr 2026 11:20:34 -0400 Subject: [PATCH 11/21] Fixed integration of dv and handling of constant C There was a typo where dv was differentiated when it should have been integrated. This commit remedies that, along with several bug fixes such as factoring out a constant of the integrand when applying IBP. The methodology for handling the constant C was also changed to be more reliable and consistent, removing the need for a recursion level counter. --- src/Multiply.cpp | 88 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 63 insertions(+), 25 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 4a6cb90a..437e1d6e 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -527,50 +527,88 @@ auto Multiply::Integrate(const Expression& integrationVariable, int // LIPET: Trigonometry is not implemented yet in Oasis - Multiply unsimplified_vdu { - *(mult->GetMostSigOp().Differentiate(integrationVariable)), - *(mult->GetLeastSigOp().Differentiate(integrationVariable))}; + auto v = mult->GetLeastSigOp().Integrate(integrationVariable); + auto du = mult->GetMostSigOp().Differentiate(integrationVariable); + + // Remove the +C from the integral + auto vPlusC = RecursiveCast>(*v); + v = vPlusC->GetMostSigOp().Copy(); // Accept(simplifyVisitor).value(); + + Real v_coefficient_local { 1 }; + Real du_coefficient_local { 1 }; + std::unique_ptr v_coefficient = v_coefficient_local.Copy(); + std::unique_ptr du_coefficient = du_coefficient_local.Copy(); + + // Make a copy of v so that the coefficients can be factored out + // when computing vdu + auto v_vdu = v->Copy(); + + // Factor out the coefficient for v + if (auto v_coeff = RecursiveCast>(*v); v_coeff != nullptr) { + v_coefficient = v_coeff->GetMostSigOp().Accept(simplifyVisitor).value(); + v_vdu = v_coeff->GetLeastSigOp().Accept(simplifyVisitor).value(); + } + + // Factor out the cofficient for du + if (auto du_coeff = RecursiveCast>(*du); du_coeff != nullptr) { + du_coefficient = du_coeff->GetMostSigOp().Accept(simplifyVisitor).value(); + du = du_coeff->GetLeastSigOp().Accept(simplifyVisitor).value(); + } + + // Attain one constant for the coefficient of vdu + Multiply vdu_coefficient { *v_coefficient, *du_coefficient }; + auto vdu_coefficient_simplified = vdu_coefficient.Accept(simplifyVisitor).value(); + + Multiply part_vdu { *v_vdu, *du }; + auto part_vdu_simplified = part_vdu.Accept(simplifyVisitor).value(); // std::move(part_vdu).Accept(simplifyVisitor).value(); + + auto simplified_vdu = part_vdu_simplified->Copy(); - auto simplified_vdu = unsimplified_vdu.Accept(simplifyVisitor).value();; std::unique_ptr integrated_vdu; // Check if simplified_vdu would be integrated by parts again. // If so, include the recursive level with the integration call // Otherwise, call the normal integration function. - if (auto multiply = RecursiveCast>(*simplified_vdu); multiply != nullptr ) { - integrated_vdu = multiply->Integrate(integrationVariable, recurseLevel + 1); - } - else { - integrated_vdu = simplified_vdu->Integrate(integrationVariable); + + // TODO: Below code may be redundant, remove after additional testing + // if (auto multiply = RecursiveCast>(*simplified_vdu); multiply != nullptr) { + // integrated_vdu = multiply->Integrate(integrationVariable, recurseLevel + 1); + // } else { + // integrated_vdu = simplified_vdu->Integrate(integrationVariable); + // } + // There is an implicit base case in that eventually, the integrand of vdu will not be a multiply instance + integrated_vdu = simplified_vdu->Integrate(integrationVariable); + + if (auto removeC = RecursiveCast>(*integrated_vdu); removeC != nullptr) { + integrated_vdu = removeC->GetMostSigOp().Copy();//Accept(simplifyVisitor).value(); } + // Apply coefficient correction for integrated_vdu + Multiply full_vdu { *vdu_coefficient_simplified, *integrated_vdu}; + + integrated_vdu = full_vdu.Copy(); + // Prevent infinite recursion // if an expression cannot be integrated, return the integral returned by the integration attempt + // TODO: Assess if necessary, given implicit recursion base case if (simplified_vdu->Equals(*integrated_vdu)) { return integrated_vdu->Copy(); } - // Correct the coefficient of the constant variable "C", if necessary - // Subtracting a negative "C" would produce a positive "C" in the final answer with an odd final recurse level - if (recurseLevel % 2 == 1) { - if (auto adjusted_vdu = RecursiveCast>(*integrated_vdu); adjusted_vdu != nullptr) { - Subtract corrected_coefficient { - adjusted_vdu->GetMostSigOp(), - adjusted_vdu->GetLeastSigOp() - }; - integrated_vdu = corrected_coefficient.Accept(simplifyVisitor).value(); - } - } - // Apply the formula: integral(udv) = uv - integral(vdu) - Subtract, Expression> subtractor { - Multiply { mult->GetMostSigOp(), mult->GetLeastSigOp() }, - *integrated_vdu + Add, Expression>, Variable> adder { + Subtract, Expression> { + Multiply { mult->GetMostSigOp(), *v }, + *integrated_vdu + }, + Variable {"C"} }; - return subtractor.Accept(simplifyVisitor).value(); + return adder.Accept(simplifyVisitor).value(); } + + } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; From 8b8eae78da984ab81b78ec431a1fc8f3fc8a3e85 Mon Sep 17 00:00:00 2001 From: bguyott Date: Sat, 11 Apr 2026 15:03:36 -0400 Subject: [PATCH 12/21] Fixed IBP tests 2 and 3 With the modified implementation of IBP in the previous commit, test cases two and three required alteration to their trees due to how SimplifyVisitor simplifies expressions. These test cases now pass. --- tests/IntegrateTests.cpp | 97 +++++++++++++++++++--------------------- 1 file changed, 45 insertions(+), 52 deletions(-) diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index e7d9562c..25183081 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -205,13 +205,12 @@ TEST_CASE("Integration By Parts: Variable and Euler's Number", "[Integrate][Vari } }; - auto ptr = integral.Accept(simplifyVisitor).value(); auto integrated = integrand.Integrate(var); REQUIRE((integrated->Equals(*ptr))); - // integrated = integrand.SwapOperands().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]") @@ -224,18 +223,18 @@ TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Expo }; Oasis::Add, - Oasis::Add, Oasis::Multiply, Oasis::Real>>>>, Oasis::Variable> integral { + Oasis::Add, Oasis::Multiply>>>, 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 {-1}, + Oasis::Real {-2}, Oasis::Add { - Oasis::Multiply { Oasis::Real { 2 }, Oasis::Variable { var.GetName() } }, - Oasis::Real { -2 } + Oasis::Variable { var.GetName() }, + Oasis::Real { -1 } } }, } @@ -248,50 +247,44 @@ TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Expo auto ptr = integral.Accept(simplifyVisitor).value(); auto integrated = integrand.Integrate(var); REQUIRE((integrated->Equals(*ptr))); -// -// integrated = integrand.SwapOperands().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::Log> integrand { + Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add, + Oasis::Add, + Oasis::Log>, Oasis::Divide>>, 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))); } -// TEST_CASE("Integration By Parts: Exponent and Logarithm", "[Integrate][Exponent][Logarithm]") -// { -// Oasis::Variable var { "x" }; -// -// Oasis::Multiply, Oasis::Log> integrand { -// Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2 } }, -// Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } -// }; -// -// Oasis::Add, -// Oasis::Exponent>, Oasis::Log>, -// Oasis::Multiply, Oasis::Exponent>>, -// Oasis::Variable> integral { -// Oasis::Add { -// Oasis::Subtract { -// Oasis::Multiply { -// Oasis::Multiply { -// Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 3 } }, -// Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } -// }, -// Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } -// }, -// Oasis::Multiply { -// Oasis::Divide { Oasis::Real { 1 }, Oasis::Real { 9 } }, -// Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3 } } -// } -// }, -// 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))); -// } -// -// -// From b1dcc6d200e97498a60f172b6699aa5e62b78b66 Mon Sep 17 00:00:00 2001 From: bguyott Date: Sat, 11 Apr 2026 16:00:24 -0400 Subject: [PATCH 13/21] Added Variable-Logarithm Integration Test --- tests/IntegrateTests.cpp | 38 ++++++++++++++++++++++++++++++++++++-- 1 file changed, 36 insertions(+), 2 deletions(-) diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index 25183081..797d581b 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -253,6 +253,41 @@ TEST_CASE("Integration By Parts: Exponent and Euler's Number", "[Integrate][Expo } +TEST_CASE("Integration By Parts: Variable and Logarithm", "[Integrate][Variable][Logarithm]") +{ + Oasis::Variable var { "x" }; + + Oasis::Multiply> integrand { + Oasis::Variable { var.GetName() }, + Oasis::Log { Oasis::EulerNumber{}, Oasis::Variable { var.GetName() } } + }; + + Oasis::Add, + Oasis::Add, + Oasis::Log>, Oasis::Divide>>, 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" }; @@ -286,5 +321,4 @@ TEST_CASE("Integration By Parts: Exponent and Logarithm", "[Integrate][Exponent] integrated = integrand.SwapOperands().Integrate(var); REQUIRE((integrated->Equals(*ptr))); -} - +} \ No newline at end of file From a910cdc7ab46bb30d1abd0c2141876fbb981f476 Mon Sep 17 00:00:00 2001 From: bguyott Date: Tue, 14 Apr 2026 17:16:37 -0400 Subject: [PATCH 14/21] Revised comments, variable names, and refactored code for simplicity --- include/Oasis/Multiply.hpp | 2 - src/Multiply.cpp | 140 +++++++++++++++---------------------- 2 files changed, 57 insertions(+), 85 deletions(-) diff --git a/include/Oasis/Multiply.hpp b/include/Oasis/Multiply.hpp index 84313d8d..fc0b1b9d 100644 --- a/include/Oasis/Multiply.hpp +++ b/include/Oasis/Multiply.hpp @@ -22,8 +22,6 @@ class Multiply : public BinaryExpression { EXPRESSION_TYPE(Multiply) EXPRESSION_CATEGORY(Associative | Commutative | BinExp) -private: - [[nodiscard]] auto Integrate(const Expression& integrationVariable, int recurseLevel) const -> std::unique_ptr; }; /// @endcond diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 437e1d6e..3a42ae02 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -459,11 +459,7 @@ auto Multiply::Simplify() const -> std::unique_ptr return simplifiedMultiply.Copy(); } -auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { - return Integrate(integrationVariable, 1); -} - -auto Multiply::Integrate(const Expression& integrationVariable, int recurseLevel) const -> std::unique_ptr +auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; // Single integration variable @@ -491,124 +487,102 @@ auto Multiply::Integrate(const Expression& integrationVariable, int } } - // TODO: Finish implementation of integration by parts + // Apply Integration By Parts else { - // In the form of mult{u, dv} - Multiply copy { - this->GetMostSigOp(), - this->GetLeastSigOp() - }; - - auto mult = std::make_unique>(copy); - - // Assume MostSigOp is u and LeastSigOp is dv - // Check if by LIPET it's necessary to swap operands + // Start by assuming u is MostSigOp and dv is LeastSigOp + auto multiplyCopy = std::make_unique>(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 (!(mult->GetMostSigOp().Is>()) - && mult->GetLeastSigOp().Is>()) { - mult = std::make_unique>(mult->SwapOperands()); + if (!(multiplyCopy->GetMostSigOp().Is>()) + && multiplyCopy->GetLeastSigOp().Is>()) { + multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); } // LIPET: Inverse trigonometry is not implemented yet in Oasis // LIPET: Swap polynomial (linear variable) and EulerNumber - if ((mult->GetLeastSigOp().Is() - && mult->GetMostSigOp().Is>())) { - mult = std::make_unique>(mult->SwapOperands()); + if ((multiplyCopy->GetLeastSigOp().Is() + && multiplyCopy->GetMostSigOp().Is>())) { + multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); } // LIPET: Swap polynomial (exponential) and EulerNumber - if (auto test1 = RecursiveCast>(mult->GetLeastSigOp()); test1 != nullptr) { - if (auto test2 = RecursiveCast>(mult->GetMostSigOp()); test2 != nullptr) { - mult = std::make_unique>(mult->SwapOperands()); + if (auto polynomial = RecursiveCast>(multiplyCopy->GetLeastSigOp()); polynomial != nullptr) { + if (auto euler = RecursiveCast>(multiplyCopy->GetMostSigOp()); euler != nullptr) { + multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); } } // LIPET: Trigonometry is not implemented yet in Oasis - auto v = mult->GetLeastSigOp().Integrate(integrationVariable); - auto du = mult->GetMostSigOp().Differentiate(integrationVariable); + // 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 the integral + // Remove the +C from v auto vPlusC = RecursiveCast>(*v); - v = vPlusC->GetMostSigOp().Copy(); // Accept(simplifyVisitor).value(); + v = vPlusC->GetMostSigOp().Copy(); - Real v_coefficient_local { 1 }; - Real du_coefficient_local { 1 }; - std::unique_ptr v_coefficient = v_coefficient_local.Copy(); - std::unique_ptr du_coefficient = du_coefficient_local.Copy(); + // Initialize the coefficients of v and du to 1 + auto vCoefficient = std::make_unique(Real {1}); + auto duCoefficient = std::make_unique(Real {1}); - // Make a copy of v so that the coefficients can be factored out - // when computing vdu - auto v_vdu = v->Copy(); + // 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 - if (auto v_coeff = RecursiveCast>(*v); v_coeff != nullptr) { - v_coefficient = v_coeff->GetMostSigOp().Accept(simplifyVisitor).value(); - v_vdu = v_coeff->GetLeastSigOp().Accept(simplifyVisitor).value(); + // Factor out the coefficient for v; v_VDU now has a coefficient of one + if (auto coefficientMultiply = RecursiveCast>(*v); coefficientMultiply != nullptr) { + vCoefficient = RecursiveCast(coefficientMultiply->GetMostSigOp()); + v_VDU = coefficientMultiply->GetLeastSigOp().Accept(simplifyVisitor).value(); } - // Factor out the cofficient for du - if (auto du_coeff = RecursiveCast>(*du); du_coeff != nullptr) { - du_coefficient = du_coeff->GetMostSigOp().Accept(simplifyVisitor).value(); - du = du_coeff->GetLeastSigOp().Accept(simplifyVisitor).value(); + // Factor out the coefficient for du; du now has a coefficient of one + if (auto coefficientMultiply = RecursiveCast>(*du); coefficientMultiply != nullptr) { + duCoefficient = RecursiveCast(coefficientMultiply->GetMostSigOp()); + du = coefficientMultiply->GetLeastSigOp().Accept(simplifyVisitor).value(); } - // Attain one constant for the coefficient of vdu - Multiply vdu_coefficient { *v_coefficient, *du_coefficient }; - auto vdu_coefficient_simplified = vdu_coefficient.Accept(simplifyVisitor).value(); - - Multiply part_vdu { *v_vdu, *du }; - auto part_vdu_simplified = part_vdu.Accept(simplifyVisitor).value(); // std::move(part_vdu).Accept(simplifyVisitor).value(); - - auto simplified_vdu = part_vdu_simplified->Copy(); - + // Attain a constant for the coefficient of vdu + auto vduCoefficient = std::make_unique>(Multiply{*vCoefficient, *duCoefficient})->Accept(simplifyVisitor).value(); - std::unique_ptr integrated_vdu; + // Multiply v and du to attain partVDU, which is vdu with a coefficient of one + auto partVDU = std::make_unique>(Multiply{*v_VDU, *du})->Accept(simplifyVisitor).value(); - // Check if simplified_vdu would be integrated by parts again. - // If so, include the recursive level with the integration call - // Otherwise, call the normal integration function. - - // TODO: Below code may be redundant, remove after additional testing - // if (auto multiply = RecursiveCast>(*simplified_vdu); multiply != nullptr) { - // integrated_vdu = multiply->Integrate(integrationVariable, recurseLevel + 1); - // } else { - // integrated_vdu = simplified_vdu->Integrate(integrationVariable); - // } - // There is an implicit base case in that eventually, the integrand of vdu will not be a multiply instance - integrated_vdu = simplified_vdu->Integrate(integrationVariable); - - if (auto removeC = RecursiveCast>(*integrated_vdu); removeC != nullptr) { - integrated_vdu = removeC->GetMostSigOp().Copy();//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{*vduCoefficient, *partVDU}); + if (multiplyCopy->Equals(*fullVDU)) { + Integral integral { *(multiplyCopy->Copy()), *(integrationVariable.Copy()) }; + return integral.Copy(); } - // Apply coefficient correction for integrated_vdu - Multiply full_vdu { *vdu_coefficient_simplified, *integrated_vdu}; - - integrated_vdu = full_vdu.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); - // Prevent infinite recursion - // if an expression cannot be integrated, return the integral returned by the integration attempt - // TODO: Assess if necessary, given implicit recursion base case - if (simplified_vdu->Equals(*integrated_vdu)) { - return integrated_vdu->Copy(); + // Remove the +C from the integration result + if (auto removeC = RecursiveCast>(*integratedPartVDU); removeC != nullptr) { + integratedPartVDU = removeC->GetMostSigOp().Copy(); } - // Apply the formula: integral(udv) = uv - integral(vdu) + // Apply coefficient correction for integratedPartVDU, thus building integratedFullVDU + auto integratedFullVDU = std::make_unique>(Multiply{*vduCoefficient, *integratedPartVDU}); + + // Apply the Integration By Parts formula, and add the +C at the end Add, Expression>, Variable> adder { Subtract, Expression> { - Multiply { mult->GetMostSigOp(), *v }, - *integrated_vdu + Multiply { multiplyCopy->GetMostSigOp(), *v }, + *integratedFullVDU }, Variable {"C"} }; return adder.Accept(simplifyVisitor).value(); } - - } Integral integral { *(this->Copy()), *(integrationVariable.Copy()) }; From c2e8418524213df66d5777e2a090f5b8335676b7 Mon Sep 17 00:00:00 2001 From: bguyott Date: Thu, 16 Apr 2026 14:10:44 -0400 Subject: [PATCH 15/21] Removed deprecated simplify method Removed a deprecated method before rebasing. --- src/Multiply.cpp | 440 ----------------------------------------------- 1 file changed, 440 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 3a42ae02..52e87edf 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -19,446 +19,6 @@ namespace Oasis { -// Mark for deletion: -// Keep in until merge to resolve other methods dependent on this one -auto Multiply::Simplify() const -> std::unique_ptr -{ - auto simplifiedMultiplicand = mostSigOp->Simplify(); - auto simplifiedMultiplier = leastSigOp->Simplify(); - - SimplifyVisitor simplifyVisitor; - - Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier }; - if (auto onezerocase = RecursiveCast>(simplifiedMultiply); onezerocase != nullptr) { - const Real& multiplicand = onezerocase->GetMostSigOp(); - const Expression& multiplier = onezerocase->GetLeastSigOp(); - if (std::abs(multiplicand.GetValue()) <= EPSILON) { - return std::make_unique(Real { 0.0 }); - } - if (multiplicand.GetValue() == 1) { - return multiplier.Simplify(); - } - } - if (auto realCase = RecursiveCast>(simplifiedMultiply); realCase != nullptr) { - const Real& multiplicand = realCase->GetMostSigOp(); - const Real& multiplier = realCase->GetLeastSigOp(); - return std::make_unique(multiplicand.GetValue() * multiplier.GetValue()); - } - - if (auto ImgCase = RecursiveCast>(simplifiedMultiply); ImgCase != nullptr) { - return std::make_unique(-1.0); - } - - if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { - auto e = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; - auto m = e.Accept(simplifyVisitor); - if (!m) { - return Divide { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); - } - return Divide { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); - } - - if (auto exprCase = RecursiveCast>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { - return std::make_unique>(exprCase->GetMostSigOp(), Real { 2.0 }); - } - } - - if (auto rMatrixCase = RecursiveCast>(simplifiedMultiply); rMatrixCase != nullptr) { - return std::make_unique(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()); - } - - if (auto matrixCase = RecursiveCast>(simplifiedMultiply); matrixCase != nullptr) { - const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp(); - const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp(); - - if (leftTerm.GetCols() == rightTerm.GetRows()) { - return std::make_unique(leftTerm.GetMatrix() * rightTerm.GetMatrix()); - } else { - // ERROR: INVALID DIMENSION - return std::make_unique>(leftTerm, rightTerm); - } - } - - // Commented out to not cause massive problems with things that need factored expressions - // // c*(a-b) - // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { - // if (negated->GetMostSigOp().GetValue()<0){ - // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, - // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); - // } else if (negated->GetMostSigOp().GetValue()>0){ - // return Subtract{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, - // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); - // } else { - // return Real{0}.Copy(); - // } - // - // } - // - // // c*(a+b) - // if (auto negated = RecursiveCast>>(simplifiedMultiply); negated != nullptr) { - // return Add{Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetMostSigOp()}, - // Multiply{negated->GetMostSigOp(), negated->GetLeastSigOp().GetLeastSigOp()}}.Simplify(); - // } - - if (auto multCase = RecursiveCast>>(simplifiedMultiply); multCase != nullptr) { - auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - if (!ms) { - return Divide { m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize(); - } - // This doesn't have to check errors because errors should have already been caught above - return Divide { *std::move(ms).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value(); - } - - if (auto multCase = RecursiveCast, Divide>>(simplifiedMultiply); multCase != nullptr) { - auto m = Multiply { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Multiply { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return this->Generalize(); - } - return Divide { *std::move(ms).value(), *std::move(ns).value() }.Accept(simplifyVisitor).value(); - } - - if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp(), e); - } - return std::make_unique>(exprCase->GetMostSigOp(), - *(std::move(s).value())); - } - } - - // x*x^n - if (auto exprCase = RecursiveCast>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp(), e); - } - return std::make_unique>(exprCase->GetMostSigOp(), *(std::move(s).value())); - } - } - - if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetLeastSigOp(), e); - } - return std::make_unique>(exprCase->GetLeastSigOp(), *(std::move(s).value())); - } - } - - // x^n*x^m - if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), e); - } - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s).value())); - } - } - - // a*x*x - if (auto exprCase = RecursiveCast, Expression>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }); - } - } - - // a*x*b*x - if (auto exprCase = RecursiveCast, Multiply>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return this->Generalize(); - } - return std::make_unique>( - *(std::move(ms).value()), - *(std::move(ns).value())); - } - } - - // a*x^n*x - if (auto exprCase = RecursiveCast>, Expression>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) { - auto e = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e }); - } - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - *(std::move(s).value()) }); - } - } - - // a*x*x^n - if (auto exprCase = RecursiveCast, Exponent>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), e }); - } - return std::make_unique>(exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - *(std::move(s).value()) }); - } - if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), - Exponent { exprCase->GetMostSigOp().GetMostSigOp(), e }); - } - return std::make_unique>(exprCase->GetMostSigOp().GetLeastSigOp(), - Exponent { exprCase->GetMostSigOp().GetMostSigOp(), - *(std::move(s).value()) }); - } - } - - // a*x^n*b*x - if (auto exprCase = RecursiveCast, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>( - m, Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), n }); - } - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - *(std::move(ns).value()) }); - } - } - - if (auto exprCase = RecursiveCast, Multiply, Expression>>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>(m, - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - n }); - } - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp(), - *(std::move(ns).value()) }); - } - } - - if (auto exprCase = RecursiveCast>, Multiply>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>(m, - Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), - n }); - } - - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetLeastSigOp().GetLeastSigOp(), - *(std::move(ns).value()) }); - } - } - - // a*x^n*x^m - if (auto exprCase = RecursiveCast>, Exponent>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) { - auto e = Add { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - return std::make_unique>( - exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - e }); - } - return std::make_unique>( - exprCase->GetMostSigOp().GetMostSigOp(), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - *(std::move(s).value()) }); - } - } - - // a*x^n*b*x^m - if (auto exprCase = RecursiveCast>, Multiply>>>(simplifiedMultiply); exprCase != nullptr) { - if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) { - auto m = Multiply { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }; - auto ms = m.Accept(simplifyVisitor); - auto n = Add { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }; - auto ns = n.Accept(simplifyVisitor); - if (!ms || !ns) { - return std::make_unique>( - m, - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - n }); - } - return std::make_unique>( - *(std::move(ms).value()), - Exponent { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), - *(std::move(ns).value()) }); - } - } - - // if (auto negate = RecursiveCast>>>(simplifiedMultiply); negate != nullptr){ - // return Add{Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetMostSigOp()}, - // Multiply{negate->GetMostSigOp(), negate->GetLeastSigOp().GetOperand().GetLeastSigOp()}}.Simplify(); - // } - - // multiply add like terms - std::vector> multiplies; - std::vector> vals; - simplifiedMultiply.Flatten(multiplies); - for (const auto& multiplicand : multiplies) { - size_t i = 0; - if (auto real = RecursiveCast(*multiplicand); real != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast(*vals[i]); valI != nullptr) { - vals[i] = Real { valI->GetValue() * real->GetValue() }.Generalize(); - break; - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - vals.push_back(multiplicand->Generalize()); - } - continue; - } - // single i - if (auto img = RecursiveCast(*multiplicand); img != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { Imaginary {}, e }.Generalize(); - } else { - vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); - } - - break; - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - vals.push_back(Exponent { Imaginary {}, Real { 1.0 } }.Generalize()); - } - continue; - } - // i^n - if (auto img = RecursiveCast>(*multiplicand); img != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - auto e = Add { valI->GetLeastSigOp(), img->GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { Imaginary {}, e }.Generalize(); - } else { - vals[i] = Exponent { Imaginary {}, *(std::move(s).value()) }.Generalize(); - } - break; - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); - vals.push_back(img->Generalize()); - } - continue; - } - // expr^n - if (auto expr = RecursiveCast>(*multiplicand); expr != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) { - auto e = Add { valI->GetLeastSigOp(), expr->GetLeastSigOp() }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); - } else { - vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); - } - break; - } - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - // vals.push_back(Multiply { img->GetMostSigOp(), Imaginary {} }.Generalize()); - vals.push_back(expr->Generalize()); - } - continue; - } - // single expr - if (auto expr = RecursiveCast(*multiplicand); expr != nullptr) { - for (; i < vals.size(); i++) { - if (auto valI = RecursiveCast>(*vals[i]); valI != nullptr) { - if (valI->GetMostSigOp().Equals(*expr)) { - auto e = Add { valI->GetLeastSigOp(), Real { 1.0 } }; - auto s = e.Accept(simplifyVisitor); - if (!s) { - vals[i] = Exponent { valI->GetMostSigOp(), e }.Generalize(); - } else { - vals[i] = Exponent { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize(); - } - break; - } - } - } - if (i >= vals.size()) { - // check to make sure it is one thing only - vals.push_back(Exponent { *expr, Real { 1.0 } }.Generalize()); - } - continue; - } - } - - // makes all expr^1 into expr - for (auto& val : vals) { - if (auto exp = RecursiveCast>(*val); exp != nullptr) { - if (exp->GetLeastSigOp().GetValue() == 1.0) { - val = exp->GetMostSigOp().Generalize(); - } - } - if (auto mul = RecursiveCast>(*val); mul != nullptr) { - if (mul->GetMostSigOp().GetValue() == 1.0) { - val = mul->GetLeastSigOp().Generalize(); - } - } - } - - return BuildFromVector(vals); - - return simplifiedMultiply.Copy(); -} - auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; From add033262d3433a99c5f361b8ea726c224eb412c Mon Sep 17 00:00:00 2001 From: bguyott Date: Thu, 16 Apr 2026 14:38:23 -0400 Subject: [PATCH 16/21] Removed comment for accuracy IBP is finished, removing the comment stating it was a work in progress. --- src/Integral.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Integral.cpp b/src/Integral.cpp index 92896833..ac901ad8 100644 --- a/src/Integral.cpp +++ b/src/Integral.cpp @@ -124,8 +124,6 @@ auto Integral::Simplify(const Expression& upper, const Expression& l // Integration by Parts // Needs Trig and Euler's Number - // Work in Progress - // Partial Fraction Decomposition // Work in Progress From 557fef4e8f80aa6edf8c678fd434a419a4f939ef Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 17 Apr 2026 16:36:47 -0400 Subject: [PATCH 17/21] Fixed Clang Format --- src/Exponent.cpp | 2 +- src/Multiply.cpp | 26 +++++++++++++------------- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/src/Exponent.cpp b/src/Exponent.cpp index 65f040c9..65e624d2 100644 --- a/src/Exponent.cpp +++ b/src/Exponent.cpp @@ -57,7 +57,7 @@ auto Exponent::Integrate(const Expression& integrationVariable) cons const Variable& expPow = eulerBase->GetLeastSigOp(); Add adder { - Exponent { EulerNumber{}, expPow }, + Exponent { EulerNumber {}, expPow }, Variable { "C" } }; diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 52e87edf..de33852f 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -50,12 +50,12 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons // Apply Integration By Parts else { // Start by assuming u is MostSigOp and dv is LeastSigOp - auto multiplyCopy = std::make_unique>(Multiply{this->GetMostSigOp(), this->GetLeastSigOp()}); + auto multiplyCopy = std::make_unique>(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>()) - && multiplyCopy->GetLeastSigOp().Is>()) { + && multiplyCopy->GetLeastSigOp().Is>()) { multiplyCopy = std::make_unique>(multiplyCopy->SwapOperands()); } @@ -85,8 +85,8 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons v = vPlusC->GetMostSigOp().Copy(); // Initialize the coefficients of v and du to 1 - auto vCoefficient = std::make_unique(Real {1}); - auto duCoefficient = std::make_unique(Real {1}); + auto vCoefficient = std::make_unique(Real { 1 }); + auto duCoefficient = std::make_unique(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 @@ -105,15 +105,15 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons } // Attain a constant for the coefficient of vdu - auto vduCoefficient = std::make_unique>(Multiply{*vCoefficient, *duCoefficient})->Accept(simplifyVisitor).value(); + auto vduCoefficient = std::make_unique>(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{*v_VDU, *du})->Accept(simplifyVisitor).value(); + auto partVDU = std::make_unique>(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{*vduCoefficient, *partVDU}); + auto fullVDU = std::make_unique>(Multiply { *vduCoefficient, *partVDU }); if (multiplyCopy->Equals(*fullVDU)) { Integral integral { *(multiplyCopy->Copy()), *(integrationVariable.Copy()) }; return integral.Copy(); @@ -130,15 +130,14 @@ auto Multiply::Integrate(const Expression& integrationVariable) cons } // Apply coefficient correction for integratedPartVDU, thus building integratedFullVDU - auto integratedFullVDU = std::make_unique>(Multiply{*vduCoefficient, *integratedPartVDU}); + auto integratedFullVDU = std::make_unique>(Multiply { *vduCoefficient, *integratedPartVDU }); // Apply the Integration By Parts formula, and add the +C at the end Add, Expression>, Variable> adder { - Subtract, Expression> { + Subtract, Expression> { Multiply { multiplyCopy->GetMostSigOp(), *v }, - *integratedFullVDU - }, - Variable {"C"} + *integratedFullVDU }, + Variable { "C" } }; return adder.Accept(simplifyVisitor).value(); @@ -183,7 +182,8 @@ auto Multiply::Differentiate(const Expression& differentiationVariab Multiply>>(Add, Multiply> { Multiply { *add, *right }, Multiply { *add2, *left } }) - ->Accept(simplifyVisitor).value(); + ->Accept(simplifyVisitor) + .value(); } } } From 73b4879e0d5ce7236d03ef262d53411b56ef9797 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 17 Apr 2026 16:59:07 -0400 Subject: [PATCH 18/21] Removed extraneous space in #includes --- src/Multiply.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index de33852f..7072f16c 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -9,10 +9,9 @@ #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/Log.hpp" #include "Oasis/Subtract.hpp" #define EPSILON 10E-6 From 7fb8777499d738a7766a6c2090cd3e59f7b56ca9 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 17 Apr 2026 17:02:51 -0400 Subject: [PATCH 19/21] Fixed codestyle formatting --- src/Multiply.cpp | 1 - tests/IntegrateTests.cpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Multiply.cpp b/src/Multiply.cpp index 7072f16c..dd522474 100644 --- a/src/Multiply.cpp +++ b/src/Multiply.cpp @@ -17,7 +17,6 @@ #define EPSILON 10E-6 namespace Oasis { - auto Multiply::Integrate(const Expression& integrationVariable) const -> std::unique_ptr { SimplifyVisitor simplifyVisitor {}; diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index 797d581b..b9a5b5fb 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -4,12 +4,12 @@ #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/SimplifyVisitor.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" -#include "Oasis/Log.hpp" inline Oasis::SimplifyVisitor simplifyVisitor{}; From 398c4bf8f7aaf0bfdb44d4e1fa54dc19d2213487 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 17 Apr 2026 17:05:18 -0400 Subject: [PATCH 20/21] Fixed format in #includes --- tests/IntegrateTests.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index b9a5b5fb..1f555a74 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -7,9 +7,9 @@ #include "Oasis/Log.hpp" #include "Oasis/Multiply.hpp" #include "Oasis/Real.hpp" -#include "Oasis/SimplifyVisitor.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" +#include "Oasis/SimplifyVisitor.hpp" inline Oasis::SimplifyVisitor simplifyVisitor{}; From ddea0294de75f52c256d6eefb7e27c6aff55b325 Mon Sep 17 00:00:00 2001 From: bguyott Date: Fri, 17 Apr 2026 17:17:56 -0400 Subject: [PATCH 21/21] Added EulerNumber integration test case --- tests/IntegrateTests.cpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp index 1f555a74..4ef48f57 100644 --- a/tests/IntegrateTests.cpp +++ b/tests/IntegrateTests.cpp @@ -185,6 +185,24 @@ TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]") 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::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" };