From 8526ca22262a23f491dd4d102eb995458d619a7c Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 23 May 2016 18:27:20 +1000 Subject: [PATCH 01/34] Make Stein gcd more generic for polynomial. --- include/boost/math/common_factor_rt.hpp | 50 ++++++++++++++++----------------- include/boost/math/tools/polynomial.hpp | 49 ++++++++++++++++++++++++++++++-- test/test_polynomial.cpp | 12 +++++++- 3 files changed, 83 insertions(+), 28 deletions(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 1dda577b4..bcb88ef0b 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -65,6 +65,11 @@ namespace boost { return a < b; } + inline static void subtract(T& a, const T& b) + { + a -= b; + } + enum method_type { method_euclid = 0, @@ -84,19 +89,6 @@ namespace boost { template struct gcd_traits : public gcd_traits_defaults {}; // - // Special handling for polynomials: - // - namespace tools { - template - class polynomial; - } - - template - struct gcd_traits > : public gcd_traits_defaults - { - static const boost::math::tools::polynomial& abs(const boost::math::tools::polynomial& val) { return val; } - }; - // // Some platforms have fast bitscan operations, that allow us to implement // make_odd much more efficiently: // @@ -277,7 +269,7 @@ namespace detail shifts = std::min(gcd_traits::make_odd(u), gcd_traits::make_odd(v)); - while(gcd_traits::less(1, v)) + while(u != v) { u %= v; v -= u; @@ -290,7 +282,8 @@ namespace detail if(gcd_traits::less(u, v)) swap(u, v); } - return (v == 1 ? v : u) << shifts; + u <<= shifts; + return u; } /** Stein gcd (aka 'binary gcd') @@ -301,25 +294,32 @@ namespace detail SteinDomain Stein_gcd(SteinDomain m, SteinDomain n) { using std::swap; - BOOST_ASSERT(m >= 0); - BOOST_ASSERT(n >= 0); - if (m == SteinDomain(0)) + // BOOST_ASSERT(m >= 0); + // BOOST_ASSERT(n >= 0); + if (!m) return n; - if (n == SteinDomain(0)) + if (!n) return m; // m > 0 && n > 0 - int d_m = gcd_traits::make_odd(m); - int d_n = gcd_traits::make_odd(n); + unsigned shifts = std::min(gcd_traits::make_odd(m), gcd_traits::make_odd(n)); // odd(m) && odd(n) while (m != n) { - if (n > m) + if (gcd_traits::less(m, n)) swap(n, m); - m -= n; - gcd_traits::make_odd(m); + // n <= m (by degree) + gcd_traits::subtract(m, n); + // gcd_traits::make_odd(m); + if (!m) + { + n <<= shifts; + return n; + } + BOOST_ASSERT(even(m)); + do m >>= 1; while (even(m)); } // m == n - m <<= std::min(d_m, d_n); + m <<= shifts; return m; } diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index df1bde374..68151c704 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -444,8 +445,7 @@ class polynomial : template polynomial& operator >>=(U const &n) { - BOOST_ASSERT(n <= m_data.size()); - m_data.erase(m_data.begin(), m_data.begin() + n); + m_data.erase(m_data.begin(), m_data.begin() + std::min(size_type(n), m_data.size())); return *this; } @@ -696,6 +696,7 @@ polynomial pow(polynomial base, int exp) return result; } + template inline std::basic_ostream& operator << (std::basic_ostream& os, const polynomial& poly) { @@ -710,6 +711,50 @@ inline std::basic_ostream& operator << (std::basic_ostream +struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defaults< boost::math::tools::polynomial > +{ + static boost::math::tools::polynomial + abs(const boost::math::tools::polynomial &val) + { + return val.size() > 0 && val.data().back() < T(0) ? -val : val; + } + + inline static int make_odd(boost::math::tools::polynomial &x) + { + unsigned r = 0; + while (even(x)) + { + x >>= 1; + r++; + } + return r; + } + + inline static bool + less(boost::math::tools::polynomial const &a, boost::math::tools::polynomial const &b) + { + return a.size() < b.size(); + } + + inline static void + subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) + { + T ratio = a.data().front() / b.data().front(); + a -= ratio * b; + // normalize coefficients so that leading coefficient is whole + if (std::floor(a.data().back()) != a.data().back()) + { + ratio = T(1) / a.data().back(); + a *= ratio; + } + } +}; + } // namespace math } // namespace boost diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index 0ce45354a..b286d127e 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -161,8 +161,18 @@ BOOST_AUTO_TEST_CASE( test_gcd ) polynomial const u(d8.begin(), d8.end()); polynomial const v(d6.begin(), d6.end()); polynomial const w(d2.begin(), d2.end()); - polynomial const d = boost::math::gcd(u, v); + polynomial d = boost::math::gcd(u, v); BOOST_CHECK_EQUAL(w, d); + d = boost::math::detail::Stein_gcd(u, v); + BOOST_CHECK_EQUAL(w, d); + boost::array const s3 = {{-2, -3, 0, 1}}; + boost::array const s2 = {{-4, 0, 1}}; + boost::array const d_expected = {{-2, 1}}; + polynomial const x(s3.begin(), s3.end()); + polynomial const y(s2.begin(), s2.end()); + polynomial const z(d_expected.begin(), d_expected.end()); + d = boost::math::detail::Stein_gcd(x, y); + BOOST_CHECK_EQUAL(z, d); } // Sanity checks to make sure I didn't break it. From 579f33248ac50e950cceffd56e8b5ee00e88612a Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 20 Jun 2016 20:49:53 +1000 Subject: [PATCH 02/34] [gcd] Stein_gcd: assert that at least one of the operands is not zero. --- include/boost/math/common_factor_rt.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index bcb88ef0b..5434391d5 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -294,6 +294,7 @@ namespace detail SteinDomain Stein_gcd(SteinDomain m, SteinDomain n) { using std::swap; + BOOST_ASSERT(m || n); // BOOST_ASSERT(m >= 0); // BOOST_ASSERT(n >= 0); if (!m) From b3596f01580240aa13811590d5449654eb14c740 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 20 Jun 2016 20:50:05 +1000 Subject: [PATCH 03/34] [gcd] Add/update comments. --- include/boost/math/common_factor_rt.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 5434391d5..5fedcc3cf 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -295,13 +295,14 @@ namespace detail { using std::swap; BOOST_ASSERT(m || n); + // TODO: What is a negative polynomial??? // BOOST_ASSERT(m >= 0); // BOOST_ASSERT(n >= 0); if (!m) return n; if (!n) return m; - // m > 0 && n > 0 + // m > 0 && n > 0 ??? unsigned shifts = std::min(gcd_traits::make_odd(m), gcd_traits::make_odd(n)); // odd(m) && odd(n) while (m != n) From 58c17ae8f11eae3168f4c643f0e9e8478fc678ac Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 23 Jun 2016 11:14:41 +1000 Subject: [PATCH 04/34] [polynomial] Make the gcd_traits subtraction more debug friendly. --- include/boost/math/tools/polynomial.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index 68151c704..1ece4e846 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -744,13 +744,13 @@ struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defau inline static void subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) { - T ratio = a.data().front() / b.data().front(); - a -= ratio * b; + T const ratio_of_constants = a.data().front() / b.data().front(); + a -= ratio_of_constants * b; // normalize coefficients so that leading coefficient is whole if (std::floor(a.data().back()) != a.data().back()) { - ratio = T(1) / a.data().back(); - a *= ratio; + T const inv_leading_coeff = T(1) / a.data().back(); + a *= inv_leading_coeff; } } }; From a6b550a6f0ffeaff00b4826fda86e4101d71c2cd Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 23 Jun 2016 11:15:27 +1000 Subject: [PATCH 05/34] [polynomial] Make separate Euclid and Stein gcd test cases. --- test/test_polynomial.cpp | 49 ++++++++++++++++++++++++++++++++++-------------- 1 file changed, 35 insertions(+), 14 deletions(-) diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index b286d127e..281fc48db 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -146,7 +146,7 @@ BOOST_AUTO_TEST_CASE( test_division_over_ufd ) } -BOOST_AUTO_TEST_CASE( test_gcd ) +BOOST_AUTO_TEST_CASE( test_Euclidean_gcd ) { /* NOTE: Euclidean gcd is not yet customized to return THE greatest * common polynomial divisor. If d is THE greatest common divisior of u and @@ -157,24 +157,45 @@ BOOST_AUTO_TEST_CASE( test_gcd ) */ boost::array const d8 = {{105, 278, -88, -56, 16}}; boost::array const d6 = {{70, 232, -44, -64, 16}}; - boost::array const d2 = {{-35, 24, -4}}; - polynomial const u(d8.begin(), d8.end()); - polynomial const v(d6.begin(), d6.end()); - polynomial const w(d2.begin(), d2.end()); - polynomial d = boost::math::gcd(u, v); - BOOST_CHECK_EQUAL(w, d); - d = boost::math::detail::Stein_gcd(u, v); - BOOST_CHECK_EQUAL(w, d); + boost::array const d2 = {{35, -24, 4}}; + polynomial u(d8.begin(), d8.end()); + polynomial v(d6.begin(), d6.end()); + polynomial expected(d2.begin(), d2.end()); + polynomial result = boost::math::detail::Euclid_gcd(u, v); + BOOST_CHECK_EQUAL(result, expected); + boost::array const s3 = {{-2, -3, 0, 1}}; + boost::array const s2 = {{-4, 0, 1}}; + boost::array const d_expected = {{-2, 1}}; + u = polynomial(s3.begin(), s3.end()); + v = polynomial(s2.begin(), s2.end()); + expected = polynomial(d_expected.begin(), d_expected.end()); + result = boost::math::detail::Euclid_gcd(u, v); + BOOST_CHECK_EQUAL(result, expected); +} + + +BOOST_AUTO_TEST_CASE( test_Stein_gcd ) +{ + // Stein gcd is not quite right for polynomials yet. + boost::array const d8 = {{105, 278, -88, -56, 16}}; + boost::array const d6 = {{70, 232, -44, -64, 16}}; + boost::array const d2 = {{35, -24, 4}}; + polynomial u(d8.begin(), d8.end()); + polynomial v(d6.begin(), d6.end()); + polynomial expected(d2.begin(), d2.end()); + polynomial result = boost::math::detail::Stein_gcd(u, v); + BOOST_CHECK_EQUAL(result, expected); boost::array const s3 = {{-2, -3, 0, 1}}; boost::array const s2 = {{-4, 0, 1}}; boost::array const d_expected = {{-2, 1}}; - polynomial const x(s3.begin(), s3.end()); - polynomial const y(s2.begin(), s2.end()); - polynomial const z(d_expected.begin(), d_expected.end()); - d = boost::math::detail::Stein_gcd(x, y); - BOOST_CHECK_EQUAL(z, d); + u = polynomial(s3.begin(), s3.end()); + v = polynomial(s2.begin(), s2.end()); + expected = polynomial(d_expected.begin(), d_expected.end()); + result = boost::math::detail::Stein_gcd(u, v); + BOOST_CHECK_EQUAL(result, expected); } + // Sanity checks to make sure I didn't break it. typedef boost::mpl::list integral_test_types; typedef boost::mpl::list non_integral_test_types; From 3cb1a25a31c0e06ab8156350583538179eb7c8bd Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 23 Jun 2016 11:22:13 +1000 Subject: [PATCH 06/34] [gcd] Improve comments, use more asserts in Stein_gcd. --- include/boost/math/common_factor_rt.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 5fedcc3cf..155d07b37 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -309,15 +309,16 @@ namespace detail { if (gcd_traits::less(m, n)) swap(n, m); - // n <= m (by degree) + BOOST_ASSERT(n.degree() <= m.degree()); gcd_traits::subtract(m, n); - // gcd_traits::make_odd(m); + BOOST_ASSERT(even(m)); + // With polynomials, it is possible that m is now zero. if (!m) { n <<= shifts; return n; } - BOOST_ASSERT(even(m)); + // gcd_traits::make_odd(m); do m >>= 1; while (even(m)); } // m == n From 3caf63531a2d124f4b537172965c4dfaeed7dea1 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sat, 25 Jun 2016 16:33:47 +1000 Subject: [PATCH 07/34] [polynomial] Simplify gcd_traits::subtract. Define constant_coefficient() for readability. --- include/boost/math/tools/polynomial.hpp | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index 1ece4e846..c0bff5ff4 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -710,6 +710,13 @@ inline std::basic_ostream& operator << (std::basic_ostream +T constant_coefficient(polynomial const &x) +{ + return x ? x.data().front() : T(0); +} + } // namespace tools // @@ -744,14 +751,13 @@ struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defau inline static void subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) { - T const ratio_of_constants = a.data().front() / b.data().front(); - a -= ratio_of_constants * b; + using std::floor; + + T const r = constant_coefficient(a) / constant_coefficient(b); + a -= r * b; // normalize coefficients so that leading coefficient is whole - if (std::floor(a.data().back()) != a.data().back()) - { - T const inv_leading_coeff = T(1) / a.data().back(); - a *= inv_leading_coeff; - } + if (floor(a.data().back()) != a.data().back()) + a /= a.data().back(); } }; From d5074742ec410b736262b9a819ead2859a6b0986 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sat, 25 Jun 2016 16:34:21 +1000 Subject: [PATCH 08/34] [polynomial] Expand unit tests of Euclid and Stein gcd. --- test/test_polynomial.cpp | 85 ++++++++++++++++++++++++++++-------------------- 1 file changed, 50 insertions(+), 35 deletions(-) diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index 281fc48db..14b2be383 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -145,8 +145,16 @@ BOOST_AUTO_TEST_CASE( test_division_over_ufd ) BOOST_CHECK_EQUAL(aa % aa, zero); } +typedef boost::mpl::list integral_test_types; +typedef boost::mpl::list pod_floating_point_types; +typedef boost::mpl::list mp_floating_point_types; +typedef boost::mpl::list mp_non_integral_test_types; +typedef boost::mpl::joint_view floating_point_types; +typedef boost::mpl::joint_view non_integral_test_types; +typedef boost::mpl::joint_view all_test_types; -BOOST_AUTO_TEST_CASE( test_Euclidean_gcd ) + +BOOST_AUTO_TEST_CASE_TEMPLATE( test_Euclidean_gcd, T, floating_point_types) { /* NOTE: Euclidean gcd is not yet customized to return THE greatest * common polynomial divisor. If d is THE greatest common divisior of u and @@ -155,52 +163,59 @@ BOOST_AUTO_TEST_CASE( test_Euclidean_gcd ) * Alpha do. * This test is an example of the fact that it returns -d. */ - boost::array const d8 = {{105, 278, -88, -56, 16}}; - boost::array const d6 = {{70, 232, -44, -64, 16}}; - boost::array const d2 = {{35, -24, 4}}; - polynomial u(d8.begin(), d8.end()); - polynomial v(d6.begin(), d6.end()); - polynomial expected(d2.begin(), d2.end()); - polynomial result = boost::math::detail::Euclid_gcd(u, v); - BOOST_CHECK_EQUAL(result, expected); - boost::array const s3 = {{-2, -3, 0, 1}}; - boost::array const s2 = {{-4, 0, 1}}; - boost::array const d_expected = {{-2, 1}}; - u = polynomial(s3.begin(), s3.end()); - v = polynomial(s2.begin(), s2.end()); - expected = polynomial(d_expected.begin(), d_expected.end()); + boost::array const d8 = {{105, 278, -88, -56, 16}}; + boost::array const d6 = {{70, 232, -44, -64, 16}}; + boost::array const d2 = {{35, -24, 4}}; + polynomial u(d8.begin(), d8.end()); + polynomial v(d6.begin(), d6.end()); + polynomial expected(d2.begin(), d2.end()); + polynomial result; result = boost::math::detail::Euclid_gcd(u, v); - BOOST_CHECK_EQUAL(result, expected); + BOOST_CHECK(result == expected || result == expected * T(-1.0)); + result = boost::math::detail::Euclid_gcd(v, u); + BOOST_CHECK(result == expected || result == expected * T(-1.0)); + + boost::array const s3 = {{-2, -3, 0, 1}}; + boost::array const s2 = {{-4, 0, 1}}; + boost::array const d_expected = {{-2, 1}}; + u = polynomial(s3.begin(), s3.end()); + v = polynomial(s2.begin(), s2.end()); + expected = polynomial(d_expected.begin(), d_expected.end()); + result = boost::math::detail::Euclid_gcd(u, v); + BOOST_CHECK(result == expected || result == expected * T(-1.0)); + result = boost::math::detail::Euclid_gcd(v, u); + BOOST_CHECK(result == expected || result == expected * T(-1.0)); } -BOOST_AUTO_TEST_CASE( test_Stein_gcd ) +BOOST_AUTO_TEST_CASE_TEMPLATE( test_Stein_gcd, T, pod_floating_point_types ) { // Stein gcd is not quite right for polynomials yet. - boost::array const d8 = {{105, 278, -88, -56, 16}}; - boost::array const d6 = {{70, 232, -44, -64, 16}}; - boost::array const d2 = {{35, -24, 4}}; - polynomial u(d8.begin(), d8.end()); - polynomial v(d6.begin(), d6.end()); - polynomial expected(d2.begin(), d2.end()); - polynomial result = boost::math::detail::Stein_gcd(u, v); + boost::array const d8 = {{105, 278, -88, -56, 16}}; + boost::array const d6 = {{70, 232, -44, -64, 16}}; + boost::array const d2 = {{35, -24, 4}}; + polynomial u(d8.begin(), d8.end()); + polynomial v(d6.begin(), d6.end()); + polynomial expected(d2.begin(), d2.end()); + polynomial result; + result = boost::math::detail::Stein_gcd(u, v); + BOOST_CHECK_EQUAL(result, expected); + result = boost::math::detail::Stein_gcd(v, u); BOOST_CHECK_EQUAL(result, expected); - boost::array const s3 = {{-2, -3, 0, 1}}; - boost::array const s2 = {{-4, 0, 1}}; - boost::array const d_expected = {{-2, 1}}; - u = polynomial(s3.begin(), s3.end()); - v = polynomial(s2.begin(), s2.end()); - expected = polynomial(d_expected.begin(), d_expected.end()); + + boost::array const s3 = {{-2, -3, 0, 1}}; + boost::array const s2 = {{-4, 0, 1}}; + boost::array const d_expected = {{-2, 1}}; + u = polynomial(s3.begin(), s3.end()); + v = polynomial(s2.begin(), s2.end()); + expected = polynomial(d_expected.begin(), d_expected.end()); result = boost::math::detail::Stein_gcd(u, v); BOOST_CHECK_EQUAL(result, expected); + result = boost::math::detail::Stein_gcd(v, u); + BOOST_CHECK_EQUAL(result, expected); } -// Sanity checks to make sure I didn't break it. -typedef boost::mpl::list integral_test_types; -typedef boost::mpl::list non_integral_test_types; -typedef boost::mpl::joint_view all_test_types; - BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types ) { polynomial const a(d3a.begin(), d3a.end()); From 7e10edd571954c48946d35bee7b3cb87995af8a7 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 27 Jun 2016 22:05:42 +1000 Subject: [PATCH 09/34] [gcd] Tweak comments and assertions. --- include/boost/math/common_factor_rt.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 155d07b37..56bc67432 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -289,37 +289,37 @@ namespace detail /** Stein gcd (aka 'binary gcd') * * From Mathematics to Generic Programming, Alexander Stepanov, Daniel Rose + * + * What can we say about the Stein Domain? + * - It must have zero. + * - It must have parity: being exclusively even or odd. + * - It must have a smallest prime (2 for integers, x for polynomials, etc). + * - It must have shift operations that add or remove factors of the smallest prime. + * - It does not work for negative integers. */ template SteinDomain Stein_gcd(SteinDomain m, SteinDomain n) { using std::swap; BOOST_ASSERT(m || n); - // TODO: What is a negative polynomial??? - // BOOST_ASSERT(m >= 0); - // BOOST_ASSERT(n >= 0); if (!m) return n; if (!n) return m; - // m > 0 && n > 0 ??? unsigned shifts = std::min(gcd_traits::make_odd(m), gcd_traits::make_odd(n)); // odd(m) && odd(n) while (m != n) { if (gcd_traits::less(m, n)) swap(n, m); - BOOST_ASSERT(n.degree() <= m.degree()); gcd_traits::subtract(m, n); - BOOST_ASSERT(even(m)); - // With polynomials, it is possible that m is now zero. + // With polynomials for example, it is possible that m is now zero. if (!m) { n <<= shifts; return n; } - // gcd_traits::make_odd(m); - do m >>= 1; while (even(m)); + gcd_traits::make_odd(m); } // m == n m <<= shifts; From 577547a13755cb32d6105c9a067c59b4296ec88a Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 27 Jun 2016 22:06:48 +1000 Subject: [PATCH 10/34] [polynomial] Add Antoine Joux's implementation of subtraction for Stein. --- include/boost/math/tools/polynomial.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index c0bff5ff4..dd6fcf0bc 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -751,6 +751,8 @@ struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defau inline static void subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) { +#if 0 + // Stepanov's implementation; suffers from floating point inaccuracy. using std::floor; T const r = constant_coefficient(a) / constant_coefficient(b); @@ -758,6 +760,12 @@ struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defau // normalize coefficients so that leading coefficient is whole if (floor(a.data().back()) != a.data().back()) a /= a.data().back(); +#else + // Antoine Joux's implementation: huge coefficients. + T const tmp = constant_coefficient(a); + a *= constant_coefficient(b); + a -= tmp * b; +#endif } }; From 0df601679ea2ce4e99ed4ecd7d21340e7dd52fee Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Wed, 29 Jun 2016 16:57:02 +1000 Subject: [PATCH 11/34] More specific gcd_traits for polynomial... --- include/boost/math/tools/polynomial.hpp | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index dd6fcf0bc..839ede1e7 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -721,14 +721,16 @@ T constant_coefficient(polynomial const &x) // // Special handling for polynomials: +// Note that gcd_traits_polynomial is templated on the coefficient type T, +// not on polynomial. // template -struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defaults< boost::math::tools::polynomial > +struct gcd_traits_polynomial : public gcd_traits_defaults< boost::math::tools::polynomial > { - static boost::math::tools::polynomial + static boost::math::tools::polynomial const & abs(const boost::math::tools::polynomial &val) { - return val.size() > 0 && val.data().back() < T(0) ? -val : val; + return val; } inline static int make_odd(boost::math::tools::polynomial &x) @@ -768,11 +770,15 @@ struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_defau #endif } }; + + +template +struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_polynomial +{ +}; + } // namespace math } // namespace boost #endif // BOOST_MATH_TOOLS_POLYNOMIAL_HPP - - - From 97d0566689c035c2f882651aca8262c021f4cc3b Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 3 Jul 2016 22:06:18 +1000 Subject: [PATCH 12/34] [gcd] Add odd/even functions for integers, mostly for debugging. --- include/boost/math/common_factor_rt.hpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 56bc67432..b7c8c4d85 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -245,6 +245,23 @@ namespace boost { }; #endif +template +inline +typename enable_if_c::is_integer, bool>::type +odd(const T& x) +{ + return x & 1u; +} + +template +inline +typename enable_if_c::is_integer, bool>::type +even(const T& x) +{ + return !odd(x); +} + + namespace detail { @@ -313,6 +330,7 @@ namespace detail if (gcd_traits::less(m, n)) swap(n, m); gcd_traits::subtract(m, n); + BOOST_ASSERT(even(m)); // With polynomials for example, it is possible that m is now zero. if (!m) { From 11f3d04d025a164ed81dec5f346ddac8b1b8fed2 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 3 Jul 2016 22:08:07 +1000 Subject: [PATCH 13/34] [polynomial] Split gcd_traits on coefficient type: float vs int. Also adds negate() method for in-place negation. --- include/boost/math/tools/polynomial.hpp | 95 +++++++++++++++++++++++++++++++-- 1 file changed, 90 insertions(+), 5 deletions(-) diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index 839ede1e7..180f9ac93 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -491,6 +492,12 @@ class polynomial : using namespace boost::lambda; m_data.erase(std::find_if(m_data.rbegin(), m_data.rend(), _1 != T(0)).base(), m_data.end()); } + + // Negate in-place. + void negate() + { + std::transform(data().begin(), data().end(), data().begin(), std::negate()); + } private: template @@ -717,6 +724,14 @@ T constant_coefficient(polynomial const &x) return x ? x.data().front() : T(0); } +// Trivial but useful convenience function referred to simply as l() in Knuth. +template +T leading_coefficient(polynomial const &x) +{ + BOOST_ASSERT(x); + return x.data().back(); +} + } // namespace tools // @@ -724,8 +739,15 @@ T constant_coefficient(polynomial const &x) // Note that gcd_traits_polynomial is templated on the coefficient type T, // not on polynomial. // -template -struct gcd_traits_polynomial : public gcd_traits_defaults< boost::math::tools::polynomial > +template +struct gcd_traits_polynomial +{}; + + + +// gcd_traits for Z[x]. +template +struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > { static boost::math::tools::polynomial const & abs(const boost::math::tools::polynomial &val) @@ -733,6 +755,22 @@ struct gcd_traits_polynomial : public gcd_traits_defaults< boost::math::tools::p return val; } + inline static + void normalize(boost::math::tools::polynomial &x) + { + using boost::lambda::_1; + + // This does assume that the coefficients are totally ordered. + if (x) + { + if (leading_coefficient(x) < T(0)) + x.negate(); + // Find the first non-zero, because we can't do gcd(0, 0). + T const d = gcd_range(find_if(x.data().begin(), x.data().end(), _1 != T(0)), x.data().end()).first; + x /= d; + } + } + inline static int make_odd(boost::math::tools::polynomial &x) { unsigned r = 0; @@ -753,20 +791,67 @@ struct gcd_traits_polynomial : public gcd_traits_defaults< boost::math::tools::p inline static void subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) { -#if 0 + // Antoine Joux's implementation: huge coefficients. + T const tmp = constant_coefficient(a); + a *= constant_coefficient(b); + a -= tmp * b; + normalize(a); + } +}; + + +template +struct gcd_traits_polynomial >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > +{ + static boost::math::tools::polynomial const & + abs(const boost::math::tools::polynomial &val) + { + return val; + } + + + inline static int make_odd(boost::math::tools::polynomial &x) + { + unsigned r = 0; + while (even(x)) + { + x >>= 1; + r++; + } + return r; + } + + inline static bool + less(boost::math::tools::polynomial const &a, boost::math::tools::polynomial const &b) + { + return a.size() < b.size(); + } + + inline static void + subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) + { +#if 1 // Stepanov's implementation; suffers from floating point inaccuracy. using std::floor; - + T const r = constant_coefficient(a) / constant_coefficient(b); a -= r * b; // normalize coefficients so that leading coefficient is whole - if (floor(a.data().back()) != a.data().back()) + if (a && floor(a.data().back()) != a.data().back()) a /= a.data().back(); #else + using boost::lambda::_1; + // Antoine Joux's implementation: huge coefficients. T const tmp = constant_coefficient(a); a *= constant_coefficient(b); a -= tmp * b; + if (std::numeric_limits::has_infinity) + { + T const &inf(std::numeric_limits::infinity()); + if (std::find_if(a.data().begin(), a.data().end(), _1 == inf || _1 == -inf) != a.data().end()) + throw std::domain_error("Floating point overflow."); + } #endif } }; From 54a8835e16cbd934b2daf0b25c95483c763d9197 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 3 Jul 2016 22:08:40 +1000 Subject: [PATCH 14/34] [polynomial] Make the Stein_gcd unit tests more organized. --- test/test_polynomial.cpp | 154 ++++++++++++++++++++++++++++++++++++----------- 1 file changed, 119 insertions(+), 35 deletions(-) diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index 14b2be383..a9334ff0f 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -15,10 +15,16 @@ #include #include #include +#include + #include +using namespace boost::math; using namespace boost::math::tools; using namespace std; +using boost::math::detail::Stein_gcd; +using boost::math::detail::Euclid_gcd; + template struct answer @@ -145,13 +151,119 @@ BOOST_AUTO_TEST_CASE( test_division_over_ufd ) BOOST_CHECK_EQUAL(aa % aa, zero); } -typedef boost::mpl::list integral_test_types; -typedef boost::mpl::list pod_floating_point_types; +typedef boost::mpl::list pod_integral_test_types; +typedef boost::mpl::list pod_floating_point_types; +typedef boost::mpl::joint_view pod_test_types; typedef boost::mpl::list mp_floating_point_types; typedef boost::mpl::list mp_non_integral_test_types; typedef boost::mpl::joint_view floating_point_types; typedef boost::mpl::joint_view non_integral_test_types; -typedef boost::mpl::joint_view all_test_types; +typedef boost::mpl::joint_view all_test_types; + +// Test data. + +typedef boost::tuple< boost::array, boost::array, boost::array > test_datum; + +std::vector make_test_data() +{ + std::vector< test_datum > test_data; + boost::array x_data = {{105, 278, -88, -56, 16}}; + boost::array y_data = {{70, 232, -44, -64, 16}}; + boost::array z_data = {{35, -24, 4}}; + test_data.push_back(boost::make_tuple(x_data, y_data, z_data)); + + return test_data; +} + +std::vector foo = make_test_data(); + +template +struct FM2GP_Ex_8_3__1 +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_Ex_8_3__1() + { + boost::array const x_data = {{105, 278, -88, -56, 16}}; + boost::array const y_data = {{70, 232, -44, -64, 16}}; + boost::array const z_data = {{35, -24, 4}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + +template +struct FM2GP_Ex_8_3__2 +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_Ex_8_3__2() + { + boost::array const x_data = {{1, -6, -8, 6, 7}}; + boost::array const y_data = {{1, -5, -2, 15, 11}}; + boost::array const z_data = {{1, 2, 1}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + + +template +struct FM2GP_trivial +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_trivial() + { + boost::array const x_data = {{-2, -3, 0, 1}}; + boost::array const y_data = {{-4, 0, 1}}; + boost::array const z_data= {{-2, 1}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_1, T, pod_test_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = Stein_gcd(fixture_type::x, fixture_type::y); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Stein_gcd(fixture_type::y, fixture_type::x); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_2, T, pod_test_types, FM2GP_Ex_8_3__2 ) +{ + typedef FM2GP_Ex_8_3__2 fixture_type; + polynomial w; + w = Stein_gcd(fixture_type::x, fixture_type::y); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Stein_gcd(fixture_type::y, fixture_type::x); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_trivial, T, pod_test_types, FM2GP_trivial ) +{ + typedef FM2GP_trivial fixture_type; + polynomial w; + w = Stein_gcd(fixture_type::x, fixture_type::y); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Stein_gcd(fixture_type::y, fixture_type::x); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} BOOST_AUTO_TEST_CASE_TEMPLATE( test_Euclidean_gcd, T, floating_point_types) @@ -170,9 +282,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( test_Euclidean_gcd, T, floating_point_types) polynomial v(d6.begin(), d6.end()); polynomial expected(d2.begin(), d2.end()); polynomial result; - result = boost::math::detail::Euclid_gcd(u, v); + result = Euclid_gcd(u, v); BOOST_CHECK(result == expected || result == expected * T(-1.0)); - result = boost::math::detail::Euclid_gcd(v, u); + result = Euclid_gcd(v, u); BOOST_CHECK(result == expected || result == expected * T(-1.0)); boost::array const s3 = {{-2, -3, 0, 1}}; @@ -181,41 +293,13 @@ BOOST_AUTO_TEST_CASE_TEMPLATE( test_Euclidean_gcd, T, floating_point_types) u = polynomial(s3.begin(), s3.end()); v = polynomial(s2.begin(), s2.end()); expected = polynomial(d_expected.begin(), d_expected.end()); - result = boost::math::detail::Euclid_gcd(u, v); + result = Euclid_gcd(u, v); BOOST_CHECK(result == expected || result == expected * T(-1.0)); - result = boost::math::detail::Euclid_gcd(v, u); + result = Euclid_gcd(v, u); BOOST_CHECK(result == expected || result == expected * T(-1.0)); } -BOOST_AUTO_TEST_CASE_TEMPLATE( test_Stein_gcd, T, pod_floating_point_types ) -{ - // Stein gcd is not quite right for polynomials yet. - boost::array const d8 = {{105, 278, -88, -56, 16}}; - boost::array const d6 = {{70, 232, -44, -64, 16}}; - boost::array const d2 = {{35, -24, 4}}; - polynomial u(d8.begin(), d8.end()); - polynomial v(d6.begin(), d6.end()); - polynomial expected(d2.begin(), d2.end()); - polynomial result; - result = boost::math::detail::Stein_gcd(u, v); - BOOST_CHECK_EQUAL(result, expected); - result = boost::math::detail::Stein_gcd(v, u); - BOOST_CHECK_EQUAL(result, expected); - - boost::array const s3 = {{-2, -3, 0, 1}}; - boost::array const s2 = {{-4, 0, 1}}; - boost::array const d_expected = {{-2, 1}}; - u = polynomial(s3.begin(), s3.end()); - v = polynomial(s2.begin(), s2.end()); - expected = polynomial(d_expected.begin(), d_expected.end()); - result = boost::math::detail::Stein_gcd(u, v); - BOOST_CHECK_EQUAL(result, expected); - result = boost::math::detail::Stein_gcd(v, u); - BOOST_CHECK_EQUAL(result, expected); -} - - BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types ) { polynomial const a(d3a.begin(), d3a.end()); From aff5c078bcc924cc57ed387ca6d4d1c7d7e1d23f Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Wed, 6 Jul 2016 11:00:58 +1000 Subject: [PATCH 15/34] [polynomial] Move gcd_traits into new file. --- include/boost/math/tools/polynomial.hpp | 130 ----------------------- include/boost/math/tools/polynomial_gcd.hpp | 159 ++++++++++++++++++++++++++++ test/test_polynomial.cpp | 1 + 3 files changed, 160 insertions(+), 130 deletions(-) create mode 100644 include/boost/math/tools/polynomial_gcd.hpp diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index 180f9ac93..d7c90729b 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -733,136 +733,6 @@ T leading_coefficient(polynomial const &x) } } // namespace tools - -// -// Special handling for polynomials: -// Note that gcd_traits_polynomial is templated on the coefficient type T, -// not on polynomial. -// -template -struct gcd_traits_polynomial -{}; - - - -// gcd_traits for Z[x]. -template -struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > -{ - static boost::math::tools::polynomial const & - abs(const boost::math::tools::polynomial &val) - { - return val; - } - - inline static - void normalize(boost::math::tools::polynomial &x) - { - using boost::lambda::_1; - - // This does assume that the coefficients are totally ordered. - if (x) - { - if (leading_coefficient(x) < T(0)) - x.negate(); - // Find the first non-zero, because we can't do gcd(0, 0). - T const d = gcd_range(find_if(x.data().begin(), x.data().end(), _1 != T(0)), x.data().end()).first; - x /= d; - } - } - - inline static int make_odd(boost::math::tools::polynomial &x) - { - unsigned r = 0; - while (even(x)) - { - x >>= 1; - r++; - } - return r; - } - - inline static bool - less(boost::math::tools::polynomial const &a, boost::math::tools::polynomial const &b) - { - return a.size() < b.size(); - } - - inline static void - subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) - { - // Antoine Joux's implementation: huge coefficients. - T const tmp = constant_coefficient(a); - a *= constant_coefficient(b); - a -= tmp * b; - normalize(a); - } -}; - - -template -struct gcd_traits_polynomial >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > -{ - static boost::math::tools::polynomial const & - abs(const boost::math::tools::polynomial &val) - { - return val; - } - - - inline static int make_odd(boost::math::tools::polynomial &x) - { - unsigned r = 0; - while (even(x)) - { - x >>= 1; - r++; - } - return r; - } - - inline static bool - less(boost::math::tools::polynomial const &a, boost::math::tools::polynomial const &b) - { - return a.size() < b.size(); - } - - inline static void - subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) - { -#if 1 - // Stepanov's implementation; suffers from floating point inaccuracy. - using std::floor; - - T const r = constant_coefficient(a) / constant_coefficient(b); - a -= r * b; - // normalize coefficients so that leading coefficient is whole - if (a && floor(a.data().back()) != a.data().back()) - a /= a.data().back(); -#else - using boost::lambda::_1; - - // Antoine Joux's implementation: huge coefficients. - T const tmp = constant_coefficient(a); - a *= constant_coefficient(b); - a -= tmp * b; - if (std::numeric_limits::has_infinity) - { - T const &inf(std::numeric_limits::infinity()); - if (std::find_if(a.data().begin(), a.data().end(), _1 == inf || _1 == -inf) != a.data().end()) - throw std::domain_error("Floating point overflow."); - } -#endif - } -}; - - -template -struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_polynomial -{ - -}; - } // namespace math } // namespace boost diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp new file mode 100644 index 000000000..52644dd48 --- /dev/null +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -0,0 +1,159 @@ +// (C) Copyright John Maddock and Jeremy William Murphy 2016. + +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP +#define BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP + +#ifdef _MSC_VER +#pragma once +#endif + +#include +#include + +namespace boost { namespace math { + +// +// Special handling for polynomials: +// Note that gcd_traits_polynomial is templated on the coefficient type T, +// not on polynomial. +// +template +struct gcd_traits_polynomial +{}; + + +// gcd_traits for Z[x]. +template +struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > +{ + typedef boost::math::tools::polynomial polynomial_type; + + // Just do a no-op for abs() as we'll normalize() later. + static polynomial_type + abs(const polynomial_type &val) + { + return leading_coefficient(val) < T(0) ? -val : val; + } + + inline static + void normalize(polynomial_type &x) + { + using boost::lambda::_1; + + // This does assume that the coefficients are totally ordered. + if (x) + { + if (leading_coefficient(x) < T(0)) + x.negate(); + // Find the first non-zero, because we can't do gcd(0, 0). + T const d = gcd_range(find_if(x.data().begin(), x.data().end(), _1 != T(0)), x.data().end()).first; + x /= d; + } + } + + inline static int make_odd(boost::math::tools::polynomial &x) + { + unsigned r = 0; + while (even(x)) + { + x >>= 1; + r++; + } + return r; + } + + inline static bool + less(boost::math::tools::polynomial const &a, boost::math::tools::polynomial const &b) + { + return a.size() < b.size(); + } + + inline static void + subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) + { + // Antoine Joux's implementation: huge coefficients. + T const tmp = constant_coefficient(a); + a *= constant_coefficient(b); + a -= tmp * b; + normalize(a); + } +}; + + +// TODO: Enable for multiprecision floating point types. +template +struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > +{ + typedef boost::math::tools::polynomial polynomial_type; + typedef gcd_traits_defaults< boost::math::tools::polynomial > parent_class; + using typename parent_class::method_type; + static const method_type method = method_type::method_euclid; + + static polynomial_type + abs(const polynomial_type &val) + { + return leading_coefficient(val) < T(0) ? -val : val; + } + + + inline static int make_odd(polynomial_type &x) + { + unsigned r = 0; + while (even(x)) + { + x >>= 1; + r++; + } + return r; + } + + inline static bool + less(polynomial_type const &a, polynomial_type const &b) + { + return a.size() < b.size(); + } + + inline static void + subtract(polynomial_type &a, polynomial_type const &b) + { +#if 1 + // Stepanov's implementation; suffers from floating point inaccuracy. + using std::floor; + + T const r = constant_coefficient(a) / constant_coefficient(b); + a -= r * b; + // normalize coefficients so that leading coefficient is whole + if (a && floor(a.data().back()) != a.data().back()) + a /= a.data().back(); +#else + using boost::lambda::_1; + + // Antoine Joux's implementation: huge coefficients. + T const tmp = constant_coefficient(a); + a *= constant_coefficient(b); + a -= tmp * b; + if (std::numeric_limits::has_infinity) + { + T const &inf(std::numeric_limits::infinity()); + if (std::find_if(a.data().begin(), a.data().end(), _1 == inf || _1 == -inf) != a.data().end()) + throw std::domain_error("Floating point overflow."); + } +#endif + } +}; + + +template +struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_polynomial +{ + +}; + +} // namespace math +} // namespace boost + +#endif // BOOST_MATH_TOOLS_POLYNOMIAL_GCD_HPP diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index a9334ff0f..d996c9814 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -7,6 +7,7 @@ #define BOOST_TEST_MAIN #include #include +#include #include #include #include From ca97a3b1249f575c3cfa255e912876b774f4e1e8 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Wed, 6 Jul 2016 11:01:41 +1000 Subject: [PATCH 16/34] Redesign gcd unit tests to use fixtures. --- test/test_polynomial.cpp | 183 +++++++++++++++++++++++++++++++++++++---------- 1 file changed, 146 insertions(+), 37 deletions(-) diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index d996c9814..3d9466e0a 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -153,16 +153,17 @@ BOOST_AUTO_TEST_CASE( test_division_over_ufd ) } typedef boost::mpl::list pod_integral_test_types; -typedef boost::mpl::list pod_floating_point_types; +typedef boost::mpl::list integral_test_types; +typedef boost::mpl::list pod_floating_point_types; typedef boost::mpl::joint_view pod_test_types; typedef boost::mpl::list mp_floating_point_types; typedef boost::mpl::list mp_non_integral_test_types; typedef boost::mpl::joint_view floating_point_types; typedef boost::mpl::joint_view non_integral_test_types; -typedef boost::mpl::joint_view all_test_types; +typedef boost::mpl::joint_view all_test_types; // Test data. - +/* typedef boost::tuple< boost::array, boost::array, boost::array > test_datum; std::vector make_test_data() @@ -177,6 +178,7 @@ std::vector make_test_data() } std::vector foo = make_test_data(); +*/ template struct FM2GP_Ex_8_3__1 @@ -216,6 +218,25 @@ struct FM2GP_Ex_8_3__2 template +struct FM2GP_mixed +{ + polynomial x; + polynomial y; + polynomial z; + + FM2GP_mixed() + { + boost::array const x_data = {{-2.2, -3.3, 0, 1}}; + boost::array const y_data = {{-4.4, 0, 1}}; + boost::array const z_data= {{-2, 1}}; + x = polynomial(x_data.begin(), x_data.end()); + y = polynomial(y_data.begin(), y_data.end()); + z = polynomial(z_data.begin(), z_data.end()); + } +}; + + +template struct FM2GP_trivial { polynomial x; @@ -233,8 +254,46 @@ struct FM2GP_trivial } }; +BOOST_AUTO_TEST_SUITE(Stein_gcd_ufd) + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_1_int, T, integral_test_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = Stein_gcd(fixture_type::x, fixture_type::y); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Stein_gcd(fixture_type::y, fixture_type::x); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_2_int, T, integral_test_types, FM2GP_Ex_8_3__2 ) +{ + typedef FM2GP_Ex_8_3__2 fixture_type; + polynomial w; + w = Stein_gcd(fixture_type::x, fixture_type::y); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Stein_gcd(fixture_type::y, fixture_type::x); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} -BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_1, T, pod_test_types, FM2GP_Ex_8_3__1 ) +BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_trivial, T, integral_test_types, FM2GP_trivial ) +{ + typedef FM2GP_trivial fixture_type; + polynomial w; + w = Stein_gcd(fixture_type::x, fixture_type::y); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Stein_gcd(fixture_type::y, fixture_type::x); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_AUTO_TEST_SUITE_END() + + +BOOST_AUTO_TEST_SUITE(Stein_gcd_field) + +BOOST_AUTO_TEST_CASE_EXPECTED_FAILURES(Ex_8_3__1, 2) + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, pod_floating_point_types, FM2GP_Ex_8_3__1 ) { typedef FM2GP_Ex_8_3__1 fixture_type; polynomial w; @@ -245,7 +304,9 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_1, T, pod_test_types, FM2GP_Ex_ } -BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_2, T, pod_test_types, FM2GP_Ex_8_3__2 ) +BOOST_AUTO_TEST_CASE_EXPECTED_FAILURES(Ex_8_3__2, 2) + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, pod_floating_point_types, FM2GP_Ex_8_3__2 ) { typedef FM2GP_Ex_8_3__2 fixture_type; polynomial w; @@ -256,7 +317,7 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_2, T, pod_test_types, FM2GP_Ex_ } -BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_trivial, T, pod_test_types, FM2GP_trivial ) +BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, pod_floating_point_types, FM2GP_trivial ) { typedef FM2GP_trivial fixture_type; polynomial w; @@ -266,41 +327,89 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_trivial, T, pod_test_types, FM2 BOOST_CHECK_EQUAL(w, fixture_type::z); } +BOOST_AUTO_TEST_SUITE_END() -BOOST_AUTO_TEST_CASE_TEMPLATE( test_Euclidean_gcd, T, floating_point_types) -{ - /* NOTE: Euclidean gcd is not yet customized to return THE greatest - * common polynomial divisor. If d is THE greatest common divisior of u and - * v, then gcd(u, v) will return d or -d according to the algorithm. - * By convention, it should return d, as for example Maxima and Wolfram - * Alpha do. - * This test is an example of the fact that it returns -d. - */ - boost::array const d8 = {{105, 278, -88, -56, 16}}; - boost::array const d6 = {{70, 232, -44, -64, 16}}; - boost::array const d2 = {{35, -24, 4}}; - polynomial u(d8.begin(), d8.end()); - polynomial v(d6.begin(), d6.end()); - polynomial expected(d2.begin(), d2.end()); - polynomial result; - result = Euclid_gcd(u, v); - BOOST_CHECK(result == expected || result == expected * T(-1.0)); - result = Euclid_gcd(v, u); - BOOST_CHECK(result == expected || result == expected * T(-1.0)); - - boost::array const s3 = {{-2, -3, 0, 1}}; - boost::array const s2 = {{-4, 0, 1}}; - boost::array const d_expected = {{-2, 1}}; - u = polynomial(s3.begin(), s3.end()); - v = polynomial(s2.begin(), s2.end()); - expected = polynomial(d_expected.begin(), d_expected.end()); - result = Euclid_gcd(u, v); - BOOST_CHECK(result == expected || result == expected * T(-1.0)); - result = Euclid_gcd(v, u); - BOOST_CHECK(result == expected || result == expected * T(-1.0)); + +BOOST_AUTO_TEST_SUITE(Euclid_gcd_field) + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, floating_point_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = gcd_traits< polynomial >::abs(Euclid_gcd(fixture_type::x, fixture_type::y)); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = gcd_traits< polynomial >::abs(Euclid_gcd(fixture_type::y, fixture_type::x)); + BOOST_CHECK_EQUAL(w, fixture_type::z); } +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, floating_point_types, FM2GP_Ex_8_3__2 ) +{ + typedef FM2GP_Ex_8_3__2 fixture_type; + polynomial w; + w = gcd_traits< polynomial >::abs(Euclid_gcd(fixture_type::x, fixture_type::y)); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = gcd_traits< polynomial >::abs(Euclid_gcd(fixture_type::y, fixture_type::x)); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, floating_point_types, FM2GP_trivial ) +{ + typedef FM2GP_trivial fixture_type; + polynomial w; + w = gcd_traits< polynomial >::abs(Euclid_gcd(fixture_type::x, fixture_type::y)); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = gcd_traits< polynomial >::abs(Euclid_gcd(fixture_type::y, fixture_type::x)); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_AUTO_TEST_SUITE_END() + + +BOOST_AUTO_TEST_SUITE(Euclid_gcd_ufd) + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, integral_test_types, FM2GP_Ex_8_3__1 ) +{ + typedef FM2GP_Ex_8_3__1 fixture_type; + polynomial w; + w = Euclid_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Euclid_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, integral_test_types, FM2GP_Ex_8_3__2 ) +{ + typedef FM2GP_Ex_8_3__2 fixture_type; + polynomial w; + w = Euclid_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Euclid_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + + +BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, integral_test_types, FM2GP_trivial ) +{ + typedef FM2GP_trivial fixture_type; + polynomial w; + w = Euclid_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); + w = Euclid_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); + BOOST_CHECK_EQUAL(w, fixture_type::z); +} + +BOOST_AUTO_TEST_SUITE_END() + + BOOST_AUTO_TEST_CASE_TEMPLATE( test_addition, T, all_test_types ) { polynomial const a(d3a.begin(), d3a.end()); From 92bc6777b61e353082ce2342b399f92aa2200873 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 5 Dec 2016 19:34:37 +1100 Subject: [PATCH 17/34] Use Stepanov or Joux implementation of subtraction as appropriate. --- include/boost/math/tools/polynomial_gcd.hpp | 49 ++++++++++++++++++----------- 1 file changed, 30 insertions(+), 19 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index 52644dd48..13ddca0bf 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -120,29 +120,40 @@ struct gcd_traits_polynomial::i inline static void subtract(polynomial_type &a, polynomial_type const &b) { -#if 1 - // Stepanov's implementation; suffers from floating point inaccuracy. - using std::floor; + using std::modf; + // We want to use Stepanov's implementation as often as possible because + // it results in the smallest coefficients, however, whole numbers take + // precedence. T const r = constant_coefficient(a) / constant_coefficient(b); - a -= r * b; - // normalize coefficients so that leading coefficient is whole - if (a && floor(a.data().back()) != a.data().back()) - a /= a.data().back(); -#else - using boost::lambda::_1; - - // Antoine Joux's implementation: huge coefficients. - T const tmp = constant_coefficient(a); - a *= constant_coefficient(b); - a -= tmp * b; - if (std::numeric_limits::has_infinity) + T r_int; + T const r_frac = modf(r, &r_int); + // Should we also consider 0.25, 0.125, etc? + if (r_frac == T(0) || r_frac == T(0.5)) { - T const &inf(std::numeric_limits::infinity()); - if (std::find_if(a.data().begin(), a.data().end(), _1 == inf || _1 == -inf) != a.data().end()) - throw std::domain_error("Floating point overflow."); - } + // Stepanov's implementation; suffers from floating point inaccuracy + // when r does not divide ___ evenly. + a -= r * b; + // normalize coefficients so that leading coefficient is whole + if (a && modf(a.data().back(), &r_int) != T(0)) + a /= a.data().back(); + } + else + { + // Antoine Joux's implementation: produces huge coefficients. + T const tmp = constant_coefficient(a); + a *= constant_coefficient(b); + a -= tmp * b; +#ifndef NDEBUG + using boost::lambda::_1; + if (std::numeric_limits::has_infinity) + { + T const &inf(std::numeric_limits::infinity()); + if (std::find_if(a.data().begin(), a.data().end(), _1 == inf || _1 == -inf) != a.data().end()) + throw std::domain_error("Floating point overflow."); + } #endif + } } }; From 5ae445e6706567356295ff1d50d62b513ae6bc6e Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 5 Dec 2016 21:57:32 +1100 Subject: [PATCH 18/34] Common gcd traits for polynomials. --- include/boost/math/tools/polynomial_gcd.hpp | 98 +++++++++++++---------------- 1 file changed, 42 insertions(+), 56 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index 13ddca0bf..d451a4292 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -16,45 +16,18 @@ namespace boost { namespace math { -// -// Special handling for polynomials: -// Note that gcd_traits_polynomial is templated on the coefficient type T, -// not on polynomial. -// -template -struct gcd_traits_polynomial -{}; - - -// gcd_traits for Z[x]. +// Common gcd traits for polynomials. template -struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > +struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math::tools::polynomial > { typedef boost::math::tools::polynomial polynomial_type; - // Just do a no-op for abs() as we'll normalize() later. static polynomial_type abs(const polynomial_type &val) { return leading_coefficient(val) < T(0) ? -val : val; } - inline static - void normalize(polynomial_type &x) - { - using boost::lambda::_1; - - // This does assume that the coefficients are totally ordered. - if (x) - { - if (leading_coefficient(x) < T(0)) - x.negate(); - // Find the first non-zero, because we can't do gcd(0, 0). - T const d = gcd_range(find_if(x.data().begin(), x.data().end(), _1 != T(0)), x.data().end()).first; - x /= d; - } - } - inline static int make_odd(boost::math::tools::polynomial &x) { unsigned r = 0; @@ -71,6 +44,43 @@ struct gcd_traits_polynomial::is { return a.size() < b.size(); } + + inline static + void normalize(polynomial_type &x) + { + using boost::lambda::_1; + + // This does assume that the coefficients are totally ordered. + if (x) + { + if (leading_coefficient(x) < T(0)) + x.negate(); + // Find the first non-zero, because we can't do gcd(0, 0). + T const d = gcd_range(find_if(x.data().begin(), x.data().end(), _1 != T(0)), x.data().end()).first; + x /= d; + } + } +}; + + + + +// +// Special handling for polynomials: +// Note that gcd_traits_polynomial is templated on the coefficient type T, +// not on polynomial. +// +template +struct gcd_traits_polynomial +{}; + + +// gcd_traits for Z[x]. +template +struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_polynomial_defaults +{ + typedef boost::math::tools::polynomial polynomial_type; + using gcd_traits_polynomial_defaults::normalize; inline static void subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) @@ -86,37 +96,13 @@ struct gcd_traits_polynomial::is // TODO: Enable for multiprecision floating point types. template -struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_defaults< boost::math::tools::polynomial > +struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_polynomial_defaults { typedef boost::math::tools::polynomial polynomial_type; typedef gcd_traits_defaults< boost::math::tools::polynomial > parent_class; using typename parent_class::method_type; - static const method_type method = method_type::method_euclid; - - static polynomial_type - abs(const polynomial_type &val) - { - return leading_coefficient(val) < T(0) ? -val : val; - } - - - inline static int make_odd(polynomial_type &x) - { - unsigned r = 0; - while (even(x)) - { - x >>= 1; - r++; - } - return r; - } - - inline static bool - less(polynomial_type const &a, polynomial_type const &b) - { - return a.size() < b.size(); - } - + // static const method_type method = method_type::method_euclid; + inline static void subtract(polynomial_type &a, polynomial_type const &b) { From 0629a9ff542b2970c4c2d35470ac3c00fbf02080 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 5 Dec 2016 23:42:49 +1100 Subject: [PATCH 19/34] A kludge to allow gcd, which makes Stein_gcd work. --- include/boost/math/common_factor_rt.hpp | 22 +++++++++++++++++++++- 1 file changed, 21 insertions(+), 1 deletion(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index d21f93d6c..2912a162c 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -70,6 +70,11 @@ namespace boost { a -= b; } + inline static void modulo(T& a, const T& b) + { + a %= b; + } + enum method_type { method_euclid = 0, @@ -244,6 +249,21 @@ namespace boost { BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val) { unsigned result = gcd_traits::find_lsb(val); val >>= result; return result; } }; #endif + + // This is a kludge to prove a point and should be replaced by something + // that covers all (PoD) floating point types. + template <> + struct gcd_traits : public gcd_traits_defaults + { + typedef double T; + + inline static void modulo(T& a, const T& b) + { + using std::fmod; + + a = fmod(a, b); + } + }; template inline @@ -356,7 +376,7 @@ namespace detail using std::swap; while (b != EuclideanDomain(0)) { - a %= b; + gcd_traits::modulo(a, b); swap(a, b); } return a; From 252474156e2f1a64f91b070ff955c0f00c96111a Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 5 Dec 2016 23:43:29 +1100 Subject: [PATCH 20/34] If gcd_traits::normalize works, then the tests pass. --- test/test_polynomial.cpp | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index 1f7565342..dc6a486ab 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -309,28 +309,28 @@ BOOST_AUTO_TEST_SUITE_END() BOOST_AUTO_TEST_SUITE(Stein_gcd_field) -BOOST_AUTO_TEST_CASE_EXPECTED_FAILURES(Ex_8_3__1, 2) - BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, pod_floating_point_types, FM2GP_Ex_8_3__1 ) { typedef FM2GP_Ex_8_3__1 fixture_type; polynomial w; w = Stein_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); w = Stein_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); } -BOOST_AUTO_TEST_CASE_EXPECTED_FAILURES(Ex_8_3__2, 2) - BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, pod_floating_point_types, FM2GP_Ex_8_3__2 ) { typedef FM2GP_Ex_8_3__2 fixture_type; polynomial w; w = Stein_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); w = Stein_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); } @@ -340,8 +340,10 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( trivial, T, pod_floating_point_types, FM2GP_tr typedef FM2GP_trivial fixture_type; polynomial w; w = Stein_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); w = Stein_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); } From fa0fac2fe449c41d73ff7a5a5ef4445b48a87735 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 5 Dec 2016 23:55:09 +1100 Subject: [PATCH 21/34] Remove gcd header from polynomial. --- include/boost/math/tools/polynomial.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/boost/math/tools/polynomial.hpp b/include/boost/math/tools/polynomial.hpp index ee013cded..73210fb48 100644 --- a/include/boost/math/tools/polynomial.hpp +++ b/include/boost/math/tools/polynomial.hpp @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include From 4df9cf560340b804dabacb7092d48bc12f4e7efe Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 5 Dec 2016 23:57:56 +1100 Subject: [PATCH 22/34] Use the convenience typedef polynomial_type. --- include/boost/math/tools/polynomial_gcd.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index d451a4292..d9afa7344 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -28,7 +28,7 @@ struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math: return leading_coefficient(val) < T(0) ? -val : val; } - inline static int make_odd(boost::math::tools::polynomial &x) + inline static int make_odd(polynomial_type &x) { unsigned r = 0; while (even(x)) @@ -40,7 +40,7 @@ struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math: } inline static bool - less(boost::math::tools::polynomial const &a, boost::math::tools::polynomial const &b) + less(polynomial_type const &a, polynomial_type const &b) { return a.size() < b.size(); } @@ -83,7 +83,7 @@ struct gcd_traits_polynomial::is using gcd_traits_polynomial_defaults::normalize; inline static void - subtract(boost::math::tools::polynomial &a, boost::math::tools::polynomial const &b) + subtract(polynomial_type &a, polynomial_type const &b) { // Antoine Joux's implementation: huge coefficients. T const tmp = constant_coefficient(a); @@ -99,7 +99,7 @@ template struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_polynomial_defaults { typedef boost::math::tools::polynomial polynomial_type; - typedef gcd_traits_defaults< boost::math::tools::polynomial > parent_class; + typedef gcd_traits_defaults< polynomial_type > parent_class; using typename parent_class::method_type; // static const method_type method = method_type::method_euclid; From cf0ff3d05bf16d73bdf4500f4ef92d4ae02bc3b9 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 09:43:42 +1100 Subject: [PATCH 23/34] Include normalize in abs_defaults; add modulo defaults; always normalize. --- include/boost/math/common_factor_rt.hpp | 55 +++++++++++++++++++-------------- test/test_polynomial.cpp | 6 ++++ 2 files changed, 38 insertions(+), 23 deletions(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 2912a162c..646aee180 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -36,7 +36,10 @@ namespace boost { struct gcd_traits_abs_defaults { inline static const T& abs(const T& val) { return val; } + + inline static void normalize(T&) {} }; + template struct gcd_traits_abs_defaults { @@ -45,10 +48,32 @@ namespace boost { using std::abs; return abs(val); } + + inline static void normalize(T& val) { val = abs(val); } + }; + + + template ::value || (std::numeric_limits::is_specialized && !std::numeric_limits::is_integer)> + struct gcd_traits_modulo_defaults + { + inline static void modulo(T& a, const T& b) + { + using std::fmod; + a = fmod(a, b); + } + }; + + template + struct gcd_traits_modulo_defaults + { + inline static void modulo(T& a, const T& b) + { + a %= b; + } }; template - struct gcd_traits_defaults : public gcd_traits_abs_defaults + struct gcd_traits_defaults : public gcd_traits_abs_defaults, gcd_traits_modulo_defaults { BOOST_FORCEINLINE static unsigned make_odd(T& val) { @@ -60,6 +85,7 @@ namespace boost { } return r; } + inline static bool less(const T& a, const T& b) { return a < b; @@ -70,11 +96,6 @@ namespace boost { a -= b; } - inline static void modulo(T& a, const T& b) - { - a %= b; - } - enum method_type { method_euclid = 0, @@ -88,12 +109,13 @@ namespace boost { boost::has_right_shift_assign::value && boost::has_left_shift_assign::value && boost::has_less::value ? method_binary : method_euclid; }; + // // Default gcd_traits just inherits from defaults: // template struct gcd_traits : public gcd_traits_defaults {}; - // + // Some platforms have fast bitscan operations, that allow us to implement // make_odd much more efficiently: // @@ -249,21 +271,6 @@ namespace boost { BOOST_FORCEINLINE static unsigned make_odd(wchar_t& val) { unsigned result = gcd_traits::find_lsb(val); val >>= result; return result; } }; #endif - - // This is a kludge to prove a point and should be replaced by something - // that covers all (PoD) floating point types. - template <> - struct gcd_traits : public gcd_traits_defaults - { - typedef double T; - - inline static void modulo(T& a, const T& b) - { - using std::fmod; - - a = fmod(a, b); - } - }; template inline @@ -421,7 +428,9 @@ namespace detail template inline Integer gcd(Integer const &a, Integer const &b) { - return detail::optimal_gcd_select(static_cast(gcd_traits::abs(a)), static_cast(gcd_traits::abs(b))); + Integer result = detail::optimal_gcd_select(static_cast(gcd_traits::abs(a)), static_cast(gcd_traits::abs(b))); + gcd_traits::normalize(result); + return result; } template diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index dc6a486ab..a3c9c7d7f 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -279,8 +279,10 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_1_int, T, integral_test_types, typedef FM2GP_Ex_8_3__1 fixture_type; polynomial w; w = Stein_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); w = Stein_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); } @@ -289,8 +291,10 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_2_int, T, integral_test_types, typedef FM2GP_Ex_8_3__2 fixture_type; polynomial w; w = Stein_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); w = Stein_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); } @@ -299,8 +303,10 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( test_Stein_gcd_trivial, T, integral_test_types typedef FM2GP_trivial fixture_type; polynomial w; w = Stein_gcd(fixture_type::x, fixture_type::y); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); w = Stein_gcd(fixture_type::y, fixture_type::x); + gcd_traits< polynomial >::normalize(w); BOOST_CHECK_EQUAL(w, fixture_type::z); } From 010de67b3669a6e2a31fd22fd7653f1f54cf91d8 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 09:45:26 +1100 Subject: [PATCH 24/34] Use the same approach for Z[x] as for R[X]. --- include/boost/math/tools/polynomial_gcd.hpp | 25 ++++++++++++++++++++----- 1 file changed, 20 insertions(+), 5 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index d9afa7344..867d893c0 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -85,11 +85,26 @@ struct gcd_traits_polynomial::is inline static void subtract(polynomial_type &a, polynomial_type const &b) { - // Antoine Joux's implementation: huge coefficients. - T const tmp = constant_coefficient(a); - a *= constant_coefficient(b); - a -= tmp * b; - normalize(a); + // We want to use Stepanov's implementation as often as possible because + // it results in the smallest coefficients, however, whole numbers take + // precedence. + T m = constant_coefficient(a); + gcd_traits::modulo(m, constant_coefficient(b)); + if (m == 0) + { + T const r = constant_coefficient(a) / constant_coefficient(b); + // Stepanov's implementation; suffers from floating point inaccuracy + // when r does not divide ___ evenly. + a -= r * b; + } + else + { + // Antoine Joux's implementation: produces huge coefficients. + T const tmp = constant_coefficient(a); + a *= constant_coefficient(b); + a -= tmp * b; + normalize(a); + } } }; From 2f141cf756ec1e5107dc30865898868ee5fba7f3 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 10:18:55 +1100 Subject: [PATCH 25/34] Consolidate the polynomial normalization gcd_trait. --- include/boost/math/tools/polynomial_gcd.hpp | 92 +++++++---------------------- 1 file changed, 22 insertions(+), 70 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index 867d893c0..32a7c9fb4 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -60,28 +60,7 @@ struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math: x /= d; } } -}; - - - - -// -// Special handling for polynomials: -// Note that gcd_traits_polynomial is templated on the coefficient type T, -// not on polynomial. -// -template -struct gcd_traits_polynomial -{}; - -// gcd_traits for Z[x]. -template -struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_polynomial_defaults -{ - typedef boost::math::tools::polynomial polynomial_type; - using gcd_traits_polynomial_defaults::normalize; - inline static void subtract(polynomial_type &a, polynomial_type const &b) { @@ -109,61 +88,34 @@ struct gcd_traits_polynomial::is }; -// TODO: Enable for multiprecision floating point types. +// +// Special handling for polynomials: +// Note that gcd_traits_polynomial is templated on the coefficient type T, +// not on polynomial. +// +template +struct gcd_traits_polynomial +{}; + + +// Note: Use these trait classes for Z[x] and R[x], etc to customize for operations +// on different kinds of polynomials. + +// gcd_traits for Z[x]. template -struct gcd_traits_polynomial::is_integer >::type > : public gcd_traits_polynomial_defaults -{ - typedef boost::math::tools::polynomial polynomial_type; - typedef gcd_traits_defaults< polynomial_type > parent_class; - using typename parent_class::method_type; - // static const method_type method = method_type::method_euclid; +struct gcd_traits_polynomial::is_integer>::type> : public gcd_traits_polynomial_defaults +{}; - inline static void - subtract(polynomial_type &a, polynomial_type const &b) - { - using std::modf; - - // We want to use Stepanov's implementation as often as possible because - // it results in the smallest coefficients, however, whole numbers take - // precedence. - T const r = constant_coefficient(a) / constant_coefficient(b); - T r_int; - T const r_frac = modf(r, &r_int); - // Should we also consider 0.25, 0.125, etc? - if (r_frac == T(0) || r_frac == T(0.5)) - { - // Stepanov's implementation; suffers from floating point inaccuracy - // when r does not divide ___ evenly. - a -= r * b; - // normalize coefficients so that leading coefficient is whole - if (a && modf(a.data().back(), &r_int) != T(0)) - a /= a.data().back(); - } - else - { - // Antoine Joux's implementation: produces huge coefficients. - T const tmp = constant_coefficient(a); - a *= constant_coefficient(b); - a -= tmp * b; -#ifndef NDEBUG - using boost::lambda::_1; - if (std::numeric_limits::has_infinity) - { - T const &inf(std::numeric_limits::infinity()); - if (std::find_if(a.data().begin(), a.data().end(), _1 == inf || _1 == -inf) != a.data().end()) - throw std::domain_error("Floating point overflow."); - } -#endif - } - } -}; + +// gcd_traits for R[x]. +template +struct gcd_traits_polynomial::value || (std::numeric_limits::is_specialized && !std::numeric_limits::is_integer) >::type> : public gcd_traits_polynomial_defaults +{}; template struct gcd_traits< boost::math::tools::polynomial > : public gcd_traits_polynomial -{ - -}; +{}; } // namespace math } // namespace boost From 32f2e83ae4849fad5ab91d407abca6ca079f4ac4 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 10:37:58 +1100 Subject: [PATCH 26/34] To implicitly cast to bool or not, that is the question. The answer is yes, for now. --- include/boost/math/tools/polynomial_gcd.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index 32a7c9fb4..893253921 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -69,7 +69,7 @@ struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math: // precedence. T m = constant_coefficient(a); gcd_traits::modulo(m, constant_coefficient(b)); - if (m == 0) + if (!m) { T const r = constant_coefficient(a) / constant_coefficient(b); // Stepanov's implementation; suffers from floating point inaccuracy From 4a0bfd8e86b7e87620fea68e704ff591e4be40e9 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 11:36:16 +1100 Subject: [PATCH 27/34] TODO: Decorate the failing Euclid test. --- test/test_polynomial.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/test/test_polynomial.cpp b/test/test_polynomial.cpp index a3c9c7d7f..7029e2353 100644 --- a/test/test_polynomial.cpp +++ b/test/test_polynomial.cpp @@ -327,6 +327,8 @@ BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__1, T, pod_floating_point_types, FM2GP_ BOOST_CHECK_EQUAL(w, fixture_type::z); } +// TODO +// BOOST_AUTO_TEST_CASE_EXPECTED_FAILURES(Ex_8_3__2, 6) BOOST_FIXTURE_TEST_CASE_TEMPLATE( Ex_8_3__2, T, pod_floating_point_types, FM2GP_Ex_8_3__2 ) { From b5d54da915d9f852d17fe959fd298e1eb7bdae03 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 13:11:09 +1100 Subject: [PATCH 28/34] Explicitly convert return value to bool. --- include/boost/math/common_factor_rt.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index 646aee180..ec3f3f08d 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -277,7 +277,7 @@ inline typename enable_if_c::is_integer, bool>::type odd(const T& x) { - return x & 1u; + return bool(x & 1u); } template From 2982b3b5cc8c71ec66bb9319aad552fbfbbb3bac Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Thu, 8 Dec 2016 13:17:45 +1100 Subject: [PATCH 29/34] Use library's mixed_binary_gcd, forbid gcd(0, 0). --- reporting/performance/test_gcd.cpp | 50 ++++++-------------------------------- 1 file changed, 8 insertions(+), 42 deletions(-) diff --git a/reporting/performance/test_gcd.cpp b/reporting/performance/test_gcd.cpp index 8e365db99..baf0924f0 100644 --- a/reporting/performance/test_gcd.cpp +++ b/reporting/performance/test_gcd.cpp @@ -169,47 +169,7 @@ T binary_textbook(T u, T v) return u + v; } -// -// The Mixed Binary Euclid Algorithm -// Sidi Mohamed Sedjelmaci -// Electronic Notes in Discrete Mathematics 35 (2009) 169–176 -// -template -T mixed_binary_gcd(T u, T v) -{ - using std::swap; - if(u < v) - swap(u, v); - - unsigned shifts = 0; - - if(!u) - return v; - if(!v) - return u; - - while(even(u) && even(v)) - { - u >>= 1u; - v >>= 1u; - ++shifts; - } - while(v > 1) - { - u %= v; - v -= u; - if(!u) - return v << shifts; - if(!v) - return u << shifts; - while(even(u)) u >>= 1u; - while(even(v)) v >>= 1u; - if(u < v) - swap(u, v); - } - return (v == 1 ? v : u) << shifts; -} template void test_type(const char* name) @@ -220,7 +180,12 @@ void test_type(const char* name) for(unsigned i = 0; i < 1000; ++i) { - data.push_back(std::make_pair(get_prime_products(), get_prime_products())); + data.push_back(pair()); + do + { + data.back() = std::make_pair(get_prime_products(), get_prime_products()); + } + while (data.back().first == 0 && data.back().second == 0); } std::string row_name("gcd<"); row_name += name; @@ -342,7 +307,8 @@ inline typename enable_if_c Date: Mon, 12 Dec 2016 01:11:36 +1100 Subject: [PATCH 30/34] Simplify Euclidean algorithm very slightly using (x) vs (x != 0). --- include/boost/math/common_factor_rt.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/boost/math/common_factor_rt.hpp b/include/boost/math/common_factor_rt.hpp index ec3f3f08d..e3a5c1c2a 100644 --- a/include/boost/math/common_factor_rt.hpp +++ b/include/boost/math/common_factor_rt.hpp @@ -381,7 +381,7 @@ namespace detail inline EuclideanDomain Euclid_gcd(EuclideanDomain a, EuclideanDomain b) { using std::swap; - while (b != EuclideanDomain(0)) + while (b) { gcd_traits::modulo(a, b); swap(a, b); From f4f94f18357c15e9e3562cc5c1927e68b57c8239 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Mon, 12 Dec 2016 01:23:08 +1100 Subject: [PATCH 31/34] Instead of normalizing after Joux's subtraction, find gcd(a(0), b(0))first. --- include/boost/math/tools/polynomial_gcd.hpp | 18 ++++++++---------- 1 file changed, 8 insertions(+), 10 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index 893253921..dd1547851 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -64,25 +64,23 @@ struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math: inline static void subtract(polynomial_type &a, polynomial_type const &b) { - // We want to use Stepanov's implementation as often as possible because - // it results in the smallest coefficients, however, whole numbers take - // precedence. + // We use Stepanov's implementation when the constant coefficients + // divide evenly, and Joux's otherwise. T m = constant_coefficient(a); gcd_traits::modulo(m, constant_coefficient(b)); if (!m) { T const r = constant_coefficient(a) / constant_coefficient(b); - // Stepanov's implementation; suffers from floating point inaccuracy - // when r does not divide ___ evenly. a -= r * b; } else { - // Antoine Joux's implementation: produces huge coefficients. - T const tmp = constant_coefficient(a); - a *= constant_coefficient(b); - a -= tmp * b; - normalize(a); + // Antoine Joux's implementation tempered by coefficient gcd. + T const a0 = constant_coefficient(a); + T const b0 = constant_coefficient(b); + T const gcd_a0b0 = gcd(a0, b0); + a *= b0 / gcd_a0b0; + a -= (a0 / gcd_a0b0) * b; } } }; From cd4f1255e492de4ebd918e86318410978c589a99 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Wed, 14 Dec 2016 19:37:41 +1100 Subject: [PATCH 32/34] Add performance reporting test for gcd of polynomials. --- reporting/performance/test_polynomial_gcd.cpp | 239 ++++++++++++++++++++++++++ 1 file changed, 239 insertions(+) create mode 100644 reporting/performance/test_polynomial_gcd.cpp diff --git a/reporting/performance/test_polynomial_gcd.cpp b/reporting/performance/test_polynomial_gcd.cpp new file mode 100644 index 000000000..3f5817871 --- /dev/null +++ b/reporting/performance/test_polynomial_gcd.cpp @@ -0,0 +1,239 @@ +// Copyright Jeremy Murphy 2016. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifdef _MSC_VER +# pragma warning (disable : 4224) +#endif + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include "fibonacci.hpp" +#include "../../test/table_type.hpp" +#include "table_helper.hpp" +#include "performance.hpp" + + +using namespace std; +using namespace boost::math::tools; + +polynomial total_sum(0); + + +template +double exec_timed_test_foo(Func f, const Table& data, double min_elapsed = 0.5) +{ + double t = 0; + unsigned repeats = 1; + typename Table::value_type::first_type sum{0}; + stopwatch w; + do + { + for(unsigned count = 0; count < repeats; ++count) + { + for(typename Table::size_type n = 0; n < data.size(); ++n) + sum += f(data[n].first, data[n].second); + } + + t = boost::chrono::duration_cast>(w.elapsed()).count(); + if(t < min_elapsed) + repeats *= 2; + } + while(t < min_elapsed); + total_sum += sum; + return t / repeats; +} + + +template +struct test_function_template +{ + vector > const & data; + const char* data_name; + + test_function_template(vector > const &data, const char* name) : data(data), data_name(name) {} + + template + void operator()(pair const &f) const + { + auto result = exec_timed_test_foo(f.first, data); + auto table_name = string("gcd method comparison with ") + compiler_name() + string(" on ") + platform_name(); + + report_execution_time(result, + table_name, + string(data_name), + string(f.second) + "\n" + boost_name()); + } +}; + +boost::random::mt19937 rng; +boost::random::uniform_int_distribution<> d_0_6(0, 6); +boost::random::uniform_int_distribution<> d_1_20(1, 20); + +template +T get_prime_products() +{ + int n_primes = d_0_6(rng); + switch(n_primes) + { + case 0: + // Generate a power of 2: + // return static_cast(1u) << d_1_20(rng); + return std::pow(2, d_1_20(rng)); + case 1: + // prime number: + return boost::math::prime(d_1_20(rng) + 3); + } + T result = 1; + for(int i = 0; i < n_primes; ++i) + result *= boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3) * boost::math::prime(d_1_20(rng) + 3); + return result; +} + +template +T get_uniform_random() +{ + static boost::random::uniform_int_distribution minmax(std::numeric_limits::min(), std::numeric_limits::max()); + return minmax(rng); +} + +template +inline bool even(T const& val) +{ + return !(val & 1u); +} + + +template +typename boost::enable_if_c::value, void>::type +gcd_algorithms(OutputIterator output) +{ + *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial >, "Stein_gcd"); +} + + +template +typename boost::enable_if_c::value, void>::type +gcd_algorithms(OutputIterator output) +{ + *output++ = std::make_pair(boost::math::detail::Euclid_gcd< polynomial >, "Euclid_gcd"); + *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial >, "Stein_gcd"); +} + + +template +void test_type(const char* name) +{ + using namespace boost::math::detail; + std::vector, polynomial > > data; + + std::cout << "Testing: " << typeid(T).name() << std::endl; + + for(unsigned i = 0; i < 10; ++i) + { + data.push_back(pair< polynomial, polynomial >()); + for (unsigned j = 0; j != 10; j++) + { + data.back().first.data().push_back(get_prime_products()); + data.back().second.data().push_back(get_prime_products()); + } + } + std::string row_name("gcd<"); + row_name += name; + row_name += "> (random prime number products)"; + + typedef pair< function(polynomial, polynomial)>, string> f_test; + vector test_functions; + gcd_algorithms(back_inserter(test_functions)); + for_each(begin(test_functions), end(test_functions), test_function_template< polynomial >(data, row_name.c_str())); + + data.clear(); +#if 0 + for(unsigned i = 0; i < 1000; ++i) + { + data.push_back(std::make_pair(get_uniform_random(), get_uniform_random())); + } + row_name.erase(); + row_name += "gcd<"; + row_name += name; + row_name += "> (uniform random numbers)"; + for_each(begin(test_functions), end(test_functions), test_function_template(data, row_name.c_str())); + + // Fibonacci number tests: + row_name.erase(); + row_name += "gcd<"; + row_name += name; + row_name += "> (adjacent Fibonacci numbers)"; + for_each(begin(test_functions), end(test_functions), test_function_template(fibonacci_numbers_permution_1(), row_name.c_str())); + + row_name.erase(); + row_name += "gcd<"; + row_name += name; + row_name += "> (permutations of Fibonacci numbers)"; + for_each(begin(test_functions), end(test_functions), test_function_template(fibonacci_numbers_permution_2(), row_name.c_str())); + + row_name.erase(); + row_name += "gcd<"; + row_name += name; + row_name += "> (Trivial cases)"; + for_each(begin(test_functions), end(test_functions), test_function_template(trivial_gcd_test_cases(), row_name.c_str())); +#endif +} + +/*******************************************************************************************************************/ + +template +T generate_random(unsigned bits_wanted) +{ + static boost::random::mt19937 gen; + typedef boost::random::mt19937::result_type random_type; + + T max_val; + unsigned digits; + if(std::numeric_limits::is_bounded && (bits_wanted == (unsigned)std::numeric_limits::digits)) + { + max_val = (std::numeric_limits::max)(); + digits = std::numeric_limits::digits; + } + else + { + max_val = T(1) << bits_wanted; + digits = bits_wanted; + } + + unsigned bits_per_r_val = std::numeric_limits::digits - 1; + while((random_type(1) << bits_per_r_val) > (gen.max)()) --bits_per_r_val; + + unsigned terms_needed = digits / bits_per_r_val + 1; + + T val = 0; + for(unsigned i = 0; i < terms_needed; ++i) + { + val *= (gen.max)(); + val += gen(); + } + val %= max_val; + return val; +} + + + +int main() +{ + test_type< int >("polynomial"); +} From 1cd5679cbc411ba8b4f660637a62f2acf6cbe239 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 22 Jan 2017 18:50:05 +1100 Subject: [PATCH 33/34] Modularize the various subtraction methods for Stein gcd. --- include/boost/math/tools/polynomial_gcd.hpp | 68 +++++++++++++++++++++-------- 1 file changed, 50 insertions(+), 18 deletions(-) diff --git a/include/boost/math/tools/polynomial_gcd.hpp b/include/boost/math/tools/polynomial_gcd.hpp index dd1547851..7ae602c47 100644 --- a/include/boost/math/tools/polynomial_gcd.hpp +++ b/include/boost/math/tools/polynomial_gcd.hpp @@ -61,27 +61,59 @@ struct gcd_traits_polynomial_defaults : public gcd_traits_defaults< boost::math: } } + /** + * Alexander Stepanov's subtraction normalizes b with the ratio of the + * constant coefficients before subtracting it from a. + * + * It does not appear to be reliable due to the floating point arithmetic + * inaccuracy in the division operation. + */ + inline static void + Stepanov_subtraction(polynomial_type &a, polynomial_type const &b) + { + using std::modf; + + T r = constant_coefficient(a) / constant_coefficient(b); + a -= r * b; + if (a && modf(a.data().back(), &r)) + a /= a.data().back(); + } + + /** + * Joux subtraction multiplies each polynomial by the other's constant + * coefficient then subtracts one from the other. + * + * Joux, Antoine. Algorithmic cryptanalysis. CRC Press, 2009. + */ + inline static void + Joux_subtraction(polynomial_type &a, polynomial_type const &b) + { + T const a0 = constant_coefficient(a); + a *= constant_coefficient(b); + a -= a0 * b; + } + + + /** + * Factored Joux subtraction first calculates the gcd of the constant + * coefficients and removes that factor before multiplying the polynomials. + * The effect is to reduce intermediate value swell. + */ + inline static void + factored_Joux_subtraction(polynomial_type &a, polynomial_type const &b) + { + T const a0 = constant_coefficient(a); + T const b0 = constant_coefficient(b); + T const gcd_a0b0 = gcd(a0, b0); + a *= b0 / gcd_a0b0; + a -= a0 / gcd_a0b0 * b; + } + + inline static void subtract(polynomial_type &a, polynomial_type const &b) { - // We use Stepanov's implementation when the constant coefficients - // divide evenly, and Joux's otherwise. - T m = constant_coefficient(a); - gcd_traits::modulo(m, constant_coefficient(b)); - if (!m) - { - T const r = constant_coefficient(a) / constant_coefficient(b); - a -= r * b; - } - else - { - // Antoine Joux's implementation tempered by coefficient gcd. - T const a0 = constant_coefficient(a); - T const b0 = constant_coefficient(b); - T const gcd_a0b0 = gcd(a0, b0); - a *= b0 / gcd_a0b0; - a -= (a0 / gcd_a0b0) * b; - } + gcd_traits< boost::math::tools::polynomial >::Joux_subtraction(a, b); } }; From 082cb40c1859ac3cc5baf951f799cbbb1fd12848 Mon Sep 17 00:00:00 2001 From: "Jeremy W. Murphy" Date: Sun, 22 Jan 2017 18:50:55 +1100 Subject: [PATCH 34/34] Update performance test to support different Stein implementations. --- reporting/performance/test_polynomial_gcd.cpp | 168 +++++++++++++++++++++++--- 1 file changed, 151 insertions(+), 17 deletions(-) diff --git a/reporting/performance/test_polynomial_gcd.cpp b/reporting/performance/test_polynomial_gcd.cpp index 3f5817871..7e50dd9f6 100644 --- a/reporting/performance/test_polynomial_gcd.cpp +++ b/reporting/performance/test_polynomial_gcd.cpp @@ -1,4 +1,4 @@ -// Copyright Jeremy Murphy 2016. +// Copyright Jeremy Murphy 2017. // Use, modification and distribution are subject to the // Boost Software License, Version 1.0. (See accompanying file // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) @@ -11,6 +11,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -23,6 +26,8 @@ #include #include #include +#include + #include "fibonacci.hpp" #include "../../test/table_type.hpp" #include "table_helper.hpp" @@ -31,8 +36,94 @@ using namespace std; using namespace boost::math::tools; +using boost::math::gcd_traits; polynomial total_sum(0); +polynomial total_sum_int(0); +typedef boost::multiprecision::cpp_bin_float_quad mp_float_type; +polynomial total_sum_mp(0); + + +template +typename boost::enable_if_c::value, void>::type +add_total(polynomial const &x) +{ + total_sum += x; +} + + +template +typename boost::enable_if_c::value, void>::type +add_total(polynomial const &x) +{ + total_sum += x; +} + +void add_total(polynomial const &x) +{ + total_sum_mp += x; +} + +template +SteinDomain Stein_gcd_Joux(SteinDomain m, SteinDomain n) +{ + using std::swap; + BOOST_ASSERT(m || n); + if (!m) + return n; + if (!n) + return m; + unsigned shifts = std::min(gcd_traits::make_odd(m), gcd_traits::make_odd(n)); + // odd(m) && odd(n) + while (m != n) + { + if (gcd_traits::less(m, n)) + swap(n, m); + gcd_traits::Joux_subtraction(m, n); + BOOST_ASSERT(even(m)); + // With polynomials for example, it is possible that m is now zero. + if (!m) + { + n <<= shifts; + return n; + } + gcd_traits::make_odd(m); + } + // m == n + m <<= shifts; + return m; +} + + +template +SteinDomain Stein_gcd_factored_Joux(SteinDomain m, SteinDomain n) +{ + using std::swap; + BOOST_ASSERT(m || n); + if (!m) + return n; + if (!n) + return m; + unsigned shifts = std::min(gcd_traits::make_odd(m), gcd_traits::make_odd(n)); + // odd(m) && odd(n) + while (m != n) + { + if (gcd_traits::less(m, n)) + swap(n, m); + gcd_traits::factored_Joux_subtraction(m, n); + BOOST_ASSERT(even(m)); + // With polynomials for example, it is possible that m is now zero. + if (!m) + { + n <<= shifts; + return n; + } + gcd_traits::make_odd(m); + } + // m == n + m <<= shifts; + return m; +} template @@ -55,7 +146,7 @@ double exec_timed_test_foo(Func f, const Table& data, double min_elapsed = 0.5) repeats *= 2; } while(t < min_elapsed); - total_sum += sum; + add_total(sum); return t / repeats; } @@ -66,13 +157,17 @@ struct test_function_template vector > const & data; const char* data_name; - test_function_template(vector > const &data, const char* name) : data(data), data_name(name) {} + test_function_template(vector > const &data, const char* name) : data(data), data_name(name) + { + cerr << "Testing: " << name << endl; + } template void operator()(pair const &f) const { + cerr << "Algorithm: " << f.second << endl; auto result = exec_timed_test_foo(f.first, data); - auto table_name = string("gcd method comparison with ") + compiler_name() + string(" on ") + platform_name(); + auto table_name = string("polynomial gcd method comparison with ") + compiler_name() + string(" on ") + platform_name(); report_execution_time(result, table_name, @@ -105,13 +200,39 @@ T get_prime_products() return result; } + +// T is integral. +template +typename boost::enable_if_c::value, T>::type +get_uniform_random() +{ + static boost::random::uniform_int_distribution minmax(0, 255); + return minmax(rng); +} + + +// T is floating point, limit to small values so that IVS does not crash test. +template +typename boost::enable_if_c::value, T>::type +get_uniform_random() +{ + using std::round; + static boost::random::uniform_real_distribution minmax(0, 3); + return round(minmax(rng)); +} + + +// T is not a POD. template -T get_uniform_random() +typename boost::enable_if_c::value, T>::type +get_uniform_random() { - static boost::random::uniform_int_distribution minmax(std::numeric_limits::min(), std::numeric_limits::max()); - return minmax(rng); + using std::round; + static boost::random::uniform_real_distribution minmax(0, 3); + return round(minmax(rng)); } + template inline bool even(T const& val) { @@ -126,32 +247,42 @@ gcd_algorithms(OutputIterator output) *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial >, "Stein_gcd"); } - template -typename boost::enable_if_c::value, void>::type +typename boost::enable_if_c::value, void>::type gcd_algorithms(OutputIterator output) { - *output++ = std::make_pair(boost::math::detail::Euclid_gcd< polynomial >, "Euclid_gcd"); *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial >, "Stein_gcd"); + // *output++ = std::make_pair(boost::math::detail::subresultant_gcd< polynomial >, "subresultant gcd"); +} + + +template +typename boost::enable_if_c::value || !std::numeric_limits::is_integer, void>::type +gcd_algorithms(OutputIterator output) +{ + // *output++ = std::make_pair(boost::math::detail::Euclid_gcd< polynomial >, "Euclid_gcd"); + *output++ = std::make_pair(boost::math::detail::Stein_gcd< polynomial >, "Stein_gcd (Stepanov-factored Joux)"); + *output++ = std::make_pair(Stein_gcd_factored_Joux< polynomial >, "Stein_gcd (factored Joux)"); + *output++ = std::make_pair(Stein_gcd_Joux< polynomial >, "Stein_gcd (Joux)"); } template -void test_type(const char* name) +void test_type(const std::string name) { using namespace boost::math::detail; std::vector, polynomial > > data; - std::cout << "Testing: " << typeid(T).name() << std::endl; - for(unsigned i = 0; i < 10; ++i) { data.push_back(pair< polynomial, polynomial >()); - for (unsigned j = 0; j != 10; j++) + for (unsigned j = 0; j != 5; j++) { - data.back().first.data().push_back(get_prime_products()); - data.back().second.data().push_back(get_prime_products()); + data.back().first.data().push_back(get_uniform_random()); + data.back().second.data().push_back(get_uniform_random()); } + data.back().first.normalize(); + data.back().second.normalize(); } std::string row_name("gcd<"); row_name += name; @@ -235,5 +366,8 @@ T generate_random(unsigned bits_wanted) int main() { - test_type< int >("polynomial"); + // test_type("polynomial"); + test_type("polynomial"); + test_type("polynomial"); + test_type("polynomial"); }