Divide and Conquer Division
---------------------------

   For divisors larger than `DC_THRESHOLD', division is done by
dividing.  Or to be precise by a recursive divide and conquer algorithm
based on work by Moenck and Borodin, Jebelean, and Burnikel and Ziegler
(see References).

   The algorithm consists essentially of recognising that a 2NxN
division can be done with the basecase division algorithm (see
Basecase Division), but using N/2 limbs as a base, not just a single
limb.  This way the multiplications that arise are (N/2)x(N/2) and can
take advantage of Karatsuba and higher multiplication algorithms (see
Multiplication Algorithms).  The "digits" of the quotient are formed
by recursive Nx(N/2) divisions.

   If the (N/2)x(N/2) multiplies are done with a basecase multiplication
then the work is about the same as a basecase division, but with more
function call overheads and with some subtractions separated from the
multiplies.  These overheads mean that it's only when N/2 is above
`KARATSUBA_MUL_THRESHOLD' that divide and conquer is of use.

   `DC_THRESHOLD' is based on the divisor size N, so it will be
somewhere above twice `KARATSUBA_MUL_THRESHOLD', but how much above
depends on the CPU.  An optimized `mpn_mul_basecase' can lower
`DC_THRESHOLD' a little by offering a ready-made advantage over
repeated `mpn_submul_1' calls.

   Divide and conquer is asymptotically O(M(N)*log(N)) where M(N) is
the time for an NxN multiplication done with FFTs.  The actual time is
a sum over multiplications of the recursed sizes, as can be seen near
the end of section 2.2 of Burnikel and Ziegler.  For example, within
the Toom-3 range, divide and conquer is 2.63*M(N).  With higher
algorithms the M(N) term improves and the multiplier tends to log(N).
In practice, at moderate to large sizes, a 2NxN division is about 2 to
4 times slower than an NxN multiplication.

   Newton's method used for division is asymptotically O(M(N)) and
should therefore be superior to divide and conquer, but it's believed
this would only be for large to very large N.

Exact Division
--------------

   A so-called exact division is when the dividend is known to be an
exact multiple of the divisor.  Jebelean's exact division algorithm
uses this knowledge to make some significant optimizations (see
References).

   The idea can be illustrated in decimal for example with 368154
divided by 543.  Because the low digit of the dividend is 4, the low
digit of the quotient must be 8.  This is arrived at from 4*7 mod 10,
using the fact 7 is the modular inverse of 3 (the low digit of the
divisor), since 3*7 == 1 mod 10.  So 8*543=4344 can be subtracted from
the dividend leaving 363810.  Notice the low digit has become zero.

   The procedure is repeated at the second digit, with the next
quotient digit 7 (7 == 1*7 mod 10), subtracting 7*543=3801, leaving
325800.  And finally at the third digit with quotient digit 6 (8*7 mod
10), subtracting 6*543=3258 leaving 0.  So the quotient is 678.

   Notice however that the multiplies and subtractions don't need to
extend past the low three digits of the dividend, since that's enough
to determine the three quotient digits.  For the last quotient digit no
subtraction is needed at all.  On a 2NxN division like this one, only
about half the work of a normal basecase division is necessary.

   For an NxM exact division producing Q=N-M quotient limbs, the saving
over a normal basecase division is in two parts.  Firstly, each of the
Q quotient limbs needs only one multiply, not a 2x1 divide and
multiply.  Secondly, the crossproducts are reduced when Q>M to
Q*M-M*(M+1)/2, or when Q<=M to Q*(Q-1)/2.  Notice the savings are
complementary.  If Q is big then many divisions are saved, or if Q is
small then the crossproducts reduce to a small number.

   The modular inverse used is calculated efficiently by
`modlimb_invert' in `gmp-impl.h'.  This does four multiplies for a
32-bit limb, or six for a 64-bit limb.  `tune/modlinv.c' has some
alternate implementations that might suit processors better at bit
twiddling than multiplying.

   The sub-quadratic exact division described by Jebelean in "Exact
Division with Karatsuba Complexity" is not currently implemented.  It
uses a rearrangement similar to the divide and conquer for normal
division (see Divide and Conquer Division), but operating from low
to high.  A further possibility not currently implemented is
"Bidirectional Exact Integer Division" by Krandick and Jebelean which
forms quotient limbs from both the high and low ends of the dividend,
and can halve once more the number of crossproducts needed in a 2NxN
division.

   A special case exact division by 3 exists in `mpn_divexact_by3',
supporting Toom-3 multiplication and `mpq' canonicalizations.  It forms
quotient digits with a multiply by the modular inverse of 3 (which is
`0xAA..AAB') and uses two comparisons to determine a borrow for the next
limb.  The multiplications don't need to be on the dependent chain, as
long as the effect of the borrows is applied.  Only a few optimized
assembler implementations currently exist.

Exact Remainder
---------------

   If the exact division algorithm is done with a full subtraction at
each stage and the dividend isn't a multiple of the divisor, then low
zero limbs are produced but with a remainder in the high limbs.  For
dividend a, divisor d, quotient q, and b = 2^mp_bits_per_limb, then
this remainder r is of the form

     a = q*d + r*b^n

   n represents the number of zero limbs produced by the subtractions,
that being the number of limbs produced for q.  r will be in the range
0<=r<d and can be viewed as a remainder, but one shifted up by a factor
of b^n.

   Carrying out full subtractions at each stage means the same number
of cross products must be done as a normal division, but there's still
some single limb divisions saved.  When d is a single limb some
simplifications arise, providing good speedups on a number of
processors.

   `mpn_bdivmod', `mpn_divexact_by3', `mpn_modexact_1_odd' and the
`redc' function in `mpz_powm' differ subtly in how they return r,
leading to some negations in the above formula, but all are essentially
the same.

   Clearly r is zero when a is a multiple of d, and this leads to
divisibility or congruence tests which are potentially more efficient
than a normal division.

   The factor of b^n on r can be ignored in a GCD when d is odd, hence
the use of `mpn_bdivmod' in `mpn_gcd', and the use of
`mpn_modexact_1_odd' by `mpn_gcd_1' and `mpz_kronecker_ui' etc (see
Greatest Common Divisor Algorithms).

   Montgomery's REDC method for modular multiplications uses operands
of the form of x*b^-n and y*b^-n and on calculating (x*b^-n)*(y*b^-n)
uses the factor of b^n in the exact remainder to reach a product in the
same form (x*y)*b^-n (see Modular Powering Algorithm).

   Notice that r generally gives no useful information about the
ordinary remainder a mod d since b^n mod d could be anything.  If
however b^n == 1 mod d, then r is the negative of the ordinary
remainder.  This occurs whenever d is a factor of b^n-1, as for example
with 3 in `mpn_divexact_by3'.  Other such factors include 5, 17 and
257, but no particular use has been found for this.

Small Quotient Division
-----------------------

   An NxM division where the number of quotient limbs Q=N-M is small
can be optimized somewhat.

   An ordinary basecase division normalizes the divisor by shifting it
to make the high bit set, shifting the dividend accordingly, and
shifting the remainder back down at the end of the calculation.  This
is wasteful if only a few quotient limbs are to be formed.  Instead a
division of just the top 2*Q limbs of the dividend by the top Q limbs
of the divisor can be used to form a trial quotient.  This requires
only those limbs normalized, not the whole of the divisor and dividend.

   A multiply and subtract then applies the trial quotient to the M-Q
unused limbs of the divisor and N-Q dividend limbs (which includes Q
limbs remaining from the trial quotient division).  The starting trial
quotient can be 1 or 2 too big, but all cases of 2 too big and most
cases of 1 too big are detected by first comparing the most significant
limbs that will arise from the subtraction.  An addback is done if the
quotient still turns out to be 1 too big.

   This whole procedure is essentially the same as one step of the
basecase algorithm done in a Q limb base, though with the trial
quotient test done only with the high limbs, not an entire Q limb
"digit" product.  The correctness of this weaker test can be
established by following the argument of Knuth section 4.3.1 exercise
20 but with the v2*q>b*r+u2 condition appropriately relaxed.

Greatest Common Divisor
=======================

Binary GCD
----------

   At small sizes GMP uses an O(N^2) binary style GCD.  This is
described in many textbooks, for example Knuth section 4.5.2 algorithm
B.  It simply consists of successively reducing operands a and b using
gcd(a,b) = gcd(min(a,b),abs(a-b)), and also that if a and b are first
made odd then abs(a-b) is even and factors of two can be discarded.

   Variants like letting a-b become negative and doing a different next
step are of interest only as far as they suit particular CPUs, since on
small operands it's machine dependent factors that determine
performance.

   The Euclidean GCD algorithm, as per Knuth algorithms E and A,
reduces using a mod b but this has so far been found to be slower
everywhere.  One reason the binary method does well is that the implied
quotient at each step is usually small, so often only one or two
subtractions are needed to get the same effect as a division.
Quotients 1, 2 and 3 for example occur 67.7% of the time, see Knuth
section 4.5.3 Theorem E.

   When the implied quotient is large, meaning b is much smaller than
a, then a division is worthwhile.  This is the basis for the initial a
mod b reductions in `mpn_gcd' and `mpn_gcd_1' (the latter for both Nx1
and 1x1 cases).  But after that initial reduction, big quotients occur
too rarely to make it worth checking for them.

Accelerated GCD
---------------

   For sizes above `GCD_ACCEL_THRESHOLD', GMP uses the Accelerated GCD
algorithm described independently by Weber and Jebelean (the latter as
the "Generalized Binary" algorithm), see References.  This
algorithm is still O(N^2), but is much faster than the binary algorithm
since it does fewer multi-precision operations.  It consists of
alternating the k-ary reduction by Sorenson, and a "dmod" exact
remainder reduction.

   For operands u and v the k-ary reduction replaces u with n*v-d*u
where n and d are single limb values chosen to give two trailing zero
limbs on that value, which can be stripped.  n and d are calculated
using an algorithm similar to half of a two limb GCD (see `find_a' in
`mpn/generic/gcd.c').

   When u and v differ in size by more than a certain number of bits, a
dmod is performed to zero out bits at the low end of the larger.  It
consists of an exact remainder style division applied to an appropriate
number of bits (see Exact Division, and Exact Remainder).
This is faster than a k-ary reduction but useful only when the operands
differ in size.  There's a dmod after each k-ary reduction, and if the
dmod leaves the operands still differing in size then it's repeated.

   The k-ary reduction step can introduce spurious factors into the GCD
calculated, and these are eliminated at the end by taking GCDs with the
original inputs gcd(u,gcd(v,g)) using the binary algorithm.  Since g is
almost always small this takes very little time.

   At small sizes the algorithm needs a good implementation of
`find_a'.  At larger sizes it's dominated by `mpn_addmul_1' applying n
and d.

Extended GCD
------------

   The extended GCD calculates gcd(a,b) and also cofactors x and y
satisfying a*x+b*y=gcd(a,b).  Lehmer's multi-step improvement of the
extended Euclidean algorithm is used.  See Knuth section 4.5.2
algorithm L, and `mpn/generic/gcdext.c'.  This is an O(N^2) algorithm.

   The multipliers at each step are found using single limb
calculations for sizes up to `GCDEXT_THRESHOLD', or double limb
calculations above that.  The single limb code is faster but doesn't
produce full-limb multipliers, hence not making full use of the
`mpn_addmul_1' calls.

   When a CPU has a data-dependent multiplier, meaning one which is
faster on operands with fewer bits, the extra work in the double-limb
calculation might only save some looping overheads, leading to a large
`GCDEXT_THRESHOLD'.

   Currently the single limb calculation doesn't optimize for the small
quotients that often occur, and this can lead to unusually low values of
`GCDEXT_THRESHOLD', depending on the CPU.

   An analysis of double-limb calculations can be found in "A
Double-Digit Lehmer-Euclid Algorithm" by Jebelean (see References).
The code in GMP was developed independently.

   It should be noted that when a double limb calculation is used, it's
used for the whole of that GCD, it doesn't fall back to single limb
part way through.  This is because as the algorithm proceeds, the
inputs a and b are reduced, but the cofactors x and y grow, so the
multipliers at each step are applied to a roughly constant total number
of limbs.

Jacobi Symbol
-------------

   `mpz_jacobi' and `mpz_kronecker' are currently implemented with a
simple binary algorithm similar to that described for the GCDs (see
Binary GCD).  They're not very fast when both inputs are large.
Lehmer's multi-step improvement or a binary based multi-step algorithm
is likely to be better.

   When one operand fits a single limb, and that includes
`mpz_kronecker_ui' and friends, an initial reduction is done with
either `mpn_mod_1' or `mpn_modexact_1_odd', followed by the binary
algorithm on a single limb.  The binary algorithm is well suited to a
single limb, and the whole calculation in this case is quite efficient.

   In all the routines sign changes for the result are accumulated
using some bit twiddling, avoiding table lookups or conditional jumps.

Powering Algorithms
===================

Normal Powering
---------------

   Normal `mpz' or `mpf' powering uses a simple binary algorithm,
successively squaring and then multiplying by the base when a 1 bit is
seen in the exponent, as per Knuth section 4.6.3.  The "left to right"
variant described there is used rather than algorithm A, since it's
just as easy and can be done with somewhat less temporary memory.

Modular Powering
----------------

   Modular powering is implemented using a 2^k-ary sliding window
algorithm, as per "Handbook of Applied Cryptography" algorithm 14.85
(see References).  k is chosen according to the size of the
exponent.  Larger exponents use larger values of k, the choice being
made to minimize the average number of multiplications that must
supplement the squaring.

   The modular multiplies and squares use either a simple division or
the REDC method by Montgomery (see References).  REDC is a little
faster, essentially saving N single limb divisions in a fashion similar
to an exact remainder (see Exact Remainder).  The current REDC has
some limitations.  It's only O(N^2) so above `POWM_THRESHOLD' division
becomes faster and is used.  It doesn't attempt to detect small bases,
but rather always uses a REDC form, which is usually a full size
operand.  And lastly it's only applied to odd moduli.

Root Extraction Algorithms
==========================

Square Root
-----------

   Square roots are taken using the "Karatsuba Square Root" algorithm
by Paul Zimmermann (see References).  This is expressed in a divide
and conquer form, but as noted in the paper it can also be viewed as a
discrete variant of Newton's method.

   In the Karatsuba multiplication range this is an O(1.5*M(N/2))
algorithm, where M(n) is the time to multiply two numbers of n limbs.
In the FFT multiplication range this grows to a bound of O(6*M(N/2)).
In practice a factor of about 1.5 to 1.8 is found in the Karatsuba and
Toom-3 ranges, growing to 2 or 3 in the FFT range.

   The algorithm does all its calculations in integers and the resulting
`mpn_sqrtrem' is used for both `mpz_sqrt' and `mpf_sqrt'.  The extended
precision given by `mpf_sqrt_ui' is obtained by padding with zero limbs.

Nth Root
--------

   Integer Nth roots are taken using Newton's method with the following
iteration, where A is the input and n is the root to be taken.

              1         A
     a[i+1] = - * ( --------- + (n-1)*a[i] )
              n     a[i]^(n-1)

   The initial approximation a[1] is generated bitwise by successively
powering a trial root with or without new 1 bits, aiming to be just
above the true root.  The iteration converges quadratically when
started from a good approximation.  When n is large more initial bits
are needed to get good convergence.  The current implementation is not
particularly well optimized.

Perfect Square
--------------

   `mpz_perfect_square_p' is able to quickly exclude most non-squares by
checking whether the input is a quadratic residue modulo some small
integers.

   The first test is modulo 256 which means simply examining the least
significant byte.  Only 44 different values occur as the low byte of a
square, so 82.8% of non-squares can be immediately excluded.  Similar
tests modulo primes from 3 to 29 exclude 99.5% of those remaining, or
if a limb is 64 bits then primes up to 53 are used, excluding 99.99%.
A single Nx1 remainder using `PP' from `gmp-impl.h' quickly gives all
these remainders.

   A square root must still be taken for any value that passes the
residue tests, to verify it's really a square and not one of the 0.086%
(or 0.000156% for 64 bits) non-squares that get through.  See Square Root Algorithm.

Perfect Power
-------------

   Detecting perfect powers is required by some factorization
algorithms.  Currently `mpz_perfect_power_p' is implemented using
repeated Nth root extractions, though naturally only prime roots need
to be considered.  (see Nth Root Algorithm.)

   If a prime divisor p with multiplicity e can be found, then only
roots which are divisors of e need to be considered, much reducing the
work necessary.  To this end divisibility by a set of small primes is
checked.

Radix Conversion
================

   Radix conversions are less important than other algorithms.  A
program dominated by conversions should probably use a different data
representation.

Binary to Radix
---------------

   Conversions from binary to a power-of-2 radix use a simple and fast
O(N) bit extraction algorithm.

   Conversions from binary to other radices use repeated divisions,
first by the biggest power of the radix that fits in a single limb,
then by the radix on the remainders.  This is an O(N^2) algorithm and
can be quite time-consuming on large inputs.

Radix to Binary
---------------

   Conversions from a power-of-2 radix into binary use a simple and fast
O(N) bitwise concatenation algorithm.

   Conversions from other radices use repeated multiplications, first
accumulating as many digits as fit in a limb, then doing an Nx1
multi-precision multiplication.  This is O(N^2) and is certainly
sub-optimal on sizes above the Karatsuba multiply threshold.

Other Algorithms
================

Factorial
---------

   Factorials n! are calculated by a simple product from 1 to n, but
arranged into certain sub-products.

   First as many factors as fit in a limb are accumulated, then two of
those multiplied to give a 2-limb product.  When two 2-limb products
are ready they're multiplied to a 4-limb product, and when two 4-limbs
are ready they're multiplied to an 8-limb product, etc.  A stack of
outstanding products is built up, with two of the same size multiplied
together when ready.

   Arranging for multiplications to have operands the same (or nearly
the same) size means the Karatsuba and higher multiplication algorithms
can be used.  And even on sizes below the Karatsuba threshold an NxN
multiply will give a basecase multiply more to work on.

   An obvious improvement not currently implemented would be to strip
factors of 2 from the products and apply them at the end with a bit
shift.  Another possibility would be to determine the prime
factorization of the result (which can be done easily), and use a
powering method, at each stage squaring then multiplying in those
primes with a 1 in their exponent at that point.  The advantage would
be some multiplies turned into squares.

Binomial Coefficients
---------------------

   Binomial coefficients C(n,k) are calculated by first arranging k <=
n/2 using C(n,k) = C(n,n-k) if necessary, and then evaluating the
following product simply from i=2 to i=k.

                           k  (n-k+i)
     C(n,k) =  (n-k+1) * prod -------
                          i=2    i

   It's easy to show that each denominator i will divide the product so
far, so the exact division algorithm is used (see Exact Division).

   The numerators n-k+i and denominators i are first accumulated into
as many fit a limb, to save multi-precision operations, though for
`mpz_bin_ui' this applies only to the divisors, since n is an `mpz_t'
and n-k+i in general won't fit in a limb at all.

   An obvious improvement would be to strip factors of 2 from each
multiplier and divisor and count them separately, to be applied with a
bit shift at the end.  Factors of 3 and perhaps 5 could even be handled
similarly.  Another possibility, if n is not too big, would be to
determine the prime factorization of the result based on the factorials
involved, and power up those primes appropriately.  This would help
most when k is near n/2.

Fibonacci Numbers
-----------------

   The Fibonacci functions `mpz_fib_ui' and `mpz_fib2_ui' are designed
for calculating isolated F[n] or F[n],F[n-1] values efficiently.

   For small n, a table of single limb values in `__gmp_fib_table' is
used.  On a 32-bit limb this goes up to F[47], or on a 64-bit limb up
to F[93].  For convenience the table starts at F[-1].

   Beyond the table, values are generated with a binary powering
algorithm, calculating a pair F[n] and F[n-1] working from high to low
across the bits of n.  The formulas used are

     F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
     F[2k-1] =   F[k]^2 + F[k-1]^2

     F[2k] = F[2k+1] - F[2k-1]

   At each step, k is the high b bits of n.  If the next bit of n is 0
then F[2k],F[2k-1] is used, or if it's a 1 then F[2k+1],F[2k] is used,
and the process repeated until all bits of n are incorporated.  Notice
these formulas require just two squares per bit of n.

   It'd be possible to handle the first few n above the single limb
table with simple additions, using the defining Fibonacci recurrence
F[k+1]=F[k]+F[k-1], but this is not done since it usually turns out to
be faster for only about 10 or 20 values of n, and including a block of
code for just those doesn't seem worthwhile.  If they really mattered
it'd be better to extend the data table.

   Using a table avoids lots of calculations on small numbers, and
makes small n go fast.  A bigger table would make more small n go fast,
it's just a question of balancing size against desired speed.  For GMP
the code is kept compact, with the emphasis primarily on a good
powering algorithm.

   `mpz_fib2_ui' returns both F[n] and F[n-1], but `mpz_fib_ui' is only
interested in F[n].  In this case the last step of the algorithm can
become one multiply instead of two squares.  One of the following two
formulas is used, according as n is odd or even.

     F[2k]   = F[k]*(F[k]+2F[k-1])

     F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k

   F[2k+1] here is the same as above, just rearranged to be a multiply.
For interest, the 2*(-1)^k term both here and above can be applied
just to the low limb of the calculation, without a carry or borrow into
further limbs, which saves some code size.  See comments with
`mpz_fib_ui' and the internal `mpn_fib2_ui' for how this is done.

Lucas Numbers
-------------

   `mpz_lucnum2_ui' derives a pair of Lucas numbers from a pair of
Fibonacci numbers with the following simple formulas.

     L[k]   =   F[k] + 2*F[k-1]
     L[k-1] = 2*F[k] -   F[k-1]

   `mpz_lucnum_ui' is only interested in L[n], and some work can be
saved.  Trailing zero bits on n can be handled with a single square
each.

     L[2k] = L[k]^2 - 2*(-1)^k

   And the lowest 1 bit can be handled with one multiply of a pair of
Fibonacci numbers, similar to what `mpz_fib_ui' does.

     L[2k+1] = 5*F[k-1]*(2*F[k]+F[k-1]) - 4*(-1)^k

Assembler Coding
================

   The assembler subroutines in GMP are the most significant source of
speed at small to moderate sizes.  At larger sizes algorithm selection
becomes more important, but of course speedups in low level routines
will still speed up everything proportionally.

   Carry handling and widening multiplies that are important for GMP
can't be easily expressed in C.  GCC `asm' blocks help a lot and are
provided in `longlong.h', but hand coding low level routines invariably
offers a speedup over generic C by a factor of anything from 2 to 10.

Code Organisation
-----------------

   The various `mpn' subdirectories contain machine-dependent code,
written in C or assembler.  The `mpn/generic' subdirectory contains
default code, used when there's no machine-specific version of a
particular file.

   Each `mpn' subdirectory is for an ISA family.  Generally 32-bit and
64-bit variants in a family cannot share code and will have separate
directories.  Within a family further subdirectories may exist for CPU
variants.

Assembler Basics
----------------

   `mpn_addmul_1' and `mpn_submul_1' are the most important routines
for overall GMP performance.  All multiplications and divisions come
down to repeated calls to these.  `mpn_add_n', `mpn_sub_n',
`mpn_lshift' and `mpn_rshift' are next most important.

   On some CPUs assembler versions of the internal functions
`mpn_mul_basecase' and `mpn_sqr_basecase' give significant speedups,
mainly through avoiding function call overheads.  They can also
potentially make better use of a wide superscalar processor.

   The restrictions on overlaps between sources and destinations (see
Low-level Functions) are designed to facilitate a variety of
implementations.  For example, knowing `mpn_add_n' won't have partly
overlapping sources and destination means reading can be done far ahead
of writing on superscalar processors, and loops can be vectorized on a
vector processor, depending on the carry handling.

Carry Propagation
-----------------

   The problem that presents most challenges in GMP is propagating
carries from one limb to the next.  In functions like `mpn_addmul_1' and
`mpn_add_n', carries are the only dependencies between limb operations.

   On processors with carry flags, a straightforward CISC style `adc' is
generally best.  AMD K6 `mpn_addmul_1' however is an example of an
unusual set of circumstances where a branch works out better.

   On RISC processors generally an add and compare for overflow is
used.  This sort of thing can be seen in `mpn/generic/aors_n.c'.  Some
carry propagation schemes require 4 instructions, meaning at least 4
cycles per limb, but other schemes may use just 1 or 2.  On wide
superscalar processors performance may be completely determined by the
number of dependent instructions between carry-in and carry-out for
each limb.

   On vector processors good use can be made of the fact that a carry
bit only very rarely propagates more than one limb.  When adding a
single bit to a limb, there's only a carry out if that limb was
`0xFF...FF' which on random data will be only 1 in 2^mp_bits_per_limb.
`mpn/cray/add_n.c' is an example of this, it adds all limbs in
parallel, adds one set of carry bits in parallel and then only rarely
needs to fall through to a loop propagating further carries.

   On the x86s, GCC (as of version 2.95.2) doesn't generate
particularly good code for the RISC style idioms that are necessary to
handle carry bits in C.  Often conditional jumps are generated where
`adc' or `sbb' forms would be better.  And so unfortunately almost any
loop involving carry bits needs to be coded in assembler for best
results.

Cache Handling
--------------

   GMP aims to perform well both on operands that fit entirely in L1
cache and those that don't.  In the assembler subroutines this means
prefetching, either always or when large enough operands are presented.

   Pre-fetching sources combines well with loop unrolling, since a
prefetch can be initiated once per unrolled loop (or more than once if
the loop processes more than one cache line).

   Pre-fetching destinations won't be necessary if the CPU has a big
enough store queue.  Older processors without a write-allocate L1
however will want destination prefetching, to avoid repeated
write-throughs, unless they can keep up with the rate at which
destination limbs are produced.

   The distance ahead to prefetch will be determined by the rate data is
processed versus the time it takes to bring a line up to L1.  Naturally
the net data rate from L2 or RAM will always limit the rate of data
processing.  Prefetch distance may also be limited by the number of
prefetches the processor can have in progress at any one time.

   If a special prefetch instruction doesn't exist then a plain load
can be used, so long as the CPU supports out-of-order loads.  But this
may mean having a second copy of a loop so that the last few limbs can
be processed without prefetching, since reading past the end of an
operand must be avoided.

Floating Point
--------------

   Floating point arithmetic is used in GMP for multiplications on CPUs
with poor integer multipliers.  Floating point generally doesn't suit
other operations like additions or shifts, due to difficulties
implementing carry handling.

   With IEEE 53-bit double precision floats, integer multiplications
producing up to 53 bits will give exact results.  Breaking a
multiplication into 16x32->48 bit pieces is convenient.  With some care
though three 21x32->53 bit products can be used to do a 64x32 multiply,
if one of those 21x32 parts uses the sign bit.

   Generally limbs want to be treated as unsigned, but on some CPUs
floating point conversions only treat integers as signed.  Copying
through a zero extended memory region or testing and adjusting for a
sign bit may be necessary.

   Currently floating point FFTs aren't used for large multiplications.
On some processors they probably have a good chance of being
worthwhile, if great care is taken with precision control.

SIMD Instructions
-----------------

   The single-instruction multiple-data support in current
microprocessors is aimed at signal processing algorithms where each
data point can be treated more or less independently.  There's
generally not much support for propagating the sort of carries that
arise in GMP.

   SIMD multiplications of say four 16x16 bit multiplies only do as much
work as one 32x32 from GMP's point of view, and need some shifts and
adds besides.  But of course if say the SIMD form is fully pipelined
and uses less instruction decoding then it may still be worthwhile.

   On the 80x86 chips, MMX has so far found a use in `mpn_rshift' and
`mpn_lshift' since it allows 64-bit operations, and is used in a special
case for 16-bit multipliers in the P55 `mpn_mul_1'.  3DNow and SSE
haven't found a use so far.

Software Pipelining
-------------------

   Software pipelining consists of scheduling instructions around the
branch point in a loop.  For example a loop taking a checksum of an
array of limbs might have a load and an add, but the load wouldn't be
for that add, rather for the one next time around the loop.  Each load
then is effectively scheduled back in the previous iteration, allowing
latency to be hidden.

   Naturally this is wanted only when doing things like loads or
multiplies that take a few cycles to complete, and only where a CPU has
multiple functional units so that other work can be done while waiting.

   A pipeline with several stages will have a data value in progress at
each stage and each loop iteration moves them along one stage.  This is
like juggling.

   Within the loop some moves between registers may be necessary to
have the right values in the right places for each iteration.  Loop
unrolling can help this, with each unrolled block able to use different
registers for different values, even if some shuffling is still needed
just before going back to the top of the loop.

Loop Unrolling
--------------

   Loop unrolling consists of replicating code so that several limbs are
processed in each loop.  At a minimum this reduces loop overheads by a
corresponding factor, but it can also allow better register usage, for
example alternately using one register combination and then another.
Judicious use of `m4' macros can help avoid lots of duplication in the
source code.

   Unrolling is commonly done to a power of 2 multiple so the number of
unrolled loops and the number of remaining limbs can be calculated with
a shift and mask.  But other multiples can be used too, just by
subtracting each N limbs processed from a counter and waiting for less
than N remaining (or offsetting the counter by N so it goes negative
when there's less than N remaining).

   The limbs not a multiple of the unrolling can be handled in various
ways, for example

   * A simple loop at the end (or the start) to process the excess.
     Care will be wanted that it isn't too much slower than the
     unrolled part.

   * A set of binary tests, for example after an 8-limb unrolling, test
     for 4 more limbs to process, then a further 2 more or not, and
     finally 1 more or not.  This will probably take more code space
     than a simple loop.

   * A `switch' statement, providing separate code for each possible
     excess, for example an 8-limb unrolling would have separate code
     for 0 remaining, 1 remaining, etc, up to 7 remaining.  This might
     take a lot of code, but may be the best way to optimize all cases
     in combination with a deep pipelined loop.

   * A computed jump into the middle of the loop, thus making the first
     iteration handle the excess.  This should make times smoothly
     increase with size, which is attractive, but setups for the jump
     and adjustments for pointers can be tricky and could become quite
     difficult in combination with deep pipelining.

   One way to write the setups and finishups for a pipelined unrolled
loop is simply to duplicate the loop at the start and the end, then
delete instructions at the start which have no valid antecedents, and
delete instructions at the end whose results are unwanted.  Sizes not a
multiple of the unrolling can then be handled as desired.

Internals
*********

   *This chapter is provided only for informational purposes and the
various internals described here may change in future GMP releases.
Applications expecting to be compatible with future releases should use
only the documented interfaces described in previous chapters.*

Integer Internals
=================

   `mpz_t' variables represent integers using sign and magnitude, in
space dynamically allocated and reallocated.  The fields are as follows.

`_mp_size'
     The number of limbs, or the negative of that when representing a
     negative integer.  Zero is represented by `_mp_size' set to zero,
     in which case the `_mp_d' data is unused.

`_mp_d'
     A pointer to an array of limbs which is the magnitude.  These are
     stored "little endian" as per the `mpn' functions, so `_mp_d[0]'
     is the least significant limb and `_mp_d[ABS(_mp_size)-1]' is the
     most significant.  Whenever `_mp_size' is non-zero, the most
     significant limb is non-zero.

     Currently there's always at least one limb allocated, so for
     instance `mpz_set_ui' never needs to reallocate, and `mpz_get_ui'
     can fetch `_mp_d[0]' unconditionally (though its value is then
     only wanted if `_mp_size' is non-zero).

`_mp_alloc'
     `_mp_alloc' is the number of limbs currently allocated at `_mp_d',
     and naturally `_mp_alloc >= ABS(_mp_size)'.  When an `mpz' routine
     is about to (or might be about to) increase `_mp_size', it checks
     `_mp_alloc' to see whether there's enough space, and reallocates
     if not.  `MPZ_REALLOC' is generally used for this.

   The various bitwise logical functions like `mpz_and' behave as if
negative values were twos complement.  But sign and magnitude is always
used internally, and necessary adjustments are made during the
calculations.  Sometimes this isn't pretty, but sign and magnitude are
best for other routines.

   Some internal temporary variables are setup with `MPZ_TMP_INIT' and
these have `_mp_d' space obtained from `TMP_ALLOC' rather than the
memory allocation functions.  Care is taken to ensure that these are
big enough that no reallocation is necessary (since it would have
unpredictable consequences).

Rational Internals
==================

   `mpq_t' variables represent rationals using an `mpz_t' numerator and
denominator (see Integer Internals).

   The canonical form adopted is denominator positive (and non-zero),
no common factors between numerator and denominator, and zero uniquely
represented as 0/1.

   It's believed that casting out common factors at each stage of a
calculation is best in general.  A GCD is an O(N^2) operation so it's
better to do a few small ones immediately than to delay and have to do
a big one later.  Knowing the numerator and denominator have no common
factors can be used for example in `mpq_mul' to make only two cross
GCDs necessary, not four.

   This general approach to common factors is badly sub-optimal in the
presence of simple factorizations or little prospect for cancellation,
but GMP has no way to know when this will occur.  As per
Efficiency, that's left to applications.  The `mpq_t' framework might
still suit, with `mpq_numref' and `mpq_denref' for direct access to the
numerator and denominator, or of course `mpz_t' variables can be used
directly.