File NumberTheory.h

↰ Parent directory (ADS/numeric)

The arithmetic of section 17.2 and of exercises 17.1-17.3, on unsigned integers of 64 bits:

Definition (ADS/numeric/NumberTheory.h)

Detailed Description

  • gcd(), of Euclid, :math:` gcd( s , m ) = gcd( m , s bmod m ) , O( log m ) steps; extendedEuclid() also gives p and q with :math: s p + m q = gcd( s , m ) `: from ( d , p’ , q’ ) of ( m , s mod m ), ( d , q’ , p’ - lfloor s / m rfloor q’ ) (the solution of exercise 17.2 in the text has the two swapped);

  • modInverse(), the p with :math:` p s bmod m = 1 ` when :math:` gcd( s , m ) = 1 ` (exercise 17.3);

  • mulMod(), :math:` a b bmod m ` exact for every modulus of 64 bits, the modular product of all of Namespace ADS: the product of 64 bits for :math:` m le 2^{32} ; mulMod61() for the Mersenne prime :math: p = 2^{61} - 1 ` of the universal hashing [see File UniversalHashing.h], four products of 32 bits folded by :math:` 2^{61} equiv 1 pmod{p} `, faster than a division of 128 bits; the product of 128 bits where the compiler has one, and the doubling of Ahmes, O( log b ) additions, where it has not (MSVC);

  • powMod(), :math:` x^v bmod m ` by repeated squaring, O( log v ) multiplications (exercise 17.1), :math:` lfloor log_2 v rfloor + 1 ` squarings and as many products, each one of Montgomery when m is odd (see Montgomery64 below): the products and the squarings are two chains that do not depend on each other (from the lowest bit up), and the bit selects the product with no branch. Below :math:` 2^{62} ` the forms are kept in :math:` [ 0 , 2 m ) ` (the lazy reduction): for :math:` a , b < 2 m , :math: T = a b < 4 m^2 < m R , so that :math: ( T - q m ) / R + m in ( 0 , 2 m ) ` and the comparison that ends each product is spared; below :math:` 2^{32} ` the forms are of :math:` R = 2^{32} `, the products of 64 bits, the form of x by one division of 64 by 32 bits, the trailing zeros of v by squarings and the power from the lowest bit 1 of the odd part, so that the form of 1 (a division more) is never needed;

  • binaryGcd(), beyond the book, the binary gcd of Stein (Algorithm 4.5.2B of Knuth): with :math:` k ` the trailing zeros of :math:` a lor b , :math: gcd( a , b ) = 2^k gcd( a’ , b’ ) ` for the odd parts, and :math:` gcd( a , b ) = gcd( min( a , b ) , | a - b | ) ` with the trailing zeros of the even difference shifted out, by one instruction; each step halves at least :math:` max( a , b ) , hence at most :math: log_2 a + log_2 b + 2 ` steps and no division, the minimum and the difference by conditional moves, with no branch to mispredict; on x86-64 the trailing zeros are those of :math:` b - a ` (equal to those of :math:` | a - b | `), so that the count starts with the difference and the step is a subtraction, a count and a shift, the two conditional moves beside them;

  • Montgomery64, beyond the book, the arithmetic modulo an odd n of 64 bits in the form :math:` bar{x} = x R bmod n , :math: R = 2^{64} : the product is :math: mathrm{REDC}( bar{a} bar{b} ) = bar{a} bar{b} R^{-1} bmod n , where for :math: T = h R + l < n R ` and :math:` q = l n^{-1} bmod R ` the low half of :math:` q n ` equals l, so that :math:` ( T - q n ) / R = h - lfloor q n / R rfloor in ( - n , n ) ` and one addition of n when negative ends it: two products of 128 bits and one of 64 instead of a division of 128 bits, exact for every odd n; :math:` n^{-1} bmod 2^{64} ` comes from :math:` x_0 = 3 n oplus 2 , right on 5 bits, by 4 steps of Newton :math: x_{i+1} = x_i ( 2 - n x_i ) , each doubling the bits that are right, and :math: R bmod n ` and :math:` R^2 bmod n ` by one division of 128 by 64 bits each;

  • isPrimeTrial(), the divisions up to :math:` sqrt{n} , :math: O( 2^{d/2} ) ` for n of d bits;

  • witness(), b is a witness that n is composite if :math:` gcd( n , b ) neq 1 ` or, with :math:` n - 1 = 2^v m ` and m odd, neither :math:` b^m bmod n = 1 ` nor :math:` b^{2^i m} bmod n = n - 1 ` for some 0 <= i < v; millerRabin() tries K random bases in [ 2 , n - 2 ], at least 3/4 of which are witnesses for a composite n, and errs with probability at most :math:` ( 1 / 4 )^K `, only by saying prime;

  • isPrime(), beyond the book, exact on 64 bits: first the divisions by the 46 odd primes up to 211, by one product each with no branch, all at once on vectors (for the odd p, :math:` p mid n iff n p^{-1} bmod 2^{64} le lfloor ( 2^{64} - 1 ) / p rfloor , since :math: x mapsto x p^{-1} ` is a bijection that sends the multiples of p onto :math:` [ 0 , lfloor ( 2^{64} - 1 ) / p rfloor ] ), which leave about :math: prod_{3 le p le 211} ( 1 - 1 / p ) approx 0.21 ` of the odd numbers; then the test of Miller and Rabin on the fewest bases known to have no strong pseudoprime below n: below :math:` 2^{32} ` one base read in a table of 1024 at the hash :math:` lfloor ( n cdot 2654435761 bmod 2^{32} ) / 2^{22} rfloor ` (the hashed bases of Forisek and Jancina), the least base, at most 67, of which no odd composite of the bucket with no prime factor up to 211 is a strong pseudoprime (computed for Namespace ADS by the exhaustive search over the :math:` 2^{31} ` odd numbers), in an arithmetic of Montgomery of 32 bits that needs only :math:` n^{-1} bmod 2^{32} ` (the form of b by one division, the power from the lowest bit of d, the comparisons after one reduction), up to the 7 bases of Sinclair for every n of 64 bits. The k bases run in lockstep, k chains of products that the processor overlaps, so that 2 cost about as much as 1; with 3 or more, the exponent goes by windows of 4 bits with a table of :math:` b^0 , dots , b^{15} , :math: lceil log_2 d rceil + lceil log_2 d / 4 rceil + 14 ` products per base instead of up to :math:` 2 lceil log_2 d rceil . The arithmetic is that of Montgomery of 32 bits below :math: 2^{32} ` ( :math:` R = 2^{32} , its constants by divisions of 32 bits), the lazy one below :math: 2^{62} `, and Montgomery64 above;

  • factorize(), beyond the book, the prime factors of a word: the divisions by the primes below 256, then, on a composite cofactor m (isPrime() decides), the rho of Pollard with the cycle detection of Brent, in the Montgomery arithmetic of m: the sequence :math:` x_{i+1} = x_i^2 + c bmod m ` falls into a cycle modulo the least prime p of m after :math:` O( sqrt{p} ) = O( m^{1/4} ) ` steps in expectation (the birthday paradox, the map taken as random), which :math:` gcd( | x_i - x_j | , m ) ` finds; the differences are multiplied together, 128 at a time, before each gcd, and the last batch is gone through again one by one if the gcd is m; with another c if the factor is m itself. O( n^{1/4} ) products expected for n of 64 bits;

  • sieve(), the sieve of Eratosthenes on a SET of the integers from 2 to n [see File BoolSet.h], beyond the book: the least item p left in the set is a prime, and its multiples are removed; since the multiples of the primes below p are already gone, the removal can start at :math:` p^2 , and the sieve stops at :math: p > sqrt{n} , in :math: O( n log log n ) `.

The contracts [see File Contract.h]: extendedEuclid() needs s , m >= 0 and gives d = p s + q m; modInverse() needs m >= 1 and gcd( s , m ) = 1, and its result r satisfies s r = 1 modulo m; the products and the powers modulo m need m >= 1, and mulMod61() factors below :math:` 2^{61} - 1 `; witness() needs an odd n >= 3 and a base 1 <= b < n, and millerRabin() K >= 0; sieve() gives the primes up to n, increasing, as isPrimeTrial() finds them (an audit); factorize() needs n >= 1 and gives primes whose product is n (an audit). The functions are defined in NumberTheory.tpp.

  1. Montresor, Algoritmi e Strutture Dati, 2025, section 17.2 and exercises 17.1-17.3, https://algoritmicamente.it/materiali/libro/capitoli/cap17.pdf

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

  1. Stein, Computational problems associated with Racah algebra, Journal of Computational Physics 1(3), 397-405, 1967, doi:10.1016/0021-9991(67)90047-2

    1. Knuth, The Art of Computer Programming, vol. 2, 3rd ed., Addison-Wesley, 1998, section 4.5.2 (Algorithm B)

    1. Miller, Riemann’s hypothesis and tests for primality, Journal of Computer and System Sciences 13(3), 300-317, 1976, doi:10.1016/S0022-0000(76)80043-8

    1. Rabin, Probabilistic algorithm for testing primality, Journal of Number Theory 12(1), 128-138, 1980, doi:10.1016/0022-314X(80)90084-0

  1. Granlund, P. L. Montgomery, Division by invariant integers using multiplication, Proceedings of the ACM SIGPLAN 1994 Conference on Programming Language Design and Implementation (PLDI), 61-72, 1994, doi:10.1145/178243.178249

  1. Jaeschke, On strong pseudoprimes to several bases, Mathematics of Computation 61(204), 915-926, 1993, doi:10.1090/S0025-5718-1993-1192971-8

  1. Forisek, J. Jancina, Fast primality testing for integers that fit into a machine word, SOFSEM 2015 Student Research Forum, CEUR-WS 1326, 20-30, 2015, https://ceur-ws.org/Vol-1326/020-Forisek.pdf

  1. Sinclair and others, the deterministic sets of bases for the test of Miller and Rabin, https://miller-rabin.appspot.com

    1. Pollard, A Monte Carlo method for factorization, BIT 15(3), 331-334, 1975, doi:10.1007/BF01933667

    1. Brent, An improved Monte Carlo factorization algorithm, BIT Numerical Mathematics 20(2), 176-184, 1980, doi:10.1007/BF01933190

Donato Meoli

Dipartimento di Informatica Universita’ di Pisa

by Donato Meoli

Includes

Included By

Namespaces

Classes

Functions

Typedefs

Variables