From cab664d086ded9d356d2828325f6d62b3aecebdd Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 26 Sep 2021 11:25:34 +0100 Subject: [PATCH] We need double_limb type for quotient calculation in Lehmers routine. Fixes https://github.com/boostorg/multiprecision/issues/370. --- include/boost/multiprecision/cpp_int/misc.hpp | 10 +++----- test/Jamfile.v2 | 1 + test/git_issue_370.cpp | 25 +++++++++++++++++++ 3 files changed, 30 insertions(+), 6 deletions(-) create mode 100644 test/git_issue_370.cpp diff --git a/include/boost/multiprecision/cpp_int/misc.hpp b/include/boost/multiprecision/cpp_int/misc.hpp index cd9f29d2..2963bede 100644 --- a/include/boost/multiprecision/cpp_int/misc.hpp +++ b/include/boost/multiprecision/cpp_int/misc.hpp @@ -732,12 +732,10 @@ void eval_gcd_lehmer(cpp_int_backend v > ~limb_type(0), and therefore -// that the result will fit into a single limb_type. -// We also know that most of the time this is called the result will be 1. +// We know that most of the time this is called the result will be 1. // For small limb counts, this almost doubles the performance of Lehmer's routine! // -BOOST_FORCEINLINE void divide_subtract(limb_type& q, double_limb_type& u, const double_limb_type& v) +BOOST_FORCEINLINE void divide_subtract(double_limb_type& q, double_limb_type& u, const double_limb_type& v) { BOOST_ASSERT(q == 1); // precondition on entry. u -= v; @@ -746,7 +744,7 @@ BOOST_FORCEINLINE void divide_subtract(limb_type& q, double_limb_type& u, const u -= v; if (++q > 30) { - limb_type t = u / v; + double_limb_type t = u / v; u -= t * v; q += t; } @@ -855,7 +853,7 @@ void eval_gcd_lehmer(cpp_int_backendgmp mpfr : no ] ] diff --git a/test/git_issue_370.cpp b/test/git_issue_370.cpp new file mode 100644 index 00000000..ae56a87b --- /dev/null +++ b/test/git_issue_370.cpp @@ -0,0 +1,25 @@ +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2019 John Maddock. Distributed under 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) + +#include +#include "test.hpp" + +int main() +{ + const boost::multiprecision::cpp_int a("500000000000000052504760255204421627079393309355027816932345132815919505535709229444276879024105562954502314530690391078574434507015318513443905076213688875017942541908041275407131568575177172639474548726709751235383681696449966404295647940685784470144122251803020020951078103818191513659921807053133698549053838430992170843235673537548059987539601671975279280846041564435631581262016246808786828637048154067265620710396778995313534536353760281048487250661054626168637371167135426013683337484254647996964562455566714879467026196409913165805735073230830136024016362543811769017875638974011487488573436"); + const boost::multiprecision::cpp_int b("1500000000000000157514280765613264881238179928065083450797035398447758516607127688332830637072316688863506943592071173235723303521045955540331715228641066625053827625724123826221394705725531517918423646180129253706151045089349899212886943822057353410432366755409060062853234311454574540979765421159386595647161515292215193506006556519037965168192736708179557957863203557666055574947146355487693991882510747766220045897624670399027877365714431356466054500731862264092476764347207739651025585146903094168986610767496468412336047796468657032646893153521091155634158263410282629846280069312485301157888001"); + const boost::multiprecision::cpp_int r = gcd(a, b); + + if (a % r) + return 1; + if (b % r) + return 1; + if (gcd(a / r, b / r) != 1) + return 1; + + return 0; +} + +