Program Listing for File BigInteger.h
↰ Return to documentation for file (ADS/numeric/BigInteger.h)
/*--------------------------------------------------------------------------*/
/*--------------------------- File BigInteger.h ----------------------------*/
/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
/*----------------------------- DEFINITIONS --------------------------------*/
/*--------------------------------------------------------------------------*/
#ifndef __ADS_BigInteger
#define __ADS_BigInteger
/* self-identification: #endif constant */
/*--------------------------------------------------------------------------*/
/*------------------------------ INCLUDES ----------------------------------*/
/*--------------------------------------------------------------------------*/
#include "ADS/Contract.h"
#include <algorithm>
#include <array>
#include <bit>
#include <cmath>
#include <compare>
#include <cstddef>
#include <cstdint>
#include <initializer_list>
#include <memory>
#include <string>
#include <utility>
#include <vector>
// the width of the product of two digits done by the compiler: 128 with
// unsigned __int128, 64 with the intrinsics of MSVC on x64, 0 for none
#ifndef ADS_BIGINT_WIDE
#if defined( ADS_BIGINT_PORTABLE )
#define ADS_BIGINT_WIDE 0
#elif defined( __SIZEOF_INT128__ )
#define ADS_BIGINT_WIDE 128
#elif defined( _MSC_VER ) && defined( _M_X64 )
#define ADS_BIGINT_WIDE 64
#else
#define ADS_BIGINT_WIDE 0
#endif
#endif
#if ADS_BIGINT_WIDE == 64
#include <intrin.h>
#endif
// the lanes of the NTT: 256 with AVX2, else 0, one at a time
#ifndef ADS_BIGINT_SIMD
#if defined( __AVX2__ )
#define ADS_BIGINT_SIMD 256
#else
#define ADS_BIGINT_SIMD 0
#endif
#endif
#if ADS_BIGINT_SIMD > 0
#include <immintrin.h>
#endif
// the NTT on doubles, with the fused product and sum, where the lanes of
// AVX2 and the products of 128 bits are there
#ifndef ADS_BIGINT_FNTT
#if ( ADS_BIGINT_SIMD >= 256 ) && defined( __FMA__ ) && \
( ADS_BIGINT_WIDE == 128 )
#define ADS_BIGINT_FNTT 1
#else
#define ADS_BIGINT_FNTT 0
#endif
#endif
#if ADS_BIGINT_FNTT
#include <cmath>
#endif
// the loops of the sums and of the products by one digit in the assembly
// of x86-64, with MULX, where GCC or Clang compile for a processor that has
// BMI2 and ADX; ADS_BIGINT_NOASM keeps them in C++
#ifndef ADS_BIGINT_ASM
#if ( ADS_BIGINT_WIDE == 128 ) && defined( __x86_64__ ) && \
defined( __GNUC__ ) && defined( __BMI2__ ) && defined( __ADX__ ) && \
! defined( ADS_BIGINT_NOASM )
#define ADS_BIGINT_ASM 1
#else
#define ADS_BIGINT_ASM 0
#endif
#endif
// the backend of GMP, with ADS_WITH_GMP: the products, the quotients, the
// conversions and the greatest common divisor go to the mpn and mpz layers
// of GMP in the ranges of lengths where those are faster, on the digits of
// BigInteger as they are; only where a limb of GMP has 64 bits, all of
// them numbers, else BigInteger stays on its own code
#if defined( ADS_WITH_GMP )
#include <cstring>
#include <gmp.h>
#if ! defined( ADS_BIGINT_GMP ) && ( GMP_LIMB_BITS == 64 ) && \
( GMP_NAIL_BITS == 0 )
#define ADS_BIGINT_GMP 1
#endif
#endif
#ifndef ADS_BIGINT_GMP
#define ADS_BIGINT_GMP 0
#endif
/*--------------------------------------------------------------------------*/
/*--------------------------- NAMESPACE ------------------------------------*/
/*--------------------------------------------------------------------------*/
namespace ADS {
namespace detail {
template< class T >
struct NoInit : std::allocator< T > {
template< class U > struct rebind { using other = NoInit< U >; };
NoInit( void ) = default;
template< class U > NoInit( const NoInit< U > & ) noexcept {}
template< class U > void construct( U * p ) noexcept {
::new( static_cast< void * >( p ) ) U;
}
template< class U , class... A > void construct( U * p , A &&... a ) {
::new( static_cast< void * >( p ) ) U( std::forward< A >( a )... );
}
};
#if defined( __GNUC__ ) || defined( __clang__ )
#define ADS_BIGINT_NOINLINE __attribute__(( noinline ))
#elif defined( _MSC_VER )
#define ADS_BIGINT_NOINLINE __declspec( noinline )
#else
#define ADS_BIGINT_NOINLINE
#endif
struct DigitPool {
static constexpr unsigned kBins = 17 , kKeep = 8;
std::uint64_t * bin[ kBins ][ kKeep ];
unsigned cnt[ kBins ];
bool armed , dead;
};
inline constinit thread_local DigitPool digitPool{};
struct DigitReaper {
~DigitReaper( void ) {
for( unsigned k = 0 ; k < DigitPool::kBins ; ++k )
for( unsigned i = 0 ; i < digitPool.cnt[ k ] ; ++i )
std::allocator< std::uint64_t >().deallocate( digitPool.bin[ k ][ i ] ,
std::size_t( 1 ) << k );
for( unsigned & c : digitPool.cnt ) c = 0;
digitPool.dead = true; // the blocks freed later go back
}
};
inline thread_local DigitReaper digitReaper;
class Digits {
public:
using T = std::uint64_t;
Digits( void ) noexcept {}
explicit Digits( std::size_t n ) { resize( n ); }
Digits( std::size_t n , T v ) { resize( n , v ); }
Digits( const T * b , const T * e ) {
resize( std::size_t( e - b ) );
std::copy( b , e , data() );
}
Digits( std::initializer_list< T > l ) {
resize( l.size() );
std::copy( l.begin() , l.end() , data() );
}
Digits( const Digits & o ) : Digits( o.data() , o.data() + o.f_n ) {}
Digits( Digits && o ) noexcept { take( o ); }
Digits & operator=( const Digits & o ) {
if( this != &o ) {
f_n = 0;
resize( o.f_n );
std::copy( o.data() , o.data() + o.f_n , data() );
}
return( *this );
}
Digits & operator=( Digits && o ) noexcept {
if( this != &o ) {
free();
take( o );
}
return( *this );
}
~Digits( void ) { free(); }
[[nodiscard]] std::size_t size( void ) const { return( f_n ); }
[[nodiscard]] bool empty( void ) const { return( f_n == 0 ); }
T * data( void ) { return( f_p ); }
const T * data( void ) const { return( f_p ); }
T * begin( void ) { return( data() ); }
T * end( void ) { return( data() + f_n ); }
const T * begin( void ) const { return( data() ); }
const T * end( void ) const { return( data() + f_n ); }
T & operator[]( std::size_t i ) { return( data()[ i ] ); }
const T & operator[]( std::size_t i ) const { return( data()[ i ] ); }
T & back( void ) { return( data()[ f_n - 1 ] ); }
const T & back( void ) const { return( data()[ f_n - 1 ] ); }
void reserve( std::size_t n ) { grow( n , n ); }
void resize( std::size_t n ) {
grow( n , n );
f_n = n;
}
void resize( std::size_t n , T v ) {
grow( n , n );
if( n > f_n ) std::fill( data() + f_n , data() + n , v );
f_n = n;
}
void assign( std::size_t n , T v ) {
f_n = 0;
resize( n , v );
}
void push_back( T v ) {
grow( f_n + 1 , 2 * f_n );
data()[ f_n++ ] = v;
}
void pop_back( void ) { --f_n; }
friend bool operator==( const Digits & a , const Digits & b ) {
return( ( a.f_n == b.f_n ) &&
std::equal( a.data() , a.data() + a.f_n , b.data() ) );
}
private:
static constexpr std::size_t kIn = 2;
void grow( std::size_t n , std::size_t c ) {
if( n > f_cap ) regrow( n , c );
}
void free( void ) {
if( f_p != f_in ) {
release( f_p , f_cap );
f_p = f_in;
f_cap = kIn;
}
}
ADS_BIGINT_NOINLINE void regrow( std::size_t n , std::size_t c ) {
c = std::max( n , c );
T * q;
if( c < ( std::size_t( 1 ) << ( DigitPool::kBins - 1 ) ) ) {
const unsigned k = unsigned( std::bit_width( std::max< std::size_t >( c ,
4 ) - 1 ) );
c = std::size_t( 1 ) << k;
DigitPool & P = digitPool;
q = P.cnt[ k ] ? P.bin[ k ][ --P.cnt[ k ] ]
: std::allocator< T >().allocate( c );
}
else q = std::allocator< T >().allocate( c );
std::copy( data() , data() + f_n , q );
free();
f_p = q;
f_cap = c;
}
ADS_BIGINT_NOINLINE static void release( T * p , std::size_t cap ) {
DigitPool & P = digitPool;
const unsigned k = unsigned( std::countr_zero( cap ) );
if( std::has_single_bit( cap ) && ( k < DigitPool::kBins ) &&
( P.cnt[ k ] < ( k < 12 ? DigitPool::kKeep : 2u ) ) && ! P.dead ) {
if( ! P.armed ) { // the reaper, once per thread
(void) &digitReaper;
P.armed = true;
}
P.bin[ k ][ P.cnt[ k ]++ ] = p;
}
else std::allocator< T >().deallocate( p , cap );
}
void take( Digits & o ) {
f_n = o.f_n;
f_cap = o.f_cap;
if( o.f_p != o.f_in ) f_p = o.f_p;
else std::copy( o.f_in , o.f_in + f_n , f_in );
o.f_p = o.f_in;
o.f_n = 0;
o.f_cap = kIn;
}
T * f_p = f_in;
std::size_t f_n = 0 , f_cap = kIn;
T f_in[ kIn ] = {};
}; // end( class Digits )
} // end( namespace detail )
/*--------------------------------------------------------------------------*/
/*--------------------------- CLASS BigInteger -----------------------------*/
/*--------------------------------------------------------------------------*/
class BigInteger {
/*--------------------------------------------------------------------------*/
public:
using Limb = std::uint64_t;
// the thresholds, in digits, from the timings of bench/bignum.cpp on a
// Xeon Gold 5120, against GMP
static constexpr std::size_t kKaratsuba = 40;
static constexpr std::size_t kToom3 = 120;
static constexpr std::size_t kToom4 = 200;
static constexpr std::size_t kToom32 = 64;
static constexpr std::size_t kFFT = 768;
static constexpr std::size_t kSqrKaratsuba = 56;
static constexpr std::size_t kSqrToom3 = 160;
static constexpr std::size_t kSqrToom4 = 300;
static constexpr std::size_t kSqrFFT = 768;
static constexpr std::size_t kDivide = 96;
static constexpr std::size_t kConvert = 24;
static constexpr std::size_t kConvertF = 64;
static constexpr std::size_t kFracLeaf = 32;
static constexpr std::size_t kFracLevel = 4;
static constexpr std::size_t kFracCyclic = 448;
static constexpr std::size_t npos = std::size_t( -1 );
static constexpr std::size_t kModFold = 16;
static constexpr std::size_t kHensel = 20;
static constexpr std::size_t kNewton = 768;
static constexpr std::size_t kNewtonKept = 512;
static constexpr std::size_t kRecip = 32;
static constexpr std::size_t kDivCut = 900;
static constexpr std::size_t kDivShort = 16;
static constexpr std::size_t kDivBlock = 40;
static constexpr std::size_t kDivBlock6 = 64;
static constexpr std::size_t kDivExact2 = 32;
static constexpr std::size_t kHgcd = 128;
static constexpr std::size_t kGcdDC = 300;
// with ADS_WITH_GMP, the lengths in digits that go to GMP, where it was
// faster in the timings of the same benchmark: the products whose shorter
// factor has [ kGmpMul , kGmpMulTop ) digits, [ kGmpMul , kGmpMulLongTop )
// when the longer one is twice as long or more, the squares below
// kGmpSqrTop; the quotients with the remainder whose divisor is below
// kGmpDivTop or whose quotient is below kGmpShort, the quotients alone
// longer than the divisor as those, the others from kGmpQuot digits of
// divisor below kGmpQuotTop, or at most half as long as the divisor; the
// exact quotients below kGmpExactTop, the conversions to a base not a
// power of 2 from kConvert digits below kGmpToStringTop, from one at any
// length; the greatest common divisors whose shorter argument has
// kGmpGcd digits or more
static constexpr std::size_t kGmpMul = 5;
static constexpr std::size_t kGmpMulTop = 1536;
static constexpr std::size_t kGmpMulLongTop = 704;
static constexpr std::size_t kGmpSqrTop = 1920;
static constexpr std::size_t kGmpDivTop = 3584;
static constexpr std::size_t kGmpShort = 512;
static constexpr std::size_t kGmpQuot = 8;
static constexpr std::size_t kGmpQuotTop = 49152;
static constexpr std::size_t kGmpExactTop = 262144;
static constexpr std::size_t kGmpToStringTop = 49152;
static constexpr std::size_t kGmpGcd = 2048;
BigInteger( long long v = 0 );
explicit BigInteger( const std::string & s , int base = 10 );
[[nodiscard]] std::string toString( int base = 10 ) const;
[[nodiscard]] std::string toStringSchool( int base = 10 ) const;
static BigInteger fromStringSchool( const std::string & s , int base = 10 );
[[nodiscard]] int sign( void ) const;
[[nodiscard]] BigInteger abs( void ) const;
[[nodiscard]] std::size_t limbs( void ) const;
BigInteger operator-( void ) const;
friend BigInteger operator+( const BigInteger & a , const BigInteger & b );
friend BigInteger operator-( const BigInteger & a , const BigInteger & b );
friend BigInteger operator*( const BigInteger & a , const BigInteger & b );
friend BigInteger operator/( const BigInteger & a , const BigInteger & b );
friend BigInteger operator%( const BigInteger & a , const BigInteger & b );
BigInteger & operator+=( const BigInteger & b );
BigInteger & operator-=( const BigInteger & b );
BigInteger & operator*=( const BigInteger & b );
friend bool operator==( const BigInteger & a , const BigInteger & b ) =
default;
friend std::strong_ordering operator<=>( const BigInteger & a ,
const BigInteger & b );
static std::pair< BigInteger , BigInteger > divMod( const BigInteger & a ,
const BigInteger & b );
static std::pair< BigInteger , BigInteger > divModSchool(
const BigInteger & a , const BigInteger & b );
static BigInteger multiplySchool( const BigInteger & a ,
const BigInteger & b );
static BigInteger multiplyKaratsuba( const BigInteger & a ,
const BigInteger & b );
static BigInteger multiplyFFT( const BigInteger & a , const BigInteger & b );
static BigInteger divExact( BigInteger a , const BigInteger & b );
static void addMul( BigInteger & r , const BigInteger & a ,
const BigInteger & b );
static BigInteger pow( const BigInteger & x , unsigned long n );
static std::pair< BigInteger , BigInteger > sqrtRem( const BigInteger & a );
static BigInteger sqrt( const BigInteger & a );
static BigInteger sqrtNewton( const BigInteger & a );
static BigInteger gcd( const BigInteger & a , const BigInteger & b );
static BigInteger gcdEuclid( BigInteger a , BigInteger b );
/*--------------------------------------------------------------------------*/
private:
using Mag = detail::Digits;
static constexpr std::size_t kNTTMax = std::size_t( 1 ) << 25;
BigInteger( bool neg , Mag d );
static void trim( Mag & d );
static int compareMag( const Mag & a , const Mag & b );
static int cmpN( const Limb * a , const Limb * b , std::size_t n );
static Limb addN( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb );
static Limb subN( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb );
static Limb addTo( Limb * r , std::size_t nr , const Limb * x ,
std::size_t nx );
static Limb subFrom( Limb * r , std::size_t nr , const Limb * x ,
std::size_t nx );
static Limb mul1( Limb * r , const Limb * a , std::size_t n , Limb m );
static Limb addMul1( Limb * r , const Limb * a , std::size_t n , Limb m );
static Limb subMul1( Limb * r , const Limb * a , std::size_t n , Limb m );
static void addMul2( Limb * r , const Limb * a , std::size_t n , Limb m0 ,
Limb m1 );
static bool absDiff( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb );
static Mag addMag( const Mag & a , const Mag & b );
static Mag subMag( const Mag & a , const Mag & b );
static BigInteger addSigned( const BigInteger & a , const BigInteger & b ,
bool bneg );
static void mulSchool( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb );
static void sqrSchool( Limb * r , const Limb * a , std::size_t n );
static void sqrTile( Limb * r , const Limb * a , std::size_t n );
static void mulN( Limb * r , const Limb * a , const Limb * b ,
std::size_t n , Limb * w , std::size_t fft ,
std::size_t toom );
static std::size_t mulScratch( std::size_t n );
static void karatsuba( Limb * r , const Limb * a , const Limb * b ,
std::size_t n , Limb * w , std::size_t fft ,
std::size_t toom );
static void toom3( Limb * r , const Limb * a , const Limb * b ,
std::size_t n , Limb * w , std::size_t fft ,
std::size_t toom );
static void toom4( Limb * r , const Limb * a , const Limb * b ,
std::size_t n , Limb * w , std::size_t fft ,
std::size_t toom );
static bool toom32( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , Limb * w ,
std::size_t fft , std::size_t toom );
static bool toom42( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , Limb * w ,
std::size_t fft , std::size_t toom );
static void divExact3( Limb * a , std::size_t n );
static void rshN( Limb * a , std::size_t n , int s );
static void divExactB1( Limb * a , std::size_t n , Limb f );
static Limb addLsh1( Limb * r , const Limb * a , const Limb * b ,
std::size_t n );
static Limb subLsh1( Limb * r , const Limb * a , const Limb * b ,
std::size_t n );
static void rsh1Add( Limb * r , const Limb * a , const Limb * b ,
std::size_t n );
static void rsh1Sub( Limb * r , const Limb * a , const Limb * b ,
std::size_t n );
static void mulNTT( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb );
template< class O >
static void mulNTT3( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , std::size_t n );
template< class O >
static void mulFNTT( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , std::size_t n ,
bool cyclic = false );
static void mulCyclic( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , std::size_t n );
static void mulP( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb ,
std::size_t fft = kFFT , std::size_t toom = kToom3 );
static Mag mulMag( const Mag & a , const Mag & b , std::size_t fft = kFFT ,
std::size_t toom = kToom3 );
static void mulRec( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , std::size_t fft ,
std::size_t toom );
static Limb divSmall( Mag & a , Limb d );
static Limb divSmall( Mag & a , Limb dn , Limb v , int s );
static Limb modSmall( const Limb * a , std::size_t n , Limb d );
static void divExactSmall( Mag & a , Limb d );
static void mulAddSmall( Mag & a , Limb m , Limb c );
static Limb divBase( Limb * a , std::size_t m , const Limb * b ,
std::size_t n , Limb * q , bool cut = false );
template< std::size_t K >
static bool divBlocks( Limb * a , std::size_t & m , const Limb * b ,
std::size_t n , Limb * q , bool cut , Limb d1 ,
Limb d0 , Limb v , std::size_t least );
static Limb divRec( Limb * a , std::size_t m , const Limb * b ,
std::size_t n , Limb * q , Limb * t );
static void recip( Limb * x , const Limb * a , std::size_t n );
static void divShort( const Mag & u , const Mag & v , Mag & q , Mag * r );
static Limb divNewton( Limb * a , std::size_t m , const Limb * b ,
std::size_t n , Limb * q , Mag * keep = nullptr ,
bool quot = false );
static std::pair< Mag , Mag > divModMag( const Mag & u , const Mag & v ,
bool school = false ,
Mag * keep = nullptr ,
bool rem = true );
struct HMat {
explicit HMat( std::size_t len );
Mag d;
Limb * p[ 4 ];
std::size_t n , alloc;
};
struct HMat1 { Limb u[ 4 ]; };
static bool hgcd2( Limb ah , Limb al , Limb bh , Limb bl , HMat1 & M );
static std::size_t hgcdInv1( const HMat1 & M , Limb * a , Limb * b ,
std::size_t n );
static void hgcdMul1( HMat & M , const HMat1 & M1 );
static void hgcdMul( HMat & M , const HMat & M1 );
static constexpr std::size_t kHgcdMul7 = 32;
static void hgcdMul7( HMat & M , const HMat & M1 );
static void hgcdUpdate( HMat & M , const Limb * q , std::size_t qn ,
unsigned c );
static std::size_t hgcdAdjust( const HMat & M , std::size_t n , Limb * a ,
Limb * b , std::size_t p );
static std::size_t gcdSubdiv( Limb * a , Limb * b , std::size_t n ,
std::size_t s , HMat * M , Mag * g );
static std::size_t hgcdStep( std::size_t n , Limb * a , Limb * b ,
std::size_t s , HMat & M );
static std::size_t hgcd( Limb * a , Limb * b , std::size_t n , HMat & M );
static Mag gcd22( Limb u1 , Limb u0 , Limb v1 , Limb v0 );
static Mag gcdMag( const Mag & x , const Mag & y );
static void figures( Mag a , std::size_t pad , int base ,
std::string & out , Mag * b = nullptr ,
std::size_t padb = 0 );
static void toFigures( const Mag & x , std::size_t pad , int base ,
const std::vector< Mag > & pw ,
const std::vector< std::size_t > & zs , int k ,
std::vector< Mag > & rc , std::vector< Mag > & rf ,
std::string & out , std::size_t j );
static void fracFigures( const Mag & x , std::size_t j , int base , int k ,
Limb P , const std::vector< Mag > & pw ,
const std::vector< std::size_t > & zs ,
std::vector< Mag > & rf , std::string & out );
static void fracTree( Mag & f , std::size_t j , int base , int k , Limb P ,
const std::vector< Mag > & pw ,
const std::vector< std::size_t > & zs , bool up ,
char * out );
static void fracLeaf( Mag & f , std::size_t c , Limb P , int k , int base ,
bool up , char * out );
static bool mulBand( Limb * r , const Limb * a , std::size_t na ,
const Limb * b , std::size_t nb , std::size_t s ,
std::size_t e );
static Mag fromChunks( const Limb * c , std::size_t n ,
const std::vector< Mag > & pw , Limb P );
static void sqrt2( Limb h , Limb l , Mag & s , Mag & r );
static void sqrtRec( const Limb * m , std::size_t k , Mag & s , Mag & r ,
bool root = false );
static std::pair< Mag , Mag > sqrtMag( const Mag & a , bool rem );
static constexpr std::size_t kSqrtSmall = 64;
static std::size_t sqrtN( const Limb * m , std::size_t k , Limb * s ,
Limb * r , bool root = false );
static std::uint64_t gcdWord( std::uint64_t a , std::uint64_t b );
static Limb gcdLimb( const Limb * a , std::size_t n , Limb d );
#if ADS_BIGINT_GMP
static mpz_srcptr gmpView( mpz_t z , const Limb * d , std::size_t n );
static Mag gmpDigits( mpz_srcptr z );
static bool gmpDivides( std::size_t u , std::size_t v , bool rem );
#endif
[[nodiscard]] bool normal( void ) const;
bool f_neg = false;
Mag f_d;
}; // end( class BigInteger )
/*--------------------------------------------------------------------------*/
} // end( namespace ADS )
/*--------------------------------------------------------------------------*/
/*-------------------------- IMPLEMENTATION --------------------------------*/
/*--------------------------------------------------------------------------*/
#include "ADS/numeric/impl/BigInteger.tpp"
/*--------------------------------------------------------------------------*/
#endif /* BigInteger.h included */
/*--------------------------------------------------------------------------*/
/*------------------------- End File BigInteger.h --------------------------*/
/*--------------------------------------------------------------------------*/