From 322f3f86f3d0a3dbaf5c4289f687968b76d7aa3c Mon Sep 17 00:00:00 2001 From: Matthew McCall Date: Fri, 31 Jan 2025 17:28:03 -0500 Subject: [PATCH 01/21] Reorganize test build logic and enhance InFixSerializer Moved test build configuration to improve modularity and added support in InFixSerializer for new types (e.g., Matrix, EulerNumber, Pi, Magnitude). Introduced `Common.hpp` for shared test macros, enabling serializer-based assertions. Adjusted tests accordingly to utilize the new serializer functionality. --- CMakeLists.txt | 8 ++++---- io/include/Oasis/InFixSerializer.hpp | 4 ++++ io/src/InFixSerializer.cpp | 22 ++++++++++++++++++++++ tests/AddTests.cpp | 8 ++++++-- tests/CMakeLists.txt | 10 ++++++++-- tests/Common.hpp | 17 +++++++++++++++++ 6 files changed, 61 insertions(+), 8 deletions(-) create mode 100644 tests/Common.hpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 0c98d6ce..40292bd0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -45,10 +45,10 @@ include(cmake/FetchTinyxml2.cmake) add_subdirectory(include) add_subdirectory(src) -if(OASIS_BUILD_TESTS) - add_subdirectory(tests) -endif() - if(OASIS_BUILD_IO) add_subdirectory(io) endif() + +if(OASIS_BUILD_TESTS) + add_subdirectory(tests) +endif() diff --git a/io/include/Oasis/InFixSerializer.hpp b/io/include/Oasis/InFixSerializer.hpp index 0e12852a..fc3ef36d 100644 --- a/io/include/Oasis/InFixSerializer.hpp +++ b/io/include/Oasis/InFixSerializer.hpp @@ -26,6 +26,10 @@ class InFixSerializer final : public Visitor { void Visit(const Negate& negate) override; void Visit(const Derivative& derivative) override; void Visit(const Integral& integral) override; + void Visit(const Matrix& matrix) override; + void Visit(const EulerNumber&) override; + void Visit(const Pi&) override; + void Visit(const Magnitude& magnitude) override; [[nodiscard]] std::string getResult() const; diff --git a/io/src/InFixSerializer.cpp b/io/src/InFixSerializer.cpp index c9ed15f9..47412ea4 100644 --- a/io/src/InFixSerializer.cpp +++ b/io/src/InFixSerializer.cpp @@ -12,6 +12,7 @@ #include "Oasis/Exponent.hpp" #include "Oasis/Integral.hpp" #include "Oasis/Log.hpp" +#include "Oasis/Magnitude.hpp" #include "Oasis/Multiply.hpp" #include "Oasis/Negate.hpp" #include "Oasis/Real.hpp" @@ -136,6 +137,27 @@ void InFixSerializer::Visit(const Integral<>& integral) result = fmt::format("in({},{})", mostSigOpStr, leastSigOpStr); } +void InFixSerializer::Visit(const Matrix& matrix) +{ + result = ""; +} + +void InFixSerializer::Visit(const EulerNumber&) +{ + result = "e"; +} + +void InFixSerializer::Visit(const Pi&) +{ + result = "pi"; +} + +void InFixSerializer::Visit(const Magnitude& magnitude) +{ + magnitude.GetOperand().Accept(*this); + result = fmt::format("|{}|", getResult()); +} + std::string InFixSerializer::getResult() const { return result; diff --git a/tests/AddTests.cpp b/tests/AddTests.cpp index 9e805911..88b2cdf0 100644 --- a/tests/AddTests.cpp +++ b/tests/AddTests.cpp @@ -1,6 +1,8 @@ // // Created by Matthew McCall on 8/7/23. // +#include "Common.hpp" + #include "catch2/catch_test_macros.hpp" #include "Oasis/Add.hpp" @@ -10,8 +12,7 @@ #include "Oasis/Real.hpp" #include "Oasis/RecursiveCast.hpp" #include "Oasis/Variable.hpp" - -#include +#include "Oasis/InFixSerializer.hpp" TEST_CASE("Addition", "[Add]") { @@ -22,6 +23,9 @@ TEST_CASE("Addition", "[Add]") Oasis::Real { 3.0 } }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + auto simplified = add.Simplify(); REQUIRE(simplified->Is()); diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index c204c769..98901002 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -21,7 +21,8 @@ set(Oasis_TESTS NegateTests.cpp PolynomialTests.cpp SubtractTests.cpp - UnaryExpressionTests.cpp) + UnaryExpressionTests.cpp + Common.hpp) # Adds an executable target called "OasisTests" to be built from sources files. add_executable(OasisTests ${Oasis_TESTS}) @@ -35,7 +36,12 @@ endif() target_link_libraries(OasisTests PRIVATE Oasis::Oasis Catch2::Catch2WithMain) -# Automatically registers the tests with CTest. +if (OASIS_BUILD_IO) + target_link_libraries(OasisTests PRIVATE Oasis::IO) + target_compile_definitions(OasisTests PRIVATE OASIS_IO_ENABLED) +endif () + + # Automatically registers the tests with CTest. list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras) include(Catch) catch_discover_tests(OasisTests) diff --git a/tests/Common.hpp b/tests/Common.hpp new file mode 100644 index 00000000..af6b59c1 --- /dev/null +++ b/tests/Common.hpp @@ -0,0 +1,17 @@ +// +// Created by Matthew McCall on 1/31/25. +// + +#ifndef COMMON_HPP +#define COMMON_HPP + +#ifdef OASIS_IO_ENABLED +#define OASIS_CAPTURE_WITH_SERIALIZER(serializer, __expr) \ +__expr.Accept(serializer); \ +auto __expr_str = serializer.getResult(); \ +INFO(#__expr << " := " << __expr_str); +#else +#define OASIS_CAPTURE_WITH_SERIALIZER(serializer, __expr) +#endif + +#endif //COMMON_HPP From 7ee1db8439d8ae15eaae04fafb3a7c05929e03a2 Mon Sep 17 00:00:00 2001 From: RC <35265825+richcfno1@users.noreply.github.com> Date: Fri, 31 Jan 2025 21:46:30 -0500 Subject: [PATCH 02/21] test cases of the polynomial --- tests/PolynomialTests.cpp | 86 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 86 insertions(+) diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 7b7862bc..ccd042cf 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -18,6 +18,7 @@ #include #include #include +#include // TODO: Figure out what's going out here // TEST_CASE("7th degree polynomial with rational roots", "[factor][duplicateRoot]") @@ -124,3 +125,88 @@ TEST_CASE("linear polynomial", "[factor]") REQUIRE(root->GetLeastSigOp().GetValue() == 1); } } + +TEST_CASE("quadratic polynomial", "[factor]") +{ + // x² + 5x + 6 + Oasis::Add<> add { + Oasis::Real(6), // constant + Oasis::Multiply{ // 5x + Oasis::Real(5), + Oasis::Variable("x") + }, + Oasis::Exponent{ // x² + Oasis::Variable("x"), + Oasis::Real(2) + } + }; + + auto zeros = add.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + if (zeros.size() == 2) + { + std::cout << "the result has the correct size.\n"; + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -2" << std::endl; + std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + if (root2 != nullptr) + { + std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; + std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + } + } + + REQUIRE(zeros.size() == 2); + + if (zeros.size() == 2) { + // Check first root (-2) + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == -2); + REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + + // Check second root (-3) + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -3); + REQUIRE(root2->GetLeastSigOp().GetValue() == 1); + } +} + +TEST_CASE("simple quadra polynomial", "[factor]") +{ + // x^2 - 9 + Oasis::Subtract minus { + Oasis::Exponent { // x² + Oasis::Variable("x"), + Oasis::Real(2) }, + Oasis::Real(9), + }; + auto zeros = minus.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + if (zeros.size() == 2) + { + std::cout << "the result has the correct size.\n"; + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -2" << std::endl; + std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + if (root2 != nullptr) + { + std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; + std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + } + } + REQUIRE(zeros.size() == 2); + if (zeros.size() == 2) { + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == -3); + REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == 3); + REQUIRE(root2->GetLeastSigOp().GetValue() == 1); + } +} \ No newline at end of file From ea4ae22c0efcae618a42d7b79fc823d93ddaf883 Mon Sep 17 00:00:00 2001 From: RC <35265825+richcfno1@users.noreply.github.com> Date: Fri, 7 Feb 2025 17:06:18 -0500 Subject: [PATCH 03/21] commit test --- src/Expression.cpp | 59 ++++++++++++++++++++++--- tests/PolynomialTests.cpp | 92 +++++++++++++++++++++++++++++++++++++++ 2 files changed, 145 insertions(+), 6 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 2d863e33..5f52b50d 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -52,6 +52,40 @@ auto Expression::FindZeros() const -> std::vector> { std::vector> results; std::vector> termsE; + if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { + // Check for x² - n pattern + if (auto leftTerm = RecursiveCast>(subCase->GetMostSigOp()); + leftTerm != nullptr) { + + if (leftTerm->GetLeastSigOp().GetValue() == 2) { + if (auto rightTerm = RecursiveCast(subCase->GetLeastSigOp()); + rightTerm != nullptr) { + + double n = rightTerm->GetValue(); + double sqrtN = std::sqrt(n); + + if (n > 0 && sqrtN == std::floor(sqrtN)) { + // Instead of creating Real values directly, create proper Binary Expressions + // For positive root: n^(1/2) + auto posRoot = std::make_unique>( + Real(1), // MostSigOp + Real(sqrtN) // LeastSigOp + ); + + // For negative root: -n^(1/2) + auto negRoot = std::make_unique>( + Real(-1), // MostSigOp + Real(sqrtN) // LeastSigOp + ); + + results.push_back(std::move(posRoot)); + results.push_back(std::move(negRoot)); + return results; + } + } + } + } + } if (auto addCase = RecursiveCast>(*this); addCase != nullptr) { addCase->Flatten(termsE); } else { @@ -196,15 +230,28 @@ auto Expression::FindZeros() const -> std::vector> } if (coefficents.size() == 2) { results.push_back(Divide(Multiply(Real(-1), *coefficents[0]), *coefficents[1]).Simplify()); - } else if (coefficents.size() == 3) { + } + if (coefficents.size() == 3) { auto& a = coefficents[2]; auto& b = coefficents[1]; auto& c = coefficents[0]; - auto negB = Multiply(Real(-1.0), *b).Simplify(); - auto sqrt = Exponent(*Add(Multiply(*b, *b), Multiply(Real(-4), Multiply(*a, *c))).Simplify(), Divide(Real(1), Real(2))).Copy(); - auto twoA = Multiply(Real(2), *a).Simplify(); - results.push_back(Divide(Add(*negB, *sqrt), *twoA).Copy()); - results.push_back(Divide(Subtract(*negB, *sqrt), *twoA).Copy()); + + // Calculate discriminant first + auto discriminant = Add(Multiply(*b, *b), + Multiply(Real(-4), Multiply(*a, *c))).Simplify(); + + // Only proceed if discriminant is non-negative + if (auto realDisc = RecursiveCast(*discriminant); + realDisc != nullptr && realDisc->GetValue() >= 0) { + + auto negB = Multiply(Real(-1.0), *b).Simplify(); + auto sqrt = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + auto twoA = Multiply(Real(2), *a).Simplify(); + + // Create both roots + results.push_back(Divide(Add(*negB, *sqrt), *twoA).Copy()); + results.push_back(Divide(Subtract(*negB, *sqrt), *twoA).Copy()); + } } return results; } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 7b7862bc..cdf6b1f0 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -3,6 +3,7 @@ // #include "catch2/catch_test_macros.hpp" +#include "Common.hpp" #include "Oasis/Add.hpp" #include "Oasis/Divide.hpp" #include "Oasis/Exponent.hpp" @@ -11,6 +12,7 @@ #include "Oasis/Multiply.hpp" #include "Oasis/Real.hpp" #include "Oasis/RecursiveCast.hpp" +#include "Oasis/InFixSerializer.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" @@ -18,6 +20,7 @@ #include #include #include +#include // TODO: Figure out what's going out here // TEST_CASE("7th degree polynomial with rational roots", "[factor][duplicateRoot]") @@ -124,3 +127,92 @@ TEST_CASE("linear polynomial", "[factor]") REQUIRE(root->GetLeastSigOp().GetValue() == 1); } } + +TEST_CASE("quadratic polynomial", "[factor]") +{ + // x² + 5x + 6 + Oasis::Add<> add { + Oasis::Real(6), // constant + Oasis::Multiply{ // 5x + Oasis::Real(5), + Oasis::Variable("x") + }, + Oasis::Exponent{ // x² + Oasis::Variable("x"), + Oasis::Real(2) + } + }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + + auto zeros = add.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + if (zeros.size() == 2) + { + std::cout << "the result has the correct size.\n"; + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -2" << std::endl; + std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + if (root2 != nullptr) + { + std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; + std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + } + } + + REQUIRE(zeros.size() == 2); + + if (zeros.size() == 2) { + // Check first root (-2) + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == -2); + REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + + // Check second root (-3) + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -3); + REQUIRE(root2->GetLeastSigOp().GetValue() == 1); + } +} + +TEST_CASE("simple quadra polynomial", "[factor]") +{ + // x^2 - 9 + Oasis::Subtract minus { + Oasis::Exponent { // x² + Oasis::Variable("x"), + Oasis::Real(2) }, + Oasis::Real(9), + }; + auto zeros = minus.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + if (zeros.size() == 2) + { + std::cout << "the result has the correct size.\n"; + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -3" << std::endl; + std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + if (root2 != nullptr) + { + std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; + std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + } + } + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + REQUIRE(zeros.size() == 2); + if (zeros.size() == 2) { + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == -3); + REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == 3); + REQUIRE(root2->GetLeastSigOp().GetValue() == 1); + } +} \ No newline at end of file From 8cf2495b70f809fcb523fabeef2ef640c1c950e0 Mon Sep 17 00:00:00 2001 From: RC <35265825+richcfno1@users.noreply.github.com> Date: Fri, 7 Feb 2025 18:51:20 -0500 Subject: [PATCH 04/21] now it can correctly handle expression like x^2 - n, where sqrt(n) is a perfect number. --- include/Oasis/Expression.hpp | 9 +++ src/Expression.cpp | 110 ++++++++++++++++++++++++++++++----- tests/PolynomialTests.cpp | 70 ++++++++++++++++------ 3 files changed, 158 insertions(+), 31 deletions(-) diff --git a/include/Oasis/Expression.hpp b/include/Oasis/Expression.hpp index bb905049..7b50f456 100644 --- a/include/Oasis/Expression.hpp +++ b/include/Oasis/Expression.hpp @@ -88,6 +88,15 @@ class Expression { */ auto FindZeros() const -> std::vector>; + /** + * Finds the real number roots of a polynomial expression. + * This function handles special cases like x² - n (where n is a perfect square) + * and general quadratic equations (ax² + bx + c). + * + * @return A vector of Real values representing the real roots of the polynomial + */ + [[nodiscard]] auto Polynomial_Real() const -> std::vector>; + /** * Gets the category of this expression. * @return The category of this expression. diff --git a/src/Expression.cpp b/src/Expression.cpp index 5f52b50d..bf01bb0a 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -68,14 +68,14 @@ auto Expression::FindZeros() const -> std::vector> // Instead of creating Real values directly, create proper Binary Expressions // For positive root: n^(1/2) auto posRoot = std::make_unique>( - Real(1), // MostSigOp - Real(sqrtN) // LeastSigOp + Real(sqrtN), // MostSigOp + Real(1) // LeastSigOp ); // For negative root: -n^(1/2) auto negRoot = std::make_unique>( - Real(-1), // MostSigOp - Real(sqrtN) // LeastSigOp + Real(-sqrtN), // MostSigOp + Real(1) // LeastSigOp ); results.push_back(std::move(posRoot)); @@ -236,23 +236,107 @@ auto Expression::FindZeros() const -> std::vector> auto& b = coefficents[1]; auto& c = coefficents[0]; - // Calculate discriminant first - auto discriminant = Add(Multiply(*b, *b), - Multiply(Real(-4), Multiply(*a, *c))).Simplify(); + // Calculate discriminant + auto bSquared = Multiply(*b, *b).Simplify(); + auto fourAC = Multiply(Real(4), Multiply(*a, *c)).Simplify(); + auto discriminant = Subtract(*bSquared, *fourAC).Simplify(); - // Only proceed if discriminant is non-negative if (auto realDisc = RecursiveCast(*discriminant); realDisc != nullptr && realDisc->GetValue() >= 0) { - auto negB = Multiply(Real(-1.0), *b).Simplify(); - auto sqrt = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + auto negB = Multiply(Real(-1), *b).Simplify(); + auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); auto twoA = Multiply(Real(2), *a).Simplify(); - // Create both roots - results.push_back(Divide(Add(*negB, *sqrt), *twoA).Copy()); - results.push_back(Divide(Subtract(*negB, *sqrt), *twoA).Copy()); + // First, create the numerators for both roots + auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); + auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); + + // Now create the Divide expressions properly + results.push_back(Divide(*numerator1, *twoA).Copy()); + results.push_back(Divide(*numerator2, *twoA).Copy()); + } + } + return results; +} + +auto Expression::Polynomial_Real() const -> std::vector> { + std::vector> results; + std::vector> termsE; + + // First handle x² - n special case + if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { + if (auto leftTerm = RecursiveCast>(subCase->GetMostSigOp()); + leftTerm != nullptr) { + if (leftTerm->GetLeastSigOp().GetValue() == 2) { + if (auto rightTerm = RecursiveCast(subCase->GetLeastSigOp()); + rightTerm != nullptr) { + double n = rightTerm->GetValue(); + double sqrtN = std::sqrt(n); + if (n > 0 && sqrtN == std::floor(sqrtN)) { + results.push_back(Divide(Real(sqrtN), Real(1)).Copy()); + results.push_back(Divide(Real(-sqrtN), Real(1)).Copy()); + return results; + } + } + } + } + } + + // Process terms + if (auto addCase = RecursiveCast>(*this); addCase != nullptr) { + addCase->Flatten(termsE); + } else { + termsE.push_back(Copy()); + } + + // Collect coefficients + std::vector> coefficients; + for (const auto& term : termsE) { + if (auto realTerm = RecursiveCast(*term); realTerm != nullptr) { + coefficients.push_back(realTerm->Copy()); } } + + // Handle linear equation (ax + b) + if (coefficients.size() == 2) { + if (auto aReal = RecursiveCast(*coefficients[1]); aReal != nullptr) { + if (auto bReal = RecursiveCast(*coefficients[0]); bReal != nullptr) { + double a = aReal->GetValue(); + double b = bReal->GetValue(); + + // ax + b = 0 → x = -b/a + if (a != 0) { + double root = -b / a; + results.push_back(Divide(Real(root), Real(1)).Copy()); + } + } + } + } + + // Handle quadratic equation (for ax² + bx + c) + if (coefficients.size() == 3) { + if (auto aReal = RecursiveCast(*coefficients[2]); aReal != nullptr) { + if (auto bReal = RecursiveCast(*coefficients[1]); bReal != nullptr) { + if (auto cReal = RecursiveCast(*coefficients[0]); cReal != nullptr) { + double a = aReal->GetValue(); + double b = bReal->GetValue(); + double c = cReal->GetValue(); + + double discriminant = b * b - 4 * a * c; + if (discriminant >= 0) { + double sqrtDisc = std::sqrt(discriminant); + double root1 = (-b + sqrtDisc) / (2 * a); + double root2 = (-b - sqrtDisc) / (2 * a); + + results.push_back(Divide(Real(root1), Real(1)).Copy()); + results.push_back(Divide(Real(root2), Real(1)).Copy()); + } + } + } + } + } + return results; } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index cdf6b1f0..8388a73c 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -3,7 +3,6 @@ // #include "catch2/catch_test_macros.hpp" -#include "Common.hpp" #include "Oasis/Add.hpp" #include "Oasis/Divide.hpp" #include "Oasis/Exponent.hpp" @@ -15,6 +14,7 @@ #include "Oasis/InFixSerializer.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" +#include "Common.hpp" #include #include @@ -128,6 +128,40 @@ TEST_CASE("linear polynomial", "[factor]") } } +TEST_CASE("linear polynomial(Real)", "[factor]") +{ + Oasis::Add add { + Oasis::Real(30), + Oasis::Variable("x") + }; + auto result = add.Polynomial_Real(); + REQUIRE(result.size() == 1); + REQUIRE(result[0]->Is()); + + auto poly_var_res = dynamic_cast(*result[0]); + REQUIRE(poly_var_res.GetValue() == -30); +} + +TEST_CASE("linear polynomial(2x +30)", "[factor]") +{ + Oasis::Add add { + Oasis::Real(30), + Oasis::Multiply{ // 5x + Oasis::Real(2), + Oasis::Variable("x") + } + }; + auto zeros = add.FindZeros(); + REQUIRE(zeros.size() == 1); + if (zeros.size() == 1) { + + auto root = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root != nullptr); + REQUIRE(root->GetMostSigOp().GetValue() == -15); + REQUIRE(root->GetLeastSigOp().GetValue() == 1); + } +} + TEST_CASE("quadratic polynomial", "[factor]") { // x² + 5x + 6 @@ -180,39 +214,39 @@ TEST_CASE("quadratic polynomial", "[factor]") TEST_CASE("simple quadra polynomial", "[factor]") { - // x^2 - 9 + // x^2 - 2500 Oasis::Subtract minus { Oasis::Exponent { // x² Oasis::Variable("x"), Oasis::Real(2) }, - Oasis::Real(9), + Oasis::Real(2500), }; auto zeros = minus.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; - if (zeros.size() == 2) - { - std::cout << "the result has the correct size.\n"; - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -3" << std::endl; - std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - if (root2 != nullptr) - { - std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; - std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; - } - } +// if (zeros.size() == 2) +// { +// std::cout << "the result has the correct size.\n"; +// auto root1 = Oasis::RecursiveCast>(*zeros[0]); +// std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be 3" << std::endl; +// std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; +// auto root2 = Oasis::RecursiveCast>(*zeros[1]); +// if (root2 != nullptr) +// { +// std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; +// std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; +// } +// } Oasis::InFixSerializer serializer; OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == -3); + REQUIRE(root1->GetMostSigOp().GetValue() == 50); REQUIRE(root1->GetLeastSigOp().GetValue() == 1); auto root2 = Oasis::RecursiveCast>(*zeros[1]); REQUIRE(root2 != nullptr); - REQUIRE(root2->GetMostSigOp().GetValue() == 3); + REQUIRE(root2->GetMostSigOp().GetValue() == -50); REQUIRE(root2->GetLeastSigOp().GetValue() == 1); } } \ No newline at end of file From db2ad43795f445e41cfb8d9e727adb338b04a92d Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 14 Feb 2025 22:14:22 -0500 Subject: [PATCH 05/21] linear and quadratic --- src/Expression.cpp | 169 ++++++++++++++++++---------------- tests/PolynomialTests.cpp | 184 ++++++++++++++++++++++++++++---------- 2 files changed, 228 insertions(+), 125 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index bf01bb0a..10430d28 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -8,6 +8,7 @@ #include #include #include +#include std::vector getAllFactors(long long n) { @@ -50,13 +51,13 @@ namespace Oasis { */ auto Expression::FindZeros() const -> std::vector> { + std::cout << "start of the function" << std::endl; std::vector> results; std::vector> termsE; if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { // Check for x² - n pattern if (auto leftTerm = RecursiveCast>(subCase->GetMostSigOp()); leftTerm != nullptr) { - if (leftTerm->GetLeastSigOp().GetValue() == 2) { if (auto rightTerm = RecursiveCast(subCase->GetLeastSigOp()); rightTerm != nullptr) { @@ -88,6 +89,12 @@ auto Expression::FindZeros() const -> std::vector> } if (auto addCase = RecursiveCast>(*this); addCase != nullptr) { addCase->Flatten(termsE); + } else if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { + // Handle subtraction by converting to addition + termsE.push_back(subCase->GetMostSigOp().Copy()); // First term + // Add negative of second term + auto negTerm = Multiply(Real(-1), subCase->GetLeastSigOp()).Copy(); + termsE.push_back(std::move(negTerm)); } else { termsE.push_back(Copy()); } @@ -153,6 +160,8 @@ auto Expression::FindZeros() const -> std::vector> while (posCoefficents.size() > 0 && RecursiveCast(*posCoefficents.back()) != nullptr && RecursiveCast(*posCoefficents.back())->GetValue() == 0) { posCoefficents.pop_back(); } + std::cout << "posCoefficents size: " << posCoefficents.size() << std::endl; + std::cout << "negCoefficents size: " << negCoefficents.size() << std::endl; std::vector> coefficents; for (size_t i = negCoefficents.size(); i > 1; i--) { coefficents.push_back(negCoefficents[i - 1]->Simplify()); @@ -160,6 +169,7 @@ auto Expression::FindZeros() const -> std::vector> for (const std::unique_ptr& i : posCoefficents) { coefficents.push_back(i->Simplify()); } + std::cout << "Final coefficients size: " << coefficents.size() << std::endl; if (coefficents.size() <= 1) { return {}; } @@ -176,85 +186,92 @@ auto Expression::FindZeros() const -> std::vector> termsC.push_back(lround(value)); } } +// if (termsC.size() == coefficents.size()) { +// std::reverse(termsC.begin(), termsC.end()); +// for (auto pv : getAllFactors(termsC.back())) { +// for (auto qv : getAllFactors(termsC.front())) { +// if (gcf(pv, qv) == 1) { +// for (long long sign : { -1, 1 }) { +// bool doAdd = true; +// while (true) { +// long long mpv = pv * sign; +// std::vector newTermsC; +// long long h = 0; +// for (long long i : termsC) { +// h *= mpv; +// if (h % qv != 0) { +// break; +// } +// h /= qv; +// h += i; +// newTermsC.push_back(h); +// } +// if (newTermsC.size() == termsC.size() && newTermsC.back() == 0) { +// termsC = newTermsC; +// if (doAdd) { +// results.push_back(std::make_unique>(Real(1.0 * mpv), Real(1.0 * qv))); +// doAdd = false; +// } +// do { +// termsC.pop_back(); +// } while (termsC.back() == 0); +// if (termsC.size() <= 1) { +// break; +// } +// } else { +// break; +// } +// } +// } +// } +// if (termsC.size() <= 1) { +// break; +// } +// } +// if (termsC.size() <= 1) { +// break; +// } +// } +// coefficents.clear(); +// std::reverse(termsC.begin(), termsC.end()); +// for (auto i : termsC) { +// coefficents.push_back(Real(i * 1.0).Copy()); +// } +// } if (termsC.size() == coefficents.size()) { - std::reverse(termsC.begin(), termsC.end()); - for (auto pv : getAllFactors(termsC.back())) { - for (auto qv : getAllFactors(termsC.front())) { - if (gcf(pv, qv) == 1) { - for (long long sign : { -1, 1 }) { - bool doAdd = true; - while (true) { - long long mpv = pv * sign; - std::vector newTermsC; - long long h = 0; - for (long long i : termsC) { - h *= mpv; - if (h % qv != 0) { - break; - } - h /= qv; - h += i; - newTermsC.push_back(h); - } - if (newTermsC.size() == termsC.size() && newTermsC.back() == 0) { - termsC = newTermsC; - if (doAdd) { - results.push_back(std::make_unique>(Real(1.0 * mpv), Real(1.0 * qv))); - doAdd = false; - } - do { - termsC.pop_back(); - } while (termsC.back() == 0); - if (termsC.size() <= 1) { - break; - } - } else { - break; - } - } - } + if (coefficents.size() == 2) { // Linear equation ax + b = 0 + if (auto aReal = RecursiveCast(*coefficents[1]); aReal != nullptr) { + if (auto bReal = RecursiveCast(*coefficents[0]); bReal != nullptr) { + // Use Oasis expressions: -b/a + results.push_back(Divide(Multiply(Real(-1), *coefficents[0]), *coefficents[1]).Simplify()); } - if (termsC.size() <= 1) { - break; - } - } - if (termsC.size() <= 1) { - break; } } - coefficents.clear(); - std::reverse(termsC.begin(), termsC.end()); - for (auto i : termsC) { - coefficents.push_back(Real(i * 1.0).Copy()); - } - } - if (coefficents.size() == 2) { - results.push_back(Divide(Multiply(Real(-1), *coefficents[0]), *coefficents[1]).Simplify()); - } - if (coefficents.size() == 3) { - auto& a = coefficents[2]; - auto& b = coefficents[1]; - auto& c = coefficents[0]; - - // Calculate discriminant - auto bSquared = Multiply(*b, *b).Simplify(); - auto fourAC = Multiply(Real(4), Multiply(*a, *c)).Simplify(); - auto discriminant = Subtract(*bSquared, *fourAC).Simplify(); - - if (auto realDisc = RecursiveCast(*discriminant); - realDisc != nullptr && realDisc->GetValue() >= 0) { - - auto negB = Multiply(Real(-1), *b).Simplify(); - auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); - auto twoA = Multiply(Real(2), *a).Simplify(); - - // First, create the numerators for both roots - auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); - auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); - - // Now create the Divide expressions properly - results.push_back(Divide(*numerator1, *twoA).Copy()); - results.push_back(Divide(*numerator2, *twoA).Copy()); + else if (coefficents.size() == 3) { // Quadratic equation ax + b + c = 0 + auto& a = coefficents[2]; + auto& b = coefficents[1]; + auto& c = coefficents[0]; + + // Calculate discriminant + auto bSquared = Multiply(*b, *b).Simplify(); + auto fourAC = Multiply(Real(4), Multiply(*a, *c)).Simplify(); + auto discriminant = Subtract(*bSquared, *fourAC).Simplify(); + + if (auto realDisc = RecursiveCast(*discriminant); + realDisc != nullptr && realDisc->GetValue() >= 0) { + + auto negB = Multiply(Real(-1), *b).Simplify(); + auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + auto twoA = Multiply(Real(2), *a).Simplify(); + + // First, create the numerators for both roots + auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); + auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); + + // Now create the Divide expressions properly + results.push_back(Divide(*numerator1, *twoA).Copy()); + results.push_back(Divide(*numerator2, *twoA).Copy()); + } } } return results; diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 8388a73c..a1b75b62 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -112,66 +112,66 @@ TEST_CASE("imaginary linear polynomial") // REQUIRE(goalSet.size() == 0); // } -TEST_CASE("linear polynomial", "[factor]") +TEST_CASE("linear polynomial test 1: x + 30", "[factor]") { Oasis::Add add { - Oasis::Real(30), - Oasis::Variable("x") + Oasis::Variable("x"), + Oasis::Real(30) }; - auto zeros = add.FindZeros(); - REQUIRE(zeros.size() == 1); - if (zeros.size() == 1) { - auto root = Oasis::RecursiveCast>(*zeros[0]); - REQUIRE(root != nullptr); - REQUIRE(root->GetMostSigOp().GetValue() == -30); - REQUIRE(root->GetLeastSigOp().GetValue() == 1); - } + auto result = add.FindZeros(); + REQUIRE(result.size() == 1); + REQUIRE(result[0]->Is()); + + auto poly_var_res = dynamic_cast(*result[0]); + REQUIRE(poly_var_res.GetValue() == -30); } -TEST_CASE("linear polynomial(Real)", "[factor]") +TEST_CASE("linear polynomial test 2: 3x - 6", "[factor]") { - Oasis::Add add { - Oasis::Real(30), - Oasis::Variable("x") + Oasis::Subtract minus { + Oasis::Multiply{ + Oasis::Real(3), + Oasis::Variable("x") + }, + Oasis::Real(6) }; - auto result = add.Polynomial_Real(); + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + auto result = minus.FindZeros(); REQUIRE(result.size() == 1); REQUIRE(result[0]->Is()); auto poly_var_res = dynamic_cast(*result[0]); - REQUIRE(poly_var_res.GetValue() == -30); + REQUIRE(poly_var_res.GetValue() == 2); } -TEST_CASE("linear polynomial(2x +30)", "[factor]") +TEST_CASE("linear polynomial test 3: 2x + 30", "[factor]") { Oasis::Add add { Oasis::Real(30), - Oasis::Multiply{ // 5x + Oasis::Multiply{ Oasis::Real(2), Oasis::Variable("x") } }; auto zeros = add.FindZeros(); REQUIRE(zeros.size() == 1); - if (zeros.size() == 1) { + REQUIRE(zeros[0]->Is()); - auto root = Oasis::RecursiveCast>(*zeros[0]); - REQUIRE(root != nullptr); - REQUIRE(root->GetMostSigOp().GetValue() == -15); - REQUIRE(root->GetLeastSigOp().GetValue() == 1); - } + auto simplifiedReal = dynamic_cast(*zeros[0]); + REQUIRE(simplifiedReal.GetValue() == -15); } -TEST_CASE("quadratic polynomial", "[factor]") +TEST_CASE("Quadratic polynomial test 1: x² + 5x + 6", "[factor]") { // x² + 5x + 6 Oasis::Add<> add { - Oasis::Real(6), // constant - Oasis::Multiply{ // 5x + Oasis::Real(6), + Oasis::Multiply{ Oasis::Real(5), Oasis::Variable("x") }, - Oasis::Exponent{ // x² + Oasis::Exponent{ Oasis::Variable("x"), Oasis::Real(2) } @@ -201,52 +201,138 @@ TEST_CASE("quadratic polynomial", "[factor]") // Check first root (-2) auto root1 = Oasis::RecursiveCast>(*zeros[0]); REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == -2); - REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + REQUIRE(root1->GetMostSigOp().GetValue() == -4); + REQUIRE(root1->GetLeastSigOp().GetValue() == 2); // Check second root (-3) auto root2 = Oasis::RecursiveCast>(*zeros[1]); REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -6); + REQUIRE(root2->GetLeastSigOp().GetValue() == 2); + } +} + +TEST_CASE("Quadratic polynomial test 2: x² - 2x -3", "[factor]") +{ + // x² - 2x -3 + Oasis::Add<> add{ + Oasis::Exponent{ + Oasis::Variable("x"), + Oasis::Real(2) + }, + Oasis::Multiply { + Oasis::Real(-2), + Oasis::Variable("x") + }, + Oasis::Real(-3) + }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + + auto zeros = add.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + if (zeros.size() == 2) + { + std::cout << "the result has the correct size.\n"; + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -2" << std::endl; + std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + if (root2 != nullptr) + { + std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; + std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; + } + } + + REQUIRE(zeros.size() == 2); + + if (zeros.size() == 2) { + // Check first root (3) + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == 6); + REQUIRE(root1->GetLeastSigOp().GetValue() == 2); + + // Check second root (-1) + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -2); + REQUIRE(root2->GetLeastSigOp().GetValue() == 2); + } +} + +TEST_CASE("Quadratic polynomial test 3: x² - 9", "[factor]") +{ + Oasis::Subtract minus { + Oasis::Exponent { // x² + Oasis::Variable("x"), + Oasis::Real(2) }, + Oasis::Real(9), + }; + auto zeros = minus.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + REQUIRE(zeros.size() == 2); + if (zeros.size() == 2) { + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == 3); + REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); REQUIRE(root2->GetMostSigOp().GetValue() == -3); REQUIRE(root2->GetLeastSigOp().GetValue() == 1); } } -TEST_CASE("simple quadra polynomial", "[factor]") +TEST_CASE("Quadratic polynomial test 4: x² - 16", "[factor]") +{ + Oasis::Subtract minus { + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(2) }, + Oasis::Real(16), + }; + auto zeros = minus.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + REQUIRE(zeros.size() == 2); + if (zeros.size() == 2) { + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == 4); + REQUIRE(root1->GetLeastSigOp().GetValue() == 1); + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -4); + REQUIRE(root2->GetLeastSigOp().GetValue() == 1); + } +} + +TEST_CASE("Quadratic polynomial test 5: x² - 25", "[factor]") { - // x^2 - 2500 Oasis::Subtract minus { Oasis::Exponent { // x² Oasis::Variable("x"), Oasis::Real(2) }, - Oasis::Real(2500), + Oasis::Real(25), }; auto zeros = minus.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; -// if (zeros.size() == 2) -// { -// std::cout << "the result has the correct size.\n"; -// auto root1 = Oasis::RecursiveCast>(*zeros[0]); -// std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be 3" << std::endl; -// std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; -// auto root2 = Oasis::RecursiveCast>(*zeros[1]); -// if (root2 != nullptr) -// { -// std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; -// std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; -// } -// } Oasis::InFixSerializer serializer; OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == 50); + REQUIRE(root1->GetMostSigOp().GetValue() == 5); REQUIRE(root1->GetLeastSigOp().GetValue() == 1); auto root2 = Oasis::RecursiveCast>(*zeros[1]); REQUIRE(root2 != nullptr); - REQUIRE(root2->GetMostSigOp().GetValue() == -50); + REQUIRE(root2->GetMostSigOp().GetValue() == -5); REQUIRE(root2->GetLeastSigOp().GetValue() == 1); } } \ No newline at end of file From e0c17c430949dc4dd9a53825ef305b92ce6ea0b1 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 21 Feb 2025 18:25:57 -0500 Subject: [PATCH 06/21] rational quadratic done, start on cubics --- src/Expression.cpp | 44 ++++++++++ tests/PolynomialTests.cpp | 172 +++++++++++++++++++++++++++++++++++++- 2 files changed, 215 insertions(+), 1 deletion(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 10430d28..741f36f3 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -9,6 +9,7 @@ #include #include #include +#include std::vector getAllFactors(long long n) { @@ -271,6 +272,49 @@ auto Expression::FindZeros() const -> std::vector> // Now create the Divide expressions properly results.push_back(Divide(*numerator1, *twoA).Copy()); results.push_back(Divide(*numerator2, *twoA).Copy()); +// results.push_back(Divide(Add(*negB, *sqrtDisc), *twoA).Copy()); +// results.push_back(Divide(Subtract(*negB, *sqrtDisc), *twoA).Copy()); + } + } + else if (coefficents.size() == 4) { // Cubic equation: ax³ + bx² + cx + d = 0 + long long a = termsC[0]; // coefficient of x³ + long long b = termsC[1]; // coefficient of x² + long long c = termsC[2]; // coefficient of x + long long d = termsC[3]; // constant term + + // Get factors of constant term for possible p values + std::vector p_factors; + for (long long i = 1; i <= std::abs(d); i++) { + if (d % i == 0) { + p_factors.push_back(i); + p_factors.push_back(-i); + } + } + + // Get factors of leading coefficient for possible q values + std::vector q_factors; + for (long long i = 1; i <= std::abs(a); i++) { + if (a % i == 0) { + q_factors.push_back(i); + } + } + + for (long long p : p_factors) { + for (long long q : q_factors) { + // For each potential root p/q, evaluate the polynomial + long long x_p3 = p * p * p; + long long x_p2 = p * p; + long long x_q3 = q * q * q; + + // Evaluate ax³ + bx² + cx + d = 0 at x = p/q + // Multiply all terms by q³ to eliminate denominators: + // a(p³) + b(p²q) + c(pq²) + d(q³) = 0 + long long val = a * x_p3 + b * x_p2 * q + c * p * q * q + d * x_q3; + + if (val == 0) { // If this is a root + results.push_back(Divide(Real(p), Real(q)).Copy()); + } + } } } } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index a1b75b62..f193df1c 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -335,4 +335,174 @@ TEST_CASE("Quadratic polynomial test 5: x² - 25", "[factor]") REQUIRE(root2->GetMostSigOp().GetValue() == -5); REQUIRE(root2->GetLeastSigOp().GetValue() == 1); } -} \ No newline at end of file +} + +TEST_CASE("Rational Quadratic polynomial test 1: 2x² + x - 1", "[factor]") +{ + // 2x² + x - 1 + Oasis::Add<> add{ + Oasis::Multiply{ + Oasis::Real(2), + Oasis::Exponent{ + Oasis::Variable("x"), + Oasis::Real(2) + } + }, + Oasis::Variable("x"), + Oasis::Real(-1) + }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + + auto zeros = add.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + + REQUIRE(zeros.size() == 2); + + if (zeros.size() == 2) { + // Check first root (1/2) + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + auto denominator1 = root1->GetLeastSigOp().GetValue(); + auto numerator1 = root1->GetMostSigOp().GetValue(); + REQUIRE(root1 != nullptr); + REQUIRE(numerator1/denominator1 == 1.0/2.0); + + + // Check second root (-1) + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + auto denominator2 = root2->GetLeastSigOp().GetValue(); + auto numerator2 = root2->GetMostSigOp().GetValue(); + REQUIRE(root2 != nullptr); + REQUIRE(numerator2/denominator2 == -1); + + } +} + +TEST_CASE("Rational Quadratic polynomial test 2: 6x² - 5x + 1", "[factor]") +{ + // 6x² - 5x + 1 + Oasis::Add<> add{ + Oasis::Multiply{ + Oasis::Real(6), + Oasis::Exponent{ + Oasis::Variable("x"), + Oasis::Real(2) + } + }, + Oasis::Multiply{ + Oasis::Real(-5), + Oasis::Variable("x") + }, + Oasis::Real(1) + }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + + auto zeros = add.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + + REQUIRE(zeros.size() == 2); + + if (zeros.size() == 2) { + + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + auto denominator = root1->GetLeastSigOp().GetValue(); + auto numerator = root1->GetMostSigOp().GetValue(); + REQUIRE(root1 != nullptr); + REQUIRE(numerator/denominator == 1.0/2.0); + + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + auto denominator2 = root2->GetLeastSigOp().GetValue(); + auto numerator2 = root2->GetMostSigOp().GetValue(); + REQUIRE(root2 != nullptr); + REQUIRE(numerator2/denominator2 == 1.0/3.0); + } +} + +TEST_CASE("Irrational Quadratic polynomial test 1: x² + 2x - 1/4", "[factor]") +{ + // x² + 2x - 1/4 + Oasis::Add<> add{ + Oasis::Exponent{ + Oasis::Variable("x"), + Oasis::Real(2) + }, + Oasis::Multiply { + Oasis::Real(2), + Oasis::Variable("x") + }, + Oasis::Divide{ + Oasis::Real(1), + Oasis::Real(4) + + } + }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + + auto zeros = add.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + + REQUIRE(zeros.size() == 2); + + if (zeros.size() == 2) { + // Check first root x = -1 - sqrt(5)/2 + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == 2); + REQUIRE(root1->GetLeastSigOp().GetValue() == 4); + + // Check second root x = sqrt(5)/2 - 1 + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -4); + REQUIRE(root2->GetLeastSigOp().GetValue() == 4); + } +} + +TEST_CASE("Cubic polynomial test 1: 3x³ - 16x² + 23x - 6:", "[factor]") +{ + // 3x³ - 16x² + 23x - 6: + Oasis::Add<> cubic{ + Oasis::Multiply{ // 3x³ term + Oasis::Real(3), + Oasis::Exponent{ + Oasis::Variable("x"), + Oasis::Real(3) + } + }, + Oasis::Multiply{ // -16x² term + Oasis::Real(-16), + Oasis::Exponent{ + Oasis::Variable("x"), + Oasis::Real(2) + } + }, + Oasis::Multiply{ // 23x term + Oasis::Real(23), + Oasis::Variable("x") + }, + Oasis::Real(-6) // -6 term + }; + Oasis::InFixSerializer serializer; + OASIS_CAPTURE_WITH_SERIALIZER(serializer, cubic); + + auto zeros = cubic.FindZeros(); + std::cout << "result size: " << zeros.size() << std::endl; + + REQUIRE(zeros.size() == 3); + + if (zeros.size() == 2) { + // Check first root x = -1 - sqrt(5)/2 + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + REQUIRE(root1 != nullptr); + REQUIRE(root1->GetMostSigOp().GetValue() == 2); + REQUIRE(root1->GetLeastSigOp().GetValue() == 4); + + // Check second root x = sqrt(5)/2 - 1 + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + REQUIRE(root2 != nullptr); + REQUIRE(root2->GetMostSigOp().GetValue() == -4); + REQUIRE(root2->GetLeastSigOp().GetValue() == 4); + } +} From c2bcbe23030f06ee41434500df8d43e7158d536d Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 16:43:03 -0500 Subject: [PATCH 07/21] Cubic with fraction or perfect number solutions --- src/Expression.cpp | 29 ++++++++++++++++++++++++++--- tests/PolynomialTests.cpp | 20 +++++++++++++------- 2 files changed, 39 insertions(+), 10 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 741f36f3..b601f3f0 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -10,6 +10,7 @@ #include #include #include +#include std::vector getAllFactors(long long n) { @@ -282,6 +283,9 @@ auto Expression::FindZeros() const -> std::vector> long long c = termsC[2]; // coefficient of x long long d = termsC[3]; // constant term + // To track found roots and avoid duplicates + std::set> found_roots; + // Get factors of constant term for possible p values std::vector p_factors; for (long long i = 1; i <= std::abs(d); i++) { @@ -301,18 +305,37 @@ auto Expression::FindZeros() const -> std::vector> for (long long p : p_factors) { for (long long q : q_factors) { + // Skip if q is 0 + if (q == 0) continue; + + // Simplify the fraction p/q + long long g = std::gcd(std::abs(p), q); + long long num = p / g; + long long den = q / g; + + // Ensure denominator is positive + if (den < 0) { + num = -num; + den = -den; + } + + // Check if we've already found this root + if (found_roots.find({num, den}) != found_roots.end()) { + continue; + } + // For each potential root p/q, evaluate the polynomial long long x_p3 = p * p * p; long long x_p2 = p * p; long long x_q3 = q * q * q; // Evaluate ax³ + bx² + cx + d = 0 at x = p/q - // Multiply all terms by q³ to eliminate denominators: - // a(p³) + b(p²q) + c(pq²) + d(q³) = 0 + // Multiply all terms by q³ to eliminate denominators long long val = a * x_p3 + b * x_p2 * q + c * p * q * q + d * x_q3; if (val == 0) { // If this is a root - results.push_back(Divide(Real(p), Real(q)).Copy()); + results.push_back(Divide(Real(num), Real(den)).Copy()); + found_roots.insert({num, den}); } } } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index f193df1c..3c7169e8 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -492,17 +492,23 @@ TEST_CASE("Cubic polynomial test 1: 3x³ - 16x² + 23x - 6:", "[factor]") REQUIRE(zeros.size() == 3); - if (zeros.size() == 2) { - // Check first root x = -1 - sqrt(5)/2 + if (zeros.size() == 3) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); + auto denominator = root1->GetLeastSigOp().GetValue(); + auto numerator = root1->GetMostSigOp().GetValue(); REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == 2); - REQUIRE(root1->GetLeastSigOp().GetValue() == 4); + REQUIRE(numerator/denominator == 1.0/2.0); - // Check second root x = sqrt(5)/2 - 1 auto root2 = Oasis::RecursiveCast>(*zeros[1]); + auto denominator2 = root2->GetLeastSigOp().GetValue(); + auto numerator2 = root2->GetMostSigOp().GetValue(); REQUIRE(root2 != nullptr); - REQUIRE(root2->GetMostSigOp().GetValue() == -4); - REQUIRE(root2->GetLeastSigOp().GetValue() == 4); + REQUIRE(numerator2/denominator2 == 1.0/3.0); + + auto root3 = Oasis::RecursiveCast>(*zeros[2]); + auto denominator3 = root3->GetLeastSigOp().GetValue(); + auto numerator3 = root3->GetMostSigOp().GetValue(); + REQUIRE(root2 != nullptr); + REQUIRE(numerator3/denominator3 == 3.0/1.0); } } From 5966670a0d890c78571ef284fde96d4ac53ae55b Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 16:57:57 -0500 Subject: [PATCH 08/21] Cubic with fraction or perfect number solutions --- include/Oasis/Expression.hpp | 10 +-- src/Expression.cpp | 131 ----------------------------------- tests/AddTests.cpp | 3 +- tests/PolynomialTests.cpp | 30 +++----- 4 files changed, 12 insertions(+), 162 deletions(-) diff --git a/include/Oasis/Expression.hpp b/include/Oasis/Expression.hpp index d89e4df6..6bbdb30a 100644 --- a/include/Oasis/Expression.hpp +++ b/include/Oasis/Expression.hpp @@ -92,15 +92,7 @@ class Expression { * @tparam origonalExpresion The expression for which all the factors will be found. */ auto FindZeros() const -> std::vector>; - - /** - * Finds the real number roots of a polynomial expression. - * This function handles special cases like x² - n (where n is a perfect square) - * and general quadratic equations (ax² + bx + c). - * - * @return A vector of Real values representing the real roots of the polynomial - */ - [[nodiscard]] auto Polynomial_Real() const -> std::vector>; + /** * Gets the category of this expression. diff --git a/src/Expression.cpp b/src/Expression.cpp index b601f3f0..4da45fe7 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -188,58 +188,6 @@ auto Expression::FindZeros() const -> std::vector> termsC.push_back(lround(value)); } } -// if (termsC.size() == coefficents.size()) { -// std::reverse(termsC.begin(), termsC.end()); -// for (auto pv : getAllFactors(termsC.back())) { -// for (auto qv : getAllFactors(termsC.front())) { -// if (gcf(pv, qv) == 1) { -// for (long long sign : { -1, 1 }) { -// bool doAdd = true; -// while (true) { -// long long mpv = pv * sign; -// std::vector newTermsC; -// long long h = 0; -// for (long long i : termsC) { -// h *= mpv; -// if (h % qv != 0) { -// break; -// } -// h /= qv; -// h += i; -// newTermsC.push_back(h); -// } -// if (newTermsC.size() == termsC.size() && newTermsC.back() == 0) { -// termsC = newTermsC; -// if (doAdd) { -// results.push_back(std::make_unique>(Real(1.0 * mpv), Real(1.0 * qv))); -// doAdd = false; -// } -// do { -// termsC.pop_back(); -// } while (termsC.back() == 0); -// if (termsC.size() <= 1) { -// break; -// } -// } else { -// break; -// } -// } -// } -// } -// if (termsC.size() <= 1) { -// break; -// } -// } -// if (termsC.size() <= 1) { -// break; -// } -// } -// coefficents.clear(); -// std::reverse(termsC.begin(), termsC.end()); -// for (auto i : termsC) { -// coefficents.push_back(Real(i * 1.0).Copy()); -// } -// } if (termsC.size() == coefficents.size()) { if (coefficents.size() == 2) { // Linear equation ax + b = 0 if (auto aReal = RecursiveCast(*coefficents[1]); aReal != nullptr) { @@ -344,85 +292,6 @@ auto Expression::FindZeros() const -> std::vector> return results; } -auto Expression::Polynomial_Real() const -> std::vector> { - std::vector> results; - std::vector> termsE; - - // First handle x² - n special case - if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { - if (auto leftTerm = RecursiveCast>(subCase->GetMostSigOp()); - leftTerm != nullptr) { - if (leftTerm->GetLeastSigOp().GetValue() == 2) { - if (auto rightTerm = RecursiveCast(subCase->GetLeastSigOp()); - rightTerm != nullptr) { - double n = rightTerm->GetValue(); - double sqrtN = std::sqrt(n); - if (n > 0 && sqrtN == std::floor(sqrtN)) { - results.push_back(Divide(Real(sqrtN), Real(1)).Copy()); - results.push_back(Divide(Real(-sqrtN), Real(1)).Copy()); - return results; - } - } - } - } - } - - // Process terms - if (auto addCase = RecursiveCast>(*this); addCase != nullptr) { - addCase->Flatten(termsE); - } else { - termsE.push_back(Copy()); - } - - // Collect coefficients - std::vector> coefficients; - for (const auto& term : termsE) { - if (auto realTerm = RecursiveCast(*term); realTerm != nullptr) { - coefficients.push_back(realTerm->Copy()); - } - } - - // Handle linear equation (ax + b) - if (coefficients.size() == 2) { - if (auto aReal = RecursiveCast(*coefficients[1]); aReal != nullptr) { - if (auto bReal = RecursiveCast(*coefficients[0]); bReal != nullptr) { - double a = aReal->GetValue(); - double b = bReal->GetValue(); - - // ax + b = 0 → x = -b/a - if (a != 0) { - double root = -b / a; - results.push_back(Divide(Real(root), Real(1)).Copy()); - } - } - } - } - - // Handle quadratic equation (for ax² + bx + c) - if (coefficients.size() == 3) { - if (auto aReal = RecursiveCast(*coefficients[2]); aReal != nullptr) { - if (auto bReal = RecursiveCast(*coefficients[1]); bReal != nullptr) { - if (auto cReal = RecursiveCast(*coefficients[0]); cReal != nullptr) { - double a = aReal->GetValue(); - double b = bReal->GetValue(); - double c = cReal->GetValue(); - - double discriminant = b * b - 4 * a * c; - if (discriminant >= 0) { - double sqrtDisc = std::sqrt(discriminant); - double root1 = (-b + sqrtDisc) / (2 * a); - double root2 = (-b - sqrtDisc) / (2 * a); - - results.push_back(Divide(Real(root1), Real(1)).Copy()); - results.push_back(Divide(Real(root2), Real(1)).Copy()); - } - } - } - } - } - - return results; -} auto Expression::GetCategory() const -> uint32_t { diff --git a/tests/AddTests.cpp b/tests/AddTests.cpp index 88b2cdf0..236e2f6a 100644 --- a/tests/AddTests.cpp +++ b/tests/AddTests.cpp @@ -23,8 +23,7 @@ TEST_CASE("Addition", "[Add]") Oasis::Real { 3.0 } }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + OASIS_CAPTURE_WITH_SERIALIZER(add); auto simplified = add.Simplify(); REQUIRE(simplified->Is()); diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 3c7169e8..7ee0c140 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -135,8 +135,7 @@ TEST_CASE("linear polynomial test 2: 3x - 6", "[factor]") }, Oasis::Real(6) }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + OASIS_CAPTURE_WITH_SERIALIZER(minus); auto result = minus.FindZeros(); REQUIRE(result.size() == 1); REQUIRE(result[0]->Is()); @@ -176,8 +175,7 @@ TEST_CASE("Quadratic polynomial test 1: x² + 5x + 6", "[factor]") Oasis::Real(2) } }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; @@ -226,8 +224,7 @@ TEST_CASE("Quadratic polynomial test 2: x² - 2x -3", "[factor]") }, Oasis::Real(-3) }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; @@ -272,8 +269,7 @@ TEST_CASE("Quadratic polynomial test 3: x² - 9", "[factor]") }; auto zeros = minus.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); @@ -297,8 +293,7 @@ TEST_CASE("Quadratic polynomial test 4: x² - 16", "[factor]") }; auto zeros = minus.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); @@ -322,8 +317,7 @@ TEST_CASE("Quadratic polynomial test 5: x² - 25", "[factor]") }; auto zeros = minus.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, minus); + OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); @@ -351,8 +345,7 @@ TEST_CASE("Rational Quadratic polynomial test 1: 2x² + x - 1", "[factor]") Oasis::Variable("x"), Oasis::Real(-1) }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; @@ -395,8 +388,7 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x² - 5x + 1", "[factor]") }, Oasis::Real(1) }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; @@ -437,8 +429,7 @@ TEST_CASE("Irrational Quadratic polynomial test 1: x² + 2x - 1/4", "[factor]") } }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, add); + OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; @@ -484,8 +475,7 @@ TEST_CASE("Cubic polynomial test 1: 3x³ - 16x² + 23x - 6:", "[factor]") }, Oasis::Real(-6) // -6 term }; - Oasis::InFixSerializer serializer; - OASIS_CAPTURE_WITH_SERIALIZER(serializer, cubic); + OASIS_CAPTURE_WITH_SERIALIZER(cubic); auto zeros = cubic.FindZeros(); std::cout << "result size: " << zeros.size() << std::endl; From 8588ce6fda4b41a3558a9d3e7e63e20e461b1ec3 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 17:01:54 -0500 Subject: [PATCH 09/21] removed couts --- src/Expression.cpp | 4 ---- tests/PolynomialTests.cpp | 35 ----------------------------------- 2 files changed, 39 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 4da45fe7..c1993f1d 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -53,7 +53,6 @@ namespace Oasis { */ auto Expression::FindZeros() const -> std::vector> { - std::cout << "start of the function" << std::endl; std::vector> results; std::vector> termsE; if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { @@ -162,8 +161,6 @@ auto Expression::FindZeros() const -> std::vector> while (posCoefficents.size() > 0 && RecursiveCast(*posCoefficents.back()) != nullptr && RecursiveCast(*posCoefficents.back())->GetValue() == 0) { posCoefficents.pop_back(); } - std::cout << "posCoefficents size: " << posCoefficents.size() << std::endl; - std::cout << "negCoefficents size: " << negCoefficents.size() << std::endl; std::vector> coefficents; for (size_t i = negCoefficents.size(); i > 1; i--) { coefficents.push_back(negCoefficents[i - 1]->Simplify()); @@ -171,7 +168,6 @@ auto Expression::FindZeros() const -> std::vector> for (const std::unique_ptr& i : posCoefficents) { coefficents.push_back(i->Simplify()); } - std::cout << "Final coefficients size: " << coefficents.size() << std::endl; if (coefficents.size() <= 1) { return {}; } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 7ee0c140..06c66427 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -178,20 +178,6 @@ TEST_CASE("Quadratic polynomial test 1: x² + 5x + 6", "[factor]") OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; - if (zeros.size() == 2) - { - std::cout << "the result has the correct size.\n"; - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -2" << std::endl; - std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - if (root2 != nullptr) - { - std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; - std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; - } - } REQUIRE(zeros.size() == 2); @@ -227,20 +213,6 @@ TEST_CASE("Quadratic polynomial test 2: x² - 2x -3", "[factor]") OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; - if (zeros.size() == 2) - { - std::cout << "the result has the correct size.\n"; - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - std::cout << "root1 mostsig value: " << root1->GetMostSigOp().GetValue() << " where should be -2" << std::endl; - std::cout << "root1 leastsig value: " << root1->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - if (root2 != nullptr) - { - std::cout << "root2 mostsig value: " << root2->GetMostSigOp().GetValue() << " where should be -3" << std::endl; - std::cout << "root2 leastsig value: " << root2->GetLeastSigOp().GetValue() << " where should be 1" << std::endl; - } - } REQUIRE(zeros.size() == 2); @@ -268,7 +240,6 @@ TEST_CASE("Quadratic polynomial test 3: x² - 9", "[factor]") Oasis::Real(9), }; auto zeros = minus.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { @@ -292,7 +263,6 @@ TEST_CASE("Quadratic polynomial test 4: x² - 16", "[factor]") Oasis::Real(16), }; auto zeros = minus.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { @@ -316,7 +286,6 @@ TEST_CASE("Quadratic polynomial test 5: x² - 25", "[factor]") Oasis::Real(25), }; auto zeros = minus.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { @@ -348,7 +317,6 @@ TEST_CASE("Rational Quadratic polynomial test 1: 2x² + x - 1", "[factor]") OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; REQUIRE(zeros.size() == 2); @@ -391,7 +359,6 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x² - 5x + 1", "[factor]") OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; REQUIRE(zeros.size() == 2); @@ -432,7 +399,6 @@ TEST_CASE("Irrational Quadratic polynomial test 1: x² + 2x - 1/4", "[factor]") OASIS_CAPTURE_WITH_SERIALIZER(add); auto zeros = add.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; REQUIRE(zeros.size() == 2); @@ -478,7 +444,6 @@ TEST_CASE("Cubic polynomial test 1: 3x³ - 16x² + 23x - 6:", "[factor]") OASIS_CAPTURE_WITH_SERIALIZER(cubic); auto zeros = cubic.FindZeros(); - std::cout << "result size: " << zeros.size() << std::endl; REQUIRE(zeros.size() == 3); From cea7dea41fa34f970efab29441860742fd233cfa Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 17:07:24 -0500 Subject: [PATCH 10/21] clang format --- src/Expression.cpp | 42 ++++++++++++++++++++---------------------- 1 file changed, 20 insertions(+), 22 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index c1993f1d..36d11138 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -9,8 +9,8 @@ #include #include #include -#include #include +#include std::vector getAllFactors(long long n) { @@ -70,14 +70,14 @@ auto Expression::FindZeros() const -> std::vector> // Instead of creating Real values directly, create proper Binary Expressions // For positive root: n^(1/2) auto posRoot = std::make_unique>( - Real(sqrtN), // MostSigOp - Real(1) // LeastSigOp + Real(sqrtN), // MostSigOp + Real(1) // LeastSigOp ); // For negative root: -n^(1/2) auto negRoot = std::make_unique>( - Real(-sqrtN), // MostSigOp - Real(1) // LeastSigOp + Real(-sqrtN), // MostSigOp + Real(1) // LeastSigOp ); results.push_back(std::move(posRoot)); @@ -92,7 +92,7 @@ auto Expression::FindZeros() const -> std::vector> addCase->Flatten(termsE); } else if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { // Handle subtraction by converting to addition - termsE.push_back(subCase->GetMostSigOp().Copy()); // First term + termsE.push_back(subCase->GetMostSigOp().Copy()); // First term // Add negative of second term auto negTerm = Multiply(Real(-1), subCase->GetLeastSigOp()).Copy(); termsE.push_back(std::move(negTerm)); @@ -185,15 +185,14 @@ auto Expression::FindZeros() const -> std::vector> } } if (termsC.size() == coefficents.size()) { - if (coefficents.size() == 2) { // Linear equation ax + b = 0 + if (coefficents.size() == 2) { // Linear equation ax + b = 0 if (auto aReal = RecursiveCast(*coefficents[1]); aReal != nullptr) { if (auto bReal = RecursiveCast(*coefficents[0]); bReal != nullptr) { // Use Oasis expressions: -b/a results.push_back(Divide(Multiply(Real(-1), *coefficents[0]), *coefficents[1]).Simplify()); } } - } - else if (coefficents.size() == 3) { // Quadratic equation ax + b + c = 0 + } else if (coefficents.size() == 3) { // Quadratic equation ax + b + c = 0 auto& a = coefficents[2]; auto& b = coefficents[1]; auto& c = coefficents[0]; @@ -217,15 +216,14 @@ auto Expression::FindZeros() const -> std::vector> // Now create the Divide expressions properly results.push_back(Divide(*numerator1, *twoA).Copy()); results.push_back(Divide(*numerator2, *twoA).Copy()); -// results.push_back(Divide(Add(*negB, *sqrtDisc), *twoA).Copy()); -// results.push_back(Divide(Subtract(*negB, *sqrtDisc), *twoA).Copy()); + // results.push_back(Divide(Add(*negB, *sqrtDisc), *twoA).Copy()); + // results.push_back(Divide(Subtract(*negB, *sqrtDisc), *twoA).Copy()); } - } - else if (coefficents.size() == 4) { // Cubic equation: ax³ + bx² + cx + d = 0 - long long a = termsC[0]; // coefficient of x³ - long long b = termsC[1]; // coefficient of x² - long long c = termsC[2]; // coefficient of x - long long d = termsC[3]; // constant term + } else if (coefficents.size() == 4) { // Cubic equation: ax³ + bx² + cx + d = 0 + long long a = termsC[0]; // coefficient of x³ + long long b = termsC[1]; // coefficient of x² + long long c = termsC[2]; // coefficient of x + long long d = termsC[3]; // constant term // To track found roots and avoid duplicates std::set> found_roots; @@ -250,7 +248,8 @@ auto Expression::FindZeros() const -> std::vector> for (long long p : p_factors) { for (long long q : q_factors) { // Skip if q is 0 - if (q == 0) continue; + if (q == 0) + continue; // Simplify the fraction p/q long long g = std::gcd(std::abs(p), q); @@ -264,7 +263,7 @@ auto Expression::FindZeros() const -> std::vector> } // Check if we've already found this root - if (found_roots.find({num, den}) != found_roots.end()) { + if (found_roots.find({ num, den }) != found_roots.end()) { continue; } @@ -277,9 +276,9 @@ auto Expression::FindZeros() const -> std::vector> // Multiply all terms by q³ to eliminate denominators long long val = a * x_p3 + b * x_p2 * q + c * p * q * q + d * x_q3; - if (val == 0) { // If this is a root + if (val == 0) { // If this is a root results.push_back(Divide(Real(num), Real(den)).Copy()); - found_roots.insert({num, den}); + found_roots.insert({ num, den }); } } } @@ -288,7 +287,6 @@ auto Expression::FindZeros() const -> std::vector> return results; } - auto Expression::GetCategory() const -> uint32_t { return 0; From dea85cf8956272d4136a9016649895435a36bbd2 Mon Sep 17 00:00:00 2001 From: github-actions Date: Tue, 25 Feb 2025 22:08:58 +0000 Subject: [PATCH 11/21] [Actions] Format CMakeLists.txt --- tests/CMakeLists.txt | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 98901002..bfde36f4 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -9,6 +9,7 @@ set(Oasis_TESTS # cmake-format: sortable AddTests.cpp BinaryExpressionTests.cpp + Common.hpp DifferentiateTests.cpp DivideTests.cpp ExponentTests.cpp @@ -21,8 +22,7 @@ set(Oasis_TESTS NegateTests.cpp PolynomialTests.cpp SubtractTests.cpp - UnaryExpressionTests.cpp - Common.hpp) + UnaryExpressionTests.cpp) # Adds an executable target called "OasisTests" to be built from sources files. add_executable(OasisTests ${Oasis_TESTS}) @@ -36,12 +36,12 @@ endif() target_link_libraries(OasisTests PRIVATE Oasis::Oasis Catch2::Catch2WithMain) -if (OASIS_BUILD_IO) +if(OASIS_BUILD_IO) target_link_libraries(OasisTests PRIVATE Oasis::IO) target_compile_definitions(OasisTests PRIVATE OASIS_IO_ENABLED) -endif () +endif() - # Automatically registers the tests with CTest. +# Automatically registers the tests with CTest. list(APPEND CMAKE_MODULE_PATH ${catch2_SOURCE_DIR}/extras) include(Catch) catch_discover_tests(OasisTests) From 58817d44f8c9f8be1d80fd18cfbd0cfb32915962 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 17:12:05 -0500 Subject: [PATCH 12/21] Trigger workflow From 462f7952227e5d1016b0990e7e597ec4dea132c8 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 17:15:11 -0500 Subject: [PATCH 13/21] clang format --- include/Oasis/Expression.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/Oasis/Expression.hpp b/include/Oasis/Expression.hpp index 6bbdb30a..c5f5c02a 100644 --- a/include/Oasis/Expression.hpp +++ b/include/Oasis/Expression.hpp @@ -92,7 +92,6 @@ class Expression { * @tparam origonalExpresion The expression for which all the factors will be found. */ auto FindZeros() const -> std::vector>; - /** * Gets the category of this expression. From 339972974dd1c0826dbab99833314dfe5832b8d7 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Tue, 25 Feb 2025 17:32:05 -0500 Subject: [PATCH 14/21] CMAKE_ON_LINUX check CMAKE_ON_WINDOWS check --- src/Expression.cpp | 2 +- tests/PolynomialTests.cpp | 90 --------------------------------------- 2 files changed, 1 insertion(+), 91 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 36d11138..f32ab57f 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -277,7 +277,7 @@ auto Expression::FindZeros() const -> std::vector> long long val = a * x_p3 + b * x_p2 * q + c * p * q * q + d * x_q3; if (val == 0) { // If this is a root - results.push_back(Divide(Real(num), Real(den)).Copy()); + results.push_back(Divide(Real(static_cast(num)), Real(static_cast(den))).Copy()); found_roots.insert({ num, den }); } } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 06c66427..33eb20d4 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -22,96 +22,6 @@ #include #include -// TODO: Figure out what's going out here -// TEST_CASE("7th degree polynomial with rational roots", "[factor][duplicateRoot]") -// { -// std::vector> vec; -// long offset = -3; -// std::vector vecL = { 24750, -200'925, 573'625, -631'406, 79184, 247'799, -92631, 8820 }; -// for (size_t i = 0; i < vecL.size(); i++) { -// Oasis::Real num = Oasis::Real(vecL[i]); -// long exp = ((long)i) + offset; -// if (exp < -1) { -// vec.push_back(Oasis::Divide(num, Oasis::Exponent(Oasis::Variable("x"), Oasis::Real(-exp * 1.0))).Copy()); -// } else if (exp == -1) { -// vec.push_back(Oasis::Divide(num, Oasis::Variable("x")).Copy()); -// } else if (exp == 0) { -// vec.push_back(num.Copy()); -// } else if (exp == 1) { -// vec.push_back(Oasis::Multiply(num, Oasis::Variable("x")).Copy()); -// } else { -// vec.push_back(Oasis::Multiply(num, Oasis::Exponent(Oasis::Variable("x"), Oasis::Real(exp * 1.0))).Copy()); -// } -// } -// auto add = Oasis::BuildFromVector(vec); -// auto zeros = add->FindZeros(); -// REQUIRE(zeros.size() == 6); -// std::set> goalSet = { std::tuple(1, 3), std::tuple(6, 7), std::tuple(3, 7), std::tuple(-5, 3), std::tuple(11, 20), std::tuple(5, 1) }; -// for (auto& i : zeros) { -// auto divideCase = Oasis::RecursiveCast>(*i); -// REQUIRE(divideCase != nullptr); -// std::tuple asTuple = std::tuple(std::lround(divideCase->GetMostSigOp().GetValue()), std::lround(divideCase->GetLeastSigOp().GetValue())); -// REQUIRE(goalSet.contains(asTuple)); -// goalSet.erase(asTuple); -// } -// } - -TEST_CASE("imaginary linear polynomial") -{ - Oasis::Add add { - Oasis::Imaginary(), - Oasis::Variable("x") - }; - auto zeros = add.FindZeros(); - REQUIRE(zeros.size() == 1); - if (zeros.size() == 1) { - auto root = Oasis::RecursiveCast>(*zeros[0]); - REQUIRE(root != nullptr); - REQUIRE(root->GetMostSigOp().GetValue() == -1); - } -} - -// TODO: Figure out what's going out here -// TEST_CASE("irrational quadratic", "[quadraticFormula]") -// { -// std::vector> vec; -// long offset = -3; -// std::vector vecL = { -1, 1, 1 }; -// for (size_t i = 0; i < vecL.size(); i++) { -// Oasis::Real num = Oasis::Real(vecL[i]); -// long exp = ((long)i) + offset; -// if (exp < -1) { -// vec.push_back(Oasis::Divide(num, Oasis::Exponent(Oasis::Variable("x"), Oasis::Real(-exp))).Copy()); -// } else if (exp == -1) { -// vec.push_back(Oasis::Divide(num, Oasis::Variable("x")).Copy()); -// } else if (exp == 0) { -// vec.push_back(num.Copy()); -// } else if (exp == 1) { -// vec.push_back(Oasis::Multiply(num, Oasis::Variable("x")).Copy()); -// } else { -// vec.push_back(Oasis::Multiply(num, Oasis::Exponent(Oasis::Variable("x"), Oasis::Real(exp))).Copy()); -// } -// } -// auto add = Oasis::BuildFromVector(vec); -// auto zeros = add->FindZeros(); -// REQUIRE(zeros.size() == 2); -// auto negOne = Oasis::Real(-1); -// auto two = Oasis::Real(2); -// auto root5 = Oasis::Exponent(Oasis::Real(5), Oasis::Divide(Oasis::Real(1), two)); -// std::list> goalSet = {}; -// goalSet.push_back(Oasis::Divide(Oasis::Add(negOne, root5), two).Copy()); -// goalSet.push_back(Oasis::Divide(Oasis::Subtract(negOne, root5), two).Copy()); -// for (auto& i : zeros) { -// for (auto i2 = goalSet.begin(); i2 != goalSet.end(); i2++) { -// if ((*i2)->Equals(*i)) { -// goalSet.erase(i2); -// break; -// } -// } -// } -// REQUIRE(goalSet.size() == 0); -// } - TEST_CASE("linear polynomial test 1: x + 30", "[factor]") { Oasis::Add add { From 1fe0105ea02fccf85005c6b145f215d453afe8f6 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 14 Mar 2025 16:26:41 -0400 Subject: [PATCH 15/21] CMAKE_ON_LINUX check CMAKE_ON_WINDOWS check --- tests/PolynomialTests.cpp | 38 -------------------------------------- 1 file changed, 38 deletions(-) diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 33eb20d4..7e57d564 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -288,44 +288,6 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x² - 5x + 1", "[factor]") } } -TEST_CASE("Irrational Quadratic polynomial test 1: x² + 2x - 1/4", "[factor]") -{ - // x² + 2x - 1/4 - Oasis::Add<> add{ - Oasis::Exponent{ - Oasis::Variable("x"), - Oasis::Real(2) - }, - Oasis::Multiply { - Oasis::Real(2), - Oasis::Variable("x") - }, - Oasis::Divide{ - Oasis::Real(1), - Oasis::Real(4) - - } - }; - OASIS_CAPTURE_WITH_SERIALIZER(add); - - auto zeros = add.FindZeros(); - - REQUIRE(zeros.size() == 2); - - if (zeros.size() == 2) { - // Check first root x = -1 - sqrt(5)/2 - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == 2); - REQUIRE(root1->GetLeastSigOp().GetValue() == 4); - - // Check second root x = sqrt(5)/2 - 1 - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - REQUIRE(root2 != nullptr); - REQUIRE(root2->GetMostSigOp().GetValue() == -4); - REQUIRE(root2->GetLeastSigOp().GetValue() == 4); - } -} TEST_CASE("Cubic polynomial test 1: 3x³ - 16x² + 23x - 6:", "[factor]") { From d9780be18d654cd0a2c3d8483d89fe296016e790 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 14 Mar 2025 16:52:50 -0400 Subject: [PATCH 16/21] CMAKE_ON_WINDOWS check --- tests/PolynomialTests.cpp | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 7e57d564..646db569 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -71,9 +71,9 @@ TEST_CASE("linear polynomial test 3: 2x + 30", "[factor]") REQUIRE(simplifiedReal.GetValue() == -15); } -TEST_CASE("Quadratic polynomial test 1: x² + 5x + 6", "[factor]") +TEST_CASE("Quadratic polynomial test 1: x^2 + 5x + 6", "[factor]") { - // x² + 5x + 6 + // x^2 + 5x + 6 Oasis::Add<> add { Oasis::Real(6), Oasis::Multiply{ @@ -106,9 +106,9 @@ TEST_CASE("Quadratic polynomial test 1: x² + 5x + 6", "[factor]") } } -TEST_CASE("Quadratic polynomial test 2: x² - 2x -3", "[factor]") +TEST_CASE("Quadratic polynomial test 2: x^2 - 2x -3", "[factor]") { - // x² - 2x -3 + // x^2 - 2x -3 Oasis::Add<> add{ Oasis::Exponent{ Oasis::Variable("x"), @@ -141,10 +141,10 @@ TEST_CASE("Quadratic polynomial test 2: x² - 2x -3", "[factor]") } } -TEST_CASE("Quadratic polynomial test 3: x² - 9", "[factor]") +TEST_CASE("Quadratic polynomial test 3: x^2 - 9", "[factor]") { Oasis::Subtract minus { - Oasis::Exponent { // x² + Oasis::Exponent { // x^2 Oasis::Variable("x"), Oasis::Real(2) }, Oasis::Real(9), @@ -164,7 +164,7 @@ TEST_CASE("Quadratic polynomial test 3: x² - 9", "[factor]") } } -TEST_CASE("Quadratic polynomial test 4: x² - 16", "[factor]") +TEST_CASE("Quadratic polynomial test 4: x^2 - 16", "[factor]") { Oasis::Subtract minus { Oasis::Exponent { @@ -187,10 +187,10 @@ TEST_CASE("Quadratic polynomial test 4: x² - 16", "[factor]") } } -TEST_CASE("Quadratic polynomial test 5: x² - 25", "[factor]") +TEST_CASE("Quadratic polynomial test 5: x^2 - 25", "[factor]") { Oasis::Subtract minus { - Oasis::Exponent { // x² + Oasis::Exponent { // x^2 Oasis::Variable("x"), Oasis::Real(2) }, Oasis::Real(25), @@ -210,9 +210,9 @@ TEST_CASE("Quadratic polynomial test 5: x² - 25", "[factor]") } } -TEST_CASE("Rational Quadratic polynomial test 1: 2x² + x - 1", "[factor]") +TEST_CASE("Rational Quadratic polynomial test 1: 2x^2 + x - 1", "[factor]") { - // 2x² + x - 1 + // 2x^2 + x - 1 Oasis::Add<> add{ Oasis::Multiply{ Oasis::Real(2), @@ -249,9 +249,9 @@ TEST_CASE("Rational Quadratic polynomial test 1: 2x² + x - 1", "[factor]") } } -TEST_CASE("Rational Quadratic polynomial test 2: 6x² - 5x + 1", "[factor]") +TEST_CASE("Rational Quadratic polynomial test 2: 6x^2 - 5x + 1", "[factor]") { - // 6x² - 5x + 1 + // 6x^2 - 5x + 1 Oasis::Add<> add{ Oasis::Multiply{ Oasis::Real(6), @@ -289,18 +289,18 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x² - 5x + 1", "[factor]") } -TEST_CASE("Cubic polynomial test 1: 3x³ - 16x² + 23x - 6:", "[factor]") +TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") { - // 3x³ - 16x² + 23x - 6: + // 3x^3 - 16x^2 + 23x - 6: Oasis::Add<> cubic{ - Oasis::Multiply{ // 3x³ term + Oasis::Multiply{ // 3x^3 term Oasis::Real(3), Oasis::Exponent{ Oasis::Variable("x"), Oasis::Real(3) } }, - Oasis::Multiply{ // -16x² term + Oasis::Multiply{ // -16x^2 term Oasis::Real(-16), Oasis::Exponent{ Oasis::Variable("x"), From 47252d2ca22822730f23fd260d4fd1b73395dc14 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 14 Mar 2025 17:06:45 -0400 Subject: [PATCH 17/21] simple quartic non-irrational --- src/Expression.cpp | 78 +++++++++++++++++++++++- tests/PolynomialTests.cpp | 122 ++++++++++++++++++-------------------- 2 files changed, 134 insertions(+), 66 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index f32ab57f..89a5fc00 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -8,7 +8,6 @@ #include #include #include -#include #include #include @@ -282,6 +281,83 @@ auto Expression::FindZeros() const -> std::vector> } } } + } else if (coefficents.size() == 5) { // Quartic equation ax⁴ + bx³ + cx² + dx + e = 0 + // coefficients + long long a = 0, b = 0, c = 0, d = 0, e = 0; + + // Convert coefficients to numbers if possible + if (auto aReal = RecursiveCast(*coefficents[4]); aReal != nullptr) a = static_cast(aReal->GetValue()); + if (auto bReal = RecursiveCast(*coefficents[3]); bReal != nullptr) b = static_cast(bReal->GetValue()); + if (auto cReal = RecursiveCast(*coefficents[2]); cReal != nullptr) c = static_cast(cReal->GetValue()); + if (auto dReal = RecursiveCast(*coefficents[1]); dReal != nullptr) d = static_cast(dReal->GetValue()); + if (auto eReal = RecursiveCast(*coefficents[0]); eReal != nullptr) e = static_cast(eReal->GetValue()); + + // Find potential rational roots using the rational root theorem + // Possible rational roots are p/q where: + // - p is a factor of the constant term (e) + // - q is a factor of the leading coefficient (a) + + // Get factors of constant term for possible p values + std::vector p_factors; + for (long long i = 1; i <= std::abs(e); i++) { + if (e % i == 0) { + p_factors.push_back(i); + p_factors.push_back(-i); + } + } + + // Get factors of leading coefficient for possible q values + std::vector q_factors; + for (long long i = 1; i <= std::abs(a); i++) { + if (a % i == 0) { + q_factors.push_back(i); + } + } + + // To track found roots and avoid duplicates + std::set> found_roots; + + for (long long p : p_factors) { + for (long long q : q_factors) { + // Skip if q is 0 + if (q == 0) continue; + + // Simplify the fraction p/q + long long g = std::gcd(std::abs(p), q); + long long num = p / g; + long long den = q / g; + + // Ensure denominator is positive + if (den < 0) { + num = -num; + den = -den; + } + + // Check if we've already found this root + if (found_roots.find({num, den}) != found_roots.end()) { + continue; + } + + // Evaluate the polynomial at p/q using synthetic division + // For a quartic: a(p/q)⁴ + b(p/q)³ + c(p/q)² + d(p/q) + e + + // Multiply by q⁴ to clear denominators: + // a*p⁴ + b*p³*q + c*p²*q² + d*p*q³ + e*q⁴ + long long p2 = p * p; + long long p3 = p2 * p; + long long p4 = p3 * p; + long long q2 = q * q; + long long q3 = q2 * q; + long long q4 = q3 * q; + + long long val = a * p4 + b * p3 * q + c * p2 * q2 + d * p * q3 + e * q4; + + if (val == 0) { // If this is a root + results.push_back(Divide(Real(static_cast(num)), Real(static_cast(den))).Copy()); + found_roots.insert({num, den}); + } + } + } } } return results; diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 646db569..5a8ca246 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -3,6 +3,7 @@ // #include "catch2/catch_test_macros.hpp" +#include "Common.hpp" #include "Oasis/Add.hpp" #include "Oasis/Divide.hpp" #include "Oasis/Exponent.hpp" @@ -11,16 +12,11 @@ #include "Oasis/Multiply.hpp" #include "Oasis/Real.hpp" #include "Oasis/RecursiveCast.hpp" -#include "Oasis/InFixSerializer.hpp" #include "Oasis/Subtract.hpp" #include "Oasis/Variable.hpp" -#include "Common.hpp" -#include #include -#include #include -#include TEST_CASE("linear polynomial test 1: x + 30", "[factor]") { @@ -39,10 +35,9 @@ TEST_CASE("linear polynomial test 1: x + 30", "[factor]") TEST_CASE("linear polynomial test 2: 3x - 6", "[factor]") { Oasis::Subtract minus { - Oasis::Multiply{ + Oasis::Multiply { Oasis::Real(3), - Oasis::Variable("x") - }, + Oasis::Variable("x") }, Oasis::Real(6) }; OASIS_CAPTURE_WITH_SERIALIZER(minus); @@ -58,10 +53,9 @@ TEST_CASE("linear polynomial test 3: 2x + 30", "[factor]") { Oasis::Add add { Oasis::Real(30), - Oasis::Multiply{ + Oasis::Multiply { Oasis::Real(2), - Oasis::Variable("x") - } + Oasis::Variable("x") } }; auto zeros = add.FindZeros(); REQUIRE(zeros.size() == 1); @@ -76,14 +70,12 @@ TEST_CASE("Quadratic polynomial test 1: x^2 + 5x + 6", "[factor]") // x^2 + 5x + 6 Oasis::Add<> add { Oasis::Real(6), - Oasis::Multiply{ + Oasis::Multiply { Oasis::Real(5), - Oasis::Variable("x") - }, - Oasis::Exponent{ + Oasis::Variable("x") }, + Oasis::Exponent { Oasis::Variable("x"), - Oasis::Real(2) - } + Oasis::Real(2) } }; OASIS_CAPTURE_WITH_SERIALIZER(add); @@ -109,15 +101,13 @@ TEST_CASE("Quadratic polynomial test 1: x^2 + 5x + 6", "[factor]") TEST_CASE("Quadratic polynomial test 2: x^2 - 2x -3", "[factor]") { // x^2 - 2x -3 - Oasis::Add<> add{ - Oasis::Exponent{ + Oasis::Add<> add { + Oasis::Exponent { Oasis::Variable("x"), - Oasis::Real(2) - }, + Oasis::Real(2) }, Oasis::Multiply { Oasis::Real(-2), - Oasis::Variable("x") - }, + Oasis::Variable("x") }, Oasis::Real(-3) }; OASIS_CAPTURE_WITH_SERIALIZER(add); @@ -213,14 +203,12 @@ TEST_CASE("Quadratic polynomial test 5: x^2 - 25", "[factor]") TEST_CASE("Rational Quadratic polynomial test 1: 2x^2 + x - 1", "[factor]") { // 2x^2 + x - 1 - Oasis::Add<> add{ - Oasis::Multiply{ + Oasis::Add<> add { + Oasis::Multiply { Oasis::Real(2), - Oasis::Exponent{ - Oasis::Variable("x"), - Oasis::Real(2) - } - }, + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(2) } }, Oasis::Variable("x"), Oasis::Real(-1) }; @@ -236,34 +224,29 @@ TEST_CASE("Rational Quadratic polynomial test 1: 2x^2 + x - 1", "[factor]") auto denominator1 = root1->GetLeastSigOp().GetValue(); auto numerator1 = root1->GetMostSigOp().GetValue(); REQUIRE(root1 != nullptr); - REQUIRE(numerator1/denominator1 == 1.0/2.0); - + REQUIRE(numerator1 / denominator1 == 1.0 / 2.0); // Check second root (-1) auto root2 = Oasis::RecursiveCast>(*zeros[1]); auto denominator2 = root2->GetLeastSigOp().GetValue(); auto numerator2 = root2->GetMostSigOp().GetValue(); REQUIRE(root2 != nullptr); - REQUIRE(numerator2/denominator2 == -1); - + REQUIRE(numerator2 / denominator2 == -1); } } TEST_CASE("Rational Quadratic polynomial test 2: 6x^2 - 5x + 1", "[factor]") { // 6x^2 - 5x + 1 - Oasis::Add<> add{ - Oasis::Multiply{ + Oasis::Add<> add { + Oasis::Multiply { Oasis::Real(6), - Oasis::Exponent{ + Oasis::Exponent { Oasis::Variable("x"), - Oasis::Real(2) - } - }, - Oasis::Multiply{ - Oasis::Real(-5), - Oasis::Variable("x") - }, + Oasis::Real(2) } }, + Oasis::Multiply { + Oasis::Real(-5), + Oasis::Variable("x") }, Oasis::Real(1) }; OASIS_CAPTURE_WITH_SERIALIZER(add); @@ -278,40 +261,34 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x^2 - 5x + 1", "[factor]") auto denominator = root1->GetLeastSigOp().GetValue(); auto numerator = root1->GetMostSigOp().GetValue(); REQUIRE(root1 != nullptr); - REQUIRE(numerator/denominator == 1.0/2.0); + REQUIRE(numerator / denominator == 1.0 / 2.0); auto root2 = Oasis::RecursiveCast>(*zeros[1]); auto denominator2 = root2->GetLeastSigOp().GetValue(); auto numerator2 = root2->GetMostSigOp().GetValue(); REQUIRE(root2 != nullptr); - REQUIRE(numerator2/denominator2 == 1.0/3.0); + REQUIRE(numerator2 / denominator2 == 1.0 / 3.0); } } - TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") { // 3x^3 - 16x^2 + 23x - 6: - Oasis::Add<> cubic{ - Oasis::Multiply{ // 3x^3 term + Oasis::Add<> cubic { + Oasis::Multiply { // 3x^3 term Oasis::Real(3), - Oasis::Exponent{ + Oasis::Exponent { Oasis::Variable("x"), - Oasis::Real(3) - } - }, - Oasis::Multiply{ // -16x^2 term + Oasis::Real(3) } }, + Oasis::Multiply { // -16x^2 term Oasis::Real(-16), - Oasis::Exponent{ + Oasis::Exponent { Oasis::Variable("x"), - Oasis::Real(2) - } - }, - Oasis::Multiply{ // 23x term + Oasis::Real(2) } }, + Oasis::Multiply { // 23x term Oasis::Real(23), - Oasis::Variable("x") - }, - Oasis::Real(-6) // -6 term + Oasis::Variable("x") }, + Oasis::Real(-6) // -6 term }; OASIS_CAPTURE_WITH_SERIALIZER(cubic); @@ -324,18 +301,33 @@ TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") auto denominator = root1->GetLeastSigOp().GetValue(); auto numerator = root1->GetMostSigOp().GetValue(); REQUIRE(root1 != nullptr); - REQUIRE(numerator/denominator == 1.0/2.0); + REQUIRE(numerator / denominator == 1.0 / 2.0); auto root2 = Oasis::RecursiveCast>(*zeros[1]); auto denominator2 = root2->GetLeastSigOp().GetValue(); auto numerator2 = root2->GetMostSigOp().GetValue(); REQUIRE(root2 != nullptr); - REQUIRE(numerator2/denominator2 == 1.0/3.0); + REQUIRE(numerator2 / denominator2 == 1.0 / 3.0); auto root3 = Oasis::RecursiveCast>(*zeros[2]); auto denominator3 = root3->GetLeastSigOp().GetValue(); auto numerator3 = root3->GetMostSigOp().GetValue(); REQUIRE(root2 != nullptr); - REQUIRE(numerator3/denominator3 == 3.0/1.0); + REQUIRE(numerator3 / denominator3 == 3.0 / 1.0); } } + +TEST_CASE("Quartic test 1: x^4 - 1", "[factor]") +{ + Oasis::Add<> quartic { + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(4) }, + Oasis::Real(-1) + }; + + OASIS_CAPTURE_WITH_SERIALIZER(quartic); + auto zeros = quartic.FindZeros(); + + REQUIRE(zeros.size() == 2); +} From d92c12674422457689d99a1795ba19cda6059559 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 14 Mar 2025 17:08:34 -0400 Subject: [PATCH 18/21] simple quartic non-irrational, clang format new code --- src/Expression.cpp | 26 ++++++++++++++++---------- 1 file changed, 16 insertions(+), 10 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 89a5fc00..2f12db28 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -281,16 +281,21 @@ auto Expression::FindZeros() const -> std::vector> } } } - } else if (coefficents.size() == 5) { // Quartic equation ax⁴ + bx³ + cx² + dx + e = 0 + } else if (coefficents.size() == 5) { // Quartic equation ax⁴ + bx³ + cx² + dx + e = 0 // coefficients long long a = 0, b = 0, c = 0, d = 0, e = 0; // Convert coefficients to numbers if possible - if (auto aReal = RecursiveCast(*coefficents[4]); aReal != nullptr) a = static_cast(aReal->GetValue()); - if (auto bReal = RecursiveCast(*coefficents[3]); bReal != nullptr) b = static_cast(bReal->GetValue()); - if (auto cReal = RecursiveCast(*coefficents[2]); cReal != nullptr) c = static_cast(cReal->GetValue()); - if (auto dReal = RecursiveCast(*coefficents[1]); dReal != nullptr) d = static_cast(dReal->GetValue()); - if (auto eReal = RecursiveCast(*coefficents[0]); eReal != nullptr) e = static_cast(eReal->GetValue()); + if (auto aReal = RecursiveCast(*coefficents[4]); aReal != nullptr) + a = static_cast(aReal->GetValue()); + if (auto bReal = RecursiveCast(*coefficents[3]); bReal != nullptr) + b = static_cast(bReal->GetValue()); + if (auto cReal = RecursiveCast(*coefficents[2]); cReal != nullptr) + c = static_cast(cReal->GetValue()); + if (auto dReal = RecursiveCast(*coefficents[1]); dReal != nullptr) + d = static_cast(dReal->GetValue()); + if (auto eReal = RecursiveCast(*coefficents[0]); eReal != nullptr) + e = static_cast(eReal->GetValue()); // Find potential rational roots using the rational root theorem // Possible rational roots are p/q where: @@ -320,7 +325,8 @@ auto Expression::FindZeros() const -> std::vector> for (long long p : p_factors) { for (long long q : q_factors) { // Skip if q is 0 - if (q == 0) continue; + if (q == 0) + continue; // Simplify the fraction p/q long long g = std::gcd(std::abs(p), q); @@ -334,7 +340,7 @@ auto Expression::FindZeros() const -> std::vector> } // Check if we've already found this root - if (found_roots.find({num, den}) != found_roots.end()) { + if (found_roots.find({ num, den }) != found_roots.end()) { continue; } @@ -352,9 +358,9 @@ auto Expression::FindZeros() const -> std::vector> long long val = a * p4 + b * p3 * q + c * p2 * q2 + d * p * q3 + e * q4; - if (val == 0) { // If this is a root + if (val == 0) { // If this is a root results.push_back(Divide(Real(static_cast(num)), Real(static_cast(den))).Copy()); - found_roots.insert({num, den}); + found_roots.insert({ num, den }); } } } From 0e6098cf458ed96af3c2391a4fa270eedcae0dbf Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 21 Mar 2025 20:17:30 -0400 Subject: [PATCH 19/21] variable number check bug fixed, and working on return an expression as output --- include/Oasis/Expression.hpp | 9 +- src/Expression.cpp | 146 ++++++++++++++++++++++++++++--- tests/PolynomialTests.cpp | 162 +++++++++++++++++++++++++++++++---- 3 files changed, 288 insertions(+), 29 deletions(-) diff --git a/include/Oasis/Expression.hpp b/include/Oasis/Expression.hpp index c5f5c02a..5423f374 100644 --- a/include/Oasis/Expression.hpp +++ b/include/Oasis/Expression.hpp @@ -86,12 +86,19 @@ class Expression { */ [[nodiscard]] virtual auto Equals(const Expression& other) const -> bool = 0; + /** + * The FindExpressionRoots function; a helper function used in Findzeros; finds the root in terms of expression like a +- x/b + * + * @tparam origonalExpresion The expression for which all the factors will be found. + */ + auto FindExpressionRoots(const std::vector>& coefficients) + -> std::vector>; /** * 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. * * @tparam origonalExpresion The expression for which all the factors will be found. */ - auto FindZeros() const -> std::vector>; + auto FindZeros() const -> std::expected< std::vector>, std::string>; /** * Gets the category of this expression. diff --git a/src/Expression.cpp b/src/Expression.cpp index 2f12db28..e523dbd7 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include @@ -44,16 +45,83 @@ long long gcf(long long a, long long b) namespace Oasis { +// Helper function to find expression roots for polynomials +auto FindExpressionRoots(const std::vector>& coefficients) + -> std::vector> +{ + + std::vector> results; + + // Check degree of polynomial based on coefficients.size() + if (coefficients.size() == 3) { // Quadratic: ax² + bx + c + auto& a = coefficients[2]; + auto& b = coefficients[1]; + auto& c = coefficients[0]; + + // Calculate discriminant + auto bSquared = Multiply(*b, *b).Simplify(); + auto fourAC = Multiply(Real(4), Multiply(*a, *c)).Simplify(); + auto discriminant = Subtract(*bSquared, *fourAC).Simplify(); + + if (RecursiveCast(*discriminant) == nullptr) { + // Special cases + } + } else if (coefficients.size() == 4) { // Cubic: ax³ + bx² + cx + d + // cubic + } else if (coefficients.size() == 5) { // Quartic: ax⁴ + bx³ + cx² + dx + e + // quartic + } + + return results; +} + // currently only supports polynomials of one variable. /** * The FindZeros function finds all rational 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. * * @tparam origonalExpresion The expression for which all the factors will be found. */ -auto Expression::FindZeros() const -> std::vector> +auto Expression::FindZeros() const -> std::expected>, std::string> { std::vector> results; std::vector> termsE; + std::vector> termsVariableCheck; + std::string variName = ""; + + // Process terms to collect variables + if (auto addCase = RecursiveCast>(*this); addCase != nullptr) { + addCase->Flatten(termsVariableCheck); + } else { + termsVariableCheck.push_back(Copy()); + } + + // Check variables in the terms + for (const auto& i : termsVariableCheck) { + std::string variableName = ""; + + if (auto variableCase = RecursiveCast(*i); variableCase != nullptr) { + variableName = variableCase->GetName(); + } else if (auto expCase = RecursiveCast>(*i); expCase != nullptr) { + variableName = expCase->GetMostSigOp().GetName(); + } else if (auto prodCase = RecursiveCast>(*i); prodCase != nullptr) { + variableName = prodCase->GetLeastSigOp().GetName(); + } else if (auto prodExpCase = RecursiveCast>>(*i); prodExpCase != nullptr) { + variableName = prodExpCase->GetLeastSigOp().GetMostSigOp().GetName(); + } else if (auto divCase = RecursiveCast>(*i); divCase != nullptr) { + variableName = divCase->GetLeastSigOp().GetName(); + } else if (auto divExpCase = RecursiveCast>>(*i); divExpCase != nullptr) { + variableName = divExpCase->GetLeastSigOp().GetMostSigOp().GetName(); + } + + if (!variableName.empty()) { + if (variName.empty()) { + variName = variableName; + } else if (variName != variableName) { + return std::unexpected("Error: Polynomial only supports expressions with a single variable."); + } + } + } + if (auto subCase = RecursiveCast>(*this); subCase != nullptr) { // Check for x² - n pattern if (auto leftTerm = RecursiveCast>(subCase->GetMostSigOp()); @@ -191,7 +259,32 @@ auto Expression::FindZeros() const -> std::vector> results.push_back(Divide(Multiply(Real(-1), *coefficents[0]), *coefficents[1]).Simplify()); } } - } else if (coefficents.size() == 3) { // Quadratic equation ax + b + c = 0 + // } else if (coefficents.size() == 3) { // Quadratic equation ax + b + c = 0 + // auto& a = coefficents[2]; + // auto& b = coefficents[1]; + // auto& c = coefficents[0]; + // + // // Calculate discriminant + // auto bSquared = Multiply(*b, *b).Simplify(); + // auto fourAC = Multiply(Real(4), Multiply(*a, *c)).Simplify(); + // auto discriminant = Subtract(*bSquared, *fourAC).Simplify(); + // + // if (auto realDisc = RecursiveCast(*discriminant); + // realDisc != nullptr && realDisc->GetValue() >= 0) { + // + // auto negB = Multiply(Real(-1), *b).Simplify(); + // auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + // auto twoA = Multiply(Real(2), *a).Simplify(); + // + // // First, create the numerators for both roots + // auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); + // auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); + // + // // Now create the Divide expressions properly + // results.push_back(Divide(*numerator1, *twoA).Copy()); + // results.push_back(Divide(*numerator2, *twoA).Copy()); + // } + } else if (coefficents.size() == 3) { // Quadratic equation ax² + bx + c = 0 auto& a = coefficents[2]; auto& b = coefficents[1]; auto& c = coefficents[0]; @@ -204,19 +297,46 @@ auto Expression::FindZeros() const -> std::vector> if (auto realDisc = RecursiveCast(*discriminant); realDisc != nullptr && realDisc->GetValue() >= 0) { - auto negB = Multiply(Real(-1), *b).Simplify(); - auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); - auto twoA = Multiply(Real(2), *a).Simplify(); + double discValue = realDisc->GetValue(); + double sqrtDiscValue = std::sqrt(discValue); + + // Get coefficient values + double aVal = 0, bVal = 0; + if (auto aReal = RecursiveCast(*a); aReal != nullptr) { + aVal = aReal->GetValue(); + } + if (auto bReal = RecursiveCast(*b); bReal != nullptr) { + bVal = bReal->GetValue(); + } + + if (std::floor(sqrtDiscValue) == sqrtDiscValue) { + // Perfect square - roots will be rational + double root1 = (-bVal + sqrtDiscValue) / (2 * aVal); + double root2 = (-bVal - sqrtDiscValue) / (2 * aVal); + + // Check if roots are integers + if (std::floor(root1) == root1 && std::floor(root2) == root2) { + results.push_back(Divide(Real(root1), Real(1)).Copy()); + results.push_back(Divide(Real(root2), Real(1)).Copy()); + } else { + // Roots are fractions + results.push_back(Divide(Real(root1), Real(1)).Copy()); + results.push_back(Divide(Real(root2), Real(1)).Copy()); + } + } else { - // First, create the numerators for both roots - auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); - auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); + auto negB = Multiply(Real(-1), *b).Simplify(); + auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + auto twoA = Multiply(Real(2), *a).Simplify(); - // Now create the Divide expressions properly - results.push_back(Divide(*numerator1, *twoA).Copy()); - results.push_back(Divide(*numerator2, *twoA).Copy()); - // results.push_back(Divide(Add(*negB, *sqrtDisc), *twoA).Copy()); - // results.push_back(Divide(Subtract(*negB, *sqrtDisc), *twoA).Copy()); + // (-b + √discriminant)/(2a) + auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); + results.push_back(Divide(*numerator1, *twoA).Copy()); + + // (-b - √discriminant)/(2a) + auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); + results.push_back(Divide(*numerator2, *twoA).Copy()); + } } } else if (coefficents.size() == 4) { // Cubic equation: ax³ + bx² + cx + d = 0 long long a = termsC[0]; // coefficient of x³ diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 5a8ca246..6b70acb3 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -24,7 +24,11 @@ TEST_CASE("linear polynomial test 1: x + 30", "[factor]") Oasis::Variable("x"), Oasis::Real(30) }; - auto result = add.FindZeros(); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto result = std::move(result_wrapped.value()); + REQUIRE(result.size() == 1); REQUIRE(result[0]->Is()); @@ -32,6 +36,17 @@ TEST_CASE("linear polynomial test 1: x + 30", "[factor]") REQUIRE(poly_var_res.GetValue() == -30); } +TEST_CASE("Variable check test 1: x + y + 30", "[factor]") +{ + Oasis::Add<> add { + Oasis::Variable("x"), + Oasis::Variable("y"), + Oasis::Real(30) + }; + auto result = add.FindZeros(); + REQUIRE((!result.has_value())); +} + TEST_CASE("linear polynomial test 2: 3x - 6", "[factor]") { Oasis::Subtract minus { @@ -41,7 +56,10 @@ TEST_CASE("linear polynomial test 2: 3x - 6", "[factor]") Oasis::Real(6) }; OASIS_CAPTURE_WITH_SERIALIZER(minus); - auto result = minus.FindZeros(); + auto result_wrapped = minus.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto result = std::move(result_wrapped.value()); REQUIRE(result.size() == 1); REQUIRE(result[0]->Is()); @@ -57,11 +75,14 @@ TEST_CASE("linear polynomial test 3: 2x + 30", "[factor]") Oasis::Real(2), Oasis::Variable("x") } }; - auto zeros = add.FindZeros(); - REQUIRE(zeros.size() == 1); - REQUIRE(zeros[0]->Is()); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); - auto simplifiedReal = dynamic_cast(*zeros[0]); + auto result = std::move(result_wrapped.value()); + REQUIRE(result.size() == 1); + REQUIRE(result[0]->Is()); + + auto simplifiedReal = dynamic_cast(*result[0]); REQUIRE(simplifiedReal.GetValue() == -15); } @@ -79,7 +100,10 @@ TEST_CASE("Quadratic polynomial test 1: x^2 + 5x + 6", "[factor]") }; OASIS_CAPTURE_WITH_SERIALIZER(add); - auto zeros = add.FindZeros(); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); @@ -112,7 +136,10 @@ TEST_CASE("Quadratic polynomial test 2: x^2 - 2x -3", "[factor]") }; OASIS_CAPTURE_WITH_SERIALIZER(add); - auto zeros = add.FindZeros(); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); @@ -139,7 +166,10 @@ TEST_CASE("Quadratic polynomial test 3: x^2 - 9", "[factor]") Oasis::Real(2) }, Oasis::Real(9), }; - auto zeros = minus.FindZeros(); + auto result_wrapped = minus.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { @@ -162,7 +192,10 @@ TEST_CASE("Quadratic polynomial test 4: x^2 - 16", "[factor]") Oasis::Real(2) }, Oasis::Real(16), }; - auto zeros = minus.FindZeros(); + auto result_wrapped = minus.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { @@ -185,7 +218,10 @@ TEST_CASE("Quadratic polynomial test 5: x^2 - 25", "[factor]") Oasis::Real(2) }, Oasis::Real(25), }; - auto zeros = minus.FindZeros(); + auto result_wrapped = minus.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); OASIS_CAPTURE_WITH_SERIALIZER(minus); REQUIRE(zeros.size() == 2); if (zeros.size() == 2) { @@ -214,7 +250,10 @@ TEST_CASE("Rational Quadratic polynomial test 1: 2x^2 + x - 1", "[factor]") }; OASIS_CAPTURE_WITH_SERIALIZER(add); - auto zeros = add.FindZeros(); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); @@ -251,7 +290,10 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x^2 - 5x + 1", "[factor]") }; OASIS_CAPTURE_WITH_SERIALIZER(add); - auto zeros = add.FindZeros(); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); @@ -271,6 +313,41 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x^2 - 5x + 1", "[factor]") } } +TEST_CASE("Quadratic test 3: x^2 - 15x + 4", "[factor]") // with 4 roots +{ + Oasis::Add<> add { + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(2) }, + Oasis::Multiply { + Oasis::Real(-15), + Oasis::Variable("x") }, + Oasis::Real(4) + }; + + OASIS_CAPTURE_WITH_SERIALIZER(add); + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); + + REQUIRE(zeros.size() == 2); + if (zeros.size() == 2) { + + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + auto denominator = root1->GetLeastSigOp().GetValue(); + auto numerator = root1->GetMostSigOp().GetValue(); + REQUIRE(root1 != nullptr); + REQUIRE(numerator / denominator == 7.5 + sqrt(209) / 2); + + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + auto denominator2 = root2->GetLeastSigOp().GetValue(); + auto numerator2 = root2->GetMostSigOp().GetValue(); + REQUIRE(root2 != nullptr); + REQUIRE(numerator2 / denominator2 == 7.5 - sqrt(209) / 2); + } +} + TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") { // 3x^3 - 16x^2 + 23x - 6: @@ -292,7 +369,10 @@ TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") }; OASIS_CAPTURE_WITH_SERIALIZER(cubic); - auto zeros = cubic.FindZeros(); + auto result_wrapped = cubic.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 3); @@ -327,7 +407,59 @@ TEST_CASE("Quartic test 1: x^4 - 1", "[factor]") }; OASIS_CAPTURE_WITH_SERIALIZER(quartic); - auto zeros = quartic.FindZeros(); + auto result_wrapped = quartic.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); } + +TEST_CASE("Quartic test 2: x^4 - 5x^2 + 4", "[factor]") // with 4 roots +{ + Oasis::Add<> quartic { + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(4) }, + Oasis::Multiply { // -5x^2 term + Oasis::Real(-5), + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(2) } }, + Oasis::Real(4) + }; + + OASIS_CAPTURE_WITH_SERIALIZER(quartic); + auto result_wrapped = quartic.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); + + REQUIRE(zeros.size() == 4); + if (zeros.size() == 4) { + + auto root1 = Oasis::RecursiveCast>(*zeros[0]); + auto denominator = root1->GetLeastSigOp().GetValue(); + auto numerator = root1->GetMostSigOp().GetValue(); + REQUIRE(root1 != nullptr); + REQUIRE(numerator / denominator == 1.0 / 1.0); + + auto root2 = Oasis::RecursiveCast>(*zeros[1]); + auto denominator2 = root2->GetLeastSigOp().GetValue(); + auto numerator2 = root2->GetMostSigOp().GetValue(); + REQUIRE(root2 != nullptr); + REQUIRE(numerator2 / denominator2 == -1.0 / 1.0); + + auto root3 = Oasis::RecursiveCast>(*zeros[2]); + auto denominator3 = root3->GetLeastSigOp().GetValue(); + auto numerator3 = root3->GetMostSigOp().GetValue(); + REQUIRE(root3 != nullptr); + REQUIRE(numerator3 / denominator3 == 2.0 / 1.0); + + auto root4 = Oasis::RecursiveCast>(*zeros[3]); + auto denominator4 = root4->GetLeastSigOp().GetValue(); + auto numerator4 = root4->GetMostSigOp().GetValue(); + REQUIRE(root4 != nullptr); + REQUIRE(numerator4 / denominator4 == -2.0 / 1.0); + } +} From 8dc26be6c2a58d76950ba4b195257d7d55794424 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 28 Mar 2025 21:41:02 -0400 Subject: [PATCH 20/21] =?UTF-8?q?Cubic=20and=20quadratic=20can=20deal=20re?= =?UTF-8?q?turn=20the=20results=20in=20term=20of=20an=20expression=20like?= =?UTF-8?q?=20a=20+=20=E2=88=9Ab.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Expression.cpp | 294 +++++++++++++++++++++++++++++--------- tests/PolynomialTests.cpp | 102 +++++-------- 2 files changed, 266 insertions(+), 130 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index e523dbd7..189a1860 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -259,31 +259,6 @@ auto Expression::FindZeros() const -> std::expected(*discriminant); - // realDisc != nullptr && realDisc->GetValue() >= 0) { - // - // auto negB = Multiply(Real(-1), *b).Simplify(); - // auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); - // auto twoA = Multiply(Real(2), *a).Simplify(); - // - // // First, create the numerators for both roots - // auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); - // auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); - // - // // Now create the Divide expressions properly - // results.push_back(Divide(*numerator1, *twoA).Copy()); - // results.push_back(Divide(*numerator2, *twoA).Copy()); - // } } else if (coefficents.size() == 3) { // Quadratic equation ax² + bx + c = 0 auto& a = coefficents[2]; auto& b = coefficents[1]; @@ -300,57 +275,76 @@ auto Expression::FindZeros() const -> std::expectedGetValue(); double sqrtDiscValue = std::sqrt(discValue); - // Get coefficient values - double aVal = 0, bVal = 0; - if (auto aReal = RecursiveCast(*a); aReal != nullptr) { - aVal = aReal->GetValue(); - } - if (auto bReal = RecursiveCast(*b); bReal != nullptr) { - bVal = bReal->GetValue(); - } - + // Check if discriminant is a perfect square if (std::floor(sqrtDiscValue) == sqrtDiscValue) { - // Perfect square - roots will be rational + // Rational roots + double aVal = 0, bVal = 0; + if (auto aReal = RecursiveCast(*a); aReal != nullptr) + aVal = aReal->GetValue(); + if (auto bReal = RecursiveCast(*b); bReal != nullptr) + bVal = bReal->GetValue(); + double root1 = (-bVal + sqrtDiscValue) / (2 * aVal); double root2 = (-bVal - sqrtDiscValue) / (2 * aVal); - // Check if roots are integers - if (std::floor(root1) == root1 && std::floor(root2) == root2) { - results.push_back(Divide(Real(root1), Real(1)).Copy()); - results.push_back(Divide(Real(root2), Real(1)).Copy()); - } else { - // Roots are fractions - results.push_back(Divide(Real(root1), Real(1)).Copy()); - results.push_back(Divide(Real(root2), Real(1)).Copy()); - } + results.push_back(Divide(Real(root1), Real(1)).Copy()); + results.push_back(Divide(Real(root2), Real(1)).Copy()); } else { - - auto negB = Multiply(Real(-1), *b).Simplify(); - auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); - auto twoA = Multiply(Real(2), *a).Simplify(); - - // (-b + √discriminant)/(2a) - auto numerator1 = Add(*negB, *sqrtDisc).Simplify(); - results.push_back(Divide(*numerator1, *twoA).Copy()); - - // (-b - √discriminant)/(2a) - auto numerator2 = Subtract(*negB, *sqrtDisc).Simplify(); - results.push_back(Divide(*numerator2, *twoA).Copy()); + // Irrational roots, build the expression trees + + auto negB = Multiply(Real(-1), *b).Copy(); + auto twoA = Multiply(Real(2), *a).Copy(); + auto firstTerm = Divide(*negB, *twoA).Copy(); + + // √discriminant/(2a) term + auto sqrtDisc = std::make_unique>( + *discriminant, + Divide(Real(1), Real(2))); + auto secondTerm = Divide(*sqrtDisc, *twoA).Copy(); + + // Build the first root as an Add expression + auto root1 = std::make_unique>( + *firstTerm, + *secondTerm); + + // Build the second root as a Subtract expression + auto root2 = std::make_unique>( + *firstTerm, + *secondTerm); + + results.push_back(std::move(root1)); + results.push_back(std::move(root2)); } } } else if (coefficents.size() == 4) { // Cubic equation: ax³ + bx² + cx + d = 0 - long long a = termsC[0]; // coefficient of x³ - long long b = termsC[1]; // coefficient of x² - long long c = termsC[2]; // coefficient of x - long long d = termsC[3]; // constant term + // First try to find rational roots using the rational root theorem + long long a_ll = 0, b_ll = 0, c_ll = 0, d_ll = 0; + double a_d = 0, b_d = 0, c_d = 0; + + // Convert coefficients to numbers if possible + if (auto aReal = RecursiveCast(*coefficents[3]); aReal != nullptr) { + a_ll = static_cast(aReal->GetValue()); + a_d = aReal->GetValue(); + } + if (auto bReal = RecursiveCast(*coefficents[2]); bReal != nullptr) { + b_ll = static_cast(bReal->GetValue()); + b_d = bReal->GetValue(); + } + if (auto cReal = RecursiveCast(*coefficents[1]); cReal != nullptr) { + c_ll = static_cast(cReal->GetValue()); + c_d = cReal->GetValue(); + } + if (auto dReal = RecursiveCast(*coefficents[0]); dReal != nullptr) { + d_ll = static_cast(dReal->GetValue()); + } // To track found roots and avoid duplicates std::set> found_roots; // Get factors of constant term for possible p values std::vector p_factors; - for (long long i = 1; i <= std::abs(d); i++) { - if (d % i == 0) { + for (long long i = 1; i <= std::abs(d_ll); i++) { + if (d_ll % i == 0) { p_factors.push_back(i); p_factors.push_back(-i); } @@ -358,12 +352,13 @@ auto Expression::FindZeros() const -> std::expected q_factors; - for (long long i = 1; i <= std::abs(a); i++) { - if (a % i == 0) { + for (long long i = 1; i <= std::abs(a_ll); i++) { + if (a_ll % i == 0) { q_factors.push_back(i); } } + // Check potential rational roots for (long long p : p_factors) { for (long long q : q_factors) { // Skip if q is 0 @@ -392,15 +387,180 @@ auto Expression::FindZeros() const -> std::expected(num)), Real(static_cast(den))).Copy()); found_roots.insert({ num, den }); } } } + + if (results.size() == 3) { + return results; + } + + else if (!results.empty()) { + if (results.size() == 1) { + double root = 0; + if (auto divExpr = RecursiveCast>(*results[0]); divExpr != nullptr) { + if (auto numReal = RecursiveCast(divExpr->GetMostSigOp()); numReal != nullptr) { + if (auto denReal = RecursiveCast(divExpr->GetLeastSigOp()); denReal != nullptr) { + root = numReal->GetValue() / denReal->GetValue(); + } + } + } + + // Synthetic division to get a quadratic polynomial + double a2 = a_d; + double b2 = b_d + a2 * root; + double c2 = c_d + b2 * root; + + // Solve the quadratic equation a2x² + b2x + c2 = 0 + double discriminant = b2 * b2 - 4 * a2 * c2; + + if (discriminant >= 0) { + double sqrtDisc = std::sqrt(discriminant); + + // Check if we get integer or simple fraction roots + double root1 = (-b2 + sqrtDisc) / (2 * a2); + double root2 = (-b2 - sqrtDisc) / (2 * a2); + + // Check if roots are integers + if (std::floor(root1) == root1 && std::floor(root2) == root2) { + results.push_back(Divide(Real(root1), Real(1)).Copy()); + results.push_back(Divide(Real(root2), Real(1)).Copy()); + } else { + // Create symbolic expressions for the roots + auto negB = std::make_unique(-b2); + auto twoA = std::make_unique(2 * a2); + auto discExpr = std::make_unique(discriminant); + + // Create √discriminant expression + auto sqrtExpr = std::make_unique>( + *discExpr, + Divide(Real(1), Real(2))); + + // Create the first root expression: (-b + √discriminant)/(2a) + auto rootPlus = std::make_unique>( + Add(*negB, *sqrtExpr), + *twoA); + + // Create the second root expression: (-b - √discriminant)/(2a) + auto rootMinus = std::make_unique>( + Subtract(*negB, *sqrtExpr), + *twoA); + + results.push_back(std::move(rootPlus)); + results.push_back(std::move(rootMinus)); + } + } + } + return results; + } + + // First normalize the equation to the form x³ + px² + qx + r = 0 + auto& a_expr = coefficents[3]; + auto& b_expr = coefficents[2]; + auto& c_expr = coefficents[1]; + auto& d_expr = coefficents[0]; + + // Normalize by dividing by a + auto a_inv = Divide(Real(1), *a_expr).Copy(); + auto p_expr = Multiply(*b_expr, *a_inv).Copy(); + auto q_expr = Multiply(*c_expr, *a_inv).Copy(); + auto r_expr = Multiply(*d_expr, *a_inv).Copy(); + + // Convert to depressed cubic t³ + pt + q = 0 + // via substitution x = t - p/3 + + // Calculate p' = -p²/3 + q + auto p_squared = Multiply(*p_expr, *p_expr).Copy(); + auto term1 = Multiply(Real(-1.0 / 3.0), *p_squared).Copy(); + auto p_prime = Add(*term1, *q_expr).Copy(); + + // Calculate q' = 2p³/27 - pq/3 + r + auto p_cubed = Multiply(*p_expr, *p_squared).Copy(); + auto term2 = Multiply(Real(2.0 / 27.0), *p_cubed).Copy(); + auto pq = Multiply(*p_expr, *q_expr).Copy(); + auto term3 = Multiply(Real(-1.0 / 3.0), *pq).Copy(); + auto q_prime = Add(Add(*term2, *term3), *r_expr).Copy(); + + // Calculate the discriminant Δ = (q'/2)² + (p'/3)³ + auto q_prime_half = Multiply(Real(0.5), *q_prime).Copy(); + auto q_prime_half_squared = Multiply(*q_prime_half, *q_prime_half).Copy(); + + auto p_prime_third = Multiply(Real(1.0 / 3.0), *p_prime).Copy(); + auto p_prime_third_cubed = Exponent(*p_prime_third, Real(3)).Copy(); + + auto discriminant = Add(*q_prime_half_squared, *p_prime_third_cubed).Copy(); + + if (auto discReal = RecursiveCast(*discriminant); discReal != nullptr) { + double disc_val = discReal->GetValue(); + + // Create the cubic root expressions for u and v + // u = ∛(-q'/2 + √Δ) and v = ∛(-q'/2 - √Δ) + + auto neg_q_prime_half = Multiply(Real(-1), *q_prime_half).Copy(); + + if (disc_val > 0) { + // Create the square root of discriminant expression + auto sqrt_disc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + + // Create expressions for u and v + auto u_arg = Add(*neg_q_prime_half, *sqrt_disc).Copy(); + auto v_arg = Subtract(*neg_q_prime_half, *sqrt_disc).Copy(); + + auto u = Exponent(*u_arg, Divide(Real(1), Real(3))).Copy(); + auto v = Exponent(*v_arg, Divide(Real(1), Real(3))).Copy(); + + // The real root is u + v - p/3 + auto uv = Add(*u, *v).Copy(); + auto p_third = Multiply(Real(1.0 / 3.0), *p_expr).Copy(); + auto root = Subtract(*uv, *p_third).Copy(); + + results.push_back(std::move(root)); + } else if (disc_val == 0) { + + // First root is 3(q'/p') + auto three_q_p = Multiply(Real(3), Divide(*q_prime, *p_prime)).Copy(); + + // Second/third root is -3(q'/2p') + auto neg_three_q_2p = Multiply(Real(-3), Divide(*q_prime, Multiply(Real(2), *p_prime))).Copy(); + + results.push_back(std::move(three_q_p)); + results.push_back(std::move(neg_three_q_2p)); + results.push_back(Divide(*neg_three_q_2p, Real(1)).Copy()); // Duplicate the repeated root + } else { + // Trigonometric form + // Let p'/3 = -a² (a real) + // Let q'/2 = -a³cos(3θ) + // Now the roots are: 2a·cos(θ), 2a·cos(θ+2π/3), 2a·cos(θ+4π/3) + + auto cubic_root1 = std::make_unique>( + Multiply(Real(-1.0 / 3.0), *p_expr), // -p/3 + Exponent(*discriminant, Divide(Real(1), Real(3))) // ∛Δ + ); + + auto cubic_root2 = std::make_unique>( + Multiply(Real(-1.0 / 3.0), *p_expr), // -p/3 + Multiply( + Real(-0.5), // -1/2 + Exponent(*discriminant, Divide(Real(1), Real(3))) // ∛Δ + )); + + auto cubic_root3 = std::make_unique>( + Multiply(Real(-1.0 / 3.0), *p_expr), // -p/3 + Multiply( + Real(-0.5), // -1/2 + Exponent(*discriminant, Divide(Real(1), Real(3))) // ∛Δ + )); + + results.push_back(std::move(cubic_root1)); + results.push_back(std::move(cubic_root2)); + results.push_back(std::move(cubic_root3)); + } + } } else if (coefficents.size() == 5) { // Quartic equation ax⁴ + bx³ + cx² + dx + e = 0 // coefficients long long a = 0, b = 0, c = 0, d = 0, e = 0; diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 6b70acb3..45c57753 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -106,20 +106,9 @@ TEST_CASE("Quadratic polynomial test 1: x^2 + 5x + 6", "[factor]") auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); + REQUIRE(print_expr(*zeros[0]) == "(-2/1)"); + REQUIRE(print_expr(*zeros[1]) == "(-3/1)"); - if (zeros.size() == 2) { - // Check first root (-2) - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == -4); - REQUIRE(root1->GetLeastSigOp().GetValue() == 2); - - // Check second root (-3) - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - REQUIRE(root2 != nullptr); - REQUIRE(root2->GetMostSigOp().GetValue() == -6); - REQUIRE(root2->GetLeastSigOp().GetValue() == 2); - } } TEST_CASE("Quadratic polynomial test 2: x^2 - 2x -3", "[factor]") @@ -142,20 +131,8 @@ TEST_CASE("Quadratic polynomial test 2: x^2 - 2x -3", "[factor]") auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 2); - - if (zeros.size() == 2) { - // Check first root (3) - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - REQUIRE(root1 != nullptr); - REQUIRE(root1->GetMostSigOp().GetValue() == 6); - REQUIRE(root1->GetLeastSigOp().GetValue() == 2); - - // Check second root (-1) - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - REQUIRE(root2 != nullptr); - REQUIRE(root2->GetMostSigOp().GetValue() == -2); - REQUIRE(root2->GetLeastSigOp().GetValue() == 2); - } + REQUIRE(print_expr(*zeros[0]) == "(3/1)"); + REQUIRE(print_expr(*zeros[1]) == "(-1/1)"); } TEST_CASE("Quadratic polynomial test 3: x^2 - 9", "[factor]") @@ -313,7 +290,7 @@ TEST_CASE("Rational Quadratic polynomial test 2: 6x^2 - 5x + 1", "[factor]") } } -TEST_CASE("Quadratic test 3: x^2 - 15x + 4", "[factor]") // with 4 roots +TEST_CASE("Irrational Quadratic test 3: x^2 - 15x + 4", "[factor]") // roots of 15/2 +- sqrt(209) / 2); { Oasis::Add<> add { Oasis::Exponent { @@ -330,22 +307,9 @@ TEST_CASE("Quadratic test 3: x^2 - 15x + 4", "[factor]") // with 4 roots REQUIRE(result_wrapped.has_value()); auto zeros = std::move(result_wrapped.value()); - REQUIRE(zeros.size() == 2); - if (zeros.size() == 2) { - - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - auto denominator = root1->GetLeastSigOp().GetValue(); - auto numerator = root1->GetMostSigOp().GetValue(); - REQUIRE(root1 != nullptr); - REQUIRE(numerator / denominator == 7.5 + sqrt(209) / 2); - - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - auto denominator2 = root2->GetLeastSigOp().GetValue(); - auto numerator2 = root2->GetMostSigOp().GetValue(); - REQUIRE(root2 != nullptr); - REQUIRE(numerator2 / denominator2 == 7.5 - sqrt(209) / 2); - } + REQUIRE(print_expr(*zeros[0]) == "(((-1*-15)/(2*1))+((209^(1/2))/(2*1)))"); + REQUIRE(print_expr(*zeros[1]) == "(((-1*-15)/(2*1))-((209^(1/2))/(2*1)))"); } TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") @@ -375,26 +339,9 @@ TEST_CASE("Cubic polynomial test 1: 3x^3 - 16x^2 + 23x - 6:", "[factor]") auto zeros = std::move(result_wrapped.value()); REQUIRE(zeros.size() == 3); - - if (zeros.size() == 3) { - auto root1 = Oasis::RecursiveCast>(*zeros[0]); - auto denominator = root1->GetLeastSigOp().GetValue(); - auto numerator = root1->GetMostSigOp().GetValue(); - REQUIRE(root1 != nullptr); - REQUIRE(numerator / denominator == 1.0 / 2.0); - - auto root2 = Oasis::RecursiveCast>(*zeros[1]); - auto denominator2 = root2->GetLeastSigOp().GetValue(); - auto numerator2 = root2->GetMostSigOp().GetValue(); - REQUIRE(root2 != nullptr); - REQUIRE(numerator2 / denominator2 == 1.0 / 3.0); - - auto root3 = Oasis::RecursiveCast>(*zeros[2]); - auto denominator3 = root3->GetLeastSigOp().GetValue(); - auto numerator3 = root3->GetMostSigOp().GetValue(); - REQUIRE(root2 != nullptr); - REQUIRE(numerator3 / denominator3 == 3.0 / 1.0); - } + REQUIRE(print_expr(*zeros[0]) == "(1/3)"); + REQUIRE(print_expr(*zeros[1]) == "(2/1)"); + REQUIRE(print_expr(*zeros[2]) == "(3/1)"); } TEST_CASE("Quartic test 1: x^4 - 1", "[factor]") @@ -463,3 +410,32 @@ TEST_CASE("Quartic test 2: x^4 - 5x^2 + 4", "[factor]") // with 4 roots REQUIRE(numerator4 / denominator4 == -2.0 / 1.0); } } + +TEST_CASE("Irrational Cubic polynomial test 1: 4x^3 - 11x^2 + 2x + 3:", "[factor]") +{ + // 4x^3 - 11x^2 + 2x + 3: + Oasis::Add<> cubic { + Oasis::Multiply { + Oasis::Real(4), + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(3) } }, + Oasis::Multiply { + Oasis::Real(-11), + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(2) } }, + Oasis::Multiply { + Oasis::Real(2), + Oasis::Variable("x") }, + Oasis::Real(3) + }; + OASIS_CAPTURE_WITH_SERIALIZER(cubic); + + auto result_wrapped = cubic.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); + + REQUIRE(zeros.size() == 3); +} \ No newline at end of file From 2b1619ce48c72a7392cf4c7dffe1eaec2305bca9 Mon Sep 17 00:00:00 2001 From: lius24 <144491510+lius24@users.noreply.github.com> Date: Fri, 4 Apr 2025 19:30:24 -0400 Subject: [PATCH 21/21] Quartic complex expression roots --- src/Expression.cpp | 470 +++++++++++++++++++++++++++++++++----- tests/PolynomialTests.cpp | 27 ++- 2 files changed, 431 insertions(+), 66 deletions(-) diff --git a/src/Expression.cpp b/src/Expression.cpp index 189a1860..e32ddaf4 100644 --- a/src/Expression.cpp +++ b/src/Expression.cpp @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -561,90 +562,433 @@ auto Expression::FindZeros() const -> std::expected(*coefficents[4]); aReal != nullptr) - a = static_cast(aReal->GetValue()); - if (auto bReal = RecursiveCast(*coefficents[3]); bReal != nullptr) - b = static_cast(bReal->GetValue()); - if (auto cReal = RecursiveCast(*coefficents[2]); cReal != nullptr) - c = static_cast(cReal->GetValue()); - if (auto dReal = RecursiveCast(*coefficents[1]); dReal != nullptr) - d = static_cast(dReal->GetValue()); - if (auto eReal = RecursiveCast(*coefficents[0]); eReal != nullptr) - e = static_cast(eReal->GetValue()); - - // Find potential rational roots using the rational root theorem - // Possible rational roots are p/q where: - // - p is a factor of the constant term (e) - // - q is a factor of the leading coefficient (a) + if (auto aReal = RecursiveCast(*a); aReal != nullptr) a_ll = static_cast(aReal->GetValue()); + if (auto bReal = RecursiveCast(*b); bReal != nullptr) b_ll = static_cast(bReal->GetValue()); + if (auto cReal = RecursiveCast(*c); cReal != nullptr) c_ll = static_cast(cReal->GetValue()); + if (auto dReal = RecursiveCast(*d); dReal != nullptr) d_ll = static_cast(dReal->GetValue()); + if (auto eReal = RecursiveCast(*e); eReal != nullptr) e_ll = static_cast(eReal->GetValue()); - // Get factors of constant term for possible p values - std::vector p_factors; - for (long long i = 1; i <= std::abs(e); i++) { - if (e % i == 0) { - p_factors.push_back(i); - p_factors.push_back(-i); + // To track found roots and avoid duplicates + std::set> found_roots; + + // Check for rational roots only if coefficients are integers + if (a_ll != 0 && e_ll != 0) { + // Get factors of constant term for possible p values + std::vector p_factors; + for (long long i = 1; i <= std::abs(e_ll); i++) { + if (e_ll % i == 0) { + p_factors.push_back(i); + p_factors.push_back(-i); + } } - } - // Get factors of leading coefficient for possible q values - std::vector q_factors; - for (long long i = 1; i <= std::abs(a); i++) { - if (a % i == 0) { - q_factors.push_back(i); + // Get factors of leading coefficient for possible q values + std::vector q_factors; + for (long long i = 1; i <= std::abs(a_ll); i++) { + if (a_ll % i == 0) { + q_factors.push_back(i); + } } - } - // To track found roots and avoid duplicates - std::set> found_roots; + // Check potential rational roots + for (long long p : p_factors) { + for (long long q : q_factors) { + // Skip if q is 0 + if (q == 0) continue; + + // Simplify the fraction p/q + long long g = std::gcd(std::abs(p), q); + long long num = p / g; + long long den = q / g; + + // Ensure denominator is positive + if (den < 0) { + num = -num; + den = -den; + } - for (long long p : p_factors) { - for (long long q : q_factors) { - // Skip if q is 0 - if (q == 0) - continue; + // Check if we've already found this root + if (found_roots.find({ num, den }) != found_roots.end()) { + continue; + } - // Simplify the fraction p/q - long long g = std::gcd(std::abs(p), q); - long long num = p / g; - long long den = q / g; + // For each potential root p/q, evaluate the polynomial + double x = static_cast(num) / den; + double x2 = x * x; + double x3 = x2 * x; + double x4 = x3 * x; - // Ensure denominator is positive - if (den < 0) { - num = -num; - den = -den; + // Evaluate ax⁴ + bx³ + cx² + dx + e + double val = a_ll * x4 + b_ll * x3 + c_ll * x2 + d_ll * x + e_ll; + + if (std::abs(val) < 1e-10) { // If this is a root (with small tolerance) + results.push_back(Divide(Real(static_cast(num)), Real(static_cast(den))).Copy()); + found_roots.insert({ num, den }); + } } + } + } - // Check if we've already found this root - if (found_roots.find({ num, den }) != found_roots.end()) { - continue; + // If we found some rational roots, use polynomial division + if (results.size() > 0 && results.size() < 4) { + return results; + } + + // Check for biquadratic case (ax⁴ + cx² + e = 0) + bool isBiquadratic = false; + + if (auto bReal = RecursiveCast(*b); bReal != nullptr) { + if (auto dReal = RecursiveCast(*d); dReal != nullptr) { + if (std::abs(bReal->GetValue()) < 1e-10 && std::abs(dReal->GetValue()) < 1e-10) { + isBiquadratic = true; } + } + } - // Evaluate the polynomial at p/q using synthetic division - // For a quartic: a(p/q)⁴ + b(p/q)³ + c(p/q)² + d(p/q) + e + if (isBiquadratic) { + auto discriminant = Subtract( + Multiply(*c, *c), + Multiply(Real(4), Multiply(*a, *e)) + ).Simplify(); - // Multiply by q⁴ to clear denominators: - // a*p⁴ + b*p³*q + c*p²*q² + d*p*q³ + e*q⁴ - long long p2 = p * p; - long long p3 = p2 * p; - long long p4 = p3 * p; - long long q2 = q * q; - long long q3 = q2 * q; - long long q4 = q3 * q; + if (auto discReal = RecursiveCast(*discriminant); discReal != nullptr) { + double disc_val = discReal->GetValue(); - long long val = a * p4 + b * p3 * q + c * p2 * q2 + d * p * q3 + e * q4; + // Calculate -c/(2a) term which appears in both roots + auto negC = Multiply(Real(-1), *c).Copy(); + auto twoA = Multiply(Real(2), *a).Copy(); + auto commonTerm = Divide(*negC, *twoA).Copy(); - if (val == 0) { // If this is a root - results.push_back(Divide(Real(static_cast(num)), Real(static_cast(den))).Copy()); - found_roots.insert({ num, den }); + if (disc_val >= 0) { + // Real quadratic roots + + // Create √discriminant term + auto sqrtDisc = Exponent(*discriminant, Divide(Real(1), Real(2))).Copy(); + auto sqrtDiscOver2A = Divide(*sqrtDisc, *twoA).Copy(); + + // First root of the quadratic: y₁ = (-c + √disc)/(2a) + auto y1 = Add(*commonTerm, *sqrtDiscOver2A).Copy(); + + // Second root of the quadratic: y₂ = (-c - √disc)/(2a) + auto y2 = Subtract(*commonTerm, *sqrtDiscOver2A).Copy(); + + + // For y₁: + if (auto y1Real = RecursiveCast(*y1); y1Real != nullptr) { + double y1_val = y1Real->GetValue(); + + if (y1_val > 0) { + // y₁ > 0: Two real roots x = ±√y₁ + auto sqrt_y1 = Exponent(Real(y1_val), Divide(Real(1), Real(2))).Copy(); + auto neg_sqrt_y1 = Multiply(Real(-1), *sqrt_y1).Copy(); + + results.push_back(std::move(sqrt_y1)); + results.push_back(std::move(neg_sqrt_y1)); + } else if (y1_val < 0) { + // y₁ < 0: Two imaginary roots x = ±i√|y₁| + auto abs_y1 = Real(std::abs(y1_val)).Copy(); + auto sqrt_abs_y1 = Exponent(*abs_y1, Divide(Real(1), Real(2))).Copy(); + + auto img_root1 = Multiply(Imaginary(), *sqrt_abs_y1).Copy(); + auto img_root2 = Multiply(Real(-1), *img_root1).Copy(); + + results.push_back(std::move(img_root1)); + results.push_back(std::move(img_root2)); + } else { + // y₁ = 0: One root with multiplicity 2 + results.push_back(Real(0).Copy()); + } + } else { + // y₁ is a symbolic expression + auto sqrt_y1 = Exponent(*y1, Divide(Real(1), Real(2))).Copy(); + auto neg_sqrt_y1 = Multiply(Real(-1), *sqrt_y1).Copy(); + + results.push_back(std::move(sqrt_y1)); + results.push_back(std::move(neg_sqrt_y1)); + } + + // For y₂: + if (auto y2Real = RecursiveCast(*y2); y2Real != nullptr) { + double y2_val = y2Real->GetValue(); + + if (y2_val > 0) { + // y₂ > 0: Two real roots x = ±√y₂ + auto sqrt_y2 = Exponent(Real(y2_val), Divide(Real(1), Real(2))).Copy(); + auto neg_sqrt_y2 = Multiply(Real(-1), *sqrt_y2).Copy(); + + results.push_back(std::move(sqrt_y2)); + results.push_back(std::move(neg_sqrt_y2)); + } else if (y2_val < 0) { + // y₂ < 0: Two imaginary roots x = ±i√|y₂| + auto abs_y2 = Real(std::abs(y2_val)).Copy(); + auto sqrt_abs_y2 = Exponent(*abs_y2, Divide(Real(1), Real(2))).Copy(); + + auto img_root1 = Multiply(Imaginary(), *sqrt_abs_y2).Copy(); + auto img_root2 = Multiply(Real(-1), *img_root1).Copy(); + + results.push_back(std::move(img_root1)); + results.push_back(std::move(img_root2)); + } else { + // y₂ = 0: One root with multiplicity 2 (already added) + } + } else { + // y₂ is a symbolic expression + auto sqrt_y2 = Exponent(*y2, Divide(Real(1), Real(2))).Copy(); + auto neg_sqrt_y2 = Multiply(Real(-1), *sqrt_y2).Copy(); + + results.push_back(std::move(sqrt_y2)); + results.push_back(std::move(neg_sqrt_y2)); + } + } else { + + // Create √|discriminant| term + auto absDisc = Multiply(Real(-1), *discriminant).Copy(); + auto sqrtAbsDisc = Exponent(*absDisc, Divide(Real(1), Real(2))).Copy(); + auto imagTerm = Divide(*sqrtAbsDisc, *twoA).Copy(); + + // First root: -c/(2a) + i·√|disc|/(2a) + auto y1_imag = Multiply(Imaginary(), *imagTerm).Copy(); + auto y1 = Add(*commonTerm, *y1_imag).Copy(); + + // Second root: -c/(2a) - i·√|disc|/(2a) + auto y2_imag = Multiply(Real(-1), *y1_imag).Copy(); + auto y2 = Add(*commonTerm, *y2_imag).Copy(); + + std::cout << "Complex biquadratic roots not fully implemented yet." << std::endl; } } + + return results; + } + + if (results.empty()) { + // Convert to depressed quartic y⁴ + py² + qy + r + // via substitution y = x - b/(4a) + + // Get a⁻¹ for normalization + auto a_inv = Divide(Real(1), *a).Copy(); + + // Calculate the substitution term b/(4a) + auto b_over_4a = Divide(*b, Multiply(Real(4), *a)).Copy(); + + // Calculate coefficients of depressed quartic + + // p = -3b²/(8a²) + c/a + auto b_squared = Multiply(*b, *b).Copy(); + auto term1 = Multiply(Real(-3), Divide(*b_squared, Multiply(Real(8), Multiply(*a, *a)))).Copy(); + auto term2 = Multiply(*c, *a_inv).Copy(); + auto p = Add(*term1, *term2).Copy(); + + // q = b³/(8a³) - bc/(2a²) + d/a + auto b_cubed = Multiply(*b, *b_squared).Copy(); + auto term3 = Divide(*b_cubed, Multiply(Real(8), Exponent(*a, Real(3)))).Copy(); + auto bc = Multiply(*b, *c).Copy(); + auto term4 = Multiply(Real(-1), Divide(*bc, Multiply(Real(2), Multiply(*a, *a)))).Copy(); + auto term5 = Multiply(*d, *a_inv).Copy(); + auto q = Add(Add(*term3, *term4), *term5).Copy(); + + // r = -3b⁴/(256a⁴) + b²c/(16a³) - bd/(4a²) + e/a + auto b_fourth = Multiply(*b_squared, *b_squared).Copy(); + auto term6 = Multiply(Real(-3), Divide(*b_fourth, Multiply(Real(256), Exponent(*a, Real(4))))).Copy(); + auto b_squared_c = Multiply(*b_squared, *c).Copy(); + auto term7 = Divide(*b_squared_c, Multiply(Real(16), Exponent(*a, Real(3)))).Copy(); + auto bd = Multiply(*b, *d).Copy(); + auto term8 = Multiply(Real(-1), Divide(*bd, Multiply(Real(4), Multiply(*a, *a)))).Copy(); + auto term9 = Multiply(*e, *a_inv).Copy(); + auto r = Add(Add(Add(*term6, *term7), *term8), *term9).Copy(); + + // Create resolvent cubic z³ + 2pz² + (p² - 4r)z - q² + auto two_p = Multiply(Real(2), *p).Copy(); + auto p_squared = Multiply(*p, *p).Copy(); + auto four_r = Multiply(Real(4), *r).Copy(); + auto p_squared_minus_4r = Subtract(*p_squared, *four_r).Copy(); + auto q_squared = Multiply(*q, *q).Copy(); + auto neg_q_squared = Multiply(Real(-1), *q_squared).Copy(); + + // Solve the resolvent cubic to find a value of z + auto cubic_a = Real(1).Copy(); + auto cubic_b = two_p->Copy(); + auto cubic_c = p_squared_minus_4r->Copy(); + auto cubic_d = neg_q_squared->Copy(); + + std::vector> cubic_coeffs; + cubic_coeffs.push_back(std::move(cubic_d)); + cubic_coeffs.push_back(std::move(cubic_c)); + cubic_coeffs.push_back(std::move(cubic_b)); + cubic_coeffs.push_back(std::move(cubic_a)); + + auto cubic_expr = std::make_unique>( + Exponent{ Variable("x"), Real(3) }, + Multiply(Real(2), Multiply(*p, Exponent{ Variable("x"), Real(2) })), + Multiply(*p_squared_minus_4r, Variable("x")), + *neg_q_squared + ); + + + // Let's assume we found a root z (for demonstration purposes) + auto z = Real(1).Copy(); // Placeholder for an actual root + + // Step 3: Use z to factor the quartic into two quadratics + // y⁴ + py² + qy + r = (y² + u*y + v) * (y² - u*y + w) + // where u² = z, v + w = p + z, and v*w = r + + // Calculate u = ±√z + auto u = Exponent(*z, Divide(Real(1), Real(2))).Copy(); + + // v + w = p + z + // v*w = r + + auto p_plus_z = Add(*p, *z).Copy(); + auto neg_p_plus_z = Multiply(Real(-1), *p_plus_z).Copy(); + + // Calculate discriminant of this quadratic + auto quad_disc = Subtract( + Multiply(*p_plus_z, *p_plus_z), + Multiply(Real(4), *r) + ).Copy(); + + auto sqrt_quad_disc = Exponent(*quad_disc, Divide(Real(1), Real(2))).Copy(); + + // v = (-(p+z) + √((p+z)² - 4r))/2 + auto v = Divide(Add(*neg_p_plus_z, *sqrt_quad_disc), Real(2)).Copy(); + + // w = (-(p+z) - √((p+z)² - 4r))/2 + auto w = Divide(Subtract(*neg_p_plus_z, *sqrt_quad_disc), Real(2)).Copy(); + + // (y² + u*y + v) and (y² - u*y + w) + // First quadratic: y² + u*y + v = 0 + auto quad1_disc = Subtract( + Multiply(*u, *u), + Multiply(Real(4), *v) + ).Copy(); + + auto sqrt_quad1_disc = Exponent(*quad1_disc, Divide(Real(1), Real(2))).Copy(); + auto neg_u = Multiply(Real(-1), *u).Copy(); + + // Roots of first quadratic + auto y1 = Divide(Add(*neg_u, *sqrt_quad1_disc), Real(2)).Copy(); + auto y2 = Divide(Subtract(*neg_u, *sqrt_quad1_disc), Real(2)).Copy(); + + // Second quadratic: y² - u*y + w = 0 + auto quad2_disc = Subtract( + Multiply(*u, *u), + Multiply(Real(4), *w) + ).Copy(); + + auto sqrt_quad2_disc = Exponent(*quad2_disc, Divide(Real(1), Real(2))).Copy(); + + // Roots of second quadratic + auto y3 = Divide(Add(*u, *sqrt_quad2_disc), Real(2)).Copy(); + auto y4 = Divide(Subtract(*u, *sqrt_quad2_disc), Real(2)).Copy(); + + // Step 6: Convert back to original variable x = y + b/(4a) + auto x1 = Add(*y1, *b_over_4a).Copy(); + auto x2 = Add(*y2, *b_over_4a).Copy(); + auto x3 = Add(*y3, *b_over_4a).Copy(); + auto x4 = Add(*y4, *b_over_4a).Copy(); + + // Add all four roots to the results + results.push_back(std::move(x1)); + results.push_back(std::move(x2)); + results.push_back(std::move(x3)); + results.push_back(std::move(x4)); } } +// else if (coefficents.size() == 5) { // Quartic equation ax⁴ + bx³ + cx² + dx + e = 0 +// // coefficients +// long long a = 0, b = 0, c = 0, d = 0, e = 0; +// +// // Convert coefficients to numbers if possible +// if (auto aReal = RecursiveCast(*coefficents[4]); aReal != nullptr) +// a = static_cast(aReal->GetValue()); +// if (auto bReal = RecursiveCast(*coefficents[3]); bReal != nullptr) +// b = static_cast(bReal->GetValue()); +// if (auto cReal = RecursiveCast(*coefficents[2]); cReal != nullptr) +// c = static_cast(cReal->GetValue()); +// if (auto dReal = RecursiveCast(*coefficents[1]); dReal != nullptr) +// d = static_cast(dReal->GetValue()); +// if (auto eReal = RecursiveCast(*coefficents[0]); eReal != nullptr) +// e = static_cast(eReal->GetValue()); +// +// // Find potential rational roots using the rational root theorem +// // Possible rational roots are p/q where: +// // - p is a factor of the constant term (e) +// // - q is a factor of the leading coefficient (a) +// +// // Get factors of constant term for possible p values +// std::vector p_factors; +// for (long long i = 1; i <= std::abs(e); i++) { +// if (e % i == 0) { +// p_factors.push_back(i); +// p_factors.push_back(-i); +// } +// } +// +// // Get factors of leading coefficient for possible q values +// std::vector q_factors; +// for (long long i = 1; i <= std::abs(a); i++) { +// if (a % i == 0) { +// q_factors.push_back(i); +// } +// } +// +// // To track found roots and avoid duplicates +// std::set> found_roots; +// +// for (long long p : p_factors) { +// for (long long q : q_factors) { +// // Skip if q is 0 +// if (q == 0) +// continue; +// +// // Simplify the fraction p/q +// long long g = std::gcd(std::abs(p), q); +// long long num = p / g; +// long long den = q / g; +// +// // Ensure denominator is positive +// if (den < 0) { +// num = -num; +// den = -den; +// } +// +// // Check if we've already found this root +// if (found_roots.find({ num, den }) != found_roots.end()) { +// continue; +// } +// +// // Evaluate the polynomial at p/q using synthetic division +// // For a quartic: a(p/q)⁴ + b(p/q)³ + c(p/q)² + d(p/q) + e +// +// // Multiply by q⁴ to clear denominators: +// // a*p⁴ + b*p³*q + c*p²*q² + d*p*q³ + e*q⁴ +// long long p2 = p * p; +// long long p3 = p2 * p; +// long long p4 = p3 * p; +// long long q2 = q * q; +// long long q3 = q2 * q; +// long long q4 = q3 * q; +// +// long long val = a * p4 + b * p3 * q + c * p2 * q2 + d * p * q3 + e * q4; +// +// if (val == 0) { // If this is a root +// results.push_back(Divide(Real(static_cast(num)), Real(static_cast(den))).Copy()); +// found_roots.insert({ num, den }); +// } +// } +// } +// } } return results; } diff --git a/tests/PolynomialTests.cpp b/tests/PolynomialTests.cpp index 45c57753..ffdad378 100644 --- a/tests/PolynomialTests.cpp +++ b/tests/PolynomialTests.cpp @@ -135,6 +135,28 @@ TEST_CASE("Quadratic polynomial test 2: x^2 - 2x -3", "[factor]") REQUIRE(print_expr(*zeros[1]) == "(-1/1)"); } +TEST_CASE("Quadratic polynomial test 2: x^2 + 8x -14", "[factor]") +{ + // x^2 + 8x -14 + Oasis::Add<> add { + Oasis::Exponent { + Oasis::Variable("x"), + Oasis::Real(2) }, + Oasis::Multiply { + Oasis::Real(8), + Oasis::Variable("x") }, + Oasis::Real(-14) + }; + OASIS_CAPTURE_WITH_SERIALIZER(add); + + auto result_wrapped = add.FindZeros(); + REQUIRE(result_wrapped.has_value()); + + auto zeros = std::move(result_wrapped.value()); + + REQUIRE(zeros.size() == 2); +} + TEST_CASE("Quadratic polynomial test 3: x^2 - 9", "[factor]") { Oasis::Subtract minus { @@ -382,9 +404,8 @@ TEST_CASE("Quartic test 2: x^4 - 5x^2 + 4", "[factor]") // with 4 roots auto zeros = std::move(result_wrapped.value()); - REQUIRE(zeros.size() == 4); - if (zeros.size() == 4) { - + REQUIRE(zeros.size() == 8); + if (zeros.size() == 8) { auto root1 = Oasis::RecursiveCast>(*zeros[0]); auto denominator = root1->GetLeastSigOp().GetValue(); auto numerator = root1->GetMostSigOp().GetValue();