diff --git a/include/Oasis/Expression.hpp b/include/Oasis/Expression.hpp index c1698801..5bdffdd8 100644 --- a/include/Oasis/Expression.hpp +++ b/include/Oasis/Expression.hpp @@ -1,7 +1,6 @@ #ifndef OASIS_EXPRESSION_HPP #define OASIS_EXPRESSION_HPP -#include #include #include #include @@ -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>; + /** + * 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; + /** * Gets the category of this expression. * @return The category of this expression. diff --git a/src/Expression.cpp b/src/Expression.cpp index e5fb0486..c882a8e8 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -218,6 +218,94 @@ auto Expression::FindZeros() const -> std::vector> return results; } +auto Expression::ApproximateZeros(const Expression& variable, const Expression& guess, int iterations) const -> std::unique_ptr +{ + // 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> 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 originalFunction = this->Copy(); + std::unique_ptr 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 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 evaluatedFunction = originalFunction->Substitute(variable, *x); + std::unique_ptr 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() || !evaluatedDerivative->Is()) { + return originalFunction; + } + + // Divide f'(a) / f(a) + + std::unique_ptr divided = Divide { *evaluatedFunction, *evaluatedDerivative }.Copy(); + + divided = divided->Accept(simplifyVisitor).value(); + + // Subtract x - (f'(a) / f(a)) to get the new guess + x = Subtract { *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; diff --git a/tests/ApproximationTests.cpp b/tests/ApproximationTests.cpp new file mode 100644 index 00000000..43c514d5 --- /dev/null +++ b/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 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 approximationNeg = polynomial.ApproximateZeros(x, guessNeg, 10); + std::unique_ptr 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 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 result = polynomial.ApproximateZeros(x, guess, 5); + + REQUIRE(polynomial.Equals(*result)); +} \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index efef44bd..2f55ea96 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -8,6 +8,7 @@ set(Oasis_TESTS # cmake-format: sortable AddTests.cpp + ApproximationTests.cpp BinaryExpressionTests.cpp Common.hpp DefiniteIntegralTests.cpp