File BigInteger.h

↰ Parent directory (ADS/numeric)

The integers of any length of Knuth (4.3.1), beyond the book: a sign and the magnitude as a vector of digits (limbs) in base :math:` b = 2^{64} `, the least significant first, with no leading zero, the numbers of one or two digits kept in place with no allocation.

Definition (ADS/numeric/BigInteger.h)

Detailed Description

The long operations work on spans of digits with the scratch space allocated once per call, and the digits are not zeroed before they are written. The product of two digits in two and the sums with carry are one instruction where the compiler offers them (unsigned __int128 in GCC and Clang, _umul128() and _addcarry_u64() in MSVC on x64), else the four products of the halves of 32 bits; ADS_BIGINT_PORTABLE forces the latter. Where GCC or Clang compile for x86-64 with BMI2 and ADX, the loops of the sums and of the products by one digit are in assembly, with MULX, which leaves the carry alone; ADS_BIGINT_NOASM keeps them in C++.

  • The sum and the difference are Algorithms 4.3.1A and 4.3.1S, a carry or a borrow per digit, :math:` O( n ) `.

  • The product is Algorithm 4.3.1M, :math:` O( m n ) , for the short factors: in the assembly by passes of k = 6 or 4 rows over the longer factor, a digit :math: a_i ` at a time, the window of the k + 1 places :math:` i , ldots , i + k ` in registers, :math:` r_i ` and the high halves of the k products :math:` a_i b_c ` added on the chain of carries of ADOX and their low halves on that of ADCX: as :math:` W + a_i B + r_i le ( b^k - 1 ) + ( b - 1 ) ( b^k - 1 ) + b - 1 = b^{k+1} - 1 ` for the window W and the k digits B of b, no carry leaves the window and both chains start again from 0 at the next digit. 31 instructions for 6 products, about 1.22 cycles per product on a Skylake, 1.4 on 32 digits; the rows left one at a time, and up to 8 digits the passes with no loop. Elsewhere two rows per pass in C++. A square :math:` a^2 = 2 T + D ` takes the products :math:` a_i a_j , :math: i < j , from 32 digits by blocks of 6: those of a block with the digits after it by the passes of 6 rows, those inside it by a kernel of 15 products in registers, all in T, and D, the squares :math: a_i^2 ` at :math:` b^{2i} `, added in the pass that doubles T; below 32 by tiles of 4 digits, the passes of 4 rows after each tile in T and the squares of the tiles side by side in D; the squares follow thresholds of their own (kSqr…). From kKaratsuba digits it is the product of Karatsuba (4.3.3A) in the subtractive form, :math:`[ a b = z_2 b^{2l} + bigl( z_0 + z_2 - ( a_0 - a_1 ) ( b_0 - b_1 )

    bigr) b^l + z_0 , qquad z_0 = a_0 b_0 , ; z_2 = a_1 b_1 , ` the sum of the middle in two passes on the two chains of carries, :math:` M( n ) = 3 M( n / 2 ) + O( n ) = O( n^{log_2 3} ) ; from kToom3 digits that of Toom-3, :math: a b = r( x ) , :math: x = b^k , from the five products :math: v_t = a( t ) b( t ) , :math: t = 0 , 1 , -1 , 2 , infty , and the interpolation of Bodrato and Zanoni, :math: r_3 = ( ( v_2 - v_{-1} ) / 3 - ( v_1 - v_0 ) ) / 2 -

    2 v_infty , :math: r_1 = ( v_1 - v_{-1} ) / 2 - r_3 , :math: r_2 = v_1 - v_0 - r_1 - r_3 - v_infty , the divisions exact, by 3 as the product by :math: -( b - 1 ) / 3 ` and a running difference, :math:` M( n ) = 5 M( n / 3 ) + O( n ) ; from kToom4 digits that of Toom-4, from the seven values at :math: 0 , pm 1 , pm 2 , 1 / 2 ` (times :math:` 2^6 ) and :math: infty , by an interpolation whose steps all stay on nonnegative numbers: with :math: E_1 = ( v_1 + v_{-1} ) / 2 , :math: O_1 = ( v_1 - v_{-1} ) / 2 , :math: E_2 = ( v_2 + v_{-2} ) / 2 ` and :math:` O_2 = ( v_2 - v_{-2} ) / 4 , :math: r_4 = ( ( E_2 - r_0 - 64 r_6 ) / 4 - ( E_1 - r_0 - r_6 ) ) / 3 , :math: r_2 = E_1 - r_0 - r_6 - r_4 , :math: H = ( v_{1/2} - 64 r_0 - 16 r_2 - 4 r_4 - r_6 ) / 2 = 16 r_1 + 4 r_3 + r_5 , :math: r_3 = ( 17 O_1 - O_2 - H ) / 9 , :math: r_5 = ( O_2 - O_1 - 3 r_3 ) / 15 ` and :math:` r_1 = O_1 - r_3 - r_5 , :math: M( n ) = 7 M( n / 4 ) + O( n ) = O( n^{log_4 7} ) ; from kToom32 digits a factor up to 7/4 as long as the other goes by Toom-3/2 (a in 3 pieces, b in 2, the values at 0, 1, -1 and :math: infty : :math: r_2 = ( v_1 + v_{-1} ) / 2 - r_0 , :math: r_1 = ( v_1 - v_{-1} ) / 2 - r_3 ), up to 5/2 by Toom-4/2 (also at 2: :math: r_3 = ( ( v_2 - r_0 - 4 r_2 - 16 r_4 ) / 2 - O ) / 3 ); the longer factor cut in pieces as long as the shorter one when their lengths differ; from kFFT digits it is a convolution of the digits by the number theoretic transform (4.3.3, 4.6.4) and the Chinese remainder theorem (4.3.2), :math: O( n log n ) , the top t digits of the longer factor left to rows of products by one digit when the transform of the rest is shorter by more than they cost, about :math: 12 L log_2 L ` cycles for L terms against :math:` 1.6 t n_b . With AVX2 and the fused product and sum, the transform runs on lanes of doubles, one digit per term, modulo three primes :math: p < 2^{50} ` with :math:` 3 cdot 2^{30} mid p - 1 , on :math: 2^k ` or :math:` 3 cdot 2^k ` terms, whichever is shorter, the residues exact integers in :math:` [ -p , p ] ` and the product of two of them reduced by the quotient rounded from the product by :math:` 1 / p , two rounds per pass over the terms: exact, since a coefficient, below :math: n , 2^{128} , stays below the product of the primes, :math: 2^{149} , up to :math: 2^{21} ` digits per factor. The coefficient i is :math:` y_1 + p_1 t_2 + p_1 p_2 t_3 ` by Garner, :math:` 0 le y_1 , t_2 , t_3 < 2^{50} : one digit each, so the product is :math: Y_1 + p_1 T_2 + p_1 p_2 T_3 ` for the numbers of those digits, three passes of products by one digit. Elsewhere and past that the digits are cut in two pieces of 32 bits each, modulo three primes :math:` p < 2^{31} ` with :math:` 2^{26} mid p - 1 , on 8 lanes of AVX2 or one at a time, the products by the roots of unity those of Shoup and the others those of Montgomery: exact up to transforms of :math: 2^{26} ` pieces, i.e., products of :math:` 2^{25} ` digits, past which Karatsuba splits the factors. multiplySchool(), multiplyKaratsuba() and multiplyFFT() are kept apart for the comparison.

  • The quotient and the remainder are Algorithm 4.3.1D for the short quotients: the divisor is shifted until its leading digit has the top bit set, and each trial quotient :math:` hat q ` is the quotient of the three leading digits of the remainder by the two of the divisor, by their reciprocal (Moeller and Granlund, Algorithms 5 and 6), which exceeds the true one by at most 1 and makes the correction after the multiplication and subtraction (“add back”) happen with probability about 2 / b; below kDivShort digits on the stack. From kDivBlock digits of divisor (kDivBlock6 for K = 6) the quotient goes K = 4 or 6 digits at a time: the K digits Q of the window, below :math:` v b^K , from its top :math: t + K ` digits by the top :math:` t = K + 2 ` of v, :math:` Q le hat Q le Q + 1 ` as :math:` hat Q - Q le 1 + Q / v_t < 2 , then the window less :math: hat Q v ` by one pass of K rows on the complements, :math:` r - a b = ( r_{top} - s ) b^n + overline{overline{r} + a b} ` with s the top digits of :math:` overline{r} + a b , and v added back while it is negative. A division of two digits by one, in the short division, is two products by a reciprocal of the divisor computed once, from a table and three steps of Newton (Moeller and Granlund, Algorithms 2 and 4). The remainder alone by a digit below :math: 2^{61} ` from kModFold digits folds 4 digits per step into a remainder of two digits, by the residues of :math:` b^k , with one division at the end; divExact() divides by one digit by its inverse modulo b (Hensel), from kDivExact2 digits both ways at once (Krandick and Jebelean): with :math: q = q_h b^m + q_l , :math: q_l d = ( a bmod b^m ) + c b^m ` for the last borrow :math:` c < d ` of Hensel’s division of the low m digits, so that :math:` lfloor a / b^m rfloor = q_h d + c ` and its quotient by d is :math:` q_h ; :math: m = n - n / 3 , as the chain of the division from the top is about twice as long per digit. From kDivide digits of quotient the division is recursive (Burnikel and Ziegler; Brent and Zimmermann, Algorithm 1.8): the upper half of the quotient from the upper half of the divisor, corrected at most twice, then the lower half, :math: O( M( n ) log n ) ` for products in :math:` M( n ) , a quotient longer than the divisor taken by blocks as long as the divisor. From kNewton digits of divisor, or kNewtonKept if its reciprocal is kept, it goes by blocks of about half the quotient, each block the product of the top digits of the remainder by a reciprocal of the top digits of the divisor (Newton’s method, Brent and Zimmermann, Algorithm 3.5; Barrett), short by a few units at most, and the product of the block by the divisor taken modulo :math: b^N - 1 ` by a cyclic convolution as long as the divisor, the digits that wrap known from the remainder and the lowest digit, :math:` O( M( n ) ) . The reciprocal :math: X = b^n + x ` of a of n digits, :math:` a X < b^{2n} le a ( X + 2 ) , is :math: X = X_h b^l + lfloor X_h t / b^{2h} rfloor ` from that of the top :math:` h = n - l ` digits, :math:` t = b^{n+h} - a X_h < 2 a : t from :math: a X_h ` modulo :math:` b^N - 1 , :math: N ge n , as t, its sum with :math: b^N - 1 ` and, :math:` X_h ` too large, its difference from it differ in the lowest digit, which is :math:` -a_0 x_0 bmod b ; and of t only its h + 1 digits from the place l are not 0, about :math: 1.8 ` to :math:` 2.1 M( n ) ` in all. The quotient alone, for operator/(), is first that of the top :math:` 2 m + 1 ` digits of the dividend by the top :math:` m + 2 ` of the divisor if it has m digits, :math:` m < n - 3 : with :math: v = v_t b^s + v_l ` and :math:` u = u_t b^s + u_l , the quotient :math: q_t ` of :math:` u_t ` by :math:` v_t ` is q or q + 1, as :math:` u_t / v_t - u_t / ( v_t + 1 ) < 2 / b , and it is q unless :math: r_t b^s + u_l < q_t v_l < b^{|q_t| + s} , which a remainder :math: r_t ` longer than :math:` q_t ` rules out with no product. Else it is that of the dividend with one digit more, so that the last digit tells when the quotient may be one off: below kDivCut digits of divisor by Algorithm D with the divisor cut to the j + 2 leading digits for the digit j, about half the work, and from kNewton digits with no remainder for the last block; the rare doubtful case is checked by a product. divModSchool() is Algorithm D alone. The quotient is truncated toward zero, as for the integers of C++.

  • The conversion to and from the bases 2 to 36 (4.4) works with the largest power P of the base in one digit, :math:` 10^{19} ` for the base

    1. In the bases that are powers of 2 it is a slice of bits per figure, :math:` O( n ) . In the others it is divide and conquer (Brent and Zimmermann, section 1.7), with the powers :math: P^{2^i} ` without their low digits 0 (about 30 per cent of the digits in the base 10) and their reciprocals kept per thread and base for the next conversions: the output splits the number by the power nearest the square root, the lower part padded with zeros, and the input joins two halves with one product per level, :math:` O( M( n ) log n ) ; below kConvert digits they divide or multiply by P once per many figures, as toStringSchool() and fromStringSchool() do at every length. On the output, a remainder :math: x < P^c , :math: c = 2^j ` chunks, below kConvertF digits goes by its fraction :math:` f = x / P^c ` (Bouvier and Zimmermann’s scaled remainder tree), the product by the reciprocal of :math:` P^c ` kept, rounded up, to the q + 2 digits of :math:` P^c ` plus 2: with :math:` f = ( X + theta ) / P^c , X the integer of the c chunks, the integer parts of :math: f P^i , :math: i le c , are those of :math: X / P^{c-i} ` while :math:` 0 le theta < 1 . The low half of the chunks is the fraction of :math: f P^{c/2} , the same :math: theta , by the band of the products that it needs (mulBand(), short by 2 units at most); the high half is the top digits of f, its :math: theta ` the fraction of the low half: f is rounded down if that fraction is at least 1/2, else up, and each half keeps its direction, so that :math:` theta ` stays below 1/2 by less than :math:` 3 b^{-2} ` per level if up and above it if down, and no run of 9 or of 0 carries wrong. From kFracLeaf chunks down the chunks are the integer parts of :math:` f P , :math: f P^2 `, …, f cut to the digits that the chunks left need, and the 19 figures of a chunk come by pairs from a table.

  • sqrtRem() is the Karatsuba square root of Zimmermann (also Brent and Zimmermann, Algorithm 1.12): m, shifted by an even number of bits to 2k digits with the top one at least b / 4, is :math:` A beta^2 + a_1 beta + a_0 ` with :math:` beta = b^{lfloor k / 2 rfloor} ; with :math: ( s’ , r’ ) ` the root and the remainder of A and :math:` ( q , u ) ` the quotient and the remainder of :math:` r’ beta + a_1 ` by :math:` 2 s’ , :math: s = s’ beta + q ` and :math:` r = u beta + a_0 - q^2 , if negative :math: r + 2 s - 1 ` and :math:` s - 1 , :math: O( M( n ) ) ` for the products of the divisions; up to kSqrtSmall digits on the stack, two digits from the root in doubles plus the nearest integer to :math:` ( a - x^2 ) / 2 x , then steps of 1. The remainder of a, shifted by 2t bits to m, is :math: ( r_m + 2 s_0 s_m - s_0^2 ) / 4^t ` for the t bits :math:` s_0 ` that the shift of the root drops. sqrt() needs no remainder: at the top :math:` l < h ` digits and the quotient alone of :math:` ( r’ beta + a_1 ) b ` by :math:` 2 s’ , :math: q b + d : :math: d ge 1 ` makes :math:` u ge 2 s’ / b ge beta , so :math: r ge 0 , with no remainder and no :math: q^2 , but with probability about :math: 1 / b `; sqrtNewton() is Newton’s method from above on the whole number.

  • pow() is the binary method of 4.6.3, from left to right, squares and products by x. gcd() is the algorithm of Lehmer (4.5.2L) on two digits: the leading 128 bits :math:` ( a_t , b_t ) ` of the two numbers run the algorithm of Euclid while the remainders keep more than 65 bits, so that by the lemma of Moeller its quotients are those of the numbers, :math:` ( a ; b ) = M ( alpha ; beta ) ` with :math:` det M = 1 ` and entries below :math:` 2^{62} , and one pass makes :math: ( alpha ; beta ) = ( u_{11} a - u_{01} b ; u_{00} b - u_{10} a ) , about 64 bits off both numbers; the pass has no subtraction, as :math: -b = overline{b} + 1 - b^n ` makes :math:` alpha equiv u_{11} a + u_{01} overline{b} + u_{01} ` modulo :math:` b^n , the two products of a digit summed below :math: 2^{128} . The steps run first on the leading digits alone down to remainders of :math: 2^{33} , then their matrix of entries below :math: 2^{31} ` reduces the 128 bits exactly and the leading 64 bits of the result run again: all in one digit, with the steps on each number written apart so that their state stays in registers; two digits only when the leading digits allow no step. From kGcdDC digits the top third goes to the half gcd of Moeller (after Knuth and Schoenhage): hgcd() of n digits reduces until |a - b| has :math:` s = lfloor n / 2 rfloor + 1 ` digits, from kHgcd digits by a recursive call on the top half, its matrix applied to the low digits, a few steps, and a second recursive call on the top of what is left, the two matrices multiplied (from kHgcdMul7 digits by the 7 products of Winograd’s form of Strassen’s, on signed numbers), :math:` O( M( n ) log n ) ; a step that the leading digits cannot make (a quotient of many digits) is a subtraction and a division. Down to two digits, where the binary algorithm ends; the number of one digit d = 2^t e, e odd, meets the long one by the folded remainder by e, or below kHensel digits by Hensel’s division, whose :math: a b^{-n} ` modulo e has the same gcd with e; gcdEuclid() is the algorithm of Euclid.

  • With ADS_WITH_GMP defined (by hand, or by the option of CMake of the same name), the products, the quotients, the conversions in the bases that are not powers of 2 and the greatest common divisor go to the mpn and mpz layers of GMP (mpn_mul(), mpn_sqr(), mpn_tdiv_qr(), mpz_tdiv_q(), mpz_divexact(), mpz_get_str(), mpn_set_str(), mpz_gcd()) in the ranges of lengths where GMP was the faster on a Xeon Gold 5120 (the thresholds kGmp…), on the digits of BigInteger as they are, read in place: for factors of :math:` n_a ge n_b ` digits, the product when :math:` k_{mul} le n_b < k_{top} , with :math: k_{top} ` kGmpMulTop if :math:` n_a < 2 n_b ` and kGmpMulLongTop otherwise, the square of n digits when :math:` n < ` kGmpSqrTop; for a divisor of n digits and a quotient of m or m + 1, the division with the remainder when :math:` n < ` kGmpDivTop or :math:` m < ` kGmpShort, the quotient alone when :math:` m > n ` and :math:` n < ` kGmpDivTop, or when :math:` m le n , :math: n ge ` kGmpQuot and either :math:` n < ` kGmpQuotTop or :math:` 2 m le n ; the exact quotient (by Jebelean’s algorithm, from the low digits) when :math: n < ` kGmpExactTop; the conversion of n digits to a base not a power of 2 when kConvert :math:` le n < ` kGmpToStringTop, and from one always; the gcd (by the subquadratic algorithm of Moeller) when the shorter argument has kGmpGcd digits or more. Elsewhere, and in the methods asked by name (multiplySchool(), multiplyKaratsuba(), multiplyFFT(), divModSchool(), toStringSchool(), fromStringSchool(), gcdEuclid()), BigInteger keeps its own code, which is all of it without the macro, when <gmp.h> is not even included; the same where a limb of GMP is not a word of 64 bits. GMP is free software under the LGPL v3 (or the GPL v2), as Namespace ADS is under the LGPL v3; note that GMP aborts where it cannot allocate, while BigInteger throws std::bad_alloc.

The contracts [see File Contract.h]: the divisor is not 0, the base is in [ 2 , 36 ] and the string a number in it; the invariant is the normal form (no leading zero, 0 not negative); divMod() ensures a = q b + r with |r| < |b| and r of the sign of a, and that it equals divModSchool(), operator/() that its quotient does, divExact() that its result times b is a; multiplyKaratsuba() and multiplyFFT() that they equal multiplySchool(); the conversions that they equal toStringSchool() and fromStringSchool(), gcd() that it equals gcdEuclid(), all audits. The operations are defined in BigInteger.tpp.

    1. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., Addison-Wesley, 1998, sections 4.3.1, 4.3.2, 4.3.3, 4.4, 4.5.2 and 4.6.3

  1. Karatsuba, Yu. Ofman, Multiplication of multidigit numbers on automata, Soviet Physics Doklady 7, 595-596, 1963

    1. Toom, The complexity of a scheme of functional elements realizing the multiplication of integers, Soviet Mathematics Doklady 3, 714-716, 1963

  1. Bodrato, A. Zanoni, Integer and polynomial multiplication: towards optimal Toom-Cook matrices, Proceedings of ISSAC 2007, 17-24, 2007, doi:10.1145/1277548.1277552

    1. Brent, P. Zimmermann, Modern Computer Arithmetic, Cambridge University Press, 2010, sections 1.3, 1.4, 1.6, 1.7, 2.4 and 3.4

  1. Burnikel, J. Ziegler, Fast recursive division, research report MPI-I-98-1-022, Max-Planck-Institut fuer Informatik, 1998

  1. Moeller, T. Granlund, Improved division by invariant integers, IEEE Transactions on Computers 60(2), 165-175, 2011, doi:10.1109/TC.2010.143

    1. Montgomery, Modular multiplication without trial division, Mathematics of Computation 44(170), 519-521, 1985, doi:10.1090/S0025-5718-1985-0777282-X

  1. Shoup, NTL: a library for doing number theory, the products by a constant with its precomputed quotient

  1. Harvey, Faster arithmetic for number-theoretic transforms, Journal of Symbolic Computation 60, 113-119, 2014, doi:10.1016/j.jsc.2013.09.002

  1. Jebelean, An algorithm for exact division, Journal of Symbolic Computation 15(2), 169-180, 1993, doi:10.1006/jsco.1993.1012

  1. Krandick, T. Jebelean, Bidirectional exact integer division, Journal of Symbolic Computation 21(4-6), 441-455, 1996, doi:10.1006/jsco.1996.0024

  1. Barrett, Implementing the Rivest Shamir and Adleman public key encryption algorithm on a standard digital signal processor, Proc. CRYPTO 1986, Lecture Notes in Computer Science 263, 311-323, 1987, doi:10.1007/3-540-47721-7_24

  1. Bouvier, P. Zimmermann, Division-free binary-to-decimal conversion, IEEE Transactions on Computers 63(8), 1895-1901, 2014, doi:10.1109/TC.2014.2315621

  1. Zimmermann, Karatsuba square root, INRIA research report RR-3805, 1999, https://hal.inria.fr/inria-00072854

    1. Lehmer, Euclid’s algorithm for large numbers, The American Mathematical Monthly 45(4), 227-233, 1938, doi:10.1080/00029890.1938.11990797

  1. Schoenhage, Schnelle Berechnung von Kettenbruchentwicklungen, Acta Informatica 1(2), 139-144, 1971, doi:10.1007/BF00289520

  1. Moeller, On Schoenhage’s algorithm and subquadratic integer gcd computation, Mathematics of Computation 77(261), 589-607, 2008, doi:10.1090/S0025-5718-07-02017-0

  1. Strassen, Gaussian elimination is not optimal, Numerische Mathematik 13(4), 354-356, 1969, doi:10.1007/BF02165411

  1. Winograd, On multiplication of 2 x 2 matrices, Linear Algebra and its Applications 4(4), 381-388, 1971, doi:10.1016/0024-3795(71)90009-7

  2. Granlund and the GMP development team, GNU MP: The GNU Multiple Precision Arithmetic Library, version 6.3.0, 2023, https://gmplib.org/gmp-man-6.3.0.pdf, the layers mpn and mpz of the backend

    1. Warren, Hacker’s Delight, 2nd ed., Addison-Wesley, 2012, sections 9.2 and 9.4

Donato Meoli

Dipartimento di Informatica Universita’ di Pisa

by Donato Meoli

Includes

Included By

Namespaces

Classes

Defines

Variables