diff --git a/README.md b/README.md
index 5ff519b9..bef3ffcf 100644
--- a/README.md
+++ b/README.md
@@ -108,17 +108,26 @@ target_link_libraries(MyApp PRIVATE Oasis::Oasis)
 ```c++
 #include <cstdlib>
 #include <print> // C++23
+#include <utility>
+#include <memory>
 
 #include <Oasis/Add.hpp>
 #include <Oasis/Real.hpp>
+#include <Oasis/SimplifyVisitor.hpp>
 
 int main() {
+    SimplifyVisitor simplifyVisitor{};
     Oasis::Add sum {
         Oasis::Real{2.0},
         Oasis::Real{3.0}
     };
 
-    auto simplified = sum.Simplify();
+    auto expr = sum.Accept(simplifiedVisitor);
+    if (!expr) { // Error Case
+        std::println(expr.error());
+        return EXIT_FAILURE;
+    }
+    auto simplified = std::move(expr).value();
     std::println("Result: {}", simplified->ToString());
 
     return EXIT_SUCCESS;
diff --git a/include/CMakeLists.txt b/include/CMakeLists.txt
index b41736c7..faa95afa 100644
--- a/include/CMakeLists.txt
+++ b/include/CMakeLists.txt
@@ -20,6 +20,8 @@ set(Oasis_HEADERS
     Oasis/Pi.hpp
     Oasis/Real.hpp
     Oasis/RecursiveCast.hpp
+    Oasis/SimplifyVisitor.hpp
+    Oasis/Sine.hpp
     Oasis/Subtract.hpp
     Oasis/UnaryExpression.hpp
     Oasis/Undefined.hpp
diff --git a/include/Oasis/Add.hpp b/include/Oasis/Add.hpp
index b61e6d7f..56274e3d 100644
--- a/include/Oasis/Add.hpp
+++ b/include/Oasis/Add.hpp
@@ -20,7 +20,7 @@ class Add<
 public:
     using BinaryExpression::BinaryExpression;
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression> final;
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> final;
diff --git a/include/Oasis/BinaryExpression.hpp b/include/Oasis/BinaryExpression.hpp
index 30af5a1b..f7b56a89 100644
--- a/include/Oasis/BinaryExpression.hpp
+++ b/include/Oasis/BinaryExpression.hpp
@@ -11,6 +11,7 @@
 #include <list>
 
 #include "Expression.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
 #include "RecursiveCast.hpp"
 #include "Visit.hpp"
 
@@ -199,7 +200,13 @@ class BinaryExpression : public Expression {
 
     [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> override
     {
-        return Generalize()->Simplify();
+        SimplifyVisitor simplifyVisitor {};
+        auto e = Generalize();
+        auto s = e->Accept(simplifyVisitor);
+        if (!s) {
+            return e;
+        }
+        return std::move(s).value();
     }
 
     [[nodiscard]] auto Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression> override
@@ -361,11 +368,17 @@ class BinaryExpression : public Expression {
 
     auto Substitute(const Expression& var, const Expression& val) -> std::unique_ptr<Expression> override
     {
+        // TODO: FIX WITH VISITOR?
         std::unique_ptr<Expression> left = ((GetMostSigOp()).Copy())->Substitute(var, val);
         std::unique_ptr<Expression> right = ((GetLeastSigOp().Copy())->Substitute(var, val));
         DerivedT<Expression, Expression> comb = DerivedT<Expression, Expression> { *left, *right };
-        auto ret = comb.Simplify();
-        return ret;
+
+        Oasis::SimplifyVisitor simplifyVisitor {};
+        auto simplified = comb.Accept(simplifyVisitor);
+        if (!simplified) {
+            return comb.Generalize();
+        }
+        return std::move(simplified.value());
     }
     /**
      * Swaps the operands of this expression.
diff --git a/include/Oasis/Derivative.hpp b/include/Oasis/Derivative.hpp
index bab38b2c..f7bc925f 100644
--- a/include/Oasis/Derivative.hpp
+++ b/include/Oasis/Derivative.hpp
@@ -21,7 +21,7 @@ class Derivative<Expression, Expression> : public BinaryExpression<Derivative> {
 
     Derivative(const Expression& Exp, const Expression& Var);
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> override;
 
diff --git a/include/Oasis/Divide.hpp b/include/Oasis/Divide.hpp
index 19a412ab..96d720af 100644
--- a/include/Oasis/Divide.hpp
+++ b/include/Oasis/Divide.hpp
@@ -21,7 +21,7 @@ class Divide<Expression, Expression> : public BinaryExpression<Divide> {
 
     Divide(const Expression& dividend, const Expression& divisor);
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> final;
 
diff --git a/include/Oasis/Exponent.hpp b/include/Oasis/Exponent.hpp
index bf9a3d1c..84ebb207 100644
--- a/include/Oasis/Exponent.hpp
+++ b/include/Oasis/Exponent.hpp
@@ -22,7 +22,7 @@ class Exponent<Expression, Expression> : public BinaryExpression<Exponent> {
 
     Exponent(const Expression& base, const Expression& power);
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> final;
 
diff --git a/include/Oasis/Expression.hpp b/include/Oasis/Expression.hpp
index 2053875b..c41bf97c 100644
--- a/include/Oasis/Expression.hpp
+++ b/include/Oasis/Expression.hpp
@@ -38,7 +38,8 @@ enum class ExpressionType {
     Matrix,
     Pi,
     EulerNumber,
-    Magnitude
+    Magnitude,
+    Sine
 };
 
 /**
@@ -210,7 +211,7 @@ template <IVisitor T>
 auto Expression::Accept(T& visitor) const -> typename T::RetT
 {
     try {
-        return boost::any_cast<typename T::RetT>(this->AcceptInternal(visitor));
+        return boost::any_cast<typename T::RetT>(this->AcceptInternal(static_cast<Visitor&>(visitor)));
     } catch (boost::bad_any_cast& e) {
         return std::unexpected { e.what() };
     }
diff --git a/include/Oasis/Integral.hpp b/include/Oasis/Integral.hpp
index 5241648a..42d013ba 100644
--- a/include/Oasis/Integral.hpp
+++ b/include/Oasis/Integral.hpp
@@ -23,8 +23,9 @@ class Integral<Expression, Expression> : public BinaryExpression<Integral> {
 
     Integral(const Expression& integrand, const Expression& differential);
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
-    [[nodiscard]] auto Simplify(const Expression& upper, const Expression& lower) const -> std::unique_ptr<Expression> /* final */;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+
+    [[deprecated]] [[nodiscard]] auto Simplify(const Expression& upper, const Expression& lower) const -> std::unique_ptr<Expression> /* final */;
 
     EXPRESSION_TYPE(Integral)
     EXPRESSION_CATEGORY(Associative | Commutative)
diff --git a/include/Oasis/Log.hpp b/include/Oasis/Log.hpp
index 7c1a1f7d..31e3cfb1 100644
--- a/include/Oasis/Log.hpp
+++ b/include/Oasis/Log.hpp
@@ -22,7 +22,7 @@ class Log<Expression, Expression> : public BinaryExpression<Log> {
 
     Log(const Expression& base, const Expression& argument);
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression> final;
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> final;
diff --git a/include/Oasis/Magnitude.hpp b/include/Oasis/Magnitude.hpp
index df542a7d..6b45cf99 100644
--- a/include/Oasis/Magnitude.hpp
+++ b/include/Oasis/Magnitude.hpp
@@ -13,6 +13,7 @@
 #include "Imaginary.hpp"
 #include "Matrix.hpp"
 #include "Multiply.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
 #include "Real.hpp"
 #include "RecursiveCast.hpp"
 #include "UnaryExpression.hpp"
@@ -24,7 +25,7 @@ namespace Oasis {
  * @tparam OperandT Type of child operand
  * This represents magnitude/absolute value
  */
-template <typename OperandT>
+template <IExpression OperandT>
 class Magnitude final : public UnaryExpression<Magnitude, OperandT> {
 public:
     Magnitude() = default;
@@ -38,8 +39,9 @@ class Magnitude final : public UnaryExpression<Magnitude, OperandT> {
     {
     }
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> override
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> override
     {
+        SimplifyVisitor simplifyVisitor {};
         auto simpOp = this->GetOperand().Simplify();
         if (auto realCase = RecursiveCast<Real>(*simpOp); realCase != nullptr) {
             double val = realCase->GetValue();
@@ -49,7 +51,12 @@ class Magnitude final : public UnaryExpression<Magnitude, OperandT> {
             return std::make_unique<Real>(1.0);
         }
         if (auto mulImgCase = RecursiveCast<Multiply<Expression, Imaginary>>(*simpOp); mulImgCase != nullptr) {
-            return Magnitude<Expression> { mulImgCase->GetMostSigOp() }.Simplify();
+            auto e = Magnitude<Expression> { mulImgCase->GetMostSigOp() };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return e.Generalize();
+            }
+            return std::move(s).value();
         }
         if (auto addCase = RecursiveCast<Add<Expression, Imaginary>>(*simpOp); addCase != nullptr) {
             return Exponent { Add<Expression> { Exponent<Expression> { addCase->GetMostSigOp(), Real { 2 } },
@@ -79,12 +86,13 @@ class Magnitude final : public UnaryExpression<Magnitude, OperandT> {
     [[nodiscard]] auto Differentiate(const Expression& var) const -> std::unique_ptr<Expression> override
     {
         // TODO: Implement
+        Oasis::SimplifyVisitor simplifyVisitor {};
 
         const std::unique_ptr<Expression> operandDerivative = this->GetOperand().Differentiate(var);
         return Magnitude<Expression> {
             *operandDerivative
         }
-            .Simplify();
+            .Generalize();
     }
 
     [[nodiscard]] auto Integrate(const Expression& integrationVar) const -> std::unique_ptr<Expression> override
@@ -94,7 +102,7 @@ class Magnitude final : public UnaryExpression<Magnitude, OperandT> {
         return Magnitude<Expression> {
             *operandDerivative
         }
-            .Simplify();
+            .Generalize();
     }
 
     EXPRESSION_TYPE(Magnitude)
diff --git a/include/Oasis/Multiply.hpp b/include/Oasis/Multiply.hpp
index 14546868..701486a0 100644
--- a/include/Oasis/Multiply.hpp
+++ b/include/Oasis/Multiply.hpp
@@ -19,7 +19,7 @@ class Multiply<Expression, Expression> : public BinaryExpression<Multiply> {
 public:
     using BinaryExpression::BinaryExpression;
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> final;
 
diff --git a/include/Oasis/Negate.hpp b/include/Oasis/Negate.hpp
index e9c715ae..930573b4 100644
--- a/include/Oasis/Negate.hpp
+++ b/include/Oasis/Negate.hpp
@@ -10,7 +10,7 @@
 
 namespace Oasis {
 
-template <typename OperandT = Expression>
+template <IExpression OperandT = Expression>
 class Negate final : public UnaryExpression<Negate, OperandT> {
 public:
     Negate() = default;
@@ -24,7 +24,7 @@ class Negate final : public UnaryExpression<Negate, OperandT> {
     {
     }
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> override
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> override
     {
         return Multiply {
             Real { -1.0 },
@@ -39,7 +39,7 @@ class Negate final : public UnaryExpression<Negate, OperandT> {
         return Negate<Expression> {
             *operandDerivative
         }
-            .Simplify();
+            .Generalize(); // TODO: FIX WITH VISITOR
     }
 
     EXPRESSION_TYPE(Negate)
diff --git a/include/Oasis/SimplifyVisitor.hpp b/include/Oasis/SimplifyVisitor.hpp
new file mode 100644
index 00000000..4b85e63a
--- /dev/null
+++ b/include/Oasis/SimplifyVisitor.hpp
@@ -0,0 +1,57 @@
+//
+// Created by Andrew Nazareth on 9/23/25.
+//
+
+#ifndef SIMPLIFYVISITOR_HPP
+#define SIMPLIFYVISITOR_HPP
+
+#include <format>
+#include <string>
+
+#include <gsl/gsl-lite.hpp>
+
+#include "Oasis/Visit.hpp"
+
+namespace Oasis {
+
+struct SimplifyOpts {
+    enum class AngleUnits {
+        RADIANS,
+        DEGREES,
+    } angleUnits
+        = AngleUnits::RADIANS;
+};
+
+class SimplifyVisitor final : public TypedVisitor<std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string>> {
+public:
+    SimplifyVisitor();
+    explicit SimplifyVisitor(SimplifyOpts& opts);
+
+    auto TypedVisit(const Real& real) -> RetT override;
+    auto TypedVisit(const Imaginary& imaginary) -> RetT override;
+    auto TypedVisit(const Variable& variable) -> RetT override;
+    auto TypedVisit(const Undefined& undefined) -> RetT override;
+    auto TypedVisit(const Add<Expression, Expression>& add) -> RetT override;
+    auto TypedVisit(const Subtract<Expression, Expression>& subtract) -> RetT override;
+    auto TypedVisit(const Multiply<Expression, Expression>& multiply) -> RetT override;
+    auto TypedVisit(const Divide<Expression, Expression>& divide) -> RetT override;
+    auto TypedVisit(const Exponent<Expression, Expression>& exponent) -> RetT override;
+    auto TypedVisit(const Log<Expression, Expression>& log) -> RetT override;
+    auto TypedVisit(const Negate<Expression>& negate) -> RetT override;
+    auto TypedVisit(const Sine<Expression>& sine) -> RetT override;
+    auto TypedVisit(const Derivative<Expression, Expression>& derivative) -> RetT override;
+    auto TypedVisit(const Integral<Expression, Expression>& integral) -> RetT override;
+    auto TypedVisit(const Matrix& matrix) -> RetT override;
+    auto TypedVisit(const EulerNumber&) -> RetT override;
+    auto TypedVisit(const Pi&) -> RetT override;
+    auto TypedVisit(const Magnitude<Expression>& magnitude) -> RetT override;
+
+    [[nodiscard]] SimplifyOpts GetOptions() const;
+
+private:
+    SimplifyOpts options;
+};
+
+} // Oasis
+
+#endif // SIMPLIFYVISITOR_HPP
diff --git a/include/Oasis/Sine.hpp b/include/Oasis/Sine.hpp
new file mode 100644
index 00000000..1fad4778
--- /dev/null
+++ b/include/Oasis/Sine.hpp
@@ -0,0 +1,41 @@
+//
+// Created by Andrew Nazareth on 9/19/25.
+//
+
+#ifndef OASIS_SINE_HPP
+#define OASIS_SINE_HPP
+
+#include "UnaryExpression.hpp"
+
+namespace Oasis {
+
+template <IExpression OperandT>
+class Sine;
+
+template <>
+class Sine<Expression> final : public UnaryExpression<Sine> {
+public:
+    Sine() = default;
+    Sine(const Sine& other)
+        : UnaryExpression<Sine>(other)
+    {
+    }
+
+    explicit Sine(const Expression& operand)
+        : UnaryExpression<Sine>(operand)
+    {
+    }
+
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> override;
+
+    [[nodiscard]] auto Differentiate(const Expression& var) const -> std::unique_ptr<Expression> override;
+
+    [[nodiscard]] auto Integrate(const Expression& var) const -> std::unique_ptr<Expression> override;
+
+    EXPRESSION_TYPE(Sine)
+    EXPRESSION_CATEGORY(UnExp)
+};
+
+} // Oasis
+
+#endif // OASIS_SINE_HPP
diff --git a/include/Oasis/Subtract.hpp b/include/Oasis/Subtract.hpp
index c869dc51..fcf85181 100644
--- a/include/Oasis/Subtract.hpp
+++ b/include/Oasis/Subtract.hpp
@@ -22,7 +22,7 @@ class Subtract<Expression, Expression> : public BinaryExpression<Subtract> {
 
     Subtract(const Expression& minuend, const Expression& subtrahend);
 
-    [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
+    [[deprecated]] [[nodiscard]] auto Simplify() const -> std::unique_ptr<Expression> final;
 
     [[nodiscard]] auto Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression> final;
 
diff --git a/include/Oasis/UnaryExpression.hpp b/include/Oasis/UnaryExpression.hpp
index 5f5395ca..b1673f64 100644
--- a/include/Oasis/UnaryExpression.hpp
+++ b/include/Oasis/UnaryExpression.hpp
@@ -10,7 +10,7 @@
 
 namespace Oasis {
 
-template <template <IExpression> class DerivedT, IExpression OperandT>
+template <template <IExpression> class DerivedT, IExpression OperandT = Expression>
 class UnaryExpression : public Expression {
 
     using DerivedSpecialized = DerivedT<OperandT>;
@@ -87,8 +87,8 @@ class UnaryExpression : public Expression {
     {
         std::unique_ptr<Expression> right = ((GetOperand().Copy())->Substitute(var, val));
         DerivedT<Expression> comb = DerivedT<Expression> { *right };
-        auto ret = comb.Simplify();
-        return ret;
+        //        auto ret = comb.Accept();
+        return comb.Generalize();
     }
 
     auto AcceptInternal(Visitor& visitor) const -> any override
diff --git a/include/Oasis/Visit.hpp b/include/Oasis/Visit.hpp
index 12fa687d..af614570 100644
--- a/include/Oasis/Visit.hpp
+++ b/include/Oasis/Visit.hpp
@@ -35,10 +35,13 @@ class Exponent;
 template <IExpression BaseT, IExpression ArgumentT>
 class Log;
 
-template <typename OperandT>
+template <IExpression OperandT>
 class Negate;
 
-template <typename OperandT>
+template <IExpression OperandT>
+class Sine;
+
+template <IExpression OperandT>
 class Magnitude;
 
 template <IExpression, IExpression>
@@ -63,6 +66,7 @@ class Visitor {
     virtual any Visit(const Exponent<Expression, Expression>& exponent) = 0;
     virtual any Visit(const Log<Expression, Expression>& log) = 0;
     virtual any Visit(const Negate<Expression>& negate) = 0;
+    virtual any Visit(const Sine<Expression>& sine) = 0;
     virtual any Visit(const Magnitude<Expression>& magnitude) = 0;
     virtual any Visit(const Derivative<Expression, Expression>& derivative) = 0;
     virtual any Visit(const Integral<Expression, Expression>& integral) = 0;
@@ -89,6 +93,7 @@ class TypedVisitor : public Visitor {
     auto Visit(const Exponent<Expression, Expression>& exponent) -> any final { return TypedVisit(exponent); }
     auto Visit(const Log<Expression, Expression>& log) -> any final { return TypedVisit(log); }
     auto Visit(const Negate<Expression>& negate) -> any final { return TypedVisit(negate); }
+    auto Visit(const Sine<Expression>& sine) -> any final { return TypedVisit(sine); }
     auto Visit(const Magnitude<Expression>& magnitude) -> any final { return TypedVisit(magnitude); }
     auto Visit(const Derivative<Expression, Expression>& derivative) -> any final { return TypedVisit(derivative); }
     auto Visit(const Integral<Expression, Expression>& integral) -> any final { return TypedVisit(integral); }
@@ -108,6 +113,7 @@ class TypedVisitor : public Visitor {
     virtual auto TypedVisit(const Exponent<Expression, Expression>& exponent) -> RetT = 0;
     virtual auto TypedVisit(const Log<Expression, Expression>& log) -> RetT = 0;
     virtual auto TypedVisit(const Negate<Expression>& negate) -> RetT = 0;
+    virtual auto TypedVisit(const Sine<Expression>& sine) -> RetT = 0;
     virtual auto TypedVisit(const Magnitude<Expression>& magnitude) -> RetT = 0;
     virtual auto TypedVisit(const Derivative<Expression, Expression>& derivative) -> RetT = 0;
     virtual auto TypedVisit(const Integral<Expression, Expression>& integral) -> RetT = 0;
diff --git a/io/include/Oasis/InFixSerializer.hpp b/io/include/Oasis/InFixSerializer.hpp
index 76c23f87..1b832a35 100644
--- a/io/include/Oasis/InFixSerializer.hpp
+++ b/io/include/Oasis/InFixSerializer.hpp
@@ -25,6 +25,7 @@ class InFixSerializer final : public TypedVisitor<std::expected<std::string, std
     auto TypedVisit(const Exponent<Expression, Expression>& exponent) -> RetT override;
     auto TypedVisit(const Log<Expression, Expression>& log) -> RetT override;
     auto TypedVisit(const Negate<Expression>& negate) -> RetT override;
+    auto TypedVisit(const Sine<Expression>& sine) -> RetT override;
     auto TypedVisit(const Derivative<Expression, Expression>& derivative) -> RetT override;
     auto TypedVisit(const Integral<Expression, Expression>& integral) -> RetT override;
     auto TypedVisit(const Matrix& matrix) -> RetT override;
diff --git a/io/include/Oasis/MathMLSerializer.hpp b/io/include/Oasis/MathMLSerializer.hpp
index e3dc54c5..7f83d165 100644
--- a/io/include/Oasis/MathMLSerializer.hpp
+++ b/io/include/Oasis/MathMLSerializer.hpp
@@ -5,7 +5,6 @@
 #ifndef MATHMLSERIALIZER_HPP
 #define MATHMLSERIALIZER_HPP
 
-#define gsl_CONFIG_NOT_NULL_EXPLICIT_CTOR 0
 #include <gsl/gsl-lite.hpp>
 #include <tinyxml2.h>
 
@@ -33,6 +32,7 @@ class MathMLSerializer final : public TypedVisitor<std::expected<gsl::not_null<t
     auto TypedVisit(const Exponent<Expression, Expression>& exponent) -> RetT override;
     auto TypedVisit(const Log<Expression, Expression>& log) -> RetT override;
     auto TypedVisit(const Negate<Expression>& negate) -> RetT override;
+    auto TypedVisit(const Sine<Expression>& sine) -> RetT override;
     auto TypedVisit(const Magnitude<Expression>& magnitude) -> RetT override;
     auto TypedVisit(const Derivative<Expression, Expression>& derivative) -> RetT override;
     auto TypedVisit(const Integral<Expression, Expression>& integral) -> RetT override;
diff --git a/io/include/Oasis/TeXSerializer.hpp b/io/include/Oasis/TeXSerializer.hpp
index 1a2409dc..31432c2f 100644
--- a/io/include/Oasis/TeXSerializer.hpp
+++ b/io/include/Oasis/TeXSerializer.hpp
@@ -87,6 +87,7 @@ class TeXSerializer final : public TypedVisitor<std::expected<std::string, std::
     auto TypedVisit(const Exponent<Expression, Expression>& exponent) -> RetT override;
     auto TypedVisit(const Log<Expression, Expression>& log) -> RetT override;
     auto TypedVisit(const Negate<Expression>& negate) -> RetT override;
+    auto TypedVisit(const Sine<Expression>& sine) -> RetT override;
     auto TypedVisit(const Magnitude<Expression>& magnitude) -> RetT override;
     auto TypedVisit(const Derivative<Expression, Expression>& derivative) -> RetT override;
     auto TypedVisit(const Integral<Expression, Expression>& integral) -> RetT override;
diff --git a/io/src/InFixSerializer.cpp b/io/src/InFixSerializer.cpp
index 181fb6a4..276af614 100644
--- a/io/src/InFixSerializer.cpp
+++ b/io/src/InFixSerializer.cpp
@@ -16,6 +16,7 @@
 #include "Oasis/Multiply.hpp"
 #include "Oasis/Negate.hpp"
 #include "Oasis/Real.hpp"
+#include "Oasis/Sine.hpp"
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Variable.hpp"
 
@@ -78,6 +79,13 @@ auto InFixSerializer::TypedVisit(const Negate<Expression>& negate) -> RetT
     });
 }
 
+auto InFixSerializer::TypedVisit(const Sine<Expression>& sine) -> RetT
+{
+    return sine.GetOperand().Accept(*this).transform([](const std::string& str) {
+        return std::format("sine({})", str);
+    });
+}
+
 auto InFixSerializer::TypedVisit(const Derivative<>& derivative) -> RetT
 {
     return SerializeFuncBinExp(derivative, "dd");
diff --git a/io/src/MathMLSerializer.cpp b/io/src/MathMLSerializer.cpp
index 7ae1a710..706715e7 100644
--- a/io/src/MathMLSerializer.cpp
+++ b/io/src/MathMLSerializer.cpp
@@ -33,14 +33,14 @@ auto MathMLSerializer::TypedVisit(const Real& real) -> RetT
 {
     tinyxml2::XMLElement* result = doc.NewElement("mn");
     result->SetText(std::format("{:.5}", real.GetValue()).c_str());
-    return result;
+    return gsl::not_null(gsl::make_not_null(result));
 }
 
 auto MathMLSerializer::TypedVisit(const Imaginary&) -> RetT
 {
     tinyxml2::XMLElement* result = doc.NewElement("mi");
     result->SetText("i");
-    return result;
+    return gsl::not_null(result);
 }
 
 auto MathMLSerializer::TypedVisit(const Matrix& matrix) -> RetT
@@ -72,35 +72,35 @@ auto MathMLSerializer::TypedVisit(const Matrix& matrix) -> RetT
     mrow->InsertEndChild(table);
     mrow->InsertEndChild(closeBrace);
 
-    return mrow;
+    return gsl::not_null(mrow);
 }
 
 auto MathMLSerializer::TypedVisit(const Oasis::Pi&) -> RetT
 {
     tinyxml2::XMLElement* result = doc.NewElement("mi");
     result->SetText("&pi;");
-    return result;
+    return gsl::not_null(result);
 }
 
 auto MathMLSerializer::TypedVisit(const Oasis::EulerNumber&) -> RetT
 {
     tinyxml2::XMLElement* result = doc.NewElement("mi");
     result->SetText("e");
-    return result;
+    return gsl::not_null(result);
 }
 
 auto MathMLSerializer::TypedVisit(const Variable& variable) -> RetT
 {
     tinyxml2::XMLElement* result = doc.NewElement("mi");
     result->SetText(variable.GetName().c_str());
-    return result;
+    return gsl::not_null(result);
 }
 
 auto MathMLSerializer::TypedVisit(const Undefined&) -> RetT
 {
     tinyxml2::XMLElement* const mtext = doc.NewElement("mtext");
     mtext->SetText("Undefined");
-    return mtext;
+    return gsl::not_null(mtext);
 }
 
 auto MathMLSerializer::TypedVisit(const Add<>& add) -> RetT
@@ -122,7 +122,7 @@ auto MathMLSerializer::TypedVisit(const Add<>& add) -> RetT
 
     // remove last mo
     mrow->DeleteChild(mrow->LastChild());
-    return mrow;
+    return gsl::not_null(mrow);
 }
 
 auto MathMLSerializer::TypedVisit(const Subtract<>& subtract) -> RetT
@@ -165,7 +165,7 @@ auto MathMLSerializer::TypedVisit(const Subtract<>& subtract) -> RetT
             mrow->InsertEndChild(rightParen);
         }
 
-        return mrow;
+        return gsl::not_null(mrow);
     });
 }
 
@@ -226,7 +226,7 @@ auto MathMLSerializer::TypedVisit(const Multiply<>& multiply) -> RetT
         }
     }
 
-    return mrow;
+    return gsl::not_null(mrow);
 }
 
 auto MathMLSerializer::TypedVisit(const Divide<>& divide) -> RetT
@@ -239,7 +239,7 @@ auto MathMLSerializer::TypedVisit(const Divide<>& divide) -> RetT
         mfrac->InsertEndChild(dividendElement);
         mfrac->InsertEndChild(divisorElement);
 
-        return mfrac;
+        return gsl::not_null(mfrac);
     });
 }
 
@@ -272,7 +272,7 @@ auto MathMLSerializer::TypedVisit(const Exponent<>& exponent) -> RetT
 
         msup->InsertEndChild(powerElement);
 
-        return msup;
+        return gsl::not_null(msup);
     });
 }
 
@@ -306,7 +306,7 @@ auto MathMLSerializer::TypedVisit(const Log<>& log) -> RetT
         mrow->InsertEndChild(argElement);
         mrow->InsertEndChild(rightParen);
 
-        return mrow;
+        return gsl::not_null(mrow);
     });
 }
 
@@ -333,7 +333,14 @@ auto MathMLSerializer::TypedVisit(const Negate<Expression>& negate) -> RetT
     rightParen->SetText(")");
     mrow->InsertEndChild(rightParen);
 
-    return mrow;
+    return gsl::not_null(mrow);
+}
+
+auto MathMLSerializer::TypedVisit(const Sine<Expression>&) -> RetT
+{
+    // mrow
+    tinyxml2::XMLElement* const mrow = doc.NewElement("mrow");
+    return gsl::not_null(mrow);
 }
 
 auto MathMLSerializer::TypedVisit(const Derivative<>& derivative) -> RetT
@@ -372,7 +379,7 @@ auto MathMLSerializer::TypedVisit(const Derivative<>& derivative) -> RetT
         mrow->InsertEndChild(expElement);
         mrow->InsertEndChild(rightParen);
 
-        return mrow;
+        return gsl::not_null(mrow);
     });
 }
 
@@ -399,7 +406,7 @@ auto MathMLSerializer::TypedVisit(const Integral<>& integral) -> RetT
         mrow->InsertEndChild(expElement);
         mrow->InsertEndChild(dVar);
 
-        return mrow;
+        return gsl::not_null(mrow);
     });
 }
 
@@ -426,7 +433,7 @@ auto MathMLSerializer::TypedVisit(const Magnitude<Expression>& magnitude) -> Ret
     rightParen->SetText("|");
     mrow->InsertEndChild(rightParen);
 
-    return mrow;
+    return gsl::not_null(mrow);
 }
 
 }
\ No newline at end of file
diff --git a/io/src/TeXSerializer.cpp b/io/src/TeXSerializer.cpp
index 4d837d28..0fdc5332 100644
--- a/io/src/TeXSerializer.cpp
+++ b/io/src/TeXSerializer.cpp
@@ -19,6 +19,7 @@
 #include "Oasis/Negate.hpp"
 #include "Oasis/Pi.hpp"
 #include "Oasis/Real.hpp"
+#include "Oasis/Sine.hpp"
 #include "Oasis/Subtract.hpp"
 #include "Oasis/TeXSerializer.hpp"
 #include "Oasis/Variable.hpp"
@@ -186,6 +187,13 @@ auto TeXSerializer::TypedVisit(const Negate<Expression>& negate) -> RetT
     });
 }
 
+auto TeXSerializer::TypedVisit(const Sine<Expression>& sine) -> RetT
+{
+    return sine.GetOperand().Accept(*this).transform([](const std::string& str) {
+        return std::format("sine\\left({}\\right)", str);
+    });
+}
+
 auto TeXSerializer::TypedVisit(const Magnitude<Expression>& magnitude) -> RetT
 {
     return magnitude.GetOperand().Accept(*this).transform([](const std::string& str) {
diff --git a/src/Add.cpp b/src/Add.cpp
index 3b1f9a5a..5cdedd39 100644
--- a/src/Add.cpp
+++ b/src/Add.cpp
@@ -108,8 +108,16 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
         if (auto img = RecursiveCast<Imaginary>(*addend); img != nullptr) {
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Multiply<Expression, Imaginary>>(*vals[i]); valI != nullptr) {
-                    vals[i] = Multiply<Expression> { *(Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } }.Simplify()), Imaginary {} }.Generalize();
-                    break;
+                    auto e = Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } };
+                    SimplifyVisitor sv {};
+                    auto s = e.Accept(sv);
+                    if (!s) {
+                        vals[i] = Multiply<Expression> { e, Imaginary {} }.Generalize();
+                        break;
+                    } else {
+                        vals[i] = Multiply<Expression> { *(std::move(s).value()), Imaginary {} }.Generalize();
+                        break;
+                    }
                 }
             }
             if (i >= vals.size()) {
@@ -122,8 +130,16 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
         if (auto img = RecursiveCast<Multiply<Expression, Imaginary>>(*addend); img != nullptr) {
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Multiply<Expression, Imaginary>>(*vals[i]); valI != nullptr) {
-                    vals[i] = Multiply<Expression> { *(Add<Expression> { valI->GetMostSigOp(), img->GetMostSigOp() }.Simplify()), Imaginary {} }.Generalize();
-                    break;
+                    auto e = Add<Expression> { valI->GetMostSigOp(), img->GetMostSigOp() };
+                    SimplifyVisitor sv {};
+                    auto s = e.Accept(sv);
+                    if (!s) {
+                        vals[i] = Multiply<Expression> { e, Imaginary {} }.Generalize();
+                        break;
+                    } else {
+                        vals[i] = Multiply<Expression> { *(std::move(s).value()), Imaginary {} }.Generalize();
+                        break;
+                    }
                 }
             }
             if (i >= vals.size()) {
@@ -137,8 +153,16 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Multiply<Expression, Variable>>(*vals[i]); valI != nullptr) {
                     if (valI->GetLeastSigOp().GetName() == var->GetName()) {
-                        vals[i] = Multiply<Expression> { *(Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } }.Simplify()), *var }.Generalize();
-                        break;
+                        auto e = Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } };
+                        SimplifyVisitor sv {};
+                        auto s = e.Accept(sv);
+                        if (!s) {
+                            vals[i] = Multiply<Expression> { e, *var }.Generalize();
+                            break;
+                        } else {
+                            vals[i] = Multiply<Expression> { *(std::move(s).value()), *var }.Generalize();
+                            break;
+                        }
                     } else
                         continue;
                 }
@@ -155,8 +179,16 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
                 if (auto valI = RecursiveCast<Multiply<Expression, Variable>>(*vals[i]); valI != nullptr) {
                     // if (auto zeroCase = RecursiveCast<Multiply<Real, Expression>>(*valI); zeroCase != nullptr) {}
                     if (valI->GetLeastSigOp().GetName() == var->GetLeastSigOp().GetName()) {
-                        vals[i] = Multiply<Expression> { *(Add<Expression> { valI->GetMostSigOp(), var->GetMostSigOp() }.Simplify()), valI->GetLeastSigOp() }.Generalize();
-                        break;
+                        auto e = Add<Expression> { valI->GetMostSigOp(), var->GetMostSigOp() };
+                        SimplifyVisitor sv {};
+                        auto s = e.Accept(sv);
+                        if (!s) {
+                            vals[i] = Multiply<Expression> { e, valI->GetLeastSigOp() }.Generalize();
+                            break;
+                        } else {
+                            vals[i] = Multiply<Expression> { *(std::move(s).value()), valI->GetLeastSigOp() }.Generalize();
+                            break;
+                        }
                     } else
                         continue;
                 }
@@ -172,8 +204,16 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(*vals[i]); valI != nullptr) {
                     if (valI->GetLeastSigOp().Equals(*exp)) {
-                        vals[i] = Multiply<Expression> { *(Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } }.Simplify()), *exp }.Generalize();
-                        break;
+                        auto e = Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } };
+                        SimplifyVisitor sv {};
+                        auto s = e.Accept(sv);
+                        if (!s) {
+                            vals[i] = Multiply<Expression> { e, *exp }.Generalize();
+                            break;
+                        } else {
+                            vals[i] = Multiply<Expression> { *(std::move(s).value()), *exp }.Generalize();
+                            break;
+                        }
                     } else
                         continue;
                 }
@@ -189,8 +229,16 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(*vals[i]); valI != nullptr) {
                     if (valI->GetLeastSigOp().Equals(exp->GetLeastSigOp())) {
-                        vals[i] = Multiply<Expression> { *(Add<Expression> { valI->GetMostSigOp(), exp->GetMostSigOp() }.Simplify()), valI->GetLeastSigOp() }.Generalize();
-                        break;
+                        auto e = Add<Expression> { valI->GetMostSigOp(), exp->GetMostSigOp() };
+                        SimplifyVisitor sv {};
+                        auto s = e.Accept(sv);
+                        if (!s) {
+                            vals[i] = Multiply<Expression> { e, valI->GetLeastSigOp() }.Generalize();
+                            break;
+                        } else {
+                            vals[i] = Multiply<Expression> { *(std::move(s).value()), valI->GetLeastSigOp() }.Generalize();
+                            break;
+                        }
                     } else
                         continue;
                 }
@@ -238,8 +286,15 @@ auto Add<Expression>::Simplify() const -> std::unique_ptr<Expression>
 auto Add<Expression>::Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression>
 {
     // Single integration variable
+    SimplifyVisitor simplifyVisitor;
     if (auto variable = RecursiveCast<Variable>(integrationVariable); variable != nullptr) {
-        auto simplifiedAdd = this->Simplify();
+        auto simplified = this->Accept(simplifyVisitor);
+
+        if (!simplified) {
+            return this->Generalize();
+        }
+
+        auto simplifiedAdd = std::move(simplified).value();
 
         // Make sure we're still adder
         if (auto adder = RecursiveCast<Add<Expression>>(*simplifiedAdd); adder != nullptr) {
@@ -260,8 +315,11 @@ auto Add<Expression>::Integrate(const Expression& integrationVariable) const ->
                     *(specializedLeft->GetMostSigOp().Copy()), *(specializedRight->GetMostSigOp().Copy()) },
                 Variable { "C" }
             };
-
-            return add.Simplify();
+            auto postSimplify = add.Accept(simplifyVisitor);
+            if (!postSimplify) {
+                return add.Generalize();
+            }
+            return std::move(postSimplify).value();
         }
         // If not, use other integration technique
         else {
@@ -278,7 +336,12 @@ auto Add<Expression>::Differentiate(const Expression& differentiationVariable) c
     if (auto variable = RecursiveCast<Variable>(differentiationVariable); variable != nullptr) {
         auto left = mostSigOp->Differentiate(differentiationVariable);
         auto right = leastSigOp->Differentiate(differentiationVariable);
-        return Add<Expression> { *left, *right }.Simplify();
+        SimplifyVisitor simplifyVisitor;
+        auto simplified = Add<Expression> { *left, *right }.Accept(simplifyVisitor);
+        if (!simplified) {
+            Add<Expression> { *left, *right };
+        }
+        return std::move(simplified).value();
     }
     return Copy();
 }
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 87cbc173..b6e4a232 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -24,6 +24,8 @@ set(Oasis_SOURCES
     Negate.cpp
     Pi.cpp
     Real.cpp
+    SimplifyVisitor.cpp
+    Sine.cpp
     Subtract.cpp
     # Summation.cpp
     Undefined.cpp
diff --git a/src/Divide.cpp b/src/Divide.cpp
index 6c4057fe..dadb6640 100644
--- a/src/Divide.cpp
+++ b/src/Divide.cpp
@@ -24,9 +24,26 @@ Divide<Expression>::Divide(const Expression& dividend, const Expression& divisor
 
 auto Divide<Expression>::Simplify() const -> std::unique_ptr<Expression>
 {
-    auto simplifiedDividend = mostSigOp->Simplify(); // numerator
-    auto simplifiedDivider = leastSigOp->Simplify(); // denominator
-    Divide simplifiedDivide { *simplifiedDividend, *simplifiedDivider };
+    SimplifyVisitor simplifyVisitor {};
+    auto ms = mostSigOp->Accept(simplifyVisitor); // numerator
+    auto ls = leastSigOp->Accept(simplifyVisitor); // denominator
+
+    std::unique_ptr<Expression> simplifiedDividend;
+    std::unique_ptr<Expression> simplifiedDivider;
+
+    if (!ms) {
+        simplifiedDividend = mostSigOp->Copy();
+    } else {
+        simplifiedDividend = std::move(ms).value();
+    }
+
+    if (!ls) {
+        simplifiedDivider = mostSigOp->Copy();
+    } else {
+        simplifiedDivider = std::move(ms).value();
+    }
+
+    Divide simplifiedDivide = Divide { *simplifiedDividend, *simplifiedDivider };
 
     if (auto realCase = RecursiveCast<Divide<Real>>(simplifiedDivide); realCase != nullptr) {
         const Real& dividend = realCase->GetMostSigOp();
@@ -38,13 +55,23 @@ auto Divide<Expression>::Simplify() const -> std::unique_ptr<Expression>
     if (auto compCase = RecursiveCast<Divide<Divide, Expression>>(simplifiedDivide)) {
         const Expression& dividend = compCase->GetMostSigOp().GetMostSigOp();
         const Expression& divisor = Multiply { compCase->GetMostSigOp().GetLeastSigOp(), compCase->GetLeastSigOp() };
-        return std::make_unique<Divide<Expression, Expression>>(dividend, divisor)->Simplify();
+        auto mid = std::make_unique<Divide<Expression, Expression>>(dividend, divisor);
+        auto s = mid->Accept(simplifyVisitor);
+        if (!s) {
+            return mid;
+        }
+        return std::move(s.value());
     }
 
     if (auto compCase2 = RecursiveCast<Divide<Expression, Divide>>(simplifiedDivide)) {
         const Expression& dividend = Multiply { compCase2->GetLeastSigOp().GetLeastSigOp(), compCase2->GetMostSigOp() };
         const Expression& divisor = compCase2->GetLeastSigOp().GetMostSigOp();
-        return std::make_unique<Divide<Expression, Expression>>(dividend, divisor)->Simplify();
+        auto e = std::make_unique<Divide<Expression, Expression>>(dividend, divisor);
+        auto s = e->Accept(simplifyVisitor);
+        if (!e) {
+            return e;
+        }
+        return std::move(s.value());
     }
 
     // log(a)/log(b)=log[b](a)
@@ -110,7 +137,14 @@ auto Divide<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < result.size(); i++) {
                 if (auto resIexp = RecursiveCast<Exponent<Variable, Expression>>(*result[i]); resIexp != nullptr) {
                     if (resIexp->GetMostSigOp().Equals(*var)) {
-                        result[i] = Exponent<Expression> { *var, *(Subtract<Expression> { resIexp->GetLeastSigOp(), Real { 1.0 } }.Simplify()) }.Generalize();
+                        auto e = Subtract<Expression> { resIexp->GetLeastSigOp(), Real { 1.0 } };
+                        auto s = e.Accept(simplifyVisitor);
+                        if (!s) {
+                            result[i] = Exponent<Expression> { *var, e }.Generalize();
+                        } else {
+                            result[i] = Exponent<Expression> { *var, *(s.value()) }.Generalize();
+                        }
+
                         break;
                     }
                 } else if (auto resI = RecursiveCast<Variable>(*result[i]); resI != nullptr) {
@@ -129,12 +163,26 @@ auto Divide<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < result.size(); i++) {
                 if (auto resIexp = RecursiveCast<Exponent<Variable, Expression>>(*result[i]); resIexp != nullptr) {
                     if (resIexp->GetMostSigOp().Equals(var->GetMostSigOp())) {
-                        result[i] = Exponent<Expression> { var->GetMostSigOp(), *(Subtract<Expression> { resIexp->GetLeastSigOp(), var->GetLeastSigOp() }.Simplify()) }.Generalize();
+                        auto e = Subtract<Expression> { resIexp->GetLeastSigOp(), var->GetLeastSigOp() };
+                        auto s = e.Accept(simplifyVisitor);
+                        if (!s) {
+                            result[i] = Exponent<Expression> { var->GetMostSigOp(), e }.Generalize();
+                        } else {
+                            result[i] = Exponent<Expression> { var->GetMostSigOp(), *(s.value()) }.Generalize();
+                        }
+
                         break;
                     }
                 } else if (auto resI = RecursiveCast<Variable>(*result[i]); resI != nullptr) {
                     if (resI->Equals(var->GetMostSigOp())) {
-                        result[i] = Exponent<Expression> { var->GetMostSigOp(), *(Subtract<Expression> { Real { 1.0 }, var->GetLeastSigOp() }.Simplify()) }.Generalize();
+                        auto e = Subtract<Expression> { Real { 1.0 }, var->GetLeastSigOp() };
+                        auto s = e.Accept(simplifyVisitor);
+                        if (!s) {
+                            result[i] = Exponent<Expression> { var->GetMostSigOp(), e }.Generalize();
+                        } else {
+                            result[i] = Exponent<Expression> { var->GetMostSigOp(), *(s.value()) }.Generalize();
+                        }
+
                         break;
                     }
                 }
@@ -148,23 +196,43 @@ auto Divide<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < result.size(); i++) {
                 if (auto resExpr = RecursiveCast<Exponent<Expression, Expression>>(*result[i]); resExpr != nullptr) {
                     if (expExpr->GetMostSigOp().Equals(resExpr->GetMostSigOp())) {
-                        result[i] = Exponent { expExpr->GetMostSigOp(), *(Subtract { resExpr->GetLeastSigOp(), expExpr->GetLeastSigOp() }.Simplify()) }.Simplify();
+                        auto e = Subtract { resExpr->GetLeastSigOp(), expExpr->GetLeastSigOp() };
+                        auto s = e.Accept(simplifyVisitor);
+                        if (!s) {
+                            result[i] = Exponent { expExpr->GetMostSigOp(), e }.Accept(simplifyVisitor).value();
+                        } else {
+                            result[i] = Exponent { expExpr->GetMostSigOp(), *(s.value()) }.Accept(simplifyVisitor).value();
+                        }
+
                         break;
                     }
                 } else if (result[i]->Equals(expExpr->GetMostSigOp())) {
-                    result[i] = Exponent { expExpr->GetMostSigOp(), *(Subtract { Real { 1.0 }, expExpr->GetLeastSigOp() }.Simplify()) }.Simplify();
+                    auto e = Subtract { Real { 1.0 }, expExpr->GetLeastSigOp() };
+                    auto s = e.Accept(simplifyVisitor);
+                    if (!s) {
+                        result[i] = Exponent { expExpr->GetMostSigOp(), e }.Accept(simplifyVisitor).value();
+                    } else {
+                        result[i] = Exponent { expExpr->GetMostSigOp(), *(s.value()) }.Accept(simplifyVisitor).value();
+                    }
                     break;
                 }
             }
             if (i >= result.size()) {
-                result.push_back(Exponent<Expression> { expExpr->GetMostSigOp(), Multiply { Real { -1.0 }, expExpr->GetLeastSigOp() } }.Simplify());
+                result.push_back(Exponent<Expression> { expExpr->GetMostSigOp(), Multiply { Real { -1.0 }, expExpr->GetLeastSigOp() } }.Accept(simplifyVisitor).value());
             }
             continue;
         }
         for (; i < result.size(); i++) {
             if (auto resExpr = RecursiveCast<Exponent<Expression, Expression>>(*result[i]); resExpr != nullptr) {
                 if (denom->Equals(resExpr->GetMostSigOp())) {
-                    result[i] = Exponent { *denom, *(Subtract { resExpr->GetLeastSigOp(), Real { 1.0 } }.Simplify()) }.Simplify();
+                    auto e = Subtract { resExpr->GetLeastSigOp(), Real { 1.0 } };
+                    auto s = e.Accept(simplifyVisitor);
+                    if (!s) {
+                        result[i] = Exponent { *denom, e }.Accept(simplifyVisitor).value();
+                    } else {
+                        result[i] = Exponent { *denom, *(s.value()) }.Accept(simplifyVisitor).value();
+                    }
+
                     break;
                 }
             } else if (result[i]->Equals(*denom)) {
@@ -239,9 +307,10 @@ auto Divide<Expression>::Simplify() const -> std::unique_ptr<Expression>
 
 auto Divide<Expression>::Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // Single integration variable
     if (auto variable = RecursiveCast<Variable>(integrationVariable); variable != nullptr) {
-        auto simplifiedDiv = this->Simplify();
+        auto simplifiedDiv = this->Accept(simplifyVisitor).value();
 
         // Constant case - Integrand over a divisor
         if (auto constant = RecursiveCast<Multiply<Expression, Real>>(*simplifiedDiv); constant != nullptr) {
@@ -256,9 +325,10 @@ auto Divide<Expression>::Integrate(const Expression& integrationVariable) const
 
 auto Divide<Expression>::Differentiate(const Oasis::Expression& differentiationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // Single differentiation variable
     if (auto variable = RecursiveCast<Variable>(differentiationVariable); variable != nullptr) {
-        auto simplifiedDiv = this->Simplify();
+        auto simplifiedDiv = this->Accept(simplifyVisitor).value();
 
         // Constant case - differentiation over a divisor
         if (auto constant = RecursiveCast<Divide<Expression, Real>>(*simplifiedDiv); constant != nullptr) {
@@ -284,11 +354,11 @@ auto Divide<Expression>::Differentiate(const Oasis::Expression& differentiationV
             auto rightexp = quotient->GetLeastSigOp().Copy();
             auto leftDiff = leftexp->Differentiate(differentiationVariable);
             auto rightDiff = rightexp->Differentiate(differentiationVariable);
-            auto mult1 = Multiply<Expression, Expression>(Multiply<Expression, Expression> { *(rightexp->Simplify()), *(leftDiff->Simplify()) }).Simplify()->Simplify();
-            auto mult2 = Multiply<Expression, Expression>(Multiply<Expression, Expression> { *(leftexp->Simplify()), *(rightDiff->Simplify()) }).Simplify()->Simplify();
-            auto numerator = Subtract<Expression, Expression>(Subtract<Expression, Expression> { *mult1, *mult2 }).Simplify();
-            auto denominator = Multiply<Expression, Expression>(Multiply<Expression, Expression> { *(rightexp->Simplify()), *(rightexp->Simplify()) }).Simplify();
-            return Divide<Expression, Expression>({ *(numerator->Simplify()), *(denominator->Simplify()) }).Simplify();
+            auto mult1 = Multiply<Expression, Expression>(Multiply<Expression, Expression> { *(rightexp->Simplify()), *(leftDiff->Simplify()) }).Accept(simplifyVisitor).value()->Simplify();
+            auto mult2 = Multiply<Expression, Expression>(Multiply<Expression, Expression> { *(leftexp->Simplify()), *(rightDiff->Simplify()) }).Accept(simplifyVisitor).value()->Simplify();
+            auto numerator = Subtract<Expression, Expression>(Subtract<Expression, Expression> { *mult1, *mult2 }).Accept(simplifyVisitor).value();
+            auto denominator = Multiply<Expression, Expression>(Multiply<Expression, Expression> { *(rightexp->Simplify()), *(rightexp->Simplify()) }).Accept(simplifyVisitor).value();
+            return Divide<Expression, Expression>({ *(numerator->Simplify()), *(denominator->Simplify()) }).Accept(simplifyVisitor).value();
         }
     }
 
diff --git a/src/Exponent.cpp b/src/Exponent.cpp
index 8edef57a..a28e63d5 100644
--- a/src/Exponent.cpp
+++ b/src/Exponent.cpp
@@ -99,9 +99,16 @@ auto Exponent<Expression>::Simplify() const -> std::unique_ptr<Expression>
                                            .Case(
                                                [](const Exponent<Exponent<>, Expression>&) -> bool { return true; },
                                                [](const Exponent<Exponent<>, Expression>& expExpCase, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
-                                                   Exponent exp { expExpCase.GetMostSigOp().GetMostSigOp(),
-                                                       *(Multiply { expExpCase.GetMostSigOp().GetLeastSigOp(), expExpCase.GetLeastSigOp() }.Simplify()) };
-                                                   return gsl::make_not_null(exp.Copy());
+                                                   SimplifyVisitor simplifyVisitor;
+                                                   auto e = Multiply { expExpCase.GetMostSigOp().GetLeastSigOp(), expExpCase.GetLeastSigOp() };
+                                                   auto s = e.Accept(simplifyVisitor);
+                                                   if (!s) {
+                                                       Exponent exp { expExpCase.GetMostSigOp().GetMostSigOp(), e };
+                                                       return gsl::make_not_null(exp.Copy());
+                                                   } else {
+                                                       Exponent exp { expExpCase.GetMostSigOp().GetMostSigOp(), *(std::move(s.value())) };
+                                                       return gsl::make_not_null(exp.Copy());
+                                                   }
                                                })
                                            .Case(
                                                [](const Exponent<Expression, Log<>>& logCase) -> bool {
@@ -121,9 +128,14 @@ auto Exponent<Expression>::Simplify() const -> std::unique_ptr<Expression>
 
 auto Exponent<Expression>::Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // variable integration
     if (auto variable = RecursiveCast<Variable>(integrationVariable); variable != nullptr) {
-        auto simplifiedExponent = this->Simplify();
+        auto simplified = this->Accept(simplifyVisitor);
+        if (!simplified) {
+            return this->Generalize();
+        }
+        auto simplifiedExponent = std::move(simplified).value();
 
         std::unique_ptr<Expression> integral;
         // Variable with a constant power
@@ -153,8 +165,13 @@ auto Exponent<Expression>::Integrate(const Expression& integrationVariable) cons
 auto Exponent<Expression>::Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression>
 {
     // variable diff
+    SimplifyVisitor simplifyVisitor {};
     if (auto variable = RecursiveCast<Variable>(differentiationVariable); variable != nullptr) {
-        auto simplifiedExponent = this->Simplify();
+        auto simplified = this->Accept(simplifyVisitor);
+        if (!simplified) {
+            return this->Generalize();
+        }
+        auto simplifiedExponent = std::move(simplified).value();
 
         std::unique_ptr<Expression> diff;
         // Variable with a constant power
diff --git a/src/Expression.cpp b/src/Expression.cpp
index 2d863e33..8fe32290 100644
--- a/src/Expression.cpp
+++ b/src/Expression.cpp
@@ -249,17 +249,41 @@ auto Expression::Simplify() const -> std::unique_ptr<Expression>
 } // namespace Oasis
 std::unique_ptr<Oasis::Expression> operator+(const std::unique_ptr<Oasis::Expression>& lhs, const std::unique_ptr<Oasis::Expression>& rhs)
 {
-    return Oasis::Add { *lhs, *rhs }.Simplify();
+    Oasis::SimplifyVisitor sV {};
+    auto e = Oasis::Add { *lhs, *rhs };
+    auto s = e.Accept(sV);
+    if (!s) {
+        return e.Generalize();
+    }
+    return std::move(s).value();
 }
 std::unique_ptr<Oasis::Expression> operator-(const std::unique_ptr<Oasis::Expression>& lhs, const std::unique_ptr<Oasis::Expression>& rhs)
 {
-    return Oasis::Subtract { *lhs, *rhs }.Simplify();
+    Oasis::SimplifyVisitor sV {};
+    auto e = Oasis::Subtract { *lhs, *rhs };
+    auto s = e.Accept(sV);
+    if (!s) {
+        return e.Generalize();
+    }
+    return std::move(s).value();
 }
 std::unique_ptr<Oasis::Expression> operator*(const std::unique_ptr<Oasis::Expression>& lhs, const std::unique_ptr<Oasis::Expression>& rhs)
 {
-    return Oasis::Multiply { *lhs, *rhs }.Simplify();
+    Oasis::SimplifyVisitor sV {};
+    auto e = Oasis::Multiply { *lhs, *rhs };
+    auto s = e.Accept(sV);
+    if (!s) {
+        return e.Generalize();
+    }
+    return std::move(s).value();
 }
 std::unique_ptr<Oasis::Expression> operator/(const std::unique_ptr<Oasis::Expression>& lhs, const std::unique_ptr<Oasis::Expression>& rhs)
 {
-    return Oasis::Divide { *lhs, *rhs }.Simplify();
+    Oasis::SimplifyVisitor sV {};
+    auto e = Oasis::Divide { *lhs, *rhs };
+    auto s = e.Accept(sV);
+    if (!s) {
+        return e.Generalize();
+    }
+    return std::move(s).value();
 }
diff --git a/src/Linear.cpp b/src/Linear.cpp
index 39f2fae3..0d5bd40e 100644
--- a/src/Linear.cpp
+++ b/src/Linear.cpp
@@ -6,14 +6,17 @@
 #include "Oasis/Add.hpp"
 #include "Oasis/Multiply.hpp"
 #include "Oasis/RecursiveCast.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
 #include "Oasis/Variable.hpp"
 
 namespace Oasis {
 
 auto SolveLinearSystems(std::vector<std::unique_ptr<Expression>>& exprs) -> std::map<std::string, double>
 {
+    Oasis::SimplifyVisitor simplifyVisitor {};
     for (auto& expr : exprs) {
-        expr = expr->Simplify();
+        auto simplified = expr->Accept(simplifyVisitor);
+        expr = simplified.value()->Generalize();
     }
     auto matrices = ConstructMatrices(exprs);
     auto A = matrices.first.first;
diff --git a/src/Log.cpp b/src/Log.cpp
index 1a3520da..d026a3f9 100644
--- a/src/Log.cpp
+++ b/src/Log.cpp
@@ -17,6 +17,7 @@
 #include "Oasis/Multiply.hpp"
 #include "Oasis/Pi.hpp"
 #include "Oasis/RecursiveCast.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Undefined.hpp"
 
@@ -78,10 +79,17 @@ auto Log<Expression>::Simplify() const -> std::unique_ptr<Expression>
 
     // log[a](b^x) = x * log[a](b)
     if (const auto expCase = RecursiveCast<Log<Expression, Exponent<>>>(simplifiedLog); expCase != nullptr) {
+        SimplifyVisitor simplifyVisitor {};
         const auto exponent = expCase->GetLeastSigOp();
         const IExpression auto& log = Log<Expression>(expCase->GetMostSigOp(), exponent.GetMostSigOp()); // might need to check that it isnt nullptr
         const IExpression auto& factor = exponent.GetLeastSigOp();
-        return Oasis::Multiply<Oasis::Expression>(factor, log).Simplify();
+        auto e = Oasis::Multiply<Oasis::Expression>(factor, log);
+        auto s = e.Accept(simplifyVisitor);
+        if (!s) {
+            return e.Generalize();
+        } else {
+            return std::move(s.value());
+        }
     }
 
     return simplifiedLog.Copy();
@@ -89,7 +97,8 @@ auto Log<Expression>::Simplify() const -> std::unique_ptr<Expression>
 
 auto Log<Expression>::Integrate(const Oasis::Expression& integrationVariable) const -> std::unique_ptr<Expression>
 {
-    // TODO: Implement
+    // TODO: Implement with integrate visitor?
+    SimplifyVisitor simplifyVisitor {};
     if (this->mostSigOp->Equals(EulerNumber {})) {
         // ln(x)
         if (leastSigOp->Is<Variable>() && RecursiveCast<Variable>(*leastSigOp)->Equals(integrationVariable)) {
@@ -97,7 +106,12 @@ auto Log<Expression>::Integrate(const Oasis::Expression& integrationVariable) co
             // alternate form
             // return Subtract<Expression> { Multiply<Expression> { integrationVariable, *this }, integrationVariable }.Simplify();
         }
-        auto derivative = Derivative<Expression> { *leastSigOp, integrationVariable }.Simplify();
+        auto der_Simp = Derivative<Expression> { *leastSigOp, integrationVariable }.Accept(simplifyVisitor);
+        if (!der_Simp) {
+            return this->Generalize();
+        }
+
+        auto derivative = std::move(der_Simp).value();
 
         if (derivative->Equals(Real { 0 }))
             return Add { Multiply { integrationVariable, *this }, Variable { "C" } }.Simplify();
@@ -138,8 +152,11 @@ auto Log<Expression>::Integrate(const Oasis::Expression& integrationVariable) co
 
     } else {
         // Use log identity Log[b](a) = Log(a)/Log(b)
+        SimplifyVisitor simplifyVisitor {};
         auto numer = Log<Expression> { EulerNumber {}, *(this->leastSigOp->Generalize()) };
         auto denom = Log<Expression> { EulerNumber {}, *(this->mostSigOp->Generalize()) };
+
+        // TODO: FIX
         if (numer.Equals(denom))
             return Add { integrationVariable, Variable { "C" } }.Generalize();
 
@@ -158,24 +175,37 @@ auto Log<Expression>::Integrate(const Oasis::Expression& integrationVariable) co
 
 auto Log<Expression>::Differentiate(const Oasis::Expression& differentiationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // d(log_e(6x))/dx = 1/6x * 6
     if (auto lnCase = RecursiveCast<EulerNumber>(*mostSigOp); lnCase != nullptr) {
         Divide derivative { Oasis::Real { 1.0 }, *leastSigOp };
         Derivative chain { *leastSigOp, differentiationVariable };
 
         Multiply result = Multiply<Expression> { derivative, *chain.Differentiate(differentiationVariable) };
-        return result.Simplify();
+        auto simp = result.Accept(simplifyVisitor);
+        if (!simp) {
+            return result.Generalize();
+        }
+        return std::move(simp).value();
     } else {
         if (auto RealCase = RecursiveCast<Real>(*mostSigOp); RealCase != nullptr) {
             Divide derivative { Oasis::Real { 1.0 }, Multiply<Expression> { *leastSigOp, Log { EulerNumber {}, *mostSigOp } } };
             Derivative chain { *leastSigOp, differentiationVariable };
             Multiply result = Multiply<Expression> { derivative, *chain.Differentiate(differentiationVariable) };
-            return result.Simplify();
+            auto simp = result.Accept(simplifyVisitor);
+            if (!simp) {
+                return result.Generalize();
+            }
+            return std::move(simp).value();
         } else { // Use log identity and Quotient rule
             Divide result { Subtract { Multiply { Log { EulerNumber {}, *mostSigOp }, Derivative { Log { EulerNumber {}, *leastSigOp }, differentiationVariable } },
                                 Multiply { Log { EulerNumber {}, *leastSigOp }, Derivative { Log { EulerNumber {}, *mostSigOp }, differentiationVariable } } },
                 Exponent { Log { EulerNumber {}, *mostSigOp }, Real { 2 } } };
-            return result.Simplify();
+            auto simp = result.Accept(simplifyVisitor);
+            if (!simp) {
+                return result.Generalize();
+            }
+            return std::move(simp).value();
         }
     }
 }
diff --git a/src/Multiply.cpp b/src/Multiply.cpp
index 607d7126..8cb9f1be 100644
--- a/src/Multiply.cpp
+++ b/src/Multiply.cpp
@@ -20,6 +20,8 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
     auto simplifiedMultiplicand = mostSigOp->Simplify();
     auto simplifiedMultiplier = leastSigOp->Simplify();
 
+    SimplifyVisitor simplifyVisitor;
+
     Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier };
     if (auto onezerocase = RecursiveCast<Multiply<Real, Expression>>(simplifiedMultiply); onezerocase != nullptr) {
         const Real& multiplicand = onezerocase->GetMostSigOp();
@@ -42,8 +44,12 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
     }
 
     if (auto multCase = RecursiveCast<Multiply<Real, Divide<Expression>>>(simplifiedMultiply); multCase != nullptr) {
-        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }.Simplify();
-        return Divide<Expression> { *m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize();
+        auto e = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() };
+        auto m = e.Accept(simplifyVisitor);
+        if (!m) {
+            return Divide<Expression> { e, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize();
+        }
+        return Divide<Expression> { *std::move(m).value(), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize();
     }
 
     if (auto exprCase = RecursiveCast<Multiply<Expression>>(simplifiedMultiply); exprCase != nullptr) {
@@ -90,43 +96,70 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
     //    }
 
     if (auto multCase = RecursiveCast<Multiply<Expression, Divide<Expression>>>(simplifiedMultiply); multCase != nullptr) {
-        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }.Simplify();
-        return Divide<Expression> { *m, (multCase->GetLeastSigOp().GetLeastSigOp()) }.Simplify();
+        auto m = Multiply<Expression> { 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<Expression> { *std::move(ms.value()), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(simplifyVisitor).value();
     }
 
     if (auto multCase = RecursiveCast<Multiply<Divide<Expression>, Divide<Expression>>>(simplifiedMultiply); multCase != nullptr) {
-        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }.Simplify();
-        auto n = Multiply<Expression> { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }.Simplify();
-        return Divide<Expression> { *m, *n }.Simplify();
+        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() };
+        auto ms = m.Accept(simplifyVisitor);
+        auto n = Multiply<Expression> { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() };
+        auto ns = n.Accept(simplifyVisitor);
+        if (!ms || !ns) {
+            return this->Generalize();
+        }
+        return Divide<Expression> { *std::move(ms.value()), *std::move(ns.value()) }.Accept(simplifyVisitor).value();
     }
 
     if (auto exprCase = RecursiveCast<Multiply<Expression, Exponent<Expression, Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto e = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(), e);
+            }
             return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(),
-                *(Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()));
+                *(std::move(s.value())));
         }
     }
 
     // x*x^n
     if (auto exprCase = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
-            return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(),
-                *(Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()));
+            auto e = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(), e);
+            }
+            return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(), *(std::move(s.value())));
         }
     }
 
     if (auto exprCase = RecursiveCast<Multiply<Exponent<Expression>, Expression>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) {
-            return std::make_unique<Exponent<Expression>>(exprCase->GetLeastSigOp(),
-                *(Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()));
+            auto e = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Exponent<Expression>>(exprCase->GetLeastSigOp(), e);
+            }
+            return std::make_unique<Exponent<Expression>>(exprCase->GetLeastSigOp(), *(std::move(s.value())));
         }
     }
 
     // x^n*x^m
     if (auto exprCase = RecursiveCast<Multiply<Exponent<Expression>, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
-            return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
-                *(Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }.Simplify()));
+            auto e = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(), e);
+            }
+            return std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(), *(std::move(s.value())));
         }
     }
 
@@ -141,80 +174,150 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
     // a*x*b*x
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Multiply<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) {
+            auto m = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() };
+            auto ms = m.Accept(simplifyVisitor);
+            auto n = Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } };
+            auto ns = n.Accept(simplifyVisitor);
+            if (!ms || !ns) {
+                return this->Generalize();
+            }
             return std::make_unique<Multiply<Expression>>(
-                *(Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }.Simplify()),
-                *(Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }.Simplify()));
+                *(std::move(ms.value())),
+                *(std::move(ns.value())));
         }
     }
 
     // a*x^n*x
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Expression>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) {
+            auto e = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
+                    Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), e });
+            }
             return std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
-                    *(Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()) });
+                    *(std::move(s.value())) });
         }
     }
 
     // a*x*x^n
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto e = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
+                    Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), e });
+            }
             return std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(),
-                    *(Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()) });
+                    *(std::move(s.value())) });
         }
         if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto e = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetLeastSigOp(),
+                    Exponent<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), e });
+            }
             return std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetLeastSigOp(),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetMostSigOp(),
-                    *(Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()) });
+                    *(std::move(s.value())) });
         }
     }
 
     // a*x^n*b*x
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Multiply<Expression, Exponent<Expression>>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) {
+            auto m = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() };
+            auto ms = m.Accept(simplifyVisitor);
+            auto n = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto ns = n.Accept(simplifyVisitor);
+            if (!ms || !ns) {
+                return std::make_unique<Multiply<Expression>>(
+                    m, Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), n });
+            }
             return std::make_unique<Multiply<Expression>>(
-                *(Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }.Simplify()),
+                *(std::move(ms.value())),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(),
-                    *(Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()) });
+                    *(std::move(ns.value())) });
         }
     }
 
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Multiply<Exponent<Expression>, Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) {
+            auto m = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() };
+            auto ms = m.Accept(simplifyVisitor);
+            auto n = Add<Expression> { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto ns = n.Accept(simplifyVisitor);
+            if (!ms || !ns) {
+                return std::make_unique<Multiply<Expression>>(m,
+                    Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(),
+                        n });
+            }
             return std::make_unique<Multiply<Expression>>(
-                *(Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }.Simplify()),
+                *(std::move(ms.value())),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(),
-                    *(Add<Expression> { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()) });
+                    *(std::move(ns.value())) });
         }
     }
 
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Multiply<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) {
+            auto m = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() };
+            auto ms = m.Accept(simplifyVisitor);
+            auto n = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } };
+            auto ns = n.Accept(simplifyVisitor);
+            if (!ms || !ns) {
+                return std::make_unique<Multiply<Expression>>(m,
+                    Exponent<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(),
+                        n });
+            }
+
             return std::make_unique<Multiply<Expression>>(
-                *(Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }.Simplify()),
+                *(std::move(ms.value())),
                 Exponent<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(),
-                    *(Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Simplify()) });
+                    *(std::move(ns.value())) });
         }
     }
 
     // a*x^n*x^m
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) {
+            auto e = Add<Expression> { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() };
+            auto s = e.Accept(simplifyVisitor);
+            if (!s) {
+                return std::make_unique<Multiply<Expression>>(
+                    exprCase->GetMostSigOp().GetMostSigOp(),
+                    Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
+                        e });
+            }
             return std::make_unique<Multiply<Expression>>(
                 exprCase->GetMostSigOp().GetMostSigOp(),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
-                    *(Add<Expression> { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }.Simplify()) });
+                    *(std::move(s.value())) });
         }
     }
 
     // a*x^n*b*x^m
     if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Multiply<Expression, Exponent<Expression>>>>(simplifiedMultiply); exprCase != nullptr) {
         if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) {
+            auto m = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() };
+            auto ms = m.Accept(simplifyVisitor);
+            auto n = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() };
+            auto ns = n.Accept(simplifyVisitor);
+            if (!ms || !ns) {
+                return std::make_unique<Multiply<Expression>>(
+                    m,
+                    Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
+                        n });
+            }
             return std::make_unique<Multiply<Expression>>(
-                *(Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }.Simplify()),
+                *(std::move(ms.value())),
                 Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
-                    *(Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }.Simplify()) });
+                    *(std::move(ns.value())) });
         }
     }
 
@@ -246,7 +349,14 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
         if (auto img = RecursiveCast<Imaginary>(*multiplicand); img != nullptr) {
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Exponent<Imaginary, Expression>>(*vals[i]); valI != nullptr) {
-                    vals[i] = Exponent<Expression> { Imaginary {}, *(Add<Expression> { valI->GetLeastSigOp(), Real { 1.0 } }.Simplify()) }.Generalize();
+                    auto e = Add<Expression> { valI->GetLeastSigOp(), Real { 1.0 } };
+                    auto s = e.Accept(simplifyVisitor);
+                    if (!s) {
+                        vals[i] = Exponent<Expression> { Imaginary {}, e }.Generalize();
+                    } else {
+                        vals[i] = Exponent<Expression> { Imaginary {}, *(std::move(s).value()) }.Generalize();
+                    }
+
                     break;
                 }
             }
@@ -260,7 +370,13 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
         if (auto img = RecursiveCast<Exponent<Imaginary, Expression>>(*multiplicand); img != nullptr) {
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Exponent<Imaginary, Expression>>(*vals[i]); valI != nullptr) {
-                    vals[i] = Exponent<Expression> { Imaginary {}, *(Add<Expression> { valI->GetLeastSigOp(), img->GetLeastSigOp() }.Simplify()) }.Generalize();
+                    auto e = Add<Expression> { valI->GetLeastSigOp(), img->GetLeastSigOp() };
+                    auto s = e.Accept(simplifyVisitor);
+                    if (!s) {
+                        vals[i] = Exponent<Expression> { Imaginary {}, e }.Generalize();
+                    } else {
+                        vals[i] = Exponent<Expression> { Imaginary {}, *(std::move(s).value()) }.Generalize();
+                    }
                     break;
                 }
             }
@@ -276,7 +392,13 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Exponent<Expression, Expression>>(*vals[i]); valI != nullptr) {
                     if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) {
-                        vals[i] = Exponent<Expression> { valI->GetMostSigOp(), *(Add<Expression> { valI->GetLeastSigOp(), expr->GetLeastSigOp() }.Simplify()) }.Generalize();
+                        auto e = Add<Expression> { valI->GetLeastSigOp(), expr->GetLeastSigOp() };
+                        auto s = e.Accept(simplifyVisitor);
+                        if (!s) {
+                            vals[i] = Exponent<Expression> { valI->GetMostSigOp(), e }.Generalize();
+                        } else {
+                            vals[i] = Exponent<Expression> { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize();
+                        }
                         break;
                     }
                 }
@@ -293,7 +415,13 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
             for (; i < vals.size(); i++) {
                 if (auto valI = RecursiveCast<Exponent<Expression, Expression>>(*vals[i]); valI != nullptr) {
                     if (valI->GetMostSigOp().Equals(*expr)) {
-                        vals[i] = Exponent<Expression> { valI->GetMostSigOp(), *(Add<Expression> { valI->GetLeastSigOp(), Real { 1.0 } }.Simplify()) }.Generalize();
+                        auto e = Add<Expression> { valI->GetLeastSigOp(), Real { 1.0 } };
+                        auto s = e.Accept(simplifyVisitor);
+                        if (!s) {
+                            vals[i] = Exponent<Expression> { valI->GetMostSigOp(), e }.Generalize();
+                        } else {
+                            vals[i] = Exponent<Expression> { valI->GetMostSigOp(), *(std::move(s).value()) }.Generalize();
+                        }
                         break;
                     }
                 }
@@ -327,9 +455,14 @@ auto Multiply<Expression>::Simplify() const -> std::unique_ptr<Expression>
 
 auto Multiply<Expression>::Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // Single integration variable
     if (auto variable = RecursiveCast<Variable>(integrationVariable); variable != nullptr) {
-        auto simplifiedMult = this->Simplify();
+        auto s = this->Accept(simplifyVisitor);
+        if (!s) {
+            return this->Generalize();
+        }
+        auto simplifiedMult = std::move(s).value();
 
         // Constant case - Constant number multiplied by integrand
         if (auto constant = RecursiveCast<Multiply<Real, Expression>>(*simplifiedMult); constant != nullptr) {
@@ -357,9 +490,14 @@ auto Multiply<Expression>::Integrate(const Expression& integrationVariable) cons
 
 auto Multiply<Expression>::Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // Single integration variable
     if (auto variable = RecursiveCast<Variable>(differentiationVariable); variable != nullptr) {
-        auto simplifiedMult = this->Simplify();
+        auto s = this->Accept(simplifyVisitor);
+        if (!s) {
+            return this->Generalize();
+        }
+        auto simplifiedMult = std::move(s).value();
 
         // Constant case - Constant number multiplied by differentiate
         if (auto constant = RecursiveCast<Multiply<Real, Expression>>(*simplifiedMult); constant != nullptr) {
diff --git a/src/SimplifyVisitor.cpp b/src/SimplifyVisitor.cpp
new file mode 100644
index 00000000..3dc4dd4e
--- /dev/null
+++ b/src/SimplifyVisitor.cpp
@@ -0,0 +1,1470 @@
+//
+// Created by Andrew Nazareth on 9/23/25.
+//
+
+#include <format>
+
+#include "Oasis/SimplifyVisitor.hpp"
+
+#include "Oasis/Add.hpp"
+#include "Oasis/Derivative.hpp"
+#include "Oasis/Divide.hpp"
+#include "Oasis/EulerNumber.hpp"
+#include "Oasis/Exponent.hpp"
+#include "Oasis/Integral.hpp"
+#include "Oasis/Log.hpp"
+#include "Oasis/Magnitude.hpp"
+#include "Oasis/MatchCast.hpp"
+#include "Oasis/Multiply.hpp"
+#include "Oasis/Negate.hpp"
+#include "Oasis/Pi.hpp"
+#include "Oasis/Real.hpp"
+#include "Oasis/Sine.hpp"
+#include "Oasis/Subtract.hpp"
+#include "Oasis/Undefined.hpp"
+#include "Oasis/Variable.hpp"
+
+#define EPSILON 1E-6
+
+namespace Oasis {
+
+SimplifyVisitor::SimplifyVisitor()
+{
+    options = SimplifyOpts {};
+}
+
+SimplifyVisitor::SimplifyVisitor(SimplifyOpts& opts)
+    : options(opts)
+{
+}
+
+SimplifyOpts SimplifyVisitor::GetOptions() const
+{
+    return options;
+}
+
+auto SimplifyVisitor::TypedVisit(const Real& real) -> RetT
+{
+    return gsl::not_null { std::make_unique<Real>(real) };
+}
+
+auto SimplifyVisitor::TypedVisit(const Imaginary&) -> RetT
+{
+    return gsl::not_null { std::make_unique<Imaginary>() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Variable& variable) -> RetT
+{
+    return gsl::not_null { std::make_unique<Variable>(variable) };
+}
+
+auto SimplifyVisitor::TypedVisit(const Undefined&) -> RetT
+{
+    return gsl::not_null { std::make_unique<Undefined>() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Add<>& add) -> RetT
+{
+    auto mostSigOp = add.GetMostSigOp().Copy();
+    auto leastSigOp = add.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedAugendResult = mostSigOp->Accept(*this);
+    auto simplifiedAddendResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedAugendResult) {
+        return std::unexpected { simplifiedAugendResult.error() };
+    }
+    if (!simplifiedAddendResult) {
+        return std::unexpected { simplifiedAddendResult.error() };
+    }
+
+    auto simplifiedAugend = std::move(simplifiedAugendResult).value();
+    auto simplifiedAddend = std::move(simplifiedAddendResult).value();
+
+    Add simplifiedAdd { *simplifiedAugend, *simplifiedAddend };
+
+    if (auto realCase = RecursiveCast<Add<Real>>(simplifiedAdd); realCase != nullptr) {
+        const Real& firstReal = realCase->GetMostSigOp();
+        const Real& secondReal = realCase->GetLeastSigOp();
+
+        return gsl::not_null { std::make_unique<Real>(firstReal.GetValue() + secondReal.GetValue()) };
+    }
+
+    if (auto zeroCase = RecursiveCast<Add<Real, Expression>>(simplifiedAdd); zeroCase != nullptr) {
+        if (zeroCase->GetMostSigOp().GetValue() == 0) {
+            return gsl::not_null { zeroCase->GetLeastSigOp().Generalize() };
+        }
+    }
+
+    if (auto likeTermsCase = RecursiveCast<Add<Multiply<Real, Expression>>>(simplifiedAdd); likeTermsCase != nullptr) {
+        const Oasis::IExpression auto& leftTerm = likeTermsCase->GetMostSigOp().GetLeastSigOp();
+        const Oasis::IExpression auto& rightTerm = likeTermsCase->GetLeastSigOp().GetLeastSigOp();
+
+        if (leftTerm.Equals(rightTerm)) {
+            const Real& coefficient1 = likeTermsCase->GetMostSigOp().GetMostSigOp();
+            const Real& coefficient2 = likeTermsCase->GetLeastSigOp().GetMostSigOp();
+
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(Real(coefficient1.GetValue() + coefficient2.GetValue()), leftTerm) };
+        }
+    }
+
+    // matrix + matrix
+    if (auto matrixCase = RecursiveCast<Add<Matrix, Matrix>>(simplifiedAdd); matrixCase != nullptr) {
+        const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp();
+        const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp();
+
+        if ((leftTerm.GetRows() == rightTerm.GetRows()) && (leftTerm.GetCols() == rightTerm.GetCols())) {
+            return gsl::not_null { std::make_unique<Matrix>(leftTerm.GetMatrix() + rightTerm.GetMatrix()) };
+        } else {
+            return gsl::not_null { std::make_unique<Add<Expression>>(leftTerm, rightTerm) };
+        }
+    }
+
+    // log(a) + log(b) = log(ab)
+    if (auto logCase = RecursiveCast<Add<Log<Expression, Expression>, Log<Expression, Expression>>>(simplifiedAdd); logCase != nullptr) {
+        if (logCase->GetMostSigOp().GetMostSigOp().Equals(logCase->GetLeastSigOp().GetMostSigOp())) {
+            const IExpression auto& base = logCase->GetMostSigOp().GetMostSigOp();
+            const IExpression auto& argument = Multiply<Expression>({ logCase->GetMostSigOp().GetLeastSigOp(), logCase->GetLeastSigOp().GetLeastSigOp() });
+            return gsl::not_null { std::make_unique<Log<Expression>>(base, argument) };
+        }
+    }
+
+    // x + x = 2x
+    if (simplifiedAugend->Equals(*simplifiedAddend)) {
+        return Multiply<Real, Expression> { Real { 2.0 }, *simplifiedAugend }.Accept(*this);
+    }
+
+    // 2x + x = 3x
+    if (const auto likeTermsCase2 = RecursiveCast<Add<Multiply<Real, Expression>, Expression>>(simplifiedAdd); likeTermsCase2 != nullptr) {
+        if (likeTermsCase2->GetMostSigOp().GetLeastSigOp().Equals(likeTermsCase2->GetLeastSigOp())) {
+            const Real& coeffiecent = likeTermsCase2->GetMostSigOp().GetMostSigOp();
+            return gsl::not_null { std::make_unique<Multiply<Real, Expression>>(Real { coeffiecent.GetValue() + 1 }, likeTermsCase2->GetMostSigOp().GetLeastSigOp()) };
+        }
+    }
+
+    // simplifies expressions and combines like terms
+    // ex: 1 + 2x + 3 + 5x = 4 + 7x (or 7x + 4)
+    std::vector<std::unique_ptr<Expression>> adds;
+    std::vector<std::unique_ptr<Expression>> vals;
+    simplifiedAdd.Flatten(adds);
+    for (const auto& addend : adds) {
+        // real
+        size_t i = 0;
+        if (auto real = RecursiveCast<Real>(*addend); real != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Real>(*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(addend->Generalize());
+            }
+            continue;
+        }
+        // single i
+        if (auto img = RecursiveCast<Imaginary>(*addend); img != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Multiply<Expression, Imaginary>>(*vals[i]); valI != nullptr) {
+                    auto addResult = Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } }.Accept(*this);
+                    if (!addResult) {
+                        return addResult;
+                    }
+                    vals[i] = Multiply<Expression> { *(addResult.value()), Imaginary {} }.Generalize();
+                    break;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(Multiply<Expression> { Real { 1.0 }, Imaginary {} }.Generalize());
+            }
+            continue;
+        }
+        // n*i
+        if (auto img = RecursiveCast<Multiply<Expression, Imaginary>>(*addend); img != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Multiply<Expression, Imaginary>>(*vals[i]); valI != nullptr) {
+                    auto addResult = Add<Expression> { valI->GetMostSigOp(), img->GetMostSigOp() }.Accept(*this);
+                    if (!addResult) {
+                        return addResult;
+                    }
+                    vals[i] = Multiply<Expression> { *(addResult.value()), Imaginary {} }.Generalize();
+                    break;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(img->Generalize());
+            }
+            continue;
+        }
+        // single variable
+        if (auto var = RecursiveCast<Variable>(*addend); var != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Multiply<Expression, Variable>>(*vals[i]); valI != nullptr) {
+                    if (valI->GetLeastSigOp().GetName() == var->GetName()) {
+                        auto addResult = Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } }.Accept(*this);
+                        if (!addResult) {
+                            return addResult;
+                        }
+                        vals[i] = Multiply<Expression> { *(addResult.value()), *var }.Generalize();
+                        break;
+                    } else
+                        continue;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(Multiply<Expression> { Real { 1.0 }, *var }.Generalize());
+            }
+            continue;
+        }
+        // n*variable
+        if (auto var = RecursiveCast<Multiply<Expression, Variable>>(*addend); var != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Multiply<Expression, Variable>>(*vals[i]); valI != nullptr) {
+                    // if (auto zeroCase = RecursiveCast<Multiply<Real, Expression>>(*valI); zeroCase != nullptr) {}
+                    if (valI->GetLeastSigOp().GetName() == var->GetLeastSigOp().GetName()) {
+                        auto addResult = Add<Expression> { valI->GetMostSigOp(), var->GetMostSigOp() }.Accept(*this);
+                        if (!addResult) {
+                            return addResult;
+                        }
+                        vals[i] = Multiply<Expression> { *(addResult.value()), valI->GetLeastSigOp() }.Generalize();
+                        break;
+                    } else
+                        continue;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(var->Generalize());
+            }
+            continue;
+        }
+        // single exponent
+        if (auto exp = RecursiveCast<Exponent<Expression>>(*addend); exp != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(*vals[i]); valI != nullptr) {
+                    if (valI->GetLeastSigOp().Equals(*exp)) {
+                        auto addResult = Add<Expression> { valI->GetMostSigOp(), Real { 1.0 } }.Accept(*this);
+                        if (!addResult) {
+                            return addResult;
+                        }
+                        vals[i] = Multiply<Expression> { *(addResult.value()), *exp }.Generalize();
+                        break;
+                    } else
+                        continue;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(Multiply<Expression> { Real { 1.0 }, *exp }.Generalize());
+            }
+            continue;
+        }
+        // n*exponent
+        if (auto exp = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(*addend); exp != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(*vals[i]); valI != nullptr) {
+                    if (valI->GetLeastSigOp().Equals(exp->GetLeastSigOp())) {
+                        auto addResult = Add<Expression> { valI->GetMostSigOp(), exp->GetMostSigOp() }.Accept(*this);
+                        if (!addResult) {
+                            return addResult;
+                        }
+                        vals[i] = Multiply<Expression> { *(addResult.value()), valI->GetLeastSigOp() }.Generalize();
+                        break;
+                    } else
+                        continue;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(exp->Generalize());
+            }
+            continue;
+        }
+    }
+    // rebuild equation after simplification.
+
+    for (auto& val : vals) {
+        if (auto mul = RecursiveCast<Multiply<Real, Expression>>(*val); mul != nullptr) {
+            if (mul->GetMostSigOp().GetValue() == 1.0) {
+                val = mul->GetLeastSigOp().Generalize();
+            }
+        }
+    }
+
+    // filter out zero-equivalent expressions
+    std::vector<std::unique_ptr<Expression>> avals;
+    for (auto& val : vals) {
+        if (auto real = RecursiveCast<Real>(*val); real != nullptr) {
+            if (std::abs(real->GetValue()) <= EPSILON) {
+                continue;
+            }
+        }
+        if (auto mul = RecursiveCast<Multiply<Real, Expression>>(*val); mul != nullptr) {
+            if (std::abs(mul->GetMostSigOp().GetValue()) <= EPSILON) {
+                continue;
+            }
+        }
+        avals.push_back(val->Generalize());
+    }
+
+    if (auto vec = BuildFromVector<Add>(avals); vec != nullptr) {
+        return gsl::not_null { std::move(vec) };
+    }
+
+    return gsl::not_null { simplifiedAdd.Copy() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Subtract<>& subtract) -> RetT
+{
+    auto mostSigOp = subtract.GetMostSigOp().Copy();
+    auto leastSigOp = subtract.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return std::unexpected { simplifiedLeastSigOpResult.error() };
+    }
+
+    auto simplifiedMinuend = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedSubtrahend = std::move(simplifiedLeastSigOpResult).value();
+
+    const Subtract simplifiedSubtract { *simplifiedMinuend, *simplifiedSubtrahend };
+
+    // 2 - 1 = 1
+    if (auto realCase = RecursiveCast<Subtract<Real>>(simplifiedSubtract); realCase != nullptr) {
+        const Real& minuend = realCase->GetMostSigOp();
+        const Real& subtrahend = realCase->GetLeastSigOp();
+
+        return gsl::not_null { std::make_unique<Real>(minuend.GetValue() - subtrahend.GetValue()) };
+    }
+
+    // x - x = 0
+    if (simplifiedMinuend->Equals(*simplifiedSubtrahend)) {
+        return gsl::not_null { std::make_unique<Real>(Real { 0.0 }) };
+    }
+
+    if (auto matrixCase = RecursiveCast<Subtract<Matrix, Matrix>>(simplifiedSubtract); matrixCase != nullptr) {
+        const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp();
+        const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp();
+
+        if ((leftTerm.GetRows() == rightTerm.GetRows()) && (leftTerm.GetCols() == leftTerm.GetCols())) {
+            return gsl::not_null { std::make_unique<Matrix>(leftTerm.GetMatrix() - rightTerm.GetMatrix()) };
+        } else {
+            return gsl::not_null { std::make_unique<Subtract<Expression>>(leftTerm, rightTerm) };
+        }
+    }
+
+    // ax - x = (a-1)x
+    if (const auto minusOneCase = RecursiveCast<Subtract<Multiply<>, Expression>>(simplifiedSubtract); minusOneCase != nullptr) {
+        if (minusOneCase->GetMostSigOp().GetLeastSigOp().Equals(minusOneCase->GetLeastSigOp())) {
+            const Subtract newCoefficient { minusOneCase->GetMostSigOp().GetMostSigOp(), Real { 1.0 } };
+            return Multiply { newCoefficient, minusOneCase->GetLeastSigOp() }.Accept(*this);
+        }
+    }
+
+    // x-ax = (1-a)x
+    if (const auto oneMinusCase = RecursiveCast<Subtract<Expression, Multiply<>>>(simplifiedSubtract); oneMinusCase != nullptr) {
+        if (oneMinusCase->GetMostSigOp().Equals(oneMinusCase->GetLeastSigOp().GetLeastSigOp())) {
+            const Subtract newCoefficient { Real { 1.0 }, oneMinusCase->GetLeastSigOp().GetMostSigOp() };
+            return Multiply { newCoefficient, oneMinusCase->GetMostSigOp() }.Accept(*this);
+        }
+    }
+
+    // ax-bx= (a-b)x
+    if (const auto coefficientCase = RecursiveCast<Subtract<Multiply<>>>(simplifiedSubtract); coefficientCase != nullptr) {
+        if (coefficientCase->GetMostSigOp().GetLeastSigOp().Equals(coefficientCase->GetLeastSigOp().GetLeastSigOp())) {
+            const Subtract newCoefficient { coefficientCase->GetMostSigOp().GetMostSigOp(), coefficientCase->GetLeastSigOp().GetMostSigOp() };
+            return Multiply { newCoefficient, coefficientCase->GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+        }
+    }
+
+    // log(a) - log(b) = log(a / b)
+    if (const auto logCase = RecursiveCast<Subtract<Log<>>>(simplifiedSubtract); logCase != nullptr) {
+        if (logCase->GetMostSigOp().GetMostSigOp().Equals(logCase->GetLeastSigOp().GetMostSigOp())) {
+            const IExpression auto& base = logCase->GetMostSigOp().GetMostSigOp();
+            const IExpression auto& argument = Divide({ logCase->GetMostSigOp().GetLeastSigOp(), logCase->GetLeastSigOp().GetLeastSigOp() });
+            return gsl::not_null { std::make_unique<Log<>>(base, argument) };
+        }
+    }
+
+    // makes subtraction into addition because it is easier to deal with
+    auto negated = Multiply<Expression> { Real { -1 }, *simplifiedSubtrahend };
+    if (auto added = RecursiveCast<Add<Expression>>(negated.GetLeastSigOp()); added != nullptr) {
+        auto msOp = Multiply<Expression> { Real { -1.0 }, added->GetMostSigOp() }.Accept(*this);
+        if (!msOp) {
+            return msOp;
+        }
+        auto lsOp = Multiply<Expression> { Real { -1.0 }, added->GetLeastSigOp() }.Accept(*this);
+        if (!lsOp) {
+            return lsOp;
+        }
+
+        auto RHS = Add { *(msOp.value()), *(lsOp.value()) }.Accept(*this);
+        if (!RHS) {
+            return RHS;
+        }
+        return Add { *simplifiedMinuend, *(RHS.value()) }.Accept(*this);
+    }
+    if (auto subtracted = RecursiveCast<Subtract<Expression>>(negated.GetLeastSigOp()); subtracted != nullptr) {
+        auto msOp = Multiply<Expression> { Real { -1.0 }, subtracted->GetMostSigOp() }.Accept(*this);
+        if (!msOp) {
+            return msOp;
+        }
+        auto lsOp = subtracted->GetLeastSigOp().Accept(*this);
+        if (!lsOp) {
+            return lsOp;
+        }
+
+        auto RHS = Add { *(msOp.value()), *(lsOp.value()) }.Accept(*this);
+        if (!RHS) {
+            return RHS;
+        }
+        return Add { *simplifiedMinuend, *(RHS.value()) }.Accept(*this);
+    }
+    return Add { *simplifiedMinuend, negated }.Accept(*this);
+}
+
+auto SimplifyVisitor::TypedVisit(const Multiply<>& multiply) -> RetT
+{
+    auto mostSigOp = multiply.GetMostSigOp().Copy();
+    auto leastSigOp = multiply.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return std::unexpected { simplifiedLeastSigOpResult.error() };
+    }
+
+    auto simplifiedMultiplicand = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedMultiplier = std::move(simplifiedLeastSigOpResult).value();
+
+    Multiply simplifiedMultiply { *simplifiedMultiplicand, *simplifiedMultiplier };
+    if (auto onezerocase = RecursiveCast<Multiply<Real, Expression>>(simplifiedMultiply); onezerocase != nullptr) {
+        const Real& multiplicand = onezerocase->GetMostSigOp();
+        const Expression& multiplier = onezerocase->GetLeastSigOp();
+        if (std::abs(multiplicand.GetValue()) <= EPSILON) {
+            return gsl::not_null { std::make_unique<Real>(Real { 0.0 }) };
+        }
+        if (multiplicand.GetValue() == 1) {
+            return multiplier.Accept(*this);
+        }
+    }
+    if (auto realCase = RecursiveCast<Multiply<Real>>(simplifiedMultiply); realCase != nullptr) {
+        const Real& multiplicand = realCase->GetMostSigOp();
+        const Real& multiplier = realCase->GetLeastSigOp();
+        return gsl::not_null { std::make_unique<Real>(multiplicand.GetValue() * multiplier.GetValue()) };
+    }
+
+    if (auto ImgCase = RecursiveCast<Multiply<Imaginary>>(simplifiedMultiply); ImgCase != nullptr) {
+        return gsl::not_null { std::make_unique<Real>(-1.0) };
+    }
+
+    if (auto multCase = RecursiveCast<Multiply<Real, Divide<Expression>>>(simplifiedMultiply); multCase != nullptr) {
+        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }.Accept(*this);
+        if (!m) {
+            return m;
+        }
+        return gsl::not_null { Divide<Expression> { *(m.value()), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Generalize() };
+    }
+
+    if (auto exprCase = RecursiveCast<Multiply<Expression>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp())) {
+            return gsl::not_null { std::make_unique<Exponent<Expression, Expression>>(exprCase->GetMostSigOp(), Real { 2.0 }) };
+        }
+    }
+
+    if (auto rMatrixCase = RecursiveCast<Multiply<Real, Matrix>>(simplifiedMultiply); rMatrixCase != nullptr) {
+        return gsl::not_null { std::make_unique<Matrix>(rMatrixCase->GetLeastSigOp().GetMatrix() * rMatrixCase->GetMostSigOp().GetValue()) };
+    }
+
+    if (auto matrixCase = RecursiveCast<Multiply<Matrix, Matrix>>(simplifiedMultiply); matrixCase != nullptr) {
+        const Oasis::IExpression auto& leftTerm = matrixCase->GetMostSigOp();
+        const Oasis::IExpression auto& rightTerm = matrixCase->GetLeastSigOp();
+
+        if (leftTerm.GetCols() == rightTerm.GetRows()) {
+            return gsl::not_null { std::make_unique<Matrix>(leftTerm.GetMatrix() * rightTerm.GetMatrix()) };
+        } else {
+            // ERROR: INVALID DIMENSION
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(leftTerm, rightTerm) };
+        }
+    }
+
+    if (auto multCase = RecursiveCast<Multiply<Expression, Divide<Expression>>>(simplifiedMultiply); multCase != nullptr) {
+        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }.Accept(*this);
+        if (!m) {
+            return m;
+        }
+        return Divide<Expression> { *(m.value()), (multCase->GetLeastSigOp().GetLeastSigOp()) }.Accept(*this);
+    }
+
+    if (auto multCase = RecursiveCast<Multiply<Divide<Expression>, Divide<Expression>>>(simplifiedMultiply); multCase != nullptr) {
+        auto m = Multiply<Expression> { multCase->GetMostSigOp(), multCase->GetLeastSigOp().GetMostSigOp() }.Accept(*this);
+        if (!m) {
+            return m;
+        }
+        auto n = Multiply<Expression> { multCase->GetMostSigOp().GetLeastSigOp(), multCase->GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+        if (!n) {
+            return n;
+        }
+        return Divide<Expression> { *(m.value()), *(n.value()) }.Accept(*this);
+    }
+
+    if (auto exprCase = RecursiveCast<Multiply<Expression, Exponent<Expression, Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(), *(lsOp.value())) };
+        }
+    }
+
+    // x*x^n
+    if (auto exprCase = RecursiveCast<Multiply<Expression, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp(), *(lsOp.value())) };
+        }
+    }
+
+    if (auto exprCase = RecursiveCast<Multiply<Exponent<Expression>, Expression>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetLeastSigOp().Equals(exprCase->GetMostSigOp().GetMostSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Exponent<Expression>>(exprCase->GetLeastSigOp(), *(lsOp.value())) };
+        }
+    }
+
+    // x^n*x^m
+    if (auto exprCase = RecursiveCast<Multiply<Exponent<Expression>, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Exponent<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(), *(lsOp.value())) };
+        }
+    }
+
+    // a*x*x
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Expression>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp())) {
+            return gsl::not_null { std::make_unique<Multiply<Expression, Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
+                Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }) };
+        }
+    }
+
+    // a*x*b*x
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Multiply<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) {
+            auto msOp = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }.Accept(*this);
+            if (!msOp) {
+                return msOp;
+            }
+            auto lsOp = Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), Real { 2.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(*(msOp.value()), *(lsOp.value())) };
+        }
+    }
+
+    // a*x^n*x
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Expression>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
+                Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(), *(lsOp.value()) }) };
+        }
+    }
+
+    // a*x*x^n
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetMostSigOp(),
+                Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), *(lsOp.value()) }) };
+        }
+        if (exprCase->GetMostSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(exprCase->GetMostSigOp().GetLeastSigOp(),
+                Exponent<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), *(lsOp.value()) }) };
+        }
+    }
+
+    // a*x^n*b*x
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Multiply<Expression, Exponent<Expression>>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) {
+            auto msOp = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }.Accept(*this);
+            if (!msOp) {
+                return msOp;
+            }
+            auto lsLsOp = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsLsOp) {
+                return lsLsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(
+                *(msOp.value()), Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), *(lsLsOp.value()) }) };
+        }
+    }
+
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression>, Multiply<Exponent<Expression>, Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().Equals(exprCase->GetLeastSigOp().GetMostSigOp().GetMostSigOp())) {
+            auto msOp = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+            if (!msOp) {
+                return msOp;
+            }
+            auto lsLsOp = Add<Expression> { exprCase->GetLeastSigOp().GetMostSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsLsOp) {
+                return lsLsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(
+                *(msOp.value()), Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp(), *(lsLsOp.value()) }) };
+        }
+    }
+
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Multiply<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) {
+            auto msOp = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+            if (!msOp) {
+                return msOp;
+            }
+            auto lsLsOp = Add<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+            if (!lsLsOp) {
+                return lsLsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(
+                *(msOp.value()),
+                Exponent<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp(), *(lsLsOp.value()) }) };
+        }
+    }
+
+    // a*x^n*x^m
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Exponent<Expression>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp())) {
+            auto lsOp = Add<Expression> { exprCase->GetLeastSigOp().GetMostSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(
+                exprCase->GetMostSigOp().GetMostSigOp(),
+                Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
+                    *(lsOp.value()) }) };
+        }
+    }
+
+    // a*x^n*b*x^m
+    if (auto exprCase = RecursiveCast<Multiply<Multiply<Expression, Exponent<Expression>>, Multiply<Expression, Exponent<Expression>>>>(simplifiedMultiply); exprCase != nullptr) {
+        if (exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp().Equals(exprCase->GetLeastSigOp().GetLeastSigOp().GetMostSigOp())) {
+            auto msOp = Multiply<Expression> { exprCase->GetMostSigOp().GetMostSigOp(), exprCase->GetLeastSigOp().GetMostSigOp() }.Accept(*this);
+            if (!msOp) {
+                return msOp;
+            }
+            auto lsOp = Add<Expression> { exprCase->GetLeastSigOp().GetLeastSigOp().GetLeastSigOp(), exprCase->GetMostSigOp().GetLeastSigOp().GetLeastSigOp() }.Accept(*this);
+            if (!lsOp) {
+                return lsOp;
+            }
+            return gsl::not_null { std::make_unique<Multiply<Expression>>(
+                *(msOp.value()),
+                Exponent<Expression> { exprCase->GetMostSigOp().GetLeastSigOp().GetMostSigOp(),
+                    *(lsOp.value()) }) };
+        }
+    }
+
+    // IDK why this is commented out but it was commented out where I sourced it from
+    //    if (auto negate = RecursiveCast<Multiply<Real, Negate<Subtract<Expression>>>>(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<std::unique_ptr<Expression>> multiplies;
+    std::vector<std::unique_ptr<Expression>> vals;
+    simplifiedMultiply.Flatten(multiplies);
+    for (const auto& multiplicand : multiplies) {
+        size_t i = 0;
+        if (auto real = RecursiveCast<Real>(*multiplicand); real != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Real>(*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<Imaginary>(*multiplicand); img != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Exponent<Imaginary, Expression>>(*vals[i]); valI != nullptr) {
+                    auto lsOp = Add<Expression> { valI->GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+                    if (!lsOp) {
+                        return lsOp;
+                    }
+                    vals[i] = Exponent<Expression> { Imaginary {}, *(lsOp.value()) }.Generalize();
+                    break;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(Exponent<Expression> { Imaginary {}, Real { 1.0 } }.Generalize());
+            }
+            continue;
+        }
+        // i^n
+        if (auto img = RecursiveCast<Exponent<Imaginary, Expression>>(*multiplicand); img != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Exponent<Imaginary, Expression>>(*vals[i]); valI != nullptr) {
+                    auto lsOp = Add<Expression> { valI->GetLeastSigOp(), img->GetLeastSigOp() }.Accept(*this);
+                    if (!lsOp) {
+                        return lsOp;
+                    }
+                    vals[i] = Exponent<Expression> { Imaginary {}, *(lsOp.value()) }.Generalize();
+                    break;
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                // vals.push_back(Multiply<Expression> { img->GetMostSigOp(), Imaginary {} }.Generalize());
+                vals.push_back(img->Generalize());
+            }
+            continue;
+        }
+        // expr^n
+        if (auto expr = RecursiveCast<Exponent<Expression, Expression>>(*multiplicand); expr != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Exponent<Expression, Expression>>(*vals[i]); valI != nullptr) {
+                    if (valI->GetMostSigOp().Equals(expr->GetMostSigOp())) {
+                        auto lsOp = Add<Expression> { valI->GetLeastSigOp(), expr->GetLeastSigOp() }.Accept(*this);
+                        if (!lsOp) {
+                            return lsOp;
+                        }
+                        vals[i] = Exponent<Expression> { valI->GetMostSigOp(), *(lsOp.value()) }.Generalize();
+                        break;
+                    }
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                // vals.push_back(Multiply<Expression> { img->GetMostSigOp(), Imaginary {} }.Generalize());
+                vals.push_back(expr->Generalize());
+            }
+            continue;
+        }
+        // single expr
+        if (auto expr = RecursiveCast<Expression>(*multiplicand); expr != nullptr) {
+            for (; i < vals.size(); i++) {
+                if (auto valI = RecursiveCast<Exponent<Expression, Expression>>(*vals[i]); valI != nullptr) {
+                    auto lsOp = Add<Expression> { valI->GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+                    if (!lsOp) {
+                        return lsOp;
+                    }
+                    if (valI->GetMostSigOp().Equals(*expr)) {
+                        vals[i] = Exponent<Expression> { valI->GetMostSigOp(), *(lsOp.value()) }.Generalize();
+                        break;
+                    }
+                }
+            }
+            if (i >= vals.size()) {
+                // check to make sure it is one thing only
+                vals.push_back(Exponent<Expression> { *expr, Real { 1.0 } }.Generalize());
+            }
+            continue;
+        }
+    }
+
+    // makes all expr^1 into expr
+    for (auto& val : vals) {
+        if (auto exp = RecursiveCast<Exponent<Expression, Real>>(*val); exp != nullptr) {
+            if (exp->GetLeastSigOp().GetValue() == 1.0) {
+                val = exp->GetMostSigOp().Generalize();
+            }
+        }
+        if (auto mul = RecursiveCast<Multiply<Real, Expression>>(*val); mul != nullptr) {
+            if (mul->GetMostSigOp().GetValue() == 1.0) {
+                val = mul->GetLeastSigOp().Generalize();
+            }
+        }
+    }
+
+    return gsl::not_null { BuildFromVector<Multiply>(vals) };
+}
+
+auto SimplifyVisitor::TypedVisit(const Divide<>& divide) -> RetT
+{
+    auto mostSigOp = divide.GetMostSigOp().Copy();
+    auto leastSigOp = divide.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return std::unexpected { simplifiedLeastSigOpResult.error() };
+    }
+
+    auto simplifiedDividend = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedDivider = std::move(simplifiedLeastSigOpResult).value();
+
+    Divide simplifiedDivide { *simplifiedDividend, *simplifiedDivider };
+
+    if (auto realCase = RecursiveCast<Divide<Real>>(simplifiedDivide); realCase != nullptr) {
+        const Real& dividend = realCase->GetMostSigOp();
+        const Real& divisor = realCase->GetLeastSigOp();
+        return gsl::not_null { std::make_unique<Real>(dividend.GetValue() / divisor.GetValue()) };
+    }
+
+    // log(a)/log(b)=log[b](a)
+    if (auto logCase = RecursiveCast<Divide<Log<Expression, Expression>, Log<Expression, Expression>>>(simplifiedDivide); logCase != nullptr) {
+        if (logCase->GetMostSigOp().GetMostSigOp().Equals(logCase->GetLeastSigOp().GetMostSigOp())) {
+            const IExpression auto& base = logCase->GetLeastSigOp().GetLeastSigOp();
+            const IExpression auto& argument = logCase->GetMostSigOp().GetLeastSigOp();
+            return gsl::not_null { std::make_unique<Log<Expression>>(base, argument) };
+        }
+    }
+    // convert the terms in numerator and denominator into a vector to make manipulations easier
+    std::vector<std::unique_ptr<Expression>> numerator;
+    std::vector<std::unique_ptr<Expression>> denominator;
+    std::vector<std::unique_ptr<Expression>> result;
+    std::vector<std::unique_ptr<Expression>> numeratorVals;
+    std::vector<std::unique_ptr<Expression>> denominatorVals;
+
+    if (auto multipliedNumerator = RecursiveCast<Multiply<Expression>>(*simplifiedDividend); multipliedNumerator != nullptr) {
+        multipliedNumerator->Flatten(numerator);
+    } else {
+        numerator.push_back(simplifiedDividend->Copy());
+    }
+
+    if (auto multipliedDenominator = RecursiveCast<Multiply<Expression>>(*simplifiedDivider); multipliedDenominator != nullptr) {
+        multipliedDenominator->Flatten(denominator);
+    } else {
+        denominator.push_back(simplifiedDivider->Copy());
+    }
+
+    // now that we have the terms in a vector, we have to cancel like terms and simplify them
+    result.reserve(numerator.size());
+    for (const auto& num : numerator) {
+        if (auto div = RecursiveCast<Divide<Expression>>(*num); div != nullptr) {
+            std::vector<std::unique_ptr<Expression>> flat_n;
+            std::vector<std::unique_ptr<Expression>> flat_d;
+            if (auto mN = RecursiveCast<Multiply<Expression>>(*div); mN != nullptr) {
+                mN->Flatten(flat_n);
+            } else {
+                result.push_back(div->GetMostSigOp().Copy());
+            }
+
+            if (auto mD = RecursiveCast<Multiply<Expression>>(*div); mD != nullptr) {
+                mD->Flatten(flat_d);
+            } else {
+                denominator.push_back(div->GetLeastSigOp().Copy());
+            }
+
+            for (const auto& v : flat_n) {
+                result.push_back(v->Copy());
+            }
+            for (const auto& v : flat_d) {
+                denominator.push_back(v->Copy());
+            }
+        } else {
+            result.push_back(num->Copy());
+        }
+    }
+
+    for (const auto& denom : denominator) {
+        size_t i = 0;
+        if (auto real = RecursiveCast<Real>(*denom); real != nullptr) {
+            for (; i < result.size(); i++) {
+                if (auto resI = RecursiveCast<Real>(*result[i]); resI != nullptr) {
+                    result[i] = Real { resI->GetValue() / real->GetValue() }.Generalize();
+                    break;
+                }
+            }
+            if (i >= result.size()) {
+                result.push_back(Real { 1 / real->GetValue() }.Generalize());
+            }
+            continue;
+        }
+        if (auto img = RecursiveCast<Imaginary>(*denom); img != nullptr) {
+            for (; i < result.size(); i++) {
+                if (auto resI = RecursiveCast<Imaginary>(*result[i]); resI != nullptr) {
+                    result[i] = Real { 1.0 }.Generalize();
+                    break;
+                }
+            }
+            if (i >= result.size()) {
+                result.push_back(Exponent<Expression> { Imaginary {}, Real { -1.0 } }.Generalize());
+            }
+            continue;
+        }
+        if (auto var = RecursiveCast<Variable>(*denom); var != nullptr) {
+            for (; i < result.size(); i++) {
+                if (auto resIexp = RecursiveCast<Exponent<Variable, Expression>>(*result[i]); resIexp != nullptr) {
+                    if (resIexp->GetMostSigOp().Equals(*var)) {
+                        auto lsOp = Subtract<Expression> { resIexp->GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+                        if (!lsOp) {
+                            return lsOp;
+                        }
+                        result[i] = Exponent<Expression> { *var, *(lsOp.value()) }.Generalize();
+                        break;
+                    }
+                } else if (auto resI = RecursiveCast<Variable>(*result[i]); resI != nullptr) {
+                    if (resI->Equals(*var)) {
+                        result[i] = Real { 1.0 }.Generalize();
+                        break;
+                    }
+                }
+            }
+            if (i >= result.size()) {
+                result.push_back(Exponent<Expression> { *var, Real { -1.0 } }.Generalize());
+            }
+            continue;
+        }
+        if (auto var = RecursiveCast<Exponent<Variable, Expression>>(*denom); var != nullptr) {
+            for (; i < result.size(); i++) {
+                if (auto resIexp = RecursiveCast<Exponent<Variable, Expression>>(*result[i]); resIexp != nullptr) {
+                    if (resIexp->GetMostSigOp().Equals(var->GetMostSigOp())) {
+                        auto lsOp = Subtract<Expression> { resIexp->GetLeastSigOp(), var->GetLeastSigOp() }.Accept(*this);
+                        if (!lsOp) {
+                            return lsOp;
+                        }
+                        result[i] = Exponent<Expression> { var->GetMostSigOp(), *(lsOp.value()) }.Generalize();
+                        break;
+                    }
+                } else if (auto resI = RecursiveCast<Variable>(*result[i]); resI != nullptr) {
+                    if (resI->Equals(*var)) {
+                        auto lsOp = Subtract<Expression> { Real { 1.0 }, var->GetLeastSigOp() }.Accept(*this);
+                        if (!lsOp) {
+                            return lsOp;
+                        }
+                        result[i] = Exponent<Expression> { var->GetMostSigOp(), *(lsOp.value()) }.Generalize();
+                    }
+                }
+            }
+            if (i >= result.size()) {
+                result.push_back(Exponent<Expression> { var->GetMostSigOp(), Multiply { Real { -1.0 }, var->GetLeastSigOp() } }.Generalize());
+            }
+            continue;
+        }
+        if (auto expExpr = RecursiveCast<Exponent<Expression>>(*denom); expExpr != nullptr) {
+            for (; i < result.size(); i++) {
+                if (auto resExpr = RecursiveCast<Exponent<Expression, Expression>>(*result[i]); resExpr != nullptr) {
+                    if (expExpr->GetMostSigOp().Equals(resExpr->GetMostSigOp())) {
+                        auto lsOp = Subtract { resExpr->GetLeastSigOp(), expExpr->GetLeastSigOp() }.Accept(*this);
+                        if (!lsOp) {
+                            return lsOp;
+                        }
+                        auto resultVisited = Exponent { expExpr->GetMostSigOp(), *(lsOp.value()) }.Accept(*this);
+                        if (!resultVisited) {
+                            return resultVisited;
+                        }
+                        result[i] = std::move(resultVisited.value());
+                        break;
+                    }
+                } else if (result[i]->Equals(expExpr->GetMostSigOp())) {
+                    auto lsOp = Subtract { Real { 1.0 }, expExpr->GetLeastSigOp() }.Accept(*this);
+                    if (!lsOp) {
+                        return lsOp;
+                    }
+                    auto resultVisited = Exponent { expExpr->GetMostSigOp(), *(lsOp.value()) }.Accept(*this);
+                    if (!resultVisited) {
+                        return resultVisited;
+                    }
+                    result[i] = std::move(resultVisited.value());
+                    break;
+                }
+            }
+            if (i >= result.size()) {
+                auto resultVisited = Exponent<Expression> { expExpr->GetMostSigOp(), Multiply { Real { -1.0 }, expExpr->GetLeastSigOp() } }.Accept(*this);
+                if (!resultVisited) {
+                    return resultVisited;
+                }
+                result.push_back(std::move(resultVisited.value()));
+            }
+            continue;
+        }
+        for (; i < result.size(); i++) {
+            if (auto resExpr = RecursiveCast<Exponent<Expression, Expression>>(*result[i]); resExpr != nullptr) {
+                if (denom->Equals(resExpr->GetMostSigOp())) {
+                    auto lsOp = Subtract { resExpr->GetLeastSigOp(), Real { 1.0 } }.Accept(*this);
+                    if (!lsOp) {
+                        return lsOp;
+                    }
+                    auto resultVisited = Exponent { *denom, *(lsOp.value()) }.Accept(*this);
+                    if (!resultVisited) {
+                        return resultVisited;
+                    }
+                    result[i] = std::move(resultVisited.value());
+                    break;
+                }
+            } else if (result[i]->Equals(*denom)) {
+                result[i] = Real { 1.0 }.Generalize();
+                break;
+            }
+        }
+        if (i >= result.size()) {
+            result.push_back(Exponent<Expression> { *denom, Real { -1.0 } }.Generalize());
+        }
+    }
+
+    // rebuild into tree
+    for (const auto& val : result) {
+        if (auto valI = RecursiveCast<Real>(*val); valI != nullptr) {
+            if (valI->GetValue() != 1.0) {
+                numeratorVals.push_back(valI->Generalize());
+            }
+        } else if (auto var = RecursiveCast<Variable>(*val); var != nullptr) {
+            numeratorVals.push_back(val->Generalize());
+        } else if (auto img = RecursiveCast<Imaginary>(*val); img != nullptr) {
+            numeratorVals.push_back(val->Generalize());
+        } else if (auto expR = RecursiveCast<Exponent<Expression, Real>>(*val); expR != nullptr) {
+            if (expR->GetLeastSigOp().GetValue() < 0.0) {
+                denominatorVals.push_back(Exponent { expR->GetMostSigOp(), Real { expR->GetLeastSigOp().GetValue() * -1.0 } }.Generalize());
+            } else {
+                numeratorVals.push_back(val->Generalize());
+            }
+        } else if (auto exp = RecursiveCast<Exponent<Expression, Multiply<Real, Expression>>>(*val); exp != nullptr) {
+            if (exp->GetLeastSigOp().GetMostSigOp().GetValue() < 0.0) {
+                denominatorVals.push_back(Exponent { exp->GetMostSigOp(), exp->GetLeastSigOp().GetLeastSigOp() }.Generalize());
+            } else {
+                numeratorVals.push_back(val->Generalize());
+            }
+        } else {
+            numeratorVals.push_back(val->Generalize());
+        }
+    }
+
+    // makes expr^1 into expr
+    for (auto& val : numeratorVals) {
+        if (auto exp = RecursiveCast<Exponent<Expression, Real>>(*val); exp != nullptr) {
+            if (exp->GetLeastSigOp().GetValue() == 1.0) {
+                val = exp->GetMostSigOp().Generalize();
+            }
+        }
+    }
+
+    for (auto& val : denominatorVals) {
+        if (auto exp = RecursiveCast<Exponent<Expression, Real>>(*val); exp != nullptr) {
+            if (exp->GetLeastSigOp().GetValue() == 1.0) {
+                val = exp->GetMostSigOp().Generalize();
+            }
+        }
+    }
+
+    auto dividend = numeratorVals.size() == 1 ? std::move(numeratorVals.front()) : BuildFromVector<Multiply>(numeratorVals);
+    auto divisor = denominatorVals.size() == 1 ? std::move(denominatorVals.front()) : BuildFromVector<Multiply>(denominatorVals);
+
+    // rebuild subtrees
+    if (!dividend && divisor)
+        return gsl::not_null { Divide { Real { 1.0 }, *divisor }.Copy() };
+
+    if (dividend && !divisor)
+        return gsl::not_null { dividend->Copy() };
+
+    if (!dividend && !divisor)
+        return gsl::not_null { Real { 1.0 }.Copy() };
+
+    return gsl::not_null { Divide { *dividend, *divisor }.Copy() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Exponent<>& exponent) -> RetT
+{
+    static auto match_cast = MatchCast<Expression>()
+                                 .Case(
+                                     [](const Exponent<Expression, Real>& zeroCase) -> bool {
+                                         // TODO: Optimize. We're calling RecursiveCast<Expression<Real, Expression> twice. Perhaps a map of expression types to a vector of lambdas?
+                                         const Real& power = zeroCase.GetLeastSigOp();
+                                         return power.GetValue() == 0.0;
+                                     },
+                                     [](const Exponent<Expression, Real>&, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         return gsl::make_not_null(std::make_unique<Real>(1.0));
+                                     })
+                                 .Case(
+                                     [](const Exponent<Real, Expression>& zeroCase) -> bool {
+                                         const Real& base = zeroCase.GetMostSigOp();
+                                         return base.GetValue() == 0.0;
+                                     },
+                                     [](const Exponent<Real, Expression>&, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         return gsl::make_not_null(std::make_unique<Real>(0.0));
+                                     })
+                                 .Case(
+                                     [](const Exponent<Real>&) -> bool { return true; },
+                                     [](const Exponent<Real>& realCase, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         const Real& base = realCase.GetMostSigOp();
+                                         const Real& power = realCase.GetLeastSigOp();
+                                         return gsl::make_not_null(std::make_unique<Real>(pow(base.GetValue(), power.GetValue())));
+                                     })
+                                 .Case(
+                                     [](const Exponent<Expression, Real>& oneCase) -> bool {
+                                         const Real& power = oneCase.GetLeastSigOp();
+                                         return power.GetValue() == 1.0;
+                                     },
+                                     [](const Exponent<Expression, Real>& oneCase, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         return gsl::make_not_null(oneCase.GetMostSigOp().Copy());
+                                     })
+                                 .Case(
+                                     [](const Exponent<Real, Expression>& oneCase) -> bool {
+                                         const Real& base = oneCase.GetMostSigOp();
+                                         return base.GetValue() == 1.0;
+                                     },
+                                     [](const Exponent<Real, Expression>&, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         return gsl::make_not_null(std::make_unique<Real>(1.0));
+                                     })
+                                 .Case(
+                                     [](const Exponent<Imaginary, Real>&) -> bool { return true; },
+                                     [](const Exponent<Imaginary, Real>& ImgCase, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         switch (const auto power = std::fmod((ImgCase.GetLeastSigOp()).GetValue(), 4); static_cast<int>(power)) {
+                                         case 0:
+                                             return gsl::make_not_null(std::make_unique<Real>(1));
+                                         case 1:
+                                             return gsl::make_not_null(std::make_unique<Imaginary>());
+                                         case 2:
+                                             return gsl::make_not_null(std::make_unique<Real>(-1));
+                                         case 3:
+                                             return gsl::make_not_null(std::make_unique<Negate<Imaginary>>(Imaginary {}));
+                                         default:
+                                             return std::unexpected { "std::fmod returned an invalid value" };
+                                         }
+                                     })
+                                 .Case(
+                                     [](const Exponent<Multiply<Real, Expression>, Real>& ImgCase) -> bool {
+                                         return ImgCase.GetMostSigOp().GetMostSigOp().GetValue() < 0 && ImgCase.GetLeastSigOp().GetValue() == 0.5;
+                                     },
+                                     [](const Exponent<Multiply<Real, Expression>, Real>& ImgCase, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         Multiply mul {
+                                             Multiply { Real { pow(std::abs(ImgCase.GetMostSigOp().GetMostSigOp().GetValue()), 0.5) },
+                                                 Exponent { ImgCase.GetMostSigOp().GetLeastSigOp(), Real { 0.5 } } },
+                                             Imaginary {}
+                                         };
+
+                                         return gsl::make_not_null(mul.Copy());
+                                     })
+                                 .Case(
+                                     [](const Exponent<Exponent<>, Expression>&) -> bool { return true; },
+                                     [](const Exponent<Exponent<>, Expression>& expExpCase, SimplifyVisitor* visitor) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         auto lsOp = Multiply { expExpCase.GetMostSigOp().GetLeastSigOp(), expExpCase.GetLeastSigOp() }.Accept(*visitor);
+                                         if (!lsOp) {
+                                             return lsOp;
+                                         }
+                                         Exponent exp { expExpCase.GetMostSigOp().GetMostSigOp(), *(lsOp.value()) };
+                                         return gsl::make_not_null(exp.Copy());
+                                     })
+                                 .Case(
+                                     [](const Exponent<Expression, Log<>>& logCase) -> bool {
+                                         return logCase.GetMostSigOp().Equals(logCase.GetLeastSigOp().GetMostSigOp());
+                                     },
+                                     [](const Exponent<Expression, Log<>>& logCase, const void*) -> std::expected<gsl::not_null<std::unique_ptr<Expression>>, std::string_view> {
+                                         return gsl::make_not_null(logCase.GetLeastSigOp().GetLeastSigOp().Copy());
+                                     });
+
+    auto mostSigOp = exponent.GetMostSigOp().Copy();
+    auto leastSigOp = exponent.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return std::unexpected { simplifiedLeastSigOpResult.error() };
+    }
+
+    auto simplifiedBase = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedPower = std::move(simplifiedLeastSigOpResult).value();
+
+    const Exponent simplifiedExponent { *simplifiedBase, *simplifiedPower };
+    auto matchResult = match_cast.Execute(simplifiedExponent, this).value();
+    return gsl::not_null { matchResult ? std::move(matchResult) : std::move(simplifiedExponent.Copy()) };
+}
+
+auto SimplifyVisitor::TypedVisit(const Log<>& logIn) -> RetT
+{
+    auto mostSigOp = logIn.GetMostSigOp().Copy();
+    auto leastSigOp = logIn.GetLeastSigOp().Copy();
+
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return std::unexpected { simplifiedLeastSigOpResult.error() };
+    }
+
+    auto simplifiedBase = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedPower = std::move(simplifiedLeastSigOpResult).value();
+
+    const Log simplifiedLog { *simplifiedBase, *simplifiedPower };
+
+    if (const auto realBaseCase = RecursiveCast<Log<Real, Expression>>(simplifiedLog); realBaseCase != nullptr) {
+        if (const Real& b = realBaseCase->GetMostSigOp(); b.GetValue() <= 0.0 || b.GetValue() == 1) {
+            return gsl::not_null { std::make_unique<Undefined>() };
+        }
+    }
+
+    if (const auto realExponentCase = RecursiveCast<Log<Expression, Real>>(simplifiedLog); realExponentCase != nullptr) {
+        const Real& argument = realExponentCase->GetLeastSigOp();
+
+        if (argument.GetValue() <= 0.0) {
+            return gsl::not_null { std::make_unique<Undefined>() };
+        }
+
+        if (argument.GetValue() == 1.0) {
+            return gsl::not_null { std::make_unique<Real>(0.0) };
+        }
+    }
+
+    if (const auto realCase = RecursiveCast<Log<Real>>(simplifiedLog); realCase != nullptr) {
+        const Real& base = realCase->GetMostSigOp();
+        const Real& argument = realCase->GetLeastSigOp();
+
+        return gsl::not_null { std::make_unique<Real>(log2(argument.GetValue()) * (1 / log2(base.GetValue()))) };
+    }
+
+    // log(a) with a < 0 log(-a)
+    if (auto negCase = RecursiveCast<Log<Expression, Real>>(simplifiedLog); negCase != nullptr) {
+        if (negCase->GetLeastSigOp().GetValue() < 0) {
+            return gsl::not_null { Add<Expression> { Log { negCase->GetMostSigOp(), Real { -1 * negCase->GetLeastSigOp().GetValue() } }, Multiply<Expression> { Imaginary {}, Pi {} } }.Generalize() };
+        }
+    }
+
+    // log[a](a) = 1
+    if (const auto sameCase = RecursiveCast<Log<Expression, Expression>>(simplifiedLog); sameCase != nullptr) {
+        if (sameCase->leastSigOp->Equals(*sameCase->mostSigOp)) {
+            return gsl::not_null { Real { 1 }.Generalize() };
+        }
+    }
+
+    // log[a](b^x) = x * log[a](b)
+    if (const auto expCase = RecursiveCast<Log<Expression, Exponent<>>>(simplifiedLog); expCase != nullptr) {
+        const auto exponent = expCase->GetLeastSigOp();
+        const IExpression auto& log = Log<Expression>(expCase->GetMostSigOp(), exponent.GetMostSigOp()); // might need to check that it isnt nullptr
+        const IExpression auto& factor = exponent.GetLeastSigOp();
+        return Oasis::Multiply<Oasis::Expression>(factor, log).Accept(*this);
+    }
+
+    return gsl::not_null { simplifiedLog.Copy() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Negate<Expression>& negate) -> RetT
+{
+    auto op = negate.GetOperand().Copy();
+    if (!op) {
+        return std::unexpected { "Missing operand." };
+    }
+
+    auto simplifiedMostSigOpResult = op->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+
+    auto simplifiedOp = std::move(simplifiedMostSigOpResult).value();
+
+    return std::move(Multiply { Real { -1 }, *simplifiedOp }.Accept(*this)).value();
+}
+
+auto SimplifyVisitor::TypedVisit(const Sine<Expression>& sine) -> RetT
+{
+    auto op = sine.GetOperand().Copy();
+    if (!op) {
+        return std::unexpected { "Missing operand." };
+    }
+
+    auto simplifiedMostSigOpResult = op->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return std::unexpected { simplifiedMostSigOpResult.error() };
+    }
+
+    auto simplifiedOp = std::move(simplifiedMostSigOpResult).value();
+
+    return simplifiedOp->Accept(*this);
+}
+
+auto SimplifyVisitor::TypedVisit(const Derivative<>& derivative) -> RetT
+{
+    auto mostSigOp = derivative.GetMostSigOp().Copy();
+    auto leastSigOp = derivative.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return simplifiedMostSigOpResult;
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return simplifiedLeastSigOpResult;
+    }
+
+    auto simplifiedExpression = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedVar = std::move(simplifiedLeastSigOpResult).value();
+    auto simplifiedDiff = simplifiedExpression->Differentiate(*simplifiedVar);
+    auto s = simplifiedDiff->Accept(*this);
+    if (!s) {
+        return s;
+    }
+    return gsl::not_null { std::move(s.value())->Accept(*this).value() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Integral<>& integral) -> RetT
+{
+    auto mostSigOp = integral.GetMostSigOp().Copy();
+    auto leastSigOp = integral.GetLeastSigOp().Copy();
+    if (!mostSigOp) {
+        return std::unexpected { "Missing most significant operand." };
+    }
+    if (!leastSigOp) {
+        return std::unexpected { "Missing least significant operand." };
+    }
+
+    auto simplifiedMostSigOpResult = mostSigOp->Accept(*this);
+    auto simplifiedLeastSigOpResult = leastSigOp->Accept(*this);
+
+    if (!simplifiedMostSigOpResult) {
+        return simplifiedMostSigOpResult;
+    }
+    if (!simplifiedLeastSigOpResult) {
+        return simplifiedLeastSigOpResult;
+    }
+
+    auto simplifiedIntegrand = std::move(simplifiedMostSigOpResult).value();
+    auto simplifiedDifferential = std::move(simplifiedLeastSigOpResult).value();
+
+    auto integrated = simplifiedIntegrand->Integrate(*simplifiedDifferential);
+    return gsl::not_null { std::move(integrated) };
+}
+
+auto SimplifyVisitor::TypedVisit(const Matrix& matrix) -> RetT
+{
+    return gsl::not_null { matrix.Copy() };
+}
+
+auto SimplifyVisitor::TypedVisit(const EulerNumber&) -> RetT
+{
+    return gsl::not_null { EulerNumber {}.Copy() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Pi&) -> RetT
+{
+    return gsl::not_null { Pi {}.Copy() };
+}
+
+auto SimplifyVisitor::TypedVisit(const Magnitude<Expression>& magnitude) -> RetT
+{
+    auto simplified = magnitude.GetOperand().Accept(*this);
+    if (!simplified) {
+        return simplified;
+    }
+    auto simpOp = std::move(simplified.value());
+    if (auto realCase = RecursiveCast<Real>(*simpOp); realCase != nullptr) {
+        double val = realCase->GetValue();
+        return gsl::not_null { val >= 0.0 ? std::make_unique<Real>(val) : std::make_unique<Real>(-val) };
+    }
+    if (auto imgCase = RecursiveCast<Imaginary>(*simpOp); imgCase != nullptr) {
+        return gsl::not_null { std::make_unique<Real>(1.0) };
+    }
+    if (auto mulImgCase = RecursiveCast<Multiply<Expression, Imaginary>>(*simpOp); mulImgCase != nullptr) {
+        return Magnitude<Expression> { mulImgCase->GetMostSigOp() }.Accept(*this);
+    }
+    if (auto addCase = RecursiveCast<Add<Expression, Imaginary>>(*simpOp); addCase != nullptr) {
+        return Exponent { Add<Expression> { Exponent<Expression> { addCase->GetMostSigOp(), Real { 2 } },
+                              Real { 1.0 } },
+            Real { 0.5 } }
+            .Accept(*this);
+    }
+    if (auto addCase = RecursiveCast<Add<Expression, Multiply<Expression, Imaginary>>>(*simpOp); addCase != nullptr) {
+        return Exponent { Add<Expression> { Exponent<Expression> { addCase->GetMostSigOp(), Real { 2 } },
+                              Exponent<Expression> { addCase->GetLeastSigOp().GetMostSigOp(), Real { 2 } } },
+            Real { 0.5 } }
+            .Accept(*this);
+    }
+    if (auto matrixCase = RecursiveCast<Matrix>(*simpOp); matrixCase != nullptr) {
+        double sum = 0;
+        for (size_t i = 0; i < matrixCase->GetRows(); i++) {
+            for (size_t j = 0; j < matrixCase->GetCols(); j++) {
+                sum += pow(matrixCase->GetMatrix()(i, j), 2);
+            }
+        }
+        return Exponent { Real { sum }, Real { 0.5 } }.Accept(*this);
+    }
+
+    return gsl::not_null { magnitude.Generalize() };
+}
+
+} // Oasis
\ No newline at end of file
diff --git a/src/Sine.cpp b/src/Sine.cpp
new file mode 100644
index 00000000..fcf41877
--- /dev/null
+++ b/src/Sine.cpp
@@ -0,0 +1,24 @@
+//
+// Created by Andrew Nazareth on 9/19/25.
+//
+
+#include "Oasis/Sine.hpp"
+#include "Oasis/Expression.hpp"
+#include "Oasis/UnaryExpression.hpp"
+
+namespace Oasis {
+auto Sine<Expression>::Simplify() const -> std::unique_ptr<Expression>
+{
+    return nullptr;
+}
+
+auto Sine<Expression>::Differentiate(const Expression&) const -> std::unique_ptr<Expression>
+{
+    return nullptr;
+}
+
+auto Sine<Expression>::Integrate(const Expression&) const -> std::unique_ptr<Expression>
+{
+    return nullptr;
+}
+} // Oasis
\ No newline at end of file
diff --git a/src/Subtract.cpp b/src/Subtract.cpp
index 81c9893d..df0a1e74 100644
--- a/src/Subtract.cpp
+++ b/src/Subtract.cpp
@@ -82,29 +82,50 @@ auto Subtract<Expression>::Simplify() const -> std::unique_ptr<Expression>
             return std::make_unique<Log<>>(base, argument);
         }
     }
-
+    SimplifyVisitor simplifyVisitor {};
     // makes subtraction into addition because it is easier to deal with
     auto negated = Multiply<Expression> { Real { -1 }, *simplifiedSubtrahend };
     if (auto added = RecursiveCast<Add<Expression>>(negated.GetLeastSigOp()); added != nullptr) {
-        auto RHS = Add { *(Multiply<Expression> { Real { -1.0 }, added->GetMostSigOp() }.Simplify()),
-            *(Multiply<Expression> { Real { -1.0 }, added->GetLeastSigOp() }.Simplify()) }
-                       .Simplify();
-        return Add { *simplifiedMinuend, *RHS }.Simplify();
+        auto m = Multiply<Expression> { Real { -1.0 }, added->GetMostSigOp() };
+        auto n = Multiply<Expression> { Real { -1.0 }, added->GetLeastSigOp() };
+        auto ms = m.Accept(simplifyVisitor);
+        auto ns = n.Accept(simplifyVisitor);
+        if (!ms || !ns) {
+            return Add { *simplifiedMinuend, Add { m, n } }.Generalize();
+        }
+        auto RHS = Add { *(std::move(ms.value())), *(std::move(ns.value())) };
+        auto s = RHS.Accept(simplifyVisitor);
+        if (!s) {
+            return Add { *simplifiedMinuend, RHS }.Generalize();
+        }
+        return Add { *simplifiedMinuend, *(std::move(s.value())) }.Accept(simplifyVisitor).value();
     }
     if (auto subtracted = RecursiveCast<Subtract<Expression>>(negated.GetLeastSigOp()); subtracted != nullptr) {
-        auto RHS = Add { *(Multiply<Expression> { Real { -1.0 }, subtracted->GetMostSigOp() }.Simplify()),
-            *(subtracted->GetLeastSigOp().Simplify()) }
-                       .Simplify();
-        return Add { *simplifiedMinuend, *RHS }.Simplify();
+        auto m = Multiply<Expression> { Real { -1.0 }, subtracted->GetMostSigOp() };
+        auto ms = m.Accept(simplifyVisitor);
+        if (!ms) {
+            return Add { *simplifiedMinuend, Add { m, *(subtracted->GetLeastSigOp().Simplify()) } }.Generalize();
+        }
+        auto RHS = Add { *(std::move(ms.value())), *(subtracted->GetLeastSigOp().Simplify()) };
+        auto s = RHS.Accept(simplifyVisitor);
+        if (!s) {
+            return Add { *simplifiedMinuend, RHS }.Generalize();
+        }
+        return Add { *simplifiedMinuend, *(std::move(s.value())) }.Accept(simplifyVisitor).value();
     }
     return Add { *simplifiedMinuend, negated }.Simplify();
 }
 
 auto Subtract<Expression>::Differentiate(const Expression& differentiationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // Single diff variable
     if (auto variable = RecursiveCast<Variable>(differentiationVariable); variable != nullptr) {
-        auto simplifiedSub = this->Simplify();
+        auto s = this->Accept(simplifyVisitor);
+        if (!s) {
+            return this->Generalize();
+        }
+        auto simplifiedSub = std::move(s).value();
 
         // Make sure we're still subtracting
         if (auto adder = RecursiveCast<Subtract<Expression>>(*simplifiedSub); adder != nullptr) {
@@ -122,7 +143,13 @@ auto Subtract<Expression>::Differentiate(const Expression& differentiationVariab
                 return Copy();
             }
 
-            return std::make_unique<Subtract<Expression, Expression>>(Subtract<Expression, Expression> { *(specializedLeft->Copy()), *(specializedRight->Copy()) })->Simplify();
+            auto us = std::make_unique<Subtract<Expression, Expression>>(Subtract<Expression, Expression> { *(specializedLeft->Copy()), *(specializedRight->Copy()) });
+            auto s1 = us->Accept(simplifyVisitor);
+            if (!s1) {
+                return us;
+            }
+            return std::move(s1).value();
+
         }
         // If not, use other differentiation technique
         else {
@@ -134,9 +161,14 @@ auto Subtract<Expression>::Differentiate(const Expression& differentiationVariab
 
 auto Subtract<Expression>::Integrate(const Expression& integrationVariable) const -> std::unique_ptr<Expression>
 {
+    SimplifyVisitor simplifyVisitor {};
     // Single integration variable
     if (auto variable = RecursiveCast<Variable>(integrationVariable); variable != nullptr) {
-        auto simplifiedSub = this->Simplify();
+        auto s = this->Accept(simplifyVisitor);
+        if (!s) {
+            return this->Generalize();
+        }
+        auto simplifiedSub = std::move(s).value();
 
         // Make sure we're still subtracting
         if (auto adder = RecursiveCast<Subtract<Expression>>(*simplifiedSub); adder != nullptr) {
diff --git a/tests/AddTests.cpp b/tests/AddTests.cpp
index 202eb14c..0788e033 100644
--- a/tests/AddTests.cpp
+++ b/tests/AddTests.cpp
@@ -12,6 +12,9 @@
 #include "Oasis/Real.hpp"
 #include "Oasis/RecursiveCast.hpp"
 #include "Oasis/Variable.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Addition", "[Add]")
 {
@@ -24,7 +27,7 @@ TEST_CASE("Addition", "[Add]")
 
     OASIS_CAPTURE_WITH_SERIALIZER(add);
 
-    auto simplified = add.Simplify();
+    auto simplified = add.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -42,7 +45,7 @@ TEST_CASE("Symbolic Addition", "[Add][Symbolic]")
                 Oasis::Variable { "x" } }
         };
 
-    auto simplified = add.Simplify();
+    auto simplified = add.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Multiply>());
 
     REQUIRE(Oasis::Multiply {
@@ -61,7 +64,7 @@ TEST_CASE("Symbolic Addition, variable case", "[Add][Symbolic]")
                 Oasis::Variable { "x" } }
         };
 
-    auto simplified = add.Simplify();
+    auto simplified = add.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Add {
         Oasis::Real { 1.0 },
@@ -79,7 +82,7 @@ TEST_CASE("Generalized Addition", "[Add][Generalized]")
         Oasis::Real { 3.0 }
     };
 
-    auto simplified = add.Simplify();
+    auto simplified = add.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -126,7 +129,7 @@ TEST_CASE("Imaginary Addition", "[Imaginary][Add]")
         Oasis::Imaginary {}
     };
 
-    auto spec1 = a1.Simplify();
+    auto spec1 = a1.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Multiply { Oasis::Real { 2.0 }, Oasis::Imaginary {} }.Equals(*spec1));
 }
@@ -146,7 +149,7 @@ TEST_CASE("Addition Associativity", "[Add][Associative]")
                     Oasis::Variable { "x" } } } }
     };
 
-    auto simplified1 = assoc1.Simplify();
+    auto simplified1 = assoc1.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Add {
         Oasis::Real { 5.0 },
@@ -177,7 +180,7 @@ TEST_CASE("Add Associativity with wider tree", "[Add][Associativity]")
                     Oasis::Variable { "x" } } } }
     };
 
-    auto simplified2 = assoc2.Simplify();
+    auto simplified2 = assoc2.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Add {
         Oasis::Add {
@@ -223,6 +226,6 @@ TEST_CASE("Variadic Add Constructor", "[Add]")
 
     const Oasis::Real expected { 40.0 };
 
-    const auto simplified = add.Simplify();
+    const auto simplified = add.Accept(simplifyVisitor).value();
     REQUIRE(expected.Equals(*simplified));
 }
\ No newline at end of file
diff --git a/tests/BinaryExpressionTests.cpp b/tests/BinaryExpressionTests.cpp
index 8d8c7bcb..87cd97fe 100644
--- a/tests/BinaryExpressionTests.cpp
+++ b/tests/BinaryExpressionTests.cpp
@@ -9,6 +9,9 @@
 #include "Oasis/Real.hpp"
 #include "Oasis/RecursiveCast.hpp"
 #include "Oasis/Variable.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Recursive Cast Considers Commutative Property", "[Symbolic]")
 {
@@ -231,5 +234,5 @@ TEST_CASE("Substitute Binary", "[Substitute]")
 
                     auto after = before.Substitute(Oasis::Variable { "x" }, Oasis::Real { 4.0 }); // after should some std::unique_ptr<Expression> such that it equals 2(4) + 3(4)
                     Oasis::Real twenty {20};
-    REQUIRE(after->Equals(*(twenty.Simplify())));
+    REQUIRE(after->Equals(*(twenty.Accept(simplifyVisitor).value())));
 }
\ No newline at end of file
diff --git a/tests/DifferentiateTests.cpp b/tests/DifferentiateTests.cpp
index 9d893935..65ab86d5 100644
--- a/tests/DifferentiateTests.cpp
+++ b/tests/DifferentiateTests.cpp
@@ -13,6 +13,9 @@
 #include "Oasis/Derivative.hpp"
 #include "Oasis/Log.hpp"
 #include "Oasis/EulerNumber.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Differentiate Nonzero number", "[Differentiate][Real][Nonzero]")
 {
@@ -62,7 +65,8 @@ TEST_CASE("Differentiate Power Rule", "[Differentiate][Exponent][Power]")
         Oasis::Multiply {
             Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 2.0f } },
             Oasis::Real { 3.0f } }};
-    auto ptr = diff1.Simplify();
+    auto ptr = diff1.Accept(simplifyVisitor).value();
+
     auto diffed = differentiate.Differentiate(var);
     REQUIRE((*diffed).Equals(*ptr));
 }
@@ -73,7 +77,7 @@ TEST_CASE("Differentiate Constant Rule Multiply", "[Differentiate][Multiply][Con
     Oasis::Multiply<Oasis::Real, Oasis::Variable> diff1 { Oasis::Real { 3.0f }, Oasis::Variable { var.GetName() } };
     Oasis::Real three {3.0};
 
-    auto ptr = three.Simplify();
+    auto ptr = three.Accept(simplifyVisitor).value();
     auto diffed = diff1.Differentiate(var);
     REQUIRE((diffed->Equals(*ptr)));
 }
@@ -84,8 +88,8 @@ TEST_CASE("Differentiate Constant Rule Divide", "[Differentiate][Divide][Constan
     Oasis::Divide<Oasis::Variable, Oasis::Real> diff1 { Oasis::Variable { var.GetName() }, Oasis::Real { 2.0f } };
     Oasis::Real half {0.5f};
     auto diffed = diff1.Differentiate(var);
-    auto simplified = diffed->Simplify();
-    REQUIRE((simplified->Equals(*(half.Simplify()))));
+    auto simplified = diffed->Accept(simplifyVisitor).value();
+    REQUIRE((simplified->Equals(*(half.Accept(simplifyVisitor).value()))));
 }
 
 TEST_CASE("Differentiate Add Rule Different Terms", "[Differentiate][Add][Different]")
@@ -95,9 +99,9 @@ TEST_CASE("Differentiate Add Rule Different Terms", "[Differentiate][Add][Differ
 
     Oasis::Real one {1};
     auto diffed = diff1.Differentiate(var);
-    auto simplified = diffed->Simplify();
+    auto simplified = diffed->Accept(simplifyVisitor).value();
 
-    REQUIRE(simplified->Equals(*(one.Simplify())));
+    REQUIRE(simplified->Equals(*(one.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Differentiate Subtract Rule Different Terms", "[Differentiate][Subtract][Different]")
@@ -106,9 +110,9 @@ TEST_CASE("Differentiate Subtract Rule Different Terms", "[Differentiate][Subtra
     Oasis::Subtract<Oasis::Variable, Oasis::Real> diff1 { Oasis::Variable { var.GetName() }, Oasis::Real { 2.0f } };
     Oasis::Real one {1};
     auto diffed = diff1.Differentiate(var);
-    auto simplified = diffed->Simplify();
+    auto simplified = diffed->Accept(simplifyVisitor).value();
 
-    REQUIRE(simplified->Equals(*(one.Simplify())));
+    REQUIRE(simplified->Equals(*(one.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Differentiate Add Rule Like Terms", "[Differentiate][Add][Like]")
@@ -119,9 +123,9 @@ TEST_CASE("Differentiate Add Rule Like Terms", "[Differentiate][Add][Like]")
     Oasis::Real two {2};
 
     auto diffed = diff1.Differentiate(var);
-    auto simplified = diffed->Simplify();
+    auto simplified = diffed->Accept(simplifyVisitor).value();
 
-    REQUIRE(simplified->Equals(*(two.Simplify())));
+    REQUIRE(simplified->Equals(*(two.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Differentiate Quotient Rule Like Terms", "[Differentiate][Divide][Like]")
@@ -147,9 +151,9 @@ TEST_CASE("Differentiate Quotient Rule Like Terms", "[Differentiate][Divide][Lik
             {Oasis::Variable {var.GetName() }, Oasis::Real
             {3}}}};
     auto diffed = diff1.Differentiate(var);
-    auto simplified = diffed->Simplify();
+    auto simplified = diffed->Accept(simplifyVisitor).value();
 
-    auto simplifiedAns = answer.Simplify();
+    auto simplifiedAns = answer.Accept(simplifyVisitor).value();
 
     REQUIRE(simplified->Equals(*(simplifiedAns)));
 }
@@ -161,8 +165,8 @@ TEST_CASE("Differentiate Multiple Variables Differentiate", "[Differentiate][Mul
     Oasis::Multiply<Oasis::Variable, Oasis::Variable> xy
     {Oasis::Variable {x.GetName()}, Oasis::Variable {y.GetName()}};
     auto diffed = xy.Differentiate(x);
-    auto simplified = diffed->Simplify();
-    REQUIRE(simplified->Equals(*(y.Simplify())));
+    auto simplified = diffed->Accept(simplifyVisitor).value();
+    REQUIRE(simplified->Equals(*(y.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Differentiate Multiple Variables + Real Differentiate", "[Differentiate][Multiply][Different]")
@@ -177,8 +181,8 @@ TEST_CASE("Differentiate Multiple Variables + Real Differentiate", "[Differentia
     Oasis::Multiply<Oasis::Real, Oasis::Variable> threex
             {Oasis::Real{3}, Oasis::Variable{x.GetName()}};
     auto diffed = threexy.Differentiate(y);
-    auto simplified = diffed->Simplify();
-    REQUIRE(simplified->Equals(*(threex.Simplify())));
+    auto simplified = diffed->Accept(simplifyVisitor).value();
+    REQUIRE(simplified->Equals(*(threex.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Add", "[Differentiate][Add][Like]")
@@ -188,8 +192,8 @@ TEST_CASE("Derivative Wrapper Class Add", "[Differentiate][Add][Like]")
     {Oasis::Add<Oasis::Variable, Oasis::Real> {Oasis::Variable {x.GetName()},
                                                Oasis::Real{2}}, Oasis::Variable {x.GetName()}};
     Oasis::Real one {1};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(one.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(one.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Subtract", "[Differentiate][Subtract][Like]")
@@ -199,8 +203,8 @@ TEST_CASE("Derivative Wrapper Class Subtract", "[Differentiate][Subtract][Like]"
             {Oasis::Subtract<Oasis::Variable, Oasis::Real> {Oasis::Variable {x.GetName()},
                                                        Oasis::Real{2}}, Oasis::Variable {x.GetName()}};
     Oasis::Real one {1};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(one.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(one.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Multiply", "[Differentiate][Multiply][Like]")
@@ -210,8 +214,8 @@ TEST_CASE("Derivative Wrapper Class Multiply", "[Differentiate][Multiply][Like]"
             {Oasis::Multiply<Oasis::Variable, Oasis::Real> {Oasis::Variable {x.GetName()},
                                                        Oasis::Real{2}}, Oasis::Variable {x.GetName()}};
     Oasis::Real two {2};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(two.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(two.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Divide", "[Differentiate][Divide][Like]")
@@ -221,8 +225,8 @@ TEST_CASE("Derivative Wrapper Class Divide", "[Differentiate][Divide][Like]")
             {Oasis::Divide<Oasis::Variable, Oasis::Real> {Oasis::Variable {x.GetName()},
                                                        Oasis::Real{2}}, Oasis::Variable {x.GetName()}};
     Oasis::Real half {0.5f};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(half.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(half.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Real", "[Differentiate][Real][Like]")
@@ -231,8 +235,8 @@ TEST_CASE("Derivative Wrapper Class Real", "[Differentiate][Real][Like]")
     Oasis::Derivative<Oasis::Real, Oasis::Variable> diff1
             {Oasis::Real {2}, Oasis::Variable {x.GetName()}};
     Oasis::Real zero {0};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(zero.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(zero.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Variable", "[Differentiate][Variable][Like]")
@@ -241,8 +245,8 @@ TEST_CASE("Derivative Wrapper Class Variable", "[Differentiate][Variable][Like]"
     Oasis::Derivative<Oasis::Variable, Oasis::Variable> diff1
             {Oasis::Variable {x.GetName()}, Oasis::Variable {x.GetName()}};
     Oasis::Real one {1};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(one.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(one.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Derivative Wrapper Class Exponent", "[Differentiate][Exponent][Like]")
@@ -252,8 +256,8 @@ TEST_CASE("Derivative Wrapper Class Exponent", "[Differentiate][Exponent][Like]"
             {Oasis::Exponent<Oasis::Variable, Oasis::Real> {Oasis::Variable {x.GetName()},
                                                        Oasis::Real{2}}, Oasis::Variable {x.GetName()}};
     Oasis::Multiply<Oasis::Variable, Oasis::Real> twox {Oasis::Variable {x.GetName()}, Oasis::Real {2}};
-    auto diffed = diff1.Simplify();
-    REQUIRE(diffed->Equals(*(twox.Simplify())));
+    auto diffed = diff1.Accept(simplifyVisitor).value();
+    REQUIRE(diffed->Equals(*(twox.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Product Rule", "[Differentiate][Product]")
@@ -261,7 +265,7 @@ TEST_CASE("Product Rule", "[Differentiate][Product]")
     Oasis::Variable f{"x"};
     Oasis::Multiply g{Oasis::Real{6.0}, Oasis::Variable{"x"}};
     Oasis::Derivative diff{Oasis::Multiply{f, g}, Oasis::Variable{"x"}};
-    auto simplified = diff.Simplify();
+    auto simplified = diff.Accept(simplifyVisitor).value();
     Oasis::Multiply expected{Oasis::Real{12.0}, Oasis::Variable{"x"}};
 
     REQUIRE(simplified->Equals(expected));
@@ -274,8 +278,8 @@ TEST_CASE("Natural Exponential Derivative", "[Derivative][Exponent][Euler's Numb
     Oasis::Derivative diffExp{exp, Oasis::Variable{"x"}};
     Oasis::Derivative diffExp2{exp2, Oasis::Variable{"x"}};
     Oasis::Multiply expected2{Oasis::Real{2.0}, exp2};
-    auto simplified = diffExp.Simplify();
-    auto simplified2 = diffExp2.Simplify();
+    auto simplified = diffExp.Accept(simplifyVisitor).value();
+    auto simplified2 = diffExp2.Accept(simplifyVisitor).value();
 
     REQUIRE(simplified->Equals(exp));
     REQUIRE(simplified2->Equals(expected2));
@@ -285,7 +289,8 @@ TEST_CASE("Natural Logarithm Derivative", "[Derivative][Logarithm][Euler's Numbe
 {
     Oasis::Derivative diffLog{
         Oasis::Log{Oasis::EulerNumber{},Oasis::Multiply{Oasis::Real{6},Oasis::Variable{"x"}}}, Oasis::Variable{"x"}};
-    auto diff = diffLog.Simplify();
+    auto raw = diffLog.Accept(simplifyVisitor);
+    auto diff = std::move(raw.value());
     Oasis::Divide expected{Oasis::Real{1.0}, Oasis::Variable{"x"}};
     REQUIRE(diff->Equals(expected));
 }
@@ -294,17 +299,17 @@ TEST_CASE("Real Base Logarithm Derivative", "[Derivative][Logarithm][Euler's Rea
 {
     Oasis::Derivative diffLog{
         Oasis::Log{Oasis::Real{10.0},Oasis::Multiply{Oasis::Real{6},Oasis::Variable{"x"}}}, Oasis::Variable{"x"}};
-    auto diff = diffLog.Simplify();
+    auto diff = diffLog.Accept(simplifyVisitor).value();
     Oasis::Divide expected{Oasis::Real{1.0}, Oasis::Multiply{Oasis::Variable{"x"},
                                                       Oasis::Log{Oasis::EulerNumber{}, Oasis::Real{10.0}}}};
-    REQUIRE(diff->Equals(*expected.Simplify()));
+    REQUIRE(diff->Equals(*expected.Accept(simplifyVisitor).value()));
 }
 
 TEST_CASE("Variable Base Logarithm Derivative", "[Derivative][Logarithm][Variable]")
 {
     Oasis::Derivative diffLog{
         Oasis::Log{Oasis::Variable{"y"},Oasis::Multiply{Oasis::Real{6},Oasis::Variable{"x"}}}, Oasis::Variable{"x"}};
-    auto diff = diffLog.Simplify();
+    auto diff = diffLog.Accept(simplifyVisitor).value();
     Oasis::Divide expected{Oasis::Real{1.0}, Oasis::Multiply{Oasis::Variable{"x"},
                                                       Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"y"}}}};
     REQUIRE(diff->Equals(expected));
@@ -315,7 +320,7 @@ TEST_CASE("Any Base Exponential Derivative", "[Derivative][Exponent]")
     Oasis::Exponent exp{Oasis::Variable{"a"}, Oasis::Multiply{Oasis::Real{2.0},Oasis::Variable{"x"}}};
     Oasis::Derivative diffExp{exp, Oasis::Variable{"x"}};
     Oasis::Multiply expected{Oasis::Multiply{Oasis::Real{2.0}, exp}, Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"a"}}};
-    auto simplified = diffExp.Simplify();
+    auto simplified = diffExp.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Equals(expected));
 }
 
@@ -325,7 +330,7 @@ TEST_CASE("Variable Base Exponential Derivative", "[Derivative][Exponent]")
     Oasis::Derivative diffExp{exp, Oasis::Variable{"x"}};
     Oasis::Multiply expected{exp , Oasis::Add{ Oasis::Multiply {Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"x"}},
                                     Oasis::Real{2.0}},Oasis::Real{2.0}},};
-    auto simplified = diffExp.Simplify();
+    auto simplified = diffExp.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Equals(expected));
 }
 
@@ -341,15 +346,17 @@ TEST_CASE("General Exponential Derivative", "[Derivative][Exponent]")
                                     Oasis::Multiply { Oasis::Real{4.0} ,
                                         Oasis::Log{Oasis::EulerNumber{},
                                             Oasis::Add { Oasis::Variable{"x"}, Oasis::Real{3.0} }}} } };
-    auto simplified = diffExp.Simplify();
+    auto simplified = diffExp.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Equals(expected));
 }
 
-TEST_CASE("Expresion Base Exponential Derivative", "[Derivative][Exponent]")
+TEST_CASE("Expression Base Exponential Derivative", "[Derivative][Exponent]")
 {
     Oasis::Exponent exp{ Oasis::Add { Oasis::Variable{"x"}, Oasis::Real{3.0} } , Oasis::Real{2.0}};
     Oasis::Derivative diffExp{exp, Oasis::Variable{"x"}};
-    Oasis::Multiply expected{ Oasis::Real{2.0} , Oasis::Add {Oasis::Variable{"x"} , Oasis::Real {3.0}} };
-    auto simplified = diffExp.Simplify();
+
+    Oasis::Multiply expected{ Oasis::Add {Oasis::Variable{"x"} , Oasis::Real {3.0} }, Oasis::Real{ 2.0}    };
+
+    auto simplified = diffExp.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Equals(expected));
 }
\ No newline at end of file
diff --git a/tests/DivideTests.cpp b/tests/DivideTests.cpp
index 3658faea..1e8a66a8 100644
--- a/tests/DivideTests.cpp
+++ b/tests/DivideTests.cpp
@@ -11,6 +11,9 @@
 #include "Oasis/Multiply.hpp"
 #include "Oasis/Exponent.hpp"
 #include "Oasis/Add.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Division", "[Divide]")
 {
@@ -21,7 +24,7 @@ TEST_CASE("Division", "[Divide]")
         Oasis::Real { 2.0 }
     };
 
-    auto simplified = subtract.Simplify();
+    auto simplified = subtract.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -55,12 +58,10 @@ TEST_CASE("Symbolic Division, equal variables", "[Division][Symbolic]")
             }
         };
 
-    auto simplified = div.Simplify();
-    REQUIRE(simplified);
+    auto simplified = div.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Equals(Oasis::Real { 2.0 }));
 
-    auto simplified2 = div2.Simplify();
-    REQUIRE(simplified2);
+    auto simplified2 = div2.Accept(simplifyVisitor).value();
     REQUIRE(simplified2->Equals(Oasis::Real { 2.0 }));
 }
 
@@ -80,7 +81,7 @@ TEST_CASE("Symbolic Division, unequal variables", "[Division][Symbolic]")
                     Oasis::Variable { "z" } }}
         };
 
-    auto simplified = div.Simplify();
+    auto simplified = div.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Divide {
         Oasis::Multiply {
@@ -114,7 +115,7 @@ TEST_CASE("Symbolic Division of Expressions", "[Division][Symbolic]")
                     Oasis::Real { 1 } }}
         };
 
-    auto simplified = div.Simplify();
+    auto simplified = div.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Multiply {
         Oasis::Real{ 2 },
@@ -142,7 +143,7 @@ TEST_CASE("Symbolic Division, unequal exponents", "[Division][Symbolic]")
                     Oasis::Variable { "z" } }}
         };
 
-    auto simplified = div.Simplify();
+    auto simplified = div.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Divide {
         Oasis::Multiply {
@@ -167,7 +168,7 @@ TEST_CASE("Symbolic Division, equal exponents", "[Division][Symbolic]")
                     Oasis::Variable { "z" } }}
         };
 
-    auto simplified = div.Simplify();
+    auto simplified = div.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Divide {
         Oasis::Multiply {
@@ -183,7 +184,7 @@ TEST_CASE("Division of equal variables", "[Division][Symbolic]")
     // x/x = 1
 
     Oasis::Divide div { Oasis::Variable { "x" }, Oasis::Variable { "x" } };
-    auto simplified = div.Simplify();
+    auto simplified = div.Accept(simplifyVisitor).value();
     REQUIRE(Oasis::Real { 1 }.Equals(*simplified));
 }
 
@@ -197,7 +198,7 @@ TEST_CASE("Generalized Division", "[Divide][Generalized]")
         Oasis::Real { 2.0 }
     };
 
-    auto simplified = subtract.Simplify();
+    auto simplified = subtract.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
diff --git a/tests/ExponentTests.cpp b/tests/ExponentTests.cpp
index b36da364..827ab282 100644
--- a/tests/ExponentTests.cpp
+++ b/tests/ExponentTests.cpp
@@ -12,6 +12,9 @@
 #include "Oasis/Real.hpp"
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Variable.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Zero Rule", "[Exponent][Zero]")
 {
@@ -20,7 +23,7 @@ TEST_CASE("Zero Rule", "[Exponent][Zero]")
         Oasis::Real { 0.0 }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -34,7 +37,7 @@ TEST_CASE("Second Zero Rule", "[Exponent][Zero]")
         Oasis::Real { 5.0 }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -48,7 +51,7 @@ TEST_CASE("Symbolic Zero Rule", "[Exponent][Zero][Symbolic]")
         Oasis::Real { 0.0 }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -62,7 +65,7 @@ TEST_CASE("Second Symbolic Zero Rule", "[Exponent][Zero][Symbolic]")
         Oasis::Variable { "x" }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -76,7 +79,7 @@ TEST_CASE("Symbolic One Rule", "[Exponent][One][Symbolic]")
         Oasis::Variable { "x" }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -90,7 +93,7 @@ TEST_CASE("Whole Number Exponentiation", "[Exponent][Numerical]")
         Oasis::Real { 2.0 }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -104,7 +107,7 @@ TEST_CASE("Fractional Exponentiation", "[Exponent][Root]")
         Oasis::Real { 0.5 }
     };
 
-    auto simplified = exponent.Simplify();
+    auto simplified = exponent.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -176,12 +179,12 @@ TEST_CASE("Addition of Exponents", "[Add][Exponent][Symbolic]")
                 Oasis::Real { 2.0 } } }
     };
 
-    auto simplified = add1.Simplify();
+    auto simplified = add1.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Multiply>());
 
-    auto simplified2 = add2.Simplify();
-    auto simplified3 = add3.Simplify();
-    auto simplified4 = add4.Simplify();
+    auto simplified2 = add2.Accept(simplifyVisitor).value();
+    auto simplified3 = add3.Accept(simplifyVisitor).value();
+    auto simplified4 = add4.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Multiply {
         Oasis::Real { 2.0 },
@@ -239,10 +242,10 @@ TEST_CASE("Variable Multiplication", "[Exponent][Variable][Multiplication]")
             Oasis::Variable { "x" } }
     };
 
-    auto simplified = expr.Simplify();
-    auto simplified2 = expr2.Simplify();
-    auto simplified3 = expr3.Simplify();
-    auto simplified4 = expr4.Simplify();
+    auto simplified = expr.Accept(simplifyVisitor).value();
+    auto simplified2 = expr2.Accept(simplifyVisitor).value();
+    auto simplified3 = expr3.Accept(simplifyVisitor).value();
+    auto simplified4 = expr4.Accept(simplifyVisitor).value();
 
     REQUIRE(simplified->Is<Oasis::Exponent>());
     REQUIRE(Oasis::Exponent<Oasis::Variable, Oasis::Real> { Oasis::Variable { "x" }, Oasis::Real { 2.0 } }.Equals(*simplified));
@@ -337,14 +340,14 @@ TEST_CASE("Variable with power Multiplication", "[Exponent][Variable][Multiplica
                 Oasis::Real { 2.0 } } }
     };
 
-    auto simplified = expr.Simplify();
-    auto simplified2 = expr2.Simplify();
-    auto simplified3 = expr3.Simplify();
-    auto simplified4 = expr4.Simplify();
-    auto simplified5 = expr5.Simplify();
-    auto simplified6 = expr6.Simplify();
-    auto simplified7 = expr7.Simplify();
-    auto simplified8 = expr8.Simplify();
+    auto simplified = expr.Accept(simplifyVisitor).value();
+    auto simplified2 = expr2.Accept(simplifyVisitor).value();
+    auto simplified3 = expr3.Accept(simplifyVisitor).value();
+    auto simplified4 = expr4.Accept(simplifyVisitor).value();
+    auto simplified5 = expr5.Accept(simplifyVisitor).value();
+    auto simplified6 = expr6.Accept(simplifyVisitor).value();
+    auto simplified7 = expr7.Accept(simplifyVisitor).value();
+    auto simplified8 = expr8.Accept(simplifyVisitor).value();
 
     REQUIRE(simplified->Is<Oasis::Exponent>());
     REQUIRE(Oasis::Exponent<Oasis::Variable, Oasis::Real> { Oasis::Variable { "x" }, Oasis::Real { 4.0 } }.Equals(*simplified));
@@ -413,17 +416,17 @@ TEST_CASE("Subtraction of Exponents", "[Subtract][Exponent][Symbolic]")
                 Oasis::Real { 2.0 } } }
     };
 
-    auto simplified = sub1.Simplify();
-    auto simplified2 = sub2.Simplify();
-    auto simplified3 = sub3.Simplify();
-    auto simplified4 = sub4.Simplify();
+    auto simplified = sub1.Accept(simplifyVisitor).value();
+    auto simplified2 = sub2.Accept(simplifyVisitor).value();
+    auto simplified3 = sub3.Accept(simplifyVisitor).value();
+    auto simplified4 = sub4.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Real { 0.0 }.Equals(*simplified));
     REQUIRE((Oasis::Multiply {
         Oasis::Real { 1.0 },
         Oasis::Exponent {
             Oasis::Variable { "x" },
-            Oasis::Real { 2.0 } } }).Simplify()->Equals(*simplified2));
+            Oasis::Real { 2.0 } } }).Accept(simplifyVisitor).value()->Equals(*simplified2));
     REQUIRE(Oasis::Multiply {
         Oasis::Real { -1.0 },
         Oasis::Exponent {
@@ -464,11 +467,11 @@ TEST_CASE("Imaginary Exponent Rule", "[Imaginary][Exponent]")
         Oasis::Real { 0.5 }
     };
 
-    auto simplified1 = i1.Simplify();
-    auto simplified2 = i2.Simplify();
-    auto simplified3 = i3.Simplify();
-    auto simplified4 = i4.Simplify();
-    auto simplifiedimg = img.Simplify();
+    auto simplified1 = i1.Accept(simplifyVisitor).value();
+    auto simplified2 = i2.Accept(simplifyVisitor).value();
+    auto simplified3 = i3.Accept(simplifyVisitor).value();
+    auto simplified4 = i4.Accept(simplifyVisitor).value();
+    auto simplifiedimg = img.Accept(simplifyVisitor).value();
 
     REQUIRE(i.Is<Oasis::Imaginary>());
     REQUIRE(Oasis::Imaginary {}.Equals(*simplified1));
diff --git a/tests/IntegrateTests.cpp b/tests/IntegrateTests.cpp
index db57d733..7e1c2ede 100644
--- a/tests/IntegrateTests.cpp
+++ b/tests/IntegrateTests.cpp
@@ -7,6 +7,9 @@
 #include "Oasis/Real.hpp"
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Variable.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Integrate Nonzero number", "[Integrate][Real][Nonzero]")
 {
@@ -44,7 +47,7 @@ TEST_CASE("Integrate Same Variable", "[Integrate][Variable][Same]")
         Oasis::Variable { "C" }
     };
 
-    auto ptr = integral.Simplify();
+    auto ptr = integral.Accept(simplifyVisitor).value();
     auto integrated = var.Integrate(var);
     REQUIRE(ptr->Equals(*integrated));
 }
@@ -73,7 +76,7 @@ TEST_CASE("Integrate Power Rule", "[Integrate][Exponent][Power]")
             Oasis::Exponent { Oasis::Variable { var.GetName() }, Oasis::Real { 3.0f } },
             Oasis::Real { 3.0f } },
         Oasis::Variable { "C" } } };
-    auto ptr = integral.Simplify();
+    auto ptr = integral.Accept(simplifyVisitor).value();
 
     auto integrated = integrand.Integrate(var);
     REQUIRE((*integrated).Equals(*ptr));
@@ -95,7 +98,7 @@ TEST_CASE("Integrate Constant Rule Multiply", "[Integrate][Multiply][Constant]")
                 "C" } }
     };
 
-    auto ptr = integral.Simplify();
+    auto ptr = integral.Accept(simplifyVisitor).value();
     auto integrated = integrand.Integrate(var);
     REQUIRE((integrated->Equals(*ptr)));
 }
@@ -115,9 +118,9 @@ TEST_CASE("Integrate Constant Rule Divide", "[Integrate][Divide][Constant]")
             Oasis::Variable { "C" } }
     };
     auto integrated = integrand.Integrate(var);
-    auto simplified = integrated->Simplify();
+    auto simplified = integrated->Accept(simplifyVisitor).value();
 
-    REQUIRE((simplified->Equals(*(integral.Simplify()))));
+    REQUIRE((simplified->Equals(*(integral.Accept(simplifyVisitor).value()))));
 }
 
 TEST_CASE("Integrate Add Rule Different Terms", "[Integrate][Add][Different]")
@@ -137,9 +140,9 @@ TEST_CASE("Integrate Add Rule Different Terms", "[Integrate][Add][Different]")
             Oasis::Variable { "C" } }
     };
     auto integrated = integrand.Integrate(var);
-    auto simplified = integrated->Simplify();
+    auto simplified = integrated->Accept(simplifyVisitor).value();
 
-    REQUIRE(simplified->Equals(*(integral.Simplify())));
+    REQUIRE(simplified->Equals(*(integral.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Integrate Subtract Rule Different Terms", "[Integrate][Subtract][Different]")
@@ -160,7 +163,7 @@ TEST_CASE("Integrate Subtract Rule Different Terms", "[Integrate][Subtract][Diff
     };
     auto integrated = integrand.Integrate(var);
 
-    REQUIRE(integrated->Equals(*(integral.Simplify())));
+    REQUIRE(integrated->Equals(*(integral.Accept(simplifyVisitor).value())));
 }
 
 TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]")
@@ -175,7 +178,7 @@ TEST_CASE("Integrate Add Rule Like Terms", "[Integrate][Add][Like]")
     };
 
     auto integrated = integrand.Integrate(var);
-    auto simplified = integrated->Simplify();
+    auto simplified = integrated->Accept(simplifyVisitor).value();
 
-    REQUIRE(simplified->Equals(*(integral.Simplify())));
+    REQUIRE(simplified->Equals(*(integral.Accept(simplifyVisitor).value())));
 }
\ No newline at end of file
diff --git a/tests/LogTests.cpp b/tests/LogTests.cpp
index ac9fe29c..54e0266e 100644
--- a/tests/LogTests.cpp
+++ b/tests/LogTests.cpp
@@ -17,9 +17,12 @@
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Undefined.hpp"
 #include "Oasis/Variable.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
 
 #define EPSILON 1E-6
 
+inline Oasis::SimplifyVisitor simplifyVisitor{};
+
 TEST_CASE("Log with invalid base", "[Log][Invalid Base]")
 {
     Oasis::Log negativeBase = Oasis::Log {
@@ -31,10 +34,10 @@ TEST_CASE("Log with invalid base", "[Log][Invalid Base]")
         Oasis::Real { 10.0 }
     };
 
-    auto negativeBaseSimplified = negativeBase.Simplify();
+    auto negativeBaseSimplified = negativeBase.Accept(simplifyVisitor).value();
     REQUIRE(negativeBaseSimplified->Is<Oasis::Undefined>());
 
-    auto baseOneSimplified = baseOne.Simplify();
+    auto baseOneSimplified = baseOne.Accept(simplifyVisitor).value();
     REQUIRE(baseOneSimplified->Is<Oasis::Undefined>());
 }
 
@@ -45,7 +48,7 @@ TEST_CASE("Log with invalid argument", "[Log][Invalid Argument]")
         Oasis::Real { -5.0 }
     };
 
-    auto negativeArgumentSimplified = negativeArgument.Simplify();
+    auto negativeArgumentSimplified = negativeArgument.Accept(simplifyVisitor).value();
     REQUIRE(negativeArgumentSimplified->Is<Oasis::Undefined>());
 }
 
@@ -57,7 +60,7 @@ TEST_CASE("Log of Real Numbers", "[Log][Real]")
         Oasis::Real { 1.0 }
     };
 
-    auto base5of1_Simplified = base5of1.Simplify();
+    auto base5of1_Simplified = base5of1.Accept(simplifyVisitor).value();
     REQUIRE(base5of1_Simplified->Is<Oasis::Real>());
 
     auto simplifiedReal5_1 = dynamic_cast<Oasis::Real&>(*base5of1_Simplified);
@@ -69,7 +72,7 @@ TEST_CASE("Log of Real Numbers", "[Log][Real]")
         Oasis::Real { 5.0 }
     };
 
-    auto base5of5_Simplified = base5of5.Simplify();
+    auto base5of5_Simplified = base5of5.Accept(simplifyVisitor).value();
     REQUIRE(base5of5_Simplified->Is<Oasis::Real>());
 
     auto simplifiedReal5_5 = dynamic_cast<Oasis::Real&>(*base5of5_Simplified);
@@ -81,7 +84,7 @@ TEST_CASE("Log of Real Numbers", "[Log][Real]")
         Oasis::Real { 25.0 }
     };
 
-    auto base5of25_Simplified = base5of25.Simplify();
+    auto base5of25_Simplified = base5of25.Accept(simplifyVisitor).value();
     REQUIRE(base5of25_Simplified->Is<Oasis::Real>());
 
     auto simplifiedReal5_25 = dynamic_cast<Oasis::Real&>(*base5of25_Simplified);
@@ -99,7 +102,7 @@ TEST_CASE("Log of Exponentiation", "[Log][Exponent]")
         }
     };
 
-    auto base5of5toX_Simplified = base5of5toX.Simplify();
+    auto base5of5toX_Simplified = base5of5toX.Accept(simplifyVisitor).value();
     auto simplifiedSpecialized = RecursiveCast<Oasis::Multiply<Oasis::Real, Oasis::Variable>>(*base5of5toX_Simplified);
 
     REQUIRE(simplifiedSpecialized != nullptr);
@@ -123,7 +126,7 @@ TEST_CASE("Log of Exponentiation", "[Log][Exponent]")
         }
     };
 
-    auto logOfExp_Simplified = logOfExp.Simplify();
+    auto logOfExp_Simplified = logOfExp.Accept(simplifyVisitor).value();
     REQUIRE(logOfExpResult.Equals(*logOfExp_Simplified));
 }
 
@@ -149,14 +152,14 @@ TEST_CASE("Sum of logs", "[Log][Add]")
     auto logOf4Simplified = dynamic_cast<Oasis::Real&>(*Oasis::Add {
         logOf2,
         logOf2
-    }.Simplify());
+    }.Accept(simplifyVisitor).value());
     auto logOf6Simplified = dynamic_cast<Oasis::Real&>(*Oasis::Add {
         logOf2,
         logOf3
-    }.Simplify());
+    }.Accept(simplifyVisitor).value());
     
-    REQUIRE_THAT(logOf4Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*logOf4.Simplify()).GetValue(), EPSILON));
-    REQUIRE_THAT(logOf6Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*logOf6.Simplify()).GetValue(), EPSILON));
+    REQUIRE_THAT(logOf4Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*logOf4.Accept(simplifyVisitor).value()).GetValue(), EPSILON));
+    REQUIRE_THAT(logOf6Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*logOf6.Accept(simplifyVisitor).value()).GetValue(), EPSILON));
 }
 
 TEST_CASE("Difference of logs", "[Log][Subtract]")
@@ -181,14 +184,14 @@ TEST_CASE("Difference of logs", "[Log][Subtract]")
     auto logOf2Simplified = Oasis::Subtract {
         logOf4,
         logOf2
-    }.Simplify();
+    }.Accept(simplifyVisitor).value();
     auto logOf3Simplified = dynamic_cast<Oasis::Real&>(*Oasis::Subtract {
         logOf6,
         logOf2
-    }.Simplify());
+    }.Accept(simplifyVisitor).value());
 
-    REQUIRE(logOf2Simplified->Equals(*(logOf2.Simplify())));
-    REQUIRE_THAT(logOf3Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*(logOf3.Simplify())).GetValue(), EPSILON));
+    REQUIRE(logOf2Simplified->Equals(*(logOf2.Accept(simplifyVisitor).value())));
+    REQUIRE_THAT(logOf3Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*(logOf3.Accept(simplifyVisitor).value())).GetValue(), EPSILON));
 }
 
 TEST_CASE("Change of Base", "[Log][Divide]")
@@ -218,28 +221,28 @@ TEST_CASE("Change of Base", "[Log][Divide]")
     auto logBase2Simplified = Oasis::Divide {
         logOf4,
         logOf2
-    }.Simplify();
+    }.Accept(simplifyVisitor).value();
     auto logBase3Simplified = dynamic_cast<Oasis::Real&>(*Oasis::Divide {
         logOf6,
         logOf3
-    }.Simplify());
+    }.Accept(simplifyVisitor).value());
 
     REQUIRE(logBase2Simplified->Equals(two));
-    REQUIRE_THAT(logBase3Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*(logBase3.Simplify())).GetValue(), EPSILON)); //PROBLEM
+    REQUIRE_THAT(logBase3Simplified.GetValue(), Catch::Matchers::WithinAbs(dynamic_cast<Oasis::Real&>(*(logBase3.Accept(simplifyVisitor).value())).GetValue(), EPSILON)); //PROBLEM
 }
 
 TEST_CASE("Log of Equal Base and Argument", "[Log][Expression]") {
     Oasis::Log log {Oasis::Variable{"x"},Oasis::Variable {"x"}};
-    auto simplifiedLog = log.Simplify();
+    auto simplifiedLog = log.Accept(simplifyVisitor);
     Oasis::Real expected {1};
-    REQUIRE(simplifiedLog->Equals(expected));
+    REQUIRE(simplifiedLog.value()->Equals(expected));
 }
 
 TEST_CASE("Undefined", "[UNDEFINED][UNDEFINED]") {
     REQUIRE(!Oasis::Undefined().Equals(Oasis::Undefined()));
-    REQUIRE(Oasis::Log<Oasis::Real, Oasis::Real>(Oasis::Real(-1.0), Oasis::Real(5.0)).Simplify()->Is<Oasis::Undefined>());
-    REQUIRE(Oasis::Log<Oasis::Real, Oasis::Real>(Oasis::Real(5.0), Oasis::Real(-5.0)).Simplify()->Is<Oasis::Undefined>());
-    REQUIRE(Oasis::Log<Oasis::Real, Oasis::Real>(Oasis::Real(1.0), Oasis::Real(5.0)).Simplify()->Is<Oasis::Undefined>());
+    REQUIRE(Oasis::Log<Oasis::Real, Oasis::Real>(Oasis::Real(-1.0), Oasis::Real(5.0)).Accept(simplifyVisitor).value()->Is<Oasis::Undefined>());
+    REQUIRE(Oasis::Log<Oasis::Real, Oasis::Real>(Oasis::Real(5.0), Oasis::Real(-5.0)).Accept(simplifyVisitor).value()->Is<Oasis::Undefined>());
+    REQUIRE(Oasis::Log<Oasis::Real, Oasis::Real>(Oasis::Real(1.0), Oasis::Real(5.0)).Accept(simplifyVisitor).value()->Is<Oasis::Undefined>());
 }
 
 TEST_CASE("Integral of Natural Log", "[Integral][Log][Euler]")
@@ -247,11 +250,14 @@ TEST_CASE("Integral of Natural Log", "[Integral][Log][Euler]")
     Oasis::Integral int1{Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"x"}}, Oasis::Variable{"x"}};
     Oasis::Integral int2{Oasis::Log{Oasis::EulerNumber{}, Oasis::Multiply{Oasis::Real{5}, Oasis::Variable{"x"}}}, Oasis::Variable{"x"}};
 
-    auto simp1 = int1.Simplify();
-    auto simp2 = int2.Simplify();
+    auto simp1 = int1.Accept(simplifyVisitor).value();
+    auto simp2 = int2.Accept(simplifyVisitor).value();
 
-    REQUIRE(simp1->Equals(*Oasis::Add {Oasis::Multiply{Oasis::Subtract{Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"x"}}, Oasis::Real{1}} , Oasis::Variable{"x"} }, Oasis::Variable{"C"} }.Simplify()));
-    REQUIRE(simp2->Equals(*Oasis::Add {Oasis::Multiply{Oasis::Subtract{Oasis::Log{Oasis::EulerNumber{}, Oasis::Multiply {Oasis::Real{5} ,Oasis::Variable{"x"}}}, Oasis::Real{1}} , Oasis::Variable{"x"} }, Oasis::Variable{"C"} }.Simplify()));
+    // TODO: FIX
+    auto expected1 = Oasis::Add {Oasis::Multiply{Oasis::Subtract{Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"x"}}, Oasis::Real{1}} , Oasis::Variable{"x"} }, Oasis::Variable{"C"} }.Accept(simplifyVisitor).value();
+    auto expected2 = Oasis::Add {Oasis::Multiply{Oasis::Subtract{Oasis::Log{Oasis::EulerNumber{}, Oasis::Multiply {Oasis::Real{5} ,Oasis::Variable{"x"}}}, Oasis::Real{1}} , Oasis::Variable{"x"} }, Oasis::Variable{"C"} }.Accept(simplifyVisitor).value();
+    REQUIRE(simp1->Equals(*expected1));
+    REQUIRE(simp2->Equals(*expected2));
     }
 
 TEST_CASE("Integral of Non-Natural variable base Log", "[Integral][Log][Variable]")
@@ -265,12 +271,12 @@ TEST_CASE("Integral of Non-Natural variable base Log", "[Integral][Log][Variable
     Oasis::Add eq2{Oasis::Divide{ Oasis::Multiply{Oasis::Variable{"x"}, Oasis::Subtract{Oasis::Log{Oasis::EulerNumber{}, Oasis::Multiply{Oasis::Real{5}, Oasis::Variable{"x"}}},Oasis::Real{1}}},
                         Oasis::Log{Oasis::EulerNumber{}, Oasis::Variable{"y"}}},Oasis::Variable{"C"}};
 
-    auto simp1 = int1.Simplify();
-    auto simp2 = int2.Simplify();
-    auto simpEq = intEq.Simplify();
+    auto simp1 = int1.Accept(simplifyVisitor).value();
+    auto simp2 = int2.Accept(simplifyVisitor).value();
+    auto simpEq = intEq.Accept(simplifyVisitor).value();
 
-    REQUIRE(simp1->Equals(*eq1.Simplify()));
-    REQUIRE(simp2->Equals(*eq2.Simplify()));
+    REQUIRE(simp1->Equals(*eq1.Accept(simplifyVisitor).value()));
+    REQUIRE(simp2->Equals(*eq2.Accept(simplifyVisitor).value()));
     REQUIRE(simpEq->Equals(Oasis::Add {Oasis::Variable{"x"}, Oasis::Variable{"C"} }));
 }
 
@@ -284,9 +290,9 @@ TEST_CASE("Integral of Non-natural real base log", "[Integral][Log][Variable]")
     Oasis::Add eq2{Oasis::Divide{ Oasis::Multiply{Oasis::Variable{"x"}, Oasis::Subtract{Oasis::Log{Oasis::EulerNumber{}, Oasis::Multiply{Oasis::Real{5}, Oasis::Variable{"x"}}}, Oasis::Real{1}}},
         Oasis::Log{Oasis::EulerNumber{}, Oasis::Real{10}}},Oasis::Variable{"C"}};
 
-    auto simp1 = int1.Simplify();
-    auto simp2 = int2.Simplify();
+    auto simp1 = int1.Accept(simplifyVisitor).value();
+    auto simp2 = int2.Accept(simplifyVisitor).value();
 
-    REQUIRE(simp1->Equals(*eq1.Simplify()));
-    REQUIRE(simp2->Equals(*eq2.Simplify()));
+    REQUIRE(simp1->Equals(*eq1.Accept(simplifyVisitor).value()));
+    REQUIRE(simp2->Equals(*eq2.Accept(simplifyVisitor).value()));
 }
\ No newline at end of file
diff --git a/tests/MagnitudeTests.cpp b/tests/MagnitudeTests.cpp
index bdcd63bb..a8b0d6d6 100644
--- a/tests/MagnitudeTests.cpp
+++ b/tests/MagnitudeTests.cpp
@@ -10,6 +10,9 @@
 #include "Oasis/Add.hpp"
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Matrix.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Magnitude of Real", "[Magnitude][Real]")
 {
@@ -18,10 +21,10 @@ TEST_CASE("Magnitude of Real", "[Magnitude][Real]")
     Oasis::Magnitude mag3{Oasis::Add{Oasis::Real{5}, Oasis::Real{6}}};
     Oasis::Magnitude mag4{Oasis::Subtract{Oasis::Real{5}, Oasis::Real{6}}};
 
-    auto simp1 = mag1.Simplify();
-    auto simp2 = mag2.Simplify();
-    auto simp3 = mag3.Simplify();
-    auto simp4 = mag4.Simplify();
+    auto simp1 = mag1.Accept(simplifyVisitor).value();
+    auto simp2 = mag2.Accept(simplifyVisitor).value();
+    auto simp3 = mag3.Accept(simplifyVisitor).value();
+    auto simp4 = mag4.Accept(simplifyVisitor).value();
 
     REQUIRE(simp1->Equals(Oasis::Real{5.0}));
     REQUIRE(simp2->Equals(Oasis::Real{75.0}));
@@ -35,9 +38,9 @@ TEST_CASE("Magnitude of Imaginary", "[Magnitude][Imaginary]")
     Oasis::Magnitude mag2{Oasis::Multiply{Oasis::Real{5.0}, Oasis::Imaginary{}}};
     Oasis::Magnitude mag3{Oasis::Multiply{Oasis::Real{-5.0}, Oasis::Imaginary{}}};
 
-    auto simp1 = mag1.Simplify();
-    auto simp2 = mag2.Simplify();
-    auto simp3 = mag3.Simplify();
+    auto simp1 = mag1.Accept(simplifyVisitor).value();
+    auto simp2 = mag2.Accept(simplifyVisitor).value();
+    auto simp3 = mag3.Accept(simplifyVisitor).value();
 
     REQUIRE(simp1->Equals(Oasis::Real{1.0}));
     REQUIRE(simp2->Equals(Oasis::Real{5}));
@@ -49,20 +52,20 @@ TEST_CASE("Magnitude of Complex", "[Magnitude][Real][Imaginary]")
     Oasis::Magnitude mag1{Oasis::Add{Oasis::Real{5.0}, Oasis::Imaginary{}}};
     Oasis::Magnitude mag2{Oasis::Subtract{Oasis::Real{5.0}, Oasis::Imaginary{}}};
 
-    auto simp1 = mag1.Simplify();
-    auto simp2 = mag2.Simplify();
+    auto simp1 = mag1.Accept(simplifyVisitor).value();
+    auto simp2 = mag2.Accept(simplifyVisitor).value();
 
-    REQUIRE(simp1->Equals(*Oasis::Exponent{Oasis::Real{26.0}, Oasis::Real{0.5}}.Simplify()));
-    REQUIRE(simp2->Equals(*Oasis::Exponent{Oasis::Real{26.0}, Oasis::Real{0.5}}.Simplify()));
+    REQUIRE(simp1->Equals(*Oasis::Exponent{Oasis::Real{26.0}, Oasis::Real{0.5}}.Accept(simplifyVisitor).value()));
+    REQUIRE(simp2->Equals(*Oasis::Exponent{Oasis::Real{26.0}, Oasis::Real{0.5}}.Accept(simplifyVisitor).value()));
 
     Oasis::Magnitude mag3{Oasis::Add{Oasis::Real{5.0}, Oasis::Multiply{Oasis::Real{4.0}, Oasis::Imaginary{}}}};
     Oasis::Magnitude mag4{Oasis::Subtract{Oasis::Real{7.0}, Oasis::Multiply{Oasis::Real{20.0}, Oasis::Imaginary{}}}};
 
-    auto simp3 = mag3.Simplify();
-    auto simp4 = mag4.Simplify();
+    auto simp3 = mag3.Accept(simplifyVisitor).value();
+    auto simp4 = mag4.Accept(simplifyVisitor).value();
 
-    REQUIRE(simp3->Equals(*Oasis::Exponent{Oasis::Real{41.0}, Oasis::Real{0.5}}.Simplify()));
-    REQUIRE(simp4->Equals(*Oasis::Exponent{Oasis::Real{449.0}, Oasis::Real{0.5}}.Simplify()));
+    REQUIRE(simp3->Equals(*Oasis::Exponent{Oasis::Real{41.0}, Oasis::Real{0.5}}.Accept(simplifyVisitor).value()));
+    REQUIRE(simp4->Equals(*Oasis::Exponent{Oasis::Real{449.0}, Oasis::Real{0.5}}.Accept(simplifyVisitor).value()));
 }
 
 TEST_CASE("Magnitude of Vector", "[Magnitude][Matrix1D]")
@@ -70,15 +73,16 @@ TEST_CASE("Magnitude of Vector", "[Magnitude][Matrix1D]")
     std::vector<double> vals = {5,6,7,8};
     Oasis::Magnitude mag1{Oasis::Matrix{4, 1, vals}};
     Oasis::Magnitude mag2{Oasis::Matrix{1, 4, vals}};
-    auto expected = Oasis::Exponent{Oasis::Real{174}, Oasis::Real{0.5}}.Simplify();
-    REQUIRE((mag1.Simplify())->Equals(*expected));
-    REQUIRE((mag2.Simplify())->Equals(*expected));
+    auto expected = Oasis::Exponent{Oasis::Real{174}, Oasis::Real{0.5}}.Accept(simplifyVisitor).value();
+    REQUIRE((mag1.Accept(simplifyVisitor).value())->Equals(*expected));
+    REQUIRE((mag2.Accept(simplifyVisitor).value())->Equals(*expected));
 }
 
 TEST_CASE("Magnitude of Matrix", "[Magnitude][MatrixXXD]")
 {
     std::vector<double> vals = {5,6,7,8};
     Oasis::Magnitude mag1{Oasis::Matrix{2, 2, vals}};
-    auto expected = Oasis::Exponent{Oasis::Real{174}, Oasis::Real{0.5}}.Simplify();
-    REQUIRE((mag1.Simplify())->Equals(*expected));
+    auto expected = Oasis::Exponent{Oasis::Real{174}, Oasis::Real{0.5}}.Accept(simplifyVisitor).value();
+    auto actual = mag1.Accept(simplifyVisitor).value();
+    REQUIRE(actual->Equals(*expected));
 }
\ No newline at end of file
diff --git a/tests/MatrixTests.cpp b/tests/MatrixTests.cpp
index 5ff1266c..ade15cb0 100644
--- a/tests/MatrixTests.cpp
+++ b/tests/MatrixTests.cpp
@@ -9,9 +9,12 @@
 #include "Oasis/Subtract.hpp"
 #include "Oasis/Multiply.hpp"
 #include "Oasis/RecursiveCast.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
 
 #define EPSILON 10E-6
 
+inline Oasis::SimplifyVisitor simplifyVisitor{};
+
 TEST_CASE("Create Empty Matrix", "[Matrix][Create][Empty]")
 {
     Oasis::Matrix a = {3,4};
@@ -67,7 +70,7 @@ TEST_CASE("Add Matrices", "[Matrix][Add]")
     Oasis::Matrix mat2{Oasis::MatrixXXD{{8,7},
                                         {6,5}}};
 
-    auto result = Oasis::Add<Oasis::Expression>{mat1, mat2}.Simplify();
+    auto result = Oasis::Add<Oasis::Expression>{mat1, mat2}.Accept(simplifyVisitor).value();
     auto spec = Oasis::RecursiveCast<Oasis::Matrix>(*result);
 
     REQUIRE(spec != nullptr);
@@ -82,7 +85,7 @@ TEST_CASE("Subtract Matrices", "[Matrix][Subtract]")
     Oasis::Matrix mat2{Oasis::MatrixXXD{{5,6},
                                         {7,8}}};
 
-    auto result = Oasis::Subtract<Oasis::Expression, Oasis::Expression>{mat2, mat1}.Simplify();
+    auto result = Oasis::Subtract<Oasis::Expression, Oasis::Expression>{mat2, mat1}.Accept(simplifyVisitor).value();
     auto spec = Oasis::RecursiveCast<Oasis::Matrix>(*result);
 
     REQUIRE(spec != nullptr);
@@ -95,8 +98,8 @@ TEST_CASE("Multiply Matrix and Real", "[Matrix][Real][Multiply]")
     Oasis::Matrix mat1{Oasis::MatrixXXD{{1,2},
                                         {3,4}}};
 
-    auto result1 = Oasis::Multiply{a, mat1}.Simplify();
-    auto result2 = Oasis::Multiply{mat1, a}.Simplify();
+    auto result1 = Oasis::Multiply{a, mat1}.Accept(simplifyVisitor).value();
+    auto result2 = Oasis::Multiply{mat1, a}.Accept(simplifyVisitor).value();
     auto spec1 = Oasis::RecursiveCast<Oasis::Matrix>(*result1);
     auto spec2 = Oasis::RecursiveCast<Oasis::Matrix>(*result2);
 
@@ -125,7 +128,7 @@ TEST_CASE("Multiply Matrices, same dimensions", "[Matrix][Multiply]")
     Oasis::Matrix mat1{Oasis::MatrixXXD{{5,6},{7,8}}};
     Oasis::Matrix mat2{Oasis::MatrixXXD{{1,2},{3,4}}};
     Oasis::Matrix expected{Oasis::MatrixXXD{{23,34},{31,46}}};
-    auto result = Oasis::Multiply{mat1, mat2}.Simplify();
+    auto result = Oasis::Multiply{mat1, mat2}.Accept(simplifyVisitor).value();
     REQUIRE(result->Equals(expected));
 }
 
@@ -142,12 +145,12 @@ TEST_CASE("Multiply Matrices, different dimensions", "[Matrix][Multiply]")
     Oasis::Matrix expected23{Oasis::MatrixXXD{{9,12,15},{19,26,33},{29,40,51}}};
     Oasis::Matrix expected32{Oasis::MatrixXXD{{22,28},{49,64}}};
 
-    auto res12 = Oasis::Multiply{mat1, mat2}.Simplify();
-    auto res21 = Oasis::Multiply{mat2, mat1}.Simplify();
-    auto res13 = Oasis::Multiply{mat1, mat3}.Simplify();
-    auto res31 = Oasis::Multiply{mat3, mat1}.Simplify();
-    auto res23 = Oasis::Multiply{mat2, mat3}.Simplify();
-    auto res32 = Oasis::Multiply{mat3, mat2}.Simplify();
+    auto res12 = Oasis::Multiply{mat1, mat2}.Accept(simplifyVisitor).value();
+    auto res21 = Oasis::Multiply{mat2, mat1}.Accept(simplifyVisitor).value();
+    auto res13 = Oasis::Multiply{mat1, mat3}.Accept(simplifyVisitor).value();
+    auto res31 = Oasis::Multiply{mat3, mat1}.Accept(simplifyVisitor).value();
+    auto res23 = Oasis::Multiply{mat2, mat3}.Accept(simplifyVisitor).value();
+    auto res32 = Oasis::Multiply{mat3, mat2}.Accept(simplifyVisitor).value();
 
     REQUIRE(res12->Equals(expected12));
     REQUIRE(res21->Equals(expected21));
diff --git a/tests/MultiplyTests.cpp b/tests/MultiplyTests.cpp
index 5993e0b2..9251797b 100644
--- a/tests/MultiplyTests.cpp
+++ b/tests/MultiplyTests.cpp
@@ -10,8 +10,10 @@
 #include "Oasis/Multiply.hpp"
 #include "Oasis/Real.hpp"
 #include "Oasis/RecursiveCast.hpp"
-
 #include <Oasis/Divide.hpp>
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Multiplication", "[Multiply]")
 {
@@ -22,7 +24,7 @@ TEST_CASE("Multiplication", "[Multiply]")
         Oasis::Real { 3.0 }
     };
 
-    auto simplified = subtract.Simplify();
+    auto simplified = subtract.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -38,7 +40,7 @@ TEST_CASE("Generalized Multiplication", "[Multiply][Generalized]")
         Oasis::Real { 3.0 }
     };
 
-    auto simplified = subtract.Simplify();
+    auto simplified = subtract.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -66,9 +68,9 @@ TEST_CASE("Imaginary Multiplication", "[Imaginary][Multiplication]")
             Oasis::Real { 2.0 } }
     };
 
-    auto simplified2 = i2.Simplify();
-    auto simplified3 = i3.Simplify();
-    auto simplified4 = i4.Simplify();
+    auto simplified2 = i2.Accept(simplifyVisitor).value();
+    auto simplified3 = i3.Accept(simplifyVisitor).value();
+    auto simplified4 = i4.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Multiply { Oasis::Real { -1 }, Oasis::Imaginary {} }.Equals(*simplified3));
     REQUIRE(Oasis::Real { -1.0 }.Equals(*simplified2));
@@ -88,7 +90,7 @@ TEST_CASE("Multiply Associativity", "[Multiply][Associativity]")
                 Oasis::Real { 4.0 } } }
     };
 
-    auto simplified1 = m1.Simplify();
+    auto simplified1 = m1.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Multiply {
         Oasis::Multiply {
@@ -126,7 +128,7 @@ TEST_CASE("Variadic Multiply Constructor", "[Multiply]")
 
     const Oasis::Real expected { 432.0 };
 
-    const auto simplified = multiply.Simplify();
+    const auto simplified = multiply.Accept(simplifyVisitor).value();
     REQUIRE(expected.Equals(*simplified));
 }
 
@@ -145,6 +147,6 @@ TEST_CASE("Multiplying Fractions", "[Multiply]")
 
     const Oasis::Real expected { 3.0 };
 
-    const auto simplified = multiply.Simplify();
+    const auto simplified = multiply.Accept(simplifyVisitor).value();
     REQUIRE(expected.Equals(*simplified));
 }
\ No newline at end of file
diff --git a/tests/NegateTests.cpp b/tests/NegateTests.cpp
index 4e1640ca..488e94b7 100644
--- a/tests/NegateTests.cpp
+++ b/tests/NegateTests.cpp
@@ -9,6 +9,9 @@
 #include "Oasis/Add.hpp"
 #include "Oasis/Variable.hpp"
 #include "Oasis/Real.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Negate", "[Negate]")
 {
@@ -16,7 +19,7 @@ TEST_CASE("Negate", "[Negate]")
         Oasis::Real { 1.0 }
     };
 
-    const auto simplified = negativeOne.Simplify();
+    const auto simplified = negativeOne.Accept(simplifyVisitor).value();
 
     const Oasis::Real expected { -1.0 };
     REQUIRE(simplified->Equals(expected));
diff --git a/tests/SubtractTests.cpp b/tests/SubtractTests.cpp
index cb2612b9..189fea8d 100644
--- a/tests/SubtractTests.cpp
+++ b/tests/SubtractTests.cpp
@@ -11,6 +11,9 @@
 #include "Oasis/Variable.hpp"
 #include "Oasis/Add.hpp"
 #include "Oasis/RecursiveCast.hpp"
+#include "Oasis/SimplifyVisitor.hpp"
+
+inline Oasis::SimplifyVisitor simplifyVisitor{};
 
 TEST_CASE("Subtraction", "[Subtract]")
 {
@@ -21,7 +24,7 @@ TEST_CASE("Subtraction", "[Subtract]")
         Oasis::Real { 3.0 }
     };
 
-    auto simplified = subtract.Simplify();
+    auto simplified = subtract.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -37,7 +40,7 @@ TEST_CASE("Generalized Subtraction", "[Subtract][Generalized]")
         Oasis::Real { 3.0 }
     };
 
-    auto simplified = subtract.Simplify();
+    auto simplified = subtract.Accept(simplifyVisitor).value();
     REQUIRE(simplified->Is<Oasis::Real>());
 
     auto simplifiedReal = dynamic_cast<Oasis::Real&>(*simplified);
@@ -53,7 +56,7 @@ TEST_CASE("Imaginary Subtration", "[Subtract][Imaginary]")
         Oasis::Imaginary {}
     };
 
-    auto spec1 = s1.Simplify();
+    auto spec1 = s1.Accept(simplifyVisitor).value();
 
     REQUIRE(Oasis::Multiply { Oasis::Real { 2.0 }, Oasis::Imaginary {} }.Equals(*spec1));
 }
@@ -83,5 +86,5 @@ TEST_CASE("Simplify Equation with subtraction", "[Subtract]") {
             Oasis::Real{-5.0}
     };
 
-    auto simplified = add1.Simplify();
+    auto simplified = add1.Accept(simplifyVisitor).value();
 }
\ No newline at end of file
diff --git a/tests/UnaryExpressionTests.cpp b/tests/UnaryExpressionTests.cpp
index bca7a1b6..8c5fa746 100644
--- a/tests/UnaryExpressionTests.cpp
+++ b/tests/UnaryExpressionTests.cpp
@@ -18,7 +18,7 @@ TEST_CASE("Substitute Unary", "[Substitute]")
         Oasis::Multiply{
             Oasis::Real { 3.0 },
             Oasis::Variable { "x" } }
-    }; // 2x+3x
+    }; // -2x+3x
 
     const auto after = before.Substitute(Oasis::Variable { "x" }, Oasis::Real { 4.0 }); // after should some std::unique_ptr<Expression> such that it equals 2(-4) + 3(4)
     REQUIRE(after->Equals(Oasis::Real { 4 }));