[Feature] Add Approximations of Zeros

include/Oasis/Expression.hpp

@@ -1,7 +1,6 @@
 #ifndef OASIS_EXPRESSION_HPP
 #define OASIS_EXPRESSION_HPP
 
-#include <cstdint>
 #include <expected>
 #include <memory>
 #include <vector>
@@ -88,12 +87,36 @@ class Expression {
     [[nodiscard]] virtual auto Equals(const Expression& other) const -> bool = 0;
 
     /**
-     * The FindZeros function finds all rational real zeros, and up to 2 irrational/complex zeros of a polynomial. Currently assumes an expression of the form a+bx+cx^2+dx^3+... where a, b, c, d are a integers.
+     * Finds the exact zeros of a polynomial when they can be derived symbolically.
      *
-     * @tparam origonalExpresion The expression for which all the factors will be found.
+     * This method is intended for polynomial expressions with integer coefficients. It can return
+     * rational roots exactly, and may also return roots from degree-1 and degree-2 cases using the
+     * usual closed-form formulas.
+     *
+     * Unlike `ApproximateZeros`, this does not iterate from a guess or produce a numeric estimate.
+     * If the expression is not in a supported polynomial form, an empty result is returned.
+     *
+     * @return A vector of expressions representing the zeros that were found exactly.
      */
     auto FindZeros() const -> std::vector<std::unique_ptr<Expression>>;
 
+    /**
+     * Approximates a zero of this expression using Newton's method.
+     *
+     * This method performs an iterative numerical search for a single root of the expression with
+     * respect to the supplied variable. The initial guess determines which root is targeted, and
+     * more iterations generally produce a better approximation.
+     *
+     * If the expression cannot be differentiated, evaluated, or converged safely from the supplied
+     * guess, the original expression is returned instead.
+     *
+     * @param variable The variable to differentiate and substitute with respect to.
+     * @param guess The starting point for the iteration.
+     * @param iterations The number of Newton iterations to perform.
+     * @return A root approximation, or a copy of the original expression if no approximation is found.
+     */
+    [[nodiscard]] auto ApproximateZeros(const Expression& variable, const Expression& guess, int iterations) const -> std::unique_ptr<Expression>;
+
     /**
      * Gets the category of this expression.
      * @return The category of this expression.

src/Expression.cpp

@@ -218,6 +218,94 @@ auto Expression::FindZeros() const -> std::vector<std::unique_ptr<Expression>>
     return results;
 }
 
+auto Expression::ApproximateZeros(const Expression& variable, const Expression& guess, int iterations) const -> std::unique_ptr<Expression>
+{
+    // Setup simplifyVisitor for later
+    SimplifyVisitor simplifyVisitor {};
+
+    // Helper variables
+
+    // List of guesses that we've taken.
+    // If any of these are duplicates, then we've run into a cycle and can't find a root.
+    std::vector<std::unique_ptr<Expression>> guessList;
+
+    // Variable representing the number 0.
+    // Useful in determining a situation where the guess a results in f(a) = 0 or f'(a) = 0
+    Real zero { 0.0f };
+
+    // Declare the function and its derivative
+    std::unique_ptr<Expression> originalFunction = this->Copy();
+    std::unique_ptr<Expression> derivative = originalFunction->Differentiate(variable);
+
+    // Function failed to differentiate, return the original function
+    if (derivative == nullptr)
+        return originalFunction;
+
+    // New guess (starts at the original guess's value)
+    std::unique_ptr<Expression> x = guess.Copy();
+
+    // Iterate for the given number of times
+    // The higher the number of iterations, the more accurate the result (in theory)
+    for (int i = 0; i < iterations; ++i) {
+        // Evaluated version of the functions
+        // f(a) and f'(a), for the guess a (stored in x)
+        std::unique_ptr<Expression> evaluatedFunction = originalFunction->Substitute(variable, *x);
+        std::unique_ptr<Expression> evaluatedDerivative = derivative->Substitute(variable, *x);
+
+        evaluatedFunction = evaluatedFunction->Accept(simplifyVisitor).value();
+        evaluatedDerivative = evaluatedDerivative->Accept(simplifyVisitor).value();
+
+        // If this is true, then we are in a divide-by-0 case
+        // (when evaluated_function = 0). We can't continue, so return the original function.
+        if (evaluatedDerivative->Equals(zero))
+            return originalFunction;
+
+        // If evaluated_function = 0 (and evaluated_derivative != 0),
+        // then we found an approximation. Return that value.
+        if (evaluatedFunction->Equals(zero)) {
+            // If our current approximation is just the original guess,
+            // then we weren't able to make a new guess (either we went in
+            // a loop, or f(guess) = 0). In that case, we didn't find
+            // anything useful, so we should return the original function instead.
+            if (x->Equals(guess))
+                return originalFunction;
+
+            // Otherwise, we found an actual approximation, so return that.
+            return x;
+        }
+
+        // We might not be able to evaluate this function as a number.
+        // In that case, we can't get to a number, so we need to return original function.
+        if (!evaluatedFunction->Is<Real>() || !evaluatedDerivative->Is<Real>()) {
+            return originalFunction;
+        }
+
+        // Divide f'(a) /  f(a)
+
+        std::unique_ptr<Expression> divided = Divide<Expression, Expression> { *evaluatedFunction, *evaluatedDerivative }.Copy();
+
+        divided = divided->Accept(simplifyVisitor).value();
+
+        // Subtract x - (f'(a) / f(a)) to get the new guess
+        x = Subtract<Expression, Expression> { *x, *divided }.Accept(simplifyVisitor).value();
+
+        // Add this to our list of guesses
+        guessList.push_back(x->Copy());
+
+        // Make sure that we haven't already seen this value.
+        // If we have, then we've gone into a cycle and won't be able to find a solution.
+        for (int j = 0; j < i; ++j) {
+            if (guessList[i]->Equals(*guessList[j])) {
+                // We found a duplicate and will be going in a loop
+                // We can't find anything from here, so we have to return the original function.
+                return originalFunction;
+            }
+        }
+    }
+
+    return x;
+}
+
 auto Expression::GetCategory() const -> uint32_t
 {
     return 0;

tests/ApproximationTests.cpp

@@ -0,0 +1,112 @@
+//
+// Created by Justin Romanelli on 3/20/26.
+//
+
+#include "Oasis/Add.hpp"
+#include "Oasis/Exponent.hpp"
+#include "Oasis/Expression.hpp"
+#include "Oasis/Multiply.hpp"
+#include "Oasis/Real.hpp"
+#include "Oasis/Subtract.hpp"
+#include "Oasis/Variable.hpp"
+#include "catch2/catch_test_macros.hpp"
+
+TEST_CASE("Approximation of a Linear function", "[Real][Approximation]")
+{
+    // Approximating f(x) = 5x - 10, looking for the root x = 2
+    Oasis::Variable x { "x" };
+    Oasis::Add<> linear {
+        Oasis::Multiply<> { Oasis::Real { 5.0f }, x },
+        Oasis::Real { -10.0f }
+    };
+    Oasis::Real guess { 5.0f }, ans { 2.0f };
+
+    // Approximate a root
+    std::unique_ptr<Oasis::Expression> result = linear.ApproximateZeros(x, guess, 10);
+
+    // Make sure the approximation worked correctly
+
+    // After a certain level of precision, the approximation becomes "exact", in that
+    // it truncates the decimal to the regular value if the result is precise enough.
+    // So, the approximations should be equal to the roots themselves.
+    REQUIRE(ans.Equals(*result));
+}
+
+TEST_CASE("Approximation of a Polynomial", "[Real][Approximation]")
+{
+    // Approximate a root for f(x) = x^2 - 4
+    // Looking for both roots, x = 2 and x = -2
+    Oasis::SimplifyVisitor sV {};
+    Oasis::Variable x { "x" };
+    Oasis::Subtract<> polynomial {
+        Oasis::Exponent<> {
+            x, Oasis::Real { 2.0f } },
+        Oasis::Real { 4.0f }
+    };
+
+    // Actual roots to the polynomial
+    Oasis::Real rootNeg { -2.0f }, rootPos { 2.0f };
+
+    // Guesses to start the approximations
+    Oasis::Real guessNeg { -5.0f }, guessPos { 5.0f };
+
+    // Approximations of the roots
+    // Make sure these are equal to the original roots
+    std::unique_ptr<Oasis::Expression> approximationNeg = polynomial.ApproximateZeros(x, guessNeg, 10);
+    std::unique_ptr<Oasis::Expression> approximationPos = polynomial.ApproximateZeros(x, guessPos, 10);
+
+    // Simplify the approximations
+    approximationNeg = *approximationNeg->Accept(sV);
+    approximationPos = *approximationPos->Accept(sV);
+
+    // Check that they equal the actual roots (root1Ans and root2Ans)
+    REQUIRE(rootNeg.Equals(*approximationNeg));
+    REQUIRE(rootPos.Equals(*approximationPos));
+}
+
+TEST_CASE("Approximation of a Higher-Degree Polynomial", "[Real][Approximation]")
+{
+    // Testing and setup variables
+    Oasis::Variable x { "x" };
+
+    // Approximating 256x^10 - 16x^6 = 0
+    // One root is x = 0.5 (looking for it here)
+    Oasis::Subtract<> polynomial {
+        Oasis::Multiply<> {
+            Oasis::Real { 256.0f },
+            Oasis::Exponent<> {
+                x, Oasis::Real { 10.0f } } },
+        Oasis::Multiply<> {
+            Oasis::Real { 16.0f },
+            Oasis::Exponent<> {
+                x, Oasis::Real { 6.0f } } }
+    };
+
+    // Initial guess and answer
+    Oasis::Real guess { 3.0f }, answer { 0.5f };
+
+    // Approximate the answer here
+    std::unique_ptr<Oasis::Expression> result = polynomial.ApproximateZeros(x, guess, 25);
+
+    // Check to make sure it's equal to the actual answer (with many iterations)
+    REQUIRE(answer.Equals(*result));
+}
+
+TEST_CASE("Impossible Approximation of a Polynomial (Stuck)", "[Approximation]")
+{
+    // Testing variables
+    // Attempt to approximate x^2 - 1 = 0
+    Oasis::Variable x { "x" };
+    Oasis::Subtract<> polynomial {
+        Oasis::Exponent<> {
+            x, Oasis::Real { 2.0f } },
+        Oasis::Real { 1.0f }
+    };
+    Oasis::Real guess { 0.0f };
+
+    // This should get "stuck" at x = 0, such that it can't find a better value to use
+    // So, it should return the original function. Make sure that's the case.
+    std::unique_ptr<Oasis::Expression> result = polynomial.ApproximateZeros(x, guess, 5);
+
+    REQUIRE(polynomial.Equals(*result));
+}

tests/CMakeLists.txt

@@ -8,6 +8,7 @@
 set(Oasis_TESTS
     # cmake-format: sortable
     AddTests.cpp
+    ApproximationTests.cpp
     BinaryExpressionTests.cpp
     Common.hpp
     DefiniteIntegralTests.cpp