We need double_limb type for quotient calculation in Lehmers routine.

Fixes https://github.com/boostorg/multiprecision/issues/370.
This commit is contained in:
jzmaddock
2021-09-26 11:25:34 +01:00
parent 7a04dcdfb6
commit cab664d086
3 changed files with 30 additions and 6 deletions

View File

@@ -732,12 +732,10 @@ void eval_gcd_lehmer(cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Al
// but that division is very slow.
//
// We begin with a specialized routine for division.
// We know that u > 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_backend<MinBits1, MaxBits1, SignType1, Checked1, Al
//
while (true)
{
limb_type q = 1;
double_limb_type q = 1;
double_limb_type tt = u;
divide_subtract(q, u, v);
std::swap(u, v);

View File

@@ -1112,6 +1112,7 @@ test-suite misc :
[ run bug12039.cpp no_eh_support ]
[ compile git_issue_30.cpp ]
[ run git_issue_167.cpp ]
[ run git_issue_370.cpp ]
[ run git_issue_175.cpp ]
[ run git_issue_248.cpp ]
[ run git_issue_265.cpp : : : [ check-target-builds ../config//has_mpfr : <source>gmp <source>mpfr : <build>no ] ]

25
test/git_issue_370.cpp Normal file
View File

@@ -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 <boost/multiprecision/cpp_int.hpp>
#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;
}