数学计算

开发平台：

Unix_Linux

gmp.texi：源码内容

@deftypefun mp_limb_t mpn_divmod (mp_limb_t *@var{r1p}, mp_limb_t *@var{rs2p}, mp_size_t @var{rs2n}, const mp_limb_t *@var{s3p}, mp_size_t @var{s3n})
[This function is obsolete. Please call @code{mpn_tdiv_qr} instead for best
performance.]
@end deftypefun
@deftypefn Macro mp_limb_t mpn_divexact_by3 (mp_limb_t *@var{rp}, mp_limb_t *@var{sp}, @w{mp_size_t @var{n}})
@deftypefnx Function mp_limb_t mpn_divexact_by3c (mp_limb_t *@var{rp}, mp_limb_t *@var{sp}, @w{mp_size_t @var{n}}, mp_limb_t @var{carry})
Divide @{@var{sp}, @var{n}@} by 3, expecting it to divide exactly, and writing
the result to @{@var{rp}, @var{n}@}. If 3 divides exactly, the return value is
zero and the result is the quotient. If not, the return value is non-zero and
the result won't be anything useful.
@code{mpn_divexact_by3c} takes an initial carry parameter, which can be the
return value from a previous call, so a large calculation can be done piece by
piece from low to high. @code{mpn_divexact_by3} is simply a macro calling
@code{mpn_divexact_by3c} with a 0 carry parameter.
These routines use a multiply-by-inverse and will be faster than
@code{mpn_divrem_1} on CPUs with fast multiplication but slow division.
The source @math{a}, result @math{q}, size @math{n}, initial carry @math{i},
and return value @math{c} satisfy @m{cb^n+a-i=3q, c*b^n + a-i = 3*q}, where
@m{b=2GMPraise{@code{GMP_NUMB_BITS}}, b=2^GMP_NUMB_BITS}. The
return @math{c} is always 0, 1 or 2, and the initial carry @math{i} must also
be 0, 1 or 2 (these are both borrows really). When @math{c=0} clearly
@math{q=(a-i)/3}. When @m{c neq 0, c!=0}, the remainder @math{(a-i) @bmod{}
3} is given by @math{3-c}, because @math{b @equiv{} 1 @bmod{} 3} (when
@code{mp_bits_per_limb} is even, which is always so currently).
@end deftypefn
@deftypefun mp_limb_t mpn_mod_1 (mp_limb_t *@var{s1p}, mp_size_t @var{s1n}, mp_limb_t @var{s2limb})
Divide @{@var{s1p}, @var{s1n}@} by @var{s2limb}, and return the remainder.
@var{s1n} can be zero.
@end deftypefun
@deftypefun mp_limb_t mpn_lshift (mp_limb_t *@var{rp}, const mp_limb_t *@var{sp}, mp_size_t @var{n}, unsigned int @var{count})
Shift @{@var{sp}, @var{n}@} left by @var{count} bits, and write the result to
@{@var{rp}, @var{n}@}. The bits shifted out at the left are returned in the
least significant @var{count} bits of the return value (the rest of the return
value is zero).
@var{count} must be in the range 1 to @nicode{mp_bits_per_limb}@minus{}1. The
regions @{@var{sp}, @var{n}@} and @{@var{rp}, @var{n}@} may overlap, provided
@math{@var{rp} @ge{} @var{sp}}.
This function is written in assembly for most CPUs.
@end deftypefun
@deftypefun mp_limb_t mpn_rshift (mp_limb_t *@var{rp}, const mp_limb_t *@var{sp}, mp_size_t @var{n}, unsigned int @var{count})
Shift @{@var{sp}, @var{n}@} right by @var{count} bits, and write the result to
@{@var{rp}, @var{n}@}. The bits shifted out at the right are returned in the
most significant @var{count} bits of the return value (the rest of the return
value is zero).
@var{count} must be in the range 1 to @nicode{mp_bits_per_limb}@minus{}1. The
regions @{@var{sp}, @var{n}@} and @{@var{rp}, @var{n}@} may overlap, provided
@math{@var{rp} @le{} @var{sp}}.
This function is written in assembly for most CPUs.
@end deftypefun
@deftypefun int mpn_cmp (const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Compare @{@var{s1p}, @var{n}@} and @{@var{s2p}, @var{n}@} and return a
positive value if @math{@var{s1} > @var{s2}}, 0 if they are equal, or a
negative value if @math{@var{s1} < @var{s2}}.
@end deftypefun
@deftypefun mp_size_t mpn_gcd (mp_limb_t *@var{rp}, mp_limb_t *@var{xp}, mp_size_t @var{xn}, mp_limb_t *@var{yp}, mp_size_t @var{yn})
Set @{@var{rp}, @var{retval}@} to the greatest common divisor of @{@var{xp},
@var{xn}@} and @{@var{yp}, @var{yn}@}. The result can be up to @var{yn} limbs,
the return value is the actual number produced. Both source operands are
destroyed.
@{@var{xp}, @var{xn}@} must have at least as many bits as @{@var{yp},
@var{yn}@}. @{@var{yp}, @var{yn}@} must be odd. Both operands must have
non-zero most significant limbs. No overlap is permitted between @{@var{xp},
@var{xn}@} and @{@var{yp}, @var{yn}@}.
@end deftypefun
@deftypefun mp_limb_t mpn_gcd_1 (const mp_limb_t *@var{xp}, mp_size_t @var{xn}, mp_limb_t @var{ylimb})
Return the greatest common divisor of @{@var{xp}, @var{xn}@} and @var{ylimb}.
Both operands must be non-zero.
@end deftypefun
@deftypefun mp_size_t mpn_gcdext (mp_limb_t *@var{gp}, mp_limb_t *@var{sp}, mp_size_t *@var{sn}, mp_limb_t *@var{xp}, mp_size_t @var{xn}, mp_limb_t *@var{yp}, mp_size_t @var{yn})
Let @m{U,@var{U}} be defined by @{@var{xp}, @var{xn}@} and let @m{V,@var{V}} be
defined by @{@var{yp}, @var{yn}@}.
Compute the greatest common divisor @math{G} of @math{U} and @math{V}. Compute
a cofactor @math{S} such that @math{G = US + VT}. The second cofactor @var{T}
is not computed but can easily be obtained from @m{(G - US) / V, (@var{G} -
@var{U}*@var{S}) / @var{V}} (the division will be exact). It is required that
@math{U @ge V > 0}.
@math{S} satisfies @math{S = 1} or @math{@GMPabs{S} < V / (2 G)}. @math{S =
0} if and only if @math{V} divides @math{U} (i.e., @math{G = V}).
Store @math{G} at @var{gp} and let the return value define its limb count.
Store @math{S} at @var{sp} and let |*@var{sn}| define its limb count. @math{S}
can be negative; when this happens *@var{sn} will be negative. The areas at
@var{gp} and @var{sp} should each have room for @math{@var{xn}+1} limbs.
The areas @{@var{xp}, @math{@var{xn}+1}@} and @{@var{yp}, @math{@var{yn}+1}@}
are destroyed (i.e.@: the input operands plus an extra limb past the end of
each).
Compatibility note: GMP 4.3.0 and 4.3.1 defined @math{S} less strictly.
Earlier as well as later GMP releases define @math{S} as described here.
@end deftypefun
@deftypefun mp_size_t mpn_sqrtrem (mp_limb_t *@var{r1p}, mp_limb_t *@var{r2p}, const mp_limb_t *@var{sp}, mp_size_t @var{n})
Compute the square root of @{@var{sp}, @var{n}@} and put the result at
@{@var{r1p}, @math{@GMPceil{@var{n}/2}}@} and the remainder at @{@var{r2p},
@var{retval}@}. @var{r2p} needs space for @var{n} limbs, but the return value
indicates how many are produced.
The most significant limb of @{@var{sp}, @var{n}@} must be non-zero. The
areas @{@var{r1p}, @math{@GMPceil{@var{n}/2}}@} and @{@var{sp}, @var{n}@} must
be completely separate. The areas @{@var{r2p}, @var{n}@} and @{@var{sp},
@var{n}@} must be either identical or completely separate.
If the remainder is not wanted then @var{r2p} can be @code{NULL}, and in this
case the return value is zero or non-zero according to whether the remainder
would have been zero or non-zero.
A return value of zero indicates a perfect square. See also
@code{mpz_perfect_square_p}.
@end deftypefun
@deftypefun mp_size_t mpn_get_str (unsigned char *@var{str}, int @var{base}, mp_limb_t *@var{s1p}, mp_size_t @var{s1n})
Convert @{@var{s1p}, @var{s1n}@} to a raw unsigned char array at @var{str} in
base @var{base}, and return the number of characters produced. There may be
leading zeros in the string. The string is not in ASCII; to convert it to
printable format, add the ASCII codes for @samp{0} or @samp{A}, depending on
the base and range. @var{base} can vary from 2 to 256.
The most significant limb of the input @{@var{s1p}, @var{s1n}@} must be
non-zero. The input @{@var{s1p}, @var{s1n}@} is clobbered, except when
@var{base} is a power of 2, in which case it's unchanged.
The area at @var{str} has to have space for the largest possible number
represented by a @var{s1n} long limb array, plus one extra character.
@end deftypefun
@deftypefun mp_size_t mpn_set_str (mp_limb_t *@var{rp}, const unsigned char *@var{str}, size_t @var{strsize}, int @var{base})
Convert bytes @{@var{str},@var{strsize}@} in the given @var{base} to limbs at
@var{rp}.
@math{@var{str}[0]} is the most significant byte and
@math{@var{str}[@var{strsize}-1]} is the least significant. Each byte should
be a value in the range 0 to @math{@var{base}-1}, not an ASCII character.
@var{base} can vary from 2 to 256.
The return value is the number of limbs written to @var{rp}. If the most
significant input byte is non-zero then the high limb at @var{rp} will be
non-zero, and only that exact number of limbs will be required there.
If the most significant input byte is zero then there may be high zero limbs
written to @var{rp} and included in the return value.
@var{strsize} must be at least 1, and no overlap is permitted between
@{@var{str},@var{strsize}@} and the result at @var{rp}.
@end deftypefun
@deftypefun {mp_bitcnt_t} mpn_scan0 (const mp_limb_t *@var{s1p}, mp_bitcnt_t @var{bit})
Scan @var{s1p} from bit position @var{bit} for the next clear bit.
It is required that there be a clear bit within the area at @var{s1p} at or
beyond bit position @var{bit}, so that the function has something to return.
@end deftypefun
@deftypefun {mp_bitcnt_t} mpn_scan1 (const mp_limb_t *@var{s1p}, mp_bitcnt_t @var{bit})
Scan @var{s1p} from bit position @var{bit} for the next set bit.
It is required that there be a set bit within the area at @var{s1p} at or
beyond bit position @var{bit}, so that the function has something to return.
@end deftypefun
@deftypefun void mpn_random (mp_limb_t *@var{r1p}, mp_size_t @var{r1n})
@deftypefunx void mpn_random2 (mp_limb_t *@var{r1p}, mp_size_t @var{r1n})
Generate a random number of length @var{r1n} and store it at @var{r1p}. The
most significant limb is always non-zero. @code{mpn_random} generates
uniformly distributed limb data, @code{mpn_random2} generates long strings of
zeros and ones in the binary representation.
@code{mpn_random2} is intended for testing the correctness of the @code{mpn}
routines.
@end deftypefun
@deftypefun {mp_bitcnt_t} mpn_popcount (const mp_limb_t *@var{s1p}, mp_size_t @var{n})
Count the number of set bits in @{@var{s1p}, @var{n}@}.
@end deftypefun
@deftypefun {mp_bitcnt_t} mpn_hamdist (const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Compute the hamming distance between @{@var{s1p}, @var{n}@} and @{@var{s2p},
@var{n}@}, which is the number of bit positions where the two operands have
different bit values.
@end deftypefun
@deftypefun int mpn_perfect_square_p (const mp_limb_t *@var{s1p}, mp_size_t @var{n})
Return non-zero iff @{@var{s1p}, @var{n}@} is a perfect square.
@end deftypefun
@deftypefun void mpn_and_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical and of @{@var{s1p}, @var{n}@} and @{@var{s2p},
@var{n}@}, and write the result to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_ior_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical inclusive or of @{@var{s1p}, @var{n}@} and
@{@var{s2p}, @var{n}@}, and write the result to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_xor_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical exclusive or of @{@var{s1p}, @var{n}@} and
@{@var{s2p}, @var{n}@}, and write the result to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_andn_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical and of @{@var{s1p}, @var{n}@} and the bitwise
complement of @{@var{s2p}, @var{n}@}, and write the result to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_iorn_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical inclusive or of @{@var{s1p}, @var{n}@} and the bitwise
complement of @{@var{s2p}, @var{n}@}, and write the result to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_nand_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical and of @{@var{s1p}, @var{n}@} and @{@var{s2p},
@var{n}@}, and write the bitwise complement of the result to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_nior_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical inclusive or of @{@var{s1p}, @var{n}@} and
@{@var{s2p}, @var{n}@}, and write the bitwise complement of the result to
@{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_xnor_n (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, const mp_limb_t *@var{s2p}, mp_size_t @var{n})
Perform the bitwise logical exclusive or of @{@var{s1p}, @var{n}@} and
@{@var{s2p}, @var{n}@}, and write the bitwise complement of the result to
@{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_com (mp_limb_t *@var{rp}, const mp_limb_t *@var{sp}, mp_size_t @var{n})
Perform the bitwise complement of @{@var{sp}, @var{n}@}, and write the result
to @{@var{rp}, @var{n}@}.
@end deftypefun
@deftypefun void mpn_copyi (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, mp_size_t @var{n})
Copy from @{@var{s1p}, @var{n}@} to @{@var{rp}, @var{n}@}, increasingly.
@end deftypefun
@deftypefun void mpn_copyd (mp_limb_t *@var{rp}, const mp_limb_t *@var{s1p}, mp_size_t @var{n})
Copy from @{@var{s1p}, @var{n}@} to @{@var{rp}, @var{n}@}, decreasingly.
@end deftypefun
@deftypefun void mpn_zero (mp_limb_t *@var{rp}, mp_size_t @var{n})
Zero @{@var{rp}, @var{n}@}.
@end deftypefun
@sp 1
@section Nails
@cindex Nails
@strong{Everything in this section is highly experimental and may disappear or
be subject to incompatible changes in a future version of GMP.}
Nails are an experimental feature whereby a few bits are left unused at the
top of each @code{mp_limb_t}. This can significantly improve carry handling
on some processors.
All the @code{mpn} functions accepting limb data will expect the nail bits to
be zero on entry, and will return data with the nails similarly all zero.
This applies both to limb vectors and to single limb arguments.
Nails can be enabled by configuring with @samp{--enable-nails}. By default
the number of bits will be chosen according to what suits the host processor,
but a particular number can be selected with @samp{--enable-nails=N}.
At the mpn level, a nail build is neither source nor binary compatible with a
non-nail build, strictly speaking. But programs acting on limbs only through
the mpn functions are likely to work equally well with either build, and
judicious use of the definitions below should make any program compatible with
either build, at the source level.
For the higher level routines, meaning @code{mpz} etc, a nail build should be
fully source and binary compatible with a non-nail build.
@defmac GMP_NAIL_BITS
@defmacx GMP_NUMB_BITS
@defmacx GMP_LIMB_BITS
@code{GMP_NAIL_BITS} is the number of nail bits, or 0 when nails are not in
use. @code{GMP_NUMB_BITS} is the number of data bits in a limb.
@code{GMP_LIMB_BITS} is the total number of bits in an @code{mp_limb_t}. In
all cases
@example
GMP_LIMB_BITS == GMP_NAIL_BITS + GMP_NUMB_BITS
@end example
@end defmac
@defmac GMP_NAIL_MASK
@defmacx GMP_NUMB_MASK
Bit masks for the nail and number parts of a limb. @code{GMP_NAIL_MASK} is 0
when nails are not in use.
@code{GMP_NAIL_MASK} is not often needed, since the nail part can be obtained
with @code{x >> GMP_NUMB_BITS}, and that means one less large constant, which
can help various RISC chips.
@end defmac
@defmac GMP_NUMB_MAX
The maximum value that can be stored in the number part of a limb. This is
the same as @code{GMP_NUMB_MASK}, but can be used for clarity when doing
comparisons rather than bit-wise operations.
@end defmac
The term ``nails'' comes from finger or toe nails, which are at the ends of a
limb (arm or leg). ``numb'' is short for number, but is also how the
developers felt after trying for a long time to come up with sensible names
for these things.
In the future (the distant future most likely) a non-zero nail might be
permitted, giving non-unique representations for numbers in a limb vector.
This would help vector processors since carries would only ever need to
propagate one or two limbs.
@node Random Number Functions, Formatted Output, Low-level Functions, Top
@chapter Random Number Functions
@cindex Random number functions
Sequences of pseudo-random numbers in GMP are generated using a variable of
type @code{gmp_randstate_t}, which holds an algorithm selection and a current
state. Such a variable must be initialized by a call to one of the
@code{gmp_randinit} functions, and can be seeded with one of the
@code{gmp_randseed} functions.
The functions actually generating random numbers are described in @ref{Integer
Random Numbers}, and @ref{Miscellaneous Float Functions}.
The older style random number functions don't accept a @code{gmp_randstate_t}
parameter but instead share a global variable of that type. They use a
default algorithm and are currently not seeded (though perhaps that will
change in the future). The new functions accepting a @code{gmp_randstate_t}
are recommended for applications that care about randomness.
@menu
* Random State Initialization::
* Random State Seeding::
* Random State Miscellaneous::
@end menu
@node Random State Initialization, Random State Seeding, Random Number Functions, Random Number Functions
@section Random State Initialization
@cindex Random number state
@cindex Initialization functions
@deftypefun void gmp_randinit_default (gmp_randstate_t @var{state})
Initialize @var{state} with a default algorithm. This will be a compromise
between speed and randomness, and is recommended for applications with no
special requirements. Currently this is @code{gmp_randinit_mt}.
@end deftypefun
@deftypefun void gmp_randinit_mt (gmp_randstate_t @var{state})
@cindex Mersenne twister random numbers
Initialize @var{state} for a Mersenne Twister algorithm. This algorithm is
fast and has good randomness properties.
@end deftypefun
@deftypefun void gmp_randinit_lc_2exp (gmp_randstate_t @var{state}, mpz_t @var{a}, @w{unsigned long @var{c}}, @w{mp_bitcnt_t @var{m2exp}})
@cindex Linear congruential random numbers
Initialize @var{state} with a linear congruential algorithm @m{X = (@var{a}X +
@var{c}) @bmod 2^{m2exp}, X = (@var{a}*X + @var{c}) mod 2^@var{m2exp}}.
The low bits of @math{X} in this algorithm are not very random. The least
significant bit will have a period no more than 2, and the second bit no more
than 4, etc. For this reason only the high half of each @math{X} is actually
used.
When a random number of more than @math{@var{m2exp}/2} bits is to be
generated, multiple iterations of the recurrence are used and the results
concatenated.
@end deftypefun
@deftypefun int gmp_randinit_lc_2exp_size (gmp_randstate_t @var{state}, mp_bitcnt_t @var{size})
@cindex Linear congruential random numbers
Initialize @var{state} for a linear congruential algorithm as per
@code{gmp_randinit_lc_2exp}. @var{a}, @var{c} and @var{m2exp} are selected
from a table, chosen so that @var{size} bits (or more) of each @math{X} will
be used, ie.@: @math{@var{m2exp}/2 @ge{} @var{size}}.
If successful the return value is non-zero. If @var{size} is bigger than the
table data provides then the return value is zero. The maximum @var{size}
currently supported is 128.
@end deftypefun
@deftypefun void gmp_randinit_set (gmp_randstate_t @var{rop}, gmp_randstate_t @var{op})
Initialize @var{rop} with a copy of the algorithm and state from @var{op}.
@end deftypefun
@c Although gmp_randinit, gmp_errno and related constants are obsolete, we
@c still put @findex entries for them, since they're still documented and
@c someone might be looking them up when perusing old application code.
@deftypefun void gmp_randinit (gmp_randstate_t @var{state}, @w{gmp_randalg_t @var{alg}}, @dots{})
@strong{This function is obsolete.}
@findex GMP_RAND_ALG_LC
@findex GMP_RAND_ALG_DEFAULT
Initialize @var{state} with an algorithm selected by @var{alg}. The only
choice is @code{GMP_RAND_ALG_LC}, which is @code{gmp_randinit_lc_2exp_size}
described above. A third parameter of type @code{unsigned long} is required,
this is the @var{size} for that function. @code{GMP_RAND_ALG_DEFAULT} or 0
are the same as @code{GMP_RAND_ALG_LC}.
@c For reference, this is the only place gmp_errno has been documented, and
@c due to being non thread safe we won't be adding to it's uses.
@findex gmp_errno
@findex GMP_ERROR_UNSUPPORTED_ARGUMENT
@findex GMP_ERROR_INVALID_ARGUMENT
@code{gmp_randinit} sets bits in the global variable @code{gmp_errno} to
indicate an error. @code{GMP_ERROR_UNSUPPORTED_ARGUMENT} if @var{alg} is
unsupported, or @code{GMP_ERROR_INVALID_ARGUMENT} if the @var{size} parameter
is too big. It may be noted this error reporting is not thread safe (a good
reason to use @code{gmp_randinit_lc_2exp_size} instead).
@end deftypefun
@deftypefun void gmp_randclear (gmp_randstate_t @var{state})
Free all memory occupied by @var{state}.
@end deftypefun
@node Random State Seeding, Random State Miscellaneous, Random State Initialization, Random Number Functions
@section Random State Seeding
@cindex Random number seeding
@cindex Seeding random numbers
@deftypefun void gmp_randseed (gmp_randstate_t @var{state}, mpz_t @var{seed})
@deftypefunx void gmp_randseed_ui (gmp_randstate_t @var{state}, @w{unsigned long int @var{seed}})
Set an initial seed value into @var{state}.
The size of a seed determines how many different sequences of random numbers
that it's possible to generate. The ``quality'' of the seed is the randomness
of a given seed compared to the previous seed used, and this affects the
randomness of separate number sequences. The method for choosing a seed is
critical if the generated numbers are to be used for important applications,
such as generating cryptographic keys.
Traditionally the system time has been used to seed, but care needs to be
taken with this. If an application seeds often and the resolution of the
system clock is low, then the same sequence of numbers might be repeated.
Also, the system time is quite easy to guess, so if unpredictability is
required then it should definitely not be the only source for the seed value.
On some systems there's a special device @file{/dev/random} which provides
random data better suited for use as a seed.
@end deftypefun
@node Random State Miscellaneous, , Random State Seeding, Random Number Functions
@section Random State Miscellaneous
@deftypefun {unsigned long} gmp_urandomb_ui (gmp_randstate_t @var{state}, unsigned long @var{n})
Return a uniformly distributed random number of @var{n} bits, ie.@: in the
range 0 to @m{2^n-1,2^@var{n}-1} inclusive. @var{n} must be less than or
equal to the number of bits in an @code{unsigned long}.
@end deftypefun
@deftypefun {unsigned long} gmp_urandomm_ui (gmp_randstate_t @var{state}, unsigned long @var{n})
Return a uniformly distributed random number in the range 0 to
@math{@var{n}-1}, inclusive.
@end deftypefun
@node Formatted Output, Formatted Input, Random Number Functions, Top
@chapter Formatted Output
@cindex Formatted output
@cindex @code{printf} formatted output
@menu
* Formatted Output Strings::
* Formatted Output Functions::
* C++ Formatted Output::
@end menu
@node Formatted Output Strings, Formatted Output Functions, Formatted Output, Formatted Output
@section Format Strings
@code{gmp_printf} and friends accept format strings similar to the standard C
@code{printf} (@pxref{Formatted Output,, Formatted Output, libc, The GNU C
Library Reference Manual}). A format specification is of the form
@example
% [flags] [width] [.[precision]] [type] conv
@end example
GMP adds types @samp{Z}, @samp{Q} and @samp{F} for @code{mpz_t}, @code{mpq_t}
and @code{mpf_t} respectively, @samp{M} for @code{mp_limb_t}, and @samp{N} for
an @code{mp_limb_t} array. @samp{Z}, @samp{Q}, @samp{M} and @samp{N} behave
like integers. @samp{Q} will print a @samp{/} and a denominator, if needed.
@samp{F} behaves like a float. For example,
@example
mpz_t z;
gmp_printf ("%s is an mpz %Zdn", "here", z);
mpq_t q;
gmp_printf ("a hex rational: %#40Qxn", q);
mpf_t f;
int n;
gmp_printf ("fixed point mpf %.*Ff with %d digitsn", n, f, n);
mp_limb_t l;
gmp_printf ("limb %Mun", l);
const mp_limb_t *ptr;
mp_size_t size;
gmp_printf ("limb array %Nxn", ptr, size);
@end example
For @samp{N} the limbs are expected least significant first, as per the
@code{mpn} functions (@pxref{Low-level Functions}). A negative size can be
given to print the value as a negative.
All the standard C @code{printf} types behave the same as the C library
@code{printf}, and can be freely intermixed with the GMP extensions. In the
current implementation the standard parts of the format string are simply
handed to @code{printf} and only the GMP extensions handled directly.
The flags accepted are as follows. GLIBC style @nisamp{'} is only for the
standard C types (not the GMP types), and only if the C library supports it.
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{0} @tab pad with zeros (rather than spaces)
@item @nicode{#} @tab show the base with @samp{0x}, @samp{0X} or @samp{0}
@item @nicode{+} @tab always show a sign
@item (space) @tab show a space or a @samp{-} sign
@item @nicode{'} @tab group digits, GLIBC style (not GMP types)
@end multitable
@end quotation
The optional width and precision can be given as a number within the format
string, or as a @samp{*} to take an extra parameter of type @code{int}, the
same as the standard @code{printf}.
The standard types accepted are as follows. @samp{h} and @samp{l} are
portable, the rest will depend on the compiler (or include files) for the type
and the C library for the output.
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{h} @tab @nicode{short}
@item @nicode{hh} @tab @nicode{char}
@item @nicode{j} @tab @nicode{intmax_t} or @nicode{uintmax_t}
@item @nicode{l} @tab @nicode{long} or @nicode{wchar_t}
@item @nicode{ll} @tab @nicode{long long}
@item @nicode{L} @tab @nicode{long double}
@item @nicode{q} @tab @nicode{quad_t} or @nicode{u_quad_t}
@item @nicode{t} @tab @nicode{ptrdiff_t}
@item @nicode{z} @tab @nicode{size_t}
@end multitable
@end quotation
@noindent
The GMP types are
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{F} @tab @nicode{mpf_t}, float conversions
@item @nicode{Q} @tab @nicode{mpq_t}, integer conversions
@item @nicode{M} @tab @nicode{mp_limb_t}, integer conversions
@item @nicode{N} @tab @nicode{mp_limb_t} array, integer conversions
@item @nicode{Z} @tab @nicode{mpz_t}, integer conversions
@end multitable
@end quotation
The conversions accepted are as follows. @samp{a} and @samp{A} are always
supported for @code{mpf_t} but depend on the C library for standard C float
types. @samp{m} and @samp{p} depend on the C library.
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{a} @nicode{A} @tab hex floats, C99 style
@item @nicode{c} @tab character
@item @nicode{d} @tab decimal integer
@item @nicode{e} @nicode{E} @tab scientific format float
@item @nicode{f} @tab fixed point float
@item @nicode{i} @tab same as @nicode{d}
@item @nicode{g} @nicode{G} @tab fixed or scientific float
@item @nicode{m} @tab @code{strerror} string, GLIBC style
@item @nicode{n} @tab store characters written so far
@item @nicode{o} @tab octal integer
@item @nicode{p} @tab pointer
@item @nicode{s} @tab string
@item @nicode{u} @tab unsigned integer
@item @nicode{x} @nicode{X} @tab hex integer
@end multitable
@end quotation
@samp{o}, @samp{x} and @samp{X} are unsigned for the standard C types, but for
types @samp{Z}, @samp{Q} and @samp{N} they are signed. @samp{u} is not
meaningful for @samp{Z}, @samp{Q} and @samp{N}.
@samp{M} is a proxy for the C library @samp{l} or @samp{L}, according to the
size of @code{mp_limb_t}. Unsigned conversions will be usual, but a signed
conversion can be used and will interpret the value as a twos complement
negative.
@samp{n} can be used with any type, even the GMP types.
Other types or conversions that might be accepted by the C library
@code{printf} cannot be used through @code{gmp_printf}, this includes for
instance extensions registered with GLIBC @code{register_printf_function}.
Also currently there's no support for POSIX @samp{$} style numbered arguments
(perhaps this will be added in the future).
The precision field has it's usual meaning for integer @samp{Z} and float
@samp{F} types, but is currently undefined for @samp{Q} and should not be used
with that.
@code{mpf_t} conversions only ever generate as many digits as can be
accurately represented by the operand, the same as @code{mpf_get_str} does.
Zeros will be used if necessary to pad to the requested precision. This
happens even for an @samp{f} conversion of an @code{mpf_t} which is an
integer, for instance @math{2^@W{1024}} in an @code{mpf_t} of 128 bits
precision will only produce about 40 digits, then pad with zeros to the
decimal point. An empty precision field like @samp{%.Fe} or @samp{%.Ff} can
be used to specifically request just the significant digits.
The decimal point character (or string) is taken from the current locale
settings on systems which provide @code{localeconv} (@pxref{Locales,, Locales
and Internationalization, libc, The GNU C Library Reference Manual}). The C
library will normally do the same for standard float output.
The format string is only interpreted as plain @code{char}s, multibyte
characters are not recognised. Perhaps this will change in the future.
@node Formatted Output Functions, C++ Formatted Output, Formatted Output Strings, Formatted Output
@section Functions
@cindex Output functions
Each of the following functions is similar to the corresponding C library
function. The basic @code{printf} forms take a variable argument list. The
@code{vprintf} forms take an argument pointer, see @ref{Variadic Functions,,
Variadic Functions, libc, The GNU C Library Reference Manual}, or @samp{man 3
va_start}.
It should be emphasised that if a format string is invalid, or the arguments
don't match what the format specifies, then the behaviour of any of these
functions will be unpredictable. GCC format string checking is not available,
since it doesn't recognise the GMP extensions.
The file based functions @code{gmp_printf} and @code{gmp_fprintf} will return
@math{-1} to indicate a write error. Output is not ``atomic'', so partial
output may be produced if a write error occurs. All the functions can return
@math{-1} if the C library @code{printf} variant in use returns @math{-1}, but
this shouldn't normally occur.
@deftypefun int gmp_printf (const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vprintf (const char *@var{fmt}, va_list @var{ap})
Print to the standard output @code{stdout}. Return the number of characters
written, or @math{-1} if an error occurred.
@end deftypefun
@deftypefun int gmp_fprintf (FILE *@var{fp}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vfprintf (FILE *@var{fp}, const char *@var{fmt}, va_list @var{ap})
Print to the stream @var{fp}. Return the number of characters written, or
@math{-1} if an error occurred.
@end deftypefun
@deftypefun int gmp_sprintf (char *@var{buf}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vsprintf (char *@var{buf}, const char *@var{fmt}, va_list @var{ap})
Form a null-terminated string in @var{buf}. Return the number of characters
written, excluding the terminating null.
No overlap is permitted between the space at @var{buf} and the string
@var{fmt}.
These functions are not recommended, since there's no protection against
exceeding the space available at @var{buf}.
@end deftypefun
@deftypefun int gmp_snprintf (char *@var{buf}, size_t @var{size}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vsnprintf (char *@var{buf}, size_t @var{size}, const char *@var{fmt}, va_list @var{ap})
Form a null-terminated string in @var{buf}. No more than @var{size} bytes
will be written. To get the full output, @var{size} must be enough for the
string and null-terminator.
The return value is the total number of characters which ought to have been
produced, excluding the terminating null. If @math{@var{retval} @ge{}
@var{size}} then the actual output has been truncated to the first
@math{@var{size}-1} characters, and a null appended.
No overlap is permitted between the region @{@var{buf},@var{size}@} and the
@var{fmt} string.
Notice the return value is in ISO C99 @code{snprintf} style. This is so even
if the C library @code{vsnprintf} is the older GLIBC 2.0.x style.
@end deftypefun
@deftypefun int gmp_asprintf (char **@var{pp}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vasprintf (char **@var{pp}, const char *@var{fmt}, va_list @var{ap})
Form a null-terminated string in a block of memory obtained from the current
memory allocation function (@pxref{Custom Allocation}). The block will be the
size of the string and null-terminator. The address of the block in stored to
*@var{pp}. The return value is the number of characters produced, excluding
the null-terminator.
Unlike the C library @code{asprintf}, @code{gmp_asprintf} doesn't return
@math{-1} if there's no more memory available, it lets the current allocation
function handle that.
@end deftypefun
@deftypefun int gmp_obstack_printf (struct obstack *@var{ob}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_obstack_vprintf (struct obstack *@var{ob}, const char *@var{fmt}, va_list @var{ap})
@cindex @code{obstack} output
Append to the current object in @var{ob}. The return value is the number of
characters written. A null-terminator is not written.
@var{fmt} cannot be within the current object in @var{ob}, since that object
might move as it grows.
These functions are available only when the C library provides the obstack
feature, which probably means only on GNU systems, see @ref{Obstacks,,
Obstacks, libc, The GNU C Library Reference Manual}.
@end deftypefun
@node C++ Formatted Output, , Formatted Output Functions, Formatted Output
@section C++ Formatted Output
@cindex C++ @code{ostream} output
@cindex @code{ostream} output
The following functions are provided in @file{libgmpxx} (@pxref{Headers and
Libraries}), which is built if C++ support is enabled (@pxref{Build Options}).
Prototypes are available from @code{<gmp.h>}.
@deftypefun ostream& operator<< (ostream& @var{stream}, mpz_t @var{op})
Print @var{op} to @var{stream}, using its @code{ios} formatting settings.
@code{ios::width} is reset to 0 after output, the same as the standard
@code{ostream operator<<} routines do.
In hex or octal, @var{op} is printed as a signed number, the same as for
decimal. This is unlike the standard @code{operator<<} routines on @code{int}
etc, which instead give twos complement.
@end deftypefun
@deftypefun ostream& operator<< (ostream& @var{stream}, mpq_t @var{op})
Print @var{op} to @var{stream}, using its @code{ios} formatting settings.
@code{ios::width} is reset to 0 after output, the same as the standard
@code{ostream operator<<} routines do.
Output will be a fraction like @samp{5/9}, or if the denominator is 1 then
just a plain integer like @samp{123}.
In hex or octal, @var{op} is printed as a signed value, the same as for
decimal. If @code{ios::showbase} is set then a base indicator is shown on
both the numerator and denominator (if the denominator is required).
@end deftypefun
@deftypefun ostream& operator<< (ostream& @var{stream}, mpf_t @var{op})
Print @var{op} to @var{stream}, using its @code{ios} formatting settings.
@code{ios::width} is reset to 0 after output, the same as the standard
@code{ostream operator<<} routines do.
The decimal point follows the standard library float @code{operator<<}, which
on recent systems means the @code{std::locale} imbued on @var{stream}.
Hex and octal are supported, unlike the standard @code{operator<<} on
@code{double}. The mantissa will be in hex or octal, the exponent will be in
decimal. For hex the exponent delimiter is an @samp{@@}. This is as per
@code{mpf_out_str}.
@code{ios::showbase} is supported, and will put a base on the mantissa, for
example hex @samp{0x1.8} or @samp{0x0.8}, or octal @samp{01.4} or @samp{00.4}.
This last form is slightly strange, but at least differentiates itself from
decimal.
@end deftypefun
These operators mean that GMP types can be printed in the usual C++ way, for
example,
@example
mpz_t z;
int n;
...
cout << "iteration " << n << " value " << z << "n";
@end example
But note that @code{ostream} output (and @code{istream} input, @pxref{C++
Formatted Input}) is the only overloading available for the GMP types and that
for instance using @code{+} with an @code{mpz_t} will have unpredictable
results. For classes with overloading, see @ref{C++ Class Interface}.
@node Formatted Input, C++ Class Interface, Formatted Output, Top
@chapter Formatted Input
@cindex Formatted input
@cindex @code{scanf} formatted input
@menu
* Formatted Input Strings::
* Formatted Input Functions::
* C++ Formatted Input::
@end menu
@node Formatted Input Strings, Formatted Input Functions, Formatted Input, Formatted Input
@section Formatted Input Strings
@code{gmp_scanf} and friends accept format strings similar to the standard C
@code{scanf} (@pxref{Formatted Input,, Formatted Input, libc, The GNU C
Library Reference Manual}). A format specification is of the form
@example
% [flags] [width] [type] conv
@end example
GMP adds types @samp{Z}, @samp{Q} and @samp{F} for @code{mpz_t}, @code{mpq_t}
and @code{mpf_t} respectively. @samp{Z} and @samp{Q} behave like integers.
@samp{Q} will read a @samp{/} and a denominator, if present. @samp{F} behaves
like a float.
GMP variables don't require an @code{&} when passed to @code{gmp_scanf}, since
they're already ``call-by-reference''. For example,
@example
/* to read say "a(5) = 1234" */
int n;
mpz_t z;
gmp_scanf ("a(%d) = %Zdn", &n, z);
mpq_t q1, q2;
gmp_sscanf ("0377 + 0x10/0x11", "%Qi + %Qi", q1, q2);
/* to read say "topleft (1.55,-2.66)" */
mpf_t x, y;
char buf[32];
gmp_scanf ("%31s (%Ff,%Ff)", buf, x, y);
@end example
All the standard C @code{scanf} types behave the same as in the C library
@code{scanf}, and can be freely intermixed with the GMP extensions. In the
current implementation the standard parts of the format string are simply
handed to @code{scanf} and only the GMP extensions handled directly.
The flags accepted are as follows. @samp{a} and @samp{'} will depend on
support from the C library, and @samp{'} cannot be used with GMP types.
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{*} @tab read but don't store
@item @nicode{a} @tab allocate a buffer (string conversions)
@item @nicode{'} @tab grouped digits, GLIBC style (not GMP types)
@end multitable
@end quotation
The standard types accepted are as follows. @samp{h} and @samp{l} are
portable, the rest will depend on the compiler (or include files) for the type
and the C library for the input.
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{h} @tab @nicode{short}
@item @nicode{hh} @tab @nicode{char}
@item @nicode{j} @tab @nicode{intmax_t} or @nicode{uintmax_t}
@item @nicode{l} @tab @nicode{long int}, @nicode{double} or @nicode{wchar_t}
@item @nicode{ll} @tab @nicode{long long}
@item @nicode{L} @tab @nicode{long double}
@item @nicode{q} @tab @nicode{quad_t} or @nicode{u_quad_t}
@item @nicode{t} @tab @nicode{ptrdiff_t}
@item @nicode{z} @tab @nicode{size_t}
@end multitable
@end quotation
@noindent
The GMP types are
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{F} @tab @nicode{mpf_t}, float conversions
@item @nicode{Q} @tab @nicode{mpq_t}, integer conversions
@item @nicode{Z} @tab @nicode{mpz_t}, integer conversions
@end multitable
@end quotation
The conversions accepted are as follows. @samp{p} and @samp{[} will depend on
support from the C library, the rest are standard.
@quotation
@multitable {(space)} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item @nicode{c} @tab character or characters
@item @nicode{d} @tab decimal integer
@item @nicode{e} @nicode{E} @nicode{f} @nicode{g} @nicode{G}
@tab float
@item @nicode{i} @tab integer with base indicator
@item @nicode{n} @tab characters read so far
@item @nicode{o} @tab octal integer
@item @nicode{p} @tab pointer
@item @nicode{s} @tab string of non-whitespace characters
@item @nicode{u} @tab decimal integer
@item @nicode{x} @nicode{X} @tab hex integer
@item @nicode{[} @tab string of characters in a set
@end multitable
@end quotation
@samp{e}, @samp{E}, @samp{f}, @samp{g} and @samp{G} are identical, they all
read either fixed point or scientific format, and either upper or lower case
@samp{e} for the exponent in scientific format.
C99 style hex float format (@code{printf %a}, @pxref{Formatted Output
Strings}) is always accepted for @code{mpf_t}, but for the standard float
types it will depend on the C library.
@samp{x} and @samp{X} are identical, both accept both upper and lower case
hexadecimal.
@samp{o}, @samp{u}, @samp{x} and @samp{X} all read positive or negative
values. For the standard C types these are described as ``unsigned''
conversions, but that merely affects certain overflow handling, negatives are
still allowed (per @code{strtoul}, @pxref{Parsing of Integers,, Parsing of
Integers, libc, The GNU C Library Reference Manual}). For GMP types there are
no overflows, so @samp{d} and @samp{u} are identical.
@samp{Q} type reads the numerator and (optional) denominator as given. If the
value might not be in canonical form then @code{mpq_canonicalize} must be
called before using it in any calculations (@pxref{Rational Number
Functions}).
@samp{Qi} will read a base specification separately for the numerator and
denominator. For example @samp{0x10/11} would be 16/11, whereas
@samp{0x10/0x11} would be 16/17.
@samp{n} can be used with any of the types above, even the GMP types.
@samp{*} to suppress assignment is allowed, though in that case it would do
nothing at all.
Other conversions or types that might be accepted by the C library
@code{scanf} cannot be used through @code{gmp_scanf}.
Whitespace is read and discarded before a field, except for @samp{c} and
@samp{[} conversions.
For float conversions, the decimal point character (or string) expected is
taken from the current locale settings on systems which provide
@code{localeconv} (@pxref{Locales,, Locales and Internationalization, libc,
The GNU C Library Reference Manual}). The C library will normally do the same
for standard float input.
The format string is only interpreted as plain @code{char}s, multibyte
characters are not recognised. Perhaps this will change in the future.
@node Formatted Input Functions, C++ Formatted Input, Formatted Input Strings, Formatted Input
@section Formatted Input Functions
@cindex Input functions
Each of the following functions is similar to the corresponding C library
function. The plain @code{scanf} forms take a variable argument list. The
@code{vscanf} forms take an argument pointer, see @ref{Variadic Functions,,
Variadic Functions, libc, The GNU C Library Reference Manual}, or @samp{man 3
va_start}.
It should be emphasised that if a format string is invalid, or the arguments
don't match what the format specifies, then the behaviour of any of these
functions will be unpredictable. GCC format string checking is not available,
since it doesn't recognise the GMP extensions.
No overlap is permitted between the @var{fmt} string and any of the results
produced.
@deftypefun int gmp_scanf (const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vscanf (const char *@var{fmt}, va_list @var{ap})
Read from the standard input @code{stdin}.
@end deftypefun
@deftypefun int gmp_fscanf (FILE *@var{fp}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vfscanf (FILE *@var{fp}, const char *@var{fmt}, va_list @var{ap})
Read from the stream @var{fp}.
@end deftypefun
@deftypefun int gmp_sscanf (const char *@var{s}, const char *@var{fmt}, @dots{})
@deftypefunx int gmp_vsscanf (const char *@var{s}, const char *@var{fmt}, va_list @var{ap})
Read from a null-terminated string @var{s}.
@end deftypefun
The return value from each of these functions is the same as the standard C99
@code{scanf}, namely the number of fields successfully parsed and stored.
@samp{%n} fields and fields read but suppressed by @samp{*} don't count
towards the return value.
If end of input (or a file error) is reached before a character for a field or
a literal, and if no previous non-suppressed fields have matched, then the
return value is @code{EOF} instead of 0. A whitespace character in the format
string is only an optional match and doesn't induce an @code{EOF} in this
fashion. Leading whitespace read and discarded for a field don't count as
characters for that field.
For the GMP types, input parsing follows C99 rules, namely one character of
lookahead is used and characters are read while they continue to meet the
format requirements. If this doesn't provide a complete number then the
function terminates, with that field not stored nor counted towards the return
value. For instance with @code{mpf_t} an input @samp{1.23e-XYZ} would be read
up to the @samp{X} and that character pushed back since it's not a digit. The
string @samp{1.23e-} would then be considered invalid since an @samp{e} must
be followed by at least one digit.
For the standard C types, in the current implementation GMP calls the C
library @code{scanf} functions, which might have looser rules about what
constitutes a valid input.
Note that @code{gmp_sscanf} is the same as @code{gmp_fscanf} and only does one
character of lookahead when parsing. Although clearly it could look at its
entire input, it is deliberately made identical to @code{gmp_fscanf}, the same
way C99 @code{sscanf} is the same as @code{fscanf}.
@node C++ Formatted Input, , Formatted Input Functions, Formatted Input
@section C++ Formatted Input
@cindex C++ @code{istream} input
@cindex @code{istream} input
The following functions are provided in @file{libgmpxx} (@pxref{Headers and
Libraries}), which is built only if C++ support is enabled (@pxref{Build
Options}). Prototypes are available from @code{<gmp.h>}.
@deftypefun istream& operator>> (istream& @var{stream}, mpz_t @var{rop})
Read @var{rop} from @var{stream}, using its @code{ios} formatting settings.
@end deftypefun
@deftypefun istream& operator>> (istream& @var{stream}, mpq_t @var{rop})
An integer like @samp{123} will be read, or a fraction like @samp{5/9}. No
whitespace is allowed around the @samp{/}. If the fraction is not in
canonical form then @code{mpq_canonicalize} must be called (@pxref{Rational
Number Functions}) before operating on it.
As per integer input, an @samp{0} or @samp{0x} base indicator is read when
none of @code{ios::dec}, @code{ios::oct} or @code{ios::hex} are set. This is
done separately for numerator and denominator, so that for instance
@samp{0x10/11} is @math{16/11} and @samp{0x10/0x11} is @math{16/17}.
@end deftypefun
@deftypefun istream& operator>> (istream& @var{stream}, mpf_t @var{rop})
Read @var{rop} from @var{stream}, using its @code{ios} formatting settings.
Hex or octal floats are not supported, but might be in the future, or perhaps
it's best to accept only what the standard float @code{operator>>} does.
@end deftypefun
Note that digit grouping specified by the @code{istream} locale is currently
not accepted. Perhaps this will change in the future.
@sp 1
These operators mean that GMP types can be read in the usual C++ way, for
example,
@example
mpz_t z;
...
cin >> z;
@end example
But note that @code{istream} input (and @code{ostream} output, @pxref{C++
Formatted Output}) is the only overloading available for the GMP types and
that for instance using @code{+} with an @code{mpz_t} will have unpredictable
results. For classes with overloading, see @ref{C++ Class Interface}.
@node C++ Class Interface, BSD Compatible Functions, Formatted Input, Top
@chapter C++ Class Interface
@cindex C++ interface
This chapter describes the C++ class based interface to GMP.
All GMP C language types and functions can be used in C++ programs, since
@file{gmp.h} has @code{extern "C"} qualifiers, but the class interface offers
overloaded functions and operators which may be more convenient.
Due to the implementation of this interface, a reasonably recent C++ compiler
is required, one supporting namespaces, partial specialization of templates
and member templates. For GCC this means version 2.91 or later.
@strong{Everything described in this chapter is to be considered preliminary
and might be subject to incompatible changes if some unforeseen difficulty
reveals itself.}
@menu
* C++ Interface General::
* C++ Interface Integers::
* C++ Interface Rationals::
* C++ Interface Floats::
* C++ Interface Random Numbers::
* C++ Interface Limitations::
@end menu
@node C++ Interface General, C++ Interface Integers, C++ Class Interface, C++ Class Interface
@section C++ Interface General
@noindent
All the C++ classes and functions are available with
@cindex @code{gmpxx.h}
@example
#include <gmpxx.h>
@end example
Programs should be linked with the @file{libgmpxx} and @file{libgmp}
libraries. For example,
@example
g++ mycxxprog.cc -lgmpxx -lgmp
@end example
@noindent
The classes defined are
@deftp Class mpz_class
@deftpx Class mpq_class
@deftpx Class mpf_class
@end deftp
The standard operators and various standard functions are overloaded to allow
arithmetic with these classes. For example,
@example
int
main (void)
@{
mpz_class a, b, c;
a = 1234;
b = "-5678";
c = a+b;
cout << "sum is " << c << "n";
cout << "absolute value is " << abs(c) << "n";
return 0;
@}
@end example
An important feature of the implementation is that an expression like
@code{a=b+c} results in a single call to the corresponding @code{mpz_add},
without using a temporary for the @code{b+c} part. Expressions which by their
nature imply intermediate values, like @code{a=b*c+d*e}, still use temporaries
though.
The classes can be freely intermixed in expressions, as can the classes and
the standard types @code{long}, @code{unsigned long} and @code{double}.
Smaller types like @code{int} or @code{float} can also be intermixed, since
C++ will promote them.
Note that @code{bool} is not accepted directly, but must be explicitly cast to
an @code{int} first. This is because C++ will automatically convert any
pointer to a @code{bool}, so if GMP accepted @code{bool} it would make all
sorts of invalid class and pointer combinations compile but almost certainly
not do anything sensible.
Conversions back from the classes to standard C++ types aren't done
automatically, instead member functions like @code{get_si} are provided (see
the following sections for details).
Also there are no automatic conversions from the classes to the corresponding
GMP C types, instead a reference to the underlying C object can be obtained
with the following functions,
@deftypefun mpz_t mpz_class::get_mpz_t ()
@deftypefunx mpq_t mpq_class::get_mpq_t ()
@deftypefunx mpf_t mpf_class::get_mpf_t ()
@end deftypefun
These can be used to call a C function which doesn't have a C++ class
interface. For example to set @code{a} to the GCD of @code{b} and @code{c},
@example
mpz_class a, b, c;
...
mpz_gcd (a.get_mpz_t(), b.get_mpz_t(), c.get_mpz_t());
@end example
In the other direction, a class can be initialized from the corresponding GMP
C type, or assigned to if an explicit constructor is used. In both cases this
makes a copy of the value, it doesn't create any sort of association. For
example,
@example
mpz_t z;
// ... init and calculate z ...
mpz_class x(z);
mpz_class y;
y = mpz_class (z);
@end example
There are no namespace setups in @file{gmpxx.h}, all types and functions are
simply put into the global namespace. This is what @file{gmp.h} has done in
the past, and continues to do for compatibility. The extras provided by
@file{gmpxx.h} follow GMP naming conventions and are unlikely to clash with
anything.
@node C++ Interface Integers, C++ Interface Rationals, C++ Interface General, C++ Class Interface
@section C++ Interface Integers
@deftypefun void mpz_class::mpz_class (type @var{n})
Construct an @code{mpz_class}. All the standard C++ types may be used, except
@code{long long} and @code{long double}, and all the GMP C++ classes can be
used. Any necessary conversion follows the corresponding C function, for
example @code{double} follows @code{mpz_set_d} (@pxref{Assigning Integers}).
@end deftypefun
@deftypefun void mpz_class::mpz_class (mpz_t @var{z})
Construct an @code{mpz_class} from an @code{mpz_t}. The value in @var{z} is
copied into the new @code{mpz_class}, there won't be any permanent association
between it and @var{z}.
@end deftypefun
@deftypefun void mpz_class::mpz_class (const char *@var{s})
@deftypefunx void mpz_class::mpz_class (const char *@var{s}, int @var{base} = 0)
@deftypefunx void mpz_class::mpz_class (const string& @var{s})
@deftypefunx void mpz_class::mpz_class (const string& @var{s}, int @var{base} = 0)
Construct an @code{mpz_class} converted from a string using @code{mpz_set_str}
(@pxref{Assigning Integers}).
If the string is not a valid integer, an @code{std::invalid_argument}
exception is thrown. The same applies to @code{operator=}.
@end deftypefun
@deftypefun mpz_class operator/ (mpz_class @var{a}, mpz_class @var{d})
@deftypefunx mpz_class operator% (mpz_class @var{a}, mpz_class @var{d})
Divisions involving @code{mpz_class} round towards zero, as per the
@code{mpz_tdiv_q} and @code{mpz_tdiv_r} functions (@pxref{Integer Division}).
This is the same as the C99 @code{/} and @code{%} operators.
The @code{mpz_fdiv@dots{}} or @code{mpz_cdiv@dots{}} functions can always be called
directly if desired. For example,
@example
mpz_class q, a, d;
...
mpz_fdiv_q (q.get_mpz_t(), a.get_mpz_t(), d.get_mpz_t());
@end example
@end deftypefun
@deftypefun mpz_class abs (mpz_class @var{op1})
@deftypefunx int cmp (mpz_class @var{op1}, type @var{op2})
@deftypefunx int cmp (type @var{op1}, mpz_class @var{op2})
@maybepagebreak
@deftypefunx bool mpz_class::fits_sint_p (void)
@deftypefunx bool mpz_class::fits_slong_p (void)
@deftypefunx bool mpz_class::fits_sshort_p (void)
@maybepagebreak
@deftypefunx bool mpz_class::fits_uint_p (void)
@deftypefunx bool mpz_class::fits_ulong_p (void)
@deftypefunx bool mpz_class::fits_ushort_p (void)
@maybepagebreak
@deftypefunx double mpz_class::get_d (void)
@deftypefunx long mpz_class::get_si (void)
@deftypefunx string mpz_class::get_str (int @var{base} = 10)
@deftypefunx {unsigned long} mpz_class::get_ui (void)
@maybepagebreak
@deftypefunx int mpz_class::set_str (const char *@var{str}, int @var{base})
@deftypefunx int mpz_class::set_str (const string& @var{str}, int @var{base})
@deftypefunx int sgn (mpz_class @var{op})
@deftypefunx mpz_class sqrt (mpz_class @var{op})
These functions provide a C++ class interface to the corresponding GMP C
routines.
@code{cmp} can be used with any of the classes or the standard C++ types,
except @code{long long} and @code{long double}.
@end deftypefun
@sp 1
Overloaded operators for combinations of @code{mpz_class} and @code{double}
are provided for completeness, but it should be noted that if the given
@code{double} is not an integer then the way any rounding is done is currently
unspecified. The rounding might take place at the start, in the middle, or at
the end of the operation, and it might change in the future.
Conversions between @code{mpz_class} and @code{double}, however, are defined
to follow the corresponding C functions @code{mpz_get_d} and @code{mpz_set_d}.
And comparisons are always made exactly, as per @code{mpz_cmp_d}.
@node C++ Interface Rationals, C++ Interface Floats, C++ Interface Integers, C++ Class Interface
@section C++ Interface Rationals
In all the following constructors, if a fraction is given then it should be in
canonical form, or if not then @code{mpq_class::canonicalize} called.
@deftypefun void mpq_class::mpq_class (type @var{op})
@deftypefunx void mpq_class::mpq_class (integer @var{num}, integer @var{den})
Construct an @code{mpq_class}. The initial value can be a single value of any
type, or a pair of integers (@code{mpz_class} or standard C++ integer types)
representing a fraction, except that @code{long long} and @code{long double}
are not supported. For example,
@example
mpq_class q (99);
mpq_class q (1.75);
mpq_class q (1, 3);
@end example
@end deftypefun
@deftypefun void mpq_class::mpq_class (mpq_t @var{q})
Construct an @code{mpq_class} from an @code{mpq_t}. The value in @var{q} is
copied into the new @code{mpq_class}, there won't be any permanent association
between it and @var{q}.
@end deftypefun
@deftypefun void mpq_class::mpq_class (const char *@var{s})
@deftypefunx void mpq_class::mpq_class (const char *@var{s}, int @var{base} = 0)
@deftypefunx void mpq_class::mpq_class (const string& @var{s})
@deftypefunx void mpq_class::mpq_class (const string& @var{s}, int @var{base} = 0)
Construct an @code{mpq_class} converted from a string using @code{mpq_set_str}
(@pxref{Initializing Rationals}).
If the string is not a valid rational, an @code{std::invalid_argument}
exception is thrown. The same applies to @code{operator=}.
@end deftypefun
@deftypefun void mpq_class::canonicalize ()
Put an @code{mpq_class} into canonical form, as per @ref{Rational Number
Functions}. All arithmetic operators require their operands in canonical
form, and will return results in canonical form.
@end deftypefun
@deftypefun mpq_class abs (mpq_class @var{op})
@deftypefunx int cmp (mpq_class @var{op1}, type @var{op2})
@deftypefunx int cmp (type @var{op1}, mpq_class @var{op2})
@maybepagebreak
@deftypefunx double mpq_class::get_d (void)
@deftypefunx string mpq_class::get_str (int @var{base} = 10)
@maybepagebreak
@deftypefunx int mpq_class::set_str (const char *@var{str}, int @var{base})
@deftypefunx int mpq_class::set_str (const string& @var{str}, int @var{base})
@deftypefunx int sgn (mpq_class @var{op})
These functions provide a C++ class interface to the corresponding GMP C
routines.
@code{cmp} can be used with any of the classes or the standard C++ types,
except @code{long long} and @code{long double}.
@end deftypefun
@deftypefun {mpz_class&} mpq_class::get_num ()
@deftypefunx {mpz_class&} mpq_class::get_den ()
Get a reference to an @code{mpz_class} which is the numerator or denominator
of an @code{mpq_class}. This can be used both for read and write access. If
the object returned is modified, it modifies the original @code{mpq_class}.
If direct manipulation might produce a non-canonical value, then
@code{mpq_class::canonicalize} must be called before further operations.
@end deftypefun
@deftypefun mpz_t mpq_class::get_num_mpz_t ()
@deftypefunx mpz_t mpq_class::get_den_mpz_t ()
Get a reference to the underlying @code{mpz_t} numerator or denominator of an
@code{mpq_class}. This can be passed to C functions expecting an
@code{mpz_t}. Any modifications made to the @code{mpz_t} will modify the
original @code{mpq_class}.
If direct manipulation might produce a non-canonical value, then
@code{mpq_class::canonicalize} must be called before further operations.
@end deftypefun
@deftypefun istream& operator>> (istream& @var{stream}, mpq_class& @var{rop});
Read @var{rop} from @var{stream}, using its @code{ios} formatting settings,
the same as @code{mpq_t operator>>} (@pxref{C++ Formatted Input}).
If the @var{rop} read might not be in canonical form then
@code{mpq_class::canonicalize} must be called.
@end deftypefun
@node C++ Interface Floats, C++ Interface Random Numbers, C++ Interface Rationals, C++ Class Interface
@section C++ Interface Floats
When an expression requires the use of temporary intermediate @code{mpf_class}
values, like @code{f=g*h+x*y}, those temporaries will have the same precision
as the destination @code{f}. Explicit constructors can be used if this
doesn't suit.
@deftypefun {} mpf_class::mpf_class (type @var{op})
@deftypefunx {} mpf_class::mpf_class (type @var{op}, unsigned long @var{prec})
Construct an @code{mpf_class}. Any standard C++ type can be used, except
@code{long long} and @code{long double}, and any of the GMP C++ classes can be
used.
If @var{prec} is given, the initial precision is that value, in bits. If
@var{prec} is not given, then the initial precision is determined by the type
of @var{op} given. An @code{mpz_class}, @code{mpq_class}, or C++
builtin type will give the default @code{mpf} precision (@pxref{Initializing
Floats}). An @code{mpf_class} or expression will give the precision of that
value. The precision of a binary expression is the higher of the two
operands.
@example
mpf_class f(1.5); // default precision
mpf_class f(1.5, 500); // 500 bits (at least)
mpf_class f(x); // precision of x
mpf_class f(abs(x)); // precision of x
mpf_class f(-g, 1000); // 1000 bits (at least)
mpf_class f(x+y); // greater of precisions of x and y
@end example
@end deftypefun
@deftypefun void mpf_class::mpf_class (const char *@var{s})
@deftypefunx void mpf_class::mpf_class (const char *@var{s}, unsigned long @var{prec}, int @var{base} = 0)
@deftypefunx void mpf_class::mpf_class (const string& @var{s})
@deftypefunx void mpf_class::mpf_class (const string& @var{s}, unsigned long @var{prec}, int @var{base} = 0)
Construct an @code{mpf_class} converted from a string using @code{mpf_set_str}
(@pxref{Assigning Floats}). If @var{prec} is given, the initial precision is
that value, in bits. If not, the default @code{mpf} precision
(@pxref{Initializing Floats}) is used.
If the string is not a valid float, an @code{std::invalid_argument} exception
is thrown. The same applies to @code{operator=}.
@end deftypefun
@deftypefun {mpf_class&} mpf_class::operator= (type @var{op})
Convert and store the given @var{op} value to an @code{mpf_class} object. The
same types are accepted as for the constructors above.
Note that @code{operator=} only stores a new value, it doesn't copy or change
the precision of the destination, instead the value is truncated if necessary.
This is the same as @code{mpf_set} etc. Note in particular this means for
@code{mpf_class} a copy constructor is not the same as a default constructor
plus assignment.
@example
mpf_class x (y); // x created with precision of y
mpf_class x; // x created with default precision
x = y; // value truncated to that precision
@end example
Applications using templated code may need to be careful about the assumptions
the code makes in this area, when working with @code{mpf_class} values of
various different or non-default precisions. For instance implementations of
the standard @code{complex} template have been seen in both styles above,
though of course @code{complex} is normally only actually specified for use
with the builtin float types.
@end deftypefun
@deftypefun mpf_class abs (mpf_class @var{op})
@deftypefunx mpf_class ceil (mpf_class @var{op})
@deftypefunx int cmp (mpf_class @var{op1}, type @var{op2})
@deftypefunx int cmp (type @var{op1}, mpf_class @var{op2})
@maybepagebreak
@deftypefunx bool mpf_class::fits_sint_p (void)
@deftypefunx bool mpf_class::fits_slong_p (void)
@deftypefunx bool mpf_class::fits_sshort_p (void)
@maybepagebreak
@deftypefunx bool mpf_class::fits_uint_p (void)
@deftypefunx bool mpf_class::fits_ulong_p (void)
@deftypefunx bool mpf_class::fits_ushort_p (void)
@maybepagebreak
@deftypefunx mpf_class floor (mpf_class @var{op})
@deftypefunx mpf_class hypot (mpf_class @var{op1}, mpf_class @var{op2})
@maybepagebreak
@deftypefunx double mpf_class::get_d (void)
@deftypefunx long mpf_class::get_si (void)
@deftypefunx string mpf_class::get_str (mp_exp_t& @var{exp}, int @var{base} = 10, size_t @var{digits} = 0)
@deftypefunx {unsigned long} mpf_class::get_ui (void)
@maybepagebreak
@deftypefunx int mpf_class::set_str (const char *@var{str}, int @var{base})
@deftypefunx int mpf_class::set_str (const string& @var{str}, int @var{base})
@deftypefunx int sgn (mpf_class @var{op})
@deftypefunx mpf_class sqrt (mpf_class @var{op})
@deftypefunx mpf_class trunc (mpf_class @var{op})
These functions provide a C++ class interface to the corresponding GMP C
routines.
@code{cmp} can be used with any of the classes or the standard C++ types,
except @code{long long} and @code{long double}.
The accuracy provided by @code{hypot} is not currently guaranteed.
@end deftypefun
@deftypefun {mp_bitcnt_t} mpf_class::get_prec ()
@deftypefunx void mpf_class::set_prec (mp_bitcnt_t @var{prec})
@deftypefunx void mpf_class::set_prec_raw (mp_bitcnt_t @var{prec})
Get or set the current precision of an @code{mpf_class}.
The restrictions described for @code{mpf_set_prec_raw} (@pxref{Initializing
Floats}) apply to @code{mpf_class::set_prec_raw}. Note in particular that the
@code{mpf_class} must be restored to it's allocated precision before being
destroyed. This must be done by application code, there's no automatic
mechanism for it.
@end deftypefun
@node C++ Interface Random Numbers, C++ Interface Limitations, C++ Interface Floats, C++ Class Interface
@section C++ Interface Random Numbers
@deftp Class gmp_randclass
The C++ class interface to the GMP random number functions uses
@code{gmp_randclass} to hold an algorithm selection and current state, as per
@code{gmp_randstate_t}.
@end deftp
@deftypefun {} gmp_randclass::gmp_randclass (void (*@var{randinit}) (gmp_randstate_t, @dots{}), @dots{})
Construct a @code{gmp_randclass}, using a call to the given @var{randinit}
function (@pxref{Random State Initialization}). The arguments expected are
the same as @var{randinit}, but with @code{mpz_class} instead of @code{mpz_t}.
For example,
@example
gmp_randclass r1 (gmp_randinit_default);
gmp_randclass r2 (gmp_randinit_lc_2exp_size, 32);
gmp_randclass r3 (gmp_randinit_lc_2exp, a, c, m2exp);
gmp_randclass r4 (gmp_randinit_mt);
@end example
@code{gmp_randinit_lc_2exp_size} will fail if the size requested is too big,
an @code{std::length_error} exception is thrown in that case.
@end deftypefun
@deftypefun {} gmp_randclass::gmp_randclass (gmp_randalg_t @var{alg}, @dots{})
Construct a @code{gmp_randclass} using the same parameters as
@code{gmp_randinit} (@pxref{Random State Initialization}). This function is
obsolete and the above @var{randinit} style should be preferred.
@end deftypefun
@deftypefun void gmp_randclass::seed (unsigned long int @var{s})
@deftypefunx void gmp_randclass::seed (mpz_class @var{s})
Seed a random number generator. See @pxref{Random Number Functions}, for how
to choose a good seed.
@end deftypefun
@deftypefun mpz_class gmp_randclass::get_z_bits (unsigned long @var{bits})
@deftypefunx mpz_class gmp_randclass::get_z_bits (mpz_class @var{bits})
Generate a random integer with a specified number of bits.
@end deftypefun
@deftypefun mpz_class gmp_randclass::get_z_range (mpz_class @var{n})
Generate a random integer in the range 0 to @math{@var{n}-1} inclusive.
@end deftypefun
@deftypefun mpf_class gmp_randclass::get_f ()
@deftypefunx mpf_class gmp_randclass::get_f (unsigned long @var{prec})
Generate a random float @var{f} in the range @math{0 <= @var{f} < 1}. @var{f}
will be to @var{prec} bits precision, or if @var{prec} is not given then to
the precision of the destination. For example,
@example
gmp_randclass r;
...
mpf_class f (0, 512); // 512 bits precision
f = r.get_f(); // random number, 512 bits
@end example
@end deftypefun
@node C++ Interface Limitations, , C++ Interface Random Numbers, C++ Class Interface
@section C++ Interface Limitations
@table @asis
@item @code{mpq_class} and Templated Reading
A generic piece of template code probably won't know that @code{mpq_class}
requires a @code{canonicalize} call if inputs read with @code{operator>>}
might be non-canonical. This can lead to incorrect results.
@code{operator>>} behaves as it does for reasons of efficiency. A
canonicalize can be quite time consuming on large operands, and is best
avoided if it's not necessary.
But this potential difficulty reduces the usefulness of @code{mpq_class}.
Perhaps a mechanism to tell @code{operator>>} what to do will be adopted in
the future, maybe a preprocessor define, a global flag, or an @code{ios} flag
pressed into service. Or maybe, at the risk of inconsistency, the
@code{mpq_class} @code{operator>>} could canonicalize and leave @code{mpq_t}
@code{operator>>} not doing so, for use on those occasions when that's
acceptable. Send feedback or alternate ideas to @email{gmp-bugs@@gmplib.org}.
@item Subclassing
Subclassing the GMP C++ classes works, but is not currently recommended.
Expressions involving subclasses resolve correctly (or seem to), but in normal
C++ fashion the subclass doesn't inherit constructors and assignments.
There's many of those in the GMP classes, and a good way to reestablish them
in a subclass is not yet provided.
@item Templated Expressions
A subtle difficulty exists when using expressions together with
application-defined template functions. Consider the following, with @code{T}
intended to be some numeric type,
@example
template <class T>
T fun (const T &, const T &);
@end example
@noindent
When used with, say, plain @code{mpz_class} variables, it works fine: @code{T}
is resolved as @code{mpz_class}.
@example
mpz_class f(1), g(2);
fun (f, g); // Good
@end example
@noindent
But when one of the arguments is an expression, it doesn't work.
@example
mpz_class f(1), g(2), h(3);
fun (f, g+h); // Bad
@end example
This is because @code{g+h} ends up being a certain expression template type
internal to @code{gmpxx.h}, which the C++ template resolution rules are unable
to automatically convert to @code{mpz_class}. The workaround is simply to add
an explicit cast.
@example
mpz_class f(1), g(2), h(3);
fun (f, mpz_class(g+h)); // Good
@end example
Similarly, within @code{fun} it may be necessary to cast an expression to type
@code{T} when calling a templated @code{fun2}.
@example
template <class T>
void fun (T f, T g)
@{
fun2 (f, f+g); // Bad
@}
template <class T>
void fun (T f, T g)
@{
fun2 (f, T(f+g)); // Good
@}
@end example
@end table
@node BSD Compatible Functions, Custom Allocation, C++ Class Interface, Top
@comment node-name, next, previous, up
@chapter Berkeley MP Compatible Functions
@cindex Berkeley MP compatible functions
@cindex BSD MP compatible functions
These functions are intended to be fully compatible with the Berkeley MP
library which is available on many BSD derived U*ix systems. The
@samp{--enable-mpbsd} option must be used when building GNU MP to make these
available (@pxref{Installing GMP}).
The original Berkeley MP library has a usage restriction: you cannot use the
same variable as both source and destination in a single function call. The
compatible functions in GNU MP do not share this restriction---inputs and
outputs may overlap.
It is not recommended that new programs are written using these functions.
Apart from the incomplete set of functions, the interface for initializing
@code{MINT} objects is more error prone, and the @code{pow} function collides
with @code{pow} in @file{libm.a}.
@cindex @code{mp.h}
@tindex MINT
Include the header @file{mp.h} to get the definition of the necessary types and
functions. If you are on a BSD derived system, make sure to include GNU
@file{mp.h} if you are going to link the GNU @file{libmp.a} to your program.
This means that you probably need to give the @samp{-I<dir>} option to the
compiler, where @samp{<dir>} is the directory where you have GNU @file{mp.h}.
@deftypefun {MINT *} itom (signed short int @var{initial_value})
Allocate an integer consisting of a @code{MINT} object and dynamic limb space.
Initialize the integer to @var{initial_value}. Return a pointer to the
@code{MINT} object.
@end deftypefun
@deftypefun {MINT *} xtom (char *@var{initial_value})
Allocate an integer consisting of a @code{MINT} object and dynamic limb space.
Initialize the integer from @var{initial_value}, a hexadecimal,
null-terminated C string. Return a pointer to the @code{MINT} object.
@end deftypefun
@deftypefun void move (MINT *@var{src}, MINT *@var{dest})
Set @var{dest} to @var{src} by copying. Both variables must be previously
initialized.
@end deftypefun
@deftypefun void madd (MINT *@var{src_1}, MINT *@var{src_2}, MINT *@var{destination})
Add @var{src_1} and @var{src_2} and put the sum in @var{destination}.
@end deftypefun
@deftypefun void msub (MINT *@var{src_1}, MINT *@var{src_2}, MINT *@var{destination})
Subtract @var{src_2} from @var{src_1} and put the difference in
@var{destination}.
@end deftypefun
@deftypefun void mult (MINT *@var{src_1}, MINT *@var{src_2}, MINT *@var{destination})
Multiply @var{src_1} and @var{src_2} and put the product in @var{destination}.
@end deftypefun
@deftypefun void mdiv (MINT *@var{dividend}, MINT *@var{divisor}, MINT *@var{quotient}, MINT *@var{remainder})
@deftypefunx void sdiv (MINT *@var{dividend}, signed short int @var{divisor}, MINT *@var{quotient}, signed short int *@var{remainder})
Set @var{quotient} to @var{dividend}/@var{divisor}, and @var{remainder} to
@var{dividend} mod @var{divisor}. The quotient is rounded towards zero; the
remainder has the same sign as the dividend unless it is zero.
Some implementations of these functions work differently---or not at all---for
negative arguments.
@end deftypefun
@deftypefun void msqrt (MINT *@var{op}, MINT *@var{root}, MINT *@var{remainder})
Set @var{root} to @m{lfloorsqrt{@var{op}}rfloor, the truncated integer part
of the square root of @var{op}}, like @code{mpz_sqrt}. Set @var{remainder} to
@m{(@var{op} - @var{root}^2), @var{op}@minus{}@var{root}*@var{root}}, i.e.
zero if @var{op} is a perfect square.
If @var{root} and @var{remainder} are the same variable, the results are
undefined.
@end deftypefun
@deftypefun void pow (MINT *@var{base}, MINT *@var{exp}, MINT *@var{mod}, MINT *@var{dest})
Set @var{dest} to (@var{base} raised to @var{exp}) modulo @var{mod}.
Note that the name @code{pow} clashes with @code{pow} from the standard C math
library (@pxref{Exponents and Logarithms,, Exponentiation and Logarithms,
libc, The GNU C Library Reference Manual}). An application will only be able
to use one or the other.
@end deftypefun
@deftypefun void rpow (MINT *@var{base}, signed short int @var{exp}, MINT *@var{dest})
Set @var{dest} to @var{base} raised to @var{exp}.
@end deftypefun
@deftypefun void gcd (MINT *@var{op1}, MINT *@var{op2}, MINT *@var{res})
Set @var{res} to the greatest common divisor of @var{op1} and @var{op2}.
@end deftypefun
@deftypefun int mcmp (MINT *@var{op1}, MINT *@var{op2})
Compare @var{op1} and @var{op2}. Return a positive value if @var{op1} >
@var{op2}, zero if @var{op1} = @var{op2}, and a negative value if @var{op1} <
@var{op2}.
@end deftypefun
@deftypefun void min (MINT *@var{dest})
Input a decimal string from @code{stdin}, and put the read integer in
@var{dest}. SPC and TAB are allowed in the number string, and are ignored.
@end deftypefun
@deftypefun void mout (MINT *@var{src})
Output @var{src} to @code{stdout}, as a decimal string. Also output a newline.
@end deftypefun
@deftypefun {char *} mtox (MINT *@var{op})
Convert @var{op} to a hexadecimal string, and return a pointer to the string.
The returned string is allocated using the default memory allocation function,
@code{malloc} by default. It will be @code{strlen(str)+1} bytes, that being
exactly enough for the string and null-terminator.
@end deftypefun
@deftypefun void mfree (MINT *@var{op})
De-allocate, the space used by @var{op}. @strong{This function should only be
passed a value returned by @code{itom} or @code{xtom}.}
@end deftypefun
@node Custom Allocation, Language Bindings, BSD Compatible Functions, Top
@comment node-name, next, previous, up
@chapter Custom Allocation
@cindex Custom allocation
@cindex Memory allocation
@cindex Allocation of memory
By default GMP uses @code{malloc}, @code{realloc} and @code{free} for memory
allocation, and if they fail GMP prints a message to the standard error output
and terminates the program.
Alternate functions can be specified, to allocate memory in a different way or
to have a different error action on running out of memory.
This feature is available in the Berkeley compatibility library (@pxref{BSD
Compatible Functions}) as well as the main GMP library.
@deftypefun void mp_set_memory_functions (@* void *(*@var{alloc_func_ptr}) (size_t), @* void *(*@var{realloc_func_ptr}) (void *, size_t, size_t), @* void (*@var{free_func_ptr}) (void *, size_t))
Replace the current allocation functions from the arguments. If an argument
is @code{NULL}, the corresponding default function is used.
These functions will be used for all memory allocation done by GMP, apart from
temporary space from @code{alloca} if that function is available and GMP is
configured to use it (@pxref{Build Options}).
@strong{Be sure to call @code{mp_set_memory_functions} only when there are no
active GMP objects allocated using the previous memory functions! Usually
that means calling it before any other GMP function.}
@end deftypefun
The functions supplied should fit the following declarations:
@deftypevr Function {void *} allocate_function (size_t @var{alloc_size})
Return a pointer to newly allocated space with at least @var{alloc_size}
bytes.
@end deftypevr
@deftypevr Function {void *} reallocate_function (void *@var{ptr}, size_t @var{old_size}, size_t @var{new_size})
Resize a previously allocated block @var{ptr} of @var{old_size} bytes to be
@var{new_size} bytes.
The block may be moved if necessary or if desired, and in that case the
smaller of @var{old_size} and @var{new_size} bytes must be copied to the new
location. The return value is a pointer to the resized block, that being the
new location if moved or just @var{ptr} if not.
@var{ptr} is never @code{NULL}, it's always a previously allocated block.
@var{new_size} may be bigger or smaller than @var{old_size}.
@end deftypevr
@deftypevr Function void free_function (void *@var{ptr}, size_t @var{size})
De-allocate the space pointed to by @var{ptr}.
@var{ptr} is never @code{NULL}, it's always a previously allocated block of
@var{size} bytes.
@end deftypevr
A @dfn{byte} here means the unit used by the @code{sizeof} operator.
The @var{old_size} parameters to @var{reallocate_function} and
@var{free_function} are passed for convenience, but of course can be ignored
if not needed. The default functions using @code{malloc} and friends for
instance don't use them.
No error return is allowed from any of these functions, if they return then
they must have performed the specified operation. In particular note that
@var{allocate_function} or @var{reallocate_function} mustn't return
@code{NULL}.
Getting a different fatal error action is a good use for custom allocation
functions, for example giving a graphical dialog rather than the default print
to @code{stderr}. How much is possible when genuinely out of memory is
another question though.
There's currently no defined way for the allocation functions to recover from
an error such as out of memory, they must terminate program execution. A
@code{longjmp} or throwing a C++ exception will have undefined results. This
may change in the future.
GMP may use allocated blocks to hold pointers to other allocated blocks. This
will limit the assumptions a conservative garbage collection scheme can make.
Since the default GMP allocation uses @code{malloc} and friends, those
functions will be linked in even if the first thing a program does is an
@code{mp_set_memory_functions}. It's necessary to change the GMP sources if
this is a problem.
@sp 1
@deftypefun void mp_get_memory_functions (@* void *(**@var{alloc_func_ptr}) (size_t), @* void *(**@var{realloc_func_ptr}) (void *, size_t, size_t), @* void (**@var{free_func_ptr}) (void *, size_t))
Get the current allocation functions, storing function pointers to the
locations given by the arguments. If an argument is @code{NULL}, that
function pointer is not stored.
@need 1000
For example, to get just the current free function,
@example
void (*freefunc) (void *, size_t);
mp_get_memory_functions (NULL, NULL, &freefunc);
@end example
@end deftypefun
@node Language Bindings, Algorithms, Custom Allocation, Top
@chapter Language Bindings
@cindex Language bindings
@cindex Other languages
The following packages and projects offer access to GMP from languages other
than C, though perhaps with varying levels of functionality and efficiency.
@c @spaceuref{U} is the same as @uref{U}, but with a couple of extra spaces
@c in tex, just to separate the URL from the preceding text a bit.
@iftex
@macro spaceuref {U}
@ @ @uref{U}
@end macro
@end iftex
@ifnottex
@macro spaceuref {U}
@uref{U}
@end macro
@end ifnottex
@sp 1
@table @asis
@item C++
@itemize @bullet
@item
GMP C++ class interface, @pxref{C++ Class Interface} @* Straightforward
interface, expression templates to eliminate temporaries.
@item
ALP @spaceuref{http://www-sop.inria.fr/saga/logiciels/ALP/} @* Linear algebra and
polynomials using templates.
@item
Arithmos @spaceuref{http://www.win.ua.ac.be/~cant/arithmos/} @* Rationals
with infinities and square roots.
@item
CLN @spaceuref{http://www.ginac.de/CLN/} @* High level classes for arithmetic.
@item
LiDIA @spaceuref{http://www.cdc.informatik.tu-darmstadt.de/TI/LiDIA/} @* A C++
library for computational number theory.
@item
Linbox @spaceuref{http://www.linalg.org/} @* Sparse vectors and matrices.
@item
NTL @spaceuref{http://www.shoup.net/ntl/} @* A C++ number theory library.
@end itemize
@c @item D
@c @itemize @bullet
@c @item
@c gmp-d @spaceuref{http://home.comcast.net/~benhinkle/gmp-d/}
@c @end itemize
@item Fortran
@itemize @bullet
@item
Omni F77 @spaceuref{http://phase.hpcc.jp/Omni/home.html} @* Arbitrary
precision floats.
@end itemize
@item Haskell
@itemize @bullet
@item
Glasgow Haskell Compiler @spaceuref{http://www.haskell.org/ghc/}
@end itemize
@item Java
@itemize @bullet
@item
Kaffe @spaceuref{http://www.kaffe.org/}
@item
Kissme @spaceuref{http://kissme.sourceforge.net/}
@end itemize
@item Lisp
@itemize @bullet
@item
GNU Common Lisp @spaceuref{http://www.gnu.org/software/gcl/gcl.html}
@item
Librep @spaceuref{http://librep.sourceforge.net/}
@item
@c FIXME: When there's a stable release with gmp support, just refer to it
@c rather than bothering to talk about betas.
XEmacs (21.5.18 beta and up) @spaceuref{http://www.xemacs.org} @* Optional
big integers, rationals and floats using GMP.
@end itemize
@item M4
@itemize @bullet
@item
@c FIXME: When there's a stable release with gmp support, just refer to it
@c rather than bothering to talk about betas.
GNU m4 betas @spaceuref{http://www.seindal.dk/rene/gnu/} @* Optionally provides
an arbitrary precision @code{mpeval}.
@end itemize
@item ML
@itemize @bullet
@item
MLton compiler @spaceuref{http://mlton.org/}
@end itemize
@item Objective Caml
@itemize @bullet
@item
MLGMP @spaceuref{http://www.di.ens.fr/~monniaux/programmes.html.en}
@item
Numerix @spaceuref{http://pauillac.inria.fr/~quercia/} @* Optionally using
GMP.
@end itemize
@item Oz
@itemize @bullet
@item
Mozart @spaceuref{http://www.mozart-oz.org/}
@end itemize
@item Pascal
@itemize @bullet
@item
GNU Pascal Compiler @spaceuref{http://www.gnu-pascal.de/} @* GMP unit.
@item
Numerix @spaceuref{http://pauillac.inria.fr/~quercia/} @* For Free Pascal,
optionally using GMP.
@end itemize
@item Perl
@itemize @bullet
@item
GMP module, see @file{demos/perl} in the GMP sources (@pxref{Demonstration
Programs}).
@item
Math::GMP @spaceuref{http://www.cpan.org/} @* Compatible with Math::BigInt, but
not as many functions as the GMP module above.
@item
Math::BigInt::GMP @spaceuref{http://www.cpan.org/} @* Plug Math::GMP into
normal Math::BigInt operations.
@end itemize
@need 1000
@item Pike
@itemize @bullet
@item
mpz module in the standard distribution, @uref{http://pike.ida.liu.se/}
@end itemize
@need 500
@item Prolog
@itemize @bullet
@item
SWI Prolog @spaceuref{http://www.swi-prolog.org/} @*
Arbitrary precision floats.
@end itemize
@item Python
@itemize @bullet
@item
mpz module in the standard distribution, @uref{http://www.python.org/}
@item
GMPY @uref{http://gmpy.sourceforge.net/}
@end itemize
@item Scheme
@itemize @bullet
@item
GNU Guile (upcoming 1.8) @spaceuref{http://www.gnu.org/software/guile/guile.html}
@item
RScheme @spaceuref{http://www.rscheme.org/}
@item
STklos @spaceuref{http://www.stklos.org/}
@c
@c For reference, MzScheme uses some of gmp, but (as of version 205) it only
@c has copies of some of the generic C code, and we don't consider that a
@c language binding to gmp.
@c
@end itemize
@item Smalltalk
@itemize @bullet
@item
GNU Smalltalk @spaceuref{http://www.smalltalk.org/versions/GNUSmalltalk.html}
@end itemize
@item Other
@itemize @bullet
@item
Axiom @uref{http://savannah.nongnu.org/projects/axiom} @* Computer algebra
using GCL.
@item
DrGenius @spaceuref{http://drgenius.seul.org/} @* Geometry system and
mathematical programming language.
@item
GiNaC @spaceuref{http://www.ginac.de/} @* C++ computer algebra using CLN.
@item
GOO @spaceuref{http://www.googoogaga.org/} @* Dynamic object oriented
language.
@item
Maxima @uref{http://www.ma.utexas.edu/users/wfs/maxima.html} @* Macsyma
computer algebra using GCL.
@item
Q @spaceuref{http://q-lang.sourceforge.net/} @* Equational programming system.
@item
Regina @spaceuref{http://regina.sourceforge.net/} @* Topological calculator.
@item
Yacas @spaceuref{http://www.xs4all.nl/~apinkus/yacas.html} @* Yet another
computer algebra system.
@end itemize
@end table
@node Algorithms, Internals, Language Bindings, Top
@chapter Algorithms
@cindex Algorithms
This chapter is an introduction to some of the algorithms used for various GMP
operations. The code is likely to be hard to understand without knowing
something about the algorithms.
Some GMP internals are mentioned, but applications that expect to be
compatible with future GMP releases should take care to use only the
documented functions.
@menu
* Multiplication Algorithms::
* Division Algorithms::
* Greatest Common Divisor Algorithms::
* Powering Algorithms::
* Root Extraction Algorithms::
* Radix Conversion Algorithms::
* Other Algorithms::
* Assembly Coding::
@end menu
@node Multiplication Algorithms, Division Algorithms, Algorithms, Algorithms
@section Multiplication
@cindex Multiplication algorithms
N@cross{}N limb multiplications and squares are done using one of five
algorithms, as the size N increases.
@quotation
@multitable {KaratsubaMMM} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item Algorithm @tab Threshold
@item Basecase @tab (none)
@item Karatsuba @tab @code{MUL_TOOM22_THRESHOLD}
@item Toom-3 @tab @code{MUL_TOOM33_THRESHOLD}
@item Toom-4 @tab @code{MUL_TOOM44_THRESHOLD}
@item FFT @tab @code{MUL_FFT_THRESHOLD}
@end multitable
@end quotation
Similarly for squaring, with the @code{SQR} thresholds.
N@cross{}M multiplications of operands with different sizes above
@code{MUL_TOOM22_THRESHOLD} are currently done by special Toom-inspired
algorithms or directly with FFT, depending on operand size (@pxref{Unbalanced
Multiplication}).
@menu
* Basecase Multiplication::
* Karatsuba Multiplication::
* Toom 3-Way Multiplication::
* Toom 4-Way Multiplication::
* FFT Multiplication::
* Other Multiplication::
* Unbalanced Multiplication::
@end menu
@node Basecase Multiplication, Karatsuba Multiplication, Multiplication Algorithms, Multiplication Algorithms
@subsection Basecase Multiplication
Basecase N@cross{}M multiplication is a straightforward rectangular set of
cross-products, the same as long multiplication done by hand and for that
reason sometimes known as the schoolbook or grammar school method. This is an
@m{O(NM),O(N*M)} algorithm. See Knuth section 4.3.1 algorithm M
(@pxref{References}), and the @file{mpn/generic/mul_basecase.c} code.
Assembly implementations of @code{mpn_mul_basecase} are essentially the same
as the generic C code, but have all the usual assembly tricks and
obscurities introduced for speed.
A square can be done in roughly half the time of a multiply, by using the fact
that the cross products above and below the diagonal are the same. A triangle
of products below the diagonal is formed, doubled (left shift by one bit), and
then the products on the diagonal added. This can be seen in
@file{mpn/generic/sqr_basecase.c}. Again the assembly implementations take
essentially the same approach.
@tex
defGMPline#1#2#3#4#5#6{%
hbox {%
vrule height 2.5ex depth 1ex
hbox to 2em {hfil{#2}hfil}%
vrule hbox to 2em {hfil{#3}hfil}%
vrule hbox to 2em {hfil{#4}hfil}%
vrule hbox to 2em {hfil{#5}hfil}%
vrule hbox to 2em {hfil{#6}hfil}%
vrule}}
GMPdisplay{
hbox{%
vbox{%
hbox to 1.5em {vrule height 2.5ex depth 1ex width 0pt}%
hbox {vrule height 2.5ex depth 1ex width 0pt u0hfil}%
hbox {vrule height 2.5ex depth 1ex width 0pt u1hfil}%
hbox {vrule height 2.5ex depth 1ex width 0pt u2hfil}%
hbox {vrule height 2.5ex depth 1ex width 0pt u3hfil}%
hbox {vrule height 2.5ex depth 1ex width 0pt u4hfil}%
vfill}%
vbox{%
hbox{%
hbox to 2em {hfil u0hfil}%
hbox to 2em {hfil u1hfil}%
hbox to 2em {hfil u2hfil}%
hbox to 2em {hfil u3hfil}%
hbox to 2em {hfil u4hfil}}%
vskip 0.7ex
hrule
GMPline{u0}{d}{}{}{}{}%
hrule
GMPline{u1}{}{d}{}{}{}%
hrule
GMPline{u2}{}{}{d}{}{}%
hrule
GMPline{u3}{}{}{}{d}{}%
hrule
GMPline{u4}{}{}{}{}{d}%
hrule}}}
@end tex
@ifnottex
@example
@group
u0 u1 u2 u3 u4
+---+---+---+---+---+
u0 | d | | | | |
+---+---+---+---+---+
u1 | | d | | | |
+---+---+---+---+---+
u2 | | | d | | |
+---+---+---+---+---+
u3 | | | | d | |
+---+---+---+---+---+
u4 | | | | | d |
+---+---+---+---+---+
@end group
@end example
@end ifnottex
In practice squaring isn't a full 2@cross{} faster than multiplying, it's
usually around 1.5@cross{}. Less than 1.5@cross{} probably indicates
@code{mpn_sqr_basecase} wants improving on that CPU.
On some CPUs @code{mpn_mul_basecase} can be faster than the generic C
@code{mpn_sqr_basecase} on some small sizes. @code{SQR_BASECASE_THRESHOLD} is
the size at which to use @code{mpn_sqr_basecase}, this will be zero if that
routine should be used always.
@node Karatsuba Multiplication, Toom 3-Way Multiplication, Basecase Multiplication, Multiplication Algorithms
@subsection Karatsuba Multiplication
@cindex Karatsuba multiplication
The Karatsuba multiplication algorithm is described in Knuth section 4.3.3
part A, and various other textbooks. A brief description is given here.
The inputs @math{x} and @math{y} are treated as each split into two parts of
equal length (or the most significant part one limb shorter if N is odd).
@tex
% GMPboxwidth used for all the multiplication pictures
globalnewdimenGMPboxwidth globalGMPboxwidth=5em
% GMPboxdepth and GMPboxheight are also used for the float pictures
globalnewdimenGMPboxdepth globalGMPboxdepth=1ex
globalnewdimenGMPboxheight globalGMPboxheight=2ex
gdefGMPvrule{vrule height GMPboxheight depth GMPboxdepth}
defGMPbox#1#2{%
vbox {%
hrule
hbox to 2GMPboxwidth{%
GMPvrule hfil $#1$hfil vrule hfil $#2$hfil vrule}%
hrule}}
GMPdisplay{%
vbox{%
hbox to 2GMPboxwidth {high hfil low}
vskip 0.7ex
GMPbox{x_1}{x_0}
vskip 0.5ex
GMPbox{y_1}{y_0}
}}
@end tex
@ifnottex
@example
@group
high low
+----------+----------+
| x1 | x0 |
+----------+----------+
+----------+----------+
| y1 | y0 |
+----------+----------+
@end group
@end example
@end ifnottex
Let @math{b} be the power of 2 where the split occurs, ie.@: if @ms{x,0} is
@math{k} limbs (@ms{y,0} the same) then
@m{b=2GMPraise{$k*$@code{mp_bits_per_limb}}, b=2^(k*mp_bits_per_limb)}.
With that @m{x=x_1b+x_0,x=x1*b+x0} and @m{y=y_1b+y_0,y=y1*b+y0}, and the
following holds,
@display
@m{xy = (b^2+b)x_1y_1 - b(x_1-x_0)(y_1-y_0) + (b+1)x_0y_0,
x*y = (b^2+b)*x1*y1 - b*(x1-x0)*(y1-y0) + (b+1)*x0*y0}
@end display
This formula means doing only three multiplies of (N/2)@cross{}(N/2) limbs,
whereas a basecase multiply of N@cross{}N limbs is equivalent to four
multiplies of (N/2)@cross{}(N/2). The factors @math{(b^2+b)} etc represent
the positions where the three products must be added.
@tex
defGMPboxA#1#2{%
vbox{%
hrule
hbox{%
GMPvrule
hbox to 2GMPboxwidth {hfilhbox{$#1$}hfil}%
vrule
hbox to 2GMPboxwidth {hfilhbox{$#2$}hfil}%
vrule}
hrule}}
defGMPboxB#1#2{%
hbox{%
raise GMPboxdepth hbox to GMPboxwidth {hfil #1hskip 0.5em}%
vbox{%
hrule
hbox{%
GMPvrule
hbox to 2GMPboxwidth {hfilhbox{$#2$}hfil}%
vrule}%
hrule}}}
GMPdisplay{%
vbox{%
hbox to 4GMPboxwidth {high hfil low}
vskip 0.7ex
GMPboxA{x_1y_1}{x_0y_0}
vskip 0.5ex
GMPboxB{$+$}{x_1y_1}
vskip 0.5ex
GMPboxB{$+$}{x_0y_0}
vskip 0.5ex
GMPboxB{$-$}{(x_1-x_0)(y_1-y_0)}
}}
@end tex
@ifnottex
@example
@group
high low
+--------+--------+ +--------+--------+
| x1*y1 | | x0*y0 |
+--------+--------+ +--------+--------+
+--------+--------+
add | x1*y1 |
+--------+--------+
+--------+--------+
add | x0*y0 |
+--------+--------+
+--------+--------+
sub | (x1-x0)*(y1-y0) |
+--------+--------+
@end group
@end example
@end ifnottex
The term @m{(x_1-x_0)(y_1-y_0),(x1-x0)*(y1-y0)} is best calculated as an
absolute value, and the sign used to choose to add or subtract. Notice the
sum @m{mathop{rm high}(x_0y_0)+mathop{rm low}(x_1y_1),
high(x0*y0)+low(x1*y1)} occurs twice, so it's possible to do @m{5k,5*k} limb
additions, rather than @m{6k,6*k}, but in GMP extra function call overheads
outweigh the saving.
Squaring is similar to multiplying, but with @math{x=y} the formula reduces to
an equivalent with three squares,
@display
@m{x^2 = (b^2+b)x_1^2 - b(x_1-x_0)^2 + (b+1)x_0^2,
x^2 = (b^2+b)*x1^2 - b*(x1-x0)^2 + (b+1)*x0^2}
@end display
The final result is accumulated from those three squares the same way as for
the three multiplies above. The middle term @m{(x_1-x_0)^2,(x1-x0)^2} is now
always positive.
A similar formula for both multiplying and squaring can be constructed with a
middle term @m{(x_1+x_0)(y_1+y_0),(x1+x0)*(y1+y0)}. But those sums can exceed
@math{k} limbs, leading to more carry handling and additions than the form
above.
Karatsuba multiplication is asymptotically an @math{O(N^@W{1.585})} algorithm,
the exponent being @m{log3/log2,log(3)/log(2)}, representing 3 multiplies
each @math{1/2} the size of the inputs. This is a big improvement over the
basecase multiply at @math{O(N^2)} and the advantage soon overcomes the extra
additions Karatsuba performs. @code{MUL_TOOM22_THRESHOLD} can be as little
as 10 limbs. The @code{SQR} threshold is usually about twice the @code{MUL}.
The basecase algorithm will take a time of the form @m{M(N) = aN^2 + bN + c,
M(N) = a*N^2 + b*N + c} and the Karatsuba algorithm @m{K(N) = 3M(N/2) + dN +
e, K(N) = 3*M(N/2) + d*N + e}, which expands to @m{K(N) = {3over4} aN^2 +
{3over2} bN + 3c + dN + e, K(N) = 3/4*a*N^2 + 3/2*b*N + 3*c + d*N + e}. The
factor @m{3over4, 3/4} for @math{a} means per-crossproduct speedups in the
basecase code will increase the threshold since they benefit @math{M(N)} more
than @math{K(N)}. And conversely the @m{3over2, 3/2} for @math{b} means
linear style speedups of @math{b} will increase the threshold since they
benefit @math{K(N)} more than @math{M(N)}. The latter can be seen for
instance when adding an optimized @code{mpn_sqr_diagonal} to
@code{mpn_sqr_basecase}. Of course all speedups reduce total time, and in
that sense the algorithm thresholds are merely of academic interest.
@node Toom 3-Way Multiplication, Toom 4-Way Multiplication, Karatsuba Multiplication, Multiplication Algorithms
@subsection Toom 3-Way Multiplication
@cindex Toom multiplication
The Karatsuba formula is the simplest case of a general approach to splitting
inputs that leads to both Toom and FFT algorithms. A description of
Toom can be found in Knuth section 4.3.3, with an example 3-way
calculation after Theorem A@. The 3-way form used in GMP is described here.
The operands are each considered split into 3 pieces of equal length (or the
most significant part 1 or 2 limbs shorter than the other two).
@tex
defGMPbox#1#2#3{%
vbox{%
hrule vfil
hbox to 3GMPboxwidth {%
GMPvrule
hfil$#1$hfil
vrule
hfil$#2$hfil
vrule
hfil$#3$hfil
vrule}%
vfil hrule
}}
GMPdisplay{%
vbox{%
hbox to 3GMPboxwidth {high hfil low}
vskip 0.7ex
GMPbox{x_2}{x_1}{x_0}
vskip 0.5ex
GMPbox{y_2}{y_1}{y_0}
vskip 0.5ex
}}
@end tex
@ifnottex
@example
@group
high low
+----------+----------+----------+
| x2 | x1 | x0 |
+----------+----------+----------+
+----------+----------+----------+
| y2 | y1 | y0 |
+----------+----------+----------+
@end group
@end example
@end ifnottex
@noindent
These parts are treated as the coefficients of two polynomials
@display
@group
@m{X(t) = x_2t^2 + x_1t + x_0,
X(t) = x2*t^2 + x1*t + x0}
@m{Y(t) = y_2t^2 + y_1t + y_0,
Y(t) = y2*t^2 + y1*t + y0}
@end group
@end display
Let @math{b} equal the power of 2 which is the size of the @ms{x,0}, @ms{x,1},
@ms{y,0} and @ms{y,1} pieces, ie.@: if they're @math{k} limbs each then
@m{b=2GMPraise{$k*$@code{mp_bits_per_limb}}, b=2^(k*mp_bits_per_limb)}.
With this @math{x=X(b)} and @math{y=Y(b)}.
Let a polynomial @m{W(t)=X(t)Y(t),W(t)=X(t)*Y(t)} and suppose its coefficients
are
@display
@m{W(t) = w_4t^4 + w_3t^3 + w_2t^2 + w_1t + w_0,
W(t) = w4*t^4 + w3*t^3 + w2*t^2 + w1*t + w0}
@end display
The @m{w_i,w[i]} are going to be determined, and when they are they'll give
the final result using @math{w=W(b)}, since
@m{xy=X(b)Y(b),x*y=X(b)*Y(b)=W(b)}. The coefficients will be roughly
@math{b^2} each, and the final @math{W(b)} will be an addition like,
@tex
defGMPbox#1#2{%
moveright #1GMPboxwidth
vbox{%
hrule
hbox{%
GMPvrule
hbox to 2GMPboxwidth {hfil$#2$hfil}%
vrule}%
hrule
}}
GMPdisplay{%
vbox{%
hbox to 6GMPboxwidth {high hfil low}%
vskip 0.7ex
GMPbox{0}{w_4}
vskip 0.5ex
GMPbox{1}{w_3}
vskip 0.5ex
GMPbox{2}{w_2}
vskip 0.5ex
GMPbox{3}{w_1}
vskip 0.5ex
GMPbox{4}{w_0}
}}
@end tex
@ifnottex
@example
@group
high low
+-------+-------+
| w4 |
+-------+-------+
+--------+-------+
| w3 |
+--------+-------+
+--------+-------+
| w2 |
+--------+-------+
+--------+-------+
| w1 |
+--------+-------+
+-------+-------+
| w0 |
+-------+-------+
@end group
@end example
@end ifnottex
The @m{w_i,w[i]} coefficients could be formed by a simple set of cross
products, like @m{w_4=x_2y_2,w4=x2*y2}, @m{w_3=x_2y_1+x_1y_2,w3=x2*y1+x1*y2},
@m{w_2=x_2y_0+x_1y_1+x_0y_2,w2=x2*y0+x1*y1+x0*y2} etc, but this would need all
nine @m{x_iy_j,x[i]*y[j]} for @math{i,j=0,1,2}, and would be equivalent merely
to a basecase multiply. Instead the following approach is used.
@math{X(t)} and @math{Y(t)} are evaluated and multiplied at 5 points, giving
values of @math{W(t)} at those points. In GMP the following points are used,
@quotation
@multitable {@m{t=infty,t=inf}M} {MMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMMM}
@item Point @tab Value
@item @math{t=0} @tab @m{x_0y_0,x0 * y0}, which gives @ms{w,0} immediately
@item @math{t=1} @tab @m{(x_2+x_1+x_0)(y_2+y_1+y_0),(x2+x1+x0) * (y2+y1+y0)}
@item @math{t=-1} @tab @m{(x_2-x_1+x_0)(y_2-y_1+y_0),(x2-x1+x0) * (y2-y1+y0)}
@item @math{t=2} @tab @m{(4x_2+2x_1+x_0)(4y_2+2y_1+y_0),(4*x2+2*x1+x0) * (4*y2+2*y1+y0)}
@item @m{t=infty,t=inf} @tab @m{x_2y_2,x2 * y2}, which gives @ms{w,4} immediately
@end multitable
@end quotation
At @math{t=-1} the values can be negative and that's handled using the
absolute values and tracking the sign separately. At @m{t=infty,t=inf} the
value is actually @m{lim_{ttoinfty} {X(t)Y(t)over t^4}, X(t)*Y(t)/t^4 in
the limit as t approaches infinity}, but it's much easier to think of as
simply @m{x_2y_2,x2*y2} giving @ms{w,4} immediately (much like
@m{x_0y_0,x0*y0} at @math{t=0} gives @ms{w,0} immediately).
Each of the points substituted into
@m{W(t)=w_4t^4+cdots+w_0,W(t)=w4*t^4+@dots{}+w0} gives a linear combination
of the @m{w_i,w[i]} coefficients, and the value of those combinations has just
been calculated.
@tex
GMPdisplay{%
$matrix{%
W(0) & = & & & & & & & & & w_0 cr
W(1) & = & w_4 & + & w_3 & + & w_2 & + & w_1 & + & w_0 cr
W(-1) & = & w_4 & - & w_3 & + & w_2 & - & w_1 & + & w_0 cr
W(2) & = & 16w_4 & + & 8w_3 & + & 4w_2 & + & 2w_1 & + & w_0 cr
W(infty) & = & w_4 cr
}$}
@end tex
@ifnottex
@example
@group
W(0) = w0
W(1) = w4 + w3 + w2 + w1 + w0
W(-1) = w4 - w3 + w2 - w1 + w0
W(2) = 16*w4 + 8*w3 + 4*w2 + 2*w1 + w0
W(inf) = w4
@end group
@end example
@end ifnottex
This is a set of five equations in five unknowns, and some elementary linear
algebra quickly isolates each @m{w_i,w[i]}. This involves adding or
subtracting one @math{W(t)} value from another, and a couple of divisions by
powers of 2 and one division by 3, the latter using the special
@code{mpn_divexact_by3} (@pxref{Exact Division}).
The conversion of @math{W(t)} values to the coefficients is interpolation. A
polynomial of degree 4 like @math{W(t)} is uniquely determined by values known
at 5 different points. The points are arbitrary and can be chosen to make the
linear equations come out with a convenient set of steps for quickly isolating
the @m{w_i,w[i]}.
Squaring follows the same procedure as multiplication, but there's only one
@math{X(t)} and it's evaluated at the 5 points, and those values squared to
give values of @math{W(t)}. The interpolation is then identical, and in fact
the same @code{toom3_interpolate} subroutine is used for both squaring and
multiplying.
Toom-3 is asymptotically @math{O(N^@W{1.465})}, the exponent being
@m{log5/log3,log(5)/log(3)}, representing 5 recursive multiplies of 1/3 the
original size each. This is an improvement over Karatsuba at
@math{O(N^@W{1.585})}, though Toom does more work in the evaluation and
interpolation and so it only realizes its advantage above a certain size.
Near the crossover between Toom-3 and Karatsuba there's generally a range of
sizes where the difference between the two is small.
@code{MUL_TOOM33_THRESHOLD} is a somewhat arbitrary point in that range and
successive runs of the tune program can give different values due to small
variations in measuring. A graph of time versus size for the two shows the
effect, see @file{tune/README}.
At the fairly small sizes where the Toom-3 thresholds occur it's worth
remembering that the asymptotic behaviour for Karatsuba and Toom-3 can't be
expected to make accurate predictions, due of course to the big influence of
all sorts of overheads, and the fact that only a few recursions of each are
being performed. Even at large sizes there's a good chance machine dependent
effects like cache architecture will mean actual performance deviates from
what might be predicted.
The formula given for the Karatsuba algorithm (@pxref{Karatsuba
Multiplication}) has an equivalent for Toom-3 involving only five multiplies,
but this would be complicated and unenlightening.
An alternate view of Toom-3 can be found in Zuras (@pxref{References}), using
a vector to represent the @math{x} and @math{y} splits and a matrix
multiplication for the evaluation and interpolation stages. The matrix
inverses are not meant to be actually used, and they have elements with values
much greater than in fact arise in the interpolation steps. The diagram shown
for the 3-way is attractive, but again doesn't have to be implemented that way
and for example with a bit of rearrangement just one division by 6 can be
done.
@node Toom 4-Way Multiplication, FFT Multiplication, Toom 3-Way Multiplication, Multiplication Algorithms
@subsection Toom 4-Way Multiplication
@cindex Toom multiplication
Karatsuba and Toom-3 split the operands into 2 and 3 coefficients,
respectively. Toom-4 analogously splits the operands into 4 coefficients.
Using the notation from the section on Toom-3 multiplication, we form two
polynomials:
@display
@group
@m{X(t) = x_3t^3 + x_2t^2 + x_1t + x_0,
X(t) = x3*t^3 + x2*t^2 + x1*t + x0}
@m{Y(t) = y_3t^3 + y_2t^2 + y_1t + y_0,
Y(t) = y3*t^3 + y2*t^2 + y1*t + y0}
@end group
@end display