boostorg/multiprecision
1
0
Fork 0
mirror of https://github.com/boostorg/multiprecision.git synced 2026-07-21 13:23:49 +00:00
Code Issues Packages Projects Releases Wiki Activity

Improved Karatsuba using aliasing.

Browse Source
This commit is contained in:
jzmaddock
2019-12-23 18:40:16 +00:00
parent 90abd9583d
commit 4e46a7a010
16 changed files with 196 additions and 70 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/boost/multiprecision/cpp_int/multiply.hpp
+38 -47
View File
@@ -70,7 +70,7 @@ inline BOOST_MP_CXX14_CONSTEXPR void resize_for_carry(cpp_int_backend<MinBits1,
if (result.size() < required)
result.resize(required, required);
}
const size_t karatsuba_cutoff = 100;
const size_t karatsuba_cutoff = 50;
template <unsigned MinBits1, unsigned MaxBits1, cpp_integer_type SignType1, cpp_int_check_type Checked1, class Allocator1, unsigned MinBits2, unsigned MaxBits2, cpp_integer_type SignType2, cpp_int_check_type Checked2, class Allocator2, unsigned MinBits3, unsigned MaxBits3, cpp_integer_type SignType3, cpp_int_check_type Checked3, class Allocator3>
inline BOOST_MP_CXX14_CONSTEXPR typename enable_if_c<!is_trivial_cpp_int<cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits2, MaxBits2, SignType2, Checked2, Allocator2> >::value && !is_trivial_cpp_int<cpp_int_backend<MinBits3, MaxBits3, SignType3, Checked3, Allocator3> >::value>::type
eval_multiply(
@@ -134,67 +134,52 @@ eval_multiply(
constexpr const unsigned limb_bits = sizeof(limb_type) * CHAR_BIT;
#endif
result.resize(as + bs, as + bs - 1);
typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pr = result.limbs();
BOOST_STATIC_ASSERT(double_limb_max - 2 * limb_max >= limb_max * limb_max);
std::fill(pr, pr + result.size(), 0);
double_limb_type carry = 0;
#if 1
if(as >= karatsuba_cutoff && bs >= karatsuba_cutoff)
{
unsigned n = (as > bs ? as : bs) / 2 + 1;
// write a, b as a = a_h * 2^n + a_l, b = b_h * 2^n + b_l
cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a_h, a_l, b_h, b_l;
unsigned sz = (std::min)(as, n);
a_l.resize(sz, sz);
std::copy(a.limbs(), a.limbs() + sz, a_l.limbs());
const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a_l(a.limbs(), 0, sz);
sz = (std::min)(bs, n);
b_l.resize(sz, sz);
std::copy(b.limbs(), b.limbs() + sz, b_l.limbs());
if(as > n){
a_h.resize(as - n, as - n);
std::copy(a.limbs() + n, a.limbs() + as, a_h.limbs());
}
else
{
a_h.resize(1,1);
a_h.limbs()[0] = 0;
}
if(bs > n)
{
b_h.resize(bs - n, bs - n);
std::copy(b.limbs() + n, b.limbs() + bs, b_h.limbs());
}
else
{
b_h.resize(1,1);
b_h.limbs()[0] = 0;
}
const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> b_l(b.limbs(), 0, sz);
limb_type zero = 0;
const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> a_h(as > n ? a.limbs() + n : &zero, 0, as > n ? as - n : 1);
const cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> b_h(bs > n ? b.limbs() + n : &zero, 0, bs > n ? bs - n : 1);
// x = a_h * b_ h
// y = a_l * b_l
// z = (a_h + a_l)*(b_h + b_l) - x - y
// a * b = x * (2 ^ (2 * n))+ z * (2 ^ n) + y
cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1, t2;
const typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::scoped_shared_storage storage(result, 9 * n);
cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> t1(storage, 0, 2 * n), t2(storage, 2 * n, 3 * n), t3(storage, 5 * n, 2 * n), t4(storage, 7 * n, 2 * n);
// result = | a_h*b_h | a_l*b_l |
// (bits) <-- 2*n -->
eval_multiply(t1, a_l, b_l);
eval_multiply(t2, a_h, b_h);
resize_for_carry(result, t1.size());
std::copy(t1.limbs(), t1.limbs() + t1.size(), result.limbs());
resize_for_carry(result, 2 * n + t2.size());
sz = (std::max<long long int>)(0, static_cast<long long int>((std::min)(result.size(), 2 * n + t2.size())) - static_cast<long long int>(2 * n));
// for fixed precision arithmetic a_h and b_h cannot be non-zero
// otherwise a_h * b_h cannot be placed in result.size() bits
std::copy(t2.limbs(), t2.limbs() + sz, result.limbs() + 2 * n);
eval_add(t1, t2); // t1 = a_l*b_l + a_h*b_h
eval_add(a_l, a_h);
eval_add(b_l, b_h);
eval_multiply(t2, a_l, b_l); // t2 = (a_h+a_l)*(b_h+b_l)
if constexpr (!boost::is_void<Allocator1>::value)
{
cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> result_low(result.limbs(), 0, 2 * n);
cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1> result_high(result.limbs(), 2 * n, result.size() - 2 * n);
eval_multiply(result_low, a_l, b_l);
for (unsigned i = result_low.size(); i < 2 * n; ++i)
result.limbs()[i] = 0;
eval_multiply(result_high, a_h, b_h);
for (unsigned i = result_high.size() + 2 * n; i < result.size(); ++i)
result.limbs()[i] = 0;
eval_add(t1, result_low, result_high); // t1 = a_l*b_l + a_h*b_h
}
else
{
eval_multiply(t1, a_l, b_l);
eval_multiply(t2, a_h, b_h);
std::copy(t1.limbs(), t1.limbs() + t1.size(), result.limbs());
std::copy(t2.limbs(), t2.limbs() + t2.size(), result.limbs() + 2 * n);
eval_add(t1, t2); // t1 = a_l*b_l + a_h*b_h
}
eval_add(t3, a_l, a_h);
eval_add(t4, b_l, b_h);
eval_multiply(t2, t3, t4); // t2 = (a_h+a_l)*(b_h+b_l)
eval_subtract(t2, t1);
eval_left_shift(t2, n * limb_bits); // t2 = 2^n * ( (a_h+a_l)*(b_h*b_l) - a_h*b_h - a_l*b_l)
eval_add(result, t2);
@@ -203,6 +188,12 @@ eval_multiply(
result.sign(a.sign() != b.sign());
return ;
}
#endif
typename cpp_int_backend<MinBits1, MaxBits1, SignType1, Checked1, Allocator1>::limb_pointer pr = result.limbs();
BOOST_STATIC_ASSERT(double_limb_max - 2 * limb_max >= limb_max * limb_max);
std::fill(pr, pr + result.size(), 0);
double_limb_type carry = 0;
for (unsigned i = 0; i < as; ++i)
{
performance/Jamfile.v2
+17 -15
View File
@@ -50,21 +50,23 @@ else
test-suite performance :
[ exe performance_test : performance_test.cpp
performance_test_files/test01.cpp performance_test_files/test16.cpp performance_test_files/test31.cpp
performance_test_files/test02.cpp performance_test_files/test17.cpp performance_test_files/test32.cpp
performance_test_files/test03.cpp performance_test_files/test18.cpp performance_test_files/test33.cpp
performance_test_files/test04.cpp performance_test_files/test19.cpp performance_test_files/test34.cpp
performance_test_files/test05.cpp performance_test_files/test20.cpp performance_test_files/test35.cpp
performance_test_files/test06.cpp performance_test_files/test21.cpp performance_test_files/test36.cpp
performance_test_files/test07.cpp performance_test_files/test22.cpp performance_test_files/test37.cpp
performance_test_files/test08.cpp performance_test_files/test23.cpp performance_test_files/test38.cpp
performance_test_files/test09.cpp performance_test_files/test24.cpp performance_test_files/test39.cpp
performance_test_files/test10.cpp performance_test_files/test25.cpp performance_test_files/test40.cpp
performance_test_files/test11.cpp performance_test_files/test26.cpp performance_test_files/test41.cpp
performance_test_files/test12.cpp performance_test_files/test27.cpp performance_test_files/test42.cpp
performance_test_files/test13.cpp performance_test_files/test28.cpp performance_test_files/test43.cpp
performance_test_files/test14.cpp performance_test_files/test29.cpp performance_test_files/test44.cpp
performance_test_files/test15.cpp performance_test_files/test30.cpp
performance_test_files/test01.cpp performance_test_files/test18.cpp performance_test_files/test35.cpp
performance_test_files/test02.cpp performance_test_files/test19.cpp performance_test_files/test36.cpp
performance_test_files/test03.cpp performance_test_files/test20.cpp performance_test_files/test37.cpp
performance_test_files/test04.cpp performance_test_files/test21.cpp performance_test_files/test38.cpp
performance_test_files/test05.cpp performance_test_files/test22.cpp performance_test_files/test39.cpp
performance_test_files/test06.cpp performance_test_files/test23.cpp performance_test_files/test40.cpp
performance_test_files/test07.cpp performance_test_files/test24.cpp performance_test_files/test41.cpp
performance_test_files/test08.cpp performance_test_files/test25.cpp performance_test_files/test42.cpp
performance_test_files/test09.cpp performance_test_files/test26.cpp performance_test_files/test43.cpp
performance_test_files/test10.cpp performance_test_files/test27.cpp performance_test_files/test44.cpp
performance_test_files/test11.cpp performance_test_files/test28.cpp performance_test_files/test45.cpp
performance_test_files/test12.cpp performance_test_files/test29.cpp performance_test_files/test46.cpp
performance_test_files/test13.cpp performance_test_files/test30.cpp performance_test_files/test47.cpp
performance_test_files/test14.cpp performance_test_files/test31.cpp performance_test_files/test48.cpp
performance_test_files/test15.cpp performance_test_files/test32.cpp performance_test_files/test49.cpp
performance_test_files/test16.cpp performance_test_files/test33.cpp performance_test_files/test50.cpp
performance_test_files/test17.cpp performance_test_files/test34.cpp performance_test_files/test51.cpp
/boost/system//boost_system
: release
[ check-target-builds ../config//has_gmp : <define>TEST_MPF <define>TEST_MPZ <define>TEST_MPQ <source>gmp : ]
performance/performance_test.cpp
+5 -2
View File
@@ -99,9 +99,12 @@ void quickbook_results()
}
}
#ifdef BOOST_HAS_UNISTD_H
#if defined(__HAS_INCLUDE)
#if __has_include(<sys/utsname.h>)
#define HAS_UTSNAME
#include <sys/utsname.h>
#endif
#endif
#ifdef _WIN32
#include <windows.h>
#endif
@@ -109,7 +112,7 @@ void quickbook_results()
void quickbook_platform_details()
{
std::cout << "[table:platform Platform Details\n[[Platform][";
#ifdef BOOST_HAS_UNISTD_H
#ifdef HAS_UTSNAME
utsname name;
uname(&name);
std::cout << name.sysname << " " << name.release << ", version " << name.version << ", " << name.machine << "]]\n";
performance/performance_test_files/test39.cpp
+1 -1
View File
@@ -11,6 +11,6 @@
void test39()
{
#ifdef TEST_MPZ
test<boost::multiprecision::mpz_int>("gmp_int", 1024*5);
test<boost::multiprecision::mpz_int>("gmp_int", 1024*2);
#endif
}
performance/performance_test_files/test40.cpp
+1 -1
View File
@@ -11,6 +11,6 @@
void test40()
{
#ifdef TEST_MPZ
test<boost::multiprecision::mpz_int>("gmp_int", 1024*10);
test<boost::multiprecision::mpz_int>("gmp_int", 1024*4);
#endif
}
performance/performance_test_files/test41.cpp
+1 -1
View File
@@ -11,6 +11,6 @@
void test41()
{
#ifdef TEST_CPP_INT
test<boost::multiprecision::cpp_int>("cpp_int", 1024 * 5);
test<boost::multiprecision::cpp_int>("cpp_int", (boost::multiprecision::backends::karatsuba_cutoff + 2) * sizeof(boost::multiprecision::limb_type) * CHAR_BIT);
#endif
}
performance/performance_test_files/test42.cpp
+1 -1
View File
@@ -11,6 +11,6 @@
void test42()
{
#ifdef TEST_CPP_INT
test<boost::multiprecision::cpp_int>("cpp_int", 1024 * 10);
test<boost::multiprecision::cpp_int>("cpp_int", 1024 * 2);
#endif
}
performance/performance_test_files/test43.cpp
+1 -1
View File
@@ -11,7 +11,7 @@
void test43()
{
#ifdef TEST_TOMMATH
test<boost::multiprecision::tom_int>("tommath_int", 1024*5);
test<boost::multiprecision::tom_int>("tommath_int", 1024*2);
/*
//
// These are actually too slow to test!!!
performance/performance_test_files/test44.cpp
+1 -1
View File
@@ -11,7 +11,7 @@
void test44()
{
#ifdef TEST_TOMMATH
test<boost::multiprecision::tom_int>("tommath_int", 1024*10);
test<boost::multiprecision::tom_int>("tommath_int", 1024*4);
/*
//
// These are actually too slow to test!!!
performance/performance_test_files/test45.cpp
+16
View File
@@ -0,0 +1,16 @@
///////////////////////////////////////////////////////////////
// Copyright 2019 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_CPP_INT)
#include <boost/multiprecision/cpp_int.hpp>
#endif
void test42()
{
#ifdef TEST_CPP_INT
test<boost::multiprecision::cpp_int>("cpp_int", 1024 * 4);
#endif
}
performance/performance_test_files/test46.cpp
+16
View File
@@ -0,0 +1,16 @@
///////////////////////////////////////////////////////////////
// Copyright 2019 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_CPP_INT)
#include <boost/multiprecision/cpp_int.hpp>
#endif
void test42()
{
#ifdef TEST_CPP_INT
test<boost::multiprecision::cpp_int>("cpp_int", 1024 * 8);
#endif
}
performance/performance_test_files/test47.cpp
+16
View File
@@ -0,0 +1,16 @@
///////////////////////////////////////////////////////////////
// Copyright 2019 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_CPP_INT)
#include <boost/multiprecision/cpp_int.hpp>
#endif
void test42()
{
#ifdef TEST_CPP_INT
test<boost::multiprecision::cpp_int>("cpp_int", 1024 * 16);
#endif
}
performance/performance_test_files/test48.cpp
+25
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 https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_TOMMATH)
#include <boost/multiprecision/tommath.hpp>
#endif
void test44()
{
#ifdef TEST_TOMMATH
test<boost::multiprecision::tom_int>("tommath_int", 1024*8);
/*
//
// These are actually too slow to test!!!
//
test<boost::multiprecision::tom_rational>("tom_rational", 128);
test<boost::multiprecision::tom_rational>("tom_rational", 256);
test<boost::multiprecision::tom_rational>("tom_rational", 512);
test<boost::multiprecision::tom_rational>("tom_rational", 1024);
*/
#endif
}
performance/performance_test_files/test49.cpp
+25
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 https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_TOMMATH)
#include <boost/multiprecision/tommath.hpp>
#endif
void test44()
{
#ifdef TEST_TOMMATH
test<boost::multiprecision::tom_int>("tommath_int", 1024*16);
/*
//
// These are actually too slow to test!!!
//
test<boost::multiprecision::tom_rational>("tom_rational", 128);
test<boost::multiprecision::tom_rational>("tom_rational", 256);
test<boost::multiprecision::tom_rational>("tom_rational", 512);
test<boost::multiprecision::tom_rational>("tom_rational", 1024);
*/
#endif
}
performance/performance_test_files/test50.cpp
+16
View File
@@ -0,0 +1,16 @@
///////////////////////////////////////////////////////////////
// Copyright 2019 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_MPZ)
#include <boost/multiprecision/gmp.hpp>
#endif
void test40()
{
#ifdef TEST_MPZ
test<boost::multiprecision::mpz_int>("gmp_int", 1024*8);
#endif
}
performance/performance_test_files/test51.cpp
+16
View File
@@ -0,0 +1,16 @@
///////////////////////////////////////////////////////////////
// Copyright 2019 John Maddock. Distributed under the Boost
// Software License, Version 1.0. (See accompanying file
// LICENSE_1_0.txt or copy at https://www.boost.org/LICENSE_1_0.txt
#include "../performance_test.hpp"
#if defined(TEST_MPZ)
#include <boost/multiprecision/gmp.hpp>
#endif
void test40()
{
#ifdef TEST_MPZ
test<boost::multiprecision::mpz_int>("gmp_int", 1024*16);
#endif
}