1 This is ../../gmp/doc/gmp.info, produced by makeinfo version 4.8 from
2 ../../gmp/doc/gmp.texi.
4 This manual describes how to install and use the GNU multiple
5 precision arithmetic library, version 5.0.1.
7 Copyright 1991, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000,
8 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 Free
9 Software Foundation, Inc.
11 Permission is granted to copy, distribute and/or modify this
12 document under the terms of the GNU Free Documentation License, Version
13 1.3 or any later version published by the Free Software Foundation;
14 with no Invariant Sections, with the Front-Cover Texts being "A GNU
15 Manual", and with the Back-Cover Texts being "You have freedom to copy
16 and modify this GNU Manual, like GNU software". A copy of the license
17 is included in *Note GNU Free Documentation License::.
19 INFO-DIR-SECTION GNU libraries
21 * gmp: (gmp). GNU Multiple Precision Arithmetic Library.
25 File: gmp.info, Node: Top, Next: Copying, Prev: (dir), Up: (dir)
30 This manual describes how to install and use the GNU multiple
31 precision arithmetic library, version 5.0.1.
33 Copyright 1991, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000,
34 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010 Free
35 Software Foundation, Inc.
37 Permission is granted to copy, distribute and/or modify this
38 document under the terms of the GNU Free Documentation License, Version
39 1.3 or any later version published by the Free Software Foundation;
40 with no Invariant Sections, with the Front-Cover Texts being "A GNU
41 Manual", and with the Back-Cover Texts being "You have freedom to copy
42 and modify this GNU Manual, like GNU software". A copy of the license
43 is included in *Note GNU Free Documentation License::.
48 * Copying:: GMP Copying Conditions (LGPL).
49 * Introduction to GMP:: Brief introduction to GNU MP.
50 * Installing GMP:: How to configure and compile the GMP library.
51 * GMP Basics:: What every GMP user should know.
52 * Reporting Bugs:: How to usefully report bugs.
53 * Integer Functions:: Functions for arithmetic on signed integers.
54 * Rational Number Functions:: Functions for arithmetic on rational numbers.
55 * Floating-point Functions:: Functions for arithmetic on floats.
56 * Low-level Functions:: Fast functions for natural numbers.
57 * Random Number Functions:: Functions for generating random numbers.
58 * Formatted Output:: `printf' style output.
59 * Formatted Input:: `scanf' style input.
60 * C++ Class Interface:: Class wrappers around GMP types.
61 * BSD Compatible Functions:: All functions found in BSD MP.
62 * Custom Allocation:: How to customize the internal allocation.
63 * Language Bindings:: Using GMP from other languages.
64 * Algorithms:: What happens behind the scenes.
65 * Internals:: How values are represented behind the scenes.
67 * Contributors:: Who brings you this library?
68 * References:: Some useful papers and books to read.
69 * GNU Free Documentation License::
74 File: gmp.info, Node: Copying, Next: Introduction to GMP, Prev: Top, Up: Top
76 GNU MP Copying Conditions
77 *************************
79 This library is "free"; this means that everyone is free to use it and
80 free to redistribute it on a free basis. The library is not in the
81 public domain; it is copyrighted and there are restrictions on its
82 distribution, but these restrictions are designed to permit everything
83 that a good cooperating citizen would want to do. What is not allowed
84 is to try to prevent others from further sharing any version of this
85 library that they might get from you.
87 Specifically, we want to make sure that you have the right to give
88 away copies of the library, that you receive source code or else can
89 get it if you want it, that you can change this library or use pieces
90 of it in new free programs, and that you know you can do these things.
92 To make sure that everyone has such rights, we have to forbid you to
93 deprive anyone else of these rights. For example, if you distribute
94 copies of the GNU MP library, you must give the recipients all the
95 rights that you have. You must make sure that they, too, receive or
96 can get the source code. And you must tell them their rights.
98 Also, for our own protection, we must make certain that everyone
99 finds out that there is no warranty for the GNU MP library. If it is
100 modified by someone else and passed on, we want their recipients to
101 know that what they have is not what we distributed, so that any
102 problems introduced by others will not reflect on our reputation.
104 The precise conditions of the license for the GNU MP library are
105 found in the Lesser General Public License version 3 that accompanies
106 the source code, see `COPYING.LIB'. Certain demonstration programs are
107 provided under the terms of the plain General Public License version 3,
111 File: gmp.info, Node: Introduction to GMP, Next: Installing GMP, Prev: Copying, Up: Top
113 1 Introduction to GNU MP
114 ************************
116 GNU MP is a portable library written in C for arbitrary precision
117 arithmetic on integers, rational numbers, and floating-point numbers.
118 It aims to provide the fastest possible arithmetic for all applications
119 that need higher precision than is directly supported by the basic C
122 Many applications use just a few hundred bits of precision; but some
123 applications may need thousands or even millions of bits. GMP is
124 designed to give good performance for both, by choosing algorithms
125 based on the sizes of the operands, and by carefully keeping the
126 overhead at a minimum.
128 The speed of GMP is achieved by using fullwords as the basic
129 arithmetic type, by using sophisticated algorithms, by including
130 carefully optimized assembly code for the most common inner loops for
131 many different CPUs, and by a general emphasis on speed (as opposed to
132 simplicity or elegance).
134 There is assembly code for these CPUs: ARM, DEC Alpha 21064, 21164,
135 and 21264, AMD 29000, AMD K6, K6-2, Athlon, and Athlon64, Hitachi
136 SuperH and SH-2, HPPA 1.0, 1.1 and 2.0, Intel Pentium, Pentium
137 Pro/II/III, Pentium 4, generic x86, Intel IA-64, i960, Motorola
138 MC68000, MC68020, MC88100, and MC88110, Motorola/IBM PowerPC 32 and 64,
139 National NS32000, IBM POWER, MIPS R3000, R4000, SPARCv7, SuperSPARC,
140 generic SPARCv8, UltraSPARC, DEC VAX, and Zilog Z8000. Some
141 optimizations also for Cray vector systems, Clipper, IBM ROMP (RT), and
144 For up-to-date information on GMP, please see the GMP web pages at
148 The latest version of the library is available at
150 `ftp://ftp.gnu.org/gnu/gmp/'
152 Many sites around the world mirror `ftp.gnu.org', please use a mirror
153 near you, see `http://www.gnu.org/order/ftp.html' for a full list.
155 There are three public mailing lists of interest. One for release
156 announcements, one for general questions and discussions about usage of
157 the GMP library and one for bug reports. For more information, see
159 `http://gmplib.org/mailman/listinfo/'.
161 The proper place for bug reports is <gmp-bugs@gmplib.org>. See
162 *Note Reporting Bugs:: for information about reporting bugs.
165 1.1 How to use this Manual
166 ==========================
168 Everyone should read *Note GMP Basics::. If you need to install the
169 library yourself, then read *Note Installing GMP::. If you have a
170 system with multiple ABIs, then read *Note ABI and ISA::, for the
171 compiler options that must be used on applications.
173 The rest of the manual can be used for later reference, although it
174 is probably a good idea to glance through it.
177 File: gmp.info, Node: Installing GMP, Next: GMP Basics, Prev: Introduction to GMP, Up: Top
182 GMP has an autoconf/automake/libtool based configuration system. On a
183 Unix-like system a basic build can be done with
188 Some self-tests can be run with
192 And you can install (under `/usr/local' by default) with
196 If you experience problems, please report them to
197 <gmp-bugs@gmplib.org>. See *Note Reporting Bugs::, for information on
198 what to include in useful bug reports.
204 * Notes for Package Builds::
205 * Notes for Particular Systems::
206 * Known Build Problems::
207 * Performance optimization::
210 File: gmp.info, Node: Build Options, Next: ABI and ISA, Prev: Installing GMP, Up: Installing GMP
215 All the usual autoconf configure options are available, run `./configure
216 --help' for a summary. The file `INSTALL.autoconf' has some generic
217 installation information too.
220 `configure' requires various Unix-like tools. See *Note Notes for
221 Particular Systems::, for some options on non-Unix systems.
223 It might be possible to build without the help of `configure',
224 certainly all the code is there, but unfortunately you'll be on
228 To compile in a separate build directory, `cd' to that directory,
229 and prefix the configure command with the path to the GMP source
230 directory. For example
233 /my/sources/gmp-5.0.1/configure
235 Not all `make' programs have the necessary features (`VPATH') to
236 support this. In particular, SunOS and Slowaris `make' have bugs
237 that make them unable to build in a separate directory. Use GNU
240 `--prefix' and `--exec-prefix'
241 The `--prefix' option can be used in the normal way to direct GMP
242 to install under a particular tree. The default is `/usr/local'.
244 `--exec-prefix' can be used to direct architecture-dependent files
245 like `libgmp.a' to a different location. This can be used to share
246 architecture-independent parts like the documentation, but
247 separate the dependent parts. Note however that `gmp.h' and
248 `mp.h' are architecture-dependent since they encode certain
249 aspects of `libgmp', so it will be necessary to ensure both
250 `$prefix/include' and `$exec_prefix/include' are available to the
253 `--disable-shared', `--disable-static'
254 By default both shared and static libraries are built (where
255 possible), but one or other can be disabled. Shared libraries
256 result in smaller executables and permit code sharing between
257 separate running processes, but on some CPUs are slightly slower,
258 having a small cost on each function call.
260 Native Compilation, `--build=CPU-VENDOR-OS'
261 For normal native compilation, the system can be specified with
262 `--build'. By default `./configure' uses the output from running
263 `./config.guess'. On some systems `./config.guess' can determine
264 the exact CPU type, on others it will be necessary to give it
265 explicitly. For example,
267 ./configure --build=ultrasparc-sun-solaris2.7
269 In all cases the `OS' part is important, since it controls how
270 libtool generates shared libraries. Running `./config.guess' is
271 the simplest way to see what it should be, if you don't know
274 Cross Compilation, `--host=CPU-VENDOR-OS'
275 When cross-compiling, the system used for compiling is given by
276 `--build' and the system where the library will run is given by
277 `--host'. For example when using a FreeBSD Athlon system to build
278 GNU/Linux m68k binaries,
280 ./configure --build=athlon-pc-freebsd3.5 --host=m68k-mac-linux-gnu
282 Compiler tools are sought first with the host system type as a
283 prefix. For example `m68k-mac-linux-gnu-ranlib' is tried, then
284 plain `ranlib'. This makes it possible for a set of
285 cross-compiling tools to co-exist with native tools. The prefix
286 is the argument to `--host', and this can be an alias, such as
287 `m68k-linux'. But note that tools don't have to be setup this
288 way, it's enough to just have a `PATH' with a suitable
289 cross-compiling `cc' etc.
291 Compiling for a different CPU in the same family as the build
292 system is a form of cross-compilation, though very possibly this
293 would merely be special options on a native compiler. In any case
294 `./configure' avoids depending on being able to run code on the
295 build system, which is important when creating binaries for a
296 newer CPU since they very possibly won't run on the build system.
298 In all cases the compiler must be able to produce an executable
299 (of whatever format) from a standard C `main'. Although only
300 object files will go to make up `libgmp', `./configure' uses
301 linking tests for various purposes, such as determining what
302 functions are available on the host system.
304 Currently a warning is given unless an explicit `--build' is used
305 when cross-compiling, because it may not be possible to correctly
306 guess the build system type if the `PATH' has only a
307 cross-compiling `cc'.
309 Note that the `--target' option is not appropriate for GMP. It's
310 for use when building compiler tools, with `--host' being where
311 they will run, and `--target' what they'll produce code for.
312 Ordinary programs or libraries like GMP are only interested in the
313 `--host' part, being where they'll run. (Some past versions of
314 GMP used `--target' incorrectly.)
317 In general, if you want a library that runs as fast as possible,
318 you should configure GMP for the exact CPU type your system uses.
319 However, this may mean the binaries won't run on older members of
320 the family, and might run slower on other members, older or newer.
321 The best idea is always to build GMP for the exact machine type
322 you intend to run it on.
324 The following CPUs have specific support. See `configure.in' for
325 details of what code and compiler options they select.
327 * Alpha: alpha, alphaev5, alphaev56, alphapca56, alphapca57,
328 alphaev6, alphaev67, alphaev68 alphaev7
330 * Cray: c90, j90, t90, sv1
332 * HPPA: hppa1.0, hppa1.1, hppa2.0, hppa2.0n, hppa2.0w, hppa64
334 * IA-64: ia64, itanium, itanium2
336 * MIPS: mips, mips3, mips64
338 * Motorola: m68k, m68000, m68010, m68020, m68030, m68040,
339 m68060, m68302, m68360, m88k, m88110
341 * POWER: power, power1, power2, power2sc
343 * PowerPC: powerpc, powerpc64, powerpc401, powerpc403,
344 powerpc405, powerpc505, powerpc601, powerpc602, powerpc603,
345 powerpc603e, powerpc604, powerpc604e, powerpc620, powerpc630,
346 powerpc740, powerpc7400, powerpc7450, powerpc750, powerpc801,
347 powerpc821, powerpc823, powerpc860, powerpc970
349 * SPARC: sparc, sparcv8, microsparc, supersparc, sparcv9,
350 ultrasparc, ultrasparc2, ultrasparc2i, ultrasparc3, sparc64
352 * x86 family: i386, i486, i586, pentium, pentiummmx, pentiumpro,
353 pentium2, pentium3, pentium4, k6, k62, k63, athlon, amd64,
356 * Other: a29k, arm, clipper, i960, ns32k, pyramid, sh, sh2, vax,
359 CPUs not listed will use generic C code.
362 If some of the assembly code causes problems, or if otherwise
363 desired, the generic C code can be selected with CPU `none'. For
366 ./configure --host=none-unknown-freebsd3.5
368 Note that this will run quite slowly, but it should be portable
369 and should at least make it possible to get something running if
372 Fat binary, `--enable-fat'
373 Using `--enable-fat' selects a "fat binary" build on x86, where
374 optimized low level subroutines are chosen at runtime according to
375 the CPU detected. This means more code, but gives good
376 performance on all x86 chips. (This option might become available
377 for more architectures in the future.)
380 On some systems GMP supports multiple ABIs (application binary
381 interfaces), meaning data type sizes and calling conventions. By
382 default GMP chooses the best ABI available, but a particular ABI
383 can be selected. For example
385 ./configure --host=mips64-sgi-irix6 ABI=n32
387 See *Note ABI and ISA::, for the available choices on relevant
388 CPUs, and what applications need to do.
391 By default the C compiler used is chosen from among some likely
392 candidates, with `gcc' normally preferred if it's present. The
393 usual `CC=whatever' can be passed to `./configure' to choose
396 For various systems, default compiler flags are set based on the
397 CPU and compiler. The usual `CFLAGS="-whatever"' can be passed to
398 `./configure' to use something different or to set good flags for
399 systems GMP doesn't otherwise know.
401 The `CC' and `CFLAGS' used are printed during `./configure', and
402 can be found in each generated `Makefile'. This is the easiest way
403 to check the defaults when considering changing or adding
406 Note that when `CC' and `CFLAGS' are specified on a system
407 supporting multiple ABIs it's important to give an explicit
408 `ABI=whatever', since GMP can't determine the ABI just from the
409 flags and won't be able to select the correct assembly code.
411 If just `CC' is selected then normal default `CFLAGS' for that
412 compiler will be used (if GMP recognises it). For example
413 `CC=gcc' can be used to force the use of GCC, with default flags
417 Any flags like `-D' defines or `-I' includes required by the
418 preprocessor should be set in `CPPFLAGS' rather than `CFLAGS'.
419 Compiling is done with both `CPPFLAGS' and `CFLAGS', but
420 preprocessing uses just `CPPFLAGS'. This distinction is because
421 most preprocessors won't accept all the flags the compiler does.
422 Preprocessing is done separately in some configure tests, and in
423 the `ansi2knr' support for K&R compilers.
426 Some build-time programs are compiled and run to generate
427 host-specific data tables. `CC_FOR_BUILD' is the compiler used
428 for this. It doesn't need to be in any particular ABI or mode, it
429 merely needs to generate executables that can run. The default is
430 to try the selected `CC' and some likely candidates such as `cc'
431 and `gcc', looking for something that works.
433 No flags are used with `CC_FOR_BUILD' because a simple invocation
434 like `cc foo.c' should be enough. If some particular options are
435 required they can be included as for instance `CC_FOR_BUILD="cc
438 C++ Support, `--enable-cxx'
439 C++ support in GMP can be enabled with `--enable-cxx', in which
440 case a C++ compiler will be required. As a convenience
441 `--enable-cxx=detect' can be used to enable C++ support only if a
442 compiler can be found. The C++ support consists of a library
443 `libgmpxx.la' and header file `gmpxx.h' (*note Headers and
446 A separate `libgmpxx.la' has been adopted rather than having C++
447 objects within `libgmp.la' in order to ensure dynamic linked C
448 programs aren't bloated by a dependency on the C++ standard
449 library, and to avoid any chance that the C++ compiler could be
450 required when linking plain C programs.
452 `libgmpxx.la' will use certain internals from `libgmp.la' and can
453 only be expected to work with `libgmp.la' from the same GMP
454 version. Future changes to the relevant internals will be
455 accompanied by renaming, so a mismatch will cause unresolved
456 symbols rather than perhaps mysterious misbehaviour.
458 In general `libgmpxx.la' will be usable only with the C++ compiler
459 that built it, since name mangling and runtime support are usually
460 incompatible between different compilers.
463 When C++ support is enabled, the C++ compiler and its flags can be
464 set with variables `CXX' and `CXXFLAGS' in the usual way. The
465 default for `CXX' is the first compiler that works from a list of
466 likely candidates, with `g++' normally preferred when available.
467 The default for `CXXFLAGS' is to try `CFLAGS', `CFLAGS' without
468 `-g', then for `g++' either `-g -O2' or `-O2', or for other
469 compilers `-g' or nothing. Trying `CFLAGS' this way is convenient
470 when using `gcc' and `g++' together, since the flags for `gcc' will
473 It's important that the C and C++ compilers match, meaning their
474 startup and runtime support routines are compatible and that they
475 generate code in the same ABI (if there's a choice of ABIs on the
476 system). `./configure' isn't currently able to check these things
477 very well itself, so for that reason `--disable-cxx' is the
478 default, to avoid a build failure due to a compiler mismatch.
479 Perhaps this will change in the future.
481 Incidentally, it's normally not good enough to set `CXX' to the
482 same as `CC'. Although `gcc' for instance recognises `foo.cc' as
483 C++ code, only `g++' will invoke the linker the right way when
484 building an executable or shared library from C++ object files.
486 Temporary Memory, `--enable-alloca=<choice>'
487 GMP allocates temporary workspace using one of the following three
488 methods, which can be selected with for instance
489 `--enable-alloca=malloc-reentrant'.
491 * `alloca' - C library or compiler builtin.
493 * `malloc-reentrant' - the heap, in a re-entrant fashion.
495 * `malloc-notreentrant' - the heap, with global variables.
497 For convenience, the following choices are also available.
498 `--disable-alloca' is the same as `no'.
500 * `yes' - a synonym for `alloca'.
502 * `no' - a synonym for `malloc-reentrant'.
504 * `reentrant' - `alloca' if available, otherwise
505 `malloc-reentrant'. This is the default.
507 * `notreentrant' - `alloca' if available, otherwise
508 `malloc-notreentrant'.
510 `alloca' is reentrant and fast, and is recommended. It actually
511 allocates just small blocks on the stack; larger ones use
514 `malloc-reentrant' is, as the name suggests, reentrant and thread
515 safe, but `malloc-notreentrant' is faster and should be used if
516 reentrancy is not required.
518 The two malloc methods in fact use the memory allocation functions
519 selected by `mp_set_memory_functions', these being `malloc' and
520 friends by default. *Note Custom Allocation::.
522 An additional choice `--enable-alloca=debug' is available, to help
523 when debugging memory related problems (*note Debugging::).
525 FFT Multiplication, `--disable-fft'
526 By default multiplications are done using Karatsuba, 3-way Toom,
527 and Fermat FFT. The FFT is only used on large to very large
528 operands and can be disabled to save code size if desired.
530 Berkeley MP, `--enable-mpbsd'
531 The Berkeley MP compatibility library (`libmp') and header file
532 (`mp.h') are built and installed only if `--enable-mpbsd' is used.
533 *Note BSD Compatible Functions::.
535 Assertion Checking, `--enable-assert'
536 This option enables some consistency checking within the library.
537 This can be of use while debugging, *note Debugging::.
539 Execution Profiling, `--enable-profiling=prof/gprof/instrument'
540 Enable profiling support, in one of various styles, *note
544 Various assembly versions of each mpn subroutines are provided.
545 For a given CPU, a search is made though a path to choose a
546 version of each. For example `sparcv8' has
548 MPN_PATH="sparc32/v8 sparc32 generic"
550 which means look first for v8 code, then plain sparc32 (which is
551 v7), and finally fall back on generic C. Knowledgeable users with
552 special requirements can specify a different path. Normally this
553 is completely unnecessary.
556 The source for the document you're now reading is `doc/gmp.texi',
557 in Texinfo format, see *Note Texinfo: (texinfo)Top.
559 Info format `doc/gmp.info' is included in the distribution. The
560 usual automake targets are available to make PostScript, DVI, PDF
561 and HTML (these will require various TeX and Texinfo tools).
563 DocBook and XML can be generated by the Texinfo `makeinfo' program
564 too, see *Note Options for `makeinfo': (texinfo)makeinfo options.
566 Some supplementary notes can also be found in the `doc'
571 File: gmp.info, Node: ABI and ISA, Next: Notes for Package Builds, Prev: Build Options, Up: Installing GMP
576 ABI (Application Binary Interface) refers to the calling conventions
577 between functions, meaning what registers are used and what sizes the
578 various C data types are. ISA (Instruction Set Architecture) refers to
579 the instructions and registers a CPU has available.
581 Some 64-bit ISA CPUs have both a 64-bit ABI and a 32-bit ABI
582 defined, the latter for compatibility with older CPUs in the family.
583 GMP supports some CPUs like this in both ABIs. In fact within GMP
584 `ABI' means a combination of chip ABI, plus how GMP chooses to use it.
585 For example in some 32-bit ABIs, GMP may support a limb as either a
586 32-bit `long' or a 64-bit `long long'.
588 By default GMP chooses the best ABI available for a given system,
589 and this generally gives significantly greater speed. But an ABI can
590 be chosen explicitly to make GMP compatible with other libraries, or
591 particular application requirements. For example,
595 In all cases it's vital that all object code used in a given program
596 is compiled for the same ABI.
598 Usually a limb is implemented as a `long'. When a `long long' limb
599 is used this is encoded in the generated `gmp.h'. This is convenient
600 for applications, but it does mean that `gmp.h' will vary, and can't be
601 just copied around. `gmp.h' remains compiler independent though, since
602 all compilers for a particular ABI will be expected to use the same
605 Currently no attempt is made to follow whatever conventions a system
606 has for installing library or header files built for a particular ABI.
607 This will probably only matter when installing multiple builds of GMP,
608 and it might be as simple as configuring with a special `libdir', or it
609 might require more than that. Note that builds for different ABIs need
610 to done separately, with a fresh `./configure' and `make' each.
614 On AMD64 systems supporting both 32-bit and 64-bit modes for
615 applications, the following ABI choices are available.
618 The 64-bit ABI uses 64-bit limbs and pointers and makes full
619 use of the chip architecture. This is the default.
620 Applications will usually not need special compiler flags,
621 but for reference the option is
626 The 32-bit ABI is the usual i386 conventions. This will be
627 slower, and is not recommended except for inter-operating
628 with other code not yet 64-bit capable. Applications must be
633 (In GCC 2.95 and earlier there's no `-m32' option, it's the
637 HPPA 2.0 (`hppa2.0*', `hppa64')
640 The 2.0w ABI uses 64-bit limbs and pointers and is available
641 on HP-UX 11 or up. Applications must be compiled with
647 The 2.0n ABI means the 32-bit HPPA 1.0 ABI and all its normal
648 calling conventions, but with 64-bit instructions permitted
649 within functions. GMP uses a 64-bit `long long' for a limb.
650 This ABI is available on hppa64 GNU/Linux and on HP-UX 10 or
651 higher. Applications must be compiled with
656 Note that current versions of GCC (eg. 3.2) don't generate
657 64-bit instructions for `long long' operations and so may be
658 slower than for 2.0w. (The GMP assembly code is the same
662 HPPA 2.0 CPUs can run all HPPA 1.0 and 1.1 code in the 32-bit
663 HPPA 1.0 ABI. No special compiler options are needed for
666 All three ABIs are available for CPU types `hppa2.0w', `hppa2.0'
667 and `hppa64', but for CPU type `hppa2.0n' only 2.0n or 1.0 are
670 Note that GCC on HP-UX has no options to choose between 2.0n and
671 2.0w modes, unlike HP `cc'. Instead it must be built for one or
672 the other ABI. GMP will detect how it was built, and skip to the
676 IA-64 under HP-UX (`ia64*-*-hpux*', `itanium*-*-hpux*')
677 HP-UX supports two ABIs for IA-64. GMP performance is the same in
681 In the 32-bit ABI, pointers, `int's and `long's are 32 bits
682 and GMP uses a 64 bit `long long' for a limb. Applications
683 can be compiled without any special flags since this ABI is
684 the default in both HP C and GCC, but for reference the flags
691 In the 64-bit ABI, `long's and pointers are 64 bits and GMP
692 uses a `long' for a limb. Applications must be compiled with
697 On other IA-64 systems, GNU/Linux for instance, `ABI=64' is the
701 MIPS under IRIX 6 (`mips*-*-irix[6789]')
702 IRIX 6 always has a 64-bit MIPS 3 or better CPU, and supports ABIs
703 o32, n32, and 64. n32 or 64 are recommended, and GMP performance
704 will be the same in each. The default is n32.
707 The o32 ABI is 32-bit pointers and integers, and no 64-bit
708 operations. GMP will be slower than in n32 or 64, this
709 option only exists to support old compilers, eg. GCC 2.7.2.
710 Applications can be compiled with no special flags on an old
711 compiler, or on a newer compiler with
717 The n32 ABI is 32-bit pointers and integers, but with a
718 64-bit limb using a `long long'. Applications must be
725 The 64-bit ABI is 64-bit pointers and integers. Applications
726 must be compiled with
731 Note that MIPS GNU/Linux, as of kernel version 2.2, doesn't have
732 the necessary support for n32 or 64 and so only gets a 32-bit limb
736 PowerPC 64 (`powerpc64', `powerpc620', `powerpc630', `powerpc970', `power4', `power5')
739 The AIX 64 ABI uses 64-bit limbs and pointers and is the
740 default on PowerPC 64 `*-*-aix*' systems. Applications must
747 The `mode64' ABI uses 64-bit limbs and pointers, and is the
748 default on 64-bit GNU/Linux, BSD, and Mac OS X/Darwin
749 systems. Applications must be compiled with
754 The `mode32' ABI uses a 64-bit `long long' limb but with the
755 chip still in 32-bit mode and using 32-bit calling
756 conventions. This is the default on for systems where the
757 true 64-bit ABIs are unavailable. No special compiler
758 options are needed for applications.
761 This is the basic 32-bit PowerPC ABI, with a 32-bit limb. No
762 special compiler options are needed for applications.
764 GMP speed is greatest in `aix64' and `mode32'. In `ABI=32' only
765 the 32-bit ISA is used and this doesn't make full use of a 64-bit
766 chip. On a suitable system we could perhaps use more of the ISA,
767 but there are no plans to do so.
770 Sparc V9 (`sparc64', `sparcv9', `ultrasparc*')
773 The 64-bit V9 ABI is available on the various BSD sparc64
774 ports, recent versions of Sparc64 GNU/Linux, and Solaris 2.7
775 and up (when the kernel is in 64-bit mode). GCC 3.2 or
776 higher, or Sun `cc' is required. On GNU/Linux, depending on
777 the default `gcc' mode, applications must be compiled with
781 On Solaris applications must be compiled with
783 gcc -m64 -mptr64 -Wa,-xarch=v9 -mcpu=v9
786 On the BSD sparc64 systems no special options are required,
787 since 64-bits is the only ABI available.
790 For the basic 32-bit ABI, GMP still uses as much of the V9
791 ISA as it can. In the Sun documentation this combination is
792 known as "v8plus". On GNU/Linux, depending on the default
793 `gcc' mode, applications may need to be compiled with
797 On Solaris, no special compiler options are required for
798 applications, though using something like the following is
799 recommended. (`gcc' 2.8 and earlier only support `-mv8'
805 GMP speed is greatest in `ABI=64', so it's the default where
806 available. The speed is partly because there are extra registers
807 available and partly because 64-bits is considered the more
808 important case and has therefore had better code written for it.
810 Don't be confused by the names of the `-m' and `-x' compiler
811 options, they're called `arch' but effectively control both ABI
814 On Solaris 2.6 and earlier, only `ABI=32' is available since the
815 kernel doesn't save all registers.
817 On Solaris 2.7 with the kernel in 32-bit mode, a normal native
818 build will reject `ABI=64' because the resulting executables won't
819 run. `ABI=64' can still be built if desired by making it look
820 like a cross-compile, for example
822 ./configure --build=none --host=sparcv9-sun-solaris2.7 ABI=64
825 File: gmp.info, Node: Notes for Package Builds, Next: Notes for Particular Systems, Prev: ABI and ISA, Up: Installing GMP
827 2.3 Notes for Package Builds
828 ============================
830 GMP should present no great difficulties for packaging in a binary
833 Libtool is used to build the library and `-version-info' is set
834 appropriately, having started from `3:0:0' in GMP 3.0 (*note Library
835 interface versions: (libtool)Versioning.).
837 The GMP 4 series will be upwardly binary compatible in each release
838 and will be upwardly binary compatible with all of the GMP 3 series.
839 Additional function interfaces may be added in each release, so on
840 systems where libtool versioning is not fully checked by the loader an
841 auxiliary mechanism may be needed to express that a dynamic linked
842 application depends on a new enough GMP.
844 An auxiliary mechanism may also be needed to express that
845 `libgmpxx.la' (from `--enable-cxx', *note Build Options::) requires
846 `libgmp.la' from the same GMP version, since this is not done by the
847 libtool versioning, nor otherwise. A mismatch will result in
848 unresolved symbols from the linker, or perhaps the loader.
850 When building a package for a CPU family, care should be taken to use
851 `--host' (or `--build') to choose the least common denominator among
852 the CPUs which might use the package. For example this might mean plain
853 `sparc' (meaning V7) for SPARCs.
855 For x86s, `--enable-fat' sets things up for a fat binary build,
856 making a runtime selection of optimized low level routines. This is a
857 good choice for packaging to run on a range of x86 chips.
859 Users who care about speed will want GMP built for their exact CPU
860 type, to make best use of the available optimizations. Providing a way
861 to suitably rebuild a package may be useful. This could be as simple
862 as making it possible for a user to omit `--build' (and `--host') so
863 `./config.guess' will detect the CPU. But a way to manually specify a
864 `--build' will be wanted for systems where `./config.guess' is inexact.
866 On systems with multiple ABIs, a packaged build will need to decide
867 which among the choices is to be provided, see *Note ABI and ISA::. A
868 given run of `./configure' etc will only build one ABI. If a second
869 ABI is also required then a second run of `./configure' etc must be
870 made, starting from a clean directory tree (`make distclean').
872 As noted under "ABI and ISA", currently no attempt is made to follow
873 system conventions for install locations that vary with ABI, such as
874 `/usr/lib/sparcv9' for `ABI=64' as opposed to `/usr/lib' for `ABI=32'.
875 A package build can override `libdir' and other standard variables as
878 Note that `gmp.h' is a generated file, and will be architecture and
879 ABI dependent. When attempting to install two ABIs simultaneously it
880 will be important that an application compile gets the correct `gmp.h'
881 for its desired ABI. If compiler include paths don't vary with ABI
882 options then it might be necessary to create a `/usr/include/gmp.h'
883 which tests preprocessor symbols and chooses the correct actual `gmp.h'.
886 File: gmp.info, Node: Notes for Particular Systems, Next: Known Build Problems, Prev: Notes for Package Builds, Up: Installing GMP
888 2.4 Notes for Particular Systems
889 ================================
892 On systems `*-*-aix[34]*' shared libraries are disabled by
893 default, since some versions of the native `ar' fail on the
894 convenience libraries used. A shared build can be attempted with
896 ./configure --enable-shared --disable-static
898 Note that the `--disable-static' is necessary because in a shared
899 build libtool makes `libgmp.a' a symlink to `libgmp.so',
900 apparently for the benefit of old versions of `ld' which only
901 recognise `.a', but unfortunately this is done even if a fully
902 functional `ld' is available.
905 On systems `arm*-*-*', versions of GCC up to and including 2.95.3
906 have a bug in unsigned division, giving wrong results for some
907 operands. GMP `./configure' will demand GCC 2.95.4 or later.
910 Compaq C++ on OSF 5.1 has two flavours of `iostream', a standard
911 one and an old pre-standard one (see `man iostream_intro'). GMP
912 can only use the standard one, which unfortunately is not the
913 default but must be selected by defining `__USE_STD_IOSTREAM'.
914 Configure with for instance
916 ./configure --enable-cxx CPPFLAGS=-D__USE_STD_IOSTREAM
919 On some systems, the hardware floating point has a control mode
920 which can set all operations to be done in a particular precision,
921 for instance single, double or extended on x86 systems (x87
922 floating point). The GMP functions involving a `double' cannot be
923 expected to operate to their full precision when the hardware is
924 in single precision mode. Of course this affects all code,
925 including application code, not just GMP.
927 MS-DOS and MS Windows
928 On an MS-DOS system DJGPP can be used to build GMP, and on an MS
929 Windows system Cygwin, DJGPP and MINGW can be used. All three are
930 excellent ports of GCC and the various GNU tools.
932 `http://www.cygwin.com/'
933 `http://www.delorie.com/djgpp/'
934 `http://www.mingw.org/'
936 Microsoft also publishes an Interix "Services for Unix" which can
937 be used to build GMP on Windows (with a normal `./configure'), but
938 it's not free software.
941 On systems `*-*-cygwin*', `*-*-mingw*' and `*-*-pw32*' by default
942 GMP builds only a static library, but a DLL can be built instead
945 ./configure --disable-static --enable-shared
947 Static and DLL libraries can't both be built, since certain export
948 directives in `gmp.h' must be different.
950 A MINGW DLL build of GMP can be used with Microsoft C. Libtool
951 doesn't install a `.lib' format import library, but it can be
952 created with MS `lib' as follows, and copied to the install
953 directory. Similarly for `libmp' and `libgmpxx'.
956 lib /def:libgmp-3.dll.def /out:libgmp-3.lib
958 MINGW uses the C runtime library `msvcrt.dll' for I/O, so
959 applications wanting to use the GMP I/O routines must be compiled
960 with `cl /MD' to do the same. If one of the other C runtime
961 library choices provided by MS C is desired then the suggestion is
962 to use the GMP string functions and confine I/O to the application.
964 Motorola 68k CPU Types
965 `m68k' is taken to mean 68000. `m68020' or higher will give a
966 performance boost on applicable CPUs. `m68360' can be used for
967 CPU32 series chips. `m68302' can be used for "Dragonball" series
968 chips, though this is merely a synonym for `m68000'.
971 `m4' in this release of OpenBSD has a bug in `eval' that makes it
972 unsuitable for `.asm' file processing. `./configure' will detect
973 the problem and either abort or choose another m4 in the `PATH'.
974 The bug is fixed in OpenBSD 2.7, so either upgrade or use GNU m4.
977 In GMP, CPU types `power*' and `powerpc*' will each use
978 instructions not available on the other, so it's important to
979 choose the right one for the CPU that will be used. Currently GMP
980 has no assembly code support for using just the common instruction
981 subset. To get executables that run on both, the current
982 suggestion is to use the generic C code (CPU `none'), possibly
983 with appropriate compiler options (like `-mcpu=common' for `gcc').
984 CPU `rs6000' (which is not a CPU but a family of workstations) is
985 accepted by `config.sub', but is currently equivalent to `none'.
988 `sparcv8' or `supersparc' on relevant systems will give a
989 significant performance increase over the V7 code selected by plain
993 The GMP assembly code for both 32-bit and 64-bit Sparc clobbers the
994 "application registers" `g2', `g3' and `g4', the same way that the
995 GCC default `-mapp-regs' does (*note SPARC Options: (gcc)SPARC
998 This makes that code unsuitable for use with the special V9
999 `-mcmodel=embmedany' (which uses `g4' as a data segment pointer),
1000 and for applications wanting to use those registers for special
1001 purposes. In these cases the only suggestion currently is to
1002 build GMP with CPU `none' to avoid the assembly code.
1005 `/usr/bin/m4' lacks various features needed to process `.asm'
1006 files, and instead `./configure' will automatically use
1007 `/usr/5bin/m4', which we believe is always available (if not then
1011 `i586', `pentium' or `pentiummmx' code is good for its intended P5
1012 Pentium chips, but quite slow when run on Intel P6 class chips
1013 (PPro, P-II, P-III). `i386' is a better choice when making
1014 binaries that must run on both.
1016 x86 MMX and SSE2 Code
1017 If the CPU selected has MMX code but the assembler doesn't support
1018 it, a warning is given and non-MMX code is used instead. This
1019 will be an inferior build, since the MMX code that's present is
1020 there because it's faster than the corresponding plain integer
1021 code. The same applies to SSE2.
1023 Old versions of `gas' don't support MMX instructions, in particular
1024 version 1.92.3 that comes with FreeBSD 2.2.8 or the more recent
1025 OpenBSD 3.1 doesn't.
1027 Solaris 2.6 and 2.7 `as' generate incorrect object code for
1028 register to register `movq' instructions, and so can't be used for
1029 MMX code. Install a recent `gas' if MMX code is wanted on these
1033 File: gmp.info, Node: Known Build Problems, Next: Performance optimization, Prev: Notes for Particular Systems, Up: Installing GMP
1035 2.5 Known Build Problems
1036 ========================
1038 You might find more up-to-date information at `http://gmplib.org/'.
1040 Compiler link options
1041 The version of libtool currently in use rather aggressively strips
1042 compiler options when linking a shared library. This will
1043 hopefully be relaxed in the future, but for now if this is a
1044 problem the suggestion is to create a little script to hide them,
1045 and for instance configure with
1047 ./configure CC=gcc-with-my-options
1049 DJGPP (`*-*-msdosdjgpp*')
1050 The DJGPP port of `bash' 2.03 is unable to run the `configure'
1051 script, it exits silently, having died writing a preamble to
1052 `config.log'. Use `bash' 2.04 or higher.
1054 `make all' was found to run out of memory during the final
1055 `libgmp.la' link on one system tested, despite having 64Mb
1056 available. Running `make libgmp.la' directly helped, perhaps
1057 recursing into the various subdirectories uses up memory.
1059 GNU binutils `strip' prior to 2.12
1060 `strip' from GNU binutils 2.11 and earlier should not be used on
1061 the static libraries `libgmp.a' and `libmp.a' since it will
1062 discard all but the last of multiple archive members with the same
1063 name, like the three versions of `init.o' in `libgmp.a'. Binutils
1064 2.12 or higher can be used successfully.
1066 The shared libraries `libgmp.so' and `libmp.so' are not affected by
1067 this and any version of `strip' can be used on them.
1070 On certain versions of SCO OpenServer 5 and IRIX 6.5 the native
1071 `make' is unable to handle the long dependencies list for
1072 `libgmp.la'. The symptom is a "syntax error" on the following
1073 line of the top-level `Makefile'.
1075 libgmp.la: $(libgmp_la_OBJECTS) $(libgmp_la_DEPENDENCIES)
1077 Either use GNU Make, or as a workaround remove
1078 `$(libgmp_la_DEPENDENCIES)' from that line (which will make the
1079 initial build work, but if any recompiling is done `libgmp.la'
1080 might not be rebuilt).
1082 MacOS X (`*-*-darwin*')
1083 Libtool currently only knows how to create shared libraries on
1084 MacOS X using the native `cc' (which is a modified GCC), not a
1085 plain GCC. A static-only build should work though
1086 (`--disable-shared').
1089 The system compiler on old versions of NeXT was a massacred and
1090 old GCC, even if it called itself `cc'. This compiler cannot be
1091 used to build GMP, you need to get a real GCC, and install that.
1092 (NeXT may have fixed this in release 3.3 of their system.)
1095 Bugs in GCC 2.7.2 (and 2.6.3) mean it can't be used to compile GMP
1096 on POWER or PowerPC. If you want to use GCC for these machines,
1097 get GCC 2.7.2.1 (or later).
1100 Use the GNU assembler instead of the system assembler, since the
1101 latter has serious bugs.
1104 The system `sed' prints an error "Output line too long" when
1105 libtool builds `libgmp.la'. This doesn't seem to cause any
1106 obvious ill effects, but GNU `sed' is recommended, to avoid any
1109 Sparc Solaris 2.7 with gcc 2.95.2 in `ABI=32'
1110 A shared library build of GMP seems to fail in this combination,
1111 it builds but then fails the tests, apparently due to some
1112 incorrect data relocations within `gmp_randinit_lc_2exp_size'.
1113 The exact cause is unknown, `--disable-shared' is recommended.
1116 File: gmp.info, Node: Performance optimization, Prev: Known Build Problems, Up: Installing GMP
1118 2.6 Performance optimization
1119 ============================
1121 For optimal performance, build GMP for the exact CPU type of the target
1122 computer, see *Note Build Options::.
1124 Unlike what is the case for most other programs, the compiler
1125 typically doesn't matter much, since GMP uses assembly language for the
1126 most critical operation.
1128 In particular for long-running GMP applications, and applications
1129 demanding extremely large numbers, building and running the `tuneup'
1130 program in the `tune' subdirectory, can be important. For example,
1136 will generate better contents for the `gmp-mparam.h' parameter file.
1138 To use the results, put the output in the file file indicated in the
1139 `Parameters for ...' header. Then recompile from scratch.
1141 The `tuneup' program takes one useful parameter, `-f NNN', which
1142 instructs the program how long to check FFT multiply parameters. If
1143 you're going to use GMP for extremely large numbers, you may want to
1144 run `tuneup' with a large NNN value.
1147 File: gmp.info, Node: GMP Basics, Next: Reporting Bugs, Prev: Installing GMP, Up: Top
1152 *Using functions, macros, data types, etc. not documented in this
1153 manual is strongly discouraged. If you do so your application is
1154 guaranteed to be incompatible with future versions of GMP.*
1158 * Headers and Libraries::
1159 * Nomenclature and Types::
1160 * Function Classes::
1161 * Variable Conventions::
1162 * Parameter Conventions::
1163 * Memory Management::
1165 * Useful Macros and Constants::
1166 * Compatibility with older versions::
1167 * Demonstration Programs::
1175 File: gmp.info, Node: Headers and Libraries, Next: Nomenclature and Types, Prev: GMP Basics, Up: GMP Basics
1177 3.1 Headers and Libraries
1178 =========================
1180 All declarations needed to use GMP are collected in the include file
1181 `gmp.h'. It is designed to work with both C and C++ compilers.
1185 Note however that prototypes for GMP functions with `FILE *'
1186 parameters are only provided if `<stdio.h>' is included too.
1191 Likewise `<stdarg.h>' (or `<varargs.h>') is required for prototypes
1192 with `va_list' parameters, such as `gmp_vprintf'. And `<obstack.h>'
1193 for prototypes with `struct obstack' parameters, such as
1194 `gmp_obstack_printf', when available.
1196 All programs using GMP must link against the `libgmp' library. On a
1197 typical Unix-like system this can be done with `-lgmp', for example
1199 gcc myprogram.c -lgmp
1201 GMP C++ functions are in a separate `libgmpxx' library. This is
1202 built and installed if C++ support has been enabled (*note Build
1203 Options::). For example,
1205 g++ mycxxprog.cc -lgmpxx -lgmp
1207 GMP is built using Libtool and an application can use that to link
1208 if desired, *note GNU Libtool: (libtool)Top.
1210 If GMP has been installed to a non-standard location then it may be
1211 necessary to use `-I' and `-L' compiler options to point to the right
1212 directories, and some sort of run-time path for a shared library.
1215 File: gmp.info, Node: Nomenclature and Types, Next: Function Classes, Prev: Headers and Libraries, Up: GMP Basics
1217 3.2 Nomenclature and Types
1218 ==========================
1220 In this manual, "integer" usually means a multiple precision integer, as
1221 defined by the GMP library. The C data type for such integers is
1222 `mpz_t'. Here are some examples of how to declare such integers:
1226 struct foo { mpz_t x, y; };
1230 "Rational number" means a multiple precision fraction. The C data
1231 type for these fractions is `mpq_t'. For example:
1235 "Floating point number" or "Float" for short, is an arbitrary
1236 precision mantissa with a limited precision exponent. The C data type
1237 for such objects is `mpf_t'. For example:
1241 The floating point functions accept and return exponents in the C
1242 type `mp_exp_t'. Currently this is usually a `long', but on some
1243 systems it's an `int' for efficiency.
1245 A "limb" means the part of a multi-precision number that fits in a
1246 single machine word. (We chose this word because a limb of the human
1247 body is analogous to a digit, only larger, and containing several
1248 digits.) Normally a limb is 32 or 64 bits. The C data type for a limb
1251 Counts of limbs of a multi-precision number represented in the C type
1252 `mp_size_t'. Currently this is normally a `long', but on some systems
1253 it's an `int' for efficiency, and on some systems it will be `long
1254 long' in the future.
1256 Counts of bits of a multi-precision number are represented in the C
1257 type `mp_bitcnt_t'. Currently this is always an `unsigned long', but on
1258 some systems it will be an `unsigned long long' in the future .
1260 "Random state" means an algorithm selection and current state data.
1261 The C data type for such objects is `gmp_randstate_t'. For example:
1263 gmp_randstate_t rstate;
1265 Also, in general `mp_bitcnt_t' is used for bit counts and ranges, and
1266 `size_t' is used for byte or character counts.
1269 File: gmp.info, Node: Function Classes, Next: Variable Conventions, Prev: Nomenclature and Types, Up: GMP Basics
1271 3.3 Function Classes
1272 ====================
1274 There are six classes of functions in the GMP library:
1276 1. Functions for signed integer arithmetic, with names beginning with
1277 `mpz_'. The associated type is `mpz_t'. There are about 150
1278 functions in this class. (*note Integer Functions::)
1280 2. Functions for rational number arithmetic, with names beginning with
1281 `mpq_'. The associated type is `mpq_t'. There are about 40
1282 functions in this class, but the integer functions can be used for
1283 arithmetic on the numerator and denominator separately. (*note
1284 Rational Number Functions::)
1286 3. Functions for floating-point arithmetic, with names beginning with
1287 `mpf_'. The associated type is `mpf_t'. There are about 60
1288 functions is this class. (*note Floating-point Functions::)
1290 4. Functions compatible with Berkeley MP, such as `itom', `madd', and
1291 `mult'. The associated type is `MINT'. (*note BSD Compatible
1294 5. Fast low-level functions that operate on natural numbers. These
1295 are used by the functions in the preceding groups, and you can
1296 also call them directly from very time-critical user programs.
1297 These functions' names begin with `mpn_'. The associated type is
1298 array of `mp_limb_t'. There are about 30 (hard-to-use) functions
1299 in this class. (*note Low-level Functions::)
1301 6. Miscellaneous functions. Functions for setting up custom
1302 allocation and functions for generating random numbers. (*note
1303 Custom Allocation::, and *note Random Number Functions::)
1306 File: gmp.info, Node: Variable Conventions, Next: Parameter Conventions, Prev: Function Classes, Up: GMP Basics
1308 3.4 Variable Conventions
1309 ========================
1311 GMP functions generally have output arguments before input arguments.
1312 This notation is by analogy with the assignment operator. The BSD MP
1313 compatibility functions are exceptions, having the output arguments
1316 GMP lets you use the same variable for both input and output in one
1317 call. For example, the main function for integer multiplication,
1318 `mpz_mul', can be used to square `x' and put the result back in `x' with
1322 Before you can assign to a GMP variable, you need to initialize it
1323 by calling one of the special initialization functions. When you're
1324 done with a variable, you need to clear it out, using one of the
1325 functions for that purpose. Which function to use depends on the type
1326 of variable. See the chapters on integer functions, rational number
1327 functions, and floating-point functions for details.
1329 A variable should only be initialized once, or at least cleared
1330 between each initialization. After a variable has been initialized, it
1331 may be assigned to any number of times.
1333 For efficiency reasons, avoid excessive initializing and clearing.
1334 In general, initialize near the start of a function and clear near the
1343 for (i = 1; i < 100; i++)
1346 mpz_fdiv_q (n, ...);
1353 File: gmp.info, Node: Parameter Conventions, Next: Memory Management, Prev: Variable Conventions, Up: GMP Basics
1355 3.5 Parameter Conventions
1356 =========================
1358 When a GMP variable is used as a function parameter, it's effectively a
1359 call-by-reference, meaning if the function stores a value there it will
1360 change the original in the caller. Parameters which are input-only can
1361 be designated `const' to provoke a compiler error or warning on
1362 attempting to modify them.
1364 When a function is going to return a GMP result, it should designate
1365 a parameter that it sets, like the library functions do. More than one
1366 value can be returned by having more than one output parameter, again
1367 like the library functions. A `return' of an `mpz_t' etc doesn't
1368 return the object, only a pointer, and this is almost certainly not
1371 Here's an example accepting an `mpz_t' parameter, doing a
1372 calculation, and storing the result to the indicated parameter.
1375 foo (mpz_t result, const mpz_t param, unsigned long n)
1378 mpz_mul_ui (result, param, n);
1379 for (i = 1; i < n; i++)
1380 mpz_add_ui (result, result, i*7);
1388 mpz_init_set_str (n, "123456", 0);
1390 gmp_printf ("%Zd\n", r);
1394 `foo' works even if the mainline passes the same variable for
1395 `param' and `result', just like the library functions. But sometimes
1396 it's tricky to make that work, and an application might not want to
1397 bother supporting that sort of thing.
1399 For interest, the GMP types `mpz_t' etc are implemented as
1400 one-element arrays of certain structures. This is why declaring a
1401 variable creates an object with the fields GMP needs, but then using it
1402 as a parameter passes a pointer to the object. Note that the actual
1403 fields in each `mpz_t' etc are for internal use only and should not be
1404 accessed directly by code that expects to be compatible with future GMP
1408 File: gmp.info, Node: Memory Management, Next: Reentrancy, Prev: Parameter Conventions, Up: GMP Basics
1410 3.6 Memory Management
1411 =====================
1413 The GMP types like `mpz_t' are small, containing only a couple of sizes,
1414 and pointers to allocated data. Once a variable is initialized, GMP
1415 takes care of all space allocation. Additional space is allocated
1416 whenever a variable doesn't have enough.
1418 `mpz_t' and `mpq_t' variables never reduce their allocated space.
1419 Normally this is the best policy, since it avoids frequent reallocation.
1420 Applications that need to return memory to the heap at some particular
1421 point can use `mpz_realloc2', or clear variables no longer needed.
1423 `mpf_t' variables, in the current implementation, use a fixed amount
1424 of space, determined by the chosen precision and allocated at
1425 initialization, so their size doesn't change.
1427 All memory is allocated using `malloc' and friends by default, but
1428 this can be changed, see *Note Custom Allocation::. Temporary memory
1429 on the stack is also used (via `alloca'), but this can be changed at
1430 build-time if desired, see *Note Build Options::.
1433 File: gmp.info, Node: Reentrancy, Next: Useful Macros and Constants, Prev: Memory Management, Up: GMP Basics
1438 GMP is reentrant and thread-safe, with some exceptions:
1440 * If configured with `--enable-alloca=malloc-notreentrant' (or with
1441 `--enable-alloca=notreentrant' when `alloca' is not available),
1442 then naturally GMP is not reentrant.
1444 * `mpf_set_default_prec' and `mpf_init' use a global variable for the
1445 selected precision. `mpf_init2' can be used instead, and in the
1446 C++ interface an explicit precision to the `mpf_class' constructor.
1448 * `mpz_random' and the other old random number functions use a global
1449 random state and are hence not reentrant. The newer random number
1450 functions that accept a `gmp_randstate_t' parameter can be used
1453 * `gmp_randinit' (obsolete) returns an error indication through a
1454 global variable, which is not thread safe. Applications are
1455 advised to use `gmp_randinit_default' or `gmp_randinit_lc_2exp'
1458 * `mp_set_memory_functions' uses global variables to store the
1459 selected memory allocation functions.
1461 * If the memory allocation functions set by a call to
1462 `mp_set_memory_functions' (or `malloc' and friends by default) are
1463 not reentrant, then GMP will not be reentrant either.
1465 * If the standard I/O functions such as `fwrite' are not reentrant
1466 then the GMP I/O functions using them will not be reentrant either.
1468 * It's safe for two threads to read from the same GMP variable
1469 simultaneously, but it's not safe for one to read while the
1470 another might be writing, nor for two threads to write
1471 simultaneously. It's not safe for two threads to generate a
1472 random number from the same `gmp_randstate_t' simultaneously,
1473 since this involves an update of that variable.
1476 File: gmp.info, Node: Useful Macros and Constants, Next: Compatibility with older versions, Prev: Reentrancy, Up: GMP Basics
1478 3.8 Useful Macros and Constants
1479 ===============================
1481 -- Global Constant: const int mp_bits_per_limb
1482 The number of bits per limb.
1484 -- Macro: __GNU_MP_VERSION
1485 -- Macro: __GNU_MP_VERSION_MINOR
1486 -- Macro: __GNU_MP_VERSION_PATCHLEVEL
1487 The major and minor GMP version, and patch level, respectively, as
1488 integers. For GMP i.j, these numbers will be i, j, and 0,
1489 respectively. For GMP i.j.k, these numbers will be i, j, and k,
1492 -- Global Constant: const char * const gmp_version
1493 The GMP version number, as a null-terminated string, in the form
1494 "i.j.k". This release is "5.0.1". Note that the format "i.j" was
1495 used when k was zero was used before version 4.3.0.
1498 -- Macro: __GMP_CFLAGS
1499 The compiler and compiler flags, respectively, used when compiling
1503 File: gmp.info, Node: Compatibility with older versions, Next: Demonstration Programs, Prev: Useful Macros and Constants, Up: GMP Basics
1505 3.9 Compatibility with older versions
1506 =====================================
1508 This version of GMP is upwardly binary compatible with all 4.x and 3.x
1509 versions, and upwardly compatible at the source level with all 2.x
1510 versions, with the following exceptions.
1512 * `mpn_gcd' had its source arguments swapped as of GMP 3.0, for
1513 consistency with other `mpn' functions.
1515 * `mpf_get_prec' counted precision slightly differently in GMP 3.0
1516 and 3.0.1, but in 3.1 reverted to the 2.x style.
1518 There are a number of compatibility issues between GMP 1 and GMP 2
1519 that of course also apply when porting applications from GMP 1 to GMP
1520 4. Please see the GMP 2 manual for details.
1522 The Berkeley MP compatibility library (*note BSD Compatible
1523 Functions::) is source and binary compatible with the standard `libmp'.
1526 File: gmp.info, Node: Demonstration Programs, Next: Efficiency, Prev: Compatibility with older versions, Up: GMP Basics
1528 3.10 Demonstration programs
1529 ===========================
1531 The `demos' subdirectory has some sample programs using GMP. These
1532 aren't built or installed, but there's a `Makefile' with rules for them.
1538 The following programs are provided
1540 * `pexpr' is an expression evaluator, the program used on the GMP
1543 * The `calc' subdirectory has a similar but simpler evaluator using
1546 * The `expr' subdirectory is yet another expression evaluator, a
1547 library designed for ease of use within a C program. See
1548 `demos/expr/README' for more information.
1550 * `factorize' is a Pollard-Rho factorization program.
1552 * `isprime' is a command-line interface to the `mpz_probab_prime_p'
1555 * `primes' counts or lists primes in an interval, using a sieve.
1557 * `qcn' is an example use of `mpz_kronecker_ui' to estimate quadratic
1560 * The `perl' subdirectory is a comprehensive perl interface to GMP.
1561 See `demos/perl/INSTALL' for more information. Documentation is
1562 in POD format in `demos/perl/GMP.pm'.
1564 As an aside, consideration has been given at various times to some
1565 sort of expression evaluation within the main GMP library. Going
1566 beyond something minimal quickly leads to matters like user-defined
1567 functions, looping, fixnums for control variables, etc, which are
1568 considered outside the scope of GMP (much closer to language
1569 interpreters or compilers, *Note Language Bindings::.) Something
1570 simple for program input convenience may yet be a possibility, a
1571 combination of the `expr' demo and the `pexpr' tree back-end perhaps.
1572 But for now the above evaluators are offered as illustrations.
1575 File: gmp.info, Node: Efficiency, Next: Debugging, Prev: Demonstration Programs, Up: GMP Basics
1581 On small operands, the time for function call overheads and memory
1582 allocation can be significant in comparison to actual calculation.
1583 This is unavoidable in a general purpose variable precision
1584 library, although GMP attempts to be as efficient as it can on
1585 both large and small operands.
1588 On some CPUs, in particular the x86s, the static `libgmp.a' should
1589 be used for maximum speed, since the PIC code in the shared
1590 `libgmp.so' will have a small overhead on each function call and
1591 global data address. For many programs this will be
1592 insignificant, but for long calculations there's a gain to be had.
1594 Initializing and Clearing
1595 Avoid excessive initializing and clearing of variables, since this
1596 can be quite time consuming, especially in comparison to otherwise
1597 fast operations like addition.
1599 A language interpreter might want to keep a free list or stack of
1600 initialized variables ready for use. It should be possible to
1601 integrate something like that with a garbage collector too.
1604 An `mpz_t' or `mpq_t' variable used to hold successively increasing
1605 values will have its memory repeatedly `realloc'ed, which could be
1606 quite slow or could fragment memory, depending on the C library.
1607 If an application can estimate the final size then `mpz_init2' or
1608 `mpz_realloc2' can be called to allocate the necessary space from
1609 the beginning (*note Initializing Integers::).
1611 It doesn't matter if a size set with `mpz_init2' or `mpz_realloc2'
1612 is too small, since all functions will do a further reallocation
1613 if necessary. Badly overestimating memory required will waste
1617 It's up to an application to call functions like `mpz_mul_2exp'
1618 when appropriate. General purpose functions like `mpz_mul' make
1619 no attempt to identify powers of two or other special forms,
1620 because such inputs will usually be very rare and testing every
1621 time would be wasteful.
1623 `ui' and `si' Functions
1624 The `ui' functions and the small number of `si' functions exist for
1625 convenience and should be used where applicable. But if for
1626 example an `mpz_t' contains a value that fits in an `unsigned
1627 long' there's no need extract it and call a `ui' function, just
1628 use the regular `mpz' function.
1631 `mpz_abs', `mpq_abs', `mpf_abs', `mpz_neg', `mpq_neg' and
1632 `mpf_neg' are fast when used for in-place operations like
1633 `mpz_abs(x,x)', since in the current implementation only a single
1634 field of `x' needs changing. On suitable compilers (GCC for
1635 instance) this is inlined too.
1637 `mpz_add_ui', `mpz_sub_ui', `mpf_add_ui' and `mpf_sub_ui' benefit
1638 from an in-place operation like `mpz_add_ui(x,x,y)', since usually
1639 only one or two limbs of `x' will need to be changed. The same
1640 applies to the full precision `mpz_add' etc if `y' is small. If
1641 `y' is big then cache locality may be helped, but that's all.
1643 `mpz_mul' is currently the opposite, a separate destination is
1644 slightly better. A call like `mpz_mul(x,x,y)' will, unless `y' is
1645 only one limb, make a temporary copy of `x' before forming the
1646 result. Normally that copying will only be a tiny fraction of the
1647 time for the multiply, so this is not a particularly important
1650 `mpz_set', `mpq_set', `mpq_set_num', `mpf_set', etc, make no
1651 attempt to recognise a copy of something to itself, so a call like
1652 `mpz_set(x,x)' will be wasteful. Naturally that would never be
1653 written deliberately, but if it might arise from two pointers to
1654 the same object then a test to avoid it might be desirable.
1659 Note that it's never worth introducing extra `mpz_set' calls just
1660 to get in-place operations. If a result should go to a particular
1661 variable then just direct it there and let GMP take care of data
1664 Divisibility Testing (Small Integers)
1665 `mpz_divisible_ui_p' and `mpz_congruent_ui_p' are the best
1666 functions for testing whether an `mpz_t' is divisible by an
1667 individual small integer. They use an algorithm which is faster
1668 than `mpz_tdiv_ui', but which gives no useful information about
1669 the actual remainder, only whether it's zero (or a particular
1672 However when testing divisibility by several small integers, it's
1673 best to take a remainder modulo their product, to save
1674 multi-precision operations. For instance to test whether a number
1675 is divisible by any of 23, 29 or 31 take a remainder modulo
1676 23*29*31 = 20677 and then test that.
1678 The division functions like `mpz_tdiv_q_ui' which give a quotient
1679 as well as a remainder are generally a little slower than the
1680 remainder-only functions like `mpz_tdiv_ui'. If the quotient is
1681 only rarely wanted then it's probably best to just take a
1682 remainder and then go back and calculate the quotient if and when
1683 it's wanted (`mpz_divexact_ui' can be used if the remainder is
1687 The `mpq' functions operate on `mpq_t' values with no common
1688 factors in the numerator and denominator. Common factors are
1689 checked-for and cast out as necessary. In general, cancelling
1690 factors every time is the best approach since it minimizes the
1691 sizes for subsequent operations.
1693 However, applications that know something about the factorization
1694 of the values they're working with might be able to avoid some of
1695 the GCDs used for canonicalization, or swap them for divisions.
1696 For example when multiplying by a prime it's enough to check for
1697 factors of it in the denominator instead of doing a full GCD. Or
1698 when forming a big product it might be known that very little
1699 cancellation will be possible, and so canonicalization can be left
1702 The `mpq_numref' and `mpq_denref' macros give access to the
1703 numerator and denominator to do things outside the scope of the
1704 supplied `mpq' functions. *Note Applying Integer Functions::.
1706 The canonical form for rationals allows mixed-type `mpq_t' and
1707 integer additions or subtractions to be done directly with
1708 multiples of the denominator. This will be somewhat faster than
1709 `mpq_add'. For example,
1712 mpz_add (mpq_numref(q), mpq_numref(q), mpq_denref(q));
1714 /* mpq += unsigned long */
1715 mpz_addmul_ui (mpq_numref(q), mpq_denref(q), 123UL);
1718 mpz_submul (mpq_numref(q), mpq_denref(q), z);
1721 Functions like `mpz_fac_ui', `mpz_fib_ui' and `mpz_bin_uiui' are
1722 designed for calculating isolated values. If a range of values is
1723 wanted it's probably best to call to get a starting point and
1727 Hexadecimal or octal are suggested for input or output in text
1728 form. Power-of-2 bases like these can be converted much more
1729 efficiently than other bases, like decimal. For big numbers
1730 there's usually nothing of particular interest to be seen in the
1731 digits, so the base doesn't matter much.
1733 Maybe we can hope octal will one day become the normal base for
1734 everyday use, as proposed by King Charles XII of Sweden and later
1738 File: gmp.info, Node: Debugging, Next: Profiling, Prev: Efficiency, Up: GMP Basics
1744 Depending on the system, a segmentation violation or bus error
1745 might be the only indication of stack overflow. See
1746 `--enable-alloca' choices in *Note Build Options::, for how to
1749 In new enough versions of GCC, `-fstack-check' may be able to
1750 ensure an overflow is recognised by the system before too much
1751 damage is done, or `-fstack-limit-symbol' or
1752 `-fstack-limit-register' may be able to add checking if the system
1753 itself doesn't do any (*note Options for Code Generation:
1754 (gcc)Code Gen Options.). These options must be added to the
1755 `CFLAGS' used in the GMP build (*note Build Options::), adding
1756 them just to an application will have no effect. Note also
1757 they're a slowdown, adding overhead to each function call and each
1761 The most likely cause of application problems with GMP is heap
1762 corruption. Failing to `init' GMP variables will have
1763 unpredictable effects, and corruption arising elsewhere in a
1764 program may well affect GMP. Initializing GMP variables more than
1765 once or failing to clear them will cause memory leaks.
1767 In all such cases a `malloc' debugger is recommended. On a GNU or
1768 BSD system the standard C library `malloc' has some diagnostic
1769 facilities, see *Note Allocation Debugging: (libc)Allocation
1770 Debugging, or `man 3 malloc'. Other possibilities, in no
1771 particular order, include
1773 `http://www.inf.ethz.ch/personal/biere/projects/ccmalloc/'
1774 `http://dmalloc.com/'
1775 `http://www.perens.com/FreeSoftware/' (electric fence)
1776 `http://packages.debian.org/stable/devel/fda'
1777 `http://www.gnupdate.org/components/leakbug/'
1778 `http://people.redhat.com/~otaylor/memprof/'
1779 `http://www.cbmamiga.demon.co.uk/mpatrol/'
1781 The GMP default allocation routines in `memory.c' also have a
1782 simple sentinel scheme which can be enabled with `#define DEBUG'
1783 in that file. This is mainly designed for detecting buffer
1784 overruns during GMP development, but might find other uses.
1787 On some systems the compiler options GMP uses by default can
1788 interfere with debugging. In particular on x86 and 68k systems
1789 `-fomit-frame-pointer' is used and this generally inhibits stack
1790 backtracing. Recompiling without such options may help while
1791 debugging, though the usual caveats about it potentially moving a
1792 memory problem or hiding a compiler bug will apply.
1794 GDB, the GNU Debugger
1795 A sample `.gdbinit' is included in the distribution, showing how
1796 to call some undocumented dump functions to print GMP variables
1797 from within GDB. Note that these functions shouldn't be used in
1798 final application code since they're undocumented and may be
1799 subject to incompatible changes in future versions of GMP.
1802 GMP has multiple source files with the same name, in different
1803 directories. For example `mpz', `mpq' and `mpf' each have an
1804 `init.c'. If the debugger can't already determine the right one
1805 it may help to build with absolute paths on each C file. One way
1806 to do that is to use a separate object directory with an absolute
1807 path to the source directory.
1810 /my/source/dir/gmp-5.0.1/configure
1812 This works via `VPATH', and might require GNU `make'. Alternately
1813 it might be possible to change the `.c.lo' rules appropriately.
1816 The build option `--enable-assert' is available to add some
1817 consistency checks to the library (see *Note Build Options::).
1818 These are likely to be of limited value to most applications.
1819 Assertion failures are just as likely to indicate memory
1820 corruption as a library or compiler bug.
1822 Applications using the low-level `mpn' functions, however, will
1823 benefit from `--enable-assert' since it adds checks on the
1824 parameters of most such functions, many of which have subtle
1825 restrictions on their usage. Note however that only the generic C
1826 code has checks, not the assembly code, so CPU `none' should be
1827 used for maximum checking.
1829 Temporary Memory Checking
1830 The build option `--enable-alloca=debug' arranges that each block
1831 of temporary memory in GMP is allocated with a separate call to
1832 `malloc' (or the allocation function set with
1833 `mp_set_memory_functions').
1835 This can help a malloc debugger detect accesses outside the
1836 intended bounds, or detect memory not released. In a normal
1837 build, on the other hand, temporary memory is allocated in blocks
1838 which GMP divides up for its own use, or may be allocated with a
1839 compiler builtin `alloca' which will go nowhere near any malloc
1842 Maximum Debuggability
1843 To summarize the above, a GMP build for maximum debuggability
1846 ./configure --disable-shared --enable-assert \
1847 --enable-alloca=debug --host=none CFLAGS=-g
1849 For C++, add `--enable-cxx CXXFLAGS=-g'.
1852 The GCC checker (`http://savannah.nongnu.org/projects/checker/')
1853 can be used with GMP. It contains a stub library which means GMP
1854 applications compiled with checker can use a normal GMP build.
1856 A build of GMP with checking within GMP itself can be made. This
1857 will run very very slowly. On GNU/Linux for example,
1859 ./configure --host=none-pc-linux-gnu CC=checkergcc
1861 `--host=none' must be used, since the GMP assembly code doesn't
1862 support the checking scheme. The GMP C++ features cannot be used,
1863 since current versions of checker (0.9.9.1) don't yet support the
1864 standard C++ library.
1867 The valgrind program (`http://valgrind.org/') is a memory checker
1868 for x86s. It translates and emulates machine instructions to do
1869 strong checks for uninitialized data (at the level of individual
1870 bits), memory accesses through bad pointers, and memory leaks.
1872 Recent versions of Valgrind are getting support for MMX and
1873 SSE/SSE2 instructions, for past versions GMP will need to be
1874 configured not to use those, ie. for an x86 without them (for
1875 instance plain `i486').
1878 Any suspected bug in GMP itself should be isolated to make sure
1879 it's not an application problem, see *Note Reporting Bugs::.
1882 File: gmp.info, Node: Profiling, Next: Autoconf, Prev: Debugging, Up: GMP Basics
1887 Running a program under a profiler is a good way to find where it's
1888 spending most time and where improvements can be best sought. The
1889 profiling choices for a GMP build are as follows.
1891 `--disable-profiling'
1892 The default is to add nothing special for profiling.
1894 It should be possible to just compile the mainline of a program
1895 with `-p' and use `prof' to get a profile consisting of
1896 timer-based sampling of the program counter. Most of the GMP
1897 assembly code has the necessary symbol information.
1899 This approach has the advantage of minimizing interference with
1900 normal program operation, but on most systems the resolution of
1901 the sampling is quite low (10 milliseconds for instance),
1902 requiring long runs to get accurate information.
1904 `--enable-profiling=prof'
1905 Build with support for the system `prof', which means `-p' added
1908 This provides call counting in addition to program counter
1909 sampling, which allows the most frequently called routines to be
1910 identified, and an average time spent in each routine to be
1913 The x86 assembly code has support for this option, but on other
1914 processors the assembly routines will be as if compiled without
1915 `-p' and therefore won't appear in the call counts.
1917 On some systems, such as GNU/Linux, `-p' in fact means `-pg' and in
1918 this case `--enable-profiling=gprof' described below should be used
1921 `--enable-profiling=gprof'
1922 Build with support for `gprof', which means `-pg' added to the
1925 This provides call graph construction in addition to call counting
1926 and program counter sampling, which makes it possible to count
1927 calls coming from different locations. For example the number of
1928 calls to `mpn_mul' from `mpz_mul' versus the number from
1929 `mpf_mul'. The program counter sampling is still flat though, so
1930 only a total time in `mpn_mul' would be accumulated, not a
1931 separate amount for each call site.
1933 The x86 assembly code has support for this option, but on other
1934 processors the assembly routines will be as if compiled without
1935 `-pg' and therefore not be included in the call counts.
1937 On x86 and m68k systems `-pg' and `-fomit-frame-pointer' are
1938 incompatible, so the latter is omitted from the default flags in
1939 that case, which might result in poorer code generation.
1941 Incidentally, it should be possible to use the `gprof' program
1942 with a plain `--enable-profiling=prof' build. But in that case
1943 only the `gprof -p' flat profile and call counts can be expected
1944 to be valid, not the `gprof -q' call graph.
1946 `--enable-profiling=instrument'
1947 Build with the GCC option `-finstrument-functions' added to the
1948 `CFLAGS' (*note Options for Code Generation: (gcc)Code Gen
1951 This inserts special instrumenting calls at the start and end of
1952 each function, allowing exact timing and full call graph
1955 This instrumenting is not normally a standard system feature and
1956 will require support from an external library, such as
1958 `http://sourceforge.net/projects/fnccheck/'
1960 This should be included in `LIBS' during the GMP configure so that
1961 test programs will link. For example,
1963 ./configure --enable-profiling=instrument LIBS=-lfc
1965 On a GNU system the C library provides dummy instrumenting
1966 functions, so programs compiled with this option will link. In
1967 this case it's only necessary to ensure the correct library is
1968 added when linking an application.
1970 The x86 assembly code supports this option, but on other
1971 processors the assembly routines will be as if compiled without
1972 `-finstrument-functions' meaning time spent in them will
1973 effectively be attributed to their caller.
1976 File: gmp.info, Node: Autoconf, Next: Emacs, Prev: Profiling, Up: GMP Basics
1981 Autoconf based applications can easily check whether GMP is installed.
1982 The only thing to be noted is that GMP library symbols from version 3
1983 onwards have prefixes like `__gmpz'. The following therefore would be
1986 AC_CHECK_LIB(gmp, __gmpz_init)
1988 This just uses the default `AC_CHECK_LIB' actions for found or not
1989 found, but an application that must have GMP would want to generate an
1990 error if not found. For example,
1992 AC_CHECK_LIB(gmp, __gmpz_init, ,
1993 [AC_MSG_ERROR([GNU MP not found, see http://gmplib.org/])])
1995 If functions added in some particular version of GMP are required,
1996 then one of those can be used when checking. For example `mpz_mul_si'
1997 was added in GMP 3.1,
1999 AC_CHECK_LIB(gmp, __gmpz_mul_si, ,
2001 [GNU MP not found, or not 3.1 or up, see http://gmplib.org/])])
2003 An alternative would be to test the version number in `gmp.h' using
2004 say `AC_EGREP_CPP'. That would make it possible to test the exact
2005 version, if some particular sub-minor release is known to be necessary.
2007 In general it's recommended that applications should simply demand a
2008 new enough GMP rather than trying to provide supplements for features
2009 not available in past versions.
2011 Occasionally an application will need or want to know the size of a
2012 type at configuration or preprocessing time, not just with `sizeof' in
2013 the code. This can be done in the normal way with `mp_limb_t' etc, but
2014 GMP 4.0 or up is best for this, since prior versions needed certain
2015 `-D' defines on systems using a `long long' limb. The following would
2016 suit Autoconf 2.50 or up,
2018 AC_CHECK_SIZEOF(mp_limb_t, , [#include <gmp.h>])
2021 File: gmp.info, Node: Emacs, Prev: Autoconf, Up: GMP Basics
2026 <C-h C-i> (`info-lookup-symbol') is a good way to find documentation on
2027 C functions while editing (*note Info Documentation Lookup: (emacs)Info
2030 The GMP manual can be included in such lookups by putting the
2031 following in your `.emacs',
2033 (eval-after-load "info-look"
2034 '(let ((mode-value (assoc 'c-mode (assoc 'symbol info-lookup-alist))))
2035 (setcar (nthcdr 3 mode-value)
2036 (cons '("(gmp)Function Index" nil "^ -.* " "\\>")
2037 (nth 3 mode-value)))))
2040 File: gmp.info, Node: Reporting Bugs, Next: Integer Functions, Prev: GMP Basics, Up: Top
2045 If you think you have found a bug in the GMP library, please
2046 investigate it and report it. We have made this library available to
2047 you, and it is not too much to ask you to report the bugs you find.
2049 Before you report a bug, check it's not already addressed in *Note
2050 Known Build Problems::, or perhaps *Note Notes for Particular
2051 Systems::. You may also want to check `http://gmplib.org/' for patches
2054 Please include the following in any report,
2056 * The GMP version number, and if pre-packaged or patched then say so.
2058 * A test program that makes it possible for us to reproduce the bug.
2059 Include instructions on how to run the program.
2061 * A description of what is wrong. If the results are incorrect, in
2062 what way. If you get a crash, say so.
2064 * If you get a crash, include a stack backtrace from the debugger if
2065 it's informative (`where' in `gdb', or `$C' in `adb').
2067 * Please do not send core dumps, executables or `strace's.
2069 * The configuration options you used when building GMP, if any.
2071 * The name of the compiler and its version. For `gcc', get the
2072 version with `gcc -v', otherwise perhaps `what `which cc`', or
2075 * The output from running `uname -a'.
2077 * The output from running `./config.guess', and from running
2078 `./configfsf.guess' (might be the same).
2080 * If the bug is related to `configure', then the compressed contents
2083 * If the bug is related to an `asm' file not assembling, then the
2084 contents of `config.m4' and the offending line or lines from the
2085 temporary `mpn/tmp-<file>.s'.
2087 Please make an effort to produce a self-contained report, with
2088 something definite that can be tested or debugged. Vague queries or
2089 piecemeal messages are difficult to act on and don't help the
2092 It is not uncommon that an observed problem is actually due to a bug
2093 in the compiler; the GMP code tends to explore interesting corners in
2096 If your bug report is good, we will do our best to help you get a
2097 corrected version of the library; if the bug report is poor, we won't
2098 do anything about it (except maybe ask you to send a better report).
2100 Send your report to: <gmp-bugs@gmplib.org>.
2102 If you think something in this manual is unclear, or downright
2103 incorrect, or if the language needs to be improved, please send a note
2104 to the same address.
2107 File: gmp.info, Node: Integer Functions, Next: Rational Number Functions, Prev: Reporting Bugs, Up: Top
2112 This chapter describes the GMP functions for performing integer
2113 arithmetic. These functions start with the prefix `mpz_'.
2115 GMP integers are stored in objects of type `mpz_t'.
2119 * Initializing Integers::
2120 * Assigning Integers::
2121 * Simultaneous Integer Init & Assign::
2122 * Converting Integers::
2123 * Integer Arithmetic::
2124 * Integer Division::
2125 * Integer Exponentiation::
2127 * Number Theoretic Functions::
2128 * Integer Comparisons::
2129 * Integer Logic and Bit Fiddling::
2131 * Integer Random Numbers::
2132 * Integer Import and Export::
2133 * Miscellaneous Integer Functions::
2134 * Integer Special Functions::
2137 File: gmp.info, Node: Initializing Integers, Next: Assigning Integers, Prev: Integer Functions, Up: Integer Functions
2139 5.1 Initialization Functions
2140 ============================
2142 The functions for integer arithmetic assume that all integer objects are
2143 initialized. You do that by calling the function `mpz_init'. For
2150 mpz_add (integ, ...);
2152 mpz_sub (integ, ...);
2154 /* Unless the program is about to exit, do ... */
2158 As you can see, you can store new values any number of times, once an
2159 object is initialized.
2161 -- Function: void mpz_init (mpz_t X)
2162 Initialize X, and set its value to 0.
2164 -- Function: void mpz_inits (mpz_t X, ...)
2165 Initialize a NULL-terminated list of `mpz_t' variables, and set
2168 -- Function: void mpz_init2 (mpz_t X, mp_bitcnt_t N)
2169 Initialize X, with space for N-bit numbers, and set its value to 0.
2170 Calling this function instead of `mpz_init' or `mpz_inits' is never
2171 necessary; reallocation is handled automatically by GMP when
2174 N is only the initial space, X will grow automatically in the
2175 normal way, if necessary, for subsequent values stored.
2176 `mpz_init2' makes it possible to avoid such reallocations if a
2177 maximum size is known in advance.
2179 -- Function: void mpz_clear (mpz_t X)
2180 Free the space occupied by X. Call this function for all `mpz_t'
2181 variables when you are done with them.
2183 -- Function: void mpz_clears (mpz_t X, ...)
2184 Free the space occupied by a NULL-terminated list of `mpz_t'
2187 -- Function: void mpz_realloc2 (mpz_t X, mp_bitcnt_t N)
2188 Change the space allocated for X to N bits. The value in X is
2189 preserved if it fits, or is set to 0 if not.
2191 Calling this function is never necessary; reallocation is handled
2192 automatically by GMP when needed. But this function can be used
2193 to increase the space for a variable in order to avoid repeated
2194 automatic reallocations, or to decrease it to give memory back to
2198 File: gmp.info, Node: Assigning Integers, Next: Simultaneous Integer Init & Assign, Prev: Initializing Integers, Up: Integer Functions
2200 5.2 Assignment Functions
2201 ========================
2203 These functions assign new values to already initialized integers
2204 (*note Initializing Integers::).
2206 -- Function: void mpz_set (mpz_t ROP, mpz_t OP)
2207 -- Function: void mpz_set_ui (mpz_t ROP, unsigned long int OP)
2208 -- Function: void mpz_set_si (mpz_t ROP, signed long int OP)
2209 -- Function: void mpz_set_d (mpz_t ROP, double OP)
2210 -- Function: void mpz_set_q (mpz_t ROP, mpq_t OP)
2211 -- Function: void mpz_set_f (mpz_t ROP, mpf_t OP)
2212 Set the value of ROP from OP.
2214 `mpz_set_d', `mpz_set_q' and `mpz_set_f' truncate OP to make it an
2217 -- Function: int mpz_set_str (mpz_t ROP, char *STR, int BASE)
2218 Set the value of ROP from STR, a null-terminated C string in base
2219 BASE. White space is allowed in the string, and is simply ignored.
2221 The BASE may vary from 2 to 62, or if BASE is 0, then the leading
2222 characters are used: `0x' and `0X' for hexadecimal, `0b' and `0B'
2223 for binary, `0' for octal, or decimal otherwise.
2225 For bases up to 36, case is ignored; upper-case and lower-case
2226 letters have the same value. For bases 37 to 62, upper-case
2227 letter represent the usual 10..35 while lower-case letter
2230 This function returns 0 if the entire string is a valid number in
2231 base BASE. Otherwise it returns -1.
2233 -- Function: void mpz_swap (mpz_t ROP1, mpz_t ROP2)
2234 Swap the values ROP1 and ROP2 efficiently.
2237 File: gmp.info, Node: Simultaneous Integer Init & Assign, Next: Converting Integers, Prev: Assigning Integers, Up: Integer Functions
2239 5.3 Combined Initialization and Assignment Functions
2240 ====================================================
2242 For convenience, GMP provides a parallel series of initialize-and-set
2243 functions which initialize the output and then store the value there.
2244 These functions' names have the form `mpz_init_set...'
2246 Here is an example of using one:
2250 mpz_init_set_str (pie, "3141592653589793238462643383279502884", 10);
2257 Once the integer has been initialized by any of the `mpz_init_set...'
2258 functions, it can be used as the source or destination operand for the
2259 ordinary integer functions. Don't use an initialize-and-set function
2260 on a variable already initialized!
2262 -- Function: void mpz_init_set (mpz_t ROP, mpz_t OP)
2263 -- Function: void mpz_init_set_ui (mpz_t ROP, unsigned long int OP)
2264 -- Function: void mpz_init_set_si (mpz_t ROP, signed long int OP)
2265 -- Function: void mpz_init_set_d (mpz_t ROP, double OP)
2266 Initialize ROP with limb space and set the initial numeric value
2269 -- Function: int mpz_init_set_str (mpz_t ROP, char *STR, int BASE)
2270 Initialize ROP and set its value like `mpz_set_str' (see its
2271 documentation above for details).
2273 If the string is a correct base BASE number, the function returns
2274 0; if an error occurs it returns -1. ROP is initialized even if
2275 an error occurs. (I.e., you have to call `mpz_clear' for it.)
2278 File: gmp.info, Node: Converting Integers, Next: Integer Arithmetic, Prev: Simultaneous Integer Init & Assign, Up: Integer Functions
2280 5.4 Conversion Functions
2281 ========================
2283 This section describes functions for converting GMP integers to
2284 standard C types. Functions for converting _to_ GMP integers are
2285 described in *Note Assigning Integers:: and *Note I/O of Integers::.
2287 -- Function: unsigned long int mpz_get_ui (mpz_t OP)
2288 Return the value of OP as an `unsigned long'.
2290 If OP is too big to fit an `unsigned long' then just the least
2291 significant bits that do fit are returned. The sign of OP is
2292 ignored, only the absolute value is used.
2294 -- Function: signed long int mpz_get_si (mpz_t OP)
2295 If OP fits into a `signed long int' return the value of OP.
2296 Otherwise return the least significant part of OP, with the same
2299 If OP is too big to fit in a `signed long int', the returned
2300 result is probably not very useful. To find out if the value will
2301 fit, use the function `mpz_fits_slong_p'.
2303 -- Function: double mpz_get_d (mpz_t OP)
2304 Convert OP to a `double', truncating if necessary (ie. rounding
2307 If the exponent from the conversion is too big, the result is
2308 system dependent. An infinity is returned where available. A
2309 hardware overflow trap may or may not occur.
2311 -- Function: double mpz_get_d_2exp (signed long int *EXP, mpz_t OP)
2312 Convert OP to a `double', truncating if necessary (ie. rounding
2313 towards zero), and returning the exponent separately.
2315 The return value is in the range 0.5<=abs(D)<1 and the exponent is
2316 stored to `*EXP'. D * 2^EXP is the (truncated) OP value. If OP
2317 is zero, the return is 0.0 and 0 is stored to `*EXP'.
2319 This is similar to the standard C `frexp' function (*note
2320 Normalization Functions: (libc)Normalization Functions.).
2322 -- Function: char * mpz_get_str (char *STR, int BASE, mpz_t OP)
2323 Convert OP to a string of digits in base BASE. The base argument
2324 may vary from 2 to 62 or from -2 to -36.
2326 For BASE in the range 2..36, digits and lower-case letters are
2327 used; for -2..-36, digits and upper-case letters are used; for
2328 37..62, digits, upper-case letters, and lower-case letters (in
2329 that significance order) are used.
2331 If STR is `NULL', the result string is allocated using the current
2332 allocation function (*note Custom Allocation::). The block will be
2333 `strlen(str)+1' bytes, that being exactly enough for the string and
2336 If STR is not `NULL', it should point to a block of storage large
2337 enough for the result, that being `mpz_sizeinbase (OP, BASE) + 2'.
2338 The two extra bytes are for a possible minus sign, and the
2341 A pointer to the result string is returned, being either the
2342 allocated block, or the given STR.
2345 File: gmp.info, Node: Integer Arithmetic, Next: Integer Division, Prev: Converting Integers, Up: Integer Functions
2347 5.5 Arithmetic Functions
2348 ========================
2350 -- Function: void mpz_add (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2351 -- Function: void mpz_add_ui (mpz_t ROP, mpz_t OP1, unsigned long int
2353 Set ROP to OP1 + OP2.
2355 -- Function: void mpz_sub (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2356 -- Function: void mpz_sub_ui (mpz_t ROP, mpz_t OP1, unsigned long int
2358 -- Function: void mpz_ui_sub (mpz_t ROP, unsigned long int OP1, mpz_t
2360 Set ROP to OP1 - OP2.
2362 -- Function: void mpz_mul (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2363 -- Function: void mpz_mul_si (mpz_t ROP, mpz_t OP1, long int OP2)
2364 -- Function: void mpz_mul_ui (mpz_t ROP, mpz_t OP1, unsigned long int
2366 Set ROP to OP1 times OP2.
2368 -- Function: void mpz_addmul (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2369 -- Function: void mpz_addmul_ui (mpz_t ROP, mpz_t OP1, unsigned long
2371 Set ROP to ROP + OP1 times OP2.
2373 -- Function: void mpz_submul (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2374 -- Function: void mpz_submul_ui (mpz_t ROP, mpz_t OP1, unsigned long
2376 Set ROP to ROP - OP1 times OP2.
2378 -- Function: void mpz_mul_2exp (mpz_t ROP, mpz_t OP1, mp_bitcnt_t OP2)
2379 Set ROP to OP1 times 2 raised to OP2. This operation can also be
2380 defined as a left shift by OP2 bits.
2382 -- Function: void mpz_neg (mpz_t ROP, mpz_t OP)
2385 -- Function: void mpz_abs (mpz_t ROP, mpz_t OP)
2386 Set ROP to the absolute value of OP.
2389 File: gmp.info, Node: Integer Division, Next: Integer Exponentiation, Prev: Integer Arithmetic, Up: Integer Functions
2391 5.6 Division Functions
2392 ======================
2394 Division is undefined if the divisor is zero. Passing a zero divisor
2395 to the division or modulo functions (including the modular powering
2396 functions `mpz_powm' and `mpz_powm_ui'), will cause an intentional
2397 division by zero. This lets a program handle arithmetic exceptions in
2398 these functions the same way as for normal C `int' arithmetic.
2400 -- Function: void mpz_cdiv_q (mpz_t Q, mpz_t N, mpz_t D)
2401 -- Function: void mpz_cdiv_r (mpz_t R, mpz_t N, mpz_t D)
2402 -- Function: void mpz_cdiv_qr (mpz_t Q, mpz_t R, mpz_t N, mpz_t D)
2403 -- Function: unsigned long int mpz_cdiv_q_ui (mpz_t Q, mpz_t N,
2404 unsigned long int D)
2405 -- Function: unsigned long int mpz_cdiv_r_ui (mpz_t R, mpz_t N,
2406 unsigned long int D)
2407 -- Function: unsigned long int mpz_cdiv_qr_ui (mpz_t Q, mpz_t R,
2408 mpz_t N, unsigned long int D)
2409 -- Function: unsigned long int mpz_cdiv_ui (mpz_t N,
2410 unsigned long int D)
2411 -- Function: void mpz_cdiv_q_2exp (mpz_t Q, mpz_t N, mp_bitcnt_t B)
2412 -- Function: void mpz_cdiv_r_2exp (mpz_t R, mpz_t N, mp_bitcnt_t B)
2414 -- Function: void mpz_fdiv_q (mpz_t Q, mpz_t N, mpz_t D)
2415 -- Function: void mpz_fdiv_r (mpz_t R, mpz_t N, mpz_t D)
2416 -- Function: void mpz_fdiv_qr (mpz_t Q, mpz_t R, mpz_t N, mpz_t D)
2417 -- Function: unsigned long int mpz_fdiv_q_ui (mpz_t Q, mpz_t N,
2418 unsigned long int D)
2419 -- Function: unsigned long int mpz_fdiv_r_ui (mpz_t R, mpz_t N,
2420 unsigned long int D)
2421 -- Function: unsigned long int mpz_fdiv_qr_ui (mpz_t Q, mpz_t R,
2422 mpz_t N, unsigned long int D)
2423 -- Function: unsigned long int mpz_fdiv_ui (mpz_t N,
2424 unsigned long int D)
2425 -- Function: void mpz_fdiv_q_2exp (mpz_t Q, mpz_t N, mp_bitcnt_t B)
2426 -- Function: void mpz_fdiv_r_2exp (mpz_t R, mpz_t N, mp_bitcnt_t B)
2428 -- Function: void mpz_tdiv_q (mpz_t Q, mpz_t N, mpz_t D)
2429 -- Function: void mpz_tdiv_r (mpz_t R, mpz_t N, mpz_t D)
2430 -- Function: void mpz_tdiv_qr (mpz_t Q, mpz_t R, mpz_t N, mpz_t D)
2431 -- Function: unsigned long int mpz_tdiv_q_ui (mpz_t Q, mpz_t N,
2432 unsigned long int D)
2433 -- Function: unsigned long int mpz_tdiv_r_ui (mpz_t R, mpz_t N,
2434 unsigned long int D)
2435 -- Function: unsigned long int mpz_tdiv_qr_ui (mpz_t Q, mpz_t R,
2436 mpz_t N, unsigned long int D)
2437 -- Function: unsigned long int mpz_tdiv_ui (mpz_t N,
2438 unsigned long int D)
2439 -- Function: void mpz_tdiv_q_2exp (mpz_t Q, mpz_t N, mp_bitcnt_t B)
2440 -- Function: void mpz_tdiv_r_2exp (mpz_t R, mpz_t N, mp_bitcnt_t B)
2442 Divide N by D, forming a quotient Q and/or remainder R. For the
2443 `2exp' functions, D=2^B. The rounding is in three styles, each
2444 suiting different applications.
2446 * `cdiv' rounds Q up towards +infinity, and R will have the
2447 opposite sign to D. The `c' stands for "ceil".
2449 * `fdiv' rounds Q down towards -infinity, and R will have the
2450 same sign as D. The `f' stands for "floor".
2452 * `tdiv' rounds Q towards zero, and R will have the same sign
2453 as N. The `t' stands for "truncate".
2455 In all cases Q and R will satisfy N=Q*D+R, and R will satisfy
2458 The `q' functions calculate only the quotient, the `r' functions
2459 only the remainder, and the `qr' functions calculate both. Note
2460 that for `qr' the same variable cannot be passed for both Q and R,
2461 or results will be unpredictable.
2463 For the `ui' variants the return value is the remainder, and in
2464 fact returning the remainder is all the `div_ui' functions do. For
2465 `tdiv' and `cdiv' the remainder can be negative, so for those the
2466 return value is the absolute value of the remainder.
2468 For the `2exp' variants the divisor is 2^B. These functions are
2469 implemented as right shifts and bit masks, but of course they
2470 round the same as the other functions.
2472 For positive N both `mpz_fdiv_q_2exp' and `mpz_tdiv_q_2exp' are
2473 simple bitwise right shifts. For negative N, `mpz_fdiv_q_2exp' is
2474 effectively an arithmetic right shift treating N as twos complement
2475 the same as the bitwise logical functions do, whereas
2476 `mpz_tdiv_q_2exp' effectively treats N as sign and magnitude.
2478 -- Function: void mpz_mod (mpz_t R, mpz_t N, mpz_t D)
2479 -- Function: unsigned long int mpz_mod_ui (mpz_t R, mpz_t N,
2480 unsigned long int D)
2481 Set R to N `mod' D. The sign of the divisor is ignored; the
2482 result is always non-negative.
2484 `mpz_mod_ui' is identical to `mpz_fdiv_r_ui' above, returning the
2485 remainder as well as setting R. See `mpz_fdiv_ui' above if only
2486 the return value is wanted.
2488 -- Function: void mpz_divexact (mpz_t Q, mpz_t N, mpz_t D)
2489 -- Function: void mpz_divexact_ui (mpz_t Q, mpz_t N, unsigned long D)
2490 Set Q to N/D. These functions produce correct results only when
2491 it is known in advance that D divides N.
2493 These routines are much faster than the other division functions,
2494 and are the best choice when exact division is known to occur, for
2495 example reducing a rational to lowest terms.
2497 -- Function: int mpz_divisible_p (mpz_t N, mpz_t D)
2498 -- Function: int mpz_divisible_ui_p (mpz_t N, unsigned long int D)
2499 -- Function: int mpz_divisible_2exp_p (mpz_t N, mp_bitcnt_t B)
2500 Return non-zero if N is exactly divisible by D, or in the case of
2501 `mpz_divisible_2exp_p' by 2^B.
2503 N is divisible by D if there exists an integer Q satisfying N =
2504 Q*D. Unlike the other division functions, D=0 is accepted and
2505 following the rule it can be seen that only 0 is considered
2508 -- Function: int mpz_congruent_p (mpz_t N, mpz_t C, mpz_t D)
2509 -- Function: int mpz_congruent_ui_p (mpz_t N, unsigned long int C,
2510 unsigned long int D)
2511 -- Function: int mpz_congruent_2exp_p (mpz_t N, mpz_t C, mp_bitcnt_t B)
2512 Return non-zero if N is congruent to C modulo D, or in the case of
2513 `mpz_congruent_2exp_p' modulo 2^B.
2515 N is congruent to C mod D if there exists an integer Q satisfying
2516 N = C + Q*D. Unlike the other division functions, D=0 is accepted
2517 and following the rule it can be seen that N and C are considered
2518 congruent mod 0 only when exactly equal.
2521 File: gmp.info, Node: Integer Exponentiation, Next: Integer Roots, Prev: Integer Division, Up: Integer Functions
2523 5.7 Exponentiation Functions
2524 ============================
2526 -- Function: void mpz_powm (mpz_t ROP, mpz_t BASE, mpz_t EXP, mpz_t
2528 -- Function: void mpz_powm_ui (mpz_t ROP, mpz_t BASE, unsigned long
2530 Set ROP to (BASE raised to EXP) modulo MOD.
2532 Negative EXP is supported if an inverse BASE^-1 mod MOD exists
2533 (see `mpz_invert' in *Note Number Theoretic Functions::). If an
2534 inverse doesn't exist then a divide by zero is raised.
2536 -- Function: void mpz_powm_sec (mpz_t ROP, mpz_t BASE, mpz_t EXP,
2538 Set ROP to (BASE raised to EXP) modulo MOD.
2540 It is required that EXP > 0 and that MOD is odd.
2542 This function is designed to take the same time and have the same
2543 cache access patterns for any two same-size arguments, assuming
2544 that function arguments are placed at the same position and that
2545 the machine state is identical upon function entry. This function
2546 is intended for cryptographic purposes, where resilience to
2547 side-channel attacks is desired.
2549 -- Function: void mpz_pow_ui (mpz_t ROP, mpz_t BASE, unsigned long int
2551 -- Function: void mpz_ui_pow_ui (mpz_t ROP, unsigned long int BASE,
2552 unsigned long int EXP)
2553 Set ROP to BASE raised to EXP. The case 0^0 yields 1.
2556 File: gmp.info, Node: Integer Roots, Next: Number Theoretic Functions, Prev: Integer Exponentiation, Up: Integer Functions
2558 5.8 Root Extraction Functions
2559 =============================
2561 -- Function: int mpz_root (mpz_t ROP, mpz_t OP, unsigned long int N)
2562 Set ROP to the truncated integer part of the Nth root of OP.
2563 Return non-zero if the computation was exact, i.e., if OP is ROP
2566 -- Function: void mpz_rootrem (mpz_t ROOT, mpz_t REM, mpz_t U,
2567 unsigned long int N)
2568 Set ROOT to the truncated integer part of the Nth root of U. Set
2569 REM to the remainder, U-ROOT**N.
2571 -- Function: void mpz_sqrt (mpz_t ROP, mpz_t OP)
2572 Set ROP to the truncated integer part of the square root of OP.
2574 -- Function: void mpz_sqrtrem (mpz_t ROP1, mpz_t ROP2, mpz_t OP)
2575 Set ROP1 to the truncated integer part of the square root of OP,
2576 like `mpz_sqrt'. Set ROP2 to the remainder OP-ROP1*ROP1, which
2577 will be zero if OP is a perfect square.
2579 If ROP1 and ROP2 are the same variable, the results are undefined.
2581 -- Function: int mpz_perfect_power_p (mpz_t OP)
2582 Return non-zero if OP is a perfect power, i.e., if there exist
2583 integers A and B, with B>1, such that OP equals A raised to the
2586 Under this definition both 0 and 1 are considered to be perfect
2587 powers. Negative values of OP are accepted, but of course can
2588 only be odd perfect powers.
2590 -- Function: int mpz_perfect_square_p (mpz_t OP)
2591 Return non-zero if OP is a perfect square, i.e., if the square
2592 root of OP is an integer. Under this definition both 0 and 1 are
2593 considered to be perfect squares.
2596 File: gmp.info, Node: Number Theoretic Functions, Next: Integer Comparisons, Prev: Integer Roots, Up: Integer Functions
2598 5.9 Number Theoretic Functions
2599 ==============================
2601 -- Function: int mpz_probab_prime_p (mpz_t N, int REPS)
2602 Determine whether N is prime. Return 2 if N is definitely prime,
2603 return 1 if N is probably prime (without being certain), or return
2604 0 if N is definitely composite.
2606 This function does some trial divisions, then some Miller-Rabin
2607 probabilistic primality tests. REPS controls how many such tests
2608 are done, 5 to 10 is a reasonable number, more will reduce the
2609 chances of a composite being returned as "probably prime".
2611 Miller-Rabin and similar tests can be more properly called
2612 compositeness tests. Numbers which fail are known to be composite
2613 but those which pass might be prime or might be composite. Only a
2614 few composites pass, hence those which pass are considered
2617 -- Function: void mpz_nextprime (mpz_t ROP, mpz_t OP)
2618 Set ROP to the next prime greater than OP.
2620 This function uses a probabilistic algorithm to identify primes.
2621 For practical purposes it's adequate, the chance of a composite
2622 passing will be extremely small.
2624 -- Function: void mpz_gcd (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2625 Set ROP to the greatest common divisor of OP1 and OP2. The result
2626 is always positive even if one or both input operands are negative.
2628 -- Function: unsigned long int mpz_gcd_ui (mpz_t ROP, mpz_t OP1,
2629 unsigned long int OP2)
2630 Compute the greatest common divisor of OP1 and OP2. If ROP is not
2631 `NULL', store the result there.
2633 If the result is small enough to fit in an `unsigned long int', it
2634 is returned. If the result does not fit, 0 is returned, and the
2635 result is equal to the argument OP1. Note that the result will
2636 always fit if OP2 is non-zero.
2638 -- Function: void mpz_gcdext (mpz_t G, mpz_t S, mpz_t T, mpz_t A,
2640 Set G to the greatest common divisor of A and B, and in addition
2641 set S and T to coefficients satisfying A*S + B*T = G. The value
2642 in G is always positive, even if one or both of A and B are
2643 negative. The values in S and T are chosen such that abs(S) <=
2644 abs(B) and abs(T) <= abs(A).
2646 If T is `NULL' then that value is not computed.
2648 -- Function: void mpz_lcm (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2649 -- Function: void mpz_lcm_ui (mpz_t ROP, mpz_t OP1, unsigned long OP2)
2650 Set ROP to the least common multiple of OP1 and OP2. ROP is
2651 always positive, irrespective of the signs of OP1 and OP2. ROP
2652 will be zero if either OP1 or OP2 is zero.
2654 -- Function: int mpz_invert (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2655 Compute the inverse of OP1 modulo OP2 and put the result in ROP.
2656 If the inverse exists, the return value is non-zero and ROP will
2657 satisfy 0 <= ROP < OP2. If an inverse doesn't exist the return
2658 value is zero and ROP is undefined.
2660 -- Function: int mpz_jacobi (mpz_t A, mpz_t B)
2661 Calculate the Jacobi symbol (A/B). This is defined only for B odd.
2663 -- Function: int mpz_legendre (mpz_t A, mpz_t P)
2664 Calculate the Legendre symbol (A/P). This is defined only for P
2665 an odd positive prime, and for such P it's identical to the Jacobi
2668 -- Function: int mpz_kronecker (mpz_t A, mpz_t B)
2669 -- Function: int mpz_kronecker_si (mpz_t A, long B)
2670 -- Function: int mpz_kronecker_ui (mpz_t A, unsigned long B)
2671 -- Function: int mpz_si_kronecker (long A, mpz_t B)
2672 -- Function: int mpz_ui_kronecker (unsigned long A, mpz_t B)
2673 Calculate the Jacobi symbol (A/B) with the Kronecker extension
2674 (a/2)=(2/a) when a odd, or (a/2)=0 when a even.
2676 When B is odd the Jacobi symbol and Kronecker symbol are
2677 identical, so `mpz_kronecker_ui' etc can be used for mixed
2678 precision Jacobi symbols too.
2680 For more information see Henri Cohen section 1.4.2 (*note
2681 References::), or any number theory textbook. See also the
2682 example program `demos/qcn.c' which uses `mpz_kronecker_ui'.
2684 -- Function: mp_bitcnt_t mpz_remove (mpz_t ROP, mpz_t OP, mpz_t F)
2685 Remove all occurrences of the factor F from OP and store the
2686 result in ROP. The return value is how many such occurrences were
2689 -- Function: void mpz_fac_ui (mpz_t ROP, unsigned long int OP)
2690 Set ROP to OP!, the factorial of OP.
2692 -- Function: void mpz_bin_ui (mpz_t ROP, mpz_t N, unsigned long int K)
2693 -- Function: void mpz_bin_uiui (mpz_t ROP, unsigned long int N,
2694 unsigned long int K)
2695 Compute the binomial coefficient N over K and store the result in
2696 ROP. Negative values of N are supported by `mpz_bin_ui', using
2697 the identity bin(-n,k) = (-1)^k * bin(n+k-1,k), see Knuth volume 1
2698 section 1.2.6 part G.
2700 -- Function: void mpz_fib_ui (mpz_t FN, unsigned long int N)
2701 -- Function: void mpz_fib2_ui (mpz_t FN, mpz_t FNSUB1, unsigned long
2703 `mpz_fib_ui' sets FN to to F[n], the N'th Fibonacci number.
2704 `mpz_fib2_ui' sets FN to F[n], and FNSUB1 to F[n-1].
2706 These functions are designed for calculating isolated Fibonacci
2707 numbers. When a sequence of values is wanted it's best to start
2708 with `mpz_fib2_ui' and iterate the defining F[n+1]=F[n]+F[n-1] or
2711 -- Function: void mpz_lucnum_ui (mpz_t LN, unsigned long int N)
2712 -- Function: void mpz_lucnum2_ui (mpz_t LN, mpz_t LNSUB1, unsigned
2714 `mpz_lucnum_ui' sets LN to to L[n], the N'th Lucas number.
2715 `mpz_lucnum2_ui' sets LN to L[n], and LNSUB1 to L[n-1].
2717 These functions are designed for calculating isolated Lucas
2718 numbers. When a sequence of values is wanted it's best to start
2719 with `mpz_lucnum2_ui' and iterate the defining L[n+1]=L[n]+L[n-1]
2722 The Fibonacci numbers and Lucas numbers are related sequences, so
2723 it's never necessary to call both `mpz_fib2_ui' and
2724 `mpz_lucnum2_ui'. The formulas for going from Fibonacci to Lucas
2725 can be found in *Note Lucas Numbers Algorithm::, the reverse is
2726 straightforward too.
2729 File: gmp.info, Node: Integer Comparisons, Next: Integer Logic and Bit Fiddling, Prev: Number Theoretic Functions, Up: Integer Functions
2731 5.10 Comparison Functions
2732 =========================
2734 -- Function: int mpz_cmp (mpz_t OP1, mpz_t OP2)
2735 -- Function: int mpz_cmp_d (mpz_t OP1, double OP2)
2736 -- Macro: int mpz_cmp_si (mpz_t OP1, signed long int OP2)
2737 -- Macro: int mpz_cmp_ui (mpz_t OP1, unsigned long int OP2)
2738 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero
2739 if OP1 = OP2, or a negative value if OP1 < OP2.
2741 `mpz_cmp_ui' and `mpz_cmp_si' are macros and will evaluate their
2742 arguments more than once. `mpz_cmp_d' can be called with an
2743 infinity, but results are undefined for a NaN.
2745 -- Function: int mpz_cmpabs (mpz_t OP1, mpz_t OP2)
2746 -- Function: int mpz_cmpabs_d (mpz_t OP1, double OP2)
2747 -- Function: int mpz_cmpabs_ui (mpz_t OP1, unsigned long int OP2)
2748 Compare the absolute values of OP1 and OP2. Return a positive
2749 value if abs(OP1) > abs(OP2), zero if abs(OP1) = abs(OP2), or a
2750 negative value if abs(OP1) < abs(OP2).
2752 `mpz_cmpabs_d' can be called with an infinity, but results are
2753 undefined for a NaN.
2755 -- Macro: int mpz_sgn (mpz_t OP)
2756 Return +1 if OP > 0, 0 if OP = 0, and -1 if OP < 0.
2758 This function is actually implemented as a macro. It evaluates
2759 its argument multiple times.
2762 File: gmp.info, Node: Integer Logic and Bit Fiddling, Next: I/O of Integers, Prev: Integer Comparisons, Up: Integer Functions
2764 5.11 Logical and Bit Manipulation Functions
2765 ===========================================
2767 These functions behave as if twos complement arithmetic were used
2768 (although sign-magnitude is the actual implementation). The least
2769 significant bit is number 0.
2771 -- Function: void mpz_and (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2772 Set ROP to OP1 bitwise-and OP2.
2774 -- Function: void mpz_ior (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2775 Set ROP to OP1 bitwise inclusive-or OP2.
2777 -- Function: void mpz_xor (mpz_t ROP, mpz_t OP1, mpz_t OP2)
2778 Set ROP to OP1 bitwise exclusive-or OP2.
2780 -- Function: void mpz_com (mpz_t ROP, mpz_t OP)
2781 Set ROP to the one's complement of OP.
2783 -- Function: mp_bitcnt_t mpz_popcount (mpz_t OP)
2784 If OP>=0, return the population count of OP, which is the number
2785 of 1 bits in the binary representation. If OP<0, the number of 1s
2786 is infinite, and the return value is the largest possible
2789 -- Function: mp_bitcnt_t mpz_hamdist (mpz_t OP1, mpz_t OP2)
2790 If OP1 and OP2 are both >=0 or both <0, return the hamming
2791 distance between the two operands, which is the number of bit
2792 positions where OP1 and OP2 have different bit values. If one
2793 operand is >=0 and the other <0 then the number of bits different
2794 is infinite, and the return value is the largest possible
2797 -- Function: mp_bitcnt_t mpz_scan0 (mpz_t OP, mp_bitcnt_t STARTING_BIT)
2798 -- Function: mp_bitcnt_t mpz_scan1 (mpz_t OP, mp_bitcnt_t STARTING_BIT)
2799 Scan OP, starting from bit STARTING_BIT, towards more significant
2800 bits, until the first 0 or 1 bit (respectively) is found. Return
2801 the index of the found bit.
2803 If the bit at STARTING_BIT is already what's sought, then
2804 STARTING_BIT is returned.
2806 If there's no bit found, then the largest possible `mp_bitcnt_t' is
2807 returned. This will happen in `mpz_scan0' past the end of a
2808 negative number, or `mpz_scan1' past the end of a nonnegative
2811 -- Function: void mpz_setbit (mpz_t ROP, mp_bitcnt_t BIT_INDEX)
2812 Set bit BIT_INDEX in ROP.
2814 -- Function: void mpz_clrbit (mpz_t ROP, mp_bitcnt_t BIT_INDEX)
2815 Clear bit BIT_INDEX in ROP.
2817 -- Function: void mpz_combit (mpz_t ROP, mp_bitcnt_t BIT_INDEX)
2818 Complement bit BIT_INDEX in ROP.
2820 -- Function: int mpz_tstbit (mpz_t OP, mp_bitcnt_t BIT_INDEX)
2821 Test bit BIT_INDEX in OP and return 0 or 1 accordingly.
2824 File: gmp.info, Node: I/O of Integers, Next: Integer Random Numbers, Prev: Integer Logic and Bit Fiddling, Up: Integer Functions
2826 5.12 Input and Output Functions
2827 ===============================
2829 Functions that perform input from a stdio stream, and functions that
2830 output to a stdio stream. Passing a `NULL' pointer for a STREAM
2831 argument to any of these functions will make them read from `stdin' and
2832 write to `stdout', respectively.
2834 When using any of these functions, it is a good idea to include
2835 `stdio.h' before `gmp.h', since that will allow `gmp.h' to define
2836 prototypes for these functions.
2838 -- Function: size_t mpz_out_str (FILE *STREAM, int BASE, mpz_t OP)
2839 Output OP on stdio stream STREAM, as a string of digits in base
2840 BASE. The base argument may vary from 2 to 62 or from -2 to -36.
2842 For BASE in the range 2..36, digits and lower-case letters are
2843 used; for -2..-36, digits and upper-case letters are used; for
2844 37..62, digits, upper-case letters, and lower-case letters (in
2845 that significance order) are used.
2847 Return the number of bytes written, or if an error occurred,
2850 -- Function: size_t mpz_inp_str (mpz_t ROP, FILE *STREAM, int BASE)
2851 Input a possibly white-space preceded string in base BASE from
2852 stdio stream STREAM, and put the read integer in ROP.
2854 The BASE may vary from 2 to 62, or if BASE is 0, then the leading
2855 characters are used: `0x' and `0X' for hexadecimal, `0b' and `0B'
2856 for binary, `0' for octal, or decimal otherwise.
2858 For bases up to 36, case is ignored; upper-case and lower-case
2859 letters have the same value. For bases 37 to 62, upper-case
2860 letter represent the usual 10..35 while lower-case letter
2863 Return the number of bytes read, or if an error occurred, return 0.
2865 -- Function: size_t mpz_out_raw (FILE *STREAM, mpz_t OP)
2866 Output OP on stdio stream STREAM, in raw binary format. The
2867 integer is written in a portable format, with 4 bytes of size
2868 information, and that many bytes of limbs. Both the size and the
2869 limbs are written in decreasing significance order (i.e., in
2872 The output can be read with `mpz_inp_raw'.
2874 Return the number of bytes written, or if an error occurred,
2877 The output of this can not be read by `mpz_inp_raw' from GMP 1,
2878 because of changes necessary for compatibility between 32-bit and
2881 -- Function: size_t mpz_inp_raw (mpz_t ROP, FILE *STREAM)
2882 Input from stdio stream STREAM in the format written by
2883 `mpz_out_raw', and put the result in ROP. Return the number of
2884 bytes read, or if an error occurred, return 0.
2886 This routine can read the output from `mpz_out_raw' also from GMP
2887 1, in spite of changes necessary for compatibility between 32-bit
2888 and 64-bit machines.
2891 File: gmp.info, Node: Integer Random Numbers, Next: Integer Import and Export, Prev: I/O of Integers, Up: Integer Functions
2893 5.13 Random Number Functions
2894 ============================
2896 The random number functions of GMP come in two groups; older function
2897 that rely on a global state, and newer functions that accept a state
2898 parameter that is read and modified. Please see the *Note Random
2899 Number Functions:: for more information on how to use and not to use
2900 random number functions.
2902 -- Function: void mpz_urandomb (mpz_t ROP, gmp_randstate_t STATE,
2904 Generate a uniformly distributed random integer in the range 0 to
2907 The variable STATE must be initialized by calling one of the
2908 `gmp_randinit' functions (*Note Random State Initialization::)
2909 before invoking this function.
2911 -- Function: void mpz_urandomm (mpz_t ROP, gmp_randstate_t STATE,
2913 Generate a uniform random integer in the range 0 to N-1, inclusive.
2915 The variable STATE must be initialized by calling one of the
2916 `gmp_randinit' functions (*Note Random State Initialization::)
2917 before invoking this function.
2919 -- Function: void mpz_rrandomb (mpz_t ROP, gmp_randstate_t STATE,
2921 Generate a random integer with long strings of zeros and ones in
2922 the binary representation. Useful for testing functions and
2923 algorithms, since this kind of random numbers have proven to be
2924 more likely to trigger corner-case bugs. The random number will
2925 be in the range 0 to 2^N-1, inclusive.
2927 The variable STATE must be initialized by calling one of the
2928 `gmp_randinit' functions (*Note Random State Initialization::)
2929 before invoking this function.
2931 -- Function: void mpz_random (mpz_t ROP, mp_size_t MAX_SIZE)
2932 Generate a random integer of at most MAX_SIZE limbs. The generated
2933 random number doesn't satisfy any particular requirements of
2934 randomness. Negative random numbers are generated when MAX_SIZE
2937 This function is obsolete. Use `mpz_urandomb' or `mpz_urandomm'
2940 -- Function: void mpz_random2 (mpz_t ROP, mp_size_t MAX_SIZE)
2941 Generate a random integer of at most MAX_SIZE limbs, with long
2942 strings of zeros and ones in the binary representation. Useful
2943 for testing functions and algorithms, since this kind of random
2944 numbers have proven to be more likely to trigger corner-case bugs.
2945 Negative random numbers are generated when MAX_SIZE is negative.
2947 This function is obsolete. Use `mpz_rrandomb' instead.
2950 File: gmp.info, Node: Integer Import and Export, Next: Miscellaneous Integer Functions, Prev: Integer Random Numbers, Up: Integer Functions
2952 5.14 Integer Import and Export
2953 ==============================
2955 `mpz_t' variables can be converted to and from arbitrary words of binary
2956 data with the following functions.
2958 -- Function: void mpz_import (mpz_t ROP, size_t COUNT, int ORDER,
2959 size_t SIZE, int ENDIAN, size_t NAILS, const void *OP)
2960 Set ROP from an array of word data at OP.
2962 The parameters specify the format of the data. COUNT many words
2963 are read, each SIZE bytes. ORDER can be 1 for most significant
2964 word first or -1 for least significant first. Within each word
2965 ENDIAN can be 1 for most significant byte first, -1 for least
2966 significant first, or 0 for the native endianness of the host CPU.
2967 The most significant NAILS bits of each word are skipped, this
2968 can be 0 to use the full words.
2970 There is no sign taken from the data, ROP will simply be a positive
2971 integer. An application can handle any sign itself, and apply it
2972 for instance with `mpz_neg'.
2974 There are no data alignment restrictions on OP, any address is
2977 Here's an example converting an array of `unsigned long' data, most
2978 significant element first, and host byte order within each value.
2980 unsigned long a[20];
2981 /* Initialize Z and A */
2982 mpz_import (z, 20, 1, sizeof(a[0]), 0, 0, a);
2984 This example assumes the full `sizeof' bytes are used for data in
2985 the given type, which is usually true, and certainly true for
2986 `unsigned long' everywhere we know of. However on Cray vector
2987 systems it may be noted that `short' and `int' are always stored
2988 in 8 bytes (and with `sizeof' indicating that) but use only 32 or
2989 46 bits. The NAILS feature can account for this, by passing for
2990 instance `8*sizeof(int)-INT_BIT'.
2992 -- Function: void * mpz_export (void *ROP, size_t *COUNTP, int ORDER,
2993 size_t SIZE, int ENDIAN, size_t NAILS, mpz_t OP)
2994 Fill ROP with word data from OP.
2996 The parameters specify the format of the data produced. Each word
2997 will be SIZE bytes and ORDER can be 1 for most significant word
2998 first or -1 for least significant first. Within each word ENDIAN
2999 can be 1 for most significant byte first, -1 for least significant
3000 first, or 0 for the native endianness of the host CPU. The most
3001 significant NAILS bits of each word are unused and set to zero,
3002 this can be 0 to produce full words.
3004 The number of words produced is written to `*COUNTP', or COUNTP
3005 can be `NULL' to discard the count. ROP must have enough space
3006 for the data, or if ROP is `NULL' then a result array of the
3007 necessary size is allocated using the current GMP allocation
3008 function (*note Custom Allocation::). In either case the return
3009 value is the destination used, either ROP or the allocated block.
3011 If OP is non-zero then the most significant word produced will be
3012 non-zero. If OP is zero then the count returned will be zero and
3013 nothing written to ROP. If ROP is `NULL' in this case, no block
3014 is allocated, just `NULL' is returned.
3016 The sign of OP is ignored, just the absolute value is exported. An
3017 application can use `mpz_sgn' to get the sign and handle it as
3018 desired. (*note Integer Comparisons::)
3020 There are no data alignment restrictions on ROP, any address is
3023 When an application is allocating space itself the required size
3024 can be determined with a calculation like the following. Since
3025 `mpz_sizeinbase' always returns at least 1, `count' here will be
3026 at least one, which avoids any portability problems with
3027 `malloc(0)', though if `z' is zero no space at all is actually
3028 needed (or written).
3030 numb = 8*size - nail;
3031 count = (mpz_sizeinbase (z, 2) + numb-1) / numb;
3032 p = malloc (count * size);
3035 File: gmp.info, Node: Miscellaneous Integer Functions, Next: Integer Special Functions, Prev: Integer Import and Export, Up: Integer Functions
3037 5.15 Miscellaneous Functions
3038 ============================
3040 -- Function: int mpz_fits_ulong_p (mpz_t OP)
3041 -- Function: int mpz_fits_slong_p (mpz_t OP)
3042 -- Function: int mpz_fits_uint_p (mpz_t OP)
3043 -- Function: int mpz_fits_sint_p (mpz_t OP)
3044 -- Function: int mpz_fits_ushort_p (mpz_t OP)
3045 -- Function: int mpz_fits_sshort_p (mpz_t OP)
3046 Return non-zero iff the value of OP fits in an `unsigned long int',
3047 `signed long int', `unsigned int', `signed int', `unsigned short
3048 int', or `signed short int', respectively. Otherwise, return zero.
3050 -- Macro: int mpz_odd_p (mpz_t OP)
3051 -- Macro: int mpz_even_p (mpz_t OP)
3052 Determine whether OP is odd or even, respectively. Return
3053 non-zero if yes, zero if no. These macros evaluate their argument
3056 -- Function: size_t mpz_sizeinbase (mpz_t OP, int BASE)
3057 Return the size of OP measured in number of digits in the given
3058 BASE. BASE can vary from 2 to 62. The sign of OP is ignored,
3059 just the absolute value is used. The result will be either exact
3060 or 1 too big. If BASE is a power of 2, the result is always
3061 exact. If OP is zero the return value is always 1.
3063 This function can be used to determine the space required when
3064 converting OP to a string. The right amount of allocation is
3065 normally two more than the value returned by `mpz_sizeinbase', one
3066 extra for a minus sign and one for the null-terminator.
3068 It will be noted that `mpz_sizeinbase(OP,2)' can be used to locate
3069 the most significant 1 bit in OP, counting from 1. (Unlike the
3070 bitwise functions which start from 0, *Note Logical and Bit
3071 Manipulation Functions: Integer Logic and Bit Fiddling.)
3074 File: gmp.info, Node: Integer Special Functions, Prev: Miscellaneous Integer Functions, Up: Integer Functions
3076 5.16 Special Functions
3077 ======================
3079 The functions in this section are for various special purposes. Most
3080 applications will not need them.
3082 -- Function: void mpz_array_init (mpz_t INTEGER_ARRAY, mp_size_t
3083 ARRAY_SIZE, mp_size_t FIXED_NUM_BITS)
3084 This is a special type of initialization. *Fixed* space of
3085 FIXED_NUM_BITS is allocated to each of the ARRAY_SIZE integers in
3086 INTEGER_ARRAY. There is no way to free the storage allocated by
3087 this function. Don't call `mpz_clear'!
3089 The INTEGER_ARRAY parameter is the first `mpz_t' in the array. For
3093 mpz_array_init (arr[0], 20000, 512);
3095 This function is only intended for programs that create a large
3096 number of integers and need to reduce memory usage by avoiding the
3097 overheads of allocating and reallocating lots of small blocks. In
3098 normal programs this function is not recommended.
3100 The space allocated to each integer by this function will not be
3101 automatically increased, unlike the normal `mpz_init', so an
3102 application must ensure it is sufficient for any value stored.
3103 The following space requirements apply to various routines,
3105 * `mpz_abs', `mpz_neg', `mpz_set', `mpz_set_si' and
3106 `mpz_set_ui' need room for the value they store.
3108 * `mpz_add', `mpz_add_ui', `mpz_sub' and `mpz_sub_ui' need room
3109 for the larger of the two operands, plus an extra
3112 * `mpz_mul', `mpz_mul_ui' and `mpz_mul_ui' need room for the sum
3113 of the number of bits in their operands, but each rounded up
3114 to a multiple of `mp_bits_per_limb'.
3116 * `mpz_swap' can be used between two array variables, but not
3117 between an array and a normal variable.
3119 For other functions, or if in doubt, the suggestion is to
3120 calculate in a regular `mpz_init' variable and copy the result to
3121 an array variable with `mpz_set'.
3123 -- Function: void * _mpz_realloc (mpz_t INTEGER, mp_size_t NEW_ALLOC)
3124 Change the space for INTEGER to NEW_ALLOC limbs. The value in
3125 INTEGER is preserved if it fits, or is set to 0 if not. The return
3126 value is not useful to applications and should be ignored.
3128 `mpz_realloc2' is the preferred way to accomplish allocation
3129 changes like this. `mpz_realloc2' and `_mpz_realloc' are the same
3130 except that `_mpz_realloc' takes its size in limbs.
3132 -- Function: mp_limb_t mpz_getlimbn (mpz_t OP, mp_size_t N)
3133 Return limb number N from OP. The sign of OP is ignored, just the
3134 absolute value is used. The least significant limb is number 0.
3136 `mpz_size' can be used to find how many limbs make up OP.
3137 `mpz_getlimbn' returns zero if N is outside the range 0 to
3140 -- Function: size_t mpz_size (mpz_t OP)
3141 Return the size of OP measured in number of limbs. If OP is zero,
3142 the returned value will be zero.
3145 File: gmp.info, Node: Rational Number Functions, Next: Floating-point Functions, Prev: Integer Functions, Up: Top
3147 6 Rational Number Functions
3148 ***************************
3150 This chapter describes the GMP functions for performing arithmetic on
3151 rational numbers. These functions start with the prefix `mpq_'.
3153 Rational numbers are stored in objects of type `mpq_t'.
3155 All rational arithmetic functions assume operands have a canonical
3156 form, and canonicalize their result. The canonical from means that the
3157 denominator and the numerator have no common factors, and that the
3158 denominator is positive. Zero has the unique representation 0/1.
3160 Pure assignment functions do not canonicalize the assigned variable.
3161 It is the responsibility of the user to canonicalize the assigned
3162 variable before any arithmetic operations are performed on that
3165 -- Function: void mpq_canonicalize (mpq_t OP)
3166 Remove any factors that are common to the numerator and
3167 denominator of OP, and make the denominator positive.
3171 * Initializing Rationals::
3172 * Rational Conversions::
3173 * Rational Arithmetic::
3174 * Comparing Rationals::
3175 * Applying Integer Functions::
3176 * I/O of Rationals::
3179 File: gmp.info, Node: Initializing Rationals, Next: Rational Conversions, Prev: Rational Number Functions, Up: Rational Number Functions
3181 6.1 Initialization and Assignment Functions
3182 ===========================================
3184 -- Function: void mpq_init (mpq_t X)
3185 Initialize X and set it to 0/1. Each variable should normally
3186 only be initialized once, or at least cleared out (using the
3187 function `mpq_clear') between each initialization.
3189 -- Function: void mpq_inits (mpq_t X, ...)
3190 Initialize a NULL-terminated list of `mpq_t' variables, and set
3191 their values to 0/1.
3193 -- Function: void mpq_clear (mpq_t X)
3194 Free the space occupied by X. Make sure to call this function for
3195 all `mpq_t' variables when you are done with them.
3197 -- Function: void mpq_clears (mpq_t X, ...)
3198 Free the space occupied by a NULL-terminated list of `mpq_t'
3201 -- Function: void mpq_set (mpq_t ROP, mpq_t OP)
3202 -- Function: void mpq_set_z (mpq_t ROP, mpz_t OP)
3205 -- Function: void mpq_set_ui (mpq_t ROP, unsigned long int OP1,
3206 unsigned long int OP2)
3207 -- Function: void mpq_set_si (mpq_t ROP, signed long int OP1, unsigned
3209 Set the value of ROP to OP1/OP2. Note that if OP1 and OP2 have
3210 common factors, ROP has to be passed to `mpq_canonicalize' before
3211 any operations are performed on ROP.
3213 -- Function: int mpq_set_str (mpq_t ROP, char *STR, int BASE)
3214 Set ROP from a null-terminated string STR in the given BASE.
3216 The string can be an integer like "41" or a fraction like
3217 "41/152". The fraction must be in canonical form (*note Rational
3218 Number Functions::), or if not then `mpq_canonicalize' must be
3221 The numerator and optional denominator are parsed the same as in
3222 `mpz_set_str' (*note Assigning Integers::). White space is
3223 allowed in the string, and is simply ignored. The BASE can vary
3224 from 2 to 62, or if BASE is 0 then the leading characters are
3225 used: `0x' or `0X' for hex, `0b' or `0B' for binary, `0' for
3226 octal, or decimal otherwise. Note that this is done separately
3227 for the numerator and denominator, so for instance `0xEF/100' is
3228 239/100, whereas `0xEF/0x100' is 239/256.
3230 The return value is 0 if the entire string is a valid number, or
3233 -- Function: void mpq_swap (mpq_t ROP1, mpq_t ROP2)
3234 Swap the values ROP1 and ROP2 efficiently.
3237 File: gmp.info, Node: Rational Conversions, Next: Rational Arithmetic, Prev: Initializing Rationals, Up: Rational Number Functions
3239 6.2 Conversion Functions
3240 ========================
3242 -- Function: double mpq_get_d (mpq_t OP)
3243 Convert OP to a `double', truncating if necessary (ie. rounding
3246 If the exponent from the conversion is too big or too small to fit
3247 a `double' then the result is system dependent. For too big an
3248 infinity is returned when available. For too small 0.0 is
3249 normally returned. Hardware overflow, underflow and denorm traps
3250 may or may not occur.
3252 -- Function: void mpq_set_d (mpq_t ROP, double OP)
3253 -- Function: void mpq_set_f (mpq_t ROP, mpf_t OP)
3254 Set ROP to the value of OP. There is no rounding, this conversion
3257 -- Function: char * mpq_get_str (char *STR, int BASE, mpq_t OP)
3258 Convert OP to a string of digits in base BASE. The base may vary
3259 from 2 to 36. The string will be of the form `num/den', or if the
3260 denominator is 1 then just `num'.
3262 If STR is `NULL', the result string is allocated using the current
3263 allocation function (*note Custom Allocation::). The block will be
3264 `strlen(str)+1' bytes, that being exactly enough for the string and
3267 If STR is not `NULL', it should point to a block of storage large
3268 enough for the result, that being
3270 mpz_sizeinbase (mpq_numref(OP), BASE)
3271 + mpz_sizeinbase (mpq_denref(OP), BASE) + 3
3273 The three extra bytes are for a possible minus sign, possible
3274 slash, and the null-terminator.
3276 A pointer to the result string is returned, being either the
3277 allocated block, or the given STR.
3280 File: gmp.info, Node: Rational Arithmetic, Next: Comparing Rationals, Prev: Rational Conversions, Up: Rational Number Functions
3282 6.3 Arithmetic Functions
3283 ========================
3285 -- Function: void mpq_add (mpq_t SUM, mpq_t ADDEND1, mpq_t ADDEND2)
3286 Set SUM to ADDEND1 + ADDEND2.
3288 -- Function: void mpq_sub (mpq_t DIFFERENCE, mpq_t MINUEND, mpq_t
3290 Set DIFFERENCE to MINUEND - SUBTRAHEND.
3292 -- Function: void mpq_mul (mpq_t PRODUCT, mpq_t MULTIPLIER, mpq_t
3294 Set PRODUCT to MULTIPLIER times MULTIPLICAND.
3296 -- Function: void mpq_mul_2exp (mpq_t ROP, mpq_t OP1, mp_bitcnt_t OP2)
3297 Set ROP to OP1 times 2 raised to OP2.
3299 -- Function: void mpq_div (mpq_t QUOTIENT, mpq_t DIVIDEND, mpq_t
3301 Set QUOTIENT to DIVIDEND/DIVISOR.
3303 -- Function: void mpq_div_2exp (mpq_t ROP, mpq_t OP1, mp_bitcnt_t OP2)
3304 Set ROP to OP1 divided by 2 raised to OP2.
3306 -- Function: void mpq_neg (mpq_t NEGATED_OPERAND, mpq_t OPERAND)
3307 Set NEGATED_OPERAND to -OPERAND.
3309 -- Function: void mpq_abs (mpq_t ROP, mpq_t OP)
3310 Set ROP to the absolute value of OP.
3312 -- Function: void mpq_inv (mpq_t INVERTED_NUMBER, mpq_t NUMBER)
3313 Set INVERTED_NUMBER to 1/NUMBER. If the new denominator is zero,
3314 this routine will divide by zero.
3317 File: gmp.info, Node: Comparing Rationals, Next: Applying Integer Functions, Prev: Rational Arithmetic, Up: Rational Number Functions
3319 6.4 Comparison Functions
3320 ========================
3322 -- Function: int mpq_cmp (mpq_t OP1, mpq_t OP2)
3323 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero
3324 if OP1 = OP2, and a negative value if OP1 < OP2.
3326 To determine if two rationals are equal, `mpq_equal' is faster than
3329 -- Macro: int mpq_cmp_ui (mpq_t OP1, unsigned long int NUM2, unsigned
3331 -- Macro: int mpq_cmp_si (mpq_t OP1, long int NUM2, unsigned long int
3333 Compare OP1 and NUM2/DEN2. Return a positive value if OP1 >
3334 NUM2/DEN2, zero if OP1 = NUM2/DEN2, and a negative value if OP1 <
3337 NUM2 and DEN2 are allowed to have common factors.
3339 These functions are implemented as a macros and evaluate their
3340 arguments multiple times.
3342 -- Macro: int mpq_sgn (mpq_t OP)
3343 Return +1 if OP > 0, 0 if OP = 0, and -1 if OP < 0.
3345 This function is actually implemented as a macro. It evaluates its
3346 arguments multiple times.
3348 -- Function: int mpq_equal (mpq_t OP1, mpq_t OP2)
3349 Return non-zero if OP1 and OP2 are equal, zero if they are
3350 non-equal. Although `mpq_cmp' can be used for the same purpose,
3351 this function is much faster.
3354 File: gmp.info, Node: Applying Integer Functions, Next: I/O of Rationals, Prev: Comparing Rationals, Up: Rational Number Functions
3356 6.5 Applying Integer Functions to Rationals
3357 ===========================================
3359 The set of `mpq' functions is quite small. In particular, there are few
3360 functions for either input or output. The following functions give
3361 direct access to the numerator and denominator of an `mpq_t'.
3363 Note that if an assignment to the numerator and/or denominator could
3364 take an `mpq_t' out of the canonical form described at the start of
3365 this chapter (*note Rational Number Functions::) then
3366 `mpq_canonicalize' must be called before any other `mpq' functions are
3367 applied to that `mpq_t'.
3369 -- Macro: mpz_t mpq_numref (mpq_t OP)
3370 -- Macro: mpz_t mpq_denref (mpq_t OP)
3371 Return a reference to the numerator and denominator of OP,
3372 respectively. The `mpz' functions can be used on the result of
3375 -- Function: void mpq_get_num (mpz_t NUMERATOR, mpq_t RATIONAL)
3376 -- Function: void mpq_get_den (mpz_t DENOMINATOR, mpq_t RATIONAL)
3377 -- Function: void mpq_set_num (mpq_t RATIONAL, mpz_t NUMERATOR)
3378 -- Function: void mpq_set_den (mpq_t RATIONAL, mpz_t DENOMINATOR)
3379 Get or set the numerator or denominator of a rational. These
3380 functions are equivalent to calling `mpz_set' with an appropriate
3381 `mpq_numref' or `mpq_denref'. Direct use of `mpq_numref' or
3382 `mpq_denref' is recommended instead of these functions.
3385 File: gmp.info, Node: I/O of Rationals, Prev: Applying Integer Functions, Up: Rational Number Functions
3387 6.6 Input and Output Functions
3388 ==============================
3390 When using any of these functions, it's a good idea to include `stdio.h'
3391 before `gmp.h', since that will allow `gmp.h' to define prototypes for
3394 Passing a `NULL' pointer for a STREAM argument to any of these
3395 functions will make them read from `stdin' and write to `stdout',
3398 -- Function: size_t mpq_out_str (FILE *STREAM, int BASE, mpq_t OP)
3399 Output OP on stdio stream STREAM, as a string of digits in base
3400 BASE. The base may vary from 2 to 36. Output is in the form
3401 `num/den' or if the denominator is 1 then just `num'.
3403 Return the number of bytes written, or if an error occurred,
3406 -- Function: size_t mpq_inp_str (mpq_t ROP, FILE *STREAM, int BASE)
3407 Read a string of digits from STREAM and convert them to a rational
3408 in ROP. Any initial white-space characters are read and
3409 discarded. Return the number of characters read (including white
3410 space), or 0 if a rational could not be read.
3412 The input can be a fraction like `17/63' or just an integer like
3413 `123'. Reading stops at the first character not in this form, and
3414 white space is not permitted within the string. If the input
3415 might not be in canonical form, then `mpq_canonicalize' must be
3416 called (*note Rational Number Functions::).
3418 The BASE can be between 2 and 36, or can be 0 in which case the
3419 leading characters of the string determine the base, `0x' or `0X'
3420 for hexadecimal, `0' for octal, or decimal otherwise. The leading
3421 characters are examined separately for the numerator and
3422 denominator of a fraction, so for instance `0x10/11' is 16/11,
3423 whereas `0x10/0x11' is 16/17.
3426 File: gmp.info, Node: Floating-point Functions, Next: Low-level Functions, Prev: Rational Number Functions, Up: Top
3428 7 Floating-point Functions
3429 **************************
3431 GMP floating point numbers are stored in objects of type `mpf_t' and
3432 functions operating on them have an `mpf_' prefix.
3434 The mantissa of each float has a user-selectable precision, limited
3435 only by available memory. Each variable has its own precision, and
3436 that can be increased or decreased at any time.
3438 The exponent of each float is a fixed precision, one machine word on
3439 most systems. In the current implementation the exponent is a count of
3440 limbs, so for example on a 32-bit system this means a range of roughly
3441 2^-68719476768 to 2^68719476736, or on a 64-bit system this will be
3442 greater. Note however `mpf_get_str' can only return an exponent which
3443 fits an `mp_exp_t' and currently `mpf_set_str' doesn't accept exponents
3444 bigger than a `long'.
3446 Each variable keeps a size for the mantissa data actually in use.
3447 This means that if a float is exactly represented in only a few bits
3448 then only those bits will be used in a calculation, even if the
3449 selected precision is high.
3451 All calculations are performed to the precision of the destination
3452 variable. Each function is defined to calculate with "infinite
3453 precision" followed by a truncation to the destination precision, but
3454 of course the work done is only what's needed to determine a result
3455 under that definition.
3457 The precision selected for a variable is a minimum value, GMP may
3458 increase it a little to facilitate efficient calculation. Currently
3459 this means rounding up to a whole limb, and then sometimes having a
3460 further partial limb, depending on the high limb of the mantissa. But
3461 applications shouldn't be concerned by such details.
3463 The mantissa in stored in binary, as might be imagined from the fact
3464 precisions are expressed in bits. One consequence of this is that
3465 decimal fractions like 0.1 cannot be represented exactly. The same is
3466 true of plain IEEE `double' floats. This makes both highly unsuitable
3467 for calculations involving money or other values that should be exact
3468 decimal fractions. (Suitably scaled integers, or perhaps rationals,
3469 are better choices.)
3471 `mpf' functions and variables have no special notion of infinity or
3472 not-a-number, and applications must take care not to overflow the
3473 exponent or results will be unpredictable. This might change in a
3476 Note that the `mpf' functions are _not_ intended as a smooth
3477 extension to IEEE P754 arithmetic. In particular results obtained on
3478 one computer often differ from the results on a computer with a
3479 different word size.
3483 * Initializing Floats::
3484 * Assigning Floats::
3485 * Simultaneous Float Init & Assign::
3486 * Converting Floats::
3487 * Float Arithmetic::
3488 * Float Comparison::
3490 * Miscellaneous Float Functions::
3493 File: gmp.info, Node: Initializing Floats, Next: Assigning Floats, Prev: Floating-point Functions, Up: Floating-point Functions
3495 7.1 Initialization Functions
3496 ============================
3498 -- Function: void mpf_set_default_prec (mp_bitcnt_t PREC)
3499 Set the default precision to be *at least* PREC bits. All
3500 subsequent calls to `mpf_init' will use this precision, but
3501 previously initialized variables are unaffected.
3503 -- Function: mp_bitcnt_t mpf_get_default_prec (void)
3504 Return the default precision actually used.
3506 An `mpf_t' object must be initialized before storing the first value
3507 in it. The functions `mpf_init' and `mpf_init2' are used for that
3510 -- Function: void mpf_init (mpf_t X)
3511 Initialize X to 0. Normally, a variable should be initialized
3512 once only or at least be cleared, using `mpf_clear', between
3513 initializations. The precision of X is undefined unless a default
3514 precision has already been established by a call to
3515 `mpf_set_default_prec'.
3517 -- Function: void mpf_init2 (mpf_t X, mp_bitcnt_t PREC)
3518 Initialize X to 0 and set its precision to be *at least* PREC
3519 bits. Normally, a variable should be initialized once only or at
3520 least be cleared, using `mpf_clear', between initializations.
3522 -- Function: void mpf_inits (mpf_t X, ...)
3523 Initialize a NULL-terminated list of `mpf_t' variables, and set
3524 their values to 0. The precision of the initialized variables is
3525 undefined unless a default precision has already been established
3526 by a call to `mpf_set_default_prec'.
3528 -- Function: void mpf_clear (mpf_t X)
3529 Free the space occupied by X. Make sure to call this function for
3530 all `mpf_t' variables when you are done with them.
3532 -- Function: void mpf_clears (mpf_t X, ...)
3533 Free the space occupied by a NULL-terminated list of `mpf_t'
3536 Here is an example on how to initialize floating-point variables:
3539 mpf_init (x); /* use default precision */
3540 mpf_init2 (y, 256); /* precision _at least_ 256 bits */
3542 /* Unless the program is about to exit, do ... */
3547 The following three functions are useful for changing the precision
3548 during a calculation. A typical use would be for adjusting the
3549 precision gradually in iterative algorithms like Newton-Raphson, making
3550 the computation precision closely match the actual accurate part of the
3553 -- Function: mp_bitcnt_t mpf_get_prec (mpf_t OP)
3554 Return the current precision of OP, in bits.
3556 -- Function: void mpf_set_prec (mpf_t ROP, mp_bitcnt_t PREC)
3557 Set the precision of ROP to be *at least* PREC bits. The value in
3558 ROP will be truncated to the new precision.
3560 This function requires a call to `realloc', and so should not be
3561 used in a tight loop.
3563 -- Function: void mpf_set_prec_raw (mpf_t ROP, mp_bitcnt_t PREC)
3564 Set the precision of ROP to be *at least* PREC bits, without
3565 changing the memory allocated.
3567 PREC must be no more than the allocated precision for ROP, that
3568 being the precision when ROP was initialized, or in the most recent
3571 The value in ROP is unchanged, and in particular if it had a higher
3572 precision than PREC it will retain that higher precision. New
3573 values written to ROP will use the new PREC.
3575 Before calling `mpf_clear' or the full `mpf_set_prec', another
3576 `mpf_set_prec_raw' call must be made to restore ROP to its original
3577 allocated precision. Failing to do so will have unpredictable
3580 `mpf_get_prec' can be used before `mpf_set_prec_raw' to get the
3581 original allocated precision. After `mpf_set_prec_raw' it
3582 reflects the PREC value set.
3584 `mpf_set_prec_raw' is an efficient way to use an `mpf_t' variable
3585 at different precisions during a calculation, perhaps to gradually
3586 increase precision in an iteration, or just to use various
3587 different precisions for different purposes during a calculation.
3590 File: gmp.info, Node: Assigning Floats, Next: Simultaneous Float Init & Assign, Prev: Initializing Floats, Up: Floating-point Functions
3592 7.2 Assignment Functions
3593 ========================
3595 These functions assign new values to already initialized floats (*note
3596 Initializing Floats::).
3598 -- Function: void mpf_set (mpf_t ROP, mpf_t OP)
3599 -- Function: void mpf_set_ui (mpf_t ROP, unsigned long int OP)
3600 -- Function: void mpf_set_si (mpf_t ROP, signed long int OP)
3601 -- Function: void mpf_set_d (mpf_t ROP, double OP)
3602 -- Function: void mpf_set_z (mpf_t ROP, mpz_t OP)
3603 -- Function: void mpf_set_q (mpf_t ROP, mpq_t OP)
3604 Set the value of ROP from OP.
3606 -- Function: int mpf_set_str (mpf_t ROP, char *STR, int BASE)
3607 Set the value of ROP from the string in STR. The string is of the
3608 form `M@N' or, if the base is 10 or less, alternatively `MeN'.
3609 `M' is the mantissa and `N' is the exponent. The mantissa is
3610 always in the specified base. The exponent is either in the
3611 specified base or, if BASE is negative, in decimal. The decimal
3612 point expected is taken from the current locale, on systems
3613 providing `localeconv'.
3615 The argument BASE may be in the ranges 2 to 62, or -62 to -2.
3616 Negative values are used to specify that the exponent is in
3619 For bases up to 36, case is ignored; upper-case and lower-case
3620 letters have the same value; for bases 37 to 62, upper-case letter
3621 represent the usual 10..35 while lower-case letter represent
3624 Unlike the corresponding `mpz' function, the base will not be
3625 determined from the leading characters of the string if BASE is 0.
3626 This is so that numbers like `0.23' are not interpreted as octal.
3628 White space is allowed in the string, and is simply ignored.
3629 [This is not really true; white-space is ignored in the beginning
3630 of the string and within the mantissa, but not in other places,
3631 such as after a minus sign or in the exponent. We are considering
3632 changing the definition of this function, making it fail when
3633 there is any white-space in the input, since that makes a lot of
3634 sense. Please tell us your opinion about this change. Do you
3635 really want it to accept "3 14" as meaning 314 as it does now?]
3637 This function returns 0 if the entire string is a valid number in
3638 base BASE. Otherwise it returns -1.
3640 -- Function: void mpf_swap (mpf_t ROP1, mpf_t ROP2)
3641 Swap ROP1 and ROP2 efficiently. Both the values and the
3642 precisions of the two variables are swapped.
3645 File: gmp.info, Node: Simultaneous Float Init & Assign, Next: Converting Floats, Prev: Assigning Floats, Up: Floating-point Functions
3647 7.3 Combined Initialization and Assignment Functions
3648 ====================================================
3650 For convenience, GMP provides a parallel series of initialize-and-set
3651 functions which initialize the output and then store the value there.
3652 These functions' names have the form `mpf_init_set...'
3654 Once the float has been initialized by any of the `mpf_init_set...'
3655 functions, it can be used as the source or destination operand for the
3656 ordinary float functions. Don't use an initialize-and-set function on
3657 a variable already initialized!
3659 -- Function: void mpf_init_set (mpf_t ROP, mpf_t OP)
3660 -- Function: void mpf_init_set_ui (mpf_t ROP, unsigned long int OP)
3661 -- Function: void mpf_init_set_si (mpf_t ROP, signed long int OP)
3662 -- Function: void mpf_init_set_d (mpf_t ROP, double OP)
3663 Initialize ROP and set its value from OP.
3665 The precision of ROP will be taken from the active default
3666 precision, as set by `mpf_set_default_prec'.
3668 -- Function: int mpf_init_set_str (mpf_t ROP, char *STR, int BASE)
3669 Initialize ROP and set its value from the string in STR. See
3670 `mpf_set_str' above for details on the assignment operation.
3672 Note that ROP is initialized even if an error occurs. (I.e., you
3673 have to call `mpf_clear' for it.)
3675 The precision of ROP will be taken from the active default
3676 precision, as set by `mpf_set_default_prec'.
3679 File: gmp.info, Node: Converting Floats, Next: Float Arithmetic, Prev: Simultaneous Float Init & Assign, Up: Floating-point Functions
3681 7.4 Conversion Functions
3682 ========================
3684 -- Function: double mpf_get_d (mpf_t OP)
3685 Convert OP to a `double', truncating if necessary (ie. rounding
3688 If the exponent in OP is too big or too small to fit a `double'
3689 then the result is system dependent. For too big an infinity is
3690 returned when available. For too small 0.0 is normally returned.
3691 Hardware overflow, underflow and denorm traps may or may not occur.
3693 -- Function: double mpf_get_d_2exp (signed long int *EXP, mpf_t OP)
3694 Convert OP to a `double', truncating if necessary (ie. rounding
3695 towards zero), and with an exponent returned separately.
3697 The return value is in the range 0.5<=abs(D)<1 and the exponent is
3698 stored to `*EXP'. D * 2^EXP is the (truncated) OP value. If OP
3699 is zero, the return is 0.0 and 0 is stored to `*EXP'.
3701 This is similar to the standard C `frexp' function (*note
3702 Normalization Functions: (libc)Normalization Functions.).
3704 -- Function: long mpf_get_si (mpf_t OP)
3705 -- Function: unsigned long mpf_get_ui (mpf_t OP)
3706 Convert OP to a `long' or `unsigned long', truncating any fraction
3707 part. If OP is too big for the return type, the result is
3710 See also `mpf_fits_slong_p' and `mpf_fits_ulong_p' (*note
3711 Miscellaneous Float Functions::).
3713 -- Function: char * mpf_get_str (char *STR, mp_exp_t *EXPPTR, int
3714 BASE, size_t N_DIGITS, mpf_t OP)
3715 Convert OP to a string of digits in base BASE. The base argument
3716 may vary from 2 to 62 or from -2 to -36. Up to N_DIGITS digits
3717 will be generated. Trailing zeros are not returned. No more
3718 digits than can be accurately represented by OP are ever
3719 generated. If N_DIGITS is 0 then that accurate maximum number of
3720 digits are generated.
3722 For BASE in the range 2..36, digits and lower-case letters are
3723 used; for -2..-36, digits and upper-case letters are used; for
3724 37..62, digits, upper-case letters, and lower-case letters (in
3725 that significance order) are used.
3727 If STR is `NULL', the result string is allocated using the current
3728 allocation function (*note Custom Allocation::). The block will be
3729 `strlen(str)+1' bytes, that being exactly enough for the string and
3732 If STR is not `NULL', it should point to a block of N_DIGITS + 2
3733 bytes, that being enough for the mantissa, a possible minus sign,
3734 and a null-terminator. When N_DIGITS is 0 to get all significant
3735 digits, an application won't be able to know the space required,
3736 and STR should be `NULL' in that case.
3738 The generated string is a fraction, with an implicit radix point
3739 immediately to the left of the first digit. The applicable
3740 exponent is written through the EXPPTR pointer. For example, the
3741 number 3.1416 would be returned as string "31416" and exponent 1.
3743 When OP is zero, an empty string is produced and the exponent
3746 A pointer to the result string is returned, being either the
3747 allocated block or the given STR.
3750 File: gmp.info, Node: Float Arithmetic, Next: Float Comparison, Prev: Converting Floats, Up: Floating-point Functions
3752 7.5 Arithmetic Functions
3753 ========================
3755 -- Function: void mpf_add (mpf_t ROP, mpf_t OP1, mpf_t OP2)
3756 -- Function: void mpf_add_ui (mpf_t ROP, mpf_t OP1, unsigned long int
3758 Set ROP to OP1 + OP2.
3760 -- Function: void mpf_sub (mpf_t ROP, mpf_t OP1, mpf_t OP2)
3761 -- Function: void mpf_ui_sub (mpf_t ROP, unsigned long int OP1, mpf_t
3763 -- Function: void mpf_sub_ui (mpf_t ROP, mpf_t OP1, unsigned long int
3765 Set ROP to OP1 - OP2.
3767 -- Function: void mpf_mul (mpf_t ROP, mpf_t OP1, mpf_t OP2)
3768 -- Function: void mpf_mul_ui (mpf_t ROP, mpf_t OP1, unsigned long int
3770 Set ROP to OP1 times OP2.
3772 Division is undefined if the divisor is zero, and passing a zero
3773 divisor to the divide functions will make these functions intentionally
3774 divide by zero. This lets the user handle arithmetic exceptions in
3775 these functions in the same manner as other arithmetic exceptions.
3777 -- Function: void mpf_div (mpf_t ROP, mpf_t OP1, mpf_t OP2)
3778 -- Function: void mpf_ui_div (mpf_t ROP, unsigned long int OP1, mpf_t
3780 -- Function: void mpf_div_ui (mpf_t ROP, mpf_t OP1, unsigned long int
3784 -- Function: void mpf_sqrt (mpf_t ROP, mpf_t OP)
3785 -- Function: void mpf_sqrt_ui (mpf_t ROP, unsigned long int OP)
3786 Set ROP to the square root of OP.
3788 -- Function: void mpf_pow_ui (mpf_t ROP, mpf_t OP1, unsigned long int
3790 Set ROP to OP1 raised to the power OP2.
3792 -- Function: void mpf_neg (mpf_t ROP, mpf_t OP)
3795 -- Function: void mpf_abs (mpf_t ROP, mpf_t OP)
3796 Set ROP to the absolute value of OP.
3798 -- Function: void mpf_mul_2exp (mpf_t ROP, mpf_t OP1, mp_bitcnt_t OP2)
3799 Set ROP to OP1 times 2 raised to OP2.
3801 -- Function: void mpf_div_2exp (mpf_t ROP, mpf_t OP1, mp_bitcnt_t OP2)
3802 Set ROP to OP1 divided by 2 raised to OP2.
3805 File: gmp.info, Node: Float Comparison, Next: I/O of Floats, Prev: Float Arithmetic, Up: Floating-point Functions
3807 7.6 Comparison Functions
3808 ========================
3810 -- Function: int mpf_cmp (mpf_t OP1, mpf_t OP2)
3811 -- Function: int mpf_cmp_d (mpf_t OP1, double OP2)
3812 -- Function: int mpf_cmp_ui (mpf_t OP1, unsigned long int OP2)
3813 -- Function: int mpf_cmp_si (mpf_t OP1, signed long int OP2)
3814 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero
3815 if OP1 = OP2, and a negative value if OP1 < OP2.
3817 `mpf_cmp_d' can be called with an infinity, but results are
3818 undefined for a NaN.
3820 -- Function: int mpf_eq (mpf_t OP1, mpf_t OP2, mp_bitcnt_t op3)
3821 Return non-zero if the first OP3 bits of OP1 and OP2 are equal,
3822 zero otherwise. I.e., test if OP1 and OP2 are approximately equal.
3824 Caution 1: All version of GMP up to version 4.2.4 compared just
3825 whole limbs, meaning sometimes more than OP3 bits, sometimes fewer.
3827 Caution 2: This function will consider XXX11...111 and XX100...000
3828 different, even if ... is replaced by a semi-infinite number of
3829 bits. Such numbers are really just one ulp off, and should be
3832 -- Function: void mpf_reldiff (mpf_t ROP, mpf_t OP1, mpf_t OP2)
3833 Compute the relative difference between OP1 and OP2 and store the
3834 result in ROP. This is abs(OP1-OP2)/OP1.
3836 -- Macro: int mpf_sgn (mpf_t OP)
3837 Return +1 if OP > 0, 0 if OP = 0, and -1 if OP < 0.
3839 This function is actually implemented as a macro. It evaluates
3840 its arguments multiple times.
3843 File: gmp.info, Node: I/O of Floats, Next: Miscellaneous Float Functions, Prev: Float Comparison, Up: Floating-point Functions
3845 7.7 Input and Output Functions
3846 ==============================
3848 Functions that perform input from a stdio stream, and functions that
3849 output to a stdio stream. Passing a `NULL' pointer for a STREAM
3850 argument to any of these functions will make them read from `stdin' and
3851 write to `stdout', respectively.
3853 When using any of these functions, it is a good idea to include
3854 `stdio.h' before `gmp.h', since that will allow `gmp.h' to define
3855 prototypes for these functions.
3857 -- Function: size_t mpf_out_str (FILE *STREAM, int BASE, size_t
3859 Print OP to STREAM, as a string of digits. Return the number of
3860 bytes written, or if an error occurred, return 0.
3862 The mantissa is prefixed with an `0.' and is in the given BASE,
3863 which may vary from 2 to 62 or from -2 to -36. An exponent is
3864 then printed, separated by an `e', or if the base is greater than
3865 10 then by an `@'. The exponent is always in decimal. The
3866 decimal point follows the current locale, on systems providing
3869 For BASE in the range 2..36, digits and lower-case letters are
3870 used; for -2..-36, digits and upper-case letters are used; for
3871 37..62, digits, upper-case letters, and lower-case letters (in
3872 that significance order) are used.
3874 Up to N_DIGITS will be printed from the mantissa, except that no
3875 more digits than are accurately representable by OP will be
3876 printed. N_DIGITS can be 0 to select that accurate maximum.
3878 -- Function: size_t mpf_inp_str (mpf_t ROP, FILE *STREAM, int BASE)
3879 Read a string in base BASE from STREAM, and put the read float in
3880 ROP. The string is of the form `M@N' or, if the base is 10 or
3881 less, alternatively `MeN'. `M' is the mantissa and `N' is the
3882 exponent. The mantissa is always in the specified base. The
3883 exponent is either in the specified base or, if BASE is negative,
3884 in decimal. The decimal point expected is taken from the current
3885 locale, on systems providing `localeconv'.
3887 The argument BASE may be in the ranges 2 to 36, or -36 to -2.
3888 Negative values are used to specify that the exponent is in
3891 Unlike the corresponding `mpz' function, the base will not be
3892 determined from the leading characters of the string if BASE is 0.
3893 This is so that numbers like `0.23' are not interpreted as octal.
3895 Return the number of bytes read, or if an error occurred, return 0.
3898 File: gmp.info, Node: Miscellaneous Float Functions, Prev: I/O of Floats, Up: Floating-point Functions
3900 7.8 Miscellaneous Functions
3901 ===========================
3903 -- Function: void mpf_ceil (mpf_t ROP, mpf_t OP)
3904 -- Function: void mpf_floor (mpf_t ROP, mpf_t OP)
3905 -- Function: void mpf_trunc (mpf_t ROP, mpf_t OP)
3906 Set ROP to OP rounded to an integer. `mpf_ceil' rounds to the
3907 next higher integer, `mpf_floor' to the next lower, and `mpf_trunc'
3908 to the integer towards zero.
3910 -- Function: int mpf_integer_p (mpf_t OP)
3911 Return non-zero if OP is an integer.
3913 -- Function: int mpf_fits_ulong_p (mpf_t OP)
3914 -- Function: int mpf_fits_slong_p (mpf_t OP)
3915 -- Function: int mpf_fits_uint_p (mpf_t OP)
3916 -- Function: int mpf_fits_sint_p (mpf_t OP)
3917 -- Function: int mpf_fits_ushort_p (mpf_t OP)
3918 -- Function: int mpf_fits_sshort_p (mpf_t OP)
3919 Return non-zero if OP would fit in the respective C data type, when
3920 truncated to an integer.
3922 -- Function: void mpf_urandomb (mpf_t ROP, gmp_randstate_t STATE,
3924 Generate a uniformly distributed random float in ROP, such that 0
3925 <= ROP < 1, with NBITS significant bits in the mantissa.
3927 The variable STATE must be initialized by calling one of the
3928 `gmp_randinit' functions (*Note Random State Initialization::)
3929 before invoking this function.
3931 -- Function: void mpf_random2 (mpf_t ROP, mp_size_t MAX_SIZE, mp_exp_t
3933 Generate a random float of at most MAX_SIZE limbs, with long
3934 strings of zeros and ones in the binary representation. The
3935 exponent of the number is in the interval -EXP to EXP (in limbs).
3936 This function is useful for testing functions and algorithms,
3937 since these kind of random numbers have proven to be more likely
3938 to trigger corner-case bugs. Negative random numbers are
3939 generated when MAX_SIZE is negative.
3942 File: gmp.info, Node: Low-level Functions, Next: Random Number Functions, Prev: Floating-point Functions, Up: Top
3944 8 Low-level Functions
3945 *********************
3947 This chapter describes low-level GMP functions, used to implement the
3948 high-level GMP functions, but also intended for time-critical user code.
3950 These functions start with the prefix `mpn_'.
3952 The `mpn' functions are designed to be as fast as possible, *not* to
3953 provide a coherent calling interface. The different functions have
3954 somewhat similar interfaces, but there are variations that make them
3955 hard to use. These functions do as little as possible apart from the
3956 real multiple precision computation, so that no time is spent on things
3957 that not all callers need.
3959 A source operand is specified by a pointer to the least significant
3960 limb and a limb count. A destination operand is specified by just a
3961 pointer. It is the responsibility of the caller to ensure that the
3962 destination has enough space for storing the result.
3964 With this way of specifying operands, it is possible to perform
3965 computations on subranges of an argument, and store the result into a
3966 subrange of a destination.
3968 A common requirement for all functions is that each source area
3969 needs at least one limb. No size argument may be zero. Unless
3970 otherwise stated, in-place operations are allowed where source and
3971 destination are the same, but not where they only partly overlap.
3973 The `mpn' functions are the base for the implementation of the
3974 `mpz_', `mpf_', and `mpq_' functions.
3976 This example adds the number beginning at S1P and the number
3977 beginning at S2P and writes the sum at DESTP. All areas have N limbs.
3979 cy = mpn_add_n (destp, s1p, s2p, n)
3981 It should be noted that the `mpn' functions make no attempt to
3982 identify high or low zero limbs on their operands, or other special
3983 forms. On random data such cases will be unlikely and it'd be wasteful
3984 for every function to check every time. An application knowing
3985 something about its data can take steps to trim or perhaps split its
3989 In the notation used below, a source operand is identified by the
3990 pointer to the least significant limb, and the limb count in braces.
3991 For example, {S1P, S1N}.
3993 -- Function: mp_limb_t mpn_add_n (mp_limb_t *RP, const mp_limb_t *S1P,
3994 const mp_limb_t *S2P, mp_size_t N)
3995 Add {S1P, N} and {S2P, N}, and write the N least significant limbs
3996 of the result to RP. Return carry, either 0 or 1.
3998 This is the lowest-level function for addition. It is the
3999 preferred function for addition, since it is written in assembly
4000 for most CPUs. For addition of a variable to itself (i.e., S1P
4001 equals S2P) use `mpn_lshift' with a count of 1 for optimal speed.
4003 -- Function: mp_limb_t mpn_add_1 (mp_limb_t *RP, const mp_limb_t *S1P,
4004 mp_size_t N, mp_limb_t S2LIMB)
4005 Add {S1P, N} and S2LIMB, and write the N least significant limbs
4006 of the result to RP. Return carry, either 0 or 1.
4008 -- Function: mp_limb_t mpn_add (mp_limb_t *RP, const mp_limb_t *S1P,
4009 mp_size_t S1N, const mp_limb_t *S2P, mp_size_t S2N)
4010 Add {S1P, S1N} and {S2P, S2N}, and write the S1N least significant
4011 limbs of the result to RP. Return carry, either 0 or 1.
4013 This function requires that S1N is greater than or equal to S2N.
4015 -- Function: mp_limb_t mpn_sub_n (mp_limb_t *RP, const mp_limb_t *S1P,
4016 const mp_limb_t *S2P, mp_size_t N)
4017 Subtract {S2P, N} from {S1P, N}, and write the N least significant
4018 limbs of the result to RP. Return borrow, either 0 or 1.
4020 This is the lowest-level function for subtraction. It is the
4021 preferred function for subtraction, since it is written in
4022 assembly for most CPUs.
4024 -- Function: mp_limb_t mpn_sub_1 (mp_limb_t *RP, const mp_limb_t *S1P,
4025 mp_size_t N, mp_limb_t S2LIMB)
4026 Subtract S2LIMB from {S1P, N}, and write the N least significant
4027 limbs of the result to RP. Return borrow, either 0 or 1.
4029 -- Function: mp_limb_t mpn_sub (mp_limb_t *RP, const mp_limb_t *S1P,
4030 mp_size_t S1N, const mp_limb_t *S2P, mp_size_t S2N)
4031 Subtract {S2P, S2N} from {S1P, S1N}, and write the S1N least
4032 significant limbs of the result to RP. Return borrow, either 0 or
4035 This function requires that S1N is greater than or equal to S2N.
4037 -- Function: void mpn_neg (mp_limb_t *RP, const mp_limb_t *SP,
4039 Perform the negation of {SP, N}, and write the result to {RP, N}.
4042 -- Function: void mpn_mul_n (mp_limb_t *RP, const mp_limb_t *S1P,
4043 const mp_limb_t *S2P, mp_size_t N)
4044 Multiply {S1P, N} and {S2P, N}, and write the 2*N-limb result to
4047 The destination has to have space for 2*N limbs, even if the
4048 product's most significant limb is zero. No overlap is permitted
4049 between the destination and either source.
4051 If the two input operands are the same, use `mpn_sqr'.
4053 -- Function: mp_limb_t mpn_mul (mp_limb_t *RP, const mp_limb_t *S1P,
4054 mp_size_t S1N, const mp_limb_t *S2P, mp_size_t S2N)
4055 Multiply {S1P, S1N} and {S2P, S2N}, and write the (S1N+S2N)-limb
4056 result to RP. Return the most significant limb of the result.
4058 The destination has to have space for S1N + S2N limbs, even if the
4059 product's most significant limb is zero. No overlap is permitted
4060 between the destination and either source.
4062 This function requires that S1N is greater than or equal to S2N.
4064 -- Function: void mpn_sqr (mp_limb_t *RP, const mp_limb_t *S1P,
4066 Compute the square of {S1P, N} and write the 2*N-limb result to RP.
4068 The destination has to have space for 2*N limbs, even if the
4069 result's most significant limb is zero. No overlap is permitted
4070 between the destination and the source.
4072 -- Function: mp_limb_t mpn_mul_1 (mp_limb_t *RP, const mp_limb_t *S1P,
4073 mp_size_t N, mp_limb_t S2LIMB)
4074 Multiply {S1P, N} by S2LIMB, and write the N least significant
4075 limbs of the product to RP. Return the most significant limb of
4076 the product. {S1P, N} and {RP, N} are allowed to overlap provided
4079 This is a low-level function that is a building block for general
4080 multiplication as well as other operations in GMP. It is written
4081 in assembly for most CPUs.
4083 Don't call this function if S2LIMB is a power of 2; use
4084 `mpn_lshift' with a count equal to the logarithm of S2LIMB
4085 instead, for optimal speed.
4087 -- Function: mp_limb_t mpn_addmul_1 (mp_limb_t *RP, const mp_limb_t
4088 *S1P, mp_size_t N, mp_limb_t S2LIMB)
4089 Multiply {S1P, N} and S2LIMB, and add the N least significant
4090 limbs of the product to {RP, N} and write the result to RP.
4091 Return the most significant limb of the product, plus carry-out
4094 This is a low-level function that is a building block for general
4095 multiplication as well as other operations in GMP. It is written
4096 in assembly for most CPUs.
4098 -- Function: mp_limb_t mpn_submul_1 (mp_limb_t *RP, const mp_limb_t
4099 *S1P, mp_size_t N, mp_limb_t S2LIMB)
4100 Multiply {S1P, N} and S2LIMB, and subtract the N least significant
4101 limbs of the product from {RP, N} and write the result to RP.
4102 Return the most significant limb of the product, plus borrow-out
4103 from the subtraction.
4105 This is a low-level function that is a building block for general
4106 multiplication and division as well as other operations in GMP.
4107 It is written in assembly for most CPUs.
4109 -- Function: void mpn_tdiv_qr (mp_limb_t *QP, mp_limb_t *RP, mp_size_t
4110 QXN, const mp_limb_t *NP, mp_size_t NN, const mp_limb_t *DP,
4112 Divide {NP, NN} by {DP, DN} and put the quotient at {QP, NN-DN+1}
4113 and the remainder at {RP, DN}. The quotient is rounded towards 0.
4115 No overlap is permitted between arguments, except that NP might
4116 equal RP. The dividend size NN must be greater than or equal to
4117 divisor size DN. The most significant limb of the divisor must be
4118 non-zero. The QXN operand must be zero.
4120 -- Function: mp_limb_t mpn_divrem (mp_limb_t *R1P, mp_size_t QXN,
4121 mp_limb_t *RS2P, mp_size_t RS2N, const mp_limb_t *S3P,
4123 [This function is obsolete. Please call `mpn_tdiv_qr' instead for
4126 Divide {RS2P, RS2N} by {S3P, S3N}, and write the quotient at R1P,
4127 with the exception of the most significant limb, which is
4128 returned. The remainder replaces the dividend at RS2P; it will be
4129 S3N limbs long (i.e., as many limbs as the divisor).
4131 In addition to an integer quotient, QXN fraction limbs are
4132 developed, and stored after the integral limbs. For most usages,
4135 It is required that RS2N is greater than or equal to S3N. It is
4136 required that the most significant bit of the divisor is set.
4138 If the quotient is not needed, pass RS2P + S3N as R1P. Aside from
4139 that special case, no overlap between arguments is permitted.
4141 Return the most significant limb of the quotient, either 0 or 1.
4143 The area at R1P needs to be RS2N - S3N + QXN limbs large.
4145 -- Function: mp_limb_t mpn_divrem_1 (mp_limb_t *R1P, mp_size_t QXN,
4146 mp_limb_t *S2P, mp_size_t S2N, mp_limb_t S3LIMB)
4147 -- Macro: mp_limb_t mpn_divmod_1 (mp_limb_t *R1P, mp_limb_t *S2P,
4148 mp_size_t S2N, mp_limb_t S3LIMB)
4149 Divide {S2P, S2N} by S3LIMB, and write the quotient at R1P.
4150 Return the remainder.
4152 The integer quotient is written to {R1P+QXN, S2N} and in addition
4153 QXN fraction limbs are developed and written to {R1P, QXN}.
4154 Either or both S2N and QXN can be zero. For most usages, QXN will
4157 `mpn_divmod_1' exists for upward source compatibility and is
4158 simply a macro calling `mpn_divrem_1' with a QXN of 0.
4160 The areas at R1P and S2P have to be identical or completely
4161 separate, not partially overlapping.
4163 -- Function: mp_limb_t mpn_divmod (mp_limb_t *R1P, mp_limb_t *RS2P,
4164 mp_size_t RS2N, const mp_limb_t *S3P, mp_size_t S3N)
4165 [This function is obsolete. Please call `mpn_tdiv_qr' instead for
4168 -- Macro: mp_limb_t mpn_divexact_by3 (mp_limb_t *RP, mp_limb_t *SP,
4170 -- Function: mp_limb_t mpn_divexact_by3c (mp_limb_t *RP, mp_limb_t
4171 *SP, mp_size_t N, mp_limb_t CARRY)
4172 Divide {SP, N} by 3, expecting it to divide exactly, and writing
4173 the result to {RP, N}. If 3 divides exactly, the return value is
4174 zero and the result is the quotient. If not, the return value is
4175 non-zero and the result won't be anything useful.
4177 `mpn_divexact_by3c' takes an initial carry parameter, which can be
4178 the return value from a previous call, so a large calculation can
4179 be done piece by piece from low to high. `mpn_divexact_by3' is
4180 simply a macro calling `mpn_divexact_by3c' with a 0 carry
4183 These routines use a multiply-by-inverse and will be faster than
4184 `mpn_divrem_1' on CPUs with fast multiplication but slow division.
4186 The source a, result q, size n, initial carry i, and return value
4187 c satisfy c*b^n + a-i = 3*q, where b=2^GMP_NUMB_BITS. The return
4188 c is always 0, 1 or 2, and the initial carry i must also be 0, 1
4189 or 2 (these are both borrows really). When c=0 clearly q=(a-i)/3.
4190 When c!=0, the remainder (a-i) mod 3 is given by 3-c, because b
4191 == 1 mod 3 (when `mp_bits_per_limb' is even, which is always so
4194 -- Function: mp_limb_t mpn_mod_1 (mp_limb_t *S1P, mp_size_t S1N,
4196 Divide {S1P, S1N} by S2LIMB, and return the remainder. S1N can be
4199 -- Function: mp_limb_t mpn_lshift (mp_limb_t *RP, const mp_limb_t *SP,
4200 mp_size_t N, unsigned int COUNT)
4201 Shift {SP, N} left by COUNT bits, and write the result to {RP, N}.
4202 The bits shifted out at the left are returned in the least
4203 significant COUNT bits of the return value (the rest of the return
4206 COUNT must be in the range 1 to mp_bits_per_limb-1. The regions
4207 {SP, N} and {RP, N} may overlap, provided RP >= SP.
4209 This function is written in assembly for most CPUs.
4211 -- Function: mp_limb_t mpn_rshift (mp_limb_t *RP, const mp_limb_t *SP,
4212 mp_size_t N, unsigned int COUNT)
4213 Shift {SP, N} right by COUNT bits, and write the result to {RP,
4214 N}. The bits shifted out at the right are returned in the most
4215 significant COUNT bits of the return value (the rest of the return
4218 COUNT must be in the range 1 to mp_bits_per_limb-1. The regions
4219 {SP, N} and {RP, N} may overlap, provided RP <= SP.
4221 This function is written in assembly for most CPUs.
4223 -- Function: int mpn_cmp (const mp_limb_t *S1P, const mp_limb_t *S2P,
4225 Compare {S1P, N} and {S2P, N} and return a positive value if S1 >
4226 S2, 0 if they are equal, or a negative value if S1 < S2.
4228 -- Function: mp_size_t mpn_gcd (mp_limb_t *RP, mp_limb_t *XP,
4229 mp_size_t XN, mp_limb_t *YP, mp_size_t YN)
4230 Set {RP, RETVAL} to the greatest common divisor of {XP, XN} and
4231 {YP, YN}. The result can be up to YN limbs, the return value is
4232 the actual number produced. Both source operands are destroyed.
4234 {XP, XN} must have at least as many bits as {YP, YN}. {YP, YN}
4235 must be odd. Both operands must have non-zero most significant
4236 limbs. No overlap is permitted between {XP, XN} and {YP, YN}.
4238 -- Function: mp_limb_t mpn_gcd_1 (const mp_limb_t *XP, mp_size_t XN,
4240 Return the greatest common divisor of {XP, XN} and YLIMB. Both
4241 operands must be non-zero.
4243 -- Function: mp_size_t mpn_gcdext (mp_limb_t *GP, mp_limb_t *SP,
4244 mp_size_t *SN, mp_limb_t *XP, mp_size_t XN, mp_limb_t *YP,
4246 Let U be defined by {XP, XN} and let V be defined by {YP, YN}.
4248 Compute the greatest common divisor G of U and V. Compute a
4249 cofactor S such that G = US + VT. The second cofactor T is not
4250 computed but can easily be obtained from (G - U*S) / V (the
4251 division will be exact). It is required that U >= V > 0.
4253 S satisfies S = 1 or abs(S) < V / (2 G). S = 0 if and only if V
4254 divides U (i.e., G = V).
4256 Store G at GP and let the return value define its limb count.
4257 Store S at SP and let |*SN| define its limb count. S can be
4258 negative; when this happens *SN will be negative. The areas at GP
4259 and SP should each have room for XN+1 limbs.
4261 The areas {XP, XN+1} and {YP, YN+1} are destroyed (i.e. the input
4262 operands plus an extra limb past the end of each).
4264 Compatibility note: GMP 4.3.0 and 4.3.1 defined S less strictly.
4265 Earlier as well as later GMP releases define S as described here.
4267 -- Function: mp_size_t mpn_sqrtrem (mp_limb_t *R1P, mp_limb_t *R2P,
4268 const mp_limb_t *SP, mp_size_t N)
4269 Compute the square root of {SP, N} and put the result at {R1P,
4270 ceil(N/2)} and the remainder at {R2P, RETVAL}. R2P needs space
4271 for N limbs, but the return value indicates how many are produced.
4273 The most significant limb of {SP, N} must be non-zero. The areas
4274 {R1P, ceil(N/2)} and {SP, N} must be completely separate. The
4275 areas {R2P, N} and {SP, N} must be either identical or completely
4278 If the remainder is not wanted then R2P can be `NULL', and in this
4279 case the return value is zero or non-zero according to whether the
4280 remainder would have been zero or non-zero.
4282 A return value of zero indicates a perfect square. See also
4283 `mpz_perfect_square_p'.
4285 -- Function: mp_size_t mpn_get_str (unsigned char *STR, int BASE,
4286 mp_limb_t *S1P, mp_size_t S1N)
4287 Convert {S1P, S1N} to a raw unsigned char array at STR in base
4288 BASE, and return the number of characters produced. There may be
4289 leading zeros in the string. The string is not in ASCII; to
4290 convert it to printable format, add the ASCII codes for `0' or
4291 `A', depending on the base and range. BASE can vary from 2 to 256.
4293 The most significant limb of the input {S1P, S1N} must be
4294 non-zero. The input {S1P, S1N} is clobbered, except when BASE is
4295 a power of 2, in which case it's unchanged.
4297 The area at STR has to have space for the largest possible number
4298 represented by a S1N long limb array, plus one extra character.
4300 -- Function: mp_size_t mpn_set_str (mp_limb_t *RP, const unsigned char
4301 *STR, size_t STRSIZE, int BASE)
4302 Convert bytes {STR,STRSIZE} in the given BASE to limbs at RP.
4304 STR[0] is the most significant byte and STR[STRSIZE-1] is the
4305 least significant. Each byte should be a value in the range 0 to
4306 BASE-1, not an ASCII character. BASE can vary from 2 to 256.
4308 The return value is the number of limbs written to RP. If the most
4309 significant input byte is non-zero then the high limb at RP will be
4310 non-zero, and only that exact number of limbs will be required
4313 If the most significant input byte is zero then there may be high
4314 zero limbs written to RP and included in the return value.
4316 STRSIZE must be at least 1, and no overlap is permitted between
4317 {STR,STRSIZE} and the result at RP.
4319 -- Function: mp_bitcnt_t mpn_scan0 (const mp_limb_t *S1P, mp_bitcnt_t
4321 Scan S1P from bit position BIT for the next clear bit.
4323 It is required that there be a clear bit within the area at S1P at
4324 or beyond bit position BIT, so that the function has something to
4327 -- Function: mp_bitcnt_t mpn_scan1 (const mp_limb_t *S1P, mp_bitcnt_t
4329 Scan S1P from bit position BIT for the next set bit.
4331 It is required that there be a set bit within the area at S1P at or
4332 beyond bit position BIT, so that the function has something to
4335 -- Function: void mpn_random (mp_limb_t *R1P, mp_size_t R1N)
4336 -- Function: void mpn_random2 (mp_limb_t *R1P, mp_size_t R1N)
4337 Generate a random number of length R1N and store it at R1P. The
4338 most significant limb is always non-zero. `mpn_random' generates
4339 uniformly distributed limb data, `mpn_random2' generates long
4340 strings of zeros and ones in the binary representation.
4342 `mpn_random2' is intended for testing the correctness of the `mpn'
4345 -- Function: mp_bitcnt_t mpn_popcount (const mp_limb_t *S1P, mp_size_t
4347 Count the number of set bits in {S1P, N}.
4349 -- Function: mp_bitcnt_t mpn_hamdist (const mp_limb_t *S1P, const
4350 mp_limb_t *S2P, mp_size_t N)
4351 Compute the hamming distance between {S1P, N} and {S2P, N}, which
4352 is the number of bit positions where the two operands have
4353 different bit values.
4355 -- Function: int mpn_perfect_square_p (const mp_limb_t *S1P, mp_size_t
4357 Return non-zero iff {S1P, N} is a perfect square.
4359 -- Function: void mpn_and_n (mp_limb_t *RP, const mp_limb_t *S1P,
4360 const mp_limb_t *S2P, mp_size_t N)
4361 Perform the bitwise logical and of {S1P, N} and {S2P, N}, and
4362 write the result to {RP, N}.
4364 -- Function: void mpn_ior_n (mp_limb_t *RP, const mp_limb_t *S1P,
4365 const mp_limb_t *S2P, mp_size_t N)
4366 Perform the bitwise logical inclusive or of {S1P, N} and {S2P, N},
4367 and write the result to {RP, N}.
4369 -- Function: void mpn_xor_n (mp_limb_t *RP, const mp_limb_t *S1P,
4370 const mp_limb_t *S2P, mp_size_t N)
4371 Perform the bitwise logical exclusive or of {S1P, N} and {S2P, N},
4372 and write the result to {RP, N}.
4374 -- Function: void mpn_andn_n (mp_limb_t *RP, const mp_limb_t *S1P,
4375 const mp_limb_t *S2P, mp_size_t N)
4376 Perform the bitwise logical and of {S1P, N} and the bitwise
4377 complement of {S2P, N}, and write the result to {RP, N}.
4379 -- Function: void mpn_iorn_n (mp_limb_t *RP, const mp_limb_t *S1P,
4380 const mp_limb_t *S2P, mp_size_t N)
4381 Perform the bitwise logical inclusive or of {S1P, N} and the
4382 bitwise complement of {S2P, N}, and write the result to {RP, N}.
4384 -- Function: void mpn_nand_n (mp_limb_t *RP, const mp_limb_t *S1P,
4385 const mp_limb_t *S2P, mp_size_t N)
4386 Perform the bitwise logical and of {S1P, N} and {S2P, N}, and
4387 write the bitwise complement of the result to {RP, N}.
4389 -- Function: void mpn_nior_n (mp_limb_t *RP, const mp_limb_t *S1P,
4390 const mp_limb_t *S2P, mp_size_t N)
4391 Perform the bitwise logical inclusive or of {S1P, N} and {S2P, N},
4392 and write the bitwise complement of the result to {RP, N}.
4394 -- Function: void mpn_xnor_n (mp_limb_t *RP, const mp_limb_t *S1P,
4395 const mp_limb_t *S2P, mp_size_t N)
4396 Perform the bitwise logical exclusive or of {S1P, N} and {S2P, N},
4397 and write the bitwise complement of the result to {RP, N}.
4399 -- Function: void mpn_com (mp_limb_t *RP, const mp_limb_t *SP,
4401 Perform the bitwise complement of {SP, N}, and write the result to
4404 -- Function: void mpn_copyi (mp_limb_t *RP, const mp_limb_t *S1P,
4406 Copy from {S1P, N} to {RP, N}, increasingly.
4408 -- Function: void mpn_copyd (mp_limb_t *RP, const mp_limb_t *S1P,
4410 Copy from {S1P, N} to {RP, N}, decreasingly.
4412 -- Function: void mpn_zero (mp_limb_t *RP, mp_size_t N)
4419 *Everything in this section is highly experimental and may disappear or
4420 be subject to incompatible changes in a future version of GMP.*
4422 Nails are an experimental feature whereby a few bits are left unused
4423 at the top of each `mp_limb_t'. This can significantly improve carry
4424 handling on some processors.
4426 All the `mpn' functions accepting limb data will expect the nail
4427 bits to be zero on entry, and will return data with the nails similarly
4428 all zero. This applies both to limb vectors and to single limb
4431 Nails can be enabled by configuring with `--enable-nails'. By
4432 default the number of bits will be chosen according to what suits the
4433 host processor, but a particular number can be selected with
4436 At the mpn level, a nail build is neither source nor binary
4437 compatible with a non-nail build, strictly speaking. But programs
4438 acting on limbs only through the mpn functions are likely to work
4439 equally well with either build, and judicious use of the definitions
4440 below should make any program compatible with either build, at the
4443 For the higher level routines, meaning `mpz' etc, a nail build
4444 should be fully source and binary compatible with a non-nail build.
4446 -- Macro: GMP_NAIL_BITS
4447 -- Macro: GMP_NUMB_BITS
4448 -- Macro: GMP_LIMB_BITS
4449 `GMP_NAIL_BITS' is the number of nail bits, or 0 when nails are
4450 not in use. `GMP_NUMB_BITS' is the number of data bits in a limb.
4451 `GMP_LIMB_BITS' is the total number of bits in an `mp_limb_t'. In
4454 GMP_LIMB_BITS == GMP_NAIL_BITS + GMP_NUMB_BITS
4456 -- Macro: GMP_NAIL_MASK
4457 -- Macro: GMP_NUMB_MASK
4458 Bit masks for the nail and number parts of a limb.
4459 `GMP_NAIL_MASK' is 0 when nails are not in use.
4461 `GMP_NAIL_MASK' is not often needed, since the nail part can be
4462 obtained with `x >> GMP_NUMB_BITS', and that means one less large
4463 constant, which can help various RISC chips.
4465 -- Macro: GMP_NUMB_MAX
4466 The maximum value that can be stored in the number part of a limb.
4467 This is the same as `GMP_NUMB_MASK', but can be used for clarity
4468 when doing comparisons rather than bit-wise operations.
4470 The term "nails" comes from finger or toe nails, which are at the
4471 ends of a limb (arm or leg). "numb" is short for number, but is also
4472 how the developers felt after trying for a long time to come up with
4473 sensible names for these things.
4475 In the future (the distant future most likely) a non-zero nail might
4476 be permitted, giving non-unique representations for numbers in a limb
4477 vector. This would help vector processors since carries would only
4478 ever need to propagate one or two limbs.
4481 File: gmp.info, Node: Random Number Functions, Next: Formatted Output, Prev: Low-level Functions, Up: Top
4483 9 Random Number Functions
4484 *************************
4486 Sequences of pseudo-random numbers in GMP are generated using a
4487 variable of type `gmp_randstate_t', which holds an algorithm selection
4488 and a current state. Such a variable must be initialized by a call to
4489 one of the `gmp_randinit' functions, and can be seeded with one of the
4490 `gmp_randseed' functions.
4492 The functions actually generating random numbers are described in
4493 *Note Integer Random Numbers::, and *Note Miscellaneous Float
4496 The older style random number functions don't accept a
4497 `gmp_randstate_t' parameter but instead share a global variable of that
4498 type. They use a default algorithm and are currently not seeded
4499 (though perhaps that will change in the future). The new functions
4500 accepting a `gmp_randstate_t' are recommended for applications that
4501 care about randomness.
4505 * Random State Initialization::
4506 * Random State Seeding::
4507 * Random State Miscellaneous::
4510 File: gmp.info, Node: Random State Initialization, Next: Random State Seeding, Prev: Random Number Functions, Up: Random Number Functions
4512 9.1 Random State Initialization
4513 ===============================
4515 -- Function: void gmp_randinit_default (gmp_randstate_t STATE)
4516 Initialize STATE with a default algorithm. This will be a
4517 compromise between speed and randomness, and is recommended for
4518 applications with no special requirements. Currently this is
4521 -- Function: void gmp_randinit_mt (gmp_randstate_t STATE)
4522 Initialize STATE for a Mersenne Twister algorithm. This algorithm
4523 is fast and has good randomness properties.
4525 -- Function: void gmp_randinit_lc_2exp (gmp_randstate_t STATE, mpz_t
4526 A, unsigned long C, mp_bitcnt_t M2EXP)
4527 Initialize STATE with a linear congruential algorithm X = (A*X +
4530 The low bits of X in this algorithm are not very random. The least
4531 significant bit will have a period no more than 2, and the second
4532 bit no more than 4, etc. For this reason only the high half of
4533 each X is actually used.
4535 When a random number of more than M2EXP/2 bits is to be generated,
4536 multiple iterations of the recurrence are used and the results
4539 -- Function: int gmp_randinit_lc_2exp_size (gmp_randstate_t STATE,
4541 Initialize STATE for a linear congruential algorithm as per
4542 `gmp_randinit_lc_2exp'. A, C and M2EXP are selected from a table,
4543 chosen so that SIZE bits (or more) of each X will be used, ie.
4546 If successful the return value is non-zero. If SIZE is bigger
4547 than the table data provides then the return value is zero. The
4548 maximum SIZE currently supported is 128.
4550 -- Function: void gmp_randinit_set (gmp_randstate_t ROP,
4552 Initialize ROP with a copy of the algorithm and state from OP.
4554 -- Function: void gmp_randinit (gmp_randstate_t STATE,
4555 gmp_randalg_t ALG, ...)
4556 *This function is obsolete.*
4558 Initialize STATE with an algorithm selected by ALG. The only
4559 choice is `GMP_RAND_ALG_LC', which is `gmp_randinit_lc_2exp_size'
4560 described above. A third parameter of type `unsigned long' is
4561 required, this is the SIZE for that function.
4562 `GMP_RAND_ALG_DEFAULT' or 0 are the same as `GMP_RAND_ALG_LC'.
4564 `gmp_randinit' sets bits in the global variable `gmp_errno' to
4565 indicate an error. `GMP_ERROR_UNSUPPORTED_ARGUMENT' if ALG is
4566 unsupported, or `GMP_ERROR_INVALID_ARGUMENT' if the SIZE parameter
4567 is too big. It may be noted this error reporting is not thread
4568 safe (a good reason to use `gmp_randinit_lc_2exp_size' instead).
4570 -- Function: void gmp_randclear (gmp_randstate_t STATE)
4571 Free all memory occupied by STATE.
4574 File: gmp.info, Node: Random State Seeding, Next: Random State Miscellaneous, Prev: Random State Initialization, Up: Random Number Functions
4576 9.2 Random State Seeding
4577 ========================
4579 -- Function: void gmp_randseed (gmp_randstate_t STATE, mpz_t SEED)
4580 -- Function: void gmp_randseed_ui (gmp_randstate_t STATE,
4581 unsigned long int SEED)
4582 Set an initial seed value into STATE.
4584 The size of a seed determines how many different sequences of
4585 random numbers that it's possible to generate. The "quality" of
4586 the seed is the randomness of a given seed compared to the
4587 previous seed used, and this affects the randomness of separate
4588 number sequences. The method for choosing a seed is critical if
4589 the generated numbers are to be used for important applications,
4590 such as generating cryptographic keys.
4592 Traditionally the system time has been used to seed, but care
4593 needs to be taken with this. If an application seeds often and
4594 the resolution of the system clock is low, then the same sequence
4595 of numbers might be repeated. Also, the system time is quite easy
4596 to guess, so if unpredictability is required then it should
4597 definitely not be the only source for the seed value. On some
4598 systems there's a special device `/dev/random' which provides
4599 random data better suited for use as a seed.
4602 File: gmp.info, Node: Random State Miscellaneous, Prev: Random State Seeding, Up: Random Number Functions
4604 9.3 Random State Miscellaneous
4605 ==============================
4607 -- Function: unsigned long gmp_urandomb_ui (gmp_randstate_t STATE,
4609 Return a uniformly distributed random number of N bits, ie. in the
4610 range 0 to 2^N-1 inclusive. N must be less than or equal to the
4611 number of bits in an `unsigned long'.
4613 -- Function: unsigned long gmp_urandomm_ui (gmp_randstate_t STATE,
4615 Return a uniformly distributed random number in the range 0 to
4619 File: gmp.info, Node: Formatted Output, Next: Formatted Input, Prev: Random Number Functions, Up: Top
4626 * Formatted Output Strings::
4627 * Formatted Output Functions::
4628 * C++ Formatted Output::
4631 File: gmp.info, Node: Formatted Output Strings, Next: Formatted Output Functions, Prev: Formatted Output, Up: Formatted Output
4636 `gmp_printf' and friends accept format strings similar to the standard C
4637 `printf' (*note Formatted Output: (libc)Formatted Output.). A format
4638 specification is of the form
4640 % [flags] [width] [.[precision]] [type] conv
4642 GMP adds types `Z', `Q' and `F' for `mpz_t', `mpq_t' and `mpf_t'
4643 respectively, `M' for `mp_limb_t', and `N' for an `mp_limb_t' array.
4644 `Z', `Q', `M' and `N' behave like integers. `Q' will print a `/' and a
4645 denominator, if needed. `F' behaves like a float. For example,
4648 gmp_printf ("%s is an mpz %Zd\n", "here", z);
4651 gmp_printf ("a hex rational: %#40Qx\n", q);
4655 gmp_printf ("fixed point mpf %.*Ff with %d digits\n", n, f, n);
4658 gmp_printf ("limb %Mu\n", l);
4660 const mp_limb_t *ptr;
4662 gmp_printf ("limb array %Nx\n", ptr, size);
4664 For `N' the limbs are expected least significant first, as per the
4665 `mpn' functions (*note Low-level Functions::). A negative size can be
4666 given to print the value as a negative.
4668 All the standard C `printf' types behave the same as the C library
4669 `printf', and can be freely intermixed with the GMP extensions. In the
4670 current implementation the standard parts of the format string are
4671 simply handed to `printf' and only the GMP extensions handled directly.
4673 The flags accepted are as follows. GLIBC style ' is only for the
4674 standard C types (not the GMP types), and only if the C library
4677 0 pad with zeros (rather than spaces)
4678 # show the base with `0x', `0X' or `0'
4679 + always show a sign
4680 (space) show a space or a `-' sign
4681 ' group digits, GLIBC style (not GMP types)
4683 The optional width and precision can be given as a number within the
4684 format string, or as a `*' to take an extra parameter of type `int', the
4685 same as the standard `printf'.
4687 The standard types accepted are as follows. `h' and `l' are
4688 portable, the rest will depend on the compiler (or include files) for
4689 the type and the C library for the output.
4693 j intmax_t or uintmax_t
4697 q quad_t or u_quad_t
4703 F mpf_t, float conversions
4704 Q mpq_t, integer conversions
4705 M mp_limb_t, integer conversions
4706 N mp_limb_t array, integer conversions
4707 Z mpz_t, integer conversions
4709 The conversions accepted are as follows. `a' and `A' are always
4710 supported for `mpf_t' but depend on the C library for standard C float
4711 types. `m' and `p' depend on the C library.
4713 a A hex floats, C99 style
4716 e E scientific format float
4719 g G fixed or scientific float
4720 m `strerror' string, GLIBC style
4721 n store characters written so far
4728 `o', `x' and `X' are unsigned for the standard C types, but for
4729 types `Z', `Q' and `N' they are signed. `u' is not meaningful for `Z',
4732 `M' is a proxy for the C library `l' or `L', according to the size
4733 of `mp_limb_t'. Unsigned conversions will be usual, but a signed
4734 conversion can be used and will interpret the value as a twos complement
4737 `n' can be used with any type, even the GMP types.
4739 Other types or conversions that might be accepted by the C library
4740 `printf' cannot be used through `gmp_printf', this includes for
4741 instance extensions registered with GLIBC `register_printf_function'.
4742 Also currently there's no support for POSIX `$' style numbered arguments
4743 (perhaps this will be added in the future).
4745 The precision field has it's usual meaning for integer `Z' and float
4746 `F' types, but is currently undefined for `Q' and should not be used
4749 `mpf_t' conversions only ever generate as many digits as can be
4750 accurately represented by the operand, the same as `mpf_get_str' does.
4751 Zeros will be used if necessary to pad to the requested precision. This
4752 happens even for an `f' conversion of an `mpf_t' which is an integer,
4753 for instance 2^1024 in an `mpf_t' of 128 bits precision will only
4754 produce about 40 digits, then pad with zeros to the decimal point. An
4755 empty precision field like `%.Fe' or `%.Ff' can be used to specifically
4756 request just the significant digits.
4758 The decimal point character (or string) is taken from the current
4759 locale settings on systems which provide `localeconv' (*note Locales
4760 and Internationalization: (libc)Locales.). The C library will normally
4761 do the same for standard float output.
4763 The format string is only interpreted as plain `char's, multibyte
4764 characters are not recognised. Perhaps this will change in the future.
4767 File: gmp.info, Node: Formatted Output Functions, Next: C++ Formatted Output, Prev: Formatted Output Strings, Up: Formatted Output
4772 Each of the following functions is similar to the corresponding C
4773 library function. The basic `printf' forms take a variable argument
4774 list. The `vprintf' forms take an argument pointer, see *Note Variadic
4775 Functions: (libc)Variadic Functions, or `man 3 va_start'.
4777 It should be emphasised that if a format string is invalid, or the
4778 arguments don't match what the format specifies, then the behaviour of
4779 any of these functions will be unpredictable. GCC format string
4780 checking is not available, since it doesn't recognise the GMP
4783 The file based functions `gmp_printf' and `gmp_fprintf' will return
4784 -1 to indicate a write error. Output is not "atomic", so partial
4785 output may be produced if a write error occurs. All the functions can
4786 return -1 if the C library `printf' variant in use returns -1, but this
4787 shouldn't normally occur.
4789 -- Function: int gmp_printf (const char *FMT, ...)
4790 -- Function: int gmp_vprintf (const char *FMT, va_list AP)
4791 Print to the standard output `stdout'. Return the number of
4792 characters written, or -1 if an error occurred.
4794 -- Function: int gmp_fprintf (FILE *FP, const char *FMT, ...)
4795 -- Function: int gmp_vfprintf (FILE *FP, const char *FMT, va_list AP)
4796 Print to the stream FP. Return the number of characters written,
4797 or -1 if an error occurred.
4799 -- Function: int gmp_sprintf (char *BUF, const char *FMT, ...)
4800 -- Function: int gmp_vsprintf (char *BUF, const char *FMT, va_list AP)
4801 Form a null-terminated string in BUF. Return the number of
4802 characters written, excluding the terminating null.
4804 No overlap is permitted between the space at BUF and the string
4807 These functions are not recommended, since there's no protection
4808 against exceeding the space available at BUF.
4810 -- Function: int gmp_snprintf (char *BUF, size_t SIZE, const char
4812 -- Function: int gmp_vsnprintf (char *BUF, size_t SIZE, const char
4814 Form a null-terminated string in BUF. No more than SIZE bytes
4815 will be written. To get the full output, SIZE must be enough for
4816 the string and null-terminator.
4818 The return value is the total number of characters which ought to
4819 have been produced, excluding the terminating null. If RETVAL >=
4820 SIZE then the actual output has been truncated to the first SIZE-1
4821 characters, and a null appended.
4823 No overlap is permitted between the region {BUF,SIZE} and the FMT
4826 Notice the return value is in ISO C99 `snprintf' style. This is
4827 so even if the C library `vsnprintf' is the older GLIBC 2.0.x
4830 -- Function: int gmp_asprintf (char **PP, const char *FMT, ...)
4831 -- Function: int gmp_vasprintf (char **PP, const char *FMT, va_list AP)
4832 Form a null-terminated string in a block of memory obtained from
4833 the current memory allocation function (*note Custom
4834 Allocation::). The block will be the size of the string and
4835 null-terminator. The address of the block in stored to *PP. The
4836 return value is the number of characters produced, excluding the
4839 Unlike the C library `asprintf', `gmp_asprintf' doesn't return -1
4840 if there's no more memory available, it lets the current allocation
4841 function handle that.
4843 -- Function: int gmp_obstack_printf (struct obstack *OB, const char
4845 -- Function: int gmp_obstack_vprintf (struct obstack *OB, const char
4847 Append to the current object in OB. The return value is the
4848 number of characters written. A null-terminator is not written.
4850 FMT cannot be within the current object in OB, since that object
4851 might move as it grows.
4853 These functions are available only when the C library provides the
4854 obstack feature, which probably means only on GNU systems, see
4855 *Note Obstacks: (libc)Obstacks.
4858 File: gmp.info, Node: C++ Formatted Output, Prev: Formatted Output Functions, Up: Formatted Output
4860 10.3 C++ Formatted Output
4861 =========================
4863 The following functions are provided in `libgmpxx' (*note Headers and
4864 Libraries::), which is built if C++ support is enabled (*note Build
4865 Options::). Prototypes are available from `<gmp.h>'.
4867 -- Function: ostream& operator<< (ostream& STREAM, mpz_t OP)
4868 Print OP to STREAM, using its `ios' formatting settings.
4869 `ios::width' is reset to 0 after output, the same as the standard
4870 `ostream operator<<' routines do.
4872 In hex or octal, OP is printed as a signed number, the same as for
4873 decimal. This is unlike the standard `operator<<' routines on
4874 `int' etc, which instead give twos complement.
4876 -- Function: ostream& operator<< (ostream& STREAM, mpq_t OP)
4877 Print OP to STREAM, using its `ios' formatting settings.
4878 `ios::width' is reset to 0 after output, the same as the standard
4879 `ostream operator<<' routines do.
4881 Output will be a fraction like `5/9', or if the denominator is 1
4882 then just a plain integer like `123'.
4884 In hex or octal, OP is printed as a signed value, the same as for
4885 decimal. If `ios::showbase' is set then a base indicator is shown
4886 on both the numerator and denominator (if the denominator is
4889 -- Function: ostream& operator<< (ostream& STREAM, mpf_t OP)
4890 Print OP to STREAM, using its `ios' formatting settings.
4891 `ios::width' is reset to 0 after output, the same as the standard
4892 `ostream operator<<' routines do.
4894 The decimal point follows the standard library float `operator<<',
4895 which on recent systems means the `std::locale' imbued on STREAM.
4897 Hex and octal are supported, unlike the standard `operator<<' on
4898 `double'. The mantissa will be in hex or octal, the exponent will
4899 be in decimal. For hex the exponent delimiter is an `@'. This is
4900 as per `mpf_out_str'.
4902 `ios::showbase' is supported, and will put a base on the mantissa,
4903 for example hex `0x1.8' or `0x0.8', or octal `01.4' or `00.4'.
4904 This last form is slightly strange, but at least differentiates
4905 itself from decimal.
4907 These operators mean that GMP types can be printed in the usual C++
4913 cout << "iteration " << n << " value " << z << "\n";
4915 But note that `ostream' output (and `istream' input, *note C++
4916 Formatted Input::) is the only overloading available for the GMP types
4917 and that for instance using `+' with an `mpz_t' will have unpredictable
4918 results. For classes with overloading, see *Note C++ Class Interface::.
4921 File: gmp.info, Node: Formatted Input, Next: C++ Class Interface, Prev: Formatted Output, Up: Top
4928 * Formatted Input Strings::
4929 * Formatted Input Functions::
4930 * C++ Formatted Input::
4933 File: gmp.info, Node: Formatted Input Strings, Next: Formatted Input Functions, Prev: Formatted Input, Up: Formatted Input
4935 11.1 Formatted Input Strings
4936 ============================
4938 `gmp_scanf' and friends accept format strings similar to the standard C
4939 `scanf' (*note Formatted Input: (libc)Formatted Input.). A format
4940 specification is of the form
4942 % [flags] [width] [type] conv
4944 GMP adds types `Z', `Q' and `F' for `mpz_t', `mpq_t' and `mpf_t'
4945 respectively. `Z' and `Q' behave like integers. `Q' will read a `/'
4946 and a denominator, if present. `F' behaves like a float.
4948 GMP variables don't require an `&' when passed to `gmp_scanf', since
4949 they're already "call-by-reference". For example,
4951 /* to read say "a(5) = 1234" */
4954 gmp_scanf ("a(%d) = %Zd\n", &n, z);
4957 gmp_sscanf ("0377 + 0x10/0x11", "%Qi + %Qi", q1, q2);
4959 /* to read say "topleft (1.55,-2.66)" */
4962 gmp_scanf ("%31s (%Ff,%Ff)", buf, x, y);
4964 All the standard C `scanf' types behave the same as in the C library
4965 `scanf', and can be freely intermixed with the GMP extensions. In the
4966 current implementation the standard parts of the format string are
4967 simply handed to `scanf' and only the GMP extensions handled directly.
4969 The flags accepted are as follows. `a' and `'' will depend on
4970 support from the C library, and `'' cannot be used with GMP types.
4972 * read but don't store
4973 a allocate a buffer (string conversions)
4974 ' grouped digits, GLIBC style (not GMP
4977 The standard types accepted are as follows. `h' and `l' are
4978 portable, the rest will depend on the compiler (or include files) for
4979 the type and the C library for the input.
4983 j intmax_t or uintmax_t
4984 l long int, double or wchar_t
4987 q quad_t or u_quad_t
4993 F mpf_t, float conversions
4994 Q mpq_t, integer conversions
4995 Z mpz_t, integer conversions
4997 The conversions accepted are as follows. `p' and `[' will depend on
4998 support from the C library, the rest are standard.
5000 c character or characters
5003 i integer with base indicator
5004 n characters read so far
5007 s string of non-whitespace characters
5010 [ string of characters in a set
5012 `e', `E', `f', `g' and `G' are identical, they all read either fixed
5013 point or scientific format, and either upper or lower case `e' for the
5014 exponent in scientific format.
5016 C99 style hex float format (`printf %a', *note Formatted Output
5017 Strings::) is always accepted for `mpf_t', but for the standard float
5018 types it will depend on the C library.
5020 `x' and `X' are identical, both accept both upper and lower case
5023 `o', `u', `x' and `X' all read positive or negative values. For the
5024 standard C types these are described as "unsigned" conversions, but
5025 that merely affects certain overflow handling, negatives are still
5026 allowed (per `strtoul', *note Parsing of Integers: (libc)Parsing of
5027 Integers.). For GMP types there are no overflows, so `d' and `u' are
5030 `Q' type reads the numerator and (optional) denominator as given.
5031 If the value might not be in canonical form then `mpq_canonicalize'
5032 must be called before using it in any calculations (*note Rational
5033 Number Functions::).
5035 `Qi' will read a base specification separately for the numerator and
5036 denominator. For example `0x10/11' would be 16/11, whereas `0x10/0x11'
5039 `n' can be used with any of the types above, even the GMP types.
5040 `*' to suppress assignment is allowed, though in that case it would do
5043 Other conversions or types that might be accepted by the C library
5044 `scanf' cannot be used through `gmp_scanf'.
5046 Whitespace is read and discarded before a field, except for `c' and
5049 For float conversions, the decimal point character (or string)
5050 expected is taken from the current locale settings on systems which
5051 provide `localeconv' (*note Locales and Internationalization:
5052 (libc)Locales.). The C library will normally do the same for standard
5055 The format string is only interpreted as plain `char's, multibyte
5056 characters are not recognised. Perhaps this will change in the future.
5059 File: gmp.info, Node: Formatted Input Functions, Next: C++ Formatted Input, Prev: Formatted Input Strings, Up: Formatted Input
5061 11.2 Formatted Input Functions
5062 ==============================
5064 Each of the following functions is similar to the corresponding C
5065 library function. The plain `scanf' forms take a variable argument
5066 list. The `vscanf' forms take an argument pointer, see *Note Variadic
5067 Functions: (libc)Variadic Functions, or `man 3 va_start'.
5069 It should be emphasised that if a format string is invalid, or the
5070 arguments don't match what the format specifies, then the behaviour of
5071 any of these functions will be unpredictable. GCC format string
5072 checking is not available, since it doesn't recognise the GMP
5075 No overlap is permitted between the FMT string and any of the results
5078 -- Function: int gmp_scanf (const char *FMT, ...)
5079 -- Function: int gmp_vscanf (const char *FMT, va_list AP)
5080 Read from the standard input `stdin'.
5082 -- Function: int gmp_fscanf (FILE *FP, const char *FMT, ...)
5083 -- Function: int gmp_vfscanf (FILE *FP, const char *FMT, va_list AP)
5084 Read from the stream FP.
5086 -- Function: int gmp_sscanf (const char *S, const char *FMT, ...)
5087 -- Function: int gmp_vsscanf (const char *S, const char *FMT, va_list
5089 Read from a null-terminated string S.
5091 The return value from each of these functions is the same as the
5092 standard C99 `scanf', namely the number of fields successfully parsed
5093 and stored. `%n' fields and fields read but suppressed by `*' don't
5094 count towards the return value.
5096 If end of input (or a file error) is reached before a character for
5097 a field or a literal, and if no previous non-suppressed fields have
5098 matched, then the return value is `EOF' instead of 0. A whitespace
5099 character in the format string is only an optional match and doesn't
5100 induce an `EOF' in this fashion. Leading whitespace read and discarded
5101 for a field don't count as characters for that field.
5103 For the GMP types, input parsing follows C99 rules, namely one
5104 character of lookahead is used and characters are read while they
5105 continue to meet the format requirements. If this doesn't provide a
5106 complete number then the function terminates, with that field not
5107 stored nor counted towards the return value. For instance with `mpf_t'
5108 an input `1.23e-XYZ' would be read up to the `X' and that character
5109 pushed back since it's not a digit. The string `1.23e-' would then be
5110 considered invalid since an `e' must be followed by at least one digit.
5112 For the standard C types, in the current implementation GMP calls
5113 the C library `scanf' functions, which might have looser rules about
5114 what constitutes a valid input.
5116 Note that `gmp_sscanf' is the same as `gmp_fscanf' and only does one
5117 character of lookahead when parsing. Although clearly it could look at
5118 its entire input, it is deliberately made identical to `gmp_fscanf',
5119 the same way C99 `sscanf' is the same as `fscanf'.
5122 File: gmp.info, Node: C++ Formatted Input, Prev: Formatted Input Functions, Up: Formatted Input
5124 11.3 C++ Formatted Input
5125 ========================
5127 The following functions are provided in `libgmpxx' (*note Headers and
5128 Libraries::), which is built only if C++ support is enabled (*note
5129 Build Options::). Prototypes are available from `<gmp.h>'.
5131 -- Function: istream& operator>> (istream& STREAM, mpz_t ROP)
5132 Read ROP from STREAM, using its `ios' formatting settings.
5134 -- Function: istream& operator>> (istream& STREAM, mpq_t ROP)
5135 An integer like `123' will be read, or a fraction like `5/9'. No
5136 whitespace is allowed around the `/'. If the fraction is not in
5137 canonical form then `mpq_canonicalize' must be called (*note
5138 Rational Number Functions::) before operating on it.
5140 As per integer input, an `0' or `0x' base indicator is read when
5141 none of `ios::dec', `ios::oct' or `ios::hex' are set. This is
5142 done separately for numerator and denominator, so that for instance
5143 `0x10/11' is 16/11 and `0x10/0x11' is 16/17.
5145 -- Function: istream& operator>> (istream& STREAM, mpf_t ROP)
5146 Read ROP from STREAM, using its `ios' formatting settings.
5148 Hex or octal floats are not supported, but might be in the future,
5149 or perhaps it's best to accept only what the standard float
5152 Note that digit grouping specified by the `istream' locale is
5153 currently not accepted. Perhaps this will change in the future.
5156 These operators mean that GMP types can be read in the usual C++
5163 But note that `istream' input (and `ostream' output, *note C++
5164 Formatted Output::) is the only overloading available for the GMP types
5165 and that for instance using `+' with an `mpz_t' will have unpredictable
5166 results. For classes with overloading, see *Note C++ Class Interface::.
5169 File: gmp.info, Node: C++ Class Interface, Next: BSD Compatible Functions, Prev: Formatted Input, Up: Top
5171 12 C++ Class Interface
5172 **********************
5174 This chapter describes the C++ class based interface to GMP.
5176 All GMP C language types and functions can be used in C++ programs,
5177 since `gmp.h' has `extern "C"' qualifiers, but the class interface
5178 offers overloaded functions and operators which may be more convenient.
5180 Due to the implementation of this interface, a reasonably recent C++
5181 compiler is required, one supporting namespaces, partial specialization
5182 of templates and member templates. For GCC this means version 2.91 or
5185 *Everything described in this chapter is to be considered preliminary
5186 and might be subject to incompatible changes if some unforeseen
5187 difficulty reveals itself.*
5191 * C++ Interface General::
5192 * C++ Interface Integers::
5193 * C++ Interface Rationals::
5194 * C++ Interface Floats::
5195 * C++ Interface Random Numbers::
5196 * C++ Interface Limitations::
5199 File: gmp.info, Node: C++ Interface General, Next: C++ Interface Integers, Prev: C++ Class Interface, Up: C++ Class Interface
5201 12.1 C++ Interface General
5202 ==========================
5204 All the C++ classes and functions are available with
5208 Programs should be linked with the `libgmpxx' and `libgmp'
5209 libraries. For example,
5211 g++ mycxxprog.cc -lgmpxx -lgmp
5213 The classes defined are
5219 The standard operators and various standard functions are overloaded
5220 to allow arithmetic with these classes. For example,
5230 cout << "sum is " << c << "\n";
5231 cout << "absolute value is " << abs(c) << "\n";
5236 An important feature of the implementation is that an expression like
5237 `a=b+c' results in a single call to the corresponding `mpz_add',
5238 without using a temporary for the `b+c' part. Expressions which by
5239 their nature imply intermediate values, like `a=b*c+d*e', still use
5242 The classes can be freely intermixed in expressions, as can the
5243 classes and the standard types `long', `unsigned long' and `double'.
5244 Smaller types like `int' or `float' can also be intermixed, since C++
5247 Note that `bool' is not accepted directly, but must be explicitly
5248 cast to an `int' first. This is because C++ will automatically convert
5249 any pointer to a `bool', so if GMP accepted `bool' it would make all
5250 sorts of invalid class and pointer combinations compile but almost
5251 certainly not do anything sensible.
5253 Conversions back from the classes to standard C++ types aren't done
5254 automatically, instead member functions like `get_si' are provided (see
5255 the following sections for details).
5257 Also there are no automatic conversions from the classes to the
5258 corresponding GMP C types, instead a reference to the underlying C
5259 object can be obtained with the following functions,
5261 -- Function: mpz_t mpz_class::get_mpz_t ()
5262 -- Function: mpq_t mpq_class::get_mpq_t ()
5263 -- Function: mpf_t mpf_class::get_mpf_t ()
5265 These can be used to call a C function which doesn't have a C++ class
5266 interface. For example to set `a' to the GCD of `b' and `c',
5270 mpz_gcd (a.get_mpz_t(), b.get_mpz_t(), c.get_mpz_t());
5272 In the other direction, a class can be initialized from the
5273 corresponding GMP C type, or assigned to if an explicit constructor is
5274 used. In both cases this makes a copy of the value, it doesn't create
5275 any sort of association. For example,
5278 // ... init and calculate z ...
5283 There are no namespace setups in `gmpxx.h', all types and functions
5284 are simply put into the global namespace. This is what `gmp.h' has
5285 done in the past, and continues to do for compatibility. The extras
5286 provided by `gmpxx.h' follow GMP naming conventions and are unlikely to
5287 clash with anything.
5290 File: gmp.info, Node: C++ Interface Integers, Next: C++ Interface Rationals, Prev: C++ Interface General, Up: C++ Class Interface
5292 12.2 C++ Interface Integers
5293 ===========================
5295 -- Function: void mpz_class::mpz_class (type N)
5296 Construct an `mpz_class'. All the standard C++ types may be used,
5297 except `long long' and `long double', and all the GMP C++ classes
5298 can be used. Any necessary conversion follows the corresponding C
5299 function, for example `double' follows `mpz_set_d' (*note
5300 Assigning Integers::).
5302 -- Function: void mpz_class::mpz_class (mpz_t Z)
5303 Construct an `mpz_class' from an `mpz_t'. The value in Z is
5304 copied into the new `mpz_class', there won't be any permanent
5305 association between it and Z.
5307 -- Function: void mpz_class::mpz_class (const char *S)
5308 -- Function: void mpz_class::mpz_class (const char *S, int BASE = 0)
5309 -- Function: void mpz_class::mpz_class (const string& S)
5310 -- Function: void mpz_class::mpz_class (const string& S, int BASE = 0)
5311 Construct an `mpz_class' converted from a string using
5312 `mpz_set_str' (*note Assigning Integers::).
5314 If the string is not a valid integer, an `std::invalid_argument'
5315 exception is thrown. The same applies to `operator='.
5317 -- Function: mpz_class operator/ (mpz_class A, mpz_class D)
5318 -- Function: mpz_class operator% (mpz_class A, mpz_class D)
5319 Divisions involving `mpz_class' round towards zero, as per the
5320 `mpz_tdiv_q' and `mpz_tdiv_r' functions (*note Integer Division::).
5321 This is the same as the C99 `/' and `%' operators.
5323 The `mpz_fdiv...' or `mpz_cdiv...' functions can always be called
5324 directly if desired. For example,
5328 mpz_fdiv_q (q.get_mpz_t(), a.get_mpz_t(), d.get_mpz_t());
5330 -- Function: mpz_class abs (mpz_class OP1)
5331 -- Function: int cmp (mpz_class OP1, type OP2)
5332 -- Function: int cmp (type OP1, mpz_class OP2)
5333 -- Function: bool mpz_class::fits_sint_p (void)
5334 -- Function: bool mpz_class::fits_slong_p (void)
5335 -- Function: bool mpz_class::fits_sshort_p (void)
5336 -- Function: bool mpz_class::fits_uint_p (void)
5337 -- Function: bool mpz_class::fits_ulong_p (void)
5338 -- Function: bool mpz_class::fits_ushort_p (void)
5339 -- Function: double mpz_class::get_d (void)
5340 -- Function: long mpz_class::get_si (void)
5341 -- Function: string mpz_class::get_str (int BASE = 10)
5342 -- Function: unsigned long mpz_class::get_ui (void)
5343 -- Function: int mpz_class::set_str (const char *STR, int BASE)
5344 -- Function: int mpz_class::set_str (const string& STR, int BASE)
5345 -- Function: int sgn (mpz_class OP)
5346 -- Function: mpz_class sqrt (mpz_class OP)
5347 These functions provide a C++ class interface to the corresponding
5350 `cmp' can be used with any of the classes or the standard C++
5351 types, except `long long' and `long double'.
5354 Overloaded operators for combinations of `mpz_class' and `double'
5355 are provided for completeness, but it should be noted that if the given
5356 `double' is not an integer then the way any rounding is done is
5357 currently unspecified. The rounding might take place at the start, in
5358 the middle, or at the end of the operation, and it might change in the
5361 Conversions between `mpz_class' and `double', however, are defined
5362 to follow the corresponding C functions `mpz_get_d' and `mpz_set_d'.
5363 And comparisons are always made exactly, as per `mpz_cmp_d'.
5366 File: gmp.info, Node: C++ Interface Rationals, Next: C++ Interface Floats, Prev: C++ Interface Integers, Up: C++ Class Interface
5368 12.3 C++ Interface Rationals
5369 ============================
5371 In all the following constructors, if a fraction is given then it
5372 should be in canonical form, or if not then `mpq_class::canonicalize'
5375 -- Function: void mpq_class::mpq_class (type OP)
5376 -- Function: void mpq_class::mpq_class (integer NUM, integer DEN)
5377 Construct an `mpq_class'. The initial value can be a single value
5378 of any type, or a pair of integers (`mpz_class' or standard C++
5379 integer types) representing a fraction, except that `long long'
5380 and `long double' are not supported. For example,
5386 -- Function: void mpq_class::mpq_class (mpq_t Q)
5387 Construct an `mpq_class' from an `mpq_t'. The value in Q is
5388 copied into the new `mpq_class', there won't be any permanent
5389 association between it and Q.
5391 -- Function: void mpq_class::mpq_class (const char *S)
5392 -- Function: void mpq_class::mpq_class (const char *S, int BASE = 0)
5393 -- Function: void mpq_class::mpq_class (const string& S)
5394 -- Function: void mpq_class::mpq_class (const string& S, int BASE = 0)
5395 Construct an `mpq_class' converted from a string using
5396 `mpq_set_str' (*note Initializing Rationals::).
5398 If the string is not a valid rational, an `std::invalid_argument'
5399 exception is thrown. The same applies to `operator='.
5401 -- Function: void mpq_class::canonicalize ()
5402 Put an `mpq_class' into canonical form, as per *Note Rational
5403 Number Functions::. All arithmetic operators require their
5404 operands in canonical form, and will return results in canonical
5407 -- Function: mpq_class abs (mpq_class OP)
5408 -- Function: int cmp (mpq_class OP1, type OP2)
5409 -- Function: int cmp (type OP1, mpq_class OP2)
5410 -- Function: double mpq_class::get_d (void)
5411 -- Function: string mpq_class::get_str (int BASE = 10)
5412 -- Function: int mpq_class::set_str (const char *STR, int BASE)
5413 -- Function: int mpq_class::set_str (const string& STR, int BASE)
5414 -- Function: int sgn (mpq_class OP)
5415 These functions provide a C++ class interface to the corresponding
5418 `cmp' can be used with any of the classes or the standard C++
5419 types, except `long long' and `long double'.
5421 -- Function: mpz_class& mpq_class::get_num ()
5422 -- Function: mpz_class& mpq_class::get_den ()
5423 Get a reference to an `mpz_class' which is the numerator or
5424 denominator of an `mpq_class'. This can be used both for read and
5425 write access. If the object returned is modified, it modifies the
5426 original `mpq_class'.
5428 If direct manipulation might produce a non-canonical value, then
5429 `mpq_class::canonicalize' must be called before further operations.
5431 -- Function: mpz_t mpq_class::get_num_mpz_t ()
5432 -- Function: mpz_t mpq_class::get_den_mpz_t ()
5433 Get a reference to the underlying `mpz_t' numerator or denominator
5434 of an `mpq_class'. This can be passed to C functions expecting an
5435 `mpz_t'. Any modifications made to the `mpz_t' will modify the
5436 original `mpq_class'.
5438 If direct manipulation might produce a non-canonical value, then
5439 `mpq_class::canonicalize' must be called before further operations.
5441 -- Function: istream& operator>> (istream& STREAM, mpq_class& ROP);
5442 Read ROP from STREAM, using its `ios' formatting settings, the
5443 same as `mpq_t operator>>' (*note C++ Formatted Input::).
5445 If the ROP read might not be in canonical form then
5446 `mpq_class::canonicalize' must be called.
5449 File: gmp.info, Node: C++ Interface Floats, Next: C++ Interface Random Numbers, Prev: C++ Interface Rationals, Up: C++ Class Interface
5451 12.4 C++ Interface Floats
5452 =========================
5454 When an expression requires the use of temporary intermediate
5455 `mpf_class' values, like `f=g*h+x*y', those temporaries will have the
5456 same precision as the destination `f'. Explicit constructors can be
5457 used if this doesn't suit.
5459 -- Function: mpf_class::mpf_class (type OP)
5460 -- Function: mpf_class::mpf_class (type OP, unsigned long PREC)
5461 Construct an `mpf_class'. Any standard C++ type can be used,
5462 except `long long' and `long double', and any of the GMP C++
5463 classes can be used.
5465 If PREC is given, the initial precision is that value, in bits. If
5466 PREC is not given, then the initial precision is determined by the
5467 type of OP given. An `mpz_class', `mpq_class', or C++ builtin
5468 type will give the default `mpf' precision (*note Initializing
5469 Floats::). An `mpf_class' or expression will give the precision
5470 of that value. The precision of a binary expression is the higher
5471 of the two operands.
5473 mpf_class f(1.5); // default precision
5474 mpf_class f(1.5, 500); // 500 bits (at least)
5475 mpf_class f(x); // precision of x
5476 mpf_class f(abs(x)); // precision of x
5477 mpf_class f(-g, 1000); // 1000 bits (at least)
5478 mpf_class f(x+y); // greater of precisions of x and y
5480 -- Function: void mpf_class::mpf_class (const char *S)
5481 -- Function: void mpf_class::mpf_class (const char *S, unsigned long
5483 -- Function: void mpf_class::mpf_class (const string& S)
5484 -- Function: void mpf_class::mpf_class (const string& S, unsigned long
5486 Construct an `mpf_class' converted from a string using
5487 `mpf_set_str' (*note Assigning Floats::). If PREC is given, the
5488 initial precision is that value, in bits. If not, the default
5489 `mpf' precision (*note Initializing Floats::) is used.
5491 If the string is not a valid float, an `std::invalid_argument'
5492 exception is thrown. The same applies to `operator='.
5494 -- Function: mpf_class& mpf_class::operator= (type OP)
5495 Convert and store the given OP value to an `mpf_class' object. The
5496 same types are accepted as for the constructors above.
5498 Note that `operator=' only stores a new value, it doesn't copy or
5499 change the precision of the destination, instead the value is
5500 truncated if necessary. This is the same as `mpf_set' etc. Note
5501 in particular this means for `mpf_class' a copy constructor is not
5502 the same as a default constructor plus assignment.
5504 mpf_class x (y); // x created with precision of y
5506 mpf_class x; // x created with default precision
5507 x = y; // value truncated to that precision
5509 Applications using templated code may need to be careful about the
5510 assumptions the code makes in this area, when working with
5511 `mpf_class' values of various different or non-default precisions.
5512 For instance implementations of the standard `complex' template
5513 have been seen in both styles above, though of course `complex' is
5514 normally only actually specified for use with the builtin float
5517 -- Function: mpf_class abs (mpf_class OP)
5518 -- Function: mpf_class ceil (mpf_class OP)
5519 -- Function: int cmp (mpf_class OP1, type OP2)
5520 -- Function: int cmp (type OP1, mpf_class OP2)
5521 -- Function: bool mpf_class::fits_sint_p (void)
5522 -- Function: bool mpf_class::fits_slong_p (void)
5523 -- Function: bool mpf_class::fits_sshort_p (void)
5524 -- Function: bool mpf_class::fits_uint_p (void)
5525 -- Function: bool mpf_class::fits_ulong_p (void)
5526 -- Function: bool mpf_class::fits_ushort_p (void)
5527 -- Function: mpf_class floor (mpf_class OP)
5528 -- Function: mpf_class hypot (mpf_class OP1, mpf_class OP2)
5529 -- Function: double mpf_class::get_d (void)
5530 -- Function: long mpf_class::get_si (void)
5531 -- Function: string mpf_class::get_str (mp_exp_t& EXP, int BASE = 10,
5533 -- Function: unsigned long mpf_class::get_ui (void)
5534 -- Function: int mpf_class::set_str (const char *STR, int BASE)
5535 -- Function: int mpf_class::set_str (const string& STR, int BASE)
5536 -- Function: int sgn (mpf_class OP)
5537 -- Function: mpf_class sqrt (mpf_class OP)
5538 -- Function: mpf_class trunc (mpf_class OP)
5539 These functions provide a C++ class interface to the corresponding
5542 `cmp' can be used with any of the classes or the standard C++
5543 types, except `long long' and `long double'.
5545 The accuracy provided by `hypot' is not currently guaranteed.
5547 -- Function: mp_bitcnt_t mpf_class::get_prec ()
5548 -- Function: void mpf_class::set_prec (mp_bitcnt_t PREC)
5549 -- Function: void mpf_class::set_prec_raw (mp_bitcnt_t PREC)
5550 Get or set the current precision of an `mpf_class'.
5552 The restrictions described for `mpf_set_prec_raw' (*note
5553 Initializing Floats::) apply to `mpf_class::set_prec_raw'. Note
5554 in particular that the `mpf_class' must be restored to it's
5555 allocated precision before being destroyed. This must be done by
5556 application code, there's no automatic mechanism for it.
5559 File: gmp.info, Node: C++ Interface Random Numbers, Next: C++ Interface Limitations, Prev: C++ Interface Floats, Up: C++ Class Interface
5561 12.5 C++ Interface Random Numbers
5562 =================================
5564 -- Class: gmp_randclass
5565 The C++ class interface to the GMP random number functions uses
5566 `gmp_randclass' to hold an algorithm selection and current state,
5567 as per `gmp_randstate_t'.
5569 -- Function: gmp_randclass::gmp_randclass (void (*RANDINIT)
5570 (gmp_randstate_t, ...), ...)
5571 Construct a `gmp_randclass', using a call to the given RANDINIT
5572 function (*note Random State Initialization::). The arguments
5573 expected are the same as RANDINIT, but with `mpz_class' instead of
5574 `mpz_t'. For example,
5576 gmp_randclass r1 (gmp_randinit_default);
5577 gmp_randclass r2 (gmp_randinit_lc_2exp_size, 32);
5578 gmp_randclass r3 (gmp_randinit_lc_2exp, a, c, m2exp);
5579 gmp_randclass r4 (gmp_randinit_mt);
5581 `gmp_randinit_lc_2exp_size' will fail if the size requested is too
5582 big, an `std::length_error' exception is thrown in that case.
5584 -- Function: gmp_randclass::gmp_randclass (gmp_randalg_t ALG, ...)
5585 Construct a `gmp_randclass' using the same parameters as
5586 `gmp_randinit' (*note Random State Initialization::). This
5587 function is obsolete and the above RANDINIT style should be
5590 -- Function: void gmp_randclass::seed (unsigned long int S)
5591 -- Function: void gmp_randclass::seed (mpz_class S)
5592 Seed a random number generator. See *note Random Number
5593 Functions::, for how to choose a good seed.
5595 -- Function: mpz_class gmp_randclass::get_z_bits (unsigned long BITS)
5596 -- Function: mpz_class gmp_randclass::get_z_bits (mpz_class BITS)
5597 Generate a random integer with a specified number of bits.
5599 -- Function: mpz_class gmp_randclass::get_z_range (mpz_class N)
5600 Generate a random integer in the range 0 to N-1 inclusive.
5602 -- Function: mpf_class gmp_randclass::get_f ()
5603 -- Function: mpf_class gmp_randclass::get_f (unsigned long PREC)
5604 Generate a random float F in the range 0 <= F < 1. F will be to
5605 PREC bits precision, or if PREC is not given then to the precision
5606 of the destination. For example,
5610 mpf_class f (0, 512); // 512 bits precision
5611 f = r.get_f(); // random number, 512 bits
5614 File: gmp.info, Node: C++ Interface Limitations, Prev: C++ Interface Random Numbers, Up: C++ Class Interface
5616 12.6 C++ Interface Limitations
5617 ==============================
5619 `mpq_class' and Templated Reading
5620 A generic piece of template code probably won't know that
5621 `mpq_class' requires a `canonicalize' call if inputs read with
5622 `operator>>' might be non-canonical. This can lead to incorrect
5625 `operator>>' behaves as it does for reasons of efficiency. A
5626 canonicalize can be quite time consuming on large operands, and is
5627 best avoided if it's not necessary.
5629 But this potential difficulty reduces the usefulness of
5630 `mpq_class'. Perhaps a mechanism to tell `operator>>' what to do
5631 will be adopted in the future, maybe a preprocessor define, a
5632 global flag, or an `ios' flag pressed into service. Or maybe, at
5633 the risk of inconsistency, the `mpq_class' `operator>>' could
5634 canonicalize and leave `mpq_t' `operator>>' not doing so, for use
5635 on those occasions when that's acceptable. Send feedback or
5636 alternate ideas to <gmp-bugs@gmplib.org>.
5639 Subclassing the GMP C++ classes works, but is not currently
5642 Expressions involving subclasses resolve correctly (or seem to),
5643 but in normal C++ fashion the subclass doesn't inherit
5644 constructors and assignments. There's many of those in the GMP
5645 classes, and a good way to reestablish them in a subclass is not
5648 Templated Expressions
5649 A subtle difficulty exists when using expressions together with
5650 application-defined template functions. Consider the following,
5651 with `T' intended to be some numeric type,
5654 T fun (const T &, const T &);
5656 When used with, say, plain `mpz_class' variables, it works fine:
5657 `T' is resolved as `mpz_class'.
5659 mpz_class f(1), g(2);
5662 But when one of the arguments is an expression, it doesn't work.
5664 mpz_class f(1), g(2), h(3);
5665 fun (f, g+h); // Bad
5667 This is because `g+h' ends up being a certain expression template
5668 type internal to `gmpxx.h', which the C++ template resolution
5669 rules are unable to automatically convert to `mpz_class'. The
5670 workaround is simply to add an explicit cast.
5672 mpz_class f(1), g(2), h(3);
5673 fun (f, mpz_class(g+h)); // Good
5675 Similarly, within `fun' it may be necessary to cast an expression
5676 to type `T' when calling a templated `fun2'.
5681 fun2 (f, f+g); // Bad
5687 fun2 (f, T(f+g)); // Good
5691 File: gmp.info, Node: BSD Compatible Functions, Next: Custom Allocation, Prev: C++ Class Interface, Up: Top
5693 13 Berkeley MP Compatible Functions
5694 ***********************************
5696 These functions are intended to be fully compatible with the Berkeley MP
5697 library which is available on many BSD derived U*ix systems. The
5698 `--enable-mpbsd' option must be used when building GNU MP to make these
5699 available (*note Installing GMP::).
5701 The original Berkeley MP library has a usage restriction: you cannot
5702 use the same variable as both source and destination in a single
5703 function call. The compatible functions in GNU MP do not share this
5704 restriction--inputs and outputs may overlap.
5706 It is not recommended that new programs are written using these
5707 functions. Apart from the incomplete set of functions, the interface
5708 for initializing `MINT' objects is more error prone, and the `pow'
5709 function collides with `pow' in `libm.a'.
5711 Include the header `mp.h' to get the definition of the necessary
5712 types and functions. If you are on a BSD derived system, make sure to
5713 include GNU `mp.h' if you are going to link the GNU `libmp.a' to your
5714 program. This means that you probably need to give the `-I<dir>'
5715 option to the compiler, where `<dir>' is the directory where you have
5718 -- Function: MINT * itom (signed short int INITIAL_VALUE)
5719 Allocate an integer consisting of a `MINT' object and dynamic limb
5720 space. Initialize the integer to INITIAL_VALUE. Return a pointer
5721 to the `MINT' object.
5723 -- Function: MINT * xtom (char *INITIAL_VALUE)
5724 Allocate an integer consisting of a `MINT' object and dynamic limb
5725 space. Initialize the integer from INITIAL_VALUE, a hexadecimal,
5726 null-terminated C string. Return a pointer to the `MINT' object.
5728 -- Function: void move (MINT *SRC, MINT *DEST)
5729 Set DEST to SRC by copying. Both variables must be previously
5732 -- Function: void madd (MINT *SRC_1, MINT *SRC_2, MINT *DESTINATION)
5733 Add SRC_1 and SRC_2 and put the sum in DESTINATION.
5735 -- Function: void msub (MINT *SRC_1, MINT *SRC_2, MINT *DESTINATION)
5736 Subtract SRC_2 from SRC_1 and put the difference in DESTINATION.
5738 -- Function: void mult (MINT *SRC_1, MINT *SRC_2, MINT *DESTINATION)
5739 Multiply SRC_1 and SRC_2 and put the product in DESTINATION.
5741 -- Function: void mdiv (MINT *DIVIDEND, MINT *DIVISOR, MINT *QUOTIENT,
5743 -- Function: void sdiv (MINT *DIVIDEND, signed short int DIVISOR, MINT
5744 *QUOTIENT, signed short int *REMAINDER)
5745 Set QUOTIENT to DIVIDEND/DIVISOR, and REMAINDER to DIVIDEND mod
5746 DIVISOR. The quotient is rounded towards zero; the remainder has
5747 the same sign as the dividend unless it is zero.
5749 Some implementations of these functions work differently--or not
5750 at all--for negative arguments.
5752 -- Function: void msqrt (MINT *OP, MINT *ROOT, MINT *REMAINDER)
5753 Set ROOT to the truncated integer part of the square root of OP,
5754 like `mpz_sqrt'. Set REMAINDER to OP-ROOT*ROOT, i.e. zero if OP
5755 is a perfect square.
5757 If ROOT and REMAINDER are the same variable, the results are
5760 -- Function: void pow (MINT *BASE, MINT *EXP, MINT *MOD, MINT *DEST)
5761 Set DEST to (BASE raised to EXP) modulo MOD.
5763 Note that the name `pow' clashes with `pow' from the standard C
5764 math library (*note Exponentiation and Logarithms: (libc)Exponents
5765 and Logarithms.). An application will only be able to use one or
5768 -- Function: void rpow (MINT *BASE, signed short int EXP, MINT *DEST)
5769 Set DEST to BASE raised to EXP.
5771 -- Function: void gcd (MINT *OP1, MINT *OP2, MINT *RES)
5772 Set RES to the greatest common divisor of OP1 and OP2.
5774 -- Function: int mcmp (MINT *OP1, MINT *OP2)
5775 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero
5776 if OP1 = OP2, and a negative value if OP1 < OP2.
5778 -- Function: void min (MINT *DEST)
5779 Input a decimal string from `stdin', and put the read integer in
5780 DEST. SPC and TAB are allowed in the number string, and are
5783 -- Function: void mout (MINT *SRC)
5784 Output SRC to `stdout', as a decimal string. Also output a
5787 -- Function: char * mtox (MINT *OP)
5788 Convert OP to a hexadecimal string, and return a pointer to the
5789 string. The returned string is allocated using the default memory
5790 allocation function, `malloc' by default. It will be
5791 `strlen(str)+1' bytes, that being exactly enough for the string
5792 and null-terminator.
5794 -- Function: void mfree (MINT *OP)
5795 De-allocate, the space used by OP. *This function should only be
5796 passed a value returned by `itom' or `xtom'.*
5799 File: gmp.info, Node: Custom Allocation, Next: Language Bindings, Prev: BSD Compatible Functions, Up: Top
5801 14 Custom Allocation
5802 ********************
5804 By default GMP uses `malloc', `realloc' and `free' for memory
5805 allocation, and if they fail GMP prints a message to the standard error
5806 output and terminates the program.
5808 Alternate functions can be specified, to allocate memory in a
5809 different way or to have a different error action on running out of
5812 This feature is available in the Berkeley compatibility library
5813 (*note BSD Compatible Functions::) as well as the main GMP library.
5815 -- Function: void mp_set_memory_functions (
5816 void *(*ALLOC_FUNC_PTR) (size_t),
5817 void *(*REALLOC_FUNC_PTR) (void *, size_t, size_t),
5818 void (*FREE_FUNC_PTR) (void *, size_t))
5819 Replace the current allocation functions from the arguments. If
5820 an argument is `NULL', the corresponding default function is used.
5822 These functions will be used for all memory allocation done by
5823 GMP, apart from temporary space from `alloca' if that function is
5824 available and GMP is configured to use it (*note Build Options::).
5826 *Be sure to call `mp_set_memory_functions' only when there are no
5827 active GMP objects allocated using the previous memory functions!
5828 Usually that means calling it before any other GMP function.*
5830 The functions supplied should fit the following declarations:
5832 -- Function: void * allocate_function (size_t ALLOC_SIZE)
5833 Return a pointer to newly allocated space with at least ALLOC_SIZE
5836 -- Function: void * reallocate_function (void *PTR, size_t OLD_SIZE,
5838 Resize a previously allocated block PTR of OLD_SIZE bytes to be
5841 The block may be moved if necessary or if desired, and in that
5842 case the smaller of OLD_SIZE and NEW_SIZE bytes must be copied to
5843 the new location. The return value is a pointer to the resized
5844 block, that being the new location if moved or just PTR if not.
5846 PTR is never `NULL', it's always a previously allocated block.
5847 NEW_SIZE may be bigger or smaller than OLD_SIZE.
5849 -- Function: void free_function (void *PTR, size_t SIZE)
5850 De-allocate the space pointed to by PTR.
5852 PTR is never `NULL', it's always a previously allocated block of
5855 A "byte" here means the unit used by the `sizeof' operator.
5857 The OLD_SIZE parameters to REALLOCATE_FUNCTION and FREE_FUNCTION are
5858 passed for convenience, but of course can be ignored if not needed.
5859 The default functions using `malloc' and friends for instance don't use
5862 No error return is allowed from any of these functions, if they
5863 return then they must have performed the specified operation. In
5864 particular note that ALLOCATE_FUNCTION or REALLOCATE_FUNCTION mustn't
5867 Getting a different fatal error action is a good use for custom
5868 allocation functions, for example giving a graphical dialog rather than
5869 the default print to `stderr'. How much is possible when genuinely out
5870 of memory is another question though.
5872 There's currently no defined way for the allocation functions to
5873 recover from an error such as out of memory, they must terminate
5874 program execution. A `longjmp' or throwing a C++ exception will have
5875 undefined results. This may change in the future.
5877 GMP may use allocated blocks to hold pointers to other allocated
5878 blocks. This will limit the assumptions a conservative garbage
5879 collection scheme can make.
5881 Since the default GMP allocation uses `malloc' and friends, those
5882 functions will be linked in even if the first thing a program does is an
5883 `mp_set_memory_functions'. It's necessary to change the GMP sources if
5887 -- Function: void mp_get_memory_functions (
5888 void *(**ALLOC_FUNC_PTR) (size_t),
5889 void *(**REALLOC_FUNC_PTR) (void *, size_t, size_t),
5890 void (**FREE_FUNC_PTR) (void *, size_t))
5891 Get the current allocation functions, storing function pointers to
5892 the locations given by the arguments. If an argument is `NULL',
5893 that function pointer is not stored.
5895 For example, to get just the current free function,
5897 void (*freefunc) (void *, size_t);
5899 mp_get_memory_functions (NULL, NULL, &freefunc);
5902 File: gmp.info, Node: Language Bindings, Next: Algorithms, Prev: Custom Allocation, Up: Top
5904 15 Language Bindings
5905 ********************
5907 The following packages and projects offer access to GMP from languages
5908 other than C, though perhaps with varying levels of functionality and
5913 * GMP C++ class interface, *note C++ Class Interface::
5914 Straightforward interface, expression templates to eliminate
5917 * ALP `http://www-sop.inria.fr/saga/logiciels/ALP/'
5918 Linear algebra and polynomials using templates.
5920 * Arithmos `http://www.win.ua.ac.be/~cant/arithmos/'
5921 Rationals with infinities and square roots.
5923 * CLN `http://www.ginac.de/CLN/'
5924 High level classes for arithmetic.
5926 * LiDIA `http://www.cdc.informatik.tu-darmstadt.de/TI/LiDIA/'
5927 A C++ library for computational number theory.
5929 * Linbox `http://www.linalg.org/'
5930 Sparse vectors and matrices.
5932 * NTL `http://www.shoup.net/ntl/'
5933 A C++ number theory library.
5936 * Omni F77 `http://phase.hpcc.jp/Omni/home.html'
5937 Arbitrary precision floats.
5940 * Glasgow Haskell Compiler `http://www.haskell.org/ghc/'
5943 * Kaffe `http://www.kaffe.org/'
5945 * Kissme `http://kissme.sourceforge.net/'
5948 * GNU Common Lisp `http://www.gnu.org/software/gcl/gcl.html'
5950 * Librep `http://librep.sourceforge.net/'
5952 * XEmacs (21.5.18 beta and up) `http://www.xemacs.org'
5953 Optional big integers, rationals and floats using GMP.
5956 * GNU m4 betas `http://www.seindal.dk/rene/gnu/'
5957 Optionally provides an arbitrary precision `mpeval'.
5960 * MLton compiler `http://mlton.org/'
5963 * MLGMP `http://www.di.ens.fr/~monniaux/programmes.html.en'
5965 * Numerix `http://pauillac.inria.fr/~quercia/'
5966 Optionally using GMP.
5969 * Mozart `http://www.mozart-oz.org/'
5972 * GNU Pascal Compiler `http://www.gnu-pascal.de/'
5975 * Numerix `http://pauillac.inria.fr/~quercia/'
5976 For Free Pascal, optionally using GMP.
5979 * GMP module, see `demos/perl' in the GMP sources (*note
5980 Demonstration Programs::).
5982 * Math::GMP `http://www.cpan.org/'
5983 Compatible with Math::BigInt, but not as many functions as
5984 the GMP module above.
5986 * Math::BigInt::GMP `http://www.cpan.org/'
5987 Plug Math::GMP into normal Math::BigInt operations.
5990 * mpz module in the standard distribution,
5991 `http://pike.ida.liu.se/'
5994 * SWI Prolog `http://www.swi-prolog.org/'
5995 Arbitrary precision floats.
5998 * mpz module in the standard distribution,
5999 `http://www.python.org/'
6001 * GMPY `http://gmpy.sourceforge.net/'
6004 * GNU Guile (upcoming 1.8)
6005 `http://www.gnu.org/software/guile/guile.html'
6007 * RScheme `http://www.rscheme.org/'
6009 * STklos `http://www.stklos.org/'
6013 `http://www.smalltalk.org/versions/GNUSmalltalk.html'
6016 * Axiom `http://savannah.nongnu.org/projects/axiom'
6017 Computer algebra using GCL.
6019 * DrGenius `http://drgenius.seul.org/'
6020 Geometry system and mathematical programming language.
6022 * GiNaC `http://www.ginac.de/'
6023 C++ computer algebra using CLN.
6025 * GOO `http://www.googoogaga.org/'
6026 Dynamic object oriented language.
6028 * Maxima `http://www.ma.utexas.edu/users/wfs/maxima.html'
6029 Macsyma computer algebra using GCL.
6031 * Q `http://q-lang.sourceforge.net/'
6032 Equational programming system.
6034 * Regina `http://regina.sourceforge.net/'
6035 Topological calculator.
6037 * Yacas `http://www.xs4all.nl/~apinkus/yacas.html'
6038 Yet another computer algebra system.
6042 File: gmp.info, Node: Algorithms, Next: Internals, Prev: Language Bindings, Up: Top
6047 This chapter is an introduction to some of the algorithms used for
6048 various GMP operations. The code is likely to be hard to understand
6049 without knowing something about the algorithms.
6051 Some GMP internals are mentioned, but applications that expect to be
6052 compatible with future GMP releases should take care to use only the
6053 documented functions.
6057 * Multiplication Algorithms::
6058 * Division Algorithms::
6059 * Greatest Common Divisor Algorithms::
6060 * Powering Algorithms::
6061 * Root Extraction Algorithms::
6062 * Radix Conversion Algorithms::
6063 * Other Algorithms::
6067 File: gmp.info, Node: Multiplication Algorithms, Next: Division Algorithms, Prev: Algorithms, Up: Algorithms
6072 NxN limb multiplications and squares are done using one of five
6073 algorithms, as the size N increases.
6077 Karatsuba `MUL_TOOM22_THRESHOLD'
6078 Toom-3 `MUL_TOOM33_THRESHOLD'
6079 Toom-4 `MUL_TOOM44_THRESHOLD'
6080 FFT `MUL_FFT_THRESHOLD'
6082 Similarly for squaring, with the `SQR' thresholds.
6084 NxM multiplications of operands with different sizes above
6085 `MUL_TOOM22_THRESHOLD' are currently done by special Toom-inspired
6086 algorithms or directly with FFT, depending on operand size (*note
6087 Unbalanced Multiplication::).
6091 * Basecase Multiplication::
6092 * Karatsuba Multiplication::
6093 * Toom 3-Way Multiplication::
6094 * Toom 4-Way Multiplication::
6095 * FFT Multiplication::
6096 * Other Multiplication::
6097 * Unbalanced Multiplication::
6100 File: gmp.info, Node: Basecase Multiplication, Next: Karatsuba Multiplication, Prev: Multiplication Algorithms, Up: Multiplication Algorithms
6102 16.1.1 Basecase Multiplication
6103 ------------------------------
6105 Basecase NxM multiplication is a straightforward rectangular set of
6106 cross-products, the same as long multiplication done by hand and for
6107 that reason sometimes known as the schoolbook or grammar school method.
6108 This is an O(N*M) algorithm. See Knuth section 4.3.1 algorithm M
6109 (*note References::), and the `mpn/generic/mul_basecase.c' code.
6111 Assembly implementations of `mpn_mul_basecase' are essentially the
6112 same as the generic C code, but have all the usual assembly tricks and
6113 obscurities introduced for speed.
6115 A square can be done in roughly half the time of a multiply, by
6116 using the fact that the cross products above and below the diagonal are
6117 the same. A triangle of products below the diagonal is formed, doubled
6118 (left shift by one bit), and then the products on the diagonal added.
6119 This can be seen in `mpn/generic/sqr_basecase.c'. Again the assembly
6120 implementations take essentially the same approach.
6123 +---+---+---+---+---+
6125 +---+---+---+---+---+
6127 +---+---+---+---+---+
6129 +---+---+---+---+---+
6131 +---+---+---+---+---+
6133 +---+---+---+---+---+
6135 In practice squaring isn't a full 2x faster than multiplying, it's
6136 usually around 1.5x. Less than 1.5x probably indicates
6137 `mpn_sqr_basecase' wants improving on that CPU.
6139 On some CPUs `mpn_mul_basecase' can be faster than the generic C
6140 `mpn_sqr_basecase' on some small sizes. `SQR_BASECASE_THRESHOLD' is
6141 the size at which to use `mpn_sqr_basecase', this will be zero if that
6142 routine should be used always.
6145 File: gmp.info, Node: Karatsuba Multiplication, Next: Toom 3-Way Multiplication, Prev: Basecase Multiplication, Up: Multiplication Algorithms
6147 16.1.2 Karatsuba Multiplication
6148 -------------------------------
6150 The Karatsuba multiplication algorithm is described in Knuth section
6151 4.3.3 part A, and various other textbooks. A brief description is
6154 The inputs x and y are treated as each split into two parts of equal
6155 length (or the most significant part one limb shorter if N is odd).
6158 +----------+----------+
6160 +----------+----------+
6162 +----------+----------+
6164 +----------+----------+
6166 Let b be the power of 2 where the split occurs, ie. if x0 is k limbs
6167 (y0 the same) then b=2^(k*mp_bits_per_limb). With that x=x1*b+x0 and
6168 y=y1*b+y0, and the following holds,
6170 x*y = (b^2+b)*x1*y1 - b*(x1-x0)*(y1-y0) + (b+1)*x0*y0
6172 This formula means doing only three multiplies of (N/2)x(N/2) limbs,
6173 whereas a basecase multiply of NxN limbs is equivalent to four
6174 multiplies of (N/2)x(N/2). The factors (b^2+b) etc represent the
6175 positions where the three products must be added.
6178 +--------+--------+ +--------+--------+
6180 +--------+--------+ +--------+--------+
6188 sub | (x1-x0)*(y1-y0) |
6191 The term (x1-x0)*(y1-y0) is best calculated as an absolute value,
6192 and the sign used to choose to add or subtract. Notice the sum
6193 high(x0*y0)+low(x1*y1) occurs twice, so it's possible to do 5*k limb
6194 additions, rather than 6*k, but in GMP extra function call overheads
6195 outweigh the saving.
6197 Squaring is similar to multiplying, but with x=y the formula reduces
6198 to an equivalent with three squares,
6200 x^2 = (b^2+b)*x1^2 - b*(x1-x0)^2 + (b+1)*x0^2
6202 The final result is accumulated from those three squares the same
6203 way as for the three multiplies above. The middle term (x1-x0)^2 is now
6206 A similar formula for both multiplying and squaring can be
6207 constructed with a middle term (x1+x0)*(y1+y0). But those sums can
6208 exceed k limbs, leading to more carry handling and additions than the
6211 Karatsuba multiplication is asymptotically an O(N^1.585) algorithm,
6212 the exponent being log(3)/log(2), representing 3 multiplies each 1/2
6213 the size of the inputs. This is a big improvement over the basecase
6214 multiply at O(N^2) and the advantage soon overcomes the extra additions
6215 Karatsuba performs. `MUL_TOOM22_THRESHOLD' can be as little as 10
6216 limbs. The `SQR' threshold is usually about twice the `MUL'.
6218 The basecase algorithm will take a time of the form M(N) = a*N^2 +
6219 b*N + c and the Karatsuba algorithm K(N) = 3*M(N/2) + d*N + e, which
6220 expands to K(N) = 3/4*a*N^2 + 3/2*b*N + 3*c + d*N + e. The factor 3/4
6221 for a means per-crossproduct speedups in the basecase code will
6222 increase the threshold since they benefit M(N) more than K(N). And
6223 conversely the 3/2 for b means linear style speedups of b will increase
6224 the threshold since they benefit K(N) more than M(N). The latter can
6225 be seen for instance when adding an optimized `mpn_sqr_diagonal' to
6226 `mpn_sqr_basecase'. Of course all speedups reduce total time, and in
6227 that sense the algorithm thresholds are merely of academic interest.
6230 File: gmp.info, Node: Toom 3-Way Multiplication, Next: Toom 4-Way Multiplication, Prev: Karatsuba Multiplication, Up: Multiplication Algorithms
6232 16.1.3 Toom 3-Way Multiplication
6233 --------------------------------
6235 The Karatsuba formula is the simplest case of a general approach to
6236 splitting inputs that leads to both Toom and FFT algorithms. A
6237 description of Toom can be found in Knuth section 4.3.3, with an
6238 example 3-way calculation after Theorem A. The 3-way form used in GMP
6241 The operands are each considered split into 3 pieces of equal length
6242 (or the most significant part 1 or 2 limbs shorter than the other two).
6245 +----------+----------+----------+
6247 +----------+----------+----------+
6249 +----------+----------+----------+
6251 +----------+----------+----------+
6253 These parts are treated as the coefficients of two polynomials
6255 X(t) = x2*t^2 + x1*t + x0
6256 Y(t) = y2*t^2 + y1*t + y0
6258 Let b equal the power of 2 which is the size of the x0, x1, y0 and
6259 y1 pieces, ie. if they're k limbs each then b=2^(k*mp_bits_per_limb).
6260 With this x=X(b) and y=Y(b).
6262 Let a polynomial W(t)=X(t)*Y(t) and suppose its coefficients are
6264 W(t) = w4*t^4 + w3*t^3 + w2*t^2 + w1*t + w0
6266 The w[i] are going to be determined, and when they are they'll give
6267 the final result using w=W(b), since x*y=X(b)*Y(b)=W(b). The
6268 coefficients will be roughly b^2 each, and the final W(b) will be an
6288 The w[i] coefficients could be formed by a simple set of cross
6289 products, like w4=x2*y2, w3=x2*y1+x1*y2, w2=x2*y0+x1*y1+x0*y2 etc, but
6290 this would need all nine x[i]*y[j] for i,j=0,1,2, and would be
6291 equivalent merely to a basecase multiply. Instead the following
6294 X(t) and Y(t) are evaluated and multiplied at 5 points, giving
6295 values of W(t) at those points. In GMP the following points are used,
6298 t=0 x0 * y0, which gives w0 immediately
6299 t=1 (x2+x1+x0) * (y2+y1+y0)
6300 t=-1 (x2-x1+x0) * (y2-y1+y0)
6301 t=2 (4*x2+2*x1+x0) * (4*y2+2*y1+y0)
6302 t=inf x2 * y2, which gives w4 immediately
6304 At t=-1 the values can be negative and that's handled using the
6305 absolute values and tracking the sign separately. At t=inf the value
6306 is actually X(t)*Y(t)/t^4 in the limit as t approaches infinity, but
6307 it's much easier to think of as simply x2*y2 giving w4 immediately
6308 (much like x0*y0 at t=0 gives w0 immediately).
6310 Each of the points substituted into W(t)=w4*t^4+...+w0 gives a
6311 linear combination of the w[i] coefficients, and the value of those
6312 combinations has just been calculated.
6315 W(1) = w4 + w3 + w2 + w1 + w0
6316 W(-1) = w4 - w3 + w2 - w1 + w0
6317 W(2) = 16*w4 + 8*w3 + 4*w2 + 2*w1 + w0
6320 This is a set of five equations in five unknowns, and some
6321 elementary linear algebra quickly isolates each w[i]. This involves
6322 adding or subtracting one W(t) value from another, and a couple of
6323 divisions by powers of 2 and one division by 3, the latter using the
6324 special `mpn_divexact_by3' (*note Exact Division::).
6326 The conversion of W(t) values to the coefficients is interpolation.
6327 A polynomial of degree 4 like W(t) is uniquely determined by values
6328 known at 5 different points. The points are arbitrary and can be
6329 chosen to make the linear equations come out with a convenient set of
6330 steps for quickly isolating the w[i].
6332 Squaring follows the same procedure as multiplication, but there's
6333 only one X(t) and it's evaluated at the 5 points, and those values
6334 squared to give values of W(t). The interpolation is then identical,
6335 and in fact the same `toom3_interpolate' subroutine is used for both
6336 squaring and multiplying.
6338 Toom-3 is asymptotically O(N^1.465), the exponent being
6339 log(5)/log(3), representing 5 recursive multiplies of 1/3 the original
6340 size each. This is an improvement over Karatsuba at O(N^1.585), though
6341 Toom does more work in the evaluation and interpolation and so it only
6342 realizes its advantage above a certain size.
6344 Near the crossover between Toom-3 and Karatsuba there's generally a
6345 range of sizes where the difference between the two is small.
6346 `MUL_TOOM33_THRESHOLD' is a somewhat arbitrary point in that range and
6347 successive runs of the tune program can give different values due to
6348 small variations in measuring. A graph of time versus size for the two
6349 shows the effect, see `tune/README'.
6351 At the fairly small sizes where the Toom-3 thresholds occur it's
6352 worth remembering that the asymptotic behaviour for Karatsuba and
6353 Toom-3 can't be expected to make accurate predictions, due of course to
6354 the big influence of all sorts of overheads, and the fact that only a
6355 few recursions of each are being performed. Even at large sizes
6356 there's a good chance machine dependent effects like cache architecture
6357 will mean actual performance deviates from what might be predicted.
6359 The formula given for the Karatsuba algorithm (*note Karatsuba
6360 Multiplication::) has an equivalent for Toom-3 involving only five
6361 multiplies, but this would be complicated and unenlightening.
6363 An alternate view of Toom-3 can be found in Zuras (*note
6364 References::), using a vector to represent the x and y splits and a
6365 matrix multiplication for the evaluation and interpolation stages. The
6366 matrix inverses are not meant to be actually used, and they have
6367 elements with values much greater than in fact arise in the
6368 interpolation steps. The diagram shown for the 3-way is attractive,
6369 but again doesn't have to be implemented that way and for example with
6370 a bit of rearrangement just one division by 6 can be done.
6373 File: gmp.info, Node: Toom 4-Way Multiplication, Next: FFT Multiplication, Prev: Toom 3-Way Multiplication, Up: Multiplication Algorithms
6375 16.1.4 Toom 4-Way Multiplication
6376 --------------------------------
6378 Karatsuba and Toom-3 split the operands into 2 and 3 coefficients,
6379 respectively. Toom-4 analogously splits the operands into 4
6380 coefficients. Using the notation from the section on Toom-3
6381 multiplication, we form two polynomials:
6383 X(t) = x3*t^3 + x2*t^2 + x1*t + x0
6384 Y(t) = y3*t^3 + y2*t^2 + y1*t + y0
6386 X(t) and Y(t) are evaluated and multiplied at 7 points, giving
6387 values of W(t) at those points. In GMP the following points are used,
6390 t=0 x0 * y0, which gives w0 immediately
6391 t=1/2 (x3+2*x2+4*x1+8*x0) * (y3+2*y2+4*y1+8*y0)
6392 t=-1/2 (-x3+2*x2-4*x1+8*x0) * (-y3+2*y2-4*y1+8*y0)
6393 t=1 (x3+x2+x1+x0) * (y3+y2+y1+y0)
6394 t=-1 (-x3+x2-x1+x0) * (-y3+y2-y1+y0)
6395 t=2 (8*x3+4*x2+2*x1+x0) * (8*y3+4*y2+2*y1+y0)
6396 t=inf x3 * y3, which gives w6 immediately
6398 The number of additions and subtractions for Toom-4 is much larger
6399 than for Toom-3. But several subexpressions occur multiple times, for
6400 example x2+x0, occurs for both t=1 and t=-1.
6402 Toom-4 is asymptotically O(N^1.404), the exponent being
6403 log(7)/log(4), representing 7 recursive multiplies of 1/4 the original
6407 File: gmp.info, Node: FFT Multiplication, Next: Other Multiplication, Prev: Toom 4-Way Multiplication, Up: Multiplication Algorithms
6409 16.1.5 FFT Multiplication
6410 -------------------------
6412 At large to very large sizes a Fermat style FFT multiplication is used,
6413 following Scho"nhage and Strassen (*note References::). Descriptions
6414 of FFTs in various forms can be found in many textbooks, for instance
6415 Knuth section 4.3.3 part C or Lipson chapter IX. A brief description
6416 of the form used in GMP is given here.
6418 The multiplication done is x*y mod 2^N+1, for a given N. A full
6419 product x*y is obtained by choosing N>=bits(x)+bits(y) and padding x
6420 and y with high zero limbs. The modular product is the native form for
6421 the algorithm, so padding to get a full product is unavoidable.
6423 The algorithm follows a split, evaluate, pointwise multiply,
6424 interpolate and combine similar to that described above for Karatsuba
6425 and Toom-3. A k parameter controls the split, with an FFT-k splitting
6426 into 2^k pieces of M=N/2^k bits each. N must be a multiple of
6427 (2^k)*mp_bits_per_limb so the split falls on limb boundaries, avoiding
6428 bit shifts in the split and combine stages.
6430 The evaluations, pointwise multiplications, and interpolation, are
6431 all done modulo 2^N'+1 where N' is 2M+k+3 rounded up to a multiple of
6432 2^k and of `mp_bits_per_limb'. The results of interpolation will be
6433 the following negacyclic convolution of the input pieces, and the
6434 choice of N' ensures these sums aren't truncated.
6438 w[n] = / (-1) * x[i] * y[j]
6443 The points used for the evaluation are g^i for i=0 to 2^k-1 where
6444 g=2^(2N'/2^k). g is a 2^k'th root of unity mod 2^N'+1, which produces
6445 necessary cancellations at the interpolation stage, and it's also a
6446 power of 2 so the fast Fourier transforms used for the evaluation and
6447 interpolation do only shifts, adds and negations.
6449 The pointwise multiplications are done modulo 2^N'+1 and either
6450 recurse into a further FFT or use a plain multiplication (Toom-3,
6451 Karatsuba or basecase), whichever is optimal at the size N'. The
6452 interpolation is an inverse fast Fourier transform. The resulting set
6453 of sums of x[i]*y[j] are added at appropriate offsets to give the final
6456 Squaring is the same, but x is the only input so it's one transform
6457 at the evaluate stage and the pointwise multiplies are squares. The
6458 interpolation is the same.
6460 For a mod 2^N+1 product, an FFT-k is an O(N^(k/(k-1))) algorithm,
6461 the exponent representing 2^k recursed modular multiplies each
6462 1/2^(k-1) the size of the original. Each successive k is an asymptotic
6463 improvement, but overheads mean each is only faster at bigger and
6464 bigger sizes. In the code, `MUL_FFT_TABLE' and `SQR_FFT_TABLE' are the
6465 thresholds where each k is used. Each new k effectively swaps some
6466 multiplying for some shifts, adds and overheads.
6468 A mod 2^N+1 product can be formed with a normal NxN->2N bit multiply
6469 plus a subtraction, so an FFT and Toom-3 etc can be compared directly.
6470 A k=4 FFT at O(N^1.333) can be expected to be the first faster than
6471 Toom-3 at O(N^1.465). In practice this is what's found, with
6472 `MUL_FFT_MODF_THRESHOLD' and `SQR_FFT_MODF_THRESHOLD' being between 300
6473 and 1000 limbs, depending on the CPU. So far it's been found that only
6474 very large FFTs recurse into pointwise multiplies above these sizes.
6476 When an FFT is to give a full product, the change of N to 2N doesn't
6477 alter the theoretical complexity for a given k, but for the purposes of
6478 considering where an FFT might be first used it can be assumed that the
6479 FFT is recursing into a normal multiply and that on that basis it's
6480 doing 2^k recursed multiplies each 1/2^(k-2) the size of the inputs,
6481 making it O(N^(k/(k-2))). This would mean k=7 at O(N^1.4) would be the
6482 first FFT faster than Toom-3. In practice `MUL_FFT_THRESHOLD' and
6483 `SQR_FFT_THRESHOLD' have been found to be in the k=8 range, somewhere
6484 between 3000 and 10000 limbs.
6486 The way N is split into 2^k pieces and then 2M+k+3 is rounded up to
6487 a multiple of 2^k and `mp_bits_per_limb' means that when
6488 2^k>=mp_bits_per_limb the effective N is a multiple of 2^(2k-1) bits.
6489 The +k+3 means some values of N just under such a multiple will be
6490 rounded to the next. The complexity calculations above assume that a
6491 favourable size is used, meaning one which isn't padded through
6492 rounding, and it's also assumed that the extra +k+3 bits are negligible
6493 at typical FFT sizes.
6495 The practical effect of the 2^(2k-1) constraint is to introduce a
6496 step-effect into measured speeds. For example k=8 will round N up to a
6497 multiple of 32768 bits, so for a 32-bit limb there'll be 512 limb
6498 groups of sizes for which `mpn_mul_n' runs at the same speed. Or for
6499 k=9 groups of 2048 limbs, k=10 groups of 8192 limbs, etc. In practice
6500 it's been found each k is used at quite small multiples of its size
6501 constraint and so the step effect is quite noticeable in a time versus
6504 The threshold determinations currently measure at the mid-points of
6505 size steps, but this is sub-optimal since at the start of a new step it
6506 can happen that it's better to go back to the previous k for a while.
6507 Something more sophisticated for `MUL_FFT_TABLE' and `SQR_FFT_TABLE'
6511 File: gmp.info, Node: Other Multiplication, Next: Unbalanced Multiplication, Prev: FFT Multiplication, Up: Multiplication Algorithms
6513 16.1.6 Other Multiplication
6514 ---------------------------
6516 The Toom algorithms described above (*note Toom 3-Way Multiplication::,
6517 *note Toom 4-Way Multiplication::) generalizes to split into an
6518 arbitrary number of pieces, as per Knuth section 4.3.3 algorithm C.
6519 This is not currently used. The notes here are merely for interest.
6521 In general a split into r+1 pieces is made, and evaluations and
6522 pointwise multiplications done at 2*r+1 points. A 4-way split does 7
6523 pointwise multiplies, 5-way does 9, etc. Asymptotically an (r+1)-way
6524 algorithm is O(N^(log(2*r+1)/log(r+1))). Only the pointwise
6525 multiplications count towards big-O complexity, but the time spent in
6526 the evaluate and interpolate stages grows with r and has a significant
6527 practical impact, with the asymptotic advantage of each r realized only
6528 at bigger and bigger sizes. The overheads grow as O(N*r), whereas in
6529 an r=2^k FFT they grow only as O(N*log(r)).
6531 Knuth algorithm C evaluates at points 0,1,2,...,2*r, but exercise 4
6532 uses -r,...,0,...,r and the latter saves some small multiplies in the
6533 evaluate stage (or rather trades them for additions), and has a further
6534 saving of nearly half the interpolate steps. The idea is to separate
6535 odd and even final coefficients and then perform algorithm C steps C7
6536 and C8 on them separately. The divisors at step C7 become j^2 and the
6537 multipliers at C8 become 2*t*j-j^2.
6539 Splitting odd and even parts through positive and negative points
6540 can be thought of as using -1 as a square root of unity. If a 4th root
6541 of unity was available then a further split and speedup would be
6542 possible, but no such root exists for plain integers. Going to complex
6543 integers with i=sqrt(-1) doesn't help, essentially because in Cartesian
6544 form it takes three real multiplies to do a complex multiply. The
6545 existence of 2^k'th roots of unity in a suitable ring or field lets the
6546 fast Fourier transform keep splitting and get to O(N*log(r)).
6548 Floating point FFTs use complex numbers approximating Nth roots of
6549 unity. Some processors have special support for such FFTs. But these
6550 are not used in GMP since it's very difficult to guarantee an exact
6551 result (to some number of bits). An occasional difference of 1 in the
6552 last bit might not matter to a typical signal processing algorithm, but
6553 is of course of vital importance to GMP.
6556 File: gmp.info, Node: Unbalanced Multiplication, Prev: Other Multiplication, Up: Multiplication Algorithms
6558 16.1.7 Unbalanced Multiplication
6559 --------------------------------
6561 Multiplication of operands with different sizes, both below
6562 `MUL_TOOM22_THRESHOLD' are done with plain schoolbook multiplication
6563 (*note Basecase Multiplication::).
6565 For really large operands, we invoke FFT directly.
6567 For operands between these sizes, we use Toom inspired algorithms
6568 suggested by Alberto Zanoni and Marco Bodrato. The idea is to split
6569 the operands into polynomials of different degree. GMP currently
6570 splits the smaller operand onto 2 coefficients, i.e., a polynomial of
6571 degree 1, but the larger operand can be split into 2, 3, or 4
6572 coefficients, i.e., a polynomial of degree 1 to 3.
6575 File: gmp.info, Node: Division Algorithms, Next: Greatest Common Divisor Algorithms, Prev: Multiplication Algorithms, Up: Algorithms
6577 16.2 Division Algorithms
6578 ========================
6582 * Single Limb Division::
6583 * Basecase Division::
6584 * Divide and Conquer Division::
6585 * Block-Wise Barrett Division::
6588 * Small Quotient Division::
6591 File: gmp.info, Node: Single Limb Division, Next: Basecase Division, Prev: Division Algorithms, Up: Division Algorithms
6593 16.2.1 Single Limb Division
6594 ---------------------------
6596 Nx1 division is implemented using repeated 2x1 divisions from high to
6597 low, either with a hardware divide instruction or a multiplication by
6598 inverse, whichever is best on a given CPU.
6600 The multiply by inverse follows "Improved division by invariant
6601 integers" by Mo"ller and Granlund (*note References::) and is
6602 implemented as `udiv_qrnnd_preinv' in `gmp-impl.h'. The idea is to
6603 have a fixed-point approximation to 1/d (see `invert_limb') and then
6604 multiply by the high limb (plus one bit) of the dividend to get a
6605 quotient q. With d normalized (high bit set), q is no more than 1 too
6606 small. Subtracting q*d from the dividend gives a remainder, and
6607 reveals whether q or q-1 is correct.
6609 The result is a division done with two multiplications and four or
6610 five arithmetic operations. On CPUs with low latency multipliers this
6611 can be much faster than a hardware divide, though the cost of
6612 calculating the inverse at the start may mean it's only better on
6613 inputs bigger than say 4 or 5 limbs.
6615 When a divisor must be normalized, either for the generic C
6616 `__udiv_qrnnd_c' or the multiply by inverse, the division performed is
6617 actually a*2^k by d*2^k where a is the dividend and k is the power
6618 necessary to have the high bit of d*2^k set. The bit shifts for the
6619 dividend are usually accomplished "on the fly" meaning by extracting
6620 the appropriate bits at each step. Done this way the quotient limbs
6621 come out aligned ready to store. When only the remainder is wanted, an
6622 alternative is to take the dividend limbs unshifted and calculate r = a
6623 mod d*2^k followed by an extra final step r*2^k mod d*2^k. This can
6624 help on CPUs with poor bit shifts or few registers.
6626 The multiply by inverse can be done two limbs at a time. The
6627 calculation is basically the same, but the inverse is two limbs and the
6628 divisor treated as if padded with a low zero limb. This means more
6629 work, since the inverse will need a 2x2 multiply, but the four 1x1s to
6630 do that are independent and can therefore be done partly or wholly in
6631 parallel. Likewise for a 2x1 calculating q*d. The net effect is to
6632 process two limbs with roughly the same two multiplies worth of latency
6633 that one limb at a time gives. This extends to 3 or 4 limbs at a time,
6634 though the extra work to apply the inverse will almost certainly soon
6635 reach the limits of multiplier throughput.
6637 A similar approach in reverse can be taken to process just half a
6638 limb at a time if the divisor is only a half limb. In this case the
6639 1x1 multiply for the inverse effectively becomes two (1/2)x1 for each
6640 limb, which can be a saving on CPUs with a fast half limb multiply, or
6641 in fact if the only multiply is a half limb, and especially if it's not
6645 File: gmp.info, Node: Basecase Division, Next: Divide and Conquer Division, Prev: Single Limb Division, Up: Division Algorithms
6647 16.2.2 Basecase Division
6648 ------------------------
6650 Basecase NxM division is like long division done by hand, but in base
6651 2^mp_bits_per_limb. See Knuth section 4.3.1 algorithm D, and
6652 `mpn/generic/sb_divrem_mn.c'.
6654 Briefly stated, while the dividend remains larger than the divisor,
6655 a high quotient limb is formed and the Nx1 product q*d subtracted at
6656 the top end of the dividend. With a normalized divisor (most
6657 significant bit set), each quotient limb can be formed with a 2x1
6658 division and a 1x1 multiplication plus some subtractions. The 2x1
6659 division is by the high limb of the divisor and is done either with a
6660 hardware divide or a multiply by inverse (the same as in *Note Single
6661 Limb Division::) whichever is faster. Such a quotient is sometimes one
6662 too big, requiring an addback of the divisor, but that happens rarely.
6664 With Q=N-M being the number of quotient limbs, this is an O(Q*M)
6665 algorithm and will run at a speed similar to a basecase QxM
6666 multiplication, differing in fact only in the extra multiply and divide
6667 for each of the Q quotient limbs.
6670 File: gmp.info, Node: Divide and Conquer Division, Next: Block-Wise Barrett Division, Prev: Basecase Division, Up: Division Algorithms
6672 16.2.3 Divide and Conquer Division
6673 ----------------------------------
6675 For divisors larger than `DC_DIV_QR_THRESHOLD', division is done by
6676 dividing. Or to be precise by a recursive divide and conquer algorithm
6677 based on work by Moenck and Borodin, Jebelean, and Burnikel and Ziegler
6678 (*note References::).
6680 The algorithm consists essentially of recognising that a 2NxN
6681 division can be done with the basecase division algorithm (*note
6682 Basecase Division::), but using N/2 limbs as a base, not just a single
6683 limb. This way the multiplications that arise are (N/2)x(N/2) and can
6684 take advantage of Karatsuba and higher multiplication algorithms (*note
6685 Multiplication Algorithms::). The two "digits" of the quotient are
6686 formed by recursive Nx(N/2) divisions.
6688 If the (N/2)x(N/2) multiplies are done with a basecase multiplication
6689 then the work is about the same as a basecase division, but with more
6690 function call overheads and with some subtractions separated from the
6691 multiplies. These overheads mean that it's only when N/2 is above
6692 `MUL_TOOM22_THRESHOLD' that divide and conquer is of use.
6694 `DC_DIV_QR_THRESHOLD' is based on the divisor size N, so it will be
6695 somewhere above twice `MUL_TOOM22_THRESHOLD', but how much above
6696 depends on the CPU. An optimized `mpn_mul_basecase' can lower
6697 `DC_DIV_QR_THRESHOLD' a little by offering a ready-made advantage over
6698 repeated `mpn_submul_1' calls.
6700 Divide and conquer is asymptotically O(M(N)*log(N)) where M(N) is
6701 the time for an NxN multiplication done with FFTs. The actual time is
6702 a sum over multiplications of the recursed sizes, as can be seen near
6703 the end of section 2.2 of Burnikel and Ziegler. For example, within
6704 the Toom-3 range, divide and conquer is 2.63*M(N). With higher
6705 algorithms the M(N) term improves and the multiplier tends to log(N).
6706 In practice, at moderate to large sizes, a 2NxN division is about 2 to
6707 4 times slower than an NxN multiplication.
6710 File: gmp.info, Node: Block-Wise Barrett Division, Next: Exact Division, Prev: Divide and Conquer Division, Up: Division Algorithms
6712 16.2.4 Block-Wise Barrett Division
6713 ----------------------------------
6715 For the largest divisions, a block-wise Barrett division algorithm is
6716 used. Here, the divisor is inverted to a precision determined by the
6717 relative size of the dividend and divisor. Blocks of quotient limbs
6718 are then generated by multiplying blocks from the dividend by the
6721 Our block-wise algorithm computes a smaller inverse than in the
6722 plain Barrett algorithm. For a 2n/n division, the inverse will be just
6726 File: gmp.info, Node: Exact Division, Next: Exact Remainder, Prev: Block-Wise Barrett Division, Up: Division Algorithms
6728 16.2.5 Exact Division
6729 ---------------------
6731 A so-called exact division is when the dividend is known to be an exact
6732 multiple of the divisor. Jebelean's exact division algorithm uses this
6733 knowledge to make some significant optimizations (*note References::).
6735 The idea can be illustrated in decimal for example with 368154
6736 divided by 543. Because the low digit of the dividend is 4, the low
6737 digit of the quotient must be 8. This is arrived at from 4*7 mod 10,
6738 using the fact 7 is the modular inverse of 3 (the low digit of the
6739 divisor), since 3*7 == 1 mod 10. So 8*543=4344 can be subtracted from
6740 the dividend leaving 363810. Notice the low digit has become zero.
6742 The procedure is repeated at the second digit, with the next
6743 quotient digit 7 (7 == 1*7 mod 10), subtracting 7*543=3801, leaving
6744 325800. And finally at the third digit with quotient digit 6 (8*7 mod
6745 10), subtracting 6*543=3258 leaving 0. So the quotient is 678.
6747 Notice however that the multiplies and subtractions don't need to
6748 extend past the low three digits of the dividend, since that's enough
6749 to determine the three quotient digits. For the last quotient digit no
6750 subtraction is needed at all. On a 2NxN division like this one, only
6751 about half the work of a normal basecase division is necessary.
6753 For an NxM exact division producing Q=N-M quotient limbs, the saving
6754 over a normal basecase division is in two parts. Firstly, each of the
6755 Q quotient limbs needs only one multiply, not a 2x1 divide and
6756 multiply. Secondly, the crossproducts are reduced when Q>M to
6757 Q*M-M*(M+1)/2, or when Q<=M to Q*(Q-1)/2. Notice the savings are
6758 complementary. If Q is big then many divisions are saved, or if Q is
6759 small then the crossproducts reduce to a small number.
6761 The modular inverse used is calculated efficiently by `binvert_limb'
6762 in `gmp-impl.h'. This does four multiplies for a 32-bit limb, or six
6763 for a 64-bit limb. `tune/modlinv.c' has some alternate implementations
6764 that might suit processors better at bit twiddling than multiplying.
6766 The sub-quadratic exact division described by Jebelean in "Exact
6767 Division with Karatsuba Complexity" is not currently implemented. It
6768 uses a rearrangement similar to the divide and conquer for normal
6769 division (*note Divide and Conquer Division::), but operating from low
6770 to high. A further possibility not currently implemented is
6771 "Bidirectional Exact Integer Division" by Krandick and Jebelean which
6772 forms quotient limbs from both the high and low ends of the dividend,
6773 and can halve once more the number of crossproducts needed in a 2NxN
6776 A special case exact division by 3 exists in `mpn_divexact_by3',
6777 supporting Toom-3 multiplication and `mpq' canonicalizations. It forms
6778 quotient digits with a multiply by the modular inverse of 3 (which is
6779 `0xAA..AAB') and uses two comparisons to determine a borrow for the next
6780 limb. The multiplications don't need to be on the dependent chain, as
6781 long as the effect of the borrows is applied, which can help chips with
6782 pipelined multipliers.
6785 File: gmp.info, Node: Exact Remainder, Next: Small Quotient Division, Prev: Exact Division, Up: Division Algorithms
6787 16.2.6 Exact Remainder
6788 ----------------------
6790 If the exact division algorithm is done with a full subtraction at each
6791 stage and the dividend isn't a multiple of the divisor, then low zero
6792 limbs are produced but with a remainder in the high limbs. For
6793 dividend a, divisor d, quotient q, and b = 2^mp_bits_per_limb, this
6794 remainder r is of the form
6798 n represents the number of zero limbs produced by the subtractions,
6799 that being the number of limbs produced for q. r will be in the range
6800 0<=r<d and can be viewed as a remainder, but one shifted up by a factor
6803 Carrying out full subtractions at each stage means the same number
6804 of cross products must be done as a normal division, but there's still
6805 some single limb divisions saved. When d is a single limb some
6806 simplifications arise, providing good speedups on a number of
6809 `mpn_divexact_by3', `mpn_modexact_1_odd' and the `mpn_redc_X'
6810 functions differ subtly in how they return r, leading to some negations
6811 in the above formula, but all are essentially the same.
6813 Clearly r is zero when a is a multiple of d, and this leads to
6814 divisibility or congruence tests which are potentially more efficient
6815 than a normal division.
6817 The factor of b^n on r can be ignored in a GCD when d is odd, hence
6818 the use of `mpn_modexact_1_odd' by `mpn_gcd_1' and `mpz_kronecker_ui'
6819 etc (*note Greatest Common Divisor Algorithms::).
6821 Montgomery's REDC method for modular multiplications uses operands
6822 of the form of x*b^-n and y*b^-n and on calculating (x*b^-n)*(y*b^-n)
6823 uses the factor of b^n in the exact remainder to reach a product in the
6824 same form (x*y)*b^-n (*note Modular Powering Algorithm::).
6826 Notice that r generally gives no useful information about the
6827 ordinary remainder a mod d since b^n mod d could be anything. If
6828 however b^n == 1 mod d, then r is the negative of the ordinary
6829 remainder. This occurs whenever d is a factor of b^n-1, as for example
6830 with 3 in `mpn_divexact_by3'. For a 32 or 64 bit limb other such
6831 factors include 5, 17 and 257, but no particular use has been found for
6835 File: gmp.info, Node: Small Quotient Division, Prev: Exact Remainder, Up: Division Algorithms
6837 16.2.7 Small Quotient Division
6838 ------------------------------
6840 An NxM division where the number of quotient limbs Q=N-M is small can
6841 be optimized somewhat.
6843 An ordinary basecase division normalizes the divisor by shifting it
6844 to make the high bit set, shifting the dividend accordingly, and
6845 shifting the remainder back down at the end of the calculation. This
6846 is wasteful if only a few quotient limbs are to be formed. Instead a
6847 division of just the top 2*Q limbs of the dividend by the top Q limbs
6848 of the divisor can be used to form a trial quotient. This requires
6849 only those limbs normalized, not the whole of the divisor and dividend.
6851 A multiply and subtract then applies the trial quotient to the M-Q
6852 unused limbs of the divisor and N-Q dividend limbs (which includes Q
6853 limbs remaining from the trial quotient division). The starting trial
6854 quotient can be 1 or 2 too big, but all cases of 2 too big and most
6855 cases of 1 too big are detected by first comparing the most significant
6856 limbs that will arise from the subtraction. An addback is done if the
6857 quotient still turns out to be 1 too big.
6859 This whole procedure is essentially the same as one step of the
6860 basecase algorithm done in a Q limb base, though with the trial
6861 quotient test done only with the high limbs, not an entire Q limb
6862 "digit" product. The correctness of this weaker test can be
6863 established by following the argument of Knuth section 4.3.1 exercise
6864 20 but with the v2*q>b*r+u2 condition appropriately relaxed.
6867 File: gmp.info, Node: Greatest Common Divisor Algorithms, Next: Powering Algorithms, Prev: Division Algorithms, Up: Algorithms
6869 16.3 Greatest Common Divisor
6870 ============================
6875 * Lehmer's Algorithm::
6876 * Subquadratic GCD::
6881 File: gmp.info, Node: Binary GCD, Next: Lehmer's Algorithm, Prev: Greatest Common Divisor Algorithms, Up: Greatest Common Divisor Algorithms
6886 At small sizes GMP uses an O(N^2) binary style GCD. This is described
6887 in many textbooks, for example Knuth section 4.5.2 algorithm B. It
6888 simply consists of successively reducing odd operands a and b using
6890 a,b = abs(a-b),min(a,b)
6891 strip factors of 2 from a
6893 The Euclidean GCD algorithm, as per Knuth algorithms E and A,
6894 repeatedly computes the quotient q = floor(a/b) and replaces a,b by v,
6895 u - q v. The binary algorithm has so far been found to be faster than
6896 the Euclidean algorithm everywhere. One reason the binary method does
6897 well is that the implied quotient at each step is usually small, so
6898 often only one or two subtractions are needed to get the same effect as
6899 a division. Quotients 1, 2 and 3 for example occur 67.7% of the time,
6900 see Knuth section 4.5.3 Theorem E.
6902 When the implied quotient is large, meaning b is much smaller than
6903 a, then a division is worthwhile. This is the basis for the initial a
6904 mod b reductions in `mpn_gcd' and `mpn_gcd_1' (the latter for both Nx1
6905 and 1x1 cases). But after that initial reduction, big quotients occur
6906 too rarely to make it worth checking for them.
6909 The final 1x1 GCD in `mpn_gcd_1' is done in the generic C code as
6910 described above. For two N-bit operands, the algorithm takes about
6911 0.68 iterations per bit. For optimum performance some attention needs
6912 to be paid to the way the factors of 2 are stripped from a.
6914 Firstly it may be noted that in twos complement the number of low
6915 zero bits on a-b is the same as b-a, so counting or testing can begin on
6916 a-b without waiting for abs(a-b) to be determined.
6918 A loop stripping low zero bits tends not to branch predict well,
6919 since the condition is data dependent. But on average there's only a
6920 few low zeros, so an option is to strip one or two bits arithmetically
6921 then loop for more (as done for AMD K6). Or use a lookup table to get
6922 a count for several bits then loop for more (as done for AMD K7). An
6923 alternative approach is to keep just one of a or b odd and iterate
6925 a,b = abs(a-b), min(a,b)
6929 This requires about 1.25 iterations per bit, but stripping of a
6930 single bit at each step avoids any branching. Repeating the bit strip
6931 reduces to about 0.9 iterations per bit, which may be a worthwhile
6934 Generally with the above approaches a speed of perhaps 6 cycles per
6935 bit can be achieved, which is still not terribly fast with for instance
6936 a 64-bit GCD taking nearly 400 cycles. It's this sort of time which
6937 means it's not usually advantageous to combine a set of divisibility
6940 Currently, the binary algorithm is used for GCD only when N < 3.
6943 File: gmp.info, Node: Lehmer's Algorithm, Next: Subquadratic GCD, Prev: Binary GCD, Up: Greatest Common Divisor Algorithms
6945 16.3.2 Lehmer's algorithm
6946 -------------------------
6948 Lehmer's improvement of the Euclidean algorithms is based on the
6949 observation that the initial part of the quotient sequence depends only
6950 on the most significant parts of the inputs. The variant of Lehmer's
6951 algorithm used in GMP splits off the most significant two limbs, as
6952 suggested, e.g., in "A Double-Digit Lehmer-Euclid Algorithm" by
6953 Jebelean (*note References::). The quotients of two double-limb inputs
6954 are collected as a 2 by 2 matrix with single-limb elements. This is
6955 done by the function `mpn_hgcd2'. The resulting matrix is applied to
6956 the inputs using `mpn_mul_1' and `mpn_submul_1'. Each iteration usually
6957 reduces the inputs by almost one limb. In the rare case of a large
6958 quotient, no progress can be made by examining just the most
6959 significant two limbs, and the quotient is computing using plain
6962 The resulting algorithm is asymptotically O(N^2), just as the
6963 Euclidean algorithm and the binary algorithm. The quadratic part of the
6964 work are the calls to `mpn_mul_1' and `mpn_submul_1'. For small sizes,
6965 the linear work is also significant. There are roughly N calls to the
6966 `mpn_hgcd2' function. This function uses a couple of important
6969 * It uses the same relaxed notion of correctness as `mpn_hgcd' (see
6970 next section). This means that when called with the most
6971 significant two limbs of two large numbers, the returned matrix
6972 does not always correspond exactly to the initial quotient
6973 sequence for the two large numbers; the final quotient may
6974 sometimes be one off.
6976 * It takes advantage of the fact the quotients are usually small.
6977 The division operator is not used, since the corresponding
6978 assembler instruction is very slow on most architectures. (This
6979 code could probably be improved further, it uses many branches
6980 that are unfriendly to prediction).
6982 * It switches from double-limb calculations to single-limb
6983 calculations half-way through, when the input numbers have been
6984 reduced in size from two limbs to one and a half.
6988 File: gmp.info, Node: Subquadratic GCD, Next: Extended GCD, Prev: Lehmer's Algorithm, Up: Greatest Common Divisor Algorithms
6990 16.3.3 Subquadratic GCD
6991 -----------------------
6993 For inputs larger than `GCD_DC_THRESHOLD', GCD is computed via the HGCD
6994 (Half GCD) function, as a generalization to Lehmer's algorithm.
6996 Let the inputs a,b be of size N limbs each. Put S = floor(N/2) + 1.
6997 Then HGCD(a,b) returns a transformation matrix T with non-negative
6998 elements, and reduced numbers (c;d) = T^-1 (a;b). The reduced numbers
6999 c,d must be larger than S limbs, while their difference abs(c-d) must
7000 fit in S limbs. The matrix elements will also be of size roughly N/2.
7002 The HGCD base case uses Lehmer's algorithm, but with the above stop
7003 condition that returns reduced numbers and the corresponding
7004 transformation matrix half-way through. For inputs larger than
7005 `HGCD_THRESHOLD', HGCD is computed recursively, using the divide and
7006 conquer algorithm in "On Scho"nhage's algorithm and subquadratic
7007 integer GCD computation" by Mo"ller (*note References::). The recursive
7008 algorithm consists of these main steps.
7010 * Call HGCD recursively, on the most significant N/2 limbs. Apply the
7011 resulting matrix T_1 to the full numbers, reducing them to a size
7014 * Perform a small number of division or subtraction steps to reduce
7015 the numbers to size below 3N/2. This is essential mainly for the
7016 unlikely case of large quotients.
7018 * Call HGCD recursively, on the most significant N/2 limbs of the
7019 reduced numbers. Apply the resulting matrix T_2 to the full
7020 numbers, reducing them to a size just above N/2.
7022 * Compute T = T_1 T_2.
7024 * Perform a small number of division and subtraction steps to
7025 satisfy the requirements, and return.
7027 GCD is then implemented as a loop around HGCD, similarly to Lehmer's
7028 algorithm. Where Lehmer repeatedly chops off the top two limbs, calls
7029 `mpn_hgcd2', and applies the resulting matrix to the full numbers, the
7030 subquadratic GCD chops off the most significant third of the limbs (the
7031 proportion is a tuning parameter, and 1/3 seems to be more efficient
7032 than, e.g, 1/2), calls `mpn_hgcd', and applies the resulting matrix.
7033 Once the input numbers are reduced to size below `GCD_DC_THRESHOLD',
7034 Lehmer's algorithm is used for the rest of the work.
7036 The asymptotic running time of both HGCD and GCD is O(M(N)*log(N)),
7037 where M(N) is the time for multiplying two N-limb numbers.
7040 File: gmp.info, Node: Extended GCD, Next: Jacobi Symbol, Prev: Subquadratic GCD, Up: Greatest Common Divisor Algorithms
7045 The extended GCD function, or GCDEXT, calculates gcd(a,b) and also
7046 cofactors x and y satisfying a*x+b*y=gcd(a,b). All the algorithms used
7047 for plain GCD are extended to handle this case. The binary algorithm is
7048 used only for single-limb GCDEXT. Lehmer's algorithm is used for sizes
7049 up to `GCDEXT_DC_THRESHOLD'. Above this threshold, GCDEXT is
7050 implemented as a loop around HGCD, but with more book-keeping to keep
7051 track of the cofactors. This gives the same asymptotic running time as
7052 for GCD and HGCD, O(M(N)*log(N))
7054 One difference to plain GCD is that while the inputs a and b are
7055 reduced as the algorithm proceeds, the cofactors x and y grow in size.
7056 This makes the tuning of the chopping-point more difficult. The current
7057 code chops off the most significant half of the inputs for the call to
7058 HGCD in the first iteration, and the most significant two thirds for
7059 the remaining calls. This strategy could surely be improved. Also the
7060 stop condition for the loop, where Lehmer's algorithm is invoked once
7061 the inputs are reduced below `GCDEXT_DC_THRESHOLD', could maybe be
7062 improved by taking into account the current size of the cofactors.
7065 File: gmp.info, Node: Jacobi Symbol, Prev: Extended GCD, Up: Greatest Common Divisor Algorithms
7067 16.3.5 Jacobi Symbol
7068 --------------------
7070 `mpz_jacobi' and `mpz_kronecker' are currently implemented with a
7071 simple binary algorithm similar to that described for the GCDs (*note
7072 Binary GCD::). They're not very fast when both inputs are large.
7073 Lehmer's multi-step improvement or a binary based multi-step algorithm
7074 is likely to be better.
7076 When one operand fits a single limb, and that includes
7077 `mpz_kronecker_ui' and friends, an initial reduction is done with
7078 either `mpn_mod_1' or `mpn_modexact_1_odd', followed by the binary
7079 algorithm on a single limb. The binary algorithm is well suited to a
7080 single limb, and the whole calculation in this case is quite efficient.
7082 In all the routines sign changes for the result are accumulated
7083 using some bit twiddling, avoiding table lookups or conditional jumps.
projects / xonotic / xonotic.git / blob