diff --git a/src/math/polynomial/algebraic_numbers.cpp b/src/math/polynomial/algebraic_numbers.cpp index bc654b7c4..ab70fc318 100644 --- a/src/math/polynomial/algebraic_numbers.cpp +++ b/src/math/polynomial/algebraic_numbers.cpp @@ -2028,6 +2028,20 @@ namespace algebraic_numbers { } IF_VERBOSE(9, verbose_stream() << "sturm 1\n"); + + // Check whether a can be separated from b's interval and vice versa + // this recognizes the case where the intervals overlap, + // but the anums do not lie in the intersection of the intervals. + mpq l_a, u_a, l_b, u_b; + to_mpq(qm(), la, l_a); + to_mpq(qm(), ua, u_a); + to_mpq(qm(), lb, l_b); + to_mpq(qm(), ub, u_b); + if (compare(cell_a, l_b) == sign_neg) return sign_neg; + if (compare(cell_a, u_b) == sign_pos) return sign_pos; + if (compare(cell_b, l_a) == sign_neg) return sign_pos; + if (compare(cell_b, u_a) == sign_pos) return sign_neg; + // // EXPENSIVE CASE // Let seq be the Sturm-Tarski sequence for diff --git a/src/test/algebraic_numbers.cpp b/src/test/algebraic_numbers.cpp index aed544714..9e4662359 100644 --- a/src/test/algebraic_numbers.cpp +++ b/src/test/algebraic_numbers.cpp @@ -104,6 +104,48 @@ void test_algebraic_comparison() { VERIFY(!am.eq(a, b)); // 2 != 3 } +void test_algebraic_comparison_edge_case() { + std::cout << "test_algebraic_comparison edge case\n"; + + // Let p1 = 1073741837 x^2 - 576460758745874510 x - 16106127555 + // Let p2 = p1 * (1073741837 x^2 - 576460759819616347 x -16106127555) + // = 1152921532524134569 x^4 - 1237940069261339757601884309 x^3 + // + 332307006992839334837849081482577900 x^2 + 18569101038920096364028264635 x + // + 259407344817930278025 + // Compare a = root(p1, 1) in (-8, 0) and b = root(p2, 2) in (-15/2^29, -7/2^28) + // The two numbers are different (a < b), but very close, and both are roots of p2 + + reslimit rl; + unsynch_mpq_manager qm; + anum_manager am(rl, qm); + manager m(rl, qm); + polynomial_ref x(m); + x = m.mk_polynomial(m.mk_var()); + + rational a0, a1, a2; + a0 = 161061; + a0 = (a0 * 100000) + 27555; + a1 = 576460758; + a1 = (a1 * 1000000000) + 745874510; + a2 = 10737; + a2 = (a2 * 100000) + 41837; + + rational b1; + b1 = 576460759; + b1 = (b1 * 1000000000) + 819616347; + + polynomial_ref p1(m); + polynomial_ref p2(m); + p1 = ((a2*x*x) - (a1*x)) - a0; + p2 = p1 * (((a2*x*x) - (b1*x)) - a0); + + scoped_anum a(am), b(am); + am.mk_root(p1, 1, a); + am.mk_root(p2, 2, b); + + VERIFY(!am.eq(a, b)); +} + void test_algebraic_degree() { std::cout << "test_algebraic_degree\n"; @@ -158,6 +200,7 @@ void test_algebraic_numbers() { test_algebraic_basic_operations(); test_algebraic_arithmetic(); test_algebraic_comparison(); + test_algebraic_comparison_edge_case(); test_algebraic_degree(); test_algebraic_signs(); }