1 This is gmp.info, produced by makeinfo version 6.7 from gmp.texi.
3 This manual describes how to install and use the GNU multiple precision
4 arithmetic library, version 6.2.1.
6 Copyright 1991, 1993-2016, 2018-2020 Free Software Foundation, Inc.
8 Permission is granted to copy, distribute and/or modify this document
9 under the terms of the GNU Free Documentation License, Version 1.3 or
10 any later version published by the Free Software Foundation; with no
11 Invariant Sections, with the Front-Cover Texts being "A GNU Manual", and
12 with the Back-Cover Texts being "You have freedom to copy and modify
13 this GNU Manual, like GNU software". A copy of the license is included
14 in *note GNU Free Documentation License::.
15 INFO-DIR-SECTION GNU libraries
17 * gmp: (gmp). GNU Multiple Precision Arithmetic Library.
21 File: gmp.info, Node: Top, Next: Copying, Prev: (dir), Up: (dir)
26 This manual describes how to install and use the GNU multiple precision
27 arithmetic library, version 6.2.1.
29 Copyright 1991, 1993-2016, 2018-2020 Free Software Foundation, Inc.
31 Permission is granted to copy, distribute and/or modify this document
32 under the terms of the GNU Free Documentation License, Version 1.3 or
33 any later version published by the Free Software Foundation; with no
34 Invariant Sections, with the Front-Cover Texts being "A GNU Manual", and
35 with the Back-Cover Texts being "You have freedom to copy and modify
36 this GNU Manual, like GNU software". A copy of the license is included
37 in *note GNU Free Documentation License::.
41 * Copying:: GMP Copying Conditions (LGPL).
42 * Introduction to GMP:: Brief introduction to GNU MP.
43 * Installing GMP:: How to configure and compile the GMP library.
44 * GMP Basics:: What every GMP user should know.
45 * Reporting Bugs:: How to usefully report bugs.
46 * Integer Functions:: Functions for arithmetic on signed integers.
47 * Rational Number Functions:: Functions for arithmetic on rational numbers.
48 * Floating-point Functions:: Functions for arithmetic on floats.
49 * Low-level Functions:: Fast functions for natural numbers.
50 * Random Number Functions:: Functions for generating random numbers.
51 * Formatted Output:: 'printf' style output.
52 * Formatted Input:: 'scanf' style input.
53 * C++ Class Interface:: Class wrappers around GMP types.
54 * Custom Allocation:: How to customize the internal allocation.
55 * Language Bindings:: Using GMP from other languages.
56 * Algorithms:: What happens behind the scenes.
57 * Internals:: How values are represented behind the scenes.
59 * Contributors:: Who brings you this library?
60 * References:: Some useful papers and books to read.
61 * GNU Free Documentation License::
66 File: gmp.info, Node: Copying, Next: Introduction to GMP, Prev: Top, Up: Top
68 GNU MP Copying Conditions
69 *************************
71 This library is "free"; this means that everyone is free to use it and
72 free to redistribute it on a free basis. The library is not in the
73 public domain; it is copyrighted and there are restrictions on its
74 distribution, but these restrictions are designed to permit everything
75 that a good cooperating citizen would want to do. What is not allowed
76 is to try to prevent others from further sharing any version of this
77 library that they might get from you.
79 Specifically, we want to make sure that you have the right to give
80 away copies of the library, that you receive source code or else can get
81 it if you want it, that you can change this library or use pieces of it
82 in new free programs, and that you know you can do these things.
84 To make sure that everyone has such rights, we have to forbid you to
85 deprive anyone else of these rights. For example, if you distribute
86 copies of the GNU MP library, you must give the recipients all the
87 rights that you have. You must make sure that they, too, receive or can
88 get the source code. And you must tell them their rights.
90 Also, for our own protection, we must make certain that everyone
91 finds out that there is no warranty for the GNU MP library. If it is
92 modified by someone else and passed on, we want their recipients to know
93 that what they have is not what we distributed, so that any problems
94 introduced by others will not reflect on our reputation.
96 More precisely, the GNU MP library is dual licensed, under the
97 conditions of the GNU Lesser General Public License version 3 (see
98 'COPYING.LESSERv3'), or the GNU General Public License version 2 (see
99 'COPYINGv2'). This is the recipient's choice, and the recipient also
100 has the additional option of applying later versions of these licenses.
101 (The reason for this dual licensing is to make it possible to use the
102 library with programs which are licensed under GPL version 2, but which
103 for historical or other reasons do not allow use under later versions of
106 Programs which are not part of the library itself, such as
107 demonstration programs and the GMP testsuite, are licensed under the
108 terms of the GNU General Public License version 3 (see 'COPYINGv3'), or
112 File: gmp.info, Node: Introduction to GMP, Next: Installing GMP, Prev: Copying, Up: Top
114 1 Introduction to GNU MP
115 ************************
117 GNU MP is a portable library written in C for arbitrary precision
118 arithmetic on integers, rational numbers, and floating-point numbers.
119 It aims to provide the fastest possible arithmetic for all applications
120 that need higher precision than is directly supported by the basic C
123 Many applications use just a few hundred bits of precision; but some
124 applications may need thousands or even millions of bits. GMP is
125 designed to give good performance for both, by choosing algorithms based
126 on the sizes of the operands, and by carefully keeping the overhead at a
129 The speed of GMP is achieved by using fullwords as the basic
130 arithmetic type, by using sophisticated algorithms, by including
131 carefully optimized assembly code for the most common inner loops for
132 many different CPUs, and by a general emphasis on speed (as opposed to
133 simplicity or elegance).
135 There is assembly code for these CPUs: ARM Cortex-A9, Cortex-A15, and
136 generic ARM, DEC Alpha 21064, 21164, and 21264, AMD K8 and K10 (sold
137 under many brands, e.g. Athlon64, Phenom, Opteron) Bulldozer, and
138 Bobcat, Intel Pentium, Pentium Pro/II/III, Pentium 4, Core2, Nehalem,
139 Sandy bridge, Haswell, generic x86, Intel IA-64, Motorola/IBM PowerPC 32
140 and 64 such as POWER970, POWER5, POWER6, and POWER7, MIPS 32-bit and
141 64-bit, SPARC 32-bit ad 64-bit with special support for all UltraSPARC
142 models. There is also assembly code for many obsolete CPUs.
144 For up-to-date information on GMP, please see the GMP web pages at
146 <https://gmplib.org/>
148 The latest version of the library is available at
150 <https://ftp.gnu.org/gnu/gmp/>
152 Many sites around the world mirror 'ftp.gnu.org', please use a mirror
153 near you, see <https://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 <https://gmplib.org/mailman/listinfo/>.
161 The proper place for bug reports is <gmp-bugs@gmplib.org>. See *note
162 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-6.2.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 separate
247 the dependent parts. Note however that 'gmp.h' is
248 architecture-dependent since it encodes certain aspects of
249 'libgmp', so it will be necessary to ensure both '$prefix/include'
250 and '$exec_prefix/include' are available to the compiler.
252 '--disable-shared', '--disable-static'
253 By default both shared and static libraries are built (where
254 possible), but one or other can be disabled. Shared libraries
255 result in smaller executables and permit code sharing between
256 separate running processes, but on some CPUs are slightly slower,
257 having a small cost on each function call.
259 Native Compilation, '--build=CPU-VENDOR-OS'
260 For normal native compilation, the system can be specified with
261 '--build'. By default './configure' uses the output from running
262 './config.guess'. On some systems './config.guess' can determine
263 the exact CPU type, on others it will be necessary to give it
264 explicitly. For example,
266 ./configure --build=ultrasparc-sun-solaris2.7
268 In all cases the 'OS' part is important, since it controls how
269 libtool generates shared libraries. Running './config.guess' is
270 the simplest way to see what it should be, if you don't know
273 Cross Compilation, '--host=CPU-VENDOR-OS'
274 When cross-compiling, the system used for compiling is given by
275 '--build' and the system where the library will run is given by
276 '--host'. For example when using a FreeBSD Athlon system to build
277 GNU/Linux m68k binaries,
279 ./configure --build=athlon-pc-freebsd3.5 --host=m68k-mac-linux-gnu
281 Compiler tools are sought first with the host system type as a
282 prefix. For example 'm68k-mac-linux-gnu-ranlib' is tried, then
283 plain 'ranlib'. This makes it possible for a set of
284 cross-compiling tools to co-exist with native tools. The prefix is
285 the argument to '--host', and this can be an alias, such as
286 'm68k-linux'. But note that tools don't have to be setup this way,
287 it's enough to just have a 'PATH' with a suitable cross-compiling
290 Compiling for a different CPU in the same family as the build
291 system is a form of cross-compilation, though very possibly this
292 would merely be special options on a native compiler. In any case
293 './configure' avoids depending on being able to run code on the
294 build system, which is important when creating binaries for a newer
295 CPU since they very possibly won't run on the build system.
297 In all cases the compiler must be able to produce an executable (of
298 whatever format) from a standard C 'main'. Although only object
299 files will go to make up 'libgmp', './configure' uses linking tests
300 for various purposes, such as determining what functions are
301 available on the host system.
303 Currently a warning is given unless an explicit '--build' is used
304 when cross-compiling, because it may not be possible to correctly
305 guess the build system type if the 'PATH' has only a
306 cross-compiling 'cc'.
308 Note that the '--target' option is not appropriate for GMP. It's
309 for use when building compiler tools, with '--host' being where
310 they will run, and '--target' what they'll produce code for.
311 Ordinary programs or libraries like GMP are only interested in the
312 '--host' part, being where they'll run. (Some past versions of GMP
313 used '--target' incorrectly.)
316 In general, if you want a library that runs as fast as possible,
317 you should configure GMP for the exact CPU type your system uses.
318 However, this may mean the binaries won't run on older members of
319 the family, and might run slower on other members, older or newer.
320 The best idea is always to build GMP for the exact machine type you
323 The following CPUs have specific support. See 'configure.ac' for
324 details of what code and compiler options they select.
326 * Alpha: alpha, alphaev5, alphaev56, alphapca56, alphapca57,
327 alphaev6, alphaev67, alphaev68 alphaev7
329 * Cray: c90, j90, t90, sv1
331 * HPPA: hppa1.0, hppa1.1, hppa2.0, hppa2.0n, hppa2.0w, hppa64
333 * IA-64: ia64, itanium, itanium2
335 * MIPS: mips, mips3, mips64
337 * Motorola: m68k, m68000, m68010, m68020, m68030, m68040,
338 m68060, m68302, m68360, m88k, m88110
340 * POWER: power, power1, power2, power2sc
342 * PowerPC: powerpc, powerpc64, powerpc401, powerpc403,
343 powerpc405, powerpc505, powerpc601, powerpc602, powerpc603,
344 powerpc603e, powerpc604, powerpc604e, powerpc620, powerpc630,
345 powerpc740, powerpc7400, powerpc7450, powerpc750, powerpc801,
346 powerpc821, powerpc823, powerpc860, powerpc970
348 * SPARC: sparc, sparcv8, microsparc, supersparc, sparcv9,
349 ultrasparc, ultrasparc2, ultrasparc2i, ultrasparc3, sparc64
351 * x86 family: i386, i486, i586, pentium, pentiummmx, pentiumpro,
352 pentium2, pentium3, pentium4, k6, k62, k63, athlon, amd64,
355 * Other: arm, sh, sh2, vax,
357 CPUs not listed will use generic C code.
360 If some of the assembly code causes problems, or if otherwise
361 desired, the generic C code can be selected with the configure
362 '--disable-assembly'.
364 Note that this will run quite slowly, but it should be portable and
365 should at least make it possible to get something running if all
368 Fat binary, '--enable-fat'
369 Using '--enable-fat' selects a "fat binary" build on x86, where
370 optimized low level subroutines are chosen at runtime according to
371 the CPU detected. This means more code, but gives good performance
372 on all x86 chips. (This option might become available for more
373 architectures in the future.)
376 On some systems GMP supports multiple ABIs (application binary
377 interfaces), meaning data type sizes and calling conventions. By
378 default GMP chooses the best ABI available, but a particular ABI
379 can be selected. For example
381 ./configure --host=mips64-sgi-irix6 ABI=n32
383 See *note ABI and ISA::, for the available choices on relevant
384 CPUs, and what applications need to do.
387 By default the C compiler used is chosen from among some likely
388 candidates, with 'gcc' normally preferred if it's present. The
389 usual 'CC=whatever' can be passed to './configure' to choose
392 For various systems, default compiler flags are set based on the
393 CPU and compiler. The usual 'CFLAGS="-whatever"' can be passed to
394 './configure' to use something different or to set good flags for
395 systems GMP doesn't otherwise know.
397 The 'CC' and 'CFLAGS' used are printed during './configure', and
398 can be found in each generated 'Makefile'. This is the easiest way
399 to check the defaults when considering changing or adding
402 Note that when 'CC' and 'CFLAGS' are specified on a system
403 supporting multiple ABIs it's important to give an explicit
404 'ABI=whatever', since GMP can't determine the ABI just from the
405 flags and won't be able to select the correct assembly code.
407 If just 'CC' is selected then normal default 'CFLAGS' for that
408 compiler will be used (if GMP recognises it). For example 'CC=gcc'
409 can be used to force the use of GCC, with default flags (and
413 Any flags like '-D' defines or '-I' includes required by the
414 preprocessor should be set in 'CPPFLAGS' rather than 'CFLAGS'.
415 Compiling is done with both 'CPPFLAGS' and 'CFLAGS', but
416 preprocessing uses just 'CPPFLAGS'. This distinction is because
417 most preprocessors won't accept all the flags the compiler does.
418 Preprocessing is done separately in some configure tests.
421 Some build-time programs are compiled and run to generate
422 host-specific data tables. 'CC_FOR_BUILD' is the compiler used for
423 this. It doesn't need to be in any particular ABI or mode, it
424 merely needs to generate executables that can run. The default is
425 to try the selected 'CC' and some likely candidates such as 'cc'
426 and 'gcc', looking for something that works.
428 No flags are used with 'CC_FOR_BUILD' because a simple invocation
429 like 'cc foo.c' should be enough. If some particular options are
430 required they can be included as for instance 'CC_FOR_BUILD="cc
433 C++ Support, '--enable-cxx'
434 C++ support in GMP can be enabled with '--enable-cxx', in which
435 case a C++ compiler will be required. As a convenience
436 '--enable-cxx=detect' can be used to enable C++ support only if a
437 compiler can be found. The C++ support consists of a library
438 'libgmpxx.la' and header file 'gmpxx.h' (*note Headers and
441 A separate 'libgmpxx.la' has been adopted rather than having C++
442 objects within 'libgmp.la' in order to ensure dynamic linked C
443 programs aren't bloated by a dependency on the C++ standard
444 library, and to avoid any chance that the C++ compiler could be
445 required when linking plain C programs.
447 'libgmpxx.la' will use certain internals from 'libgmp.la' and can
448 only be expected to work with 'libgmp.la' from the same GMP
449 version. Future changes to the relevant internals will be
450 accompanied by renaming, so a mismatch will cause unresolved
451 symbols rather than perhaps mysterious misbehaviour.
453 In general 'libgmpxx.la' will be usable only with the C++ compiler
454 that built it, since name mangling and runtime support are usually
455 incompatible between different compilers.
458 When C++ support is enabled, the C++ compiler and its flags can be
459 set with variables 'CXX' and 'CXXFLAGS' in the usual way. The
460 default for 'CXX' is the first compiler that works from a list of
461 likely candidates, with 'g++' normally preferred when available.
462 The default for 'CXXFLAGS' is to try 'CFLAGS', 'CFLAGS' without
463 '-g', then for 'g++' either '-g -O2' or '-O2', or for other
464 compilers '-g' or nothing. Trying 'CFLAGS' this way is convenient
465 when using 'gcc' and 'g++' together, since the flags for 'gcc' will
468 It's important that the C and C++ compilers match, meaning their
469 startup and runtime support routines are compatible and that they
470 generate code in the same ABI (if there's a choice of ABIs on the
471 system). './configure' isn't currently able to check these things
472 very well itself, so for that reason '--disable-cxx' is the
473 default, to avoid a build failure due to a compiler mismatch.
474 Perhaps this will change in the future.
476 Incidentally, it's normally not good enough to set 'CXX' to the
477 same as 'CC'. Although 'gcc' for instance recognises 'foo.cc' as
478 C++ code, only 'g++' will invoke the linker the right way when
479 building an executable or shared library from C++ object files.
481 Temporary Memory, '--enable-alloca=<choice>'
482 GMP allocates temporary workspace using one of the following three
483 methods, which can be selected with for instance
484 '--enable-alloca=malloc-reentrant'.
486 * 'alloca' - C library or compiler builtin.
487 * 'malloc-reentrant' - the heap, in a re-entrant fashion.
488 * 'malloc-notreentrant' - the heap, with global variables.
490 For convenience, the following choices are also available.
491 '--disable-alloca' is the same as 'no'.
493 * 'yes' - a synonym for 'alloca'.
494 * 'no' - a synonym for 'malloc-reentrant'.
495 * 'reentrant' - 'alloca' if available, otherwise
496 'malloc-reentrant'. This is the default.
497 * 'notreentrant' - 'alloca' if available, otherwise
498 'malloc-notreentrant'.
500 'alloca' is reentrant and fast, and is recommended. It actually
501 allocates just small blocks on the stack; larger ones use
504 'malloc-reentrant' is, as the name suggests, reentrant and thread
505 safe, but 'malloc-notreentrant' is faster and should be used if
506 reentrancy is not required.
508 The two malloc methods in fact use the memory allocation functions
509 selected by 'mp_set_memory_functions', these being 'malloc' and
510 friends by default. *Note Custom Allocation::.
512 An additional choice '--enable-alloca=debug' is available, to help
513 when debugging memory related problems (*note Debugging::).
515 FFT Multiplication, '--disable-fft'
516 By default multiplications are done using Karatsuba, 3-way Toom,
517 higher degree Toom, and Fermat FFT. The FFT is only used on large
518 to very large operands and can be disabled to save code size if
521 Assertion Checking, '--enable-assert'
522 This option enables some consistency checking within the library.
523 This can be of use while debugging, *note Debugging::.
525 Execution Profiling, '--enable-profiling=prof/gprof/instrument'
526 Enable profiling support, in one of various styles, *note
530 Various assembly versions of each mpn subroutines are provided.
531 For a given CPU, a search is made though a path to choose a version
532 of each. For example 'sparcv8' has
534 MPN_PATH="sparc32/v8 sparc32 generic"
536 which means look first for v8 code, then plain sparc32 (which is
537 v7), and finally fall back on generic C. Knowledgeable users with
538 special requirements can specify a different path. Normally this
539 is completely unnecessary.
542 The source for the document you're now reading is 'doc/gmp.texi',
543 in Texinfo format, see *note Texinfo: (texinfo)Top.
545 Info format 'doc/gmp.info' is included in the distribution. The
546 usual automake targets are available to make PostScript, DVI, PDF
547 and HTML (these will require various TeX and Texinfo tools).
549 DocBook and XML can be generated by the Texinfo 'makeinfo' program
550 too, see *note Options for 'makeinfo': (texinfo)makeinfo options.
552 Some supplementary notes can also be found in the 'doc'
556 File: gmp.info, Node: ABI and ISA, Next: Notes for Package Builds, Prev: Build Options, Up: Installing GMP
561 ABI (Application Binary Interface) refers to the calling conventions
562 between functions, meaning what registers are used and what sizes the
563 various C data types are. ISA (Instruction Set Architecture) refers to
564 the instructions and registers a CPU has available.
566 Some 64-bit ISA CPUs have both a 64-bit ABI and a 32-bit ABI defined,
567 the latter for compatibility with older CPUs in the family. GMP
568 supports some CPUs like this in both ABIs. In fact within GMP 'ABI'
569 means a combination of chip ABI, plus how GMP chooses to use it. For
570 example in some 32-bit ABIs, GMP may support a limb as either a 32-bit
571 'long' or a 64-bit 'long long'.
573 By default GMP chooses the best ABI available for a given system, and
574 this generally gives significantly greater speed. But an ABI can be
575 chosen explicitly to make GMP compatible with other libraries, or
576 particular application requirements. For example,
580 In all cases it's vital that all object code used in a given program
581 is compiled for the same ABI.
583 Usually a limb is implemented as a 'long'. When a 'long long' limb
584 is used this is encoded in the generated 'gmp.h'. This is convenient
585 for applications, but it does mean that 'gmp.h' will vary, and can't be
586 just copied around. 'gmp.h' remains compiler independent though, since
587 all compilers for a particular ABI will be expected to use the same limb
590 Currently no attempt is made to follow whatever conventions a system
591 has for installing library or header files built for a particular ABI.
592 This will probably only matter when installing multiple builds of GMP,
593 and it might be as simple as configuring with a special 'libdir', or it
594 might require more than that. Note that builds for different ABIs need
595 to done separately, with a fresh './configure' and 'make' each.
599 On AMD64 systems supporting both 32-bit and 64-bit modes for
600 applications, the following ABI choices are available.
603 The 64-bit ABI uses 64-bit limbs and pointers and makes full
604 use of the chip architecture. This is the default.
605 Applications will usually not need special compiler flags, but
606 for reference the option is
611 The 32-bit ABI is the usual i386 conventions. This will be
612 slower, and is not recommended except for inter-operating with
613 other code not yet 64-bit capable. Applications must be
618 (In GCC 2.95 and earlier there's no '-m32' option, it's the
622 The x32 ABI uses 64-bit limbs but 32-bit pointers. Like the
623 64-bit ABI, it makes full use of the chip's arithmetic
624 capabilities. This ABI is not supported by all operating
630 HPPA 2.0 ('hppa2.0*', 'hppa64')
632 The 2.0w ABI uses 64-bit limbs and pointers and is available
633 on HP-UX 11 or up. Applications must be compiled with
639 The 2.0n ABI means the 32-bit HPPA 1.0 ABI and all its normal
640 calling conventions, but with 64-bit instructions permitted
641 within functions. GMP uses a 64-bit 'long long' for a limb.
642 This ABI is available on hppa64 GNU/Linux and on HP-UX 10 or
643 higher. Applications must be compiled with
648 Note that current versions of GCC (eg. 3.2) don't generate
649 64-bit instructions for 'long long' operations and so may be
650 slower than for 2.0w. (The GMP assembly code is the same
654 HPPA 2.0 CPUs can run all HPPA 1.0 and 1.1 code in the 32-bit
655 HPPA 1.0 ABI. No special compiler options are needed for
658 All three ABIs are available for CPU types 'hppa2.0w', 'hppa2.0'
659 and 'hppa64', but for CPU type 'hppa2.0n' only 2.0n or 1.0 are
662 Note that GCC on HP-UX has no options to choose between 2.0n and
663 2.0w modes, unlike HP 'cc'. Instead it must be built for one or
664 the other ABI. GMP will detect how it was built, and skip to the
668 IA-64 under HP-UX ('ia64*-*-hpux*', 'itanium*-*-hpux*')
669 HP-UX supports two ABIs for IA-64. GMP performance is the same in
673 In the 32-bit ABI, pointers, 'int's and 'long's are 32 bits
674 and GMP uses a 64 bit 'long long' for a limb. Applications
675 can be compiled without any special flags since this ABI is
676 the default in both HP C and GCC, but for reference the flags
683 In the 64-bit ABI, 'long's and pointers are 64 bits and GMP
684 uses a 'long' for a limb. Applications must be compiled with
689 On other IA-64 systems, GNU/Linux for instance, 'ABI=64' is the
693 MIPS under IRIX 6 ('mips*-*-irix[6789]')
694 IRIX 6 always has a 64-bit MIPS 3 or better CPU, and supports ABIs
695 o32, n32, and 64. n32 or 64 are recommended, and GMP performance
696 will be the same in each. The default is n32.
699 The o32 ABI is 32-bit pointers and integers, and no 64-bit
700 operations. GMP will be slower than in n32 or 64, this option
701 only exists to support old compilers, eg. GCC 2.7.2.
702 Applications can be compiled with no special flags on an old
703 compiler, or on a newer compiler with
709 The n32 ABI is 32-bit pointers and integers, but with a 64-bit
710 limb using a 'long long'. Applications must be compiled with
716 The 64-bit ABI is 64-bit pointers and integers. Applications
717 must be compiled with
722 Note that MIPS GNU/Linux, as of kernel version 2.2, doesn't have
723 the necessary support for n32 or 64 and so only gets a 32-bit limb
727 PowerPC 64 ('powerpc64', 'powerpc620', 'powerpc630', 'powerpc970', 'power4', 'power5')
729 The AIX 64 ABI uses 64-bit limbs and pointers and is the
730 default on PowerPC 64 '*-*-aix*' systems. Applications must
736 On 64-bit GNU/Linux, BSD, and Mac OS X/Darwin systems, the
737 applications must be compiled with
742 The 'mode32' ABI uses a 64-bit 'long long' limb but with the
743 chip still in 32-bit mode and using 32-bit calling
744 conventions. This is the default for systems where the true
745 64-bit ABI is unavailable. No special compiler options are
746 typically needed for applications. This ABI is not available
750 This is the basic 32-bit PowerPC ABI, with a 32-bit limb. No
751 special compiler options are needed for applications.
753 GMP's speed is greatest for the 'mode64' ABI, the 'mode32' ABI is
754 2nd best. In 'ABI=32' only the 32-bit ISA is used and this doesn't
755 make full use of a 64-bit chip.
758 Sparc V9 ('sparc64', 'sparcv9', 'ultrasparc*')
760 The 64-bit V9 ABI is available on the various BSD sparc64
761 ports, recent versions of Sparc64 GNU/Linux, and Solaris 2.7
762 and up (when the kernel is in 64-bit mode). GCC 3.2 or
763 higher, or Sun 'cc' is required. On GNU/Linux, depending on
764 the default 'gcc' mode, applications must be compiled with
768 On Solaris applications must be compiled with
770 gcc -m64 -mptr64 -Wa,-xarch=v9 -mcpu=v9
773 On the BSD sparc64 systems no special options are required,
774 since 64-bits is the only ABI available.
777 For the basic 32-bit ABI, GMP still uses as much of the V9 ISA
778 as it can. In the Sun documentation this combination is known
779 as "v8plus". On GNU/Linux, depending on the default 'gcc'
780 mode, applications may need to be compiled with
784 On Solaris, no special compiler options are required for
785 applications, though using something like the following is
786 recommended. ('gcc' 2.8 and earlier only support '-mv8'
792 GMP speed is greatest in 'ABI=64', so it's the default where
793 available. The speed is partly because there are extra registers
794 available and partly because 64-bits is considered the more
795 important case and has therefore had better code written for it.
797 Don't be confused by the names of the '-m' and '-x' compiler
798 options, they're called 'arch' but effectively control both ABI and
801 On Solaris 2.6 and earlier, only 'ABI=32' is available since the
802 kernel doesn't save all registers.
804 On Solaris 2.7 with the kernel in 32-bit mode, a normal native
805 build will reject 'ABI=64' because the resulting executables won't
806 run. 'ABI=64' can still be built if desired by making it look like
807 a cross-compile, for example
809 ./configure --build=none --host=sparcv9-sun-solaris2.7 ABI=64
812 File: gmp.info, Node: Notes for Package Builds, Next: Notes for Particular Systems, Prev: ABI and ISA, Up: Installing GMP
814 2.3 Notes for Package Builds
815 ============================
817 GMP should present no great difficulties for packaging in a binary
820 Libtool is used to build the library and '-version-info' is set
821 appropriately, having started from '3:0:0' in GMP 3.0 (*note Library
822 interface versions: (libtool)Versioning.).
824 The GMP 4 series will be upwardly binary compatible in each release
825 and will be upwardly binary compatible with all of the GMP 3 series.
826 Additional function interfaces may be added in each release, so on
827 systems where libtool versioning is not fully checked by the loader an
828 auxiliary mechanism may be needed to express that a dynamic linked
829 application depends on a new enough GMP.
831 An auxiliary mechanism may also be needed to express that
832 'libgmpxx.la' (from '--enable-cxx', *note Build Options::) requires
833 'libgmp.la' from the same GMP version, since this is not done by the
834 libtool versioning, nor otherwise. A mismatch will result in unresolved
835 symbols from the linker, or perhaps the loader.
837 When building a package for a CPU family, care should be taken to use
838 '--host' (or '--build') to choose the least common denominator among the
839 CPUs which might use the package. For example this might mean plain
840 'sparc' (meaning V7) for SPARCs.
842 For x86s, '--enable-fat' sets things up for a fat binary build,
843 making a runtime selection of optimized low level routines. This is a
844 good choice for packaging to run on a range of x86 chips.
846 Users who care about speed will want GMP built for their exact CPU
847 type, to make best use of the available optimizations. Providing a way
848 to suitably rebuild a package may be useful. This could be as simple as
849 making it possible for a user to omit '--build' (and '--host') so
850 './config.guess' will detect the CPU. But a way to manually specify a
851 '--build' will be wanted for systems where './config.guess' is inexact.
853 On systems with multiple ABIs, a packaged build will need to decide
854 which among the choices is to be provided, see *note ABI and ISA::. A
855 given run of './configure' etc will only build one ABI. If a second ABI
856 is also required then a second run of './configure' etc must be made,
857 starting from a clean directory tree ('make distclean').
859 As noted under "ABI and ISA", currently no attempt is made to follow
860 system conventions for install locations that vary with ABI, such as
861 '/usr/lib/sparcv9' for 'ABI=64' as opposed to '/usr/lib' for 'ABI=32'.
862 A package build can override 'libdir' and other standard variables as
865 Note that 'gmp.h' is a generated file, and will be architecture and
866 ABI dependent. When attempting to install two ABIs simultaneously it
867 will be important that an application compile gets the correct 'gmp.h'
868 for its desired ABI. If compiler include paths don't vary with ABI
869 options then it might be necessary to create a '/usr/include/gmp.h'
870 which tests preprocessor symbols and chooses the correct actual 'gmp.h'.
873 File: gmp.info, Node: Notes for Particular Systems, Next: Known Build Problems, Prev: Notes for Package Builds, Up: Installing GMP
875 2.4 Notes for Particular Systems
876 ================================
879 On systems '*-*-aix[34]*' shared libraries are disabled by default,
880 since some versions of the native 'ar' fail on the convenience
881 libraries used. A shared build can be attempted with
883 ./configure --enable-shared --disable-static
885 Note that the '--disable-static' is necessary because in a shared
886 build libtool makes 'libgmp.a' a symlink to 'libgmp.so', apparently
887 for the benefit of old versions of 'ld' which only recognise '.a',
888 but unfortunately this is done even if a fully functional 'ld' is
892 On systems 'arm*-*-*', versions of GCC up to and including 2.95.3
893 have a bug in unsigned division, giving wrong results for some
894 operands. GMP './configure' will demand GCC 2.95.4 or later.
897 Compaq C++ on OSF 5.1 has two flavours of 'iostream', a standard
898 one and an old pre-standard one (see 'man iostream_intro'). GMP
899 can only use the standard one, which unfortunately is not the
900 default but must be selected by defining '__USE_STD_IOSTREAM'.
901 Configure with for instance
903 ./configure --enable-cxx CPPFLAGS=-D__USE_STD_IOSTREAM
906 On some systems, the hardware floating point has a control mode
907 which can set all operations to be done in a particular precision,
908 for instance single, double or extended on x86 systems (x87
909 floating point). The GMP functions involving a 'double' cannot be
910 expected to operate to their full precision when the hardware is in
911 single precision mode. Of course this affects all code, including
912 application code, not just GMP.
914 FreeBSD 7.x, 8.x, 9.0, 9.1, 9.2
915 'm4' in these releases of FreeBSD has an eval function which
916 ignores its 2nd and 3rd arguments, which makes it unsuitable for
917 '.asm' file processing. './configure' will detect the problem and
918 either abort or choose another m4 in the 'PATH'. The bug is fixed
919 in FreeBSD 9.3 and 10.0, so either upgrade or use GNU m4. Note
920 that the FreeBSD package system installs GNU m4 under the name
921 'gm4', which GMP cannot guess.
923 FreeBSD 7.x, 8.x, 9.x
924 GMP releases starting with 6.0 do not support 'ABI=32' on
925 FreeBSD/amd64 prior to release 10.0 of the system. The cause is a
926 broken 'limits.h', which GMP no longer works around.
928 MS-DOS and MS Windows
929 On an MS-DOS system DJGPP can be used to build GMP, and on an MS
930 Windows system Cygwin, DJGPP and MINGW can be used. All three are
931 excellent ports of GCC and the various GNU tools.
933 <https://www.cygwin.com/>
934 <http://www.delorie.com/djgpp/>
935 <http://www.mingw.org/>
937 Microsoft also publishes an Interix "Services for Unix" which can
938 be used to build GMP on Windows (with a normal './configure'), but
939 it's not free software.
942 On systems '*-*-cygwin*', '*-*-mingw*' and '*-*-pw32*' by default
943 GMP builds only a static library, but a DLL can be built instead
946 ./configure --disable-static --enable-shared
948 Static and DLL libraries can't both be built, since certain export
949 directives in 'gmp.h' must be different.
951 A MINGW DLL build of GMP can be used with Microsoft C. Libtool
952 doesn't install a '.lib' format import library, but it can be
953 created with MS 'lib' as follows, and copied to the install
954 directory. Similarly for 'libmp' and 'libgmpxx'.
957 lib /def:libgmp-3.dll.def /out:libgmp-3.lib
959 MINGW uses the C runtime library 'msvcrt.dll' for I/O, so
960 applications wanting to use the GMP I/O routines must be compiled
961 with 'cl /MD' to do the same. If one of the other C runtime
962 library choices provided by MS C is desired then the suggestion is
963 to use the GMP string functions and confine I/O to the application.
965 Motorola 68k CPU Types
966 'm68k' is taken to mean 68000. 'm68020' or higher will give a
967 performance boost on applicable CPUs. 'm68360' can be used for
968 CPU32 series chips. 'm68302' can be used for "Dragonball" series
969 chips, though this is merely a synonym for 'm68000'.
972 'm4' in these releases of NetBSD has an eval function which ignores
973 its 2nd and 3rd arguments, which makes it unsuitable for '.asm'
974 file processing. './configure' will detect the problem and either
975 abort or choose another m4 in the 'PATH'. The bug is fixed in
976 NetBSD 6, so either upgrade or use GNU m4. Note that the NetBSD
977 package system installs GNU m4 under the name 'gm4', which GMP
981 'm4' in this release of OpenBSD has a bug in 'eval' that makes it
982 unsuitable for '.asm' file processing. './configure' will detect
983 the problem and either abort or choose another m4 in the 'PATH'.
984 The bug is fixed in OpenBSD 2.7, so either upgrade or use GNU m4.
987 In GMP, CPU types 'power*' and 'powerpc*' will each use
988 instructions not available on the other, so it's important to
989 choose the right one for the CPU that will be used. Currently GMP
990 has no assembly code support for using just the common instruction
991 subset. To get executables that run on both, the current
992 suggestion is to use the generic C code ('--disable-assembly'),
993 possibly with appropriate compiler options (like '-mcpu=common' for
994 'gcc'). CPU 'rs6000' (which is not a CPU but a family of
995 workstations) is accepted by 'config.sub', but is currently
996 equivalent to '--disable-assembly'.
999 'sparcv8' or 'supersparc' on relevant systems will give a
1000 significant performance increase over the V7 code selected by plain
1004 The GMP assembly code for both 32-bit and 64-bit Sparc clobbers the
1005 "application registers" 'g2', 'g3' and 'g4', the same way that the
1006 GCC default '-mapp-regs' does (*note SPARC Options: (gcc)SPARC
1009 This makes that code unsuitable for use with the special V9
1010 '-mcmodel=embmedany' (which uses 'g4' as a data segment pointer),
1011 and for applications wanting to use those registers for special
1012 purposes. In these cases the only suggestion currently is to build
1013 GMP with '--disable-assembly' to avoid the assembly code.
1016 '/usr/bin/m4' lacks various features needed to process '.asm'
1017 files, and instead './configure' will automatically use
1018 '/usr/5bin/m4', which we believe is always available (if not then
1022 'i586', 'pentium' or 'pentiummmx' code is good for its intended P5
1023 Pentium chips, but quite slow when run on Intel P6 class chips
1024 (PPro, P-II, P-III). 'i386' is a better choice when making
1025 binaries that must run on both.
1027 x86 MMX and SSE2 Code
1028 If the CPU selected has MMX code but the assembler doesn't support
1029 it, a warning is given and non-MMX code is used instead. This will
1030 be an inferior build, since the MMX code that's present is there
1031 because it's faster than the corresponding plain integer code. The
1032 same applies to SSE2.
1034 Old versions of 'gas' don't support MMX instructions, in particular
1035 version 1.92.3 that comes with FreeBSD 2.2.8 or the more recent
1036 OpenBSD 3.1 doesn't.
1038 Solaris 2.6 and 2.7 'as' generate incorrect object code for
1039 register to register 'movq' instructions, and so can't be used for
1040 MMX code. Install a recent 'gas' if MMX code is wanted on these
1044 File: gmp.info, Node: Known Build Problems, Next: Performance optimization, Prev: Notes for Particular Systems, Up: Installing GMP
1046 2.5 Known Build Problems
1047 ========================
1049 You might find more up-to-date information at <https://gmplib.org/>.
1051 Compiler link options
1052 The version of libtool currently in use rather aggressively strips
1053 compiler options when linking a shared library. This will
1054 hopefully be relaxed in the future, but for now if this is a
1055 problem the suggestion is to create a little script to hide them,
1056 and for instance configure with
1058 ./configure CC=gcc-with-my-options
1060 DJGPP ('*-*-msdosdjgpp*')
1061 The DJGPP port of 'bash' 2.03 is unable to run the 'configure'
1062 script, it exits silently, having died writing a preamble to
1063 'config.log'. Use 'bash' 2.04 or higher.
1065 'make all' was found to run out of memory during the final
1066 'libgmp.la' link on one system tested, despite having 64Mb
1067 available. Running 'make libgmp.la' directly helped, perhaps
1068 recursing into the various subdirectories uses up memory.
1070 GNU binutils 'strip' prior to 2.12
1071 'strip' from GNU binutils 2.11 and earlier should not be used on
1072 the static libraries 'libgmp.a' and 'libmp.a' since it will discard
1073 all but the last of multiple archive members with the same name,
1074 like the three versions of 'init.o' in 'libgmp.a'. Binutils 2.12
1075 or higher can be used successfully.
1077 The shared libraries 'libgmp.so' and 'libmp.so' are not affected by
1078 this and any version of 'strip' can be used on them.
1081 On certain versions of SCO OpenServer 5 and IRIX 6.5 the native
1082 'make' is unable to handle the long dependencies list for
1083 'libgmp.la'. The symptom is a "syntax error" on the following line
1084 of the top-level 'Makefile'.
1086 libgmp.la: $(libgmp_la_OBJECTS) $(libgmp_la_DEPENDENCIES)
1088 Either use GNU Make, or as a workaround remove
1089 '$(libgmp_la_DEPENDENCIES)' from that line (which will make the
1090 initial build work, but if any recompiling is done 'libgmp.la'
1091 might not be rebuilt).
1093 MacOS X ('*-*-darwin*')
1094 Libtool currently only knows how to create shared libraries on
1095 MacOS X using the native 'cc' (which is a modified GCC), not a
1096 plain GCC. A static-only build should work though
1097 ('--disable-shared').
1100 The system compiler on old versions of NeXT was a massacred and old
1101 GCC, even if it called itself 'cc'. This compiler cannot be used
1102 to build GMP, you need to get a real GCC, and install that. (NeXT
1103 may have fixed this in release 3.3 of their system.)
1106 Bugs in GCC 2.7.2 (and 2.6.3) mean it can't be used to compile GMP
1107 on POWER or PowerPC. If you want to use GCC for these machines,
1108 get GCC 2.7.2.1 (or later).
1111 Use the GNU assembler instead of the system assembler, since the
1112 latter has serious bugs.
1115 The system 'sed' prints an error "Output line too long" when
1116 libtool builds 'libgmp.la'. This doesn't seem to cause any obvious
1117 ill effects, but GNU 'sed' is recommended, to avoid any doubt.
1119 Sparc Solaris 2.7 with gcc 2.95.2 in 'ABI=32'
1120 A shared library build of GMP seems to fail in this combination, it
1121 builds but then fails the tests, apparently due to some incorrect
1122 data relocations within 'gmp_randinit_lc_2exp_size'. The exact
1123 cause is unknown, '--disable-shared' is recommended.
1126 File: gmp.info, Node: Performance optimization, Prev: Known Build Problems, Up: Installing GMP
1128 2.6 Performance optimization
1129 ============================
1131 For optimal performance, build GMP for the exact CPU type of the target
1132 computer, see *note Build Options::.
1134 Unlike what is the case for most other programs, the compiler
1135 typically doesn't matter much, since GMP uses assembly language for the
1136 most critical operation.
1138 In particular for long-running GMP applications, and applications
1139 demanding extremely large numbers, building and running the 'tuneup'
1140 program in the 'tune' subdirectory, can be important. For example,
1146 will generate better contents for the 'gmp-mparam.h' parameter file.
1148 To use the results, put the output in the file indicated in the
1149 'Parameters for ...' header. Then recompile from scratch.
1151 The 'tuneup' program takes one useful parameter, '-f NNN', which
1152 instructs the program how long to check FFT multiply parameters. If
1153 you're going to use GMP for extremely large numbers, you may want to run
1154 'tuneup' with a large NNN value.
1157 File: gmp.info, Node: GMP Basics, Next: Reporting Bugs, Prev: Installing GMP, Up: Top
1162 *Using functions, macros, data types, etc. not documented in this manual
1163 is strongly discouraged. If you do so your application is guaranteed to
1164 be incompatible with future versions of GMP.*
1168 * Headers and Libraries::
1169 * Nomenclature and Types::
1170 * Function Classes::
1171 * Variable Conventions::
1172 * Parameter Conventions::
1173 * Memory Management::
1175 * Useful Macros and Constants::
1176 * Compatibility with older versions::
1177 * Demonstration Programs::
1185 File: gmp.info, Node: Headers and Libraries, Next: Nomenclature and Types, Prev: GMP Basics, Up: GMP Basics
1187 3.1 Headers and Libraries
1188 =========================
1190 All declarations needed to use GMP are collected in the include file
1191 'gmp.h'. It is designed to work with both C and C++ compilers.
1195 Note however that prototypes for GMP functions with 'FILE *'
1196 parameters are only provided if '<stdio.h>' is included too.
1201 Likewise '<stdarg.h>' is required for prototypes with 'va_list'
1202 parameters, such as 'gmp_vprintf'. And '<obstack.h>' for prototypes
1203 with 'struct obstack' parameters, such as 'gmp_obstack_printf', when
1206 All programs using GMP must link against the 'libgmp' library. On a
1207 typical Unix-like system this can be done with '-lgmp', for example
1209 gcc myprogram.c -lgmp
1211 GMP C++ functions are in a separate 'libgmpxx' library. This is
1212 built and installed if C++ support has been enabled (*note Build
1213 Options::). For example,
1215 g++ mycxxprog.cc -lgmpxx -lgmp
1217 GMP is built using Libtool and an application can use that to link if
1218 desired, *note GNU Libtool: (libtool)Top.
1220 If GMP has been installed to a non-standard location then it may be
1221 necessary to use '-I' and '-L' compiler options to point to the right
1222 directories, and some sort of run-time path for a shared library.
1225 File: gmp.info, Node: Nomenclature and Types, Next: Function Classes, Prev: Headers and Libraries, Up: GMP Basics
1227 3.2 Nomenclature and Types
1228 ==========================
1230 In this manual, "integer" usually means a multiple precision integer, as
1231 defined by the GMP library. The C data type for such integers is
1232 'mpz_t'. Here are some examples of how to declare such integers:
1236 struct foo { mpz_t x, y; };
1240 "Rational number" means a multiple precision fraction. The C data
1241 type for these fractions is 'mpq_t'. For example:
1245 "Floating point number" or "Float" for short, is an arbitrary
1246 precision mantissa with a limited precision exponent. The C data type
1247 for such objects is 'mpf_t'. For example:
1251 The floating point functions accept and return exponents in the C
1252 type 'mp_exp_t'. Currently this is usually a 'long', but on some
1253 systems it's an 'int' for efficiency.
1255 A "limb" means the part of a multi-precision number that fits in a
1256 single machine word. (We chose this word because a limb of the human
1257 body is analogous to a digit, only larger, and containing several
1258 digits.) Normally a limb is 32 or 64 bits. The C data type for a limb
1261 Counts of limbs of a multi-precision number represented in the C type
1262 'mp_size_t'. Currently this is normally a 'long', but on some systems
1263 it's an 'int' for efficiency, and on some systems it will be 'long long'
1266 Counts of bits of a multi-precision number are represented in the C
1267 type 'mp_bitcnt_t'. Currently this is always an 'unsigned long', but on
1268 some systems it will be an 'unsigned long long' in the future.
1270 "Random state" means an algorithm selection and current state data.
1271 The C data type for such objects is 'gmp_randstate_t'. For example:
1273 gmp_randstate_t rstate;
1275 Also, in general 'mp_bitcnt_t' is used for bit counts and ranges, and
1276 'size_t' is used for byte or character counts.
1279 File: gmp.info, Node: Function Classes, Next: Variable Conventions, Prev: Nomenclature and Types, Up: GMP Basics
1281 3.3 Function Classes
1282 ====================
1284 There are six classes of functions in the GMP library:
1286 1. Functions for signed integer arithmetic, with names beginning with
1287 'mpz_'. The associated type is 'mpz_t'. There are about 150
1288 functions in this class. (*note Integer Functions::)
1290 2. Functions for rational number arithmetic, with names beginning with
1291 'mpq_'. The associated type is 'mpq_t'. There are about 35
1292 functions in this class, but the integer functions can be used for
1293 arithmetic on the numerator and denominator separately. (*note
1294 Rational Number Functions::)
1296 3. Functions for floating-point arithmetic, with names beginning with
1297 'mpf_'. The associated type is 'mpf_t'. There are about 70
1298 functions is this class. (*note Floating-point Functions::)
1300 4. Fast low-level functions that operate on natural numbers. These
1301 are used by the functions in the preceding groups, and you can also
1302 call them directly from very time-critical user programs. These
1303 functions' names begin with 'mpn_'. The associated type is array
1304 of 'mp_limb_t'. There are about 60 (hard-to-use) functions in this
1305 class. (*note Low-level Functions::)
1307 5. Miscellaneous functions. Functions for setting up custom
1308 allocation and functions for generating random numbers. (*note
1309 Custom Allocation::, and *note Random Number Functions::)
1312 File: gmp.info, Node: Variable Conventions, Next: Parameter Conventions, Prev: Function Classes, Up: GMP Basics
1314 3.4 Variable Conventions
1315 ========================
1317 GMP functions generally have output arguments before input arguments.
1318 This notation is by analogy with the assignment operator.
1320 GMP lets you use the same variable for both input and output in one
1321 call. For example, the main function for integer multiplication,
1322 'mpz_mul', can be used to square 'x' and put the result back in 'x' with
1326 Before you can assign to a GMP variable, you need to initialize it by
1327 calling one of the special initialization functions. When you're done
1328 with a variable, you need to clear it out, using one of the functions
1329 for that purpose. Which function to use depends on the type of
1330 variable. See the chapters on integer functions, rational number
1331 functions, and floating-point functions for details.
1333 A variable should only be initialized once, or at least cleared
1334 between each initialization. After a variable has been initialized, it
1335 may be assigned to any number of times.
1337 For efficiency reasons, avoid excessive initializing and clearing.
1338 In general, initialize near the start of a function and clear near the
1347 for (i = 1; i < 100; i++)
1350 mpz_fdiv_q (n, ...);
1356 GMP types like 'mpz_t' are implemented as one-element arrays of
1357 certain structures. Declaring a variable creates an object with the
1358 fields GMP needs, but variables are normally manipulated by using the
1359 pointer to the object. For both behavior and efficiency reasons, it is
1360 discouraged to make copies of the GMP object itself (either directly or
1361 via aggregate objects containing such GMP objects). If copies are done,
1362 all of them must be used read-only; using a copy as the output of some
1363 function will invalidate all the other copies. Note that the actual
1364 fields in each 'mpz_t' etc are for internal use only and should not be
1365 accessed directly by code that expects to be compatible with future GMP
1369 File: gmp.info, Node: Parameter Conventions, Next: Memory Management, Prev: Variable Conventions, Up: GMP Basics
1371 3.5 Parameter Conventions
1372 =========================
1374 When a GMP variable is used as a function parameter, it's effectively a
1375 call-by-reference, meaning that when the function stores a value there
1376 it will change the original in the caller. Parameters which are
1377 input-only can be designated 'const' to provoke a compiler error or
1378 warning on attempting to modify them.
1380 When a function is going to return a GMP result, it should designate
1381 a parameter that it sets, like the library functions do. More than one
1382 value can be returned by having more than one output parameter, again
1383 like the library functions. A 'return' of an 'mpz_t' etc doesn't return
1384 the object, only a pointer, and this is almost certainly not what's
1387 Here's an example accepting an 'mpz_t' parameter, doing a
1388 calculation, and storing the result to the indicated parameter.
1391 foo (mpz_t result, const mpz_t param, unsigned long n)
1394 mpz_mul_ui (result, param, n);
1395 for (i = 1; i < n; i++)
1396 mpz_add_ui (result, result, i*7);
1404 mpz_init_set_str (n, "123456", 0);
1406 gmp_printf ("%Zd\n", r);
1410 Our function 'foo' works even if its caller passes the same variable
1411 for 'param' and 'result', just like the library functions. But
1412 sometimes it's tricky to make that work, and an application might not
1413 want to bother supporting that sort of thing.
1415 Since GMP types are implemented as one-element arrays, using a GMP
1416 variable as a parameter passes a pointer to the object. Hence the
1420 File: gmp.info, Node: Memory Management, Next: Reentrancy, Prev: Parameter Conventions, Up: GMP Basics
1422 3.6 Memory Management
1423 =====================
1425 The GMP types like 'mpz_t' are small, containing only a couple of sizes,
1426 and pointers to allocated data. Once a variable is initialized, GMP
1427 takes care of all space allocation. Additional space is allocated
1428 whenever a variable doesn't have enough.
1430 'mpz_t' and 'mpq_t' variables never reduce their allocated space.
1431 Normally this is the best policy, since it avoids frequent reallocation.
1432 Applications that need to return memory to the heap at some particular
1433 point can use 'mpz_realloc2', or clear variables no longer needed.
1435 'mpf_t' variables, in the current implementation, use a fixed amount
1436 of space, determined by the chosen precision and allocated at
1437 initialization, so their size doesn't change.
1439 All memory is allocated using 'malloc' and friends by default, but
1440 this can be changed, see *note Custom Allocation::. Temporary memory on
1441 the stack is also used (via 'alloca'), but this can be changed at
1442 build-time if desired, see *note Build Options::.
1445 File: gmp.info, Node: Reentrancy, Next: Useful Macros and Constants, Prev: Memory Management, Up: GMP Basics
1450 GMP is reentrant and thread-safe, with some exceptions:
1452 * If configured with '--enable-alloca=malloc-notreentrant' (or with
1453 '--enable-alloca=notreentrant' when 'alloca' is not available),
1454 then naturally GMP is not reentrant.
1456 * 'mpf_set_default_prec' and 'mpf_init' use a global variable for the
1457 selected precision. 'mpf_init2' can be used instead, and in the
1458 C++ interface an explicit precision to the 'mpf_class' constructor.
1460 * 'mpz_random' and the other old random number functions use a global
1461 random state and are hence not reentrant. The newer random number
1462 functions that accept a 'gmp_randstate_t' parameter can be used
1465 * 'gmp_randinit' (obsolete) returns an error indication through a
1466 global variable, which is not thread safe. Applications are
1467 advised to use 'gmp_randinit_default' or 'gmp_randinit_lc_2exp'
1470 * 'mp_set_memory_functions' uses global variables to store the
1471 selected memory allocation functions.
1473 * If the memory allocation functions set by a call to
1474 'mp_set_memory_functions' (or 'malloc' and friends by default) are
1475 not reentrant, then GMP will not be reentrant either.
1477 * If the standard I/O functions such as 'fwrite' are not reentrant
1478 then the GMP I/O functions using them will not be reentrant either.
1480 * It's safe for two threads to read from the same GMP variable
1481 simultaneously, but it's not safe for one to read while another
1482 might be writing, nor for two threads to write simultaneously.
1483 It's not safe for two threads to generate a random number from the
1484 same 'gmp_randstate_t' simultaneously, since this involves an
1485 update of that variable.
1488 File: gmp.info, Node: Useful Macros and Constants, Next: Compatibility with older versions, Prev: Reentrancy, Up: GMP Basics
1490 3.8 Useful Macros and Constants
1491 ===============================
1493 -- Global Constant: const int mp_bits_per_limb
1494 The number of bits per limb.
1496 -- Macro: __GNU_MP_VERSION
1497 -- Macro: __GNU_MP_VERSION_MINOR
1498 -- Macro: __GNU_MP_VERSION_PATCHLEVEL
1499 The major and minor GMP version, and patch level, respectively, as
1500 integers. For GMP i.j, these numbers will be i, j, and 0,
1501 respectively. For GMP i.j.k, these numbers will be i, j, and k,
1504 -- Global Constant: const char * const gmp_version
1505 The GMP version number, as a null-terminated string, in the form
1506 "i.j.k". This release is "6.2.1". Note that the format "i.j" was
1507 used, before version 4.3.0, when k was zero.
1510 -- Macro: __GMP_CFLAGS
1511 The compiler and compiler flags, respectively, used when compiling
1515 File: gmp.info, Node: Compatibility with older versions, Next: Demonstration Programs, Prev: Useful Macros and Constants, Up: GMP Basics
1517 3.9 Compatibility with older versions
1518 =====================================
1520 This version of GMP is upwardly binary compatible with all 5.x, 4.x, and
1521 3.x versions, and upwardly compatible at the source level with all 2.x
1522 versions, with the following exceptions.
1524 * 'mpn_gcd' had its source arguments swapped as of GMP 3.0, for
1525 consistency with other 'mpn' functions.
1527 * 'mpf_get_prec' counted precision slightly differently in GMP 3.0
1528 and 3.0.1, but in 3.1 reverted to the 2.x style.
1530 * 'mpn_bdivmod', documented as preliminary in GMP 4, has been
1533 There are a number of compatibility issues between GMP 1 and GMP 2
1534 that of course also apply when porting applications from GMP 1 to GMP 5.
1535 Please see the GMP 2 manual for details.
1538 File: gmp.info, Node: Demonstration Programs, Next: Efficiency, Prev: Compatibility with older versions, Up: GMP Basics
1540 3.10 Demonstration programs
1541 ===========================
1543 The 'demos' subdirectory has some sample programs using GMP. These
1544 aren't built or installed, but there's a 'Makefile' with rules for them.
1550 The following programs are provided
1552 * 'pexpr' is an expression evaluator, the program used on the GMP web
1554 * The 'calc' subdirectory has a similar but simpler evaluator using
1556 * The 'expr' subdirectory is yet another expression evaluator, a
1557 library designed for ease of use within a C program. See
1558 'demos/expr/README' for more information.
1559 * 'factorize' is a Pollard-Rho factorization program.
1560 * 'isprime' is a command-line interface to the 'mpz_probab_prime_p'
1562 * 'primes' counts or lists primes in an interval, using a sieve.
1563 * 'qcn' is an example use of 'mpz_kronecker_ui' to estimate quadratic
1565 * The 'perl' subdirectory is a comprehensive perl interface to GMP.
1566 See 'demos/perl/INSTALL' for more information. Documentation is in
1567 POD format in 'demos/perl/GMP.pm'.
1569 As an aside, consideration has been given at various times to some
1570 sort of expression evaluation within the main GMP library. Going beyond
1571 something minimal quickly leads to matters like user-defined functions,
1572 looping, fixnums for control variables, etc, which are considered
1573 outside the scope of GMP (much closer to language interpreters or
1574 compilers, *Note Language Bindings::.) Something simple for program
1575 input convenience may yet be a possibility, a combination of the 'expr'
1576 demo and the 'pexpr' tree back-end perhaps. But for now the above
1577 evaluators are offered as illustrations.
1580 File: gmp.info, Node: Efficiency, Next: Debugging, Prev: Demonstration Programs, Up: GMP Basics
1586 On small operands, the time for function call overheads and memory
1587 allocation can be significant in comparison to actual calculation.
1588 This is unavoidable in a general purpose variable precision
1589 library, although GMP attempts to be as efficient as it can on both
1590 large and small operands.
1593 On some CPUs, in particular the x86s, the static 'libgmp.a' should
1594 be used for maximum speed, since the PIC code in the shared
1595 'libgmp.so' will have a small overhead on each function call and
1596 global data address. For many programs this will be insignificant,
1597 but for long calculations there's a gain to be had.
1599 Initializing and Clearing
1600 Avoid excessive initializing and clearing of variables, since this
1601 can be quite time consuming, especially in comparison to otherwise
1602 fast operations like addition.
1604 A language interpreter might want to keep a free list or stack of
1605 initialized variables ready for use. It should be possible to
1606 integrate something like that with a garbage collector too.
1609 An 'mpz_t' or 'mpq_t' variable used to hold successively increasing
1610 values will have its memory repeatedly 'realloc'ed, which could be
1611 quite slow or could fragment memory, depending on the C library.
1612 If an application can estimate the final size then 'mpz_init2' or
1613 'mpz_realloc2' can be called to allocate the necessary space from
1614 the beginning (*note Initializing Integers::).
1616 It doesn't matter if a size set with 'mpz_init2' or 'mpz_realloc2'
1617 is too small, since all functions will do a further reallocation if
1618 necessary. Badly overestimating memory required will waste space
1622 It's up to an application to call functions like 'mpz_mul_2exp'
1623 when appropriate. General purpose functions like 'mpz_mul' make no
1624 attempt to identify powers of two or other special forms, because
1625 such inputs will usually be very rare and testing every time would
1628 'ui' and 'si' Functions
1629 The 'ui' functions and the small number of 'si' functions exist for
1630 convenience and should be used where applicable. But if for
1631 example an 'mpz_t' contains a value that fits in an 'unsigned long'
1632 there's no need extract it and call a 'ui' function, just use the
1633 regular 'mpz' function.
1636 'mpz_abs', 'mpq_abs', 'mpf_abs', 'mpz_neg', 'mpq_neg' and 'mpf_neg'
1637 are fast when used for in-place operations like 'mpz_abs(x,x)',
1638 since in the current implementation only a single field of 'x'
1639 needs changing. On suitable compilers (GCC for instance) this is
1642 'mpz_add_ui', 'mpz_sub_ui', 'mpf_add_ui' and 'mpf_sub_ui' benefit
1643 from an in-place operation like 'mpz_add_ui(x,x,y)', since usually
1644 only one or two limbs of 'x' will need to be changed. The same
1645 applies to the full precision 'mpz_add' etc if 'y' is small. If
1646 'y' is big then cache locality may be helped, but that's all.
1648 'mpz_mul' is currently the opposite, a separate destination is
1649 slightly better. A call like 'mpz_mul(x,x,y)' will, unless 'y' is
1650 only one limb, make a temporary copy of 'x' before forming the
1651 result. Normally that copying will only be a tiny fraction of the
1652 time for the multiply, so this is not a particularly important
1655 'mpz_set', 'mpq_set', 'mpq_set_num', 'mpf_set', etc, make no
1656 attempt to recognise a copy of something to itself, so a call like
1657 'mpz_set(x,x)' will be wasteful. Naturally that would never be
1658 written deliberately, but if it might arise from two pointers to
1659 the same object then a test to avoid it might be desirable.
1664 Note that it's never worth introducing extra 'mpz_set' calls just
1665 to get in-place operations. If a result should go to a particular
1666 variable then just direct it there and let GMP take care of data
1669 Divisibility Testing (Small Integers)
1670 'mpz_divisible_ui_p' and 'mpz_congruent_ui_p' are the best
1671 functions for testing whether an 'mpz_t' is divisible by an
1672 individual small integer. They use an algorithm which is faster
1673 than 'mpz_tdiv_ui', but which gives no useful information about the
1674 actual remainder, only whether it's zero (or a particular value).
1676 However when testing divisibility by several small integers, it's
1677 best to take a remainder modulo their product, to save
1678 multi-precision operations. For instance to test whether a number
1679 is divisible by any of 23, 29 or 31 take a remainder modulo
1680 23*29*31 = 20677 and then test that.
1682 The division functions like 'mpz_tdiv_q_ui' which give a quotient
1683 as well as a remainder are generally a little slower than the
1684 remainder-only functions like 'mpz_tdiv_ui'. If the quotient is
1685 only rarely wanted then it's probably best to just take a remainder
1686 and then go back and calculate the quotient if and when it's wanted
1687 ('mpz_divexact_ui' can be used if the remainder is zero).
1690 The 'mpq' functions operate on 'mpq_t' values with no common
1691 factors in the numerator and denominator. Common factors are
1692 checked-for and cast out as necessary. In general, cancelling
1693 factors every time is the best approach since it minimizes the
1694 sizes for subsequent operations.
1696 However, applications that know something about the factorization
1697 of the values they're working with might be able to avoid some of
1698 the GCDs used for canonicalization, or swap them for divisions.
1699 For example when multiplying by a prime it's enough to check for
1700 factors of it in the denominator instead of doing a full GCD. Or
1701 when forming a big product it might be known that very little
1702 cancellation will be possible, and so canonicalization can be left
1705 The 'mpq_numref' and 'mpq_denref' macros give access to the
1706 numerator and denominator to do things outside the scope of the
1707 supplied 'mpq' functions. *Note Applying Integer Functions::.
1709 The canonical form for rationals allows mixed-type 'mpq_t' and
1710 integer additions or subtractions to be done directly with
1711 multiples of the denominator. This will be somewhat faster than
1712 'mpq_add'. For example,
1715 mpz_add (mpq_numref(q), mpq_numref(q), mpq_denref(q));
1717 /* mpq += unsigned long */
1718 mpz_addmul_ui (mpq_numref(q), mpq_denref(q), 123UL);
1721 mpz_submul (mpq_numref(q), mpq_denref(q), z);
1724 Functions like 'mpz_fac_ui', 'mpz_fib_ui' and 'mpz_bin_uiui' are
1725 designed for calculating isolated values. If a range of values is
1726 wanted it's probably best to get a starting point and iterate from
1730 Hexadecimal or octal are suggested for input or output in text
1731 form. Power-of-2 bases like these can be converted much more
1732 efficiently than other bases, like decimal. For big numbers
1733 there's usually nothing of particular interest to be seen in the
1734 digits, so the base doesn't matter much.
1736 Maybe we can hope octal will one day become the normal base for
1737 everyday use, as proposed by King Charles XII of Sweden and later
1741 File: gmp.info, Node: Debugging, Next: Profiling, Prev: Efficiency, Up: GMP Basics
1747 Depending on the system, a segmentation violation or bus error
1748 might be the only indication of stack overflow. See
1749 '--enable-alloca' choices in *note Build Options::, for how to
1752 In new enough versions of GCC, '-fstack-check' may be able to
1753 ensure an overflow is recognised by the system before too much
1754 damage is done, or '-fstack-limit-symbol' or
1755 '-fstack-limit-register' may be able to add checking if the system
1756 itself doesn't do any (*note Options for Code Generation: (gcc)Code
1757 Gen Options.). These options must be added to the 'CFLAGS' used in
1758 the GMP build (*note Build Options::), adding them just to an
1759 application will have no effect. Note also they're a slowdown,
1760 adding overhead to each function call and each stack allocation.
1763 The most likely cause of application problems with GMP is heap
1764 corruption. Failing to 'init' GMP variables will have
1765 unpredictable effects, and corruption arising elsewhere in a
1766 program may well affect GMP. Initializing GMP variables more than
1767 once or failing to clear them will cause memory leaks.
1769 In all such cases a 'malloc' debugger is recommended. On a GNU or
1770 BSD system the standard C library 'malloc' has some diagnostic
1771 facilities, see *note Allocation Debugging: (libc)Allocation
1772 Debugging, or 'man 3 malloc'. Other possibilities, in no
1773 particular order, include
1775 <http://cs.ecs.baylor.edu/~donahoo/tools/ccmalloc/>
1776 <http://dmalloc.com/>
1777 <https://wiki.gnome.org/Apps/MemProf>
1779 The GMP default allocation routines in 'memory.c' also have a
1780 simple sentinel scheme which can be enabled with '#define DEBUG' in
1781 that file. This is mainly designed for detecting buffer overruns
1782 during GMP development, but might find other uses.
1785 On some systems the compiler options GMP uses by default can
1786 interfere with debugging. In particular on x86 and 68k systems
1787 '-fomit-frame-pointer' is used and this generally inhibits stack
1788 backtracing. Recompiling without such options may help while
1789 debugging, though the usual caveats about it potentially moving a
1790 memory problem or hiding a compiler bug will apply.
1792 GDB, the GNU Debugger
1793 A sample '.gdbinit' is included in the distribution, showing how to
1794 call some undocumented dump functions to print GMP variables from
1795 within GDB. Note that these functions shouldn't be used in final
1796 application code since they're undocumented and may be subject to
1797 incompatible changes in future versions of GMP.
1800 GMP has multiple source files with the same name, in different
1801 directories. For example 'mpz', 'mpq' and 'mpf' each have an
1802 'init.c'. If the debugger can't already determine the right one it
1803 may help to build with absolute paths on each C file. One way to
1804 do that is to use a separate object directory with an absolute path
1805 to the source directory.
1808 /my/source/dir/gmp-6.2.1/configure
1810 This works via 'VPATH', and might require GNU 'make'. Alternately
1811 it might be possible to change the '.c.lo' rules appropriately.
1814 The build option '--enable-assert' is available to add some
1815 consistency checks to the library (see *note Build Options::).
1816 These are likely to be of limited value to most applications.
1817 Assertion failures are just as likely to indicate memory corruption
1818 as a library or compiler bug.
1820 Applications using the low-level 'mpn' functions, however, will
1821 benefit from '--enable-assert' since it adds checks on the
1822 parameters of most such functions, many of which have subtle
1823 restrictions on their usage. Note however that only the generic C
1824 code has checks, not the assembly code, so '--disable-assembly'
1825 should be used for maximum checking.
1827 Temporary Memory Checking
1828 The build option '--enable-alloca=debug' arranges that each block
1829 of temporary memory in GMP is allocated with a separate call to
1830 'malloc' (or the allocation function set with
1831 'mp_set_memory_functions').
1833 This can help a malloc debugger detect accesses outside the
1834 intended bounds, or detect memory not released. In a normal build,
1835 on the other hand, temporary memory is allocated in blocks which
1836 GMP divides up for its own use, or may be allocated with a compiler
1837 builtin 'alloca' which will go nowhere near any malloc debugger
1840 Maximum Debuggability
1841 To summarize the above, a GMP build for maximum debuggability would
1844 ./configure --disable-shared --enable-assert \
1845 --enable-alloca=debug --disable-assembly CFLAGS=-g
1847 For C++, add '--enable-cxx CXXFLAGS=-g'.
1850 The GCC checker (<https://savannah.nongnu.org/projects/checker/>)
1851 can be used with GMP. It contains a stub library which means GMP
1852 applications compiled with checker can use a normal GMP build.
1854 A build of GMP with checking within GMP itself can be made. This
1855 will run very very slowly. On GNU/Linux for example,
1857 ./configure --disable-assembly CC=checkergcc
1859 '--disable-assembly' must be used, since the GMP assembly code
1860 doesn't support the checking scheme. The GMP C++ features cannot
1861 be used, since current versions of checker (0.9.9.1) don't yet
1862 support the standard C++ library.
1865 Valgrind (<http://valgrind.org/>) is a memory checker for x86, ARM,
1866 MIPS, PowerPC, and S/390. It translates and emulates machine
1867 instructions to do strong checks for uninitialized data (at the
1868 level of individual bits), memory accesses through bad pointers,
1871 Valgrind does not always support every possible instruction, in
1872 particular ones recently added to an ISA. Valgrind might therefore
1873 be incompatible with a recent GMP or even a less recent GMP which
1874 is compiled using a recent GCC.
1876 GMP's assembly code sometimes promotes a read of the limbs to some
1877 larger size, for efficiency. GMP will do this even at the start
1878 and end of a multilimb operand, using naturally aligned operations
1879 on the larger type. This may lead to benign reads outside of
1880 allocated areas, triggering complaints from Valgrind. Valgrind's
1881 option '--partial-loads-ok=yes' should help.
1884 Any suspected bug in GMP itself should be isolated to make sure
1885 it's not an application problem, see *note Reporting Bugs::.
1888 File: gmp.info, Node: Profiling, Next: Autoconf, Prev: Debugging, Up: GMP Basics
1893 Running a program under a profiler is a good way to find where it's
1894 spending most time and where improvements can be best sought. The
1895 profiling choices for a GMP build are as follows.
1897 '--disable-profiling'
1898 The default is to add nothing special for profiling.
1900 It should be possible to just compile the mainline of a program
1901 with '-p' and use 'prof' to get a profile consisting of timer-based
1902 sampling of the program counter. Most of the GMP assembly code has
1903 the necessary symbol information.
1905 This approach has the advantage of minimizing interference with
1906 normal program operation, but on most systems the resolution of the
1907 sampling is quite low (10 milliseconds for instance), requiring
1908 long runs to get accurate information.
1910 '--enable-profiling=prof'
1911 Build with support for the system 'prof', which means '-p' added to
1914 This provides call counting in addition to program counter
1915 sampling, which allows the most frequently called routines to be
1916 identified, and an average time spent in each routine to be
1919 The x86 assembly code has support for this option, but on other
1920 processors the assembly routines will be as if compiled without
1921 '-p' and therefore won't appear in the call counts.
1923 On some systems, such as GNU/Linux, '-p' in fact means '-pg' and in
1924 this case '--enable-profiling=gprof' described below should be used
1927 '--enable-profiling=gprof'
1928 Build with support for 'gprof', which means '-pg' added to the
1931 This provides call graph construction in addition to call counting
1932 and program counter sampling, which makes it possible to count
1933 calls coming from different locations. For example the number of
1934 calls to 'mpn_mul' from 'mpz_mul' versus the number from 'mpf_mul'.
1935 The program counter sampling is still flat though, so only a total
1936 time in 'mpn_mul' would be accumulated, not a separate amount for
1939 The x86 assembly code has support for this option, but on other
1940 processors the assembly routines will be as if compiled without
1941 '-pg' and therefore not be included in the call counts.
1943 On x86 and m68k systems '-pg' and '-fomit-frame-pointer' are
1944 incompatible, so the latter is omitted from the default flags in
1945 that case, which might result in poorer code generation.
1947 Incidentally, it should be possible to use the 'gprof' program with
1948 a plain '--enable-profiling=prof' build. But in that case only the
1949 'gprof -p' flat profile and call counts can be expected to be
1950 valid, not the 'gprof -q' call graph.
1952 '--enable-profiling=instrument'
1953 Build with the GCC option '-finstrument-functions' added to the
1954 'CFLAGS' (*note Options for Code Generation: (gcc)Code Gen
1957 This inserts special instrumenting calls at the start and end of
1958 each function, allowing exact timing and full call graph
1961 This instrumenting is not normally a standard system feature and
1962 will require support from an external library, such as
1964 <https://sourceforge.net/projects/fnccheck/>
1966 This should be included in 'LIBS' during the GMP configure so that
1967 test programs will link. For example,
1969 ./configure --enable-profiling=instrument LIBS=-lfc
1971 On a GNU system the C library provides dummy instrumenting
1972 functions, so programs compiled with this option will link. In
1973 this case it's only necessary to ensure the correct library is
1974 added when linking an application.
1976 The x86 assembly code supports this option, but on other processors
1977 the assembly routines will be as if compiled without
1978 '-finstrument-functions' meaning time spent in them will
1979 effectively be attributed to their caller.
1982 File: gmp.info, Node: Autoconf, Next: Emacs, Prev: Profiling, Up: GMP Basics
1987 Autoconf based applications can easily check whether GMP is installed.
1988 The only thing to be noted is that GMP library symbols from version 3
1989 onwards have prefixes like '__gmpz'. The following therefore would be a
1992 AC_CHECK_LIB(gmp, __gmpz_init)
1994 This just uses the default 'AC_CHECK_LIB' actions for found or not
1995 found, but an application that must have GMP would want to generate an
1996 error if not found. For example,
1998 AC_CHECK_LIB(gmp, __gmpz_init, ,
1999 [AC_MSG_ERROR([GNU MP not found, see https://gmplib.org/])])
2001 If functions added in some particular version of GMP are required,
2002 then one of those can be used when checking. For example 'mpz_mul_si'
2003 was added in GMP 3.1,
2005 AC_CHECK_LIB(gmp, __gmpz_mul_si, ,
2007 [GNU MP not found, or not 3.1 or up, see https://gmplib.org/])])
2009 An alternative would be to test the version number in 'gmp.h' using
2010 say 'AC_EGREP_CPP'. That would make it possible to test the exact
2011 version, if some particular sub-minor release is known to be necessary.
2013 In general it's recommended that applications should simply demand a
2014 new enough GMP rather than trying to provide supplements for features
2015 not available in past versions.
2017 Occasionally an application will need or want to know the size of a
2018 type at configuration or preprocessing time, not just with 'sizeof' in
2019 the code. This can be done in the normal way with 'mp_limb_t' etc, but
2020 GMP 4.0 or up is best for this, since prior versions needed certain '-D'
2021 defines on systems using a 'long long' limb. The following would suit
2022 Autoconf 2.50 or up,
2024 AC_CHECK_SIZEOF(mp_limb_t, , [#include <gmp.h>])
2027 File: gmp.info, Node: Emacs, Prev: Autoconf, Up: GMP Basics
2032 <C-h C-i> ('info-lookup-symbol') is a good way to find documentation on
2033 C functions while editing (*note Info Documentation Lookup: (emacs)Info
2036 The GMP manual can be included in such lookups by putting the
2037 following in your '.emacs',
2039 (eval-after-load "info-look"
2040 '(let ((mode-value (assoc 'c-mode (assoc 'symbol info-lookup-alist))))
2041 (setcar (nthcdr 3 mode-value)
2042 (cons '("(gmp)Function Index" nil "^ -.* " "\\>")
2043 (nth 3 mode-value)))))
2046 File: gmp.info, Node: Reporting Bugs, Next: Integer Functions, Prev: GMP Basics, Up: Top
2051 If you think you have found a bug in the GMP library, please investigate
2052 it and report it. We have made this library available to you, and it is
2053 not too much to ask you to report the bugs you find.
2055 Before you report a bug, check it's not already addressed in *note
2056 Known Build Problems::, or perhaps *note Notes for Particular Systems::.
2057 You may also want to check <https://gmplib.org/> for patches for this
2060 Please include the following in any report,
2062 * The GMP version number, and if pre-packaged or patched then say so.
2064 * A test program that makes it possible for us to reproduce the bug.
2065 Include instructions on how to run the program.
2067 * A description of what is wrong. If the results are incorrect, in
2068 what way. If you get a crash, say so.
2070 * If you get a crash, include a stack backtrace from the debugger if
2071 it's informative ('where' in 'gdb', or '$C' in 'adb').
2073 * Please do not send core dumps, executables or 'strace's.
2075 * The 'configure' options you used when building GMP, if any.
2077 * The output from 'configure', as printed to stdout, with any options
2080 * The name of the compiler and its version. For 'gcc', get the
2081 version with 'gcc -v', otherwise perhaps 'what `which cc`', or
2084 * The output from running 'uname -a'.
2086 * The output from running './config.guess', and from running
2087 './configfsf.guess' (might be the same).
2089 * If the bug is related to 'configure', then the compressed contents
2092 * If the bug is related to an 'asm' file not assembling, then the
2093 contents of 'config.m4' and the offending line or lines from the
2094 temporary 'mpn/tmp-<file>.s'.
2096 Please make an effort to produce a self-contained report, with
2097 something definite that can be tested or debugged. Vague queries or
2098 piecemeal messages are difficult to act on and don't help the
2101 It is not uncommon that an observed problem is actually due to a bug
2102 in the compiler; the GMP code tends to explore interesting corners in
2105 If your bug report is good, we will do our best to help you get a
2106 corrected version of the library; if the bug report is poor, we won't do
2107 anything about it (except maybe ask you to send a better report).
2109 Send your report to: <gmp-bugs@gmplib.org>.
2111 If you think something in this manual is unclear, or downright
2112 incorrect, or if the language needs to be improved, please send a note
2113 to the same address.
2116 File: gmp.info, Node: Integer Functions, Next: Rational Number Functions, Prev: Reporting Bugs, Up: Top
2121 This chapter describes the GMP functions for performing integer
2122 arithmetic. These functions start with the prefix 'mpz_'.
2124 GMP integers are stored in objects of type 'mpz_t'.
2128 * Initializing Integers::
2129 * Assigning Integers::
2130 * Simultaneous Integer Init & Assign::
2131 * Converting Integers::
2132 * Integer Arithmetic::
2133 * Integer Division::
2134 * Integer Exponentiation::
2136 * Number Theoretic Functions::
2137 * Integer Comparisons::
2138 * Integer Logic and Bit Fiddling::
2140 * Integer Random Numbers::
2141 * Integer Import and Export::
2142 * Miscellaneous Integer Functions::
2143 * Integer Special Functions::
2146 File: gmp.info, Node: Initializing Integers, Next: Assigning Integers, Prev: Integer Functions, Up: Integer Functions
2148 5.1 Initialization Functions
2149 ============================
2151 The functions for integer arithmetic assume that all integer objects are
2152 initialized. You do that by calling the function 'mpz_init'. For
2159 mpz_add (integ, ...);
2161 mpz_sub (integ, ...);
2163 /* Unless the program is about to exit, do ... */
2167 As you can see, you can store new values any number of times, once an
2168 object is initialized.
2170 -- Function: void mpz_init (mpz_t X)
2171 Initialize X, and set its value to 0.
2173 -- Function: void mpz_inits (mpz_t X, ...)
2174 Initialize a NULL-terminated list of 'mpz_t' variables, and set
2177 -- Function: void mpz_init2 (mpz_t X, mp_bitcnt_t N)
2178 Initialize X, with space for N-bit numbers, and set its value to 0.
2179 Calling this function instead of 'mpz_init' or 'mpz_inits' is never
2180 necessary; reallocation is handled automatically by GMP when
2183 While N defines the initial space, X will grow automatically in the
2184 normal way, if necessary, for subsequent values stored.
2185 'mpz_init2' makes it possible to avoid such reallocations if a
2186 maximum size is known in advance.
2188 In preparation for an operation, GMP often allocates one limb more
2189 than ultimately needed. To make sure GMP will not perform
2190 reallocation for X, you need to add the number of bits in
2193 -- Function: void mpz_clear (mpz_t X)
2194 Free the space occupied by X. Call this function for all 'mpz_t'
2195 variables when you are done with them.
2197 -- Function: void mpz_clears (mpz_t X, ...)
2198 Free the space occupied by a NULL-terminated list of 'mpz_t'
2201 -- Function: void mpz_realloc2 (mpz_t X, mp_bitcnt_t N)
2202 Change the space allocated for X to N bits. The value in X is
2203 preserved if it fits, or is set to 0 if not.
2205 Calling this function is never necessary; reallocation is handled
2206 automatically by GMP when needed. But this function can be used to
2207 increase the space for a variable in order to avoid repeated
2208 automatic reallocations, or to decrease it to give memory back to
2212 File: gmp.info, Node: Assigning Integers, Next: Simultaneous Integer Init & Assign, Prev: Initializing Integers, Up: Integer Functions
2214 5.2 Assignment Functions
2215 ========================
2217 These functions assign new values to already initialized integers (*note
2218 Initializing Integers::).
2220 -- Function: void mpz_set (mpz_t ROP, const mpz_t OP)
2221 -- Function: void mpz_set_ui (mpz_t ROP, unsigned long int OP)
2222 -- Function: void mpz_set_si (mpz_t ROP, signed long int OP)
2223 -- Function: void mpz_set_d (mpz_t ROP, double OP)
2224 -- Function: void mpz_set_q (mpz_t ROP, const mpq_t OP)
2225 -- Function: void mpz_set_f (mpz_t ROP, const mpf_t OP)
2226 Set the value of ROP from OP.
2228 'mpz_set_d', 'mpz_set_q' and 'mpz_set_f' truncate OP to make it an
2231 -- Function: int mpz_set_str (mpz_t ROP, const char *STR, int BASE)
2232 Set the value of ROP from STR, a null-terminated C string in base
2233 BASE. White space is allowed in the string, and is simply ignored.
2235 The BASE may vary from 2 to 62, or if BASE is 0, then the leading
2236 characters are used: '0x' and '0X' for hexadecimal, '0b' and '0B'
2237 for binary, '0' for octal, or decimal otherwise.
2239 For bases up to 36, case is ignored; upper-case and lower-case
2240 letters have the same value. For bases 37 to 62, upper-case letter
2241 represent the usual 10..35 while lower-case letter represent
2244 This function returns 0 if the entire string is a valid number in
2245 base BASE. Otherwise it returns -1.
2247 -- Function: void mpz_swap (mpz_t ROP1, mpz_t ROP2)
2248 Swap the values ROP1 and ROP2 efficiently.
2251 File: gmp.info, Node: Simultaneous Integer Init & Assign, Next: Converting Integers, Prev: Assigning Integers, Up: Integer Functions
2253 5.3 Combined Initialization and Assignment Functions
2254 ====================================================
2256 For convenience, GMP provides a parallel series of initialize-and-set
2257 functions which initialize the output and then store the value there.
2258 These functions' names have the form 'mpz_init_set...'
2260 Here is an example of using one:
2264 mpz_init_set_str (pie, "3141592653589793238462643383279502884", 10);
2271 Once the integer has been initialized by any of the 'mpz_init_set...'
2272 functions, it can be used as the source or destination operand for the
2273 ordinary integer functions. Don't use an initialize-and-set function on
2274 a variable already initialized!
2276 -- Function: void mpz_init_set (mpz_t ROP, const mpz_t OP)
2277 -- Function: void mpz_init_set_ui (mpz_t ROP, unsigned long int OP)
2278 -- Function: void mpz_init_set_si (mpz_t ROP, signed long int OP)
2279 -- Function: void mpz_init_set_d (mpz_t ROP, double OP)
2280 Initialize ROP with limb space and set the initial numeric value
2283 -- Function: int mpz_init_set_str (mpz_t ROP, const char *STR, int
2285 Initialize ROP and set its value like 'mpz_set_str' (see its
2286 documentation above for details).
2288 If the string is a correct base BASE number, the function returns
2289 0; if an error occurs it returns -1. ROP is initialized even if an
2290 error occurs. (I.e., you have to call 'mpz_clear' for it.)
2293 File: gmp.info, Node: Converting Integers, Next: Integer Arithmetic, Prev: Simultaneous Integer Init & Assign, Up: Integer Functions
2295 5.4 Conversion Functions
2296 ========================
2298 This section describes functions for converting GMP integers to standard
2299 C types. Functions for converting _to_ GMP integers are described in
2300 *note Assigning Integers:: and *note I/O of Integers::.
2302 -- Function: unsigned long int mpz_get_ui (const mpz_t OP)
2303 Return the value of OP as an 'unsigned long'.
2305 If OP is too big to fit an 'unsigned long' then just the least
2306 significant bits that do fit are returned. The sign of OP is
2307 ignored, only the absolute value is used.
2309 -- Function: signed long int mpz_get_si (const mpz_t OP)
2310 If OP fits into a 'signed long int' return the value of OP.
2311 Otherwise return the least significant part of OP, with the same
2314 If OP is too big to fit in a 'signed long int', the returned result
2315 is probably not very useful. To find out if the value will fit,
2316 use the function 'mpz_fits_slong_p'.
2318 -- Function: double mpz_get_d (const mpz_t OP)
2319 Convert OP to a 'double', truncating if necessary (i.e. rounding
2322 If the exponent from the conversion is too big, the result is
2323 system dependent. An infinity is returned where available. A
2324 hardware overflow trap may or may not occur.
2326 -- Function: double mpz_get_d_2exp (signed long int *EXP, const mpz_t
2328 Convert OP to a 'double', truncating if necessary (i.e. rounding
2329 towards zero), and returning the exponent separately.
2331 The return value is in the range 0.5<=abs(D)<1 and the exponent is
2332 stored to '*EXP'. D * 2^EXP is the (truncated) OP value. If OP is
2333 zero, the return is 0.0 and 0 is stored to '*EXP'.
2335 This is similar to the standard C 'frexp' function (*note
2336 (libc)Normalization Functions::).
2338 -- Function: char * mpz_get_str (char *STR, int BASE, const mpz_t OP)
2339 Convert OP to a string of digits in base BASE. The base argument
2340 may vary from 2 to 62 or from -2 to -36.
2342 For BASE in the range 2..36, digits and lower-case letters are
2343 used; for -2..-36, digits and upper-case letters are used; for
2344 37..62, digits, upper-case letters, and lower-case letters (in that
2345 significance order) are used.
2347 If STR is 'NULL', the result string is allocated using the current
2348 allocation function (*note Custom Allocation::). The block will be
2349 'strlen(str)+1' bytes, that being exactly enough for the string and
2352 If STR is not 'NULL', it should point to a block of storage large
2353 enough for the result, that being 'mpz_sizeinbase (OP, BASE) + 2'.
2354 The two extra bytes are for a possible minus sign, and the
2357 A pointer to the result string is returned, being either the
2358 allocated block, or the given STR.
2361 File: gmp.info, Node: Integer Arithmetic, Next: Integer Division, Prev: Converting Integers, Up: Integer Functions
2363 5.5 Arithmetic Functions
2364 ========================
2366 -- Function: void mpz_add (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2367 -- Function: void mpz_add_ui (mpz_t ROP, const mpz_t OP1, unsigned long
2369 Set ROP to OP1 + OP2.
2371 -- Function: void mpz_sub (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2372 -- Function: void mpz_sub_ui (mpz_t ROP, const mpz_t OP1, unsigned long
2374 -- Function: void mpz_ui_sub (mpz_t ROP, unsigned long int OP1, const
2376 Set ROP to OP1 - OP2.
2378 -- Function: void mpz_mul (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2379 -- Function: void mpz_mul_si (mpz_t ROP, const mpz_t OP1, long int OP2)
2380 -- Function: void mpz_mul_ui (mpz_t ROP, const mpz_t OP1, unsigned long
2382 Set ROP to OP1 times OP2.
2384 -- Function: void mpz_addmul (mpz_t ROP, const mpz_t OP1, const mpz_t
2386 -- Function: void mpz_addmul_ui (mpz_t ROP, const mpz_t OP1, unsigned
2388 Set ROP to ROP + OP1 times OP2.
2390 -- Function: void mpz_submul (mpz_t ROP, const mpz_t OP1, const mpz_t
2392 -- Function: void mpz_submul_ui (mpz_t ROP, const mpz_t OP1, unsigned
2394 Set ROP to ROP - OP1 times OP2.
2396 -- Function: void mpz_mul_2exp (mpz_t ROP, const mpz_t OP1, mp_bitcnt_t
2398 Set ROP to OP1 times 2 raised to OP2. This operation can also be
2399 defined as a left shift by OP2 bits.
2401 -- Function: void mpz_neg (mpz_t ROP, const mpz_t OP)
2404 -- Function: void mpz_abs (mpz_t ROP, const mpz_t OP)
2405 Set ROP to the absolute value of OP.
2408 File: gmp.info, Node: Integer Division, Next: Integer Exponentiation, Prev: Integer Arithmetic, Up: Integer Functions
2410 5.6 Division Functions
2411 ======================
2413 Division is undefined if the divisor is zero. Passing a zero divisor to
2414 the division or modulo functions (including the modular powering
2415 functions 'mpz_powm' and 'mpz_powm_ui'), will cause an intentional
2416 division by zero. This lets a program handle arithmetic exceptions in
2417 these functions the same way as for normal C 'int' arithmetic.
2419 -- Function: void mpz_cdiv_q (mpz_t Q, const mpz_t N, const mpz_t D)
2420 -- Function: void mpz_cdiv_r (mpz_t R, const mpz_t N, const mpz_t D)
2421 -- Function: void mpz_cdiv_qr (mpz_t Q, mpz_t R, const mpz_t N, const
2424 -- Function: unsigned long int mpz_cdiv_q_ui (mpz_t Q, const mpz_t N,
2425 unsigned long int D)
2426 -- Function: unsigned long int mpz_cdiv_r_ui (mpz_t R, const mpz_t N,
2427 unsigned long int D)
2428 -- Function: unsigned long int mpz_cdiv_qr_ui (mpz_t Q, mpz_t R,
2429 const mpz_t N, unsigned long int D)
2430 -- Function: unsigned long int mpz_cdiv_ui (const mpz_t N,
2431 unsigned long int D)
2433 -- Function: void mpz_cdiv_q_2exp (mpz_t Q, const mpz_t N,
2435 -- Function: void mpz_cdiv_r_2exp (mpz_t R, const mpz_t N,
2438 -- Function: void mpz_fdiv_q (mpz_t Q, const mpz_t N, const mpz_t D)
2439 -- Function: void mpz_fdiv_r (mpz_t R, const mpz_t N, const mpz_t D)
2440 -- Function: void mpz_fdiv_qr (mpz_t Q, mpz_t R, const mpz_t N, const
2443 -- Function: unsigned long int mpz_fdiv_q_ui (mpz_t Q, const mpz_t N,
2444 unsigned long int D)
2445 -- Function: unsigned long int mpz_fdiv_r_ui (mpz_t R, const mpz_t N,
2446 unsigned long int D)
2447 -- Function: unsigned long int mpz_fdiv_qr_ui (mpz_t Q, mpz_t R,
2448 const mpz_t N, unsigned long int D)
2449 -- Function: unsigned long int mpz_fdiv_ui (const mpz_t N,
2450 unsigned long int D)
2452 -- Function: void mpz_fdiv_q_2exp (mpz_t Q, const mpz_t N,
2454 -- Function: void mpz_fdiv_r_2exp (mpz_t R, const mpz_t N,
2457 -- Function: void mpz_tdiv_q (mpz_t Q, const mpz_t N, const mpz_t D)
2458 -- Function: void mpz_tdiv_r (mpz_t R, const mpz_t N, const mpz_t D)
2459 -- Function: void mpz_tdiv_qr (mpz_t Q, mpz_t R, const mpz_t N, const
2462 -- Function: unsigned long int mpz_tdiv_q_ui (mpz_t Q, const mpz_t N,
2463 unsigned long int D)
2464 -- Function: unsigned long int mpz_tdiv_r_ui (mpz_t R, const mpz_t N,
2465 unsigned long int D)
2466 -- Function: unsigned long int mpz_tdiv_qr_ui (mpz_t Q, mpz_t R,
2467 const mpz_t N, unsigned long int D)
2468 -- Function: unsigned long int mpz_tdiv_ui (const mpz_t N,
2469 unsigned long int D)
2471 -- Function: void mpz_tdiv_q_2exp (mpz_t Q, const mpz_t N,
2473 -- Function: void mpz_tdiv_r_2exp (mpz_t R, const mpz_t N,
2477 Divide N by D, forming a quotient Q and/or remainder R. For the
2478 '2exp' functions, D=2^B. The rounding is in three styles, each
2479 suiting different applications.
2481 * 'cdiv' rounds Q up towards +infinity, and R will have the
2482 opposite sign to D. The 'c' stands for "ceil".
2484 * 'fdiv' rounds Q down towards -infinity, and R will have the
2485 same sign as D. The 'f' stands for "floor".
2487 * 'tdiv' rounds Q towards zero, and R will have the same sign as
2488 N. The 't' stands for "truncate".
2490 In all cases Q and R will satisfy N=Q*D+R, and R will satisfy
2493 The 'q' functions calculate only the quotient, the 'r' functions
2494 only the remainder, and the 'qr' functions calculate both. Note
2495 that for 'qr' the same variable cannot be passed for both Q and R,
2496 or results will be unpredictable.
2498 For the 'ui' variants the return value is the remainder, and in
2499 fact returning the remainder is all the 'div_ui' functions do. For
2500 'tdiv' and 'cdiv' the remainder can be negative, so for those the
2501 return value is the absolute value of the remainder.
2503 For the '2exp' variants the divisor is 2^B. These functions are
2504 implemented as right shifts and bit masks, but of course they round
2505 the same as the other functions.
2507 For positive N both 'mpz_fdiv_q_2exp' and 'mpz_tdiv_q_2exp' are
2508 simple bitwise right shifts. For negative N, 'mpz_fdiv_q_2exp' is
2509 effectively an arithmetic right shift treating N as twos complement
2510 the same as the bitwise logical functions do, whereas
2511 'mpz_tdiv_q_2exp' effectively treats N as sign and magnitude.
2513 -- Function: void mpz_mod (mpz_t R, const mpz_t N, const mpz_t D)
2514 -- Function: unsigned long int mpz_mod_ui (mpz_t R, const mpz_t N,
2515 unsigned long int D)
2516 Set R to N 'mod' D. The sign of the divisor is ignored; the result
2517 is always non-negative.
2519 'mpz_mod_ui' is identical to 'mpz_fdiv_r_ui' above, returning the
2520 remainder as well as setting R. See 'mpz_fdiv_ui' above if only
2521 the return value is wanted.
2523 -- Function: void mpz_divexact (mpz_t Q, const mpz_t N, const mpz_t D)
2524 -- Function: void mpz_divexact_ui (mpz_t Q, const mpz_t N, unsigned
2526 Set Q to N/D. These functions produce correct results only when it
2527 is known in advance that D divides N.
2529 These routines are much faster than the other division functions,
2530 and are the best choice when exact division is known to occur, for
2531 example reducing a rational to lowest terms.
2533 -- Function: int mpz_divisible_p (const mpz_t N, const mpz_t D)
2534 -- Function: int mpz_divisible_ui_p (const mpz_t N, unsigned long int
2536 -- Function: int mpz_divisible_2exp_p (const mpz_t N, mp_bitcnt_t B)
2537 Return non-zero if N is exactly divisible by D, or in the case of
2538 'mpz_divisible_2exp_p' by 2^B.
2540 N is divisible by D if there exists an integer Q satisfying N =
2541 Q*D. Unlike the other division functions, D=0 is accepted and
2542 following the rule it can be seen that only 0 is considered
2545 -- Function: int mpz_congruent_p (const mpz_t N, const mpz_t C, const
2547 -- Function: int mpz_congruent_ui_p (const mpz_t N, unsigned long int
2548 C, unsigned long int D)
2549 -- Function: int mpz_congruent_2exp_p (const mpz_t N, const mpz_t C,
2551 Return non-zero if N is congruent to C modulo D, or in the case of
2552 'mpz_congruent_2exp_p' modulo 2^B.
2554 N is congruent to C mod D if there exists an integer Q satisfying N
2555 = C + Q*D. Unlike the other division functions, D=0 is accepted
2556 and following the rule it can be seen that N and C are considered
2557 congruent mod 0 only when exactly equal.
2560 File: gmp.info, Node: Integer Exponentiation, Next: Integer Roots, Prev: Integer Division, Up: Integer Functions
2562 5.7 Exponentiation Functions
2563 ============================
2565 -- Function: void mpz_powm (mpz_t ROP, const mpz_t BASE, const mpz_t
2566 EXP, const mpz_t MOD)
2567 -- Function: void mpz_powm_ui (mpz_t ROP, const mpz_t BASE, unsigned
2568 long int EXP, const mpz_t MOD)
2569 Set ROP to (BASE raised to EXP) modulo MOD.
2571 Negative EXP is supported if the inverse BASE^(-1) mod MOD exists
2572 (see 'mpz_invert' in *note Number Theoretic Functions::). If an
2573 inverse doesn't exist then a divide by zero is raised.
2575 -- Function: void mpz_powm_sec (mpz_t ROP, const mpz_t BASE, const
2576 mpz_t EXP, const mpz_t MOD)
2577 Set ROP to (BASE raised to EXP) modulo MOD.
2579 It is required that EXP > 0 and that MOD is odd.
2581 This function is designed to take the same time and have the same
2582 cache access patterns for any two same-size arguments, assuming
2583 that function arguments are placed at the same position and that
2584 the machine state is identical upon function entry. This function
2585 is intended for cryptographic purposes, where resilience to
2586 side-channel attacks is desired.
2588 -- Function: void mpz_pow_ui (mpz_t ROP, const mpz_t BASE, unsigned
2590 -- Function: void mpz_ui_pow_ui (mpz_t ROP, unsigned long int BASE,
2591 unsigned long int EXP)
2592 Set ROP to BASE raised to EXP. The case 0^0 yields 1.
2595 File: gmp.info, Node: Integer Roots, Next: Number Theoretic Functions, Prev: Integer Exponentiation, Up: Integer Functions
2597 5.8 Root Extraction Functions
2598 =============================
2600 -- Function: int mpz_root (mpz_t ROP, const mpz_t OP, unsigned long int
2602 Set ROP to the truncated integer part of the Nth root of OP.
2603 Return non-zero if the computation was exact, i.e., if OP is ROP to
2606 -- Function: void mpz_rootrem (mpz_t ROOT, mpz_t REM, const mpz_t U,
2607 unsigned long int N)
2608 Set ROOT to the truncated integer part of the Nth root of U. Set
2609 REM to the remainder, U-ROOT**N.
2611 -- Function: void mpz_sqrt (mpz_t ROP, const mpz_t OP)
2612 Set ROP to the truncated integer part of the square root of OP.
2614 -- Function: void mpz_sqrtrem (mpz_t ROP1, mpz_t ROP2, const mpz_t OP)
2615 Set ROP1 to the truncated integer part of the square root of OP,
2616 like 'mpz_sqrt'. Set ROP2 to the remainder OP-ROP1*ROP1, which
2617 will be zero if OP is a perfect square.
2619 If ROP1 and ROP2 are the same variable, the results are undefined.
2621 -- Function: int mpz_perfect_power_p (const mpz_t OP)
2622 Return non-zero if OP is a perfect power, i.e., if there exist
2623 integers A and B, with B>1, such that OP equals A raised to the
2626 Under this definition both 0 and 1 are considered to be perfect
2627 powers. Negative values of OP are accepted, but of course can only
2628 be odd perfect powers.
2630 -- Function: int mpz_perfect_square_p (const mpz_t OP)
2631 Return non-zero if OP is a perfect square, i.e., if the square root
2632 of OP is an integer. Under this definition both 0 and 1 are
2633 considered to be perfect squares.
2636 File: gmp.info, Node: Number Theoretic Functions, Next: Integer Comparisons, Prev: Integer Roots, Up: Integer Functions
2638 5.9 Number Theoretic Functions
2639 ==============================
2641 -- Function: int mpz_probab_prime_p (const mpz_t N, int REPS)
2642 Determine whether N is prime. Return 2 if N is definitely prime,
2643 return 1 if N is probably prime (without being certain), or return
2644 0 if N is definitely non-prime.
2646 This function performs some trial divisions, a Baillie-PSW probable
2647 prime test, then REPS-24 Miller-Rabin probabilistic primality
2648 tests. A higher REPS value will reduce the chances of a non-prime
2649 being identified as "probably prime". A composite number will be
2650 identified as a prime with an asymptotic probability of less than
2651 4^(-REPS). Reasonable values of REPS are between 15 and 50.
2653 GMP versions up to and including 6.1.2 did not use the Baillie-PSW
2654 primality test. In those older versions of GMP, this function
2655 performed REPS Miller-Rabin tests.
2657 -- Function: void mpz_nextprime (mpz_t ROP, const mpz_t OP)
2658 Set ROP to the next prime greater than OP.
2660 This function uses a probabilistic algorithm to identify primes.
2661 For practical purposes it's adequate, the chance of a composite
2662 passing will be extremely small.
2664 -- Function: void mpz_gcd (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2665 Set ROP to the greatest common divisor of OP1 and OP2. The result
2666 is always positive even if one or both input operands are negative.
2667 Except if both inputs are zero; then this function defines gcd(0,0)
2670 -- Function: unsigned long int mpz_gcd_ui (mpz_t ROP, const mpz_t OP1,
2671 unsigned long int OP2)
2672 Compute the greatest common divisor of OP1 and OP2. If ROP is not
2673 'NULL', store the result there.
2675 If the result is small enough to fit in an 'unsigned long int', it
2676 is returned. If the result does not fit, 0 is returned, and the
2677 result is equal to the argument OP1. Note that the result will
2678 always fit if OP2 is non-zero.
2680 -- Function: void mpz_gcdext (mpz_t G, mpz_t S, mpz_t T, const mpz_t A,
2682 Set G to the greatest common divisor of A and B, and in addition
2683 set S and T to coefficients satisfying A*S + B*T = G. The value in
2684 G is always positive, even if one or both of A and B are negative
2685 (or zero if both inputs are zero). The values in S and T are
2686 chosen such that normally, abs(S) < abs(B) / (2 G) and abs(T) <
2687 abs(A) / (2 G), and these relations define S and T uniquely. There
2688 are a few exceptional cases:
2690 If abs(A) = abs(B), then S = 0, T = sgn(B).
2692 Otherwise, S = sgn(A) if B = 0 or abs(B) = 2 G, and T = sgn(B) if A
2693 = 0 or abs(A) = 2 G.
2695 In all cases, S = 0 if and only if G = abs(B), i.e., if B divides A
2698 If T or G is 'NULL' then that value is not computed.
2700 -- Function: void mpz_lcm (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2701 -- Function: void mpz_lcm_ui (mpz_t ROP, const mpz_t OP1, unsigned long
2703 Set ROP to the least common multiple of OP1 and OP2. ROP is always
2704 positive, irrespective of the signs of OP1 and OP2. ROP will be
2705 zero if either OP1 or OP2 is zero.
2707 -- Function: int mpz_invert (mpz_t ROP, const mpz_t OP1, const mpz_t
2709 Compute the inverse of OP1 modulo OP2 and put the result in ROP.
2710 If the inverse exists, the return value is non-zero and ROP will
2711 satisfy 0 <= ROP < abs(OP2) (with ROP = 0 possible only when
2712 abs(OP2) = 1, i.e., in the somewhat degenerate zero ring). If an
2713 inverse doesn't exist the return value is zero and ROP is
2714 undefined. The behaviour of this function is undefined when OP2 is
2717 -- Function: int mpz_jacobi (const mpz_t A, const mpz_t B)
2718 Calculate the Jacobi symbol (A/B). This is defined only for B odd.
2720 -- Function: int mpz_legendre (const mpz_t A, const mpz_t P)
2721 Calculate the Legendre symbol (A/P). This is defined only for P an
2722 odd positive prime, and for such P it's identical to the Jacobi
2725 -- Function: int mpz_kronecker (const mpz_t A, const mpz_t B)
2726 -- Function: int mpz_kronecker_si (const mpz_t A, long B)
2727 -- Function: int mpz_kronecker_ui (const mpz_t A, unsigned long B)
2728 -- Function: int mpz_si_kronecker (long A, const mpz_t B)
2729 -- Function: int mpz_ui_kronecker (unsigned long A, const mpz_t B)
2730 Calculate the Jacobi symbol (A/B) with the Kronecker extension
2731 (a/2)=(2/a) when a odd, or (a/2)=0 when a even.
2733 When B is odd the Jacobi symbol and Kronecker symbol are identical,
2734 so 'mpz_kronecker_ui' etc can be used for mixed precision Jacobi
2737 For more information see Henri Cohen section 1.4.2 (*note
2738 References::), or any number theory textbook. See also the example
2739 program 'demos/qcn.c' which uses 'mpz_kronecker_ui'.
2741 -- Function: mp_bitcnt_t mpz_remove (mpz_t ROP, const mpz_t OP, const
2743 Remove all occurrences of the factor F from OP and store the result
2744 in ROP. The return value is how many such occurrences were
2747 -- Function: void mpz_fac_ui (mpz_t ROP, unsigned long int N)
2748 -- Function: void mpz_2fac_ui (mpz_t ROP, unsigned long int N)
2749 -- Function: void mpz_mfac_uiui (mpz_t ROP, unsigned long int N,
2750 unsigned long int M)
2751 Set ROP to the factorial of N: 'mpz_fac_ui' computes the plain
2752 factorial N!, 'mpz_2fac_ui' computes the double-factorial N!!, and
2753 'mpz_mfac_uiui' the M-multi-factorial N!^(M).
2755 -- Function: void mpz_primorial_ui (mpz_t ROP, unsigned long int N)
2756 Set ROP to the primorial of N, i.e. the product of all positive
2759 -- Function: void mpz_bin_ui (mpz_t ROP, const mpz_t N, unsigned long
2761 -- Function: void mpz_bin_uiui (mpz_t ROP, unsigned long int N,
2762 unsigned long int K)
2763 Compute the binomial coefficient N over K and store the result in
2764 ROP. Negative values of N are supported by 'mpz_bin_ui', using the
2765 identity bin(-n,k) = (-1)^k * bin(n+k-1,k), see Knuth volume 1
2766 section 1.2.6 part G.
2768 -- Function: void mpz_fib_ui (mpz_t FN, unsigned long int N)
2769 -- Function: void mpz_fib2_ui (mpz_t FN, mpz_t FNSUB1, unsigned long
2771 'mpz_fib_ui' sets FN to to F[n], the N'th Fibonacci number.
2772 'mpz_fib2_ui' sets FN to F[n], and FNSUB1 to F[n-1].
2774 These functions are designed for calculating isolated Fibonacci
2775 numbers. When a sequence of values is wanted it's best to start
2776 with 'mpz_fib2_ui' and iterate the defining F[n+1]=F[n]+F[n-1] or
2779 -- Function: void mpz_lucnum_ui (mpz_t LN, unsigned long int N)
2780 -- Function: void mpz_lucnum2_ui (mpz_t LN, mpz_t LNSUB1, unsigned long
2782 'mpz_lucnum_ui' sets LN to to L[n], the N'th Lucas number.
2783 'mpz_lucnum2_ui' sets LN to L[n], and LNSUB1 to L[n-1].
2785 These functions are designed for calculating isolated Lucas
2786 numbers. When a sequence of values is wanted it's best to start
2787 with 'mpz_lucnum2_ui' and iterate the defining L[n+1]=L[n]+L[n-1]
2790 The Fibonacci numbers and Lucas numbers are related sequences, so
2791 it's never necessary to call both 'mpz_fib2_ui' and
2792 'mpz_lucnum2_ui'. The formulas for going from Fibonacci to Lucas
2793 can be found in *note Lucas Numbers Algorithm::, the reverse is
2794 straightforward too.
2797 File: gmp.info, Node: Integer Comparisons, Next: Integer Logic and Bit Fiddling, Prev: Number Theoretic Functions, Up: Integer Functions
2799 5.10 Comparison Functions
2800 =========================
2802 -- Function: int mpz_cmp (const mpz_t OP1, const mpz_t OP2)
2803 -- Function: int mpz_cmp_d (const mpz_t OP1, double OP2)
2804 -- Macro: int mpz_cmp_si (const mpz_t OP1, signed long int OP2)
2805 -- Macro: int mpz_cmp_ui (const mpz_t OP1, unsigned long int OP2)
2806 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero if
2807 OP1 = OP2, or a negative value if OP1 < OP2.
2809 'mpz_cmp_ui' and 'mpz_cmp_si' are macros and will evaluate their
2810 arguments more than once. 'mpz_cmp_d' can be called with an
2811 infinity, but results are undefined for a NaN.
2813 -- Function: int mpz_cmpabs (const mpz_t OP1, const mpz_t OP2)
2814 -- Function: int mpz_cmpabs_d (const mpz_t OP1, double OP2)
2815 -- Function: int mpz_cmpabs_ui (const mpz_t OP1, unsigned long int OP2)
2816 Compare the absolute values of OP1 and OP2. Return a positive
2817 value if abs(OP1) > abs(OP2), zero if abs(OP1) = abs(OP2), or a
2818 negative value if abs(OP1) < abs(OP2).
2820 'mpz_cmpabs_d' can be called with an infinity, but results are
2821 undefined for a NaN.
2823 -- Macro: int mpz_sgn (const mpz_t OP)
2824 Return +1 if OP > 0, 0 if OP = 0, and -1 if OP < 0.
2826 This function is actually implemented as a macro. It evaluates its
2827 argument multiple times.
2830 File: gmp.info, Node: Integer Logic and Bit Fiddling, Next: I/O of Integers, Prev: Integer Comparisons, Up: Integer Functions
2832 5.11 Logical and Bit Manipulation Functions
2833 ===========================================
2835 These functions behave as if twos complement arithmetic were used
2836 (although sign-magnitude is the actual implementation). The least
2837 significant bit is number 0.
2839 -- Function: void mpz_and (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2840 Set ROP to OP1 bitwise-and OP2.
2842 -- Function: void mpz_ior (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2843 Set ROP to OP1 bitwise inclusive-or OP2.
2845 -- Function: void mpz_xor (mpz_t ROP, const mpz_t OP1, const mpz_t OP2)
2846 Set ROP to OP1 bitwise exclusive-or OP2.
2848 -- Function: void mpz_com (mpz_t ROP, const mpz_t OP)
2849 Set ROP to the one's complement of OP.
2851 -- Function: mp_bitcnt_t mpz_popcount (const mpz_t OP)
2852 If OP>=0, return the population count of OP, which is the number of
2853 1 bits in the binary representation. If OP<0, the number of 1s is
2854 infinite, and the return value is the largest possible
2857 -- Function: mp_bitcnt_t mpz_hamdist (const mpz_t OP1, const mpz_t OP2)
2858 If OP1 and OP2 are both >=0 or both <0, return the hamming distance
2859 between the two operands, which is the number of bit positions
2860 where OP1 and OP2 have different bit values. If one operand is >=0
2861 and the other <0 then the number of bits different is infinite, and
2862 the return value is the largest possible 'mp_bitcnt_t'.
2864 -- Function: mp_bitcnt_t mpz_scan0 (const mpz_t OP, mp_bitcnt_t
2866 -- Function: mp_bitcnt_t mpz_scan1 (const mpz_t OP, mp_bitcnt_t
2868 Scan OP, starting from bit STARTING_BIT, towards more significant
2869 bits, until the first 0 or 1 bit (respectively) is found. Return
2870 the index of the found bit.
2872 If the bit at STARTING_BIT is already what's sought, then
2873 STARTING_BIT is returned.
2875 If there's no bit found, then the largest possible 'mp_bitcnt_t' is
2876 returned. This will happen in 'mpz_scan0' past the end of a
2877 negative number, or 'mpz_scan1' past the end of a nonnegative
2880 -- Function: void mpz_setbit (mpz_t ROP, mp_bitcnt_t BIT_INDEX)
2881 Set bit BIT_INDEX in ROP.
2883 -- Function: void mpz_clrbit (mpz_t ROP, mp_bitcnt_t BIT_INDEX)
2884 Clear bit BIT_INDEX in ROP.
2886 -- Function: void mpz_combit (mpz_t ROP, mp_bitcnt_t BIT_INDEX)
2887 Complement bit BIT_INDEX in ROP.
2889 -- Function: int mpz_tstbit (const mpz_t OP, mp_bitcnt_t BIT_INDEX)
2890 Test bit BIT_INDEX in OP and return 0 or 1 accordingly.
2893 File: gmp.info, Node: I/O of Integers, Next: Integer Random Numbers, Prev: Integer Logic and Bit Fiddling, Up: Integer Functions
2895 5.12 Input and Output Functions
2896 ===============================
2898 Functions that perform input from a stdio stream, and functions that
2899 output to a stdio stream, of 'mpz' numbers. Passing a 'NULL' pointer
2900 for a STREAM argument to any of these functions will make them read from
2901 'stdin' and write to 'stdout', respectively.
2903 When using any of these functions, it is a good idea to include
2904 'stdio.h' before 'gmp.h', since that will allow 'gmp.h' to define
2905 prototypes for these functions.
2907 See also *note Formatted Output:: and *note Formatted Input::.
2909 -- Function: size_t mpz_out_str (FILE *STREAM, int BASE, const mpz_t
2911 Output OP on stdio stream STREAM, as a string of digits in base
2912 BASE. The base argument may vary from 2 to 62 or from -2 to -36.
2914 For BASE in the range 2..36, digits and lower-case letters are
2915 used; for -2..-36, digits and upper-case letters are used; for
2916 37..62, digits, upper-case letters, and lower-case letters (in that
2917 significance order) are used.
2919 Return the number of bytes written, or if an error occurred, return
2922 -- Function: size_t mpz_inp_str (mpz_t ROP, FILE *STREAM, int BASE)
2923 Input a possibly white-space preceded string in base BASE from
2924 stdio stream STREAM, and put the read integer in ROP.
2926 The BASE may vary from 2 to 62, or if BASE is 0, then the leading
2927 characters are used: '0x' and '0X' for hexadecimal, '0b' and '0B'
2928 for binary, '0' for octal, or decimal otherwise.
2930 For bases up to 36, case is ignored; upper-case and lower-case
2931 letters have the same value. For bases 37 to 62, upper-case letter
2932 represent the usual 10..35 while lower-case letter represent
2935 Return the number of bytes read, or if an error occurred, return 0.
2937 -- Function: size_t mpz_out_raw (FILE *STREAM, const mpz_t OP)
2938 Output OP on stdio stream STREAM, in raw binary format. The
2939 integer is written in a portable format, with 4 bytes of size
2940 information, and that many bytes of limbs. Both the size and the
2941 limbs are written in decreasing significance order (i.e., in
2944 The output can be read with 'mpz_inp_raw'.
2946 Return the number of bytes written, or if an error occurred, return
2949 The output of this can not be read by 'mpz_inp_raw' from GMP 1,
2950 because of changes necessary for compatibility between 32-bit and
2953 -- Function: size_t mpz_inp_raw (mpz_t ROP, FILE *STREAM)
2954 Input from stdio stream STREAM in the format written by
2955 'mpz_out_raw', and put the result in ROP. Return the number of
2956 bytes read, or if an error occurred, return 0.
2958 This routine can read the output from 'mpz_out_raw' also from GMP
2959 1, in spite of changes necessary for compatibility between 32-bit
2960 and 64-bit machines.
2963 File: gmp.info, Node: Integer Random Numbers, Next: Integer Import and Export, Prev: I/O of Integers, Up: Integer Functions
2965 5.13 Random Number Functions
2966 ============================
2968 The random number functions of GMP come in two groups; older function
2969 that rely on a global state, and newer functions that accept a state
2970 parameter that is read and modified. Please see the *note Random Number
2971 Functions:: for more information on how to use and not to use random
2974 -- Function: void mpz_urandomb (mpz_t ROP, gmp_randstate_t STATE,
2976 Generate a uniformly distributed random integer in the range 0 to
2979 The variable STATE must be initialized by calling one of the
2980 'gmp_randinit' functions (*note Random State Initialization::)
2981 before invoking this function.
2983 -- Function: void mpz_urandomm (mpz_t ROP, gmp_randstate_t STATE, const
2985 Generate a uniform random integer in the range 0 to N-1, inclusive.
2987 The variable STATE must be initialized by calling one of the
2988 'gmp_randinit' functions (*note Random State Initialization::)
2989 before invoking this function.
2991 -- Function: void mpz_rrandomb (mpz_t ROP, gmp_randstate_t STATE,
2993 Generate a random integer with long strings of zeros and ones in
2994 the binary representation. Useful for testing functions and
2995 algorithms, since this kind of random numbers have proven to be
2996 more likely to trigger corner-case bugs. The random number will be
2997 in the range 2^(N-1) to 2^N-1, inclusive.
2999 The variable STATE must be initialized by calling one of the
3000 'gmp_randinit' functions (*note Random State Initialization::)
3001 before invoking this function.
3003 -- Function: void mpz_random (mpz_t ROP, mp_size_t MAX_SIZE)
3004 Generate a random integer of at most MAX_SIZE limbs. The generated
3005 random number doesn't satisfy any particular requirements of
3006 randomness. Negative random numbers are generated when MAX_SIZE is
3009 This function is obsolete. Use 'mpz_urandomb' or 'mpz_urandomm'
3012 -- Function: void mpz_random2 (mpz_t ROP, mp_size_t MAX_SIZE)
3013 Generate a random integer of at most MAX_SIZE limbs, with long
3014 strings of zeros and ones in the binary representation. Useful for
3015 testing functions and algorithms, since this kind of random numbers
3016 have proven to be more likely to trigger corner-case bugs.
3017 Negative random numbers are generated when MAX_SIZE is negative.
3019 This function is obsolete. Use 'mpz_rrandomb' instead.
3022 File: gmp.info, Node: Integer Import and Export, Next: Miscellaneous Integer Functions, Prev: Integer Random Numbers, Up: Integer Functions
3024 5.14 Integer Import and Export
3025 ==============================
3027 'mpz_t' variables can be converted to and from arbitrary words of binary
3028 data with the following functions.
3030 -- Function: void mpz_import (mpz_t ROP, size_t COUNT, int ORDER,
3031 size_t SIZE, int ENDIAN, size_t NAILS, const void *OP)
3032 Set ROP from an array of word data at OP.
3034 The parameters specify the format of the data. COUNT many words
3035 are read, each SIZE bytes. ORDER can be 1 for most significant
3036 word first or -1 for least significant first. Within each word
3037 ENDIAN can be 1 for most significant byte first, -1 for least
3038 significant first, or 0 for the native endianness of the host CPU.
3039 The most significant NAILS bits of each word are skipped, this can
3040 be 0 to use the full words.
3042 There is no sign taken from the data, ROP will simply be a positive
3043 integer. An application can handle any sign itself, and apply it
3044 for instance with 'mpz_neg'.
3046 There are no data alignment restrictions on OP, any address is
3049 Here's an example converting an array of 'unsigned long' data, most
3050 significant element first, and host byte order within each value.
3052 unsigned long a[20];
3053 /* Initialize Z and A */
3054 mpz_import (z, 20, 1, sizeof(a[0]), 0, 0, a);
3056 This example assumes the full 'sizeof' bytes are used for data in
3057 the given type, which is usually true, and certainly true for
3058 'unsigned long' everywhere we know of. However on Cray vector
3059 systems it may be noted that 'short' and 'int' are always stored in
3060 8 bytes (and with 'sizeof' indicating that) but use only 32 or 46
3061 bits. The NAILS feature can account for this, by passing for
3062 instance '8*sizeof(int)-INT_BIT'.
3064 -- Function: void * mpz_export (void *ROP, size_t *COUNTP, int ORDER,
3065 size_t SIZE, int ENDIAN, size_t NAILS, const mpz_t OP)
3066 Fill ROP with word data from OP.
3068 The parameters specify the format of the data produced. Each word
3069 will be SIZE bytes and ORDER can be 1 for most significant word
3070 first or -1 for least significant first. Within each word ENDIAN
3071 can be 1 for most significant byte first, -1 for least significant
3072 first, or 0 for the native endianness of the host CPU. The most
3073 significant NAILS bits of each word are unused and set to zero,
3074 this can be 0 to produce full words.
3076 The number of words produced is written to '*COUNTP', or COUNTP can
3077 be 'NULL' to discard the count. ROP must have enough space for the
3078 data, or if ROP is 'NULL' then a result array of the necessary size
3079 is allocated using the current GMP allocation function (*note
3080 Custom Allocation::). In either case the return value is the
3081 destination used, either ROP or the allocated block.
3083 If OP is non-zero then the most significant word produced will be
3084 non-zero. If OP is zero then the count returned will be zero and
3085 nothing written to ROP. If ROP is 'NULL' in this case, no block is
3086 allocated, just 'NULL' is returned.
3088 The sign of OP is ignored, just the absolute value is exported. An
3089 application can use 'mpz_sgn' to get the sign and handle it as
3090 desired. (*note Integer Comparisons::)
3092 There are no data alignment restrictions on ROP, any address is
3095 When an application is allocating space itself the required size
3096 can be determined with a calculation like the following. Since
3097 'mpz_sizeinbase' always returns at least 1, 'count' here will be at
3098 least one, which avoids any portability problems with 'malloc(0)',
3099 though if 'z' is zero no space at all is actually needed (or
3102 numb = 8*size - nail;
3103 count = (mpz_sizeinbase (z, 2) + numb-1) / numb;
3104 p = malloc (count * size);
3107 File: gmp.info, Node: Miscellaneous Integer Functions, Next: Integer Special Functions, Prev: Integer Import and Export, Up: Integer Functions
3109 5.15 Miscellaneous Functions
3110 ============================
3112 -- Function: int mpz_fits_ulong_p (const mpz_t OP)
3113 -- Function: int mpz_fits_slong_p (const mpz_t OP)
3114 -- Function: int mpz_fits_uint_p (const mpz_t OP)
3115 -- Function: int mpz_fits_sint_p (const mpz_t OP)
3116 -- Function: int mpz_fits_ushort_p (const mpz_t OP)
3117 -- Function: int mpz_fits_sshort_p (const mpz_t OP)
3118 Return non-zero iff the value of OP fits in an 'unsigned long int',
3119 'signed long int', 'unsigned int', 'signed int', 'unsigned short
3120 int', or 'signed short int', respectively. Otherwise, return zero.
3122 -- Macro: int mpz_odd_p (const mpz_t OP)
3123 -- Macro: int mpz_even_p (const mpz_t OP)
3124 Determine whether OP is odd or even, respectively. Return non-zero
3125 if yes, zero if no. These macros evaluate their argument more than
3128 -- Function: size_t mpz_sizeinbase (const mpz_t OP, int BASE)
3129 Return the size of OP measured in number of digits in the given
3130 BASE. BASE can vary from 2 to 62. The sign of OP is ignored, just
3131 the absolute value is used. The result will be either exact or 1
3132 too big. If BASE is a power of 2, the result is always exact. If
3133 OP is zero the return value is always 1.
3135 This function can be used to determine the space required when
3136 converting OP to a string. The right amount of allocation is
3137 normally two more than the value returned by 'mpz_sizeinbase', one
3138 extra for a minus sign and one for the null-terminator.
3140 It will be noted that 'mpz_sizeinbase(OP,2)' can be used to locate
3141 the most significant 1 bit in OP, counting from 1. (Unlike the
3142 bitwise functions which start from 0, *Note Logical and Bit
3143 Manipulation Functions: Integer Logic and Bit Fiddling.)
3146 File: gmp.info, Node: Integer Special Functions, Prev: Miscellaneous Integer Functions, Up: Integer Functions
3148 5.16 Special Functions
3149 ======================
3151 The functions in this section are for various special purposes. Most
3152 applications will not need them.
3154 -- Function: void mpz_array_init (mpz_t INTEGER_ARRAY, mp_size_t
3155 ARRAY_SIZE, mp_size_t FIXED_NUM_BITS)
3156 *This is an obsolete function. Do not use it.*
3158 -- Function: void * _mpz_realloc (mpz_t INTEGER, mp_size_t NEW_ALLOC)
3159 Change the space for INTEGER to NEW_ALLOC limbs. The value in
3160 INTEGER is preserved if it fits, or is set to 0 if not. The return
3161 value is not useful to applications and should be ignored.
3163 'mpz_realloc2' is the preferred way to accomplish allocation
3164 changes like this. 'mpz_realloc2' and '_mpz_realloc' are the same
3165 except that '_mpz_realloc' takes its size in limbs.
3167 -- Function: mp_limb_t mpz_getlimbn (const mpz_t OP, mp_size_t N)
3168 Return limb number N from OP. The sign of OP is ignored, just the
3169 absolute value is used. The least significant limb is number 0.
3171 'mpz_size' can be used to find how many limbs make up OP.
3172 'mpz_getlimbn' returns zero if N is outside the range 0 to
3175 -- Function: size_t mpz_size (const mpz_t OP)
3176 Return the size of OP measured in number of limbs. If OP is zero,
3177 the returned value will be zero.
3179 -- Function: const mp_limb_t * mpz_limbs_read (const mpz_t X)
3180 Return a pointer to the limb array representing the absolute value
3181 of X. The size of the array is 'mpz_size(X)'. Intended for read
3184 -- Function: mp_limb_t * mpz_limbs_write (mpz_t X, mp_size_t N)
3185 -- Function: mp_limb_t * mpz_limbs_modify (mpz_t X, mp_size_t N)
3186 Return a pointer to the limb array, intended for write access. The
3187 array is reallocated as needed, to make room for N limbs. Requires
3188 N > 0. The 'mpz_limbs_modify' function returns an array that holds
3189 the old absolute value of X, while 'mpz_limbs_write' may destroy
3190 the old value and return an array with unspecified contents.
3192 -- Function: void mpz_limbs_finish (mpz_t X, mp_size_t S)
3193 Updates the internal size field of X. Used after writing to the
3194 limb array pointer returned by 'mpz_limbs_write' or
3195 'mpz_limbs_modify' is completed. The array should contain abs(S)
3196 valid limbs, representing the new absolute value for X, and the
3197 sign of X is taken from the sign of S. This function never
3198 reallocates X, so the limb pointer remains valid.
3206 xp = mpz_limbs_modify (x, 2*n);
3207 for (i = 0; i < n; i++)
3208 xp[n+i] = xp[n-1-i];
3209 mpz_limbs_finish (x, mpz_sgn (x) < 0 ? - 2*n : 2*n);
3212 -- Function: mpz_srcptr mpz_roinit_n (mpz_t X, const mp_limb_t *XP,
3214 Special initialization of X, using the given limb array and size.
3215 X should be treated as read-only: it can be passed safely as input
3216 to any mpz function, but not as an output. The array XP must point
3217 to at least a readable limb, its size is abs(XS), and the sign of X
3218 is the sign of XS. For convenience, the function returns X, but
3219 cast to a const pointer type.
3223 static const mp_limb_t y[3] = { 0x1, 0x2, 0x3 };
3225 mpz_add (x, x, mpz_roinit_n (tmp, y, 3));
3228 -- Macro: mpz_t MPZ_ROINIT_N (mp_limb_t *XP, mp_size_t XS)
3229 This macro expands to an initializer which can be assigned to an
3230 mpz_t variable. The limb array XP must point to at least a
3231 readable limb, moreover, unlike the 'mpz_roinit_n' function, the
3232 array must be normalized: if XS is non-zero, then 'XP[abs(XS)-1]'
3233 must be non-zero. Intended primarily for constant values. Using
3234 it for non-constant values requires a C compiler supporting C99.
3238 static const mp_limb_t ya[3] = { 0x1, 0x2, 0x3 };
3239 static const mpz_t y = MPZ_ROINIT_N ((mp_limb_t *) ya, 3);
3245 File: gmp.info, Node: Rational Number Functions, Next: Floating-point Functions, Prev: Integer Functions, Up: Top
3247 6 Rational Number Functions
3248 ***************************
3250 This chapter describes the GMP functions for performing arithmetic on
3251 rational numbers. These functions start with the prefix 'mpq_'.
3253 Rational numbers are stored in objects of type 'mpq_t'.
3255 All rational arithmetic functions assume operands have a canonical
3256 form, and canonicalize their result. The canonical form means that the
3257 denominator and the numerator have no common factors, and that the
3258 denominator is positive. Zero has the unique representation 0/1.
3260 Pure assignment functions do not canonicalize the assigned variable.
3261 It is the responsibility of the user to canonicalize the assigned
3262 variable before any arithmetic operations are performed on that
3265 -- Function: void mpq_canonicalize (mpq_t OP)
3266 Remove any factors that are common to the numerator and denominator
3267 of OP, and make the denominator positive.
3271 * Initializing Rationals::
3272 * Rational Conversions::
3273 * Rational Arithmetic::
3274 * Comparing Rationals::
3275 * Applying Integer Functions::
3276 * I/O of Rationals::
3279 File: gmp.info, Node: Initializing Rationals, Next: Rational Conversions, Prev: Rational Number Functions, Up: Rational Number Functions
3281 6.1 Initialization and Assignment Functions
3282 ===========================================
3284 -- Function: void mpq_init (mpq_t X)
3285 Initialize X and set it to 0/1. Each variable should normally only
3286 be initialized once, or at least cleared out (using the function
3287 'mpq_clear') between each initialization.
3289 -- Function: void mpq_inits (mpq_t X, ...)
3290 Initialize a NULL-terminated list of 'mpq_t' variables, and set
3291 their values to 0/1.
3293 -- Function: void mpq_clear (mpq_t X)
3294 Free the space occupied by X. Make sure to call this function for
3295 all 'mpq_t' variables when you are done with them.
3297 -- Function: void mpq_clears (mpq_t X, ...)
3298 Free the space occupied by a NULL-terminated list of 'mpq_t'
3301 -- Function: void mpq_set (mpq_t ROP, const mpq_t OP)
3302 -- Function: void mpq_set_z (mpq_t ROP, const mpz_t OP)
3305 -- Function: void mpq_set_ui (mpq_t ROP, unsigned long int OP1,
3306 unsigned long int OP2)
3307 -- Function: void mpq_set_si (mpq_t ROP, signed long int OP1, unsigned
3309 Set the value of ROP to OP1/OP2. Note that if OP1 and OP2 have
3310 common factors, ROP has to be passed to 'mpq_canonicalize' before
3311 any operations are performed on ROP.
3313 -- Function: int mpq_set_str (mpq_t ROP, const char *STR, int BASE)
3314 Set ROP from a null-terminated string STR in the given BASE.
3316 The string can be an integer like "41" or a fraction like "41/152".
3317 The fraction must be in canonical form (*note Rational Number
3318 Functions::), or if not then 'mpq_canonicalize' must be called.
3320 The numerator and optional denominator are parsed the same as in
3321 'mpz_set_str' (*note Assigning Integers::). White space is allowed
3322 in the string, and is simply ignored. The BASE can vary from 2 to
3323 62, or if BASE is 0 then the leading characters are used: '0x' or
3324 '0X' for hex, '0b' or '0B' for binary, '0' for octal, or decimal
3325 otherwise. Note that this is done separately for the numerator and
3326 denominator, so for instance '0xEF/100' is 239/100, whereas
3327 '0xEF/0x100' is 239/256.
3329 The return value is 0 if the entire string is a valid number, or -1
3332 -- Function: void mpq_swap (mpq_t ROP1, mpq_t ROP2)
3333 Swap the values ROP1 and ROP2 efficiently.
3336 File: gmp.info, Node: Rational Conversions, Next: Rational Arithmetic, Prev: Initializing Rationals, Up: Rational Number Functions
3338 6.2 Conversion Functions
3339 ========================
3341 -- Function: double mpq_get_d (const mpq_t OP)
3342 Convert OP to a 'double', truncating if necessary (i.e. rounding
3345 If the exponent from the conversion is too big or too small to fit
3346 a 'double' then the result is system dependent. For too big an
3347 infinity is returned when available. For too small 0.0 is normally
3348 returned. Hardware overflow, underflow and denorm traps may or may
3351 -- Function: void mpq_set_d (mpq_t ROP, double OP)
3352 -- Function: void mpq_set_f (mpq_t ROP, const mpf_t OP)
3353 Set ROP to the value of OP. There is no rounding, this conversion
3356 -- Function: char * mpq_get_str (char *STR, int BASE, const mpq_t OP)
3357 Convert OP to a string of digits in base BASE. The base argument
3358 may vary from 2 to 62 or from -2 to -36. The string will be of the
3359 form 'num/den', or if the denominator is 1 then just 'num'.
3361 For BASE in the range 2..36, digits and lower-case letters are
3362 used; for -2..-36, digits and upper-case letters are used; for
3363 37..62, digits, upper-case letters, and lower-case letters (in that
3364 significance order) are used.
3366 If STR is 'NULL', the result string is allocated using the current
3367 allocation function (*note Custom Allocation::). The block will be
3368 'strlen(str)+1' bytes, that being exactly enough for the string and
3371 If STR is not 'NULL', it should point to a block of storage large
3372 enough for the result, that being
3374 mpz_sizeinbase (mpq_numref(OP), BASE)
3375 + mpz_sizeinbase (mpq_denref(OP), BASE) + 3
3377 The three extra bytes are for a possible minus sign, possible
3378 slash, and the null-terminator.
3380 A pointer to the result string is returned, being either the
3381 allocated block, or the given STR.
3384 File: gmp.info, Node: Rational Arithmetic, Next: Comparing Rationals, Prev: Rational Conversions, Up: Rational Number Functions
3386 6.3 Arithmetic Functions
3387 ========================
3389 -- Function: void mpq_add (mpq_t SUM, const mpq_t ADDEND1, const mpq_t
3391 Set SUM to ADDEND1 + ADDEND2.
3393 -- Function: void mpq_sub (mpq_t DIFFERENCE, const mpq_t MINUEND, const
3395 Set DIFFERENCE to MINUEND - SUBTRAHEND.
3397 -- Function: void mpq_mul (mpq_t PRODUCT, const mpq_t MULTIPLIER, const
3399 Set PRODUCT to MULTIPLIER times MULTIPLICAND.
3401 -- Function: void mpq_mul_2exp (mpq_t ROP, const mpq_t OP1, mp_bitcnt_t
3403 Set ROP to OP1 times 2 raised to OP2.
3405 -- Function: void mpq_div (mpq_t QUOTIENT, const mpq_t DIVIDEND, const
3407 Set QUOTIENT to DIVIDEND/DIVISOR.
3409 -- Function: void mpq_div_2exp (mpq_t ROP, const mpq_t OP1, mp_bitcnt_t
3411 Set ROP to OP1 divided by 2 raised to OP2.
3413 -- Function: void mpq_neg (mpq_t NEGATED_OPERAND, const mpq_t OPERAND)
3414 Set NEGATED_OPERAND to -OPERAND.
3416 -- Function: void mpq_abs (mpq_t ROP, const mpq_t OP)
3417 Set ROP to the absolute value of OP.
3419 -- Function: void mpq_inv (mpq_t INVERTED_NUMBER, const mpq_t NUMBER)
3420 Set INVERTED_NUMBER to 1/NUMBER. If the new denominator is zero,
3421 this routine will divide by zero.
3424 File: gmp.info, Node: Comparing Rationals, Next: Applying Integer Functions, Prev: Rational Arithmetic, Up: Rational Number Functions
3426 6.4 Comparison Functions
3427 ========================
3429 -- Function: int mpq_cmp (const mpq_t OP1, const mpq_t OP2)
3430 -- Function: int mpq_cmp_z (const mpq_t OP1, const mpz_t OP2)
3431 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero if
3432 OP1 = OP2, and a negative value if OP1 < OP2.
3434 To determine if two rationals are equal, 'mpq_equal' is faster than
3437 -- Macro: int mpq_cmp_ui (const mpq_t OP1, unsigned long int NUM2,
3438 unsigned long int DEN2)
3439 -- Macro: int mpq_cmp_si (const mpq_t OP1, long int NUM2, unsigned long
3441 Compare OP1 and NUM2/DEN2. Return a positive value if OP1 >
3442 NUM2/DEN2, zero if OP1 = NUM2/DEN2, and a negative value if OP1 <
3445 NUM2 and DEN2 are allowed to have common factors.
3447 These functions are implemented as a macros and evaluate their
3448 arguments multiple times.
3450 -- Macro: int mpq_sgn (const mpq_t OP)
3451 Return +1 if OP > 0, 0 if OP = 0, and -1 if OP < 0.
3453 This function is actually implemented as a macro. It evaluates its
3454 argument multiple times.
3456 -- Function: int mpq_equal (const mpq_t OP1, const mpq_t OP2)
3457 Return non-zero if OP1 and OP2 are equal, zero if they are
3458 non-equal. Although 'mpq_cmp' can be used for the same purpose,
3459 this function is much faster.
3462 File: gmp.info, Node: Applying Integer Functions, Next: I/O of Rationals, Prev: Comparing Rationals, Up: Rational Number Functions
3464 6.5 Applying Integer Functions to Rationals
3465 ===========================================
3467 The set of 'mpq' functions is quite small. In particular, there are few
3468 functions for either input or output. The following functions give
3469 direct access to the numerator and denominator of an 'mpq_t'.
3471 Note that if an assignment to the numerator and/or denominator could
3472 take an 'mpq_t' out of the canonical form described at the start of this
3473 chapter (*note Rational Number Functions::) then 'mpq_canonicalize' must
3474 be called before any other 'mpq' functions are applied to that 'mpq_t'.
3476 -- Macro: mpz_t mpq_numref (const mpq_t OP)
3477 -- Macro: mpz_t mpq_denref (const mpq_t OP)
3478 Return a reference to the numerator and denominator of OP,
3479 respectively. The 'mpz' functions can be used on the result of
3482 -- Function: void mpq_get_num (mpz_t NUMERATOR, const mpq_t RATIONAL)
3483 -- Function: void mpq_get_den (mpz_t DENOMINATOR, const mpq_t RATIONAL)
3484 -- Function: void mpq_set_num (mpq_t RATIONAL, const mpz_t NUMERATOR)
3485 -- Function: void mpq_set_den (mpq_t RATIONAL, const mpz_t DENOMINATOR)
3486 Get or set the numerator or denominator of a rational. These
3487 functions are equivalent to calling 'mpz_set' with an appropriate
3488 'mpq_numref' or 'mpq_denref'. Direct use of 'mpq_numref' or
3489 'mpq_denref' is recommended instead of these functions.
3492 File: gmp.info, Node: I/O of Rationals, Prev: Applying Integer Functions, Up: Rational Number Functions
3494 6.6 Input and Output Functions
3495 ==============================
3497 Functions that perform input from a stdio stream, and functions that
3498 output to a stdio stream, of 'mpq' numbers. Passing a 'NULL' pointer
3499 for a STREAM argument to any of these functions will make them read from
3500 'stdin' and write to 'stdout', respectively.
3502 When using any of these functions, it is a good idea to include
3503 'stdio.h' before 'gmp.h', since that will allow 'gmp.h' to define
3504 prototypes for these functions.
3506 See also *note Formatted Output:: and *note Formatted Input::.
3508 -- Function: size_t mpq_out_str (FILE *STREAM, int BASE, const mpq_t
3510 Output OP on stdio stream STREAM, as a string of digits in base
3511 BASE. The base argument may vary from 2 to 62 or from -2 to -36.
3512 Output is in the form 'num/den' or if the denominator is 1 then
3515 For BASE in the range 2..36, digits and lower-case letters are
3516 used; for -2..-36, digits and upper-case letters are used; for
3517 37..62, digits, upper-case letters, and lower-case letters (in that
3518 significance order) are used.
3520 Return the number of bytes written, or if an error occurred, return
3523 -- Function: size_t mpq_inp_str (mpq_t ROP, FILE *STREAM, int BASE)
3524 Read a string of digits from STREAM and convert them to a rational
3525 in ROP. Any initial white-space characters are read and discarded.
3526 Return the number of characters read (including white space), or 0
3527 if a rational could not be read.
3529 The input can be a fraction like '17/63' or just an integer like
3530 '123'. Reading stops at the first character not in this form, and
3531 white space is not permitted within the string. If the input might
3532 not be in canonical form, then 'mpq_canonicalize' must be called
3533 (*note Rational Number Functions::).
3535 The BASE can be between 2 and 62, or can be 0 in which case the
3536 leading characters of the string determine the base, '0x' or '0X'
3537 for hexadecimal, '0b' and '0B' for binary, '0' for octal, or
3538 decimal otherwise. The leading characters are examined separately
3539 for the numerator and denominator of a fraction, so for instance
3540 '0x10/11' is 16/11, whereas '0x10/0x11' is 16/17.
3543 File: gmp.info, Node: Floating-point Functions, Next: Low-level Functions, Prev: Rational Number Functions, Up: Top
3545 7 Floating-point Functions
3546 **************************
3548 GMP floating point numbers are stored in objects of type 'mpf_t' and
3549 functions operating on them have an 'mpf_' prefix.
3551 The mantissa of each float has a user-selectable precision, in
3552 practice only limited by available memory. Each variable has its own
3553 precision, and that can be increased or decreased at any time. This
3554 selectable precision is a minimum value, GMP rounds it up to a whole
3557 The accuracy of a calculation is determined by the priorly set
3558 precision of the destination variable and the numeric values of the
3559 input variables. Input variables' set precisions do not affect
3560 calculations (except indirectly as their values might have been affected
3561 when they were assigned).
3563 The exponent of each float has fixed precision, one machine word on
3564 most systems. In the current implementation the exponent is a count of
3565 limbs, so for example on a 32-bit system this means a range of roughly
3566 2^-68719476768 to 2^68719476736, or on a 64-bit system this will be much
3567 greater. Note however that 'mpf_get_str' can only return an exponent
3568 which fits an 'mp_exp_t' and currently 'mpf_set_str' doesn't accept
3569 exponents bigger than a 'long'.
3571 Each variable keeps track of the mantissa data actually in use. This
3572 means that if a float is exactly represented in only a few bits then
3573 only those bits will be used in a calculation, even if the variable's
3574 selected precision is high. This is a performance optimization; it does
3575 not affect the numeric results.
3577 Internally, GMP sometimes calculates with higher precision than that
3578 of the destination variable in order to limit errors. Final results are
3579 always truncated to the destination variable's precision.
3581 The mantissa is stored in binary. One consequence of this is that
3582 decimal fractions like 0.1 cannot be represented exactly. The same is
3583 true of plain IEEE 'double' floats. This makes both highly unsuitable
3584 for calculations involving money or other values that should be exact
3585 decimal fractions. (Suitably scaled integers, or perhaps rationals, are
3588 The 'mpf' functions and variables have no special notion of infinity
3589 or not-a-number, and applications must take care not to overflow the
3590 exponent or results will be unpredictable.
3592 Note that the 'mpf' functions are _not_ intended as a smooth
3593 extension to IEEE P754 arithmetic. In particular results obtained on
3594 one computer often differ from the results on a computer with a
3595 different word size.
3597 New projects should consider using the GMP extension library MPFR
3598 (<http://mpfr.org>) instead. MPFR provides well-defined precision and
3599 accurate rounding, and thereby naturally extends IEEE P754.
3603 * Initializing Floats::
3604 * Assigning Floats::
3605 * Simultaneous Float Init & Assign::
3606 * Converting Floats::
3607 * Float Arithmetic::
3608 * Float Comparison::
3610 * Miscellaneous Float Functions::
3613 File: gmp.info, Node: Initializing Floats, Next: Assigning Floats, Prev: Floating-point Functions, Up: Floating-point Functions
3615 7.1 Initialization Functions
3616 ============================
3618 -- Function: void mpf_set_default_prec (mp_bitcnt_t PREC)
3619 Set the default precision to be *at least* PREC bits. All
3620 subsequent calls to 'mpf_init' will use this precision, but
3621 previously initialized variables are unaffected.
3623 -- Function: mp_bitcnt_t mpf_get_default_prec (void)
3624 Return the default precision actually used.
3626 An 'mpf_t' object must be initialized before storing the first value
3627 in it. The functions 'mpf_init' and 'mpf_init2' are used for that
3630 -- Function: void mpf_init (mpf_t X)
3631 Initialize X to 0. Normally, a variable should be initialized once
3632 only or at least be cleared, using 'mpf_clear', between
3633 initializations. The precision of X is undefined unless a default
3634 precision has already been established by a call to
3635 'mpf_set_default_prec'.
3637 -- Function: void mpf_init2 (mpf_t X, mp_bitcnt_t PREC)
3638 Initialize X to 0 and set its precision to be *at least* PREC bits.
3639 Normally, a variable should be initialized once only or at least be
3640 cleared, using 'mpf_clear', between initializations.
3642 -- Function: void mpf_inits (mpf_t X, ...)
3643 Initialize a NULL-terminated list of 'mpf_t' variables, and set
3644 their values to 0. The precision of the initialized variables is
3645 undefined unless a default precision has already been established
3646 by a call to 'mpf_set_default_prec'.
3648 -- Function: void mpf_clear (mpf_t X)
3649 Free the space occupied by X. Make sure to call this function for
3650 all 'mpf_t' variables when you are done with them.
3652 -- Function: void mpf_clears (mpf_t X, ...)
3653 Free the space occupied by a NULL-terminated list of 'mpf_t'
3656 Here is an example on how to initialize floating-point variables:
3659 mpf_init (x); /* use default precision */
3660 mpf_init2 (y, 256); /* precision _at least_ 256 bits */
3662 /* Unless the program is about to exit, do ... */
3667 The following three functions are useful for changing the precision
3668 during a calculation. A typical use would be for adjusting the
3669 precision gradually in iterative algorithms like Newton-Raphson, making
3670 the computation precision closely match the actual accurate part of the
3673 -- Function: mp_bitcnt_t mpf_get_prec (const mpf_t OP)
3674 Return the current precision of OP, in bits.
3676 -- Function: void mpf_set_prec (mpf_t ROP, mp_bitcnt_t PREC)
3677 Set the precision of ROP to be *at least* PREC bits. The value in
3678 ROP will be truncated to the new precision.
3680 This function requires a call to 'realloc', and so should not be
3681 used in a tight loop.
3683 -- Function: void mpf_set_prec_raw (mpf_t ROP, mp_bitcnt_t PREC)
3684 Set the precision of ROP to be *at least* PREC bits, without
3685 changing the memory allocated.
3687 PREC must be no more than the allocated precision for ROP, that
3688 being the precision when ROP was initialized, or in the most recent
3691 The value in ROP is unchanged, and in particular if it had a higher
3692 precision than PREC it will retain that higher precision. New
3693 values written to ROP will use the new PREC.
3695 Before calling 'mpf_clear' or the full 'mpf_set_prec', another
3696 'mpf_set_prec_raw' call must be made to restore ROP to its original
3697 allocated precision. Failing to do so will have unpredictable
3700 'mpf_get_prec' can be used before 'mpf_set_prec_raw' to get the
3701 original allocated precision. After 'mpf_set_prec_raw' it reflects
3704 'mpf_set_prec_raw' is an efficient way to use an 'mpf_t' variable
3705 at different precisions during a calculation, perhaps to gradually
3706 increase precision in an iteration, or just to use various
3707 different precisions for different purposes during a calculation.
3710 File: gmp.info, Node: Assigning Floats, Next: Simultaneous Float Init & Assign, Prev: Initializing Floats, Up: Floating-point Functions
3712 7.2 Assignment Functions
3713 ========================
3715 These functions assign new values to already initialized floats (*note
3716 Initializing Floats::).
3718 -- Function: void mpf_set (mpf_t ROP, const mpf_t OP)
3719 -- Function: void mpf_set_ui (mpf_t ROP, unsigned long int OP)
3720 -- Function: void mpf_set_si (mpf_t ROP, signed long int OP)
3721 -- Function: void mpf_set_d (mpf_t ROP, double OP)
3722 -- Function: void mpf_set_z (mpf_t ROP, const mpz_t OP)
3723 -- Function: void mpf_set_q (mpf_t ROP, const mpq_t OP)
3724 Set the value of ROP from OP.
3726 -- Function: int mpf_set_str (mpf_t ROP, const char *STR, int BASE)
3727 Set the value of ROP from the string in STR. The string is of the
3728 form 'M@N' or, if the base is 10 or less, alternatively 'MeN'. 'M'
3729 is the mantissa and 'N' is the exponent. The mantissa is always in
3730 the specified base. The exponent is either in the specified base
3731 or, if BASE is negative, in decimal. The decimal point expected is
3732 taken from the current locale, on systems providing 'localeconv'.
3734 The argument BASE may be in the ranges 2 to 62, or -62 to -2.
3735 Negative values are used to specify that the exponent is in
3738 For bases up to 36, case is ignored; upper-case and lower-case
3739 letters have the same value; for bases 37 to 62, upper-case letter
3740 represent the usual 10..35 while lower-case letter represent
3743 Unlike the corresponding 'mpz' function, the base will not be
3744 determined from the leading characters of the string if BASE is 0.
3745 This is so that numbers like '0.23' are not interpreted as octal.
3747 White space is allowed in the string, and is simply ignored. [This
3748 is not really true; white-space is ignored in the beginning of the
3749 string and within the mantissa, but not in other places, such as
3750 after a minus sign or in the exponent. We are considering changing
3751 the definition of this function, making it fail when there is any
3752 white-space in the input, since that makes a lot of sense. Please
3753 tell us your opinion about this change. Do you really want it to
3754 accept "3 14" as meaning 314 as it does now?]
3756 This function returns 0 if the entire string is a valid number in
3757 base BASE. Otherwise it returns -1.
3759 -- Function: void mpf_swap (mpf_t ROP1, mpf_t ROP2)
3760 Swap ROP1 and ROP2 efficiently. Both the values and the precisions
3761 of the two variables are swapped.
3764 File: gmp.info, Node: Simultaneous Float Init & Assign, Next: Converting Floats, Prev: Assigning Floats, Up: Floating-point Functions
3766 7.3 Combined Initialization and Assignment Functions
3767 ====================================================
3769 For convenience, GMP provides a parallel series of initialize-and-set
3770 functions which initialize the output and then store the value there.
3771 These functions' names have the form 'mpf_init_set...'
3773 Once the float has been initialized by any of the 'mpf_init_set...'
3774 functions, it can be used as the source or destination operand for the
3775 ordinary float functions. Don't use an initialize-and-set function on a
3776 variable already initialized!
3778 -- Function: void mpf_init_set (mpf_t ROP, const mpf_t OP)
3779 -- Function: void mpf_init_set_ui (mpf_t ROP, unsigned long int OP)
3780 -- Function: void mpf_init_set_si (mpf_t ROP, signed long int OP)
3781 -- Function: void mpf_init_set_d (mpf_t ROP, double OP)
3782 Initialize ROP and set its value from OP.
3784 The precision of ROP will be taken from the active default
3785 precision, as set by 'mpf_set_default_prec'.
3787 -- Function: int mpf_init_set_str (mpf_t ROP, const char *STR, int
3789 Initialize ROP and set its value from the string in STR. See
3790 'mpf_set_str' above for details on the assignment operation.
3792 Note that ROP is initialized even if an error occurs. (I.e., you
3793 have to call 'mpf_clear' for it.)
3795 The precision of ROP will be taken from the active default
3796 precision, as set by 'mpf_set_default_prec'.
3799 File: gmp.info, Node: Converting Floats, Next: Float Arithmetic, Prev: Simultaneous Float Init & Assign, Up: Floating-point Functions
3801 7.4 Conversion Functions
3802 ========================
3804 -- Function: double mpf_get_d (const mpf_t OP)
3805 Convert OP to a 'double', truncating if necessary (i.e. rounding
3808 If the exponent in OP is too big or too small to fit a 'double'
3809 then the result is system dependent. For too big an infinity is
3810 returned when available. For too small 0.0 is normally returned.
3811 Hardware overflow, underflow and denorm traps may or may not occur.
3813 -- Function: double mpf_get_d_2exp (signed long int *EXP, const mpf_t
3815 Convert OP to a 'double', truncating if necessary (i.e. rounding
3816 towards zero), and with an exponent returned separately.
3818 The return value is in the range 0.5<=abs(D)<1 and the exponent is
3819 stored to '*EXP'. D * 2^EXP is the (truncated) OP value. If OP is
3820 zero, the return is 0.0 and 0 is stored to '*EXP'.
3822 This is similar to the standard C 'frexp' function (*note
3823 (libc)Normalization Functions::).
3825 -- Function: long mpf_get_si (const mpf_t OP)
3826 -- Function: unsigned long mpf_get_ui (const mpf_t OP)
3827 Convert OP to a 'long' or 'unsigned long', truncating any fraction
3828 part. If OP is too big for the return type, the result is
3831 See also 'mpf_fits_slong_p' and 'mpf_fits_ulong_p' (*note
3832 Miscellaneous Float Functions::).
3834 -- Function: char * mpf_get_str (char *STR, mp_exp_t *EXPPTR, int BASE,
3835 size_t N_DIGITS, const mpf_t OP)
3836 Convert OP to a string of digits in base BASE. The base argument
3837 may vary from 2 to 62 or from -2 to -36. Up to N_DIGITS digits
3838 will be generated. Trailing zeros are not returned. No more
3839 digits than can be accurately represented by OP are ever generated.
3840 If N_DIGITS is 0 then that accurate maximum number of digits are
3843 For BASE in the range 2..36, digits and lower-case letters are
3844 used; for -2..-36, digits and upper-case letters are used; for
3845 37..62, digits, upper-case letters, and lower-case letters (in that
3846 significance order) are used.
3848 If STR is 'NULL', the result string is allocated using the current
3849 allocation function (*note Custom Allocation::). The block will be
3850 'strlen(str)+1' bytes, that being exactly enough for the string and
3853 If STR is not 'NULL', it should point to a block of N_DIGITS + 2
3854 bytes, that being enough for the mantissa, a possible minus sign,
3855 and a null-terminator. When N_DIGITS is 0 to get all significant
3856 digits, an application won't be able to know the space required,
3857 and STR should be 'NULL' in that case.
3859 The generated string is a fraction, with an implicit radix point
3860 immediately to the left of the first digit. The applicable
3861 exponent is written through the EXPPTR pointer. For example, the
3862 number 3.1416 would be returned as string "31416" and exponent 1.
3864 When OP is zero, an empty string is produced and the exponent
3867 A pointer to the result string is returned, being either the
3868 allocated block or the given STR.
3871 File: gmp.info, Node: Float Arithmetic, Next: Float Comparison, Prev: Converting Floats, Up: Floating-point Functions
3873 7.5 Arithmetic Functions
3874 ========================
3876 -- Function: void mpf_add (mpf_t ROP, const mpf_t OP1, const mpf_t OP2)
3877 -- Function: void mpf_add_ui (mpf_t ROP, const mpf_t OP1, unsigned long
3879 Set ROP to OP1 + OP2.
3881 -- Function: void mpf_sub (mpf_t ROP, const mpf_t OP1, const mpf_t OP2)
3882 -- Function: void mpf_ui_sub (mpf_t ROP, unsigned long int OP1, const
3884 -- Function: void mpf_sub_ui (mpf_t ROP, const mpf_t OP1, unsigned long
3886 Set ROP to OP1 - OP2.
3888 -- Function: void mpf_mul (mpf_t ROP, const mpf_t OP1, const mpf_t OP2)
3889 -- Function: void mpf_mul_ui (mpf_t ROP, const mpf_t OP1, unsigned long
3891 Set ROP to OP1 times OP2.
3893 Division is undefined if the divisor is zero, and passing a zero
3894 divisor to the divide functions will make these functions intentionally
3895 divide by zero. This lets the user handle arithmetic exceptions in
3896 these functions in the same manner as other arithmetic exceptions.
3898 -- Function: void mpf_div (mpf_t ROP, const mpf_t OP1, const mpf_t OP2)
3899 -- Function: void mpf_ui_div (mpf_t ROP, unsigned long int OP1, const
3901 -- Function: void mpf_div_ui (mpf_t ROP, const mpf_t OP1, unsigned long
3905 -- Function: void mpf_sqrt (mpf_t ROP, const mpf_t OP)
3906 -- Function: void mpf_sqrt_ui (mpf_t ROP, unsigned long int OP)
3907 Set ROP to the square root of OP.
3909 -- Function: void mpf_pow_ui (mpf_t ROP, const mpf_t OP1, unsigned long
3911 Set ROP to OP1 raised to the power OP2.
3913 -- Function: void mpf_neg (mpf_t ROP, const mpf_t OP)
3916 -- Function: void mpf_abs (mpf_t ROP, const mpf_t OP)
3917 Set ROP to the absolute value of OP.
3919 -- Function: void mpf_mul_2exp (mpf_t ROP, const mpf_t OP1, mp_bitcnt_t
3921 Set ROP to OP1 times 2 raised to OP2.
3923 -- Function: void mpf_div_2exp (mpf_t ROP, const mpf_t OP1, mp_bitcnt_t
3925 Set ROP to OP1 divided by 2 raised to OP2.
3928 File: gmp.info, Node: Float Comparison, Next: I/O of Floats, Prev: Float Arithmetic, Up: Floating-point Functions
3930 7.6 Comparison Functions
3931 ========================
3933 -- Function: int mpf_cmp (const mpf_t OP1, const mpf_t OP2)
3934 -- Function: int mpf_cmp_z (const mpf_t OP1, const mpz_t OP2)
3935 -- Function: int mpf_cmp_d (const mpf_t OP1, double OP2)
3936 -- Function: int mpf_cmp_ui (const mpf_t OP1, unsigned long int OP2)
3937 -- Function: int mpf_cmp_si (const mpf_t OP1, signed long int OP2)
3938 Compare OP1 and OP2. Return a positive value if OP1 > OP2, zero if
3939 OP1 = OP2, and a negative value if OP1 < OP2.
3941 'mpf_cmp_d' can be called with an infinity, but results are
3942 undefined for a NaN.
3944 -- Function: int mpf_eq (const mpf_t OP1, const mpf_t OP2, mp_bitcnt_t
3946 *This function is mathematically ill-defined and should not be
3949 Return non-zero if the first OP3 bits of OP1 and OP2 are equal,
3950 zero otherwise. Note that numbers like e.g., 256 (binary
3951 100000000) and 255 (binary 11111111) will never be equal by this
3952 function's measure, and furthermore that 0 will only be equal to
3955 -- Function: void mpf_reldiff (mpf_t ROP, const mpf_t OP1, const mpf_t
3957 Compute the relative difference between OP1 and OP2 and store the
3958 result in ROP. This is abs(OP1-OP2)/OP1.
3960 -- Macro: int mpf_sgn (const mpf_t OP)
3961 Return +1 if OP > 0, 0 if OP = 0, and -1 if OP < 0.
3963 This function is actually implemented as a macro. It evaluates its
3964 argument multiple times.
3967 File: gmp.info, Node: I/O of Floats, Next: Miscellaneous Float Functions, Prev: Float Comparison, Up: Floating-point Functions
3969 7.7 Input and Output Functions
3970 ==============================
3972 Functions that perform input from a stdio stream, and functions that
3973 output to a stdio stream, of 'mpf' numbers. Passing a 'NULL' pointer
3974 for a STREAM argument to any of these functions will make them read from
3975 'stdin' and write to 'stdout', respectively.
3977 When using any of these functions, it is a good idea to include
3978 'stdio.h' before 'gmp.h', since that will allow 'gmp.h' to define
3979 prototypes for these functions.
3981 See also *note Formatted Output:: and *note Formatted Input::.
3983 -- Function: size_t mpf_out_str (FILE *STREAM, int BASE, size_t
3984 N_DIGITS, const mpf_t OP)
3985 Print OP to STREAM, as a string of digits. Return the number of
3986 bytes written, or if an error occurred, return 0.
3988 The mantissa is prefixed with an '0.' and is in the given BASE,
3989 which may vary from 2 to 62 or from -2 to -36. An exponent is then
3990 printed, separated by an 'e', or if the base is greater than 10
3991 then by an '@'. The exponent is always in decimal. The decimal
3992 point follows the current locale, on systems providing
3995 For BASE in the range 2..36, digits and lower-case letters are
3996 used; for -2..-36, digits and upper-case letters are used; for
3997 37..62, digits, upper-case letters, and lower-case letters (in that
3998 significance order) are used.
4000 Up to N_DIGITS will be printed from the mantissa, except that no
4001 more digits than are accurately representable by OP will be
4002 printed. N_DIGITS can be 0 to select that accurate maximum.
4004 -- Function: size_t mpf_inp_str (mpf_t ROP, FILE *STREAM, int BASE)
4005 Read a string in base BASE from STREAM, and put the read float in
4006 ROP. The string is of the form 'M@N' or, if the base is 10 or
4007 less, alternatively 'MeN'. 'M' is the mantissa and 'N' is the
4008 exponent. The mantissa is always in the specified base. The
4009 exponent is either in the specified base or, if BASE is negative,
4010 in decimal. The decimal point expected is taken from the current
4011 locale, on systems providing 'localeconv'.
4013 The argument BASE may be in the ranges 2 to 36, or -36 to -2.
4014 Negative values are used to specify that the exponent is in
4017 Unlike the corresponding 'mpz' function, the base will not be
4018 determined from the leading characters of the string if BASE is 0.
4019 This is so that numbers like '0.23' are not interpreted as octal.
4021 Return the number of bytes read, or if an error occurred, return 0.
4024 File: gmp.info, Node: Miscellaneous Float Functions, Prev: I/O of Floats, Up: Floating-point Functions
4026 7.8 Miscellaneous Functions
4027 ===========================
4029 -- Function: void mpf_ceil (mpf_t ROP, const mpf_t OP)
4030 -- Function: void mpf_floor (mpf_t ROP, const mpf_t OP)
4031 -- Function: void mpf_trunc (mpf_t ROP, const mpf_t OP)
4032 Set ROP to OP rounded to an integer. 'mpf_ceil' rounds to the next
4033 higher integer, 'mpf_floor' to the next lower, and 'mpf_trunc' to
4034 the integer towards zero.
4036 -- Function: int mpf_integer_p (const mpf_t OP)
4037 Return non-zero if OP is an integer.
4039 -- Function: int mpf_fits_ulong_p (const mpf_t OP)
4040 -- Function: int mpf_fits_slong_p (const mpf_t OP)
4041 -- Function: int mpf_fits_uint_p (const mpf_t OP)
4042 -- Function: int mpf_fits_sint_p (const mpf_t OP)
4043 -- Function: int mpf_fits_ushort_p (const mpf_t OP)
4044 -- Function: int mpf_fits_sshort_p (const mpf_t OP)
4045 Return non-zero if OP would fit in the respective C data type, when
4046 truncated to an integer.
4048 -- Function: void mpf_urandomb (mpf_t ROP, gmp_randstate_t STATE,
4050 Generate a uniformly distributed random float in ROP, such that 0
4051 <= ROP < 1, with NBITS significant bits in the mantissa or less if
4052 the precision of ROP is smaller.
4054 The variable STATE must be initialized by calling one of the
4055 'gmp_randinit' functions (*note Random State Initialization::)
4056 before invoking this function.
4058 -- Function: void mpf_random2 (mpf_t ROP, mp_size_t MAX_SIZE, mp_exp_t
4060 Generate a random float of at most MAX_SIZE limbs, with long
4061 strings of zeros and ones in the binary representation. The
4062 exponent of the number is in the interval -EXP to EXP (in limbs).
4063 This function is useful for testing functions and algorithms, since
4064 these kind of random numbers have proven to be more likely to
4065 trigger corner-case bugs. Negative random numbers are generated
4066 when MAX_SIZE is negative.
4069 File: gmp.info, Node: Low-level Functions, Next: Random Number Functions, Prev: Floating-point Functions, Up: Top
4071 8 Low-level Functions
4072 *********************
4074 This chapter describes low-level GMP functions, used to implement the
4075 high-level GMP functions, but also intended for time-critical user code.
4077 These functions start with the prefix 'mpn_'.
4079 The 'mpn' functions are designed to be as fast as possible, *not* to
4080 provide a coherent calling interface. The different functions have
4081 somewhat similar interfaces, but there are variations that make them
4082 hard to use. These functions do as little as possible apart from the
4083 real multiple precision computation, so that no time is spent on things
4084 that not all callers need.
4086 A source operand is specified by a pointer to the least significant
4087 limb and a limb count. A destination operand is specified by just a
4088 pointer. It is the responsibility of the caller to ensure that the
4089 destination has enough space for storing the result.
4091 With this way of specifying operands, it is possible to perform
4092 computations on subranges of an argument, and store the result into a
4093 subrange of a destination.
4095 A common requirement for all functions is that each source area needs
4096 at least one limb. No size argument may be zero. Unless otherwise
4097 stated, in-place operations are allowed where source and destination are
4098 the same, but not where they only partly overlap.
4100 The 'mpn' functions are the base for the implementation of the
4101 'mpz_', 'mpf_', and 'mpq_' functions.
4103 This example adds the number beginning at S1P and the number
4104 beginning at S2P and writes the sum at DESTP. All areas have N limbs.
4106 cy = mpn_add_n (destp, s1p, s2p, n)
4108 It should be noted that the 'mpn' functions make no attempt to
4109 identify high or low zero limbs on their operands, or other special
4110 forms. On random data such cases will be unlikely and it'd be wasteful
4111 for every function to check every time. An application knowing
4112 something about its data can take steps to trim or perhaps split its
4116 In the notation used below, a source operand is identified by the
4117 pointer to the least significant limb, and the limb count in braces.
4118 For example, {S1P, S1N}.
4120 -- Function: mp_limb_t mpn_add_n (mp_limb_t *RP, const mp_limb_t *S1P,
4121 const mp_limb_t *S2P, mp_size_t N)
4122 Add {S1P, N} and {S2P, N}, and write the N least significant limbs
4123 of the result to RP. Return carry, either 0 or 1.
4125 This is the lowest-level function for addition. It is the
4126 preferred function for addition, since it is written in assembly
4127 for most CPUs. For addition of a variable to itself (i.e., S1P
4128 equals S2P) use 'mpn_lshift' with a count of 1 for optimal speed.
4130 -- Function: mp_limb_t mpn_add_1 (mp_limb_t *RP, const mp_limb_t *S1P,
4131 mp_size_t N, mp_limb_t S2LIMB)
4132 Add {S1P, N} and S2LIMB, and write the N least significant limbs of
4133 the result to RP. Return carry, either 0 or 1.
4135 -- Function: mp_limb_t mpn_add (mp_limb_t *RP, const mp_limb_t *S1P,
4136 mp_size_t S1N, const mp_limb_t *S2P, mp_size_t S2N)
4137 Add {S1P, S1N} and {S2P, S2N}, and write the S1N least significant
4138 limbs of the result to RP. Return carry, either 0 or 1.
4140 This function requires that S1N is greater than or equal to S2N.
4142 -- Function: mp_limb_t mpn_sub_n (mp_limb_t *RP, const mp_limb_t *S1P,
4143 const mp_limb_t *S2P, mp_size_t N)
4144 Subtract {S2P, N} from {S1P, N}, and write the N least significant
4145 limbs of the result to RP. Return borrow, either 0 or 1.
4147 This is the lowest-level function for subtraction. It is the
4148 preferred function for subtraction, since it is written in assembly
4151 -- Function: mp_limb_t mpn_sub_1 (mp_limb_t *RP, const mp_limb_t *S1P,
4152 mp_size_t N, mp_limb_t S2LIMB)
4153 Subtract S2LIMB from {S1P, N}, and write the N least significant
4154 limbs of the result to RP. Return borrow, either 0 or 1.
4156 -- Function: mp_limb_t mpn_sub (mp_limb_t *RP, const mp_limb_t *S1P,
4157 mp_size_t S1N, const mp_limb_t *S2P, mp_size_t S2N)
4158 Subtract {S2P, S2N} from {S1P, S1N}, and write the S1N least
4159 significant limbs of the result to RP. Return borrow, either 0 or
4162 This function requires that S1N is greater than or equal to S2N.
4164 -- Function: mp_limb_t mpn_neg (mp_limb_t *RP, const mp_limb_t *SP,
4166 Perform the negation of {SP, N}, and write the result to {RP, N}.
4167 This is equivalent to calling 'mpn_sub_n' with a N-limb zero
4168 minuend and passing {SP, N} as subtrahend. Return borrow, either 0
4171 -- Function: void mpn_mul_n (mp_limb_t *RP, const mp_limb_t *S1P, const
4172 mp_limb_t *S2P, mp_size_t N)
4173 Multiply {S1P, N} and {S2P, N}, and write the 2*N-limb result to
4176 The destination has to have space for 2*N limbs, even if the
4177 product's most significant limb is zero. No overlap is permitted
4178 between the destination and either source.
4180 If the two input operands are the same, use 'mpn_sqr'.
4182 -- Function: mp_limb_t mpn_mul (mp_limb_t *RP, const mp_limb_t *S1P,
4183 mp_size_t S1N, const mp_limb_t *S2P, mp_size_t S2N)
4184 Multiply {S1P, S1N} and {S2P, S2N}, and write the (S1N+S2N)-limb
4185 result to RP. Return the most significant limb of the result.
4187 The destination has to have space for S1N + S2N limbs, even if the
4188 product's most significant limb is zero. No overlap is permitted
4189 between the destination and either source.
4191 This function requires that S1N is greater than or equal to S2N.
4193 -- Function: void mpn_sqr (mp_limb_t *RP, const mp_limb_t *S1P,
4195 Compute the square of {S1P, N} and write the 2*N-limb result to RP.
4197 The destination has to have space for 2N limbs, even if the
4198 result's most significant limb is zero. No overlap is permitted
4199 between the destination and the source.
4201 -- Function: mp_limb_t mpn_mul_1 (mp_limb_t *RP, const mp_limb_t *S1P,
4202 mp_size_t N, mp_limb_t S2LIMB)
4203 Multiply {S1P, N} by S2LIMB, and write the N least significant
4204 limbs of the product to RP. Return the most significant limb of
4205 the product. {S1P, N} and {RP, N} are allowed to overlap provided
4208 This is a low-level function that is a building block for general
4209 multiplication as well as other operations in GMP. It is written
4210 in assembly for most CPUs.
4212 Don't call this function if S2LIMB is a power of 2; use
4213 'mpn_lshift' with a count equal to the logarithm of S2LIMB instead,
4216 -- Function: mp_limb_t mpn_addmul_1 (mp_limb_t *RP, const mp_limb_t
4217 *S1P, mp_size_t N, mp_limb_t S2LIMB)
4218 Multiply {S1P, N} and S2LIMB, and add the N least significant limbs
4219 of the product to {RP, N} and write the result to RP. Return the
4220 most significant limb of the product, plus carry-out from the
4221 addition. {S1P, N} and {RP, N} are allowed to overlap provided RP
4224 This is a low-level function that is a building block for general
4225 multiplication as well as other operations in GMP. It is written
4226 in assembly for most CPUs.
4228 -- Function: mp_limb_t mpn_submul_1 (mp_limb_t *RP, const mp_limb_t
4229 *S1P, mp_size_t N, mp_limb_t S2LIMB)
4230 Multiply {S1P, N} and S2LIMB, and subtract the N least significant
4231 limbs of the product from {RP, N} and write the result to RP.
4232 Return the most significant limb of the product, plus borrow-out
4233 from the subtraction. {S1P, N} and {RP, N} are allowed to overlap
4236 This is a low-level function that is a building block for general
4237 multiplication and division as well as other operations in GMP. It
4238 is written in assembly for most CPUs.
4240 -- Function: void mpn_tdiv_qr (mp_limb_t *QP, mp_limb_t *RP, mp_size_t
4241 QXN, const mp_limb_t *NP, mp_size_t NN, const mp_limb_t *DP,
4243 Divide {NP, NN} by {DP, DN} and put the quotient at {QP, NN-DN+1}
4244 and the remainder at {RP, DN}. The quotient is rounded towards 0.
4246 No overlap is permitted between arguments, except that NP might
4247 equal RP. The dividend size NN must be greater than or equal to
4248 divisor size DN. The most significant limb of the divisor must be
4249 non-zero. The QXN operand must be zero.
4251 -- Function: mp_limb_t mpn_divrem (mp_limb_t *R1P, mp_size_t QXN,
4252 mp_limb_t *RS2P, mp_size_t RS2N, const mp_limb_t *S3P,
4254 [This function is obsolete. Please call 'mpn_tdiv_qr' instead for
4257 Divide {RS2P, RS2N} by {S3P, S3N}, and write the quotient at R1P,
4258 with the exception of the most significant limb, which is returned.
4259 The remainder replaces the dividend at RS2P; it will be S3N limbs
4260 long (i.e., as many limbs as the divisor).
4262 In addition to an integer quotient, QXN fraction limbs are
4263 developed, and stored after the integral limbs. For most usages,
4266 It is required that RS2N is greater than or equal to S3N. It is
4267 required that the most significant bit of the divisor is set.
4269 If the quotient is not needed, pass RS2P + S3N as R1P. Aside from
4270 that special case, no overlap between arguments is permitted.
4272 Return the most significant limb of the quotient, either 0 or 1.
4274 The area at R1P needs to be RS2N - S3N + QXN limbs large.
4276 -- Function: mp_limb_t mpn_divrem_1 (mp_limb_t *R1P, mp_size_t QXN,
4277 mp_limb_t *S2P, mp_size_t S2N, mp_limb_t S3LIMB)
4278 -- Macro: mp_limb_t mpn_divmod_1 (mp_limb_t *R1P, mp_limb_t *S2P,
4279 mp_size_t S2N, mp_limb_t S3LIMB)
4280 Divide {S2P, S2N} by S3LIMB, and write the quotient at R1P. Return
4283 The integer quotient is written to {R1P+QXN, S2N} and in addition
4284 QXN fraction limbs are developed and written to {R1P, QXN}. Either
4285 or both S2N and QXN can be zero. For most usages, QXN will be
4288 'mpn_divmod_1' exists for upward source compatibility and is simply
4289 a macro calling 'mpn_divrem_1' with a QXN of 0.
4291 The areas at R1P and S2P have to be identical or completely
4292 separate, not partially overlapping.
4294 -- Function: mp_limb_t mpn_divmod (mp_limb_t *R1P, mp_limb_t *RS2P,
4295 mp_size_t RS2N, const mp_limb_t *S3P, mp_size_t S3N)
4296 [This function is obsolete. Please call 'mpn_tdiv_qr' instead for
4299 -- Function: void mpn_divexact_1 (mp_limb_t * RP, const mp_limb_t * SP,
4300 mp_size_t N, mp_limb_t D)
4301 Divide {SP, N} by D, expecting it to divide exactly, and writing
4302 the result to {RP, N}. If D doesn't divide exactly, the value
4303 written to {RP, N} is undefined. The areas at RP and SP have to be
4304 identical or completely separate, not partially overlapping.
4306 -- Macro: mp_limb_t mpn_divexact_by3 (mp_limb_t *RP, mp_limb_t *SP,
4308 -- Function: mp_limb_t mpn_divexact_by3c (mp_limb_t *RP, mp_limb_t *SP,
4309 mp_size_t N, mp_limb_t CARRY)
4310 Divide {SP, N} by 3, expecting it to divide exactly, and writing
4311 the result to {RP, N}. If 3 divides exactly, the return value is
4312 zero and the result is the quotient. If not, the return value is
4313 non-zero and the result won't be anything useful.
4315 'mpn_divexact_by3c' takes an initial carry parameter, which can be
4316 the return value from a previous call, so a large calculation can
4317 be done piece by piece from low to high. 'mpn_divexact_by3' is
4318 simply a macro calling 'mpn_divexact_by3c' with a 0 carry
4321 These routines use a multiply-by-inverse and will be faster than
4322 'mpn_divrem_1' on CPUs with fast multiplication but slow division.
4324 The source a, result q, size n, initial carry i, and return value c
4325 satisfy c*b^n + a-i = 3*q, where b=2^GMP_NUMB_BITS. The return c is
4326 always 0, 1 or 2, and the initial carry i must also be 0, 1 or 2
4327 (these are both borrows really). When c=0 clearly q=(a-i)/3. When
4328 c!=0, the remainder (a-i) mod 3 is given by 3-c, because b == 1 mod
4329 3 (when 'mp_bits_per_limb' is even, which is always so currently).
4331 -- Function: mp_limb_t mpn_mod_1 (const mp_limb_t *S1P, mp_size_t S1N,
4333 Divide {S1P, S1N} by S2LIMB, and return the remainder. S1N can be
4336 -- Function: mp_limb_t mpn_lshift (mp_limb_t *RP, const mp_limb_t *SP,
4337 mp_size_t N, unsigned int COUNT)
4338 Shift {SP, N} left by COUNT bits, and write the result to {RP, N}.
4339 The bits shifted out at the left are returned in the least
4340 significant COUNT bits of the return value (the rest of the return
4343 COUNT must be in the range 1 to mp_bits_per_limb-1. The regions
4344 {SP, N} and {RP, N} may overlap, provided RP >= SP.
4346 This function is written in assembly for most CPUs.
4348 -- Function: mp_limb_t mpn_rshift (mp_limb_t *RP, const mp_limb_t *SP,
4349 mp_size_t N, unsigned int COUNT)
4350 Shift {SP, N} right by COUNT bits, and write the result to {RP, N}.
4351 The bits shifted out at the right are returned in the most
4352 significant COUNT bits of the return value (the rest of the return
4355 COUNT must be in the range 1 to mp_bits_per_limb-1. The regions
4356 {SP, N} and {RP, N} may overlap, provided RP <= SP.
4358 This function is written in assembly for most CPUs.
4360 -- Function: int mpn_cmp (const mp_limb_t *S1P, const mp_limb_t *S2P,
4362 Compare {S1P, N} and {S2P, N} and return a positive value if S1 >
4363 S2, 0 if they are equal, or a negative value if S1 < S2.
4365 -- Function: int mpn_zero_p (const mp_limb_t *SP, mp_size_t N)
4366 Test {SP, N} and return 1 if the operand is zero, 0 otherwise.
4368 -- Function: mp_size_t mpn_gcd (mp_limb_t *RP, mp_limb_t *XP, mp_size_t
4369 XN, mp_limb_t *YP, mp_size_t YN)
4370 Set {RP, RETVAL} to the greatest common divisor of {XP, XN} and
4371 {YP, YN}. The result can be up to YN limbs, the return value is
4372 the actual number produced. Both source operands are destroyed.
4374 It is required that XN >= YN > 0, the most significant limb of {YP,
4375 YN} must be non-zero, and at least one of the two operands must be
4376 odd. No overlap is permitted between {XP, XN} and {YP, YN}.
4378 -- Function: mp_limb_t mpn_gcd_1 (const mp_limb_t *XP, mp_size_t XN,
4380 Return the greatest common divisor of {XP, XN} and YLIMB. Both
4381 operands must be non-zero.
4383 -- Function: mp_size_t mpn_gcdext (mp_limb_t *GP, mp_limb_t *SP,
4384 mp_size_t *SN, mp_limb_t *UP, mp_size_t UN, mp_limb_t *VP,
4386 Let U be defined by {UP, UN} and let V be defined by {VP, VN}.
4388 Compute the greatest common divisor G of U and V. Compute a
4389 cofactor S such that G = US + VT. The second cofactor T is not
4390 computed but can easily be obtained from (G - U*S) / V (the
4391 division will be exact). It is required that UN >= VN > 0, and the
4392 most significant limb of {VP, VN} must be non-zero.
4394 S satisfies S = 1 or abs(S) < V / (2 G). S = 0 if and only if V
4395 divides U (i.e., G = V).
4397 Store G at GP and let the return value define its limb count.
4398 Store S at SP and let |*SN| define its limb count. S can be
4399 negative; when this happens *SN will be negative. The area at GP
4400 should have room for VN limbs and the area at SP should have room
4403 Both source operands are destroyed.
4405 Compatibility notes: GMP 4.3.0 and 4.3.1 defined S less strictly.
4406 Earlier as well as later GMP releases define S as described here.
4407 GMP releases before GMP 4.3.0 required additional space for both
4408 input and output areas. More precisely, the areas {UP, UN+1} and
4409 {VP, VN+1} were destroyed (i.e. the operands plus an extra limb
4410 past the end of each), and the areas pointed to by GP and SP should
4411 each have room for UN+1 limbs.
4413 -- Function: mp_size_t mpn_sqrtrem (mp_limb_t *R1P, mp_limb_t *R2P,
4414 const mp_limb_t *SP, mp_size_t N)
4415 Compute the square root of {SP, N} and put the result at {R1P,
4416 ceil(N/2)} and the remainder at {R2P, RETVAL}. R2P needs space for
4417 N limbs, but the return value indicates how many are produced.
4419 The most significant limb of {SP, N} must be non-zero. The areas
4420 {R1P, ceil(N/2)} and {SP, N} must be completely separate. The
4421 areas {R2P, N} and {SP, N} must be either identical or completely
4424 If the remainder is not wanted then R2P can be 'NULL', and in this
4425 case the return value is zero or non-zero according to whether the
4426 remainder would have been zero or non-zero.
4428 A return value of zero indicates a perfect square. See also
4429 'mpn_perfect_square_p'.
4431 -- Function: size_t mpn_sizeinbase (const mp_limb_t *XP, mp_size_t N,
4433 Return the size of {XP,N} measured in number of digits in the given
4434 BASE. BASE can vary from 2 to 62. Requires N > 0 and XP[N-1] > 0.
4435 The result will be either exact or 1 too big. If BASE is a power
4436 of 2, the result is always exact.
4438 -- Function: mp_size_t mpn_get_str (unsigned char *STR, int BASE,
4439 mp_limb_t *S1P, mp_size_t S1N)
4440 Convert {S1P, S1N} to a raw unsigned char array at STR in base
4441 BASE, and return the number of characters produced. There may be
4442 leading zeros in the string. The string is not in ASCII; to
4443 convert it to printable format, add the ASCII codes for '0' or 'A',
4444 depending on the base and range. BASE can vary from 2 to 256.
4446 The most significant limb of the input {S1P, S1N} must be non-zero.
4447 The input {S1P, S1N} is clobbered, except when BASE is a power of
4448 2, in which case it's unchanged.
4450 The area at STR has to have space for the largest possible number
4451 represented by a S1N long limb array, plus one extra character.
4453 -- Function: mp_size_t mpn_set_str (mp_limb_t *RP, const unsigned char
4454 *STR, size_t STRSIZE, int BASE)
4455 Convert bytes {STR,STRSIZE} in the given BASE to limbs at RP.
4457 STR[0] is the most significant input byte and STR[STRSIZE-1] is the
4458 least significant input byte. Each byte should be a value in the
4459 range 0 to BASE-1, not an ASCII character. BASE can vary from 2 to
4462 The converted value is {RP,RN} where RN is the return value. If
4463 the most significant input byte STR[0] is non-zero, then RP[RN-1]
4464 will be non-zero, else RP[RN-1] and some number of subsequent limbs
4467 The area at RP has to have space for the largest possible number
4468 with STRSIZE digits in the chosen base, plus one extra limb.
4470 The input must have at least one byte, and no overlap is permitted
4471 between {STR,STRSIZE} and the result at RP.
4473 -- Function: mp_bitcnt_t mpn_scan0 (const mp_limb_t *S1P, mp_bitcnt_t
4475 Scan S1P from bit position BIT for the next clear bit.
4477 It is required that there be a clear bit within the area at S1P at
4478 or beyond bit position BIT, so that the function has something to
4481 -- Function: mp_bitcnt_t mpn_scan1 (const mp_limb_t *S1P, mp_bitcnt_t
4483 Scan S1P from bit position BIT for the next set bit.
4485 It is required that there be a set bit within the area at S1P at or
4486 beyond bit position BIT, so that the function has something to
4489 -- Function: void mpn_random (mp_limb_t *R1P, mp_size_t R1N)
4490 -- Function: void mpn_random2 (mp_limb_t *R1P, mp_size_t R1N)
4491 Generate a random number of length R1N and store it at R1P. The
4492 most significant limb is always non-zero. 'mpn_random' generates
4493 uniformly distributed limb data, 'mpn_random2' generates long
4494 strings of zeros and ones in the binary representation.
4496 'mpn_random2' is intended for testing the correctness of the 'mpn'
4499 -- Function: mp_bitcnt_t mpn_popcount (const mp_limb_t *S1P, mp_size_t
4501 Count the number of set bits in {S1P, N}.
4503 -- Function: mp_bitcnt_t mpn_hamdist (const mp_limb_t *S1P, const
4504 mp_limb_t *S2P, mp_size_t N)
4505 Compute the hamming distance between {S1P, N} and {S2P, N}, which
4506 is the number of bit positions where the two operands have
4507 different bit values.
4509 -- Function: int mpn_perfect_square_p (const mp_limb_t *S1P, mp_size_t
4511 Return non-zero iff {S1P, N} is a perfect square. The most
4512 significant limb of the input {S1P, N} must be non-zero.
4514 -- Function: void mpn_and_n (mp_limb_t *RP, const mp_limb_t *S1P, const
4515 mp_limb_t *S2P, mp_size_t N)
4516 Perform the bitwise logical and of {S1P, N} and {S2P, N}, and write
4517 the result to {RP, N}.
4519 -- Function: void mpn_ior_n (mp_limb_t *RP, const mp_limb_t *S1P, const
4520 mp_limb_t *S2P, mp_size_t N)
4521 Perform the bitwise logical inclusive or of {S1P, N} and {S2P, N},
4522 and write the result to {RP, N}.
4524 -- Function: void mpn_xor_n (mp_limb_t *RP, const mp_limb_t *S1P, const
4525 mp_limb_t *S2P, mp_size_t N)
4526 Perform the bitwise logical exclusive or of {S1P, N} and {S2P, N},
4527 and write the result to {RP, N}.
4529 -- Function: void mpn_andn_n (mp_limb_t *RP, const mp_limb_t *S1P,
4530 const mp_limb_t *S2P, mp_size_t N)
4531 Perform the bitwise logical and of {S1P, N} and the bitwise
4532 complement of {S2P, N}, and write the result to {RP, N}.
4534 -- Function: void mpn_iorn_n (mp_limb_t *RP, const mp_limb_t *S1P,
4535 const mp_limb_t *S2P, mp_size_t N)
4536 Perform the bitwise logical inclusive or of {S1P, N} and the
4537 bitwise complement of {S2P, N}, and write the result to {RP, N}.
4539 -- Function: void mpn_nand_n (mp_limb_t *RP, const mp_limb_t *S1P,
4540 const mp_limb_t *S2P, mp_size_t N)
4541 Perform the bitwise logical and of {S1P, N} and {S2P, N}, and write
4542 the bitwise complement of the result to {RP, N}.
4544 -- Function: void mpn_nior_n (mp_limb_t *RP, const mp_limb_t *S1P,
4545 const mp_limb_t *S2P, mp_size_t N)
4546 Perform the bitwise logical inclusive or of {S1P, N} and {S2P, N},
4547 and write the bitwise complement of the result to {RP, N}.
4549 -- Function: void mpn_xnor_n (mp_limb_t *RP, const mp_limb_t *S1P,
4550 const mp_limb_t *S2P, mp_size_t N)
4551 Perform the bitwise logical exclusive or of {S1P, N} and {S2P, N},
4552 and write the bitwise complement of the result to {RP, N}.
4554 -- Function: void mpn_com (mp_limb_t *RP, const mp_limb_t *SP,
4556 Perform the bitwise complement of {SP, N}, and write the result to
4559 -- Function: void mpn_copyi (mp_limb_t *RP, const mp_limb_t *S1P,
4561 Copy from {S1P, N} to {RP, N}, increasingly.
4563 -- Function: void mpn_copyd (mp_limb_t *RP, const mp_limb_t *S1P,
4565 Copy from {S1P, N} to {RP, N}, decreasingly.
4567 -- Function: void mpn_zero (mp_limb_t *RP, mp_size_t N)
4571 8.1 Low-level functions for cryptography
4572 ========================================
4574 The functions prefixed with 'mpn_sec_' and 'mpn_cnd_' are designed to
4575 perform the exact same low-level operations and have the same cache
4576 access patterns for any two same-size arguments, assuming that function
4577 arguments are placed at the same position and that the machine state is
4578 identical upon function entry. These functions are intended for
4579 cryptographic purposes, where resilience to side-channel attacks is
4582 These functions are less efficient than their "leaky" counterparts;
4583 their performance for operands of the sizes typically used for
4584 cryptographic applications is between 15% and 100% worse. For larger
4585 operands, these functions might be inadequate, since they rely on
4586 asymptotically elementary algorithms.
4588 These functions do not make any explicit allocations. Those of these
4589 functions that need scratch space accept a scratch space operand. This
4590 convention allows callers to keep sensitive data in designated memory
4591 areas. Note however that compilers may choose to spill scalar values
4592 used within these functions to their stack frame and that such scalars
4593 may contain sensitive data.
4595 In addition to these specially crafted functions, the following 'mpn'
4596 functions are naturally side-channel resistant: 'mpn_add_n',
4597 'mpn_sub_n', 'mpn_lshift', 'mpn_rshift', 'mpn_zero', 'mpn_copyi',
4598 'mpn_copyd', 'mpn_com', and the logical function ('mpn_and_n', etc).
4600 There are some exceptions from the side-channel resilience: (1) Some
4601 assembly implementations of 'mpn_lshift' identify shift-by-one as a
4602 special case. This is a problem iff the shift count is a function of
4603 sensitive data. (2) Alpha ev6 and Pentium4 using 64-bit limbs have
4604 leaky 'mpn_add_n' and 'mpn_sub_n'. (3) Alpha ev6 has a leaky
4605 'mpn_mul_1' which also makes 'mpn_sec_mul' on those systems unsafe.
4607 -- Function: mp_limb_t mpn_cnd_add_n (mp_limb_t CND, mp_limb_t *RP,
4608 const mp_limb_t *S1P, const mp_limb_t *S2P, mp_size_t N)
4609 -- Function: mp_limb_t mpn_cnd_sub_n (mp_limb_t CND, mp_limb_t *RP,
4610 const mp_limb_t *S1P, const mp_limb_t *S2P, mp_size_t N)
4611 These functions do conditional addition and subtraction. If CND is
4612 non-zero, they produce the same result as a regular 'mpn_add_n' or
4613 'mpn_sub_n', and if CND is zero, they copy {S1P,N} to the result
4614 area and return zero. The functions are designed to have timing
4615 and memory access patterns depending only on size and location of
4616 the data areas, but independent of the condition CND. Like for
4617 'mpn_add_n' and 'mpn_sub_n', on most machines, the timing will also
4618 be independent of the actual limb values.
4620 -- Function: mp_limb_t mpn_sec_add_1 (mp_limb_t *RP, const mp_limb_t
4621 *AP, mp_size_t N, mp_limb_t B, mp_limb_t *TP)
4622 -- Function: mp_limb_t mpn_sec_sub_1 (mp_limb_t *RP, const mp_limb_t
4623 *AP, mp_size_t N, mp_limb_t B, mp_limb_t *TP)
4624 Set R to A + B or A - B, respectively, where R = {RP,N}, A =
4625 {AP,N}, and B is a single limb. Returns carry.
4627 These functions take O(N) time, unlike the leaky functions
4628 'mpn_add_1' which are O(1) on average. They require scratch space
4629 of 'mpn_sec_add_1_itch(N)' and 'mpn_sec_sub_1_itch(N)' limbs,
4630 respectively, to be passed in the TP parameter. The scratch space
4631 requirements are guaranteed to be at most N limbs, and increase
4632 monotonously in the operand size.
4634 -- Function: void mpn_cnd_swap (mp_limb_t CND, volatile mp_limb_t *AP,
4635 volatile mp_limb_t *BP, mp_size_t N)
4636 If CND is non-zero, swaps the contents of the areas {AP,N} and
4637 {BP,N}. Otherwise, the areas are left unmodified. Implemented
4638 using logical operations on the limbs, with the same memory
4639 accesses independent of the value of CND.
4641 -- Function: void mpn_sec_mul (mp_limb_t *RP, const mp_limb_t *AP,
4642 mp_size_t AN, const mp_limb_t *BP, mp_size_t BN, mp_limb_t
4644 -- Function: mp_size_t mpn_sec_mul_itch (mp_size_t AN, mp_size_t BN)
4645 Set R to A * B, where A = {AP,AN}, B = {BP,BN}, and R = {RP,AN+BN}.
4647 It is required that AN >= BN > 0.
4649 No overlapping between R and the input operands is allowed. For A
4650 = B, use 'mpn_sec_sqr' for optimal performance.
4652 This function requires scratch space of 'mpn_sec_mul_itch(AN, BN)'
4653 limbs to be passed in the TP parameter. The scratch space
4654 requirements are guaranteed to increase monotonously in the operand
4657 -- Function: void mpn_sec_sqr (mp_limb_t *RP, const mp_limb_t *AP,
4658 mp_size_t AN, mp_limb_t *TP)
4659 -- Function: mp_size_t mpn_sec_sqr_itch (mp_size_t AN)
4660 Set R to A^2, where A = {AP,AN}, and R = {RP,2AN}.
4662 It is required that AN > 0.
4664 No overlapping between R and the input operands is allowed.
4666 This function requires scratch space of 'mpn_sec_sqr_itch(AN)'
4667 limbs to be passed in the TP parameter. The scratch space
4668 requirements are guaranteed to increase monotonously in the operand
4671 -- Function: void mpn_sec_powm (mp_limb_t *RP, const mp_limb_t *BP,
4672 mp_size_t BN, const mp_limb_t *EP, mp_bitcnt_t ENB, const
4673 mp_limb_t *MP, mp_size_t N, mp_limb_t *TP)
4674 -- Function: mp_size_t mpn_sec_powm_itch (mp_size_t BN, mp_bitcnt_t
4676 Set R to (B raised to E) modulo M, where R = {RP,N}, M = {MP,N},
4677 and E = {EP,ceil(ENB / 'GMP\_NUMB\_BITS')}.
4679 It is required that B > 0, that M > 0 is odd, and that E < 2^ENB,
4682 No overlapping between R and the input operands is allowed.
4684 This function requires scratch space of 'mpn_sec_powm_itch(BN, ENB,
4685 N)' limbs to be passed in the TP parameter. The scratch space
4686 requirements are guaranteed to increase monotonously in the operand
4689 -- Function: void mpn_sec_tabselect (mp_limb_t *RP, const mp_limb_t
4690 *TAB, mp_size_t N, mp_size_t NENTS, mp_size_t WHICH)
4691 Select entry WHICH from table TAB, which has NENTS entries, each N
4692 limbs. Store the selected entry at RP.
4694 This function reads the entire table to avoid side-channel
4697 -- Function: mp_limb_t mpn_sec_div_qr (mp_limb_t *QP, mp_limb_t *NP,
4698 mp_size_t NN, const mp_limb_t *DP, mp_size_t DN, mp_limb_t
4700 -- Function: mp_size_t mpn_sec_div_qr_itch (mp_size_t NN, mp_size_t DN)
4702 Set Q to the truncated quotient N / D and R to N modulo D, where N
4703 = {NP,NN}, D = {DP,DN}, Q's most significant limb is the function
4704 return value and the remaining limbs are {QP,NN-DN}, and R =
4707 It is required that NN >= DN >= 1, and that DP[DN-1] != 0. This
4708 does not imply that N >= D since N might be zero-padded.
4710 Note the overlapping between N and R. No other operand overlapping
4711 is allowed. The entire space occupied by N is overwritten.
4713 This function requires scratch space of 'mpn_sec_div_qr_itch(NN,
4714 DN)' limbs to be passed in the TP parameter.
4716 -- Function: void mpn_sec_div_r (mp_limb_t *NP, mp_size_t NN, const
4717 mp_limb_t *DP, mp_size_t DN, mp_limb_t *TP)
4718 -- Function: mp_size_t mpn_sec_div_r_itch (mp_size_t NN, mp_size_t DN)
4720 Set R to N modulo D, where N = {NP,NN}, D = {DP,DN}, and R =
4723 It is required that NN >= DN >= 1, and that DP[DN-1] != 0. This
4724 does not imply that N >= D since N might be zero-padded.
4726 Note the overlapping between N and R. No other operand overlapping
4727 is allowed. The entire space occupied by N is overwritten.
4729 This function requires scratch space of 'mpn_sec_div_r_itch(NN,
4730 DN)' limbs to be passed in the TP parameter.
4732 -- Function: int mpn_sec_invert (mp_limb_t *RP, mp_limb_t *AP, const
4733 mp_limb_t *MP, mp_size_t N, mp_bitcnt_t NBCNT, mp_limb_t *TP)
4734 -- Function: mp_size_t mpn_sec_invert_itch (mp_size_t N)
4735 Set R to the inverse of A modulo M, where R = {RP,N}, A = {AP,N},
4736 and M = {MP,N}. *This function's interface is preliminary.*
4738 If an inverse exists, return 1, otherwise return 0 and leave R
4739 undefined. In either case, the input A is destroyed.
4741 It is required that M is odd, and that NBCNT >= ceil(\log(A+1)) +
4742 ceil(\log(M+1)). A safe choice is NBCNT = 2 * N * GMP_NUMB_BITS,
4743 but a smaller value might improve performance if M or A are known
4744 to have leading zero bits.
4746 This function requires scratch space of 'mpn_sec_invert_itch(N)'
4747 limbs to be passed in the TP parameter.
4753 *Everything in this section is highly experimental and may disappear or
4754 be subject to incompatible changes in a future version of GMP.*
4756 Nails are an experimental feature whereby a few bits are left unused
4757 at the top of each 'mp_limb_t'. This can significantly improve carry
4758 handling on some processors.
4760 All the 'mpn' functions accepting limb data will expect the nail bits
4761 to be zero on entry, and will return data with the nails similarly all
4762 zero. This applies both to limb vectors and to single limb arguments.
4764 Nails can be enabled by configuring with '--enable-nails'. By
4765 default the number of bits will be chosen according to what suits the
4766 host processor, but a particular number can be selected with
4769 At the mpn level, a nail build is neither source nor binary
4770 compatible with a non-nail build, strictly speaking. But programs
4771 acting on limbs only through the mpn functions are likely to work
4772 equally well with either build, and judicious use of the definitions
4773 below should make any program compatible with either build, at the
4776 For the higher level routines, meaning 'mpz' etc, a nail build should
4777 be fully source and binary compatible with a non-nail build.
4779 -- Macro: GMP_NAIL_BITS
4780 -- Macro: GMP_NUMB_BITS
4781 -- Macro: GMP_LIMB_BITS
4782 'GMP_NAIL_BITS' is the number of nail bits, or 0 when nails are not
4783 in use. 'GMP_NUMB_BITS' is the number of data bits in a limb.
4784 'GMP_LIMB_BITS' is the total number of bits in an 'mp_limb_t'. In
4787 GMP_LIMB_BITS == GMP_NAIL_BITS + GMP_NUMB_BITS
4789 -- Macro: GMP_NAIL_MASK
4790 -- Macro: GMP_NUMB_MASK
4791 Bit masks for the nail and number parts of a limb. 'GMP_NAIL_MASK'
4792 is 0 when nails are not in use.
4794 'GMP_NAIL_MASK' is not often needed, since the nail part can be
4795 obtained with 'x >> GMP_NUMB_BITS', and that means one less large
4796 constant, which can help various RISC chips.
4798 -- Macro: GMP_NUMB_MAX
4799 The maximum value that can be stored in the number part of a limb.
4800 This is the same as 'GMP_NUMB_MASK', but can be used for clarity
4801 when doing comparisons rather than bit-wise operations.
4803 The term "nails" comes from finger or toe nails, which are at the
4804 ends of a limb (arm or leg). "numb" is short for number, but is also
4805 how the developers felt after trying for a long time to come up with
4806 sensible names for these things.
4808 In the future (the distant future most likely) a non-zero nail might
4809 be permitted, giving non-unique representations for numbers in a limb
4810 vector. This would help vector processors since carries would only ever
4811 need to propagate one or two limbs.
4814 File: gmp.info, Node: Random Number Functions, Next: Formatted Output, Prev: Low-level Functions, Up: Top
4816 9 Random Number Functions
4817 *************************
4819 Sequences of pseudo-random numbers in GMP are generated using a variable
4820 of type 'gmp_randstate_t', which holds an algorithm selection and a
4821 current state. Such a variable must be initialized by a call to one of
4822 the 'gmp_randinit' functions, and can be seeded with one of the
4823 'gmp_randseed' functions.
4825 The functions actually generating random numbers are described in
4826 *note Integer Random Numbers::, and *note Miscellaneous Float
4829 The older style random number functions don't accept a
4830 'gmp_randstate_t' parameter but instead share a global variable of that
4831 type. They use a default algorithm and are currently not seeded (though
4832 perhaps that will change in the future). The new functions accepting a
4833 'gmp_randstate_t' are recommended for applications that care about
4838 * Random State Initialization::
4839 * Random State Seeding::
4840 * Random State Miscellaneous::
4843 File: gmp.info, Node: Random State Initialization, Next: Random State Seeding, Prev: Random Number Functions, Up: Random Number Functions
4845 9.1 Random State Initialization
4846 ===============================
4848 -- Function: void gmp_randinit_default (gmp_randstate_t STATE)
4849 Initialize STATE with a default algorithm. This will be a
4850 compromise between speed and randomness, and is recommended for
4851 applications with no special requirements. Currently this is
4854 -- Function: void gmp_randinit_mt (gmp_randstate_t STATE)
4855 Initialize STATE for a Mersenne Twister algorithm. This algorithm
4856 is fast and has good randomness properties.
4858 -- Function: void gmp_randinit_lc_2exp (gmp_randstate_t STATE, const
4859 mpz_t A, unsigned long C, mp_bitcnt_t M2EXP)
4860 Initialize STATE with a linear congruential algorithm X = (A*X + C)
4863 The low bits of X in this algorithm are not very random. The least
4864 significant bit will have a period no more than 2, and the second
4865 bit no more than 4, etc. For this reason only the high half of
4866 each X is actually used.
4868 When a random number of more than M2EXP/2 bits is to be generated,
4869 multiple iterations of the recurrence are used and the results
4872 -- Function: int gmp_randinit_lc_2exp_size (gmp_randstate_t STATE,
4874 Initialize STATE for a linear congruential algorithm as per
4875 'gmp_randinit_lc_2exp'. A, C and M2EXP are selected from a table,
4876 chosen so that SIZE bits (or more) of each X will be used, i.e.
4879 If successful the return value is non-zero. If SIZE is bigger than
4880 the table data provides then the return value is zero. The maximum
4881 SIZE currently supported is 128.
4883 -- Function: void gmp_randinit_set (gmp_randstate_t ROP,
4885 Initialize ROP with a copy of the algorithm and state from OP.
4887 -- Function: void gmp_randinit (gmp_randstate_t STATE,
4888 gmp_randalg_t ALG, ...)
4889 *This function is obsolete.*
4891 Initialize STATE with an algorithm selected by ALG. The only
4892 choice is 'GMP_RAND_ALG_LC', which is 'gmp_randinit_lc_2exp_size'
4893 described above. A third parameter of type 'unsigned long' is
4894 required, this is the SIZE for that function.
4895 'GMP_RAND_ALG_DEFAULT' or 0 are the same as 'GMP_RAND_ALG_LC'.
4897 'gmp_randinit' sets bits in the global variable 'gmp_errno' to
4898 indicate an error. 'GMP_ERROR_UNSUPPORTED_ARGUMENT' if ALG is
4899 unsupported, or 'GMP_ERROR_INVALID_ARGUMENT' if the SIZE parameter
4900 is too big. It may be noted this error reporting is not thread
4901 safe (a good reason to use 'gmp_randinit_lc_2exp_size' instead).
4903 -- Function: void gmp_randclear (gmp_randstate_t STATE)
4904 Free all memory occupied by STATE.
4907 File: gmp.info, Node: Random State Seeding, Next: Random State Miscellaneous, Prev: Random State Initialization, Up: Random Number Functions
4909 9.2 Random State Seeding
4910 ========================
4912 -- Function: void gmp_randseed (gmp_randstate_t STATE, const mpz_t
4914 -- Function: void gmp_randseed_ui (gmp_randstate_t STATE,
4915 unsigned long int SEED)
4916 Set an initial seed value into STATE.
4918 The size of a seed determines how many different sequences of
4919 random numbers that it's possible to generate. The "quality" of
4920 the seed is the randomness of a given seed compared to the previous
4921 seed used, and this affects the randomness of separate number
4922 sequences. The method for choosing a seed is critical if the
4923 generated numbers are to be used for important applications, such
4924 as generating cryptographic keys.
4926 Traditionally the system time has been used to seed, but care needs
4927 to be taken with this. If an application seeds often and the
4928 resolution of the system clock is low, then the same sequence of
4929 numbers might be repeated. Also, the system time is quite easy to
4930 guess, so if unpredictability is required then it should definitely
4931 not be the only source for the seed value. On some systems there's
4932 a special device '/dev/random' which provides random data better
4933 suited for use as a seed.
4936 File: gmp.info, Node: Random State Miscellaneous, Prev: Random State Seeding, Up: Random Number Functions
4938 9.3 Random State Miscellaneous
4939 ==============================
4941 -- Function: unsigned long gmp_urandomb_ui (gmp_randstate_t STATE,
4943 Return a uniformly distributed random number of N bits, i.e. in the
4944 range 0 to 2^N-1 inclusive. N must be less than or equal to the
4945 number of bits in an 'unsigned long'.
4947 -- Function: unsigned long gmp_urandomm_ui (gmp_randstate_t STATE,
4949 Return a uniformly distributed random number in the range 0 to N-1,
4953 File: gmp.info, Node: Formatted Output, Next: Formatted Input, Prev: Random Number Functions, Up: Top
4960 * Formatted Output Strings::
4961 * Formatted Output Functions::
4962 * C++ Formatted Output::
4965 File: gmp.info, Node: Formatted Output Strings, Next: Formatted Output Functions, Prev: Formatted Output, Up: Formatted Output
4970 'gmp_printf' and friends accept format strings similar to the standard C
4971 'printf' (*note Formatted Output: (libc)Formatted Output.). A format
4972 specification is of the form
4974 % [flags] [width] [.[precision]] [type] conv
4976 GMP adds types 'Z', 'Q' and 'F' for 'mpz_t', 'mpq_t' and 'mpf_t'
4977 respectively, 'M' for 'mp_limb_t', and 'N' for an 'mp_limb_t' array.
4978 'Z', 'Q', 'M' and 'N' behave like integers. 'Q' will print a '/' and a
4979 denominator, if needed. 'F' behaves like a float. For example,
4982 gmp_printf ("%s is an mpz %Zd\n", "here", z);
4985 gmp_printf ("a hex rational: %#40Qx\n", q);
4989 gmp_printf ("fixed point mpf %.*Ff with %d digits\n", n, f, n);
4992 gmp_printf ("limb %Mu\n", l);
4994 const mp_limb_t *ptr;
4996 gmp_printf ("limb array %Nx\n", ptr, size);
4998 For 'N' the limbs are expected least significant first, as per the
4999 'mpn' functions (*note Low-level Functions::). A negative size can be
5000 given to print the value as a negative.
5002 All the standard C 'printf' types behave the same as the C library
5003 'printf', and can be freely intermixed with the GMP extensions. In the
5004 current implementation the standard parts of the format string are
5005 simply handed to 'printf' and only the GMP extensions handled directly.
5007 The flags accepted are as follows. GLIBC style ' is only for the
5008 standard C types (not the GMP types), and only if the C library supports
5011 0 pad with zeros (rather than spaces)
5012 # show the base with '0x', '0X' or '0'
5013 + always show a sign
5014 (space) show a space or a '-' sign
5015 ' group digits, GLIBC style (not GMP
5018 The optional width and precision can be given as a number within the
5019 format string, or as a '*' to take an extra parameter of type 'int', the
5020 same as the standard 'printf'.
5022 The standard types accepted are as follows. 'h' and 'l' are
5023 portable, the rest will depend on the compiler (or include files) for
5024 the type and the C library for the output.
5028 j intmax_t or uintmax_t
5032 q quad_t or u_quad_t
5038 F mpf_t, float conversions
5039 Q mpq_t, integer conversions
5040 M mp_limb_t, integer conversions
5041 N mp_limb_t array, integer conversions
5042 Z mpz_t, integer conversions
5044 The conversions accepted are as follows. 'a' and 'A' are always
5045 supported for 'mpf_t' but depend on the C library for standard C float
5046 types. 'm' and 'p' depend on the C library.
5048 a A hex floats, C99 style
5051 e E scientific format float
5054 g G fixed or scientific float
5055 m 'strerror' string, GLIBC style
5056 n store characters written so far
5063 'o', 'x' and 'X' are unsigned for the standard C types, but for types
5064 'Z', 'Q' and 'N' they are signed. 'u' is not meaningful for 'Z', 'Q'
5067 'M' is a proxy for the C library 'l' or 'L', according to the size of
5068 'mp_limb_t'. Unsigned conversions will be usual, but a signed
5069 conversion can be used and will interpret the value as a twos complement
5072 'n' can be used with any type, even the GMP types.
5074 Other types or conversions that might be accepted by the C library
5075 'printf' cannot be used through 'gmp_printf', this includes for instance
5076 extensions registered with GLIBC 'register_printf_function'. Also
5077 currently there's no support for POSIX '$' style numbered arguments
5078 (perhaps this will be added in the future).
5080 The precision field has its usual meaning for integer 'Z' and float
5081 'F' types, but is currently undefined for 'Q' and should not be used
5084 'mpf_t' conversions only ever generate as many digits as can be
5085 accurately represented by the operand, the same as 'mpf_get_str' does.
5086 Zeros will be used if necessary to pad to the requested precision. This
5087 happens even for an 'f' conversion of an 'mpf_t' which is an integer,
5088 for instance 2^1024 in an 'mpf_t' of 128 bits precision will only
5089 produce about 40 digits, then pad with zeros to the decimal point. An
5090 empty precision field like '%.Fe' or '%.Ff' can be used to specifically
5091 request just the significant digits. Without any dot and thus no
5092 precision field, a precision value of 6 will be used. Note that these
5093 rules mean that '%Ff', '%.Ff', and '%.0Ff' will all be different.
5095 The decimal point character (or string) is taken from the current
5096 locale settings on systems which provide 'localeconv' (*note Locales and
5097 Internationalization: (libc)Locales.). The C library will normally do
5098 the same for standard float output.
5100 The format string is only interpreted as plain 'char's, multibyte
5101 characters are not recognised. Perhaps this will change in the future.
5104 File: gmp.info, Node: Formatted Output Functions, Next: C++ Formatted Output, Prev: Formatted Output Strings, Up: Formatted Output
5109 Each of the following functions is similar to the corresponding C
5110 library function. The basic 'printf' forms take a variable argument
5111 list. The 'vprintf' forms take an argument pointer, see *note Variadic
5112 Functions: (libc)Variadic Functions, or 'man 3 va_start'.
5114 It should be emphasised that if a format string is invalid, or the
5115 arguments don't match what the format specifies, then the behaviour of
5116 any of these functions will be unpredictable. GCC format string
5117 checking is not available, since it doesn't recognise the GMP
5120 The file based functions 'gmp_printf' and 'gmp_fprintf' will return
5121 -1 to indicate a write error. Output is not "atomic", so partial output
5122 may be produced if a write error occurs. All the functions can return
5123 -1 if the C library 'printf' variant in use returns -1, but this
5124 shouldn't normally occur.
5126 -- Function: int gmp_printf (const char *FMT, ...)
5127 -- Function: int gmp_vprintf (const char *FMT, va_list AP)
5128 Print to the standard output 'stdout'. Return the number of
5129 characters written, or -1 if an error occurred.
5131 -- Function: int gmp_fprintf (FILE *FP, const char *FMT, ...)
5132 -- Function: int gmp_vfprintf (FILE *FP, const char *FMT, va_list AP)
5133 Print to the stream FP. Return the number of characters written,
5134 or -1 if an error occurred.
5136 -- Function: int gmp_sprintf (char *BUF, const char *FMT, ...)
5137 -- Function: int gmp_vsprintf (char *BUF, const char *FMT, va_list AP)
5138 Form a null-terminated string in BUF. Return the number of
5139 characters written, excluding the terminating null.
5141 No overlap is permitted between the space at BUF and the string
5144 These functions are not recommended, since there's no protection
5145 against exceeding the space available at BUF.
5147 -- Function: int gmp_snprintf (char *BUF, size_t SIZE, const char *FMT,
5149 -- Function: int gmp_vsnprintf (char *BUF, size_t SIZE, const char
5151 Form a null-terminated string in BUF. No more than SIZE bytes will
5152 be written. To get the full output, SIZE must be enough for the
5153 string and null-terminator.
5155 The return value is the total number of characters which ought to
5156 have been produced, excluding the terminating null. If RETVAL >=
5157 SIZE then the actual output has been truncated to the first SIZE-1
5158 characters, and a null appended.
5160 No overlap is permitted between the region {BUF,SIZE} and the FMT
5163 Notice the return value is in ISO C99 'snprintf' style. This is so
5164 even if the C library 'vsnprintf' is the older GLIBC 2.0.x style.
5166 -- Function: int gmp_asprintf (char **PP, const char *FMT, ...)
5167 -- Function: int gmp_vasprintf (char **PP, const char *FMT, va_list AP)
5168 Form a null-terminated string in a block of memory obtained from
5169 the current memory allocation function (*note Custom Allocation::).
5170 The block will be the size of the string and null-terminator. The
5171 address of the block in stored to *PP. The return value is the
5172 number of characters produced, excluding the null-terminator.
5174 Unlike the C library 'asprintf', 'gmp_asprintf' doesn't return -1
5175 if there's no more memory available, it lets the current allocation
5176 function handle that.
5178 -- Function: int gmp_obstack_printf (struct obstack *OB, const char
5180 -- Function: int gmp_obstack_vprintf (struct obstack *OB, const char
5182 Append to the current object in OB. The return value is the number
5183 of characters written. A null-terminator is not written.
5185 FMT cannot be within the current object in OB, since that object
5186 might move as it grows.
5188 These functions are available only when the C library provides the
5189 obstack feature, which probably means only on GNU systems, see
5190 *note Obstacks: (libc)Obstacks.
5193 File: gmp.info, Node: C++ Formatted Output, Prev: Formatted Output Functions, Up: Formatted Output
5195 10.3 C++ Formatted Output
5196 =========================
5198 The following functions are provided in 'libgmpxx' (*note Headers and
5199 Libraries::), which is built if C++ support is enabled (*note Build
5200 Options::). Prototypes are available from '<gmp.h>'.
5202 -- Function: ostream& operator<< (ostream& STREAM, const mpz_t OP)
5203 Print OP to STREAM, using its 'ios' formatting settings.
5204 'ios::width' is reset to 0 after output, the same as the standard
5205 'ostream operator<<' routines do.
5207 In hex or octal, OP is printed as a signed number, the same as for
5208 decimal. This is unlike the standard 'operator<<' routines on
5209 'int' etc, which instead give twos complement.
5211 -- Function: ostream& operator<< (ostream& STREAM, const mpq_t OP)
5212 Print OP to STREAM, using its 'ios' formatting settings.
5213 'ios::width' is reset to 0 after output, the same as the standard
5214 'ostream operator<<' routines do.
5216 Output will be a fraction like '5/9', or if the denominator is 1
5217 then just a plain integer like '123'.
5219 In hex or octal, OP is printed as a signed value, the same as for
5220 decimal. If 'ios::showbase' is set then a base indicator is shown
5221 on both the numerator and denominator (if the denominator is
5224 -- Function: ostream& operator<< (ostream& STREAM, const mpf_t OP)
5225 Print OP to STREAM, using its 'ios' formatting settings.
5226 'ios::width' is reset to 0 after output, the same as the standard
5227 'ostream operator<<' routines do.
5229 The decimal point follows the standard library float 'operator<<',
5230 which on recent systems means the 'std::locale' imbued on STREAM.
5232 Hex and octal are supported, unlike the standard 'operator<<' on
5233 'double'. The mantissa will be in hex or octal, the exponent will
5234 be in decimal. For hex the exponent delimiter is an '@'. This is
5235 as per 'mpf_out_str'.
5237 'ios::showbase' is supported, and will put a base on the mantissa,
5238 for example hex '0x1.8' or '0x0.8', or octal '01.4' or '00.4'.
5239 This last form is slightly strange, but at least differentiates
5240 itself from decimal.
5242 These operators mean that GMP types can be printed in the usual C++
5248 cout << "iteration " << n << " value " << z << "\n";
5250 But note that 'ostream' output (and 'istream' input, *note C++
5251 Formatted Input::) is the only overloading available for the GMP types
5252 and that for instance using '+' with an 'mpz_t' will have unpredictable
5253 results. For classes with overloading, see *note C++ Class Interface::.
5256 File: gmp.info, Node: Formatted Input, Next: C++ Class Interface, Prev: Formatted Output, Up: Top
5263 * Formatted Input Strings::
5264 * Formatted Input Functions::
5265 * C++ Formatted Input::
5268 File: gmp.info, Node: Formatted Input Strings, Next: Formatted Input Functions, Prev: Formatted Input, Up: Formatted Input
5270 11.1 Formatted Input Strings
5271 ============================
5273 'gmp_scanf' and friends accept format strings similar to the standard C
5274 'scanf' (*note Formatted Input: (libc)Formatted Input.). A format
5275 specification is of the form
5277 % [flags] [width] [type] conv
5279 GMP adds types 'Z', 'Q' and 'F' for 'mpz_t', 'mpq_t' and 'mpf_t'
5280 respectively. 'Z' and 'Q' behave like integers. 'Q' will read a '/'
5281 and a denominator, if present. 'F' behaves like a float.
5283 GMP variables don't require an '&' when passed to 'gmp_scanf', since
5284 they're already "call-by-reference". For example,
5286 /* to read say "a(5) = 1234" */
5289 gmp_scanf ("a(%d) = %Zd\n", &n, z);
5292 gmp_sscanf ("0377 + 0x10/0x11", "%Qi + %Qi", q1, q2);
5294 /* to read say "topleft (1.55,-2.66)" */
5297 gmp_scanf ("%31s (%Ff,%Ff)", buf, x, y);
5299 All the standard C 'scanf' types behave the same as in the C library
5300 'scanf', and can be freely intermixed with the GMP extensions. In the
5301 current implementation the standard parts of the format string are
5302 simply handed to 'scanf' and only the GMP extensions handled directly.
5304 The flags accepted are as follows. 'a' and ''' will depend on
5305 support from the C library, and ''' cannot be used with GMP types.
5307 * read but don't store
5308 a allocate a buffer (string conversions)
5309 ' grouped digits, GLIBC style (not GMP
5312 The standard types accepted are as follows. 'h' and 'l' are
5313 portable, the rest will depend on the compiler (or include files) for
5314 the type and the C library for the input.
5318 j intmax_t or uintmax_t
5319 l long int, double or wchar_t
5322 q quad_t or u_quad_t
5328 F mpf_t, float conversions
5329 Q mpq_t, integer conversions
5330 Z mpz_t, integer conversions
5332 The conversions accepted are as follows. 'p' and '[' will depend on
5333 support from the C library, the rest are standard.
5335 c character or characters
5339 i integer with base indicator
5340 n characters read so far
5343 s string of non-whitespace characters
5346 [ string of characters in a set
5348 'e', 'E', 'f', 'g' and 'G' are identical, they all read either fixed
5349 point or scientific format, and either upper or lower case 'e' for the
5350 exponent in scientific format.
5352 C99 style hex float format ('printf %a', *note Formatted Output
5353 Strings::) is always accepted for 'mpf_t', but for the standard float
5354 types it will depend on the C library.
5356 'x' and 'X' are identical, both accept both upper and lower case
5359 'o', 'u', 'x' and 'X' all read positive or negative values. For the
5360 standard C types these are described as "unsigned" conversions, but that
5361 merely affects certain overflow handling, negatives are still allowed
5362 (per 'strtoul', *note Parsing of Integers: (libc)Parsing of Integers.).
5363 For GMP types there are no overflows, so 'd' and 'u' are identical.
5365 'Q' type reads the numerator and (optional) denominator as given. If
5366 the value might not be in canonical form then 'mpq_canonicalize' must be
5367 called before using it in any calculations (*note Rational Number
5370 'Qi' will read a base specification separately for the numerator and
5371 denominator. For example '0x10/11' would be 16/11, whereas '0x10/0x11'
5374 'n' can be used with any of the types above, even the GMP types. '*'
5375 to suppress assignment is allowed, though in that case it would do
5378 Other conversions or types that might be accepted by the C library
5379 'scanf' cannot be used through 'gmp_scanf'.
5381 Whitespace is read and discarded before a field, except for 'c' and
5384 For float conversions, the decimal point character (or string)
5385 expected is taken from the current locale settings on systems which
5386 provide 'localeconv' (*note Locales and Internationalization:
5387 (libc)Locales.). The C library will normally do the same for standard
5390 The format string is only interpreted as plain 'char's, multibyte
5391 characters are not recognised. Perhaps this will change in the future.
5394 File: gmp.info, Node: Formatted Input Functions, Next: C++ Formatted Input, Prev: Formatted Input Strings, Up: Formatted Input
5396 11.2 Formatted Input Functions
5397 ==============================
5399 Each of the following functions is similar to the corresponding C
5400 library function. The plain 'scanf' forms take a variable argument
5401 list. The 'vscanf' forms take an argument pointer, see *note Variadic
5402 Functions: (libc)Variadic Functions, or 'man 3 va_start'.
5404 It should be emphasised that if a format string is invalid, or the
5405 arguments don't match what the format specifies, then the behaviour of
5406 any of these functions will be unpredictable. GCC format string
5407 checking is not available, since it doesn't recognise the GMP
5410 No overlap is permitted between the FMT string and any of the results
5413 -- Function: int gmp_scanf (const char *FMT, ...)
5414 -- Function: int gmp_vscanf (const char *FMT, va_list AP)
5415 Read from the standard input 'stdin'.
5417 -- Function: int gmp_fscanf (FILE *FP, const char *FMT, ...)
5418 -- Function: int gmp_vfscanf (FILE *FP, const char *FMT, va_list AP)
5419 Read from the stream FP.
5421 -- Function: int gmp_sscanf (const char *S, const char *FMT, ...)
5422 -- Function: int gmp_vsscanf (const char *S, const char *FMT, va_list
5424 Read from a null-terminated string S.
5426 The return value from each of these functions is the same as the
5427 standard C99 'scanf', namely the number of fields successfully parsed
5428 and stored. '%n' fields and fields read but suppressed by '*' don't
5429 count towards the return value.
5431 If end of input (or a file error) is reached before a character for a
5432 field or a literal, and if no previous non-suppressed fields have
5433 matched, then the return value is 'EOF' instead of 0. A whitespace
5434 character in the format string is only an optional match and doesn't
5435 induce an 'EOF' in this fashion. Leading whitespace read and discarded
5436 for a field don't count as characters for that field.
5438 For the GMP types, input parsing follows C99 rules, namely one
5439 character of lookahead is used and characters are read while they
5440 continue to meet the format requirements. If this doesn't provide a
5441 complete number then the function terminates, with that field not stored
5442 nor counted towards the return value. For instance with 'mpf_t' an
5443 input '1.23e-XYZ' would be read up to the 'X' and that character pushed
5444 back since it's not a digit. The string '1.23e-' would then be
5445 considered invalid since an 'e' must be followed by at least one digit.
5447 For the standard C types, in the current implementation GMP calls the
5448 C library 'scanf' functions, which might have looser rules about what
5449 constitutes a valid input.
5451 Note that 'gmp_sscanf' is the same as 'gmp_fscanf' and only does one
5452 character of lookahead when parsing. Although clearly it could look at
5453 its entire input, it is deliberately made identical to 'gmp_fscanf', the
5454 same way C99 'sscanf' is the same as 'fscanf'.
5457 File: gmp.info, Node: C++ Formatted Input, Prev: Formatted Input Functions, Up: Formatted Input
5459 11.3 C++ Formatted Input
5460 ========================
5462 The following functions are provided in 'libgmpxx' (*note Headers and
5463 Libraries::), which is built only if C++ support is enabled (*note Build
5464 Options::). Prototypes are available from '<gmp.h>'.
5466 -- Function: istream& operator>> (istream& STREAM, mpz_t ROP)
5467 Read ROP from STREAM, using its 'ios' formatting settings.
5469 -- Function: istream& operator>> (istream& STREAM, mpq_t ROP)
5470 An integer like '123' will be read, or a fraction like '5/9'. No
5471 whitespace is allowed around the '/'. If the fraction is not in
5472 canonical form then 'mpq_canonicalize' must be called (*note
5473 Rational Number Functions::) before operating on it.
5475 As per integer input, an '0' or '0x' base indicator is read when
5476 none of 'ios::dec', 'ios::oct' or 'ios::hex' are set. This is done
5477 separately for numerator and denominator, so that for instance
5478 '0x10/11' is 16/11 and '0x10/0x11' is 16/17.
5480 -- Function: istream& operator>> (istream& STREAM, mpf_t ROP)
5481 Read ROP from STREAM, using its 'ios' formatting settings.
5483 Hex or octal floats are not supported, but might be in the future,
5484 or perhaps it's best to accept only what the standard float
5487 Note that digit grouping specified by the 'istream' locale is
5488 currently not accepted. Perhaps this will change in the future.
5491 These operators mean that GMP types can be read in the usual C++ way,
5498 But note that 'istream' input (and 'ostream' output, *note C++
5499 Formatted Output::) is the only overloading available for the GMP types
5500 and that for instance using '+' with an 'mpz_t' will have unpredictable
5501 results. For classes with overloading, see *note C++ Class Interface::.
5504 File: gmp.info, Node: C++ Class Interface, Next: Custom Allocation, Prev: Formatted Input, Up: Top
5506 12 C++ Class Interface
5507 **********************
5509 This chapter describes the C++ class based interface to GMP.
5511 All GMP C language types and functions can be used in C++ programs,
5512 since 'gmp.h' has 'extern "C"' qualifiers, but the class interface
5513 offers overloaded functions and operators which may be more convenient.
5515 Due to the implementation of this interface, a reasonably recent C++
5516 compiler is required, one supporting namespaces, partial specialization
5517 of templates and member templates.
5519 *Everything described in this chapter is to be considered preliminary
5520 and might be subject to incompatible changes if some unforeseen
5521 difficulty reveals itself.*
5525 * C++ Interface General::
5526 * C++ Interface Integers::
5527 * C++ Interface Rationals::
5528 * C++ Interface Floats::
5529 * C++ Interface Random Numbers::
5530 * C++ Interface Limitations::
5533 File: gmp.info, Node: C++ Interface General, Next: C++ Interface Integers, Prev: C++ Class Interface, Up: C++ Class Interface
5535 12.1 C++ Interface General
5536 ==========================
5538 All the C++ classes and functions are available with
5542 Programs should be linked with the 'libgmpxx' and 'libgmp' libraries.
5545 g++ mycxxprog.cc -lgmpxx -lgmp
5547 The classes defined are
5553 The standard operators and various standard functions are overloaded
5554 to allow arithmetic with these classes. For example,
5564 cout << "sum is " << c << "\n";
5565 cout << "absolute value is " << abs(c) << "\n";
5570 An important feature of the implementation is that an expression like
5571 'a=b+c' results in a single call to the corresponding 'mpz_add', without
5572 using a temporary for the 'b+c' part. Expressions which by their nature
5573 imply intermediate values, like 'a=b*c+d*e', still use temporaries
5576 The classes can be freely intermixed in expressions, as can the
5577 classes and the standard types 'long', 'unsigned long' and 'double'.
5578 Smaller types like 'int' or 'float' can also be intermixed, since C++
5581 Note that 'bool' is not accepted directly, but must be explicitly
5582 cast to an 'int' first. This is because C++ will automatically convert
5583 any pointer to a 'bool', so if GMP accepted 'bool' it would make all
5584 sorts of invalid class and pointer combinations compile but almost
5585 certainly not do anything sensible.
5587 Conversions back from the classes to standard C++ types aren't done
5588 automatically, instead member functions like 'get_si' are provided (see
5589 the following sections for details).
5591 Also there are no automatic conversions from the classes to the
5592 corresponding GMP C types, instead a reference to the underlying C
5593 object can be obtained with the following functions,
5595 -- Function: mpz_t mpz_class::get_mpz_t ()
5596 -- Function: mpq_t mpq_class::get_mpq_t ()
5597 -- Function: mpf_t mpf_class::get_mpf_t ()
5599 These can be used to call a C function which doesn't have a C++ class
5600 interface. For example to set 'a' to the GCD of 'b' and 'c',
5604 mpz_gcd (a.get_mpz_t(), b.get_mpz_t(), c.get_mpz_t());
5606 In the other direction, a class can be initialized from the
5607 corresponding GMP C type, or assigned to if an explicit constructor is
5608 used. In both cases this makes a copy of the value, it doesn't create
5609 any sort of association. For example,
5612 // ... init and calculate z ...
5617 There are no namespace setups in 'gmpxx.h', all types and functions
5618 are simply put into the global namespace. This is what 'gmp.h' has done
5619 in the past, and continues to do for compatibility. The extras provided
5620 by 'gmpxx.h' follow GMP naming conventions and are unlikely to clash
5624 File: gmp.info, Node: C++ Interface Integers, Next: C++ Interface Rationals, Prev: C++ Interface General, Up: C++ Class Interface
5626 12.2 C++ Interface Integers
5627 ===========================
5629 -- Function: mpz_class::mpz_class (type N)
5630 Construct an 'mpz_class'. All the standard C++ types may be used,
5631 except 'long long' and 'long double', and all the GMP C++ classes
5632 can be used, although conversions from 'mpq_class' and 'mpf_class'
5633 are 'explicit'. Any necessary conversion follows the corresponding
5634 C function, for example 'double' follows 'mpz_set_d' (*note
5635 Assigning Integers::).
5637 -- Function: explicit mpz_class::mpz_class (const mpz_t Z)
5638 Construct an 'mpz_class' from an 'mpz_t'. The value in Z is copied
5639 into the new 'mpz_class', there won't be any permanent association
5642 -- Function: explicit mpz_class::mpz_class (const char *S, int BASE =
5644 -- Function: explicit mpz_class::mpz_class (const string& S, int BASE =
5646 Construct an 'mpz_class' converted from a string using
5647 'mpz_set_str' (*note Assigning Integers::).
5649 If the string is not a valid integer, an 'std::invalid_argument'
5650 exception is thrown. The same applies to 'operator='.
5652 -- Function: mpz_class operator"" _mpz (const char *STR)
5653 With C++11 compilers, integers can be constructed with the syntax
5654 '123_mpz' which is equivalent to 'mpz_class("123")'.
5656 -- Function: mpz_class operator/ (mpz_class A, mpz_class D)
5657 -- Function: mpz_class operator% (mpz_class A, mpz_class D)
5658 Divisions involving 'mpz_class' round towards zero, as per the
5659 'mpz_tdiv_q' and 'mpz_tdiv_r' functions (*note Integer Division::).
5660 This is the same as the C99 '/' and '%' operators.
5662 The 'mpz_fdiv...' or 'mpz_cdiv...' functions can always be called
5663 directly if desired. For example,
5667 mpz_fdiv_q (q.get_mpz_t(), a.get_mpz_t(), d.get_mpz_t());
5669 -- Function: mpz_class abs (mpz_class OP)
5670 -- Function: int cmp (mpz_class OP1, type OP2)
5671 -- Function: int cmp (type OP1, mpz_class OP2)
5673 -- Function: bool mpz_class::fits_sint_p (void)
5674 -- Function: bool mpz_class::fits_slong_p (void)
5675 -- Function: bool mpz_class::fits_sshort_p (void)
5677 -- Function: bool mpz_class::fits_uint_p (void)
5678 -- Function: bool mpz_class::fits_ulong_p (void)
5679 -- Function: bool mpz_class::fits_ushort_p (void)
5681 -- Function: double mpz_class::get_d (void)
5682 -- Function: long mpz_class::get_si (void)
5683 -- Function: string mpz_class::get_str (int BASE = 10)
5684 -- Function: unsigned long mpz_class::get_ui (void)
5686 -- Function: int mpz_class::set_str (const char *STR, int BASE)
5687 -- Function: int mpz_class::set_str (const string& STR, int BASE)
5688 -- Function: int sgn (mpz_class OP)
5689 -- Function: mpz_class sqrt (mpz_class OP)
5691 -- Function: mpz_class gcd (mpz_class OP1, mpz_class OP2)
5692 -- Function: mpz_class lcm (mpz_class OP1, mpz_class OP2)
5693 -- Function: mpz_class mpz_class::factorial (type OP)
5694 -- Function: mpz_class factorial (mpz_class OP)
5695 -- Function: mpz_class mpz_class::primorial (type OP)
5696 -- Function: mpz_class primorial (mpz_class OP)
5697 -- Function: mpz_class mpz_class::fibonacci (type OP)
5698 -- Function: mpz_class fibonacci (mpz_class OP)
5700 -- Function: void mpz_class::swap (mpz_class& OP)
5701 -- Function: void swap (mpz_class& OP1, mpz_class& OP2)
5702 These functions provide a C++ class interface to the corresponding
5703 GMP C routines. Calling 'factorial' or 'primorial' on a negative
5704 number is undefined.
5706 'cmp' can be used with any of the classes or the standard C++
5707 types, except 'long long' and 'long double'.
5710 Overloaded operators for combinations of 'mpz_class' and 'double' are
5711 provided for completeness, but it should be noted that if the given
5712 'double' is not an integer then the way any rounding is done is
5713 currently unspecified. The rounding might take place at the start, in
5714 the middle, or at the end of the operation, and it might change in the
5717 Conversions between 'mpz_class' and 'double', however, are defined to
5718 follow the corresponding C functions 'mpz_get_d' and 'mpz_set_d'. And
5719 comparisons are always made exactly, as per 'mpz_cmp_d'.
5722 File: gmp.info, Node: C++ Interface Rationals, Next: C++ Interface Floats, Prev: C++ Interface Integers, Up: C++ Class Interface
5724 12.3 C++ Interface Rationals
5725 ============================
5727 In all the following constructors, if a fraction is given then it should
5728 be in canonical form, or if not then 'mpq_class::canonicalize' called.
5730 -- Function: mpq_class::mpq_class (type OP)
5731 -- Function: mpq_class::mpq_class (integer NUM, integer DEN)
5732 Construct an 'mpq_class'. The initial value can be a single value
5733 of any type (conversion from 'mpf_class' is 'explicit'), or a pair
5734 of integers ('mpz_class' or standard C++ integer types)
5735 representing a fraction, except that 'long long' and 'long double'
5736 are not supported. For example,
5742 -- Function: explicit mpq_class::mpq_class (const mpq_t Q)
5743 Construct an 'mpq_class' from an 'mpq_t'. The value in Q is copied
5744 into the new 'mpq_class', there won't be any permanent association
5747 -- Function: explicit mpq_class::mpq_class (const char *S, int BASE =
5749 -- Function: explicit mpq_class::mpq_class (const string& S, int BASE =
5751 Construct an 'mpq_class' converted from a string using
5752 'mpq_set_str' (*note Initializing Rationals::).
5754 If the string is not a valid rational, an 'std::invalid_argument'
5755 exception is thrown. The same applies to 'operator='.
5757 -- Function: mpq_class operator"" _mpq (const char *STR)
5758 With C++11 compilers, integral rationals can be constructed with
5759 the syntax '123_mpq' which is equivalent to 'mpq_class(123_mpz)'.
5760 Other rationals can be built as '-1_mpq/2' or '0xb_mpq/123456_mpz'.
5762 -- Function: void mpq_class::canonicalize ()
5763 Put an 'mpq_class' into canonical form, as per *note Rational
5764 Number Functions::. All arithmetic operators require their
5765 operands in canonical form, and will return results in canonical
5768 -- Function: mpq_class abs (mpq_class OP)
5769 -- Function: int cmp (mpq_class OP1, type OP2)
5770 -- Function: int cmp (type OP1, mpq_class OP2)
5772 -- Function: double mpq_class::get_d (void)
5773 -- Function: string mpq_class::get_str (int BASE = 10)
5775 -- Function: int mpq_class::set_str (const char *STR, int BASE)
5776 -- Function: int mpq_class::set_str (const string& STR, int BASE)
5777 -- Function: int sgn (mpq_class OP)
5779 -- Function: void mpq_class::swap (mpq_class& OP)
5780 -- Function: void swap (mpq_class& OP1, mpq_class& OP2)
5781 These functions provide a C++ class interface to the corresponding
5784 'cmp' can be used with any of the classes or the standard C++
5785 types, except 'long long' and 'long double'.
5787 -- Function: mpz_class& mpq_class::get_num ()
5788 -- Function: mpz_class& mpq_class::get_den ()
5789 Get a reference to an 'mpz_class' which is the numerator or
5790 denominator of an 'mpq_class'. This can be used both for read and
5791 write access. If the object returned is modified, it modifies the
5792 original 'mpq_class'.
5794 If direct manipulation might produce a non-canonical value, then
5795 'mpq_class::canonicalize' must be called before further operations.
5797 -- Function: mpz_t mpq_class::get_num_mpz_t ()
5798 -- Function: mpz_t mpq_class::get_den_mpz_t ()
5799 Get a reference to the underlying 'mpz_t' numerator or denominator
5800 of an 'mpq_class'. This can be passed to C functions expecting an
5801 'mpz_t'. Any modifications made to the 'mpz_t' will modify the
5802 original 'mpq_class'.
5804 If direct manipulation might produce a non-canonical value, then
5805 'mpq_class::canonicalize' must be called before further operations.
5807 -- Function: istream& operator>> (istream& STREAM, mpq_class& ROP);
5808 Read ROP from STREAM, using its 'ios' formatting settings, the same
5809 as 'mpq_t operator>>' (*note C++ Formatted Input::).
5811 If the ROP read might not be in canonical form then
5812 'mpq_class::canonicalize' must be called.
5815 File: gmp.info, Node: C++ Interface Floats, Next: C++ Interface Random Numbers, Prev: C++ Interface Rationals, Up: C++ Class Interface
5817 12.4 C++ Interface Floats
5818 =========================
5820 When an expression requires the use of temporary intermediate
5821 'mpf_class' values, like 'f=g*h+x*y', those temporaries will have the
5822 same precision as the destination 'f'. Explicit constructors can be
5823 used if this doesn't suit.
5825 -- Function: mpf_class::mpf_class (type OP)
5826 -- Function: mpf_class::mpf_class (type OP, mp_bitcnt_t PREC)
5827 Construct an 'mpf_class'. Any standard C++ type can be used,
5828 except 'long long' and 'long double', and any of the GMP C++
5829 classes can be used.
5831 If PREC is given, the initial precision is that value, in bits. If
5832 PREC is not given, then the initial precision is determined by the
5833 type of OP given. An 'mpz_class', 'mpq_class', or C++ builtin type
5834 will give the default 'mpf' precision (*note Initializing
5835 Floats::). An 'mpf_class' or expression will give the precision of
5836 that value. The precision of a binary expression is the higher of
5839 mpf_class f(1.5); // default precision
5840 mpf_class f(1.5, 500); // 500 bits (at least)
5841 mpf_class f(x); // precision of x
5842 mpf_class f(abs(x)); // precision of x
5843 mpf_class f(-g, 1000); // 1000 bits (at least)
5844 mpf_class f(x+y); // greater of precisions of x and y
5846 -- Function: explicit mpf_class::mpf_class (const mpf_t F)
5847 -- Function: mpf_class::mpf_class (const mpf_t F, mp_bitcnt_t PREC)
5848 Construct an 'mpf_class' from an 'mpf_t'. The value in F is copied
5849 into the new 'mpf_class', there won't be any permanent association
5852 If PREC is given, the initial precision is that value, in bits. If
5853 PREC is not given, then the initial precision is that of F.
5855 -- Function: explicit mpf_class::mpf_class (const char *S)
5856 -- Function: mpf_class::mpf_class (const char *S, mp_bitcnt_t PREC, int
5858 -- Function: explicit mpf_class::mpf_class (const string& S)
5859 -- Function: mpf_class::mpf_class (const string& S, mp_bitcnt_t PREC,
5861 Construct an 'mpf_class' converted from a string using
5862 'mpf_set_str' (*note Assigning Floats::). If PREC is given, the
5863 initial precision is that value, in bits. If not, the default
5864 'mpf' precision (*note Initializing Floats::) is used.
5866 If the string is not a valid float, an 'std::invalid_argument'
5867 exception is thrown. The same applies to 'operator='.
5869 -- Function: mpf_class operator"" _mpf (const char *STR)
5870 With C++11 compilers, floats can be constructed with the syntax
5871 '1.23e-1_mpf' which is equivalent to 'mpf_class("1.23e-1")'.
5873 -- Function: mpf_class& mpf_class::operator= (type OP)
5874 Convert and store the given OP value to an 'mpf_class' object. The
5875 same types are accepted as for the constructors above.
5877 Note that 'operator=' only stores a new value, it doesn't copy or
5878 change the precision of the destination, instead the value is
5879 truncated if necessary. This is the same as 'mpf_set' etc. Note
5880 in particular this means for 'mpf_class' a copy constructor is not
5881 the same as a default constructor plus assignment.
5883 mpf_class x (y); // x created with precision of y
5885 mpf_class x; // x created with default precision
5886 x = y; // value truncated to that precision
5888 Applications using templated code may need to be careful about the
5889 assumptions the code makes in this area, when working with
5890 'mpf_class' values of various different or non-default precisions.
5891 For instance implementations of the standard 'complex' template
5892 have been seen in both styles above, though of course 'complex' is
5893 normally only actually specified for use with the builtin float
5896 -- Function: mpf_class abs (mpf_class OP)
5897 -- Function: mpf_class ceil (mpf_class OP)
5898 -- Function: int cmp (mpf_class OP1, type OP2)
5899 -- Function: int cmp (type OP1, mpf_class OP2)
5901 -- Function: bool mpf_class::fits_sint_p (void)
5902 -- Function: bool mpf_class::fits_slong_p (void)
5903 -- Function: bool mpf_class::fits_sshort_p (void)
5905 -- Function: bool mpf_class::fits_uint_p (void)
5906 -- Function: bool mpf_class::fits_ulong_p (void)
5907 -- Function: bool mpf_class::fits_ushort_p (void)
5909 -- Function: mpf_class floor (mpf_class OP)
5910 -- Function: mpf_class hypot (mpf_class OP1, mpf_class OP2)
5912 -- Function: double mpf_class::get_d (void)
5913 -- Function: long mpf_class::get_si (void)
5914 -- Function: string mpf_class::get_str (mp_exp_t& EXP, int BASE = 10,
5916 -- Function: unsigned long mpf_class::get_ui (void)
5918 -- Function: int mpf_class::set_str (const char *STR, int BASE)
5919 -- Function: int mpf_class::set_str (const string& STR, int BASE)
5920 -- Function: int sgn (mpf_class OP)
5921 -- Function: mpf_class sqrt (mpf_class OP)
5923 -- Function: void mpf_class::swap (mpf_class& OP)
5924 -- Function: void swap (mpf_class& OP1, mpf_class& OP2)
5925 -- Function: mpf_class trunc (mpf_class OP)
5926 These functions provide a C++ class interface to the corresponding
5929 'cmp' can be used with any of the classes or the standard C++
5930 types, except 'long long' and 'long double'.
5932 The accuracy provided by 'hypot' is not currently guaranteed.
5934 -- Function: mp_bitcnt_t mpf_class::get_prec ()
5935 -- Function: void mpf_class::set_prec (mp_bitcnt_t PREC)
5936 -- Function: void mpf_class::set_prec_raw (mp_bitcnt_t PREC)
5937 Get or set the current precision of an 'mpf_class'.
5939 The restrictions described for 'mpf_set_prec_raw' (*note
5940 Initializing Floats::) apply to 'mpf_class::set_prec_raw'. Note in
5941 particular that the 'mpf_class' must be restored to it's allocated
5942 precision before being destroyed. This must be done by application
5943 code, there's no automatic mechanism for it.
5946 File: gmp.info, Node: C++ Interface Random Numbers, Next: C++ Interface Limitations, Prev: C++ Interface Floats, Up: C++ Class Interface
5948 12.5 C++ Interface Random Numbers
5949 =================================
5951 -- Class: gmp_randclass
5952 The C++ class interface to the GMP random number functions uses
5953 'gmp_randclass' to hold an algorithm selection and current state,
5954 as per 'gmp_randstate_t'.
5956 -- Function: gmp_randclass::gmp_randclass (void (*RANDINIT)
5957 (gmp_randstate_t, ...), ...)
5958 Construct a 'gmp_randclass', using a call to the given RANDINIT
5959 function (*note Random State Initialization::). The arguments
5960 expected are the same as RANDINIT, but with 'mpz_class' instead of
5961 'mpz_t'. For example,
5963 gmp_randclass r1 (gmp_randinit_default);
5964 gmp_randclass r2 (gmp_randinit_lc_2exp_size, 32);
5965 gmp_randclass r3 (gmp_randinit_lc_2exp, a, c, m2exp);
5966 gmp_randclass r4 (gmp_randinit_mt);
5968 'gmp_randinit_lc_2exp_size' will fail if the size requested is too
5969 big, an 'std::length_error' exception is thrown in that case.
5971 -- Function: gmp_randclass::gmp_randclass (gmp_randalg_t ALG, ...)
5972 Construct a 'gmp_randclass' using the same parameters as
5973 'gmp_randinit' (*note Random State Initialization::). This
5974 function is obsolete and the above RANDINIT style should be
5977 -- Function: void gmp_randclass::seed (unsigned long int S)
5978 -- Function: void gmp_randclass::seed (mpz_class S)
5979 Seed a random number generator. See *note Random Number
5980 Functions::, for how to choose a good seed.
5982 -- Function: mpz_class gmp_randclass::get_z_bits (mp_bitcnt_t BITS)
5983 -- Function: mpz_class gmp_randclass::get_z_bits (mpz_class BITS)
5984 Generate a random integer with a specified number of bits.
5986 -- Function: mpz_class gmp_randclass::get_z_range (mpz_class N)
5987 Generate a random integer in the range 0 to N-1 inclusive.
5989 -- Function: mpf_class gmp_randclass::get_f ()
5990 -- Function: mpf_class gmp_randclass::get_f (mp_bitcnt_t PREC)
5991 Generate a random float F in the range 0 <= F < 1. F will be to
5992 PREC bits precision, or if PREC is not given then to the precision
5993 of the destination. For example,
5997 mpf_class f (0, 512); // 512 bits precision
5998 f = r.get_f(); // random number, 512 bits
6001 File: gmp.info, Node: C++ Interface Limitations, Prev: C++ Interface Random Numbers, Up: C++ Class Interface
6003 12.6 C++ Interface Limitations
6004 ==============================
6006 'mpq_class' and Templated Reading
6007 A generic piece of template code probably won't know that
6008 'mpq_class' requires a 'canonicalize' call if inputs read with
6009 'operator>>' might be non-canonical. This can lead to incorrect
6012 'operator>>' behaves as it does for reasons of efficiency. A
6013 canonicalize can be quite time consuming on large operands, and is
6014 best avoided if it's not necessary.
6016 But this potential difficulty reduces the usefulness of
6017 'mpq_class'. Perhaps a mechanism to tell 'operator>>' what to do
6018 will be adopted in the future, maybe a preprocessor define, a
6019 global flag, or an 'ios' flag pressed into service. Or maybe, at
6020 the risk of inconsistency, the 'mpq_class' 'operator>>' could
6021 canonicalize and leave 'mpq_t' 'operator>>' not doing so, for use
6022 on those occasions when that's acceptable. Send feedback or
6023 alternate ideas to <gmp-bugs@gmplib.org>.
6026 Subclassing the GMP C++ classes works, but is not currently
6029 Expressions involving subclasses resolve correctly (or seem to),
6030 but in normal C++ fashion the subclass doesn't inherit constructors
6031 and assignments. There's many of those in the GMP classes, and a
6032 good way to reestablish them in a subclass is not yet provided.
6034 Templated Expressions
6035 A subtle difficulty exists when using expressions together with
6036 application-defined template functions. Consider the following,
6037 with 'T' intended to be some numeric type,
6040 T fun (const T &, const T &);
6042 When used with, say, plain 'mpz_class' variables, it works fine:
6043 'T' is resolved as 'mpz_class'.
6045 mpz_class f(1), g(2);
6048 But when one of the arguments is an expression, it doesn't work.
6050 mpz_class f(1), g(2), h(3);
6051 fun (f, g+h); // Bad
6053 This is because 'g+h' ends up being a certain expression template
6054 type internal to 'gmpxx.h', which the C++ template resolution rules
6055 are unable to automatically convert to 'mpz_class'. The workaround
6056 is simply to add an explicit cast.
6058 mpz_class f(1), g(2), h(3);
6059 fun (f, mpz_class(g+h)); // Good
6061 Similarly, within 'fun' it may be necessary to cast an expression
6062 to type 'T' when calling a templated 'fun2'.
6067 fun2 (f, f+g); // Bad
6073 fun2 (f, T(f+g)); // Good
6077 C++11 provides several new ways in which types can be inferred:
6078 'auto', 'decltype', etc. While they can be very convenient, they
6079 don't mix well with expression templates. In this example, the
6080 addition is performed twice, as if we had defined 'sum' as a macro.
6084 mpz_class prod = sum * sum;
6086 This other example may crash, though some compilers might make it
6087 look like it is working, because the expression 'z+z' goes out of
6088 scope before it is evaluated.
6091 auto sum = z + z + z;
6092 mpz_class prod = sum * 2;
6094 It is thus strongly recommended to avoid 'auto' anywhere a GMP C++
6095 expression may appear.
6098 File: gmp.info, Node: Custom Allocation, Next: Language Bindings, Prev: C++ Class Interface, Up: Top
6100 13 Custom Allocation
6101 ********************
6103 By default GMP uses 'malloc', 'realloc' and 'free' for memory
6104 allocation, and if they fail GMP prints a message to the standard error
6105 output and terminates the program.
6107 Alternate functions can be specified, to allocate memory in a
6108 different way or to have a different error action on running out of
6111 -- Function: void mp_set_memory_functions (
6112 void *(*ALLOC_FUNC_PTR) (size_t),
6113 void *(*REALLOC_FUNC_PTR) (void *, size_t, size_t),
6114 void (*FREE_FUNC_PTR) (void *, size_t))
6115 Replace the current allocation functions from the arguments. If an
6116 argument is 'NULL', the corresponding default function is used.
6118 These functions will be used for all memory allocation done by GMP,
6119 apart from temporary space from 'alloca' if that function is
6120 available and GMP is configured to use it (*note Build Options::).
6122 *Be sure to call 'mp_set_memory_functions' only when there are no
6123 active GMP objects allocated using the previous memory functions!
6124 Usually that means calling it before any other GMP function.*
6126 The functions supplied should fit the following declarations:
6128 -- Function: void * allocate_function (size_t ALLOC_SIZE)
6129 Return a pointer to newly allocated space with at least ALLOC_SIZE
6132 -- Function: void * reallocate_function (void *PTR, size_t OLD_SIZE,
6134 Resize a previously allocated block PTR of OLD_SIZE bytes to be
6137 The block may be moved if necessary or if desired, and in that case
6138 the smaller of OLD_SIZE and NEW_SIZE bytes must be copied to the
6139 new location. The return value is a pointer to the resized block,
6140 that being the new location if moved or just PTR if not.
6142 PTR is never 'NULL', it's always a previously allocated block.
6143 NEW_SIZE may be bigger or smaller than OLD_SIZE.
6145 -- Function: void free_function (void *PTR, size_t SIZE)
6146 De-allocate the space pointed to by PTR.
6148 PTR is never 'NULL', it's always a previously allocated block of
6151 A "byte" here means the unit used by the 'sizeof' operator.
6153 The REALLOCATE_FUNCTION parameter OLD_SIZE and the FREE_FUNCTION
6154 parameter SIZE are passed for convenience, but of course they can be
6155 ignored if not needed by an implementation. The default functions using
6156 'malloc' and friends for instance don't use them.
6158 No error return is allowed from any of these functions, if they
6159 return then they must have performed the specified operation. In
6160 particular note that ALLOCATE_FUNCTION or REALLOCATE_FUNCTION mustn't
6163 Getting a different fatal error action is a good use for custom
6164 allocation functions, for example giving a graphical dialog rather than
6165 the default print to 'stderr'. How much is possible when genuinely out
6166 of memory is another question though.
6168 There's currently no defined way for the allocation functions to
6169 recover from an error such as out of memory, they must terminate program
6170 execution. A 'longjmp' or throwing a C++ exception will have undefined
6171 results. This may change in the future.
6173 GMP may use allocated blocks to hold pointers to other allocated
6174 blocks. This will limit the assumptions a conservative garbage
6175 collection scheme can make.
6177 Since the default GMP allocation uses 'malloc' and friends, those
6178 functions will be linked in even if the first thing a program does is an
6179 'mp_set_memory_functions'. It's necessary to change the GMP sources if
6183 -- Function: void mp_get_memory_functions (
6184 void *(**ALLOC_FUNC_PTR) (size_t),
6185 void *(**REALLOC_FUNC_PTR) (void *, size_t, size_t),
6186 void (**FREE_FUNC_PTR) (void *, size_t))
6187 Get the current allocation functions, storing function pointers to
6188 the locations given by the arguments. If an argument is 'NULL',
6189 that function pointer is not stored.
6191 For example, to get just the current free function,
6193 void (*freefunc) (void *, size_t);
6195 mp_get_memory_functions (NULL, NULL, &freefunc);
6198 File: gmp.info, Node: Language Bindings, Next: Algorithms, Prev: Custom Allocation, Up: Top
6200 14 Language Bindings
6201 ********************
6203 The following packages and projects offer access to GMP from languages
6204 other than C, though perhaps with varying levels of functionality and
6209 * GMP C++ class interface, *note C++ Class Interface::
6210 Straightforward interface, expression templates to eliminate
6212 * ALP <https://www-sop.inria.fr/saga/logiciels/ALP/>
6213 Linear algebra and polynomials using templates.
6214 * CLN <https://www.ginac.de/CLN/>
6215 High level classes for arithmetic.
6216 * Linbox <http://www.linalg.org/>
6217 Sparse vectors and matrices.
6218 * NTL <http://www.shoup.net/ntl/>
6219 A C++ number theory library.
6222 * Eiffelroom <http://www.eiffelroom.org/node/442>
6225 * Glasgow Haskell Compiler <https://www.haskell.org/ghc/>
6228 * Kaffe <https://github.com/kaffe/kaffe>
6231 * GNU Common Lisp <https://www.gnu.org/software/gcl/gcl.html>
6232 * Librep <http://librep.sourceforge.net/>
6233 * XEmacs (21.5.18 beta and up) <https://www.xemacs.org>
6234 Optional big integers, rationals and floats using GMP.
6237 * MLton compiler <http://mlton.org/>
6240 * MLGMP <https://opam.ocaml.org/packages/mlgmp/>
6241 * Numerix <http://pauillac.inria.fr/~quercia/>
6242 Optionally using GMP.
6245 * Mozart <https://mozart.github.io/>
6248 * GNU Pascal Compiler <http://www.gnu-pascal.de/>
6250 * Numerix <http://pauillac.inria.fr/~quercia/>
6251 For Free Pascal, optionally using GMP.
6254 * GMP module, see 'demos/perl' in the GMP sources (*note
6255 Demonstration Programs::).
6256 * Math::GMP <https://www.cpan.org/>
6257 Compatible with Math::BigInt, but not as many functions as the
6259 * Math::BigInt::GMP <https://www.cpan.org/>
6260 Plug Math::GMP into normal Math::BigInt operations.
6263 * pikempz module in the standard distribution,
6264 <https://pike.lysator.liu.se/>
6267 * SWI Prolog <http://www.swi-prolog.org/>
6268 Arbitrary precision floats.
6271 * GMPY <https://code.google.com/p/gmpy/>
6274 * <https://rubygems.org/gems/gmp>
6277 * GNU Guile <https://www.gnu.org/software/guile/guile.html>
6278 * RScheme <https://www.rscheme.org/>
6279 * STklos <http://www.stklos.net/>
6282 * GNU Smalltalk <http://smalltalk.gnu.org/>
6285 * Axiom <https://savannah.nongnu.org/projects/axiom>
6286 Computer algebra using GCL.
6287 * DrGenius <http://drgenius.seul.org/>
6288 Geometry system and mathematical programming language.
6289 * GiNaC <httsp://www.ginac.de/>
6290 C++ computer algebra using CLN.
6291 * GOO <https://www.eecs.berkeley.edu/~jrb/goo/>
6292 Dynamic object oriented language.
6293 * Maxima <https://www.ma.utexas.edu/users/wfs/maxima.html>
6294 Macsyma computer algebra using GCL.
6295 * Regina <http://regina.sourceforge.net/>
6296 Topological calculator.
6297 * Yacas <http://yacas.sourceforge.net>
6298 Yet another computer algebra system.
6301 File: gmp.info, Node: Algorithms, Next: Internals, Prev: Language Bindings, Up: Top
6306 This chapter is an introduction to some of the algorithms used for
6307 various GMP operations. The code is likely to be hard to understand
6308 without knowing something about the algorithms.
6310 Some GMP internals are mentioned, but applications that expect to be
6311 compatible with future GMP releases should take care to use only the
6312 documented functions.
6316 * Multiplication Algorithms::
6317 * Division Algorithms::
6318 * Greatest Common Divisor Algorithms::
6319 * Powering Algorithms::
6320 * Root Extraction Algorithms::
6321 * Radix Conversion Algorithms::
6322 * Other Algorithms::
6326 File: gmp.info, Node: Multiplication Algorithms, Next: Division Algorithms, Prev: Algorithms, Up: Algorithms
6331 NxN limb multiplications and squares are done using one of seven
6332 algorithms, as the size N increases.
6336 Karatsuba 'MUL_TOOM22_THRESHOLD'
6337 Toom-3 'MUL_TOOM33_THRESHOLD'
6338 Toom-4 'MUL_TOOM44_THRESHOLD'
6339 Toom-6.5 'MUL_TOOM6H_THRESHOLD'
6340 Toom-8.5 'MUL_TOOM8H_THRESHOLD'
6341 FFT 'MUL_FFT_THRESHOLD'
6343 Similarly for squaring, with the 'SQR' thresholds.
6345 NxM multiplications of operands with different sizes above
6346 'MUL_TOOM22_THRESHOLD' are currently done by special Toom-inspired
6347 algorithms or directly with FFT, depending on operand size (*note
6348 Unbalanced Multiplication::).
6352 * Basecase Multiplication::
6353 * Karatsuba Multiplication::
6354 * Toom 3-Way Multiplication::
6355 * Toom 4-Way Multiplication::
6356 * Higher degree Toom'n'half::
6357 * FFT Multiplication::
6358 * Other Multiplication::
6359 * Unbalanced Multiplication::
6362 File: gmp.info, Node: Basecase Multiplication, Next: Karatsuba Multiplication, Prev: Multiplication Algorithms, Up: Multiplication Algorithms
6364 15.1.1 Basecase Multiplication
6365 ------------------------------
6367 Basecase NxM multiplication is a straightforward rectangular set of
6368 cross-products, the same as long multiplication done by hand and for
6369 that reason sometimes known as the schoolbook or grammar school method.
6370 This is an O(N*M) algorithm. See Knuth section 4.3.1 algorithm M (*note
6371 References::), and the 'mpn/generic/mul_basecase.c' code.
6373 Assembly implementations of 'mpn_mul_basecase' are essentially the
6374 same as the generic C code, but have all the usual assembly tricks and
6375 obscurities introduced for speed.
6377 A square can be done in roughly half the time of a multiply, by using
6378 the fact that the cross products above and below the diagonal are the
6379 same. A triangle of products below the diagonal is formed, doubled
6380 (left shift by one bit), and then the products on the diagonal added.
6381 This can be seen in 'mpn/generic/sqr_basecase.c'. Again the assembly
6382 implementations take essentially the same approach.
6385 +---+---+---+---+---+
6387 +---+---+---+---+---+
6389 +---+---+---+---+---+
6391 +---+---+---+---+---+
6393 +---+---+---+---+---+
6395 +---+---+---+---+---+
6397 In practice squaring isn't a full 2x faster than multiplying, it's
6398 usually around 1.5x. Less than 1.5x probably indicates
6399 'mpn_sqr_basecase' wants improving on that CPU.
6401 On some CPUs 'mpn_mul_basecase' can be faster than the generic C
6402 'mpn_sqr_basecase' on some small sizes. 'SQR_BASECASE_THRESHOLD' is the
6403 size at which to use 'mpn_sqr_basecase', this will be zero if that
6404 routine should be used always.
6407 File: gmp.info, Node: Karatsuba Multiplication, Next: Toom 3-Way Multiplication, Prev: Basecase Multiplication, Up: Multiplication Algorithms
6409 15.1.2 Karatsuba Multiplication
6410 -------------------------------
6412 The Karatsuba multiplication algorithm is described in Knuth section
6413 4.3.3 part A, and various other textbooks. A brief description is given
6416 The inputs x and y are treated as each split into two parts of equal
6417 length (or the most significant part one limb shorter if N is odd).
6420 +----------+----------+
6422 +----------+----------+
6424 +----------+----------+
6426 +----------+----------+
6428 Let b be the power of 2 where the split occurs, i.e. if x0 is k limbs
6429 (y0 the same) then b=2^(k*mp_bits_per_limb). With that x=x1*b+x0 and
6430 y=y1*b+y0, and the following holds,
6432 x*y = (b^2+b)*x1*y1 - b*(x1-x0)*(y1-y0) + (b+1)*x0*y0
6434 This formula means doing only three multiplies of (N/2)x(N/2) limbs,
6435 whereas a basecase multiply of NxN limbs is equivalent to four
6436 multiplies of (N/2)x(N/2). The factors (b^2+b) etc represent the
6437 positions where the three products must be added.
6440 +--------+--------+ +--------+--------+
6442 +--------+--------+ +--------+--------+
6450 sub | (x1-x0)*(y1-y0) |
6453 The term (x1-x0)*(y1-y0) is best calculated as an absolute value, and
6454 the sign used to choose to add or subtract. Notice the sum
6455 high(x0*y0)+low(x1*y1) occurs twice, so it's possible to do 5*k limb
6456 additions, rather than 6*k, but in GMP extra function call overheads
6457 outweigh the saving.
6459 Squaring is similar to multiplying, but with x=y the formula reduces
6460 to an equivalent with three squares,
6462 x^2 = (b^2+b)*x1^2 - b*(x1-x0)^2 + (b+1)*x0^2
6464 The final result is accumulated from those three squares the same way
6465 as for the three multiplies above. The middle term (x1-x0)^2 is now
6468 A similar formula for both multiplying and squaring can be
6469 constructed with a middle term (x1+x0)*(y1+y0). But those sums can
6470 exceed k limbs, leading to more carry handling and additions than the
6473 Karatsuba multiplication is asymptotically an O(N^1.585) algorithm,
6474 the exponent being log(3)/log(2), representing 3 multiplies each 1/2 the
6475 size of the inputs. This is a big improvement over the basecase
6476 multiply at O(N^2) and the advantage soon overcomes the extra additions
6477 Karatsuba performs. 'MUL_TOOM22_THRESHOLD' can be as little as 10
6478 limbs. The 'SQR' threshold is usually about twice the 'MUL'.
6480 The basecase algorithm will take a time of the form M(N) = a*N^2 +
6481 b*N + c and the Karatsuba algorithm K(N) = 3*M(N/2) + d*N + e, which
6482 expands to K(N) = 3/4*a*N^2 + 3/2*b*N + 3*c + d*N + e. The factor 3/4
6483 for a means per-crossproduct speedups in the basecase code will increase
6484 the threshold since they benefit M(N) more than K(N). And conversely the
6485 3/2 for b means linear style speedups of b will increase the threshold
6486 since they benefit K(N) more than M(N). The latter can be seen for
6487 instance when adding an optimized 'mpn_sqr_diagonal' to
6488 'mpn_sqr_basecase'. Of course all speedups reduce total time, and in
6489 that sense the algorithm thresholds are merely of academic interest.
6492 File: gmp.info, Node: Toom 3-Way Multiplication, Next: Toom 4-Way Multiplication, Prev: Karatsuba Multiplication, Up: Multiplication Algorithms
6494 15.1.3 Toom 3-Way Multiplication
6495 --------------------------------
6497 The Karatsuba formula is the simplest case of a general approach to
6498 splitting inputs that leads to both Toom and FFT algorithms. A
6499 description of Toom can be found in Knuth section 4.3.3, with an example
6500 3-way calculation after Theorem A. The 3-way form used in GMP is
6503 The operands are each considered split into 3 pieces of equal length
6504 (or the most significant part 1 or 2 limbs shorter than the other two).
6507 +----------+----------+----------+
6509 +----------+----------+----------+
6511 +----------+----------+----------+
6513 +----------+----------+----------+
6515 These parts are treated as the coefficients of two polynomials
6517 X(t) = x2*t^2 + x1*t + x0
6518 Y(t) = y2*t^2 + y1*t + y0
6520 Let b equal the power of 2 which is the size of the x0, x1, y0 and y1
6521 pieces, i.e. if they're k limbs each then b=2^(k*mp_bits_per_limb).
6522 With this x=X(b) and y=Y(b).
6524 Let a polynomial W(t)=X(t)*Y(t) and suppose its coefficients are
6526 W(t) = w4*t^4 + w3*t^3 + w2*t^2 + w1*t + w0
6528 The w[i] are going to be determined, and when they are they'll give
6529 the final result using w=W(b), since x*y=X(b)*Y(b)=W(b). The
6530 coefficients will be roughly b^2 each, and the final W(b) will be an
6550 The w[i] coefficients could be formed by a simple set of cross
6551 products, like w4=x2*y2, w3=x2*y1+x1*y2, w2=x2*y0+x1*y1+x0*y2 etc, but
6552 this would need all nine x[i]*y[j] for i,j=0,1,2, and would be
6553 equivalent merely to a basecase multiply. Instead the following
6556 X(t) and Y(t) are evaluated and multiplied at 5 points, giving values
6557 of W(t) at those points. In GMP the following points are used,
6560 t=0 x0 * y0, which gives w0 immediately
6561 t=1 (x2+x1+x0) * (y2+y1+y0)
6562 t=-1 (x2-x1+x0) * (y2-y1+y0)
6563 t=2 (4*x2+2*x1+x0) * (4*y2+2*y1+y0)
6564 t=inf x2 * y2, which gives w4 immediately
6566 At t=-1 the values can be negative and that's handled using the
6567 absolute values and tracking the sign separately. At t=inf the value is
6568 actually X(t)*Y(t)/t^4 in the limit as t approaches infinity, but it's
6569 much easier to think of as simply x2*y2 giving w4 immediately (much like
6570 x0*y0 at t=0 gives w0 immediately).
6572 Each of the points substituted into W(t)=w4*t^4+...+w0 gives a linear
6573 combination of the w[i] coefficients, and the value of those
6574 combinations has just been calculated.
6577 W(1) = w4 + w3 + w2 + w1 + w0
6578 W(-1) = w4 - w3 + w2 - w1 + w0
6579 W(2) = 16*w4 + 8*w3 + 4*w2 + 2*w1 + w0
6582 This is a set of five equations in five unknowns, and some elementary
6583 linear algebra quickly isolates each w[i]. This involves adding or
6584 subtracting one W(t) value from another, and a couple of divisions by
6585 powers of 2 and one division by 3, the latter using the special
6586 'mpn_divexact_by3' (*note Exact Division::).
6588 The conversion of W(t) values to the coefficients is interpolation.
6589 A polynomial of degree 4 like W(t) is uniquely determined by values
6590 known at 5 different points. The points are arbitrary and can be chosen
6591 to make the linear equations come out with a convenient set of steps for
6592 quickly isolating the w[i].
6594 Squaring follows the same procedure as multiplication, but there's
6595 only one X(t) and it's evaluated at the 5 points, and those values
6596 squared to give values of W(t). The interpolation is then identical,
6597 and in fact the same 'toom_interpolate_5pts' subroutine is used for both
6598 squaring and multiplying.
6600 Toom-3 is asymptotically O(N^1.465), the exponent being
6601 log(5)/log(3), representing 5 recursive multiplies of 1/3 the original
6602 size each. This is an improvement over Karatsuba at O(N^1.585), though
6603 Toom does more work in the evaluation and interpolation and so it only
6604 realizes its advantage above a certain size.
6606 Near the crossover between Toom-3 and Karatsuba there's generally a
6607 range of sizes where the difference between the two is small.
6608 'MUL_TOOM33_THRESHOLD' is a somewhat arbitrary point in that range and
6609 successive runs of the tune program can give different values due to
6610 small variations in measuring. A graph of time versus size for the two
6611 shows the effect, see 'tune/README'.
6613 At the fairly small sizes where the Toom-3 thresholds occur it's
6614 worth remembering that the asymptotic behaviour for Karatsuba and Toom-3
6615 can't be expected to make accurate predictions, due of course to the big
6616 influence of all sorts of overheads, and the fact that only a few
6617 recursions of each are being performed. Even at large sizes there's a
6618 good chance machine dependent effects like cache architecture will mean
6619 actual performance deviates from what might be predicted.
6621 The formula given for the Karatsuba algorithm (*note Karatsuba
6622 Multiplication::) has an equivalent for Toom-3 involving only five
6623 multiplies, but this would be complicated and unenlightening.
6625 An alternate view of Toom-3 can be found in Zuras (*note
6626 References::), using a vector to represent the x and y splits and a
6627 matrix multiplication for the evaluation and interpolation stages. The
6628 matrix inverses are not meant to be actually used, and they have
6629 elements with values much greater than in fact arise in the
6630 interpolation steps. The diagram shown for the 3-way is attractive, but
6631 again doesn't have to be implemented that way and for example with a bit
6632 of rearrangement just one division by 6 can be done.
6635 File: gmp.info, Node: Toom 4-Way Multiplication, Next: Higher degree Toom'n'half, Prev: Toom 3-Way Multiplication, Up: Multiplication Algorithms
6637 15.1.4 Toom 4-Way Multiplication
6638 --------------------------------
6640 Karatsuba and Toom-3 split the operands into 2 and 3 coefficients,
6641 respectively. Toom-4 analogously splits the operands into 4
6642 coefficients. Using the notation from the section on Toom-3
6643 multiplication, we form two polynomials:
6645 X(t) = x3*t^3 + x2*t^2 + x1*t + x0
6646 Y(t) = y3*t^3 + y2*t^2 + y1*t + y0
6648 X(t) and Y(t) are evaluated and multiplied at 7 points, giving values
6649 of W(t) at those points. In GMP the following points are used,
6652 t=0 x0 * y0, which gives w0 immediately
6653 t=1/2 (x3+2*x2+4*x1+8*x0) * (y3+2*y2+4*y1+8*y0)
6654 t=-1/2 (-x3+2*x2-4*x1+8*x0) * (-y3+2*y2-4*y1+8*y0)
6655 t=1 (x3+x2+x1+x0) * (y3+y2+y1+y0)
6656 t=-1 (-x3+x2-x1+x0) * (-y3+y2-y1+y0)
6657 t=2 (8*x3+4*x2+2*x1+x0) * (8*y3+4*y2+2*y1+y0)
6658 t=inf x3 * y3, which gives w6 immediately
6660 The number of additions and subtractions for Toom-4 is much larger
6661 than for Toom-3. But several subexpressions occur multiple times, for
6662 example x2+x0, occurs for both t=1 and t=-1.
6664 Toom-4 is asymptotically O(N^1.404), the exponent being
6665 log(7)/log(4), representing 7 recursive multiplies of 1/4 the original
6669 File: gmp.info, Node: Higher degree Toom'n'half, Next: FFT Multiplication, Prev: Toom 4-Way Multiplication, Up: Multiplication Algorithms
6671 15.1.5 Higher degree Toom'n'half
6672 --------------------------------
6674 The Toom algorithms described above (*note Toom 3-Way Multiplication::,
6675 *note Toom 4-Way Multiplication::) generalizes to split into an
6676 arbitrary number of pieces. In general a split of two equally long
6677 operands into r pieces leads to evaluations and pointwise
6678 multiplications done at 2*r-1 points. To fully exploit symmetries it
6679 would be better to have a multiple of 4 points, that's why for higher
6680 degree Toom'n'half is used.
6682 Toom'n'half means that the existence of one more piece is considered
6683 for a single operand. It can be virtual, i.e. zero, or real, when the
6684 two operand are not exactly balanced. By choosing an even r, Toom-r+1/2
6685 requires 2r points, a multiple of four.
6687 The quadruplets of points include 0, inf, +1, -1 and +-2^i, +-2^-i .
6688 Each of them giving shortcuts for the evaluation phase and for some
6689 steps in the interpolation phase. Further tricks are used to reduce the
6690 memory footprint of the whole multiplication algorithm to a memory
6691 buffer equal in size to the result of the product.
6693 Current GMP uses both Toom-6'n'half and Toom-8'n'half.
6696 File: gmp.info, Node: FFT Multiplication, Next: Other Multiplication, Prev: Higher degree Toom'n'half, Up: Multiplication Algorithms
6698 15.1.6 FFT Multiplication
6699 -------------------------
6701 At large to very large sizes a Fermat style FFT multiplication is used,
6702 following Schönhage and Strassen (*note References::). Descriptions of
6703 FFTs in various forms can be found in many textbooks, for instance Knuth
6704 section 4.3.3 part C or Lipson chapter IX. A brief description of the
6705 form used in GMP is given here.
6707 The multiplication done is x*y mod 2^N+1, for a given N. A full
6708 product x*y is obtained by choosing N>=bits(x)+bits(y) and padding x and
6709 y with high zero limbs. The modular product is the native form for the
6710 algorithm, so padding to get a full product is unavoidable.
6712 The algorithm follows a split, evaluate, pointwise multiply,
6713 interpolate and combine similar to that described above for Karatsuba
6714 and Toom-3. A k parameter controls the split, with an FFT-k splitting
6715 into 2^k pieces of M=N/2^k bits each. N must be a multiple of
6716 (2^k)*mp_bits_per_limb so the split falls on limb boundaries, avoiding
6717 bit shifts in the split and combine stages.
6719 The evaluations, pointwise multiplications, and interpolation, are
6720 all done modulo 2^N'+1 where N' is 2M+k+3 rounded up to a multiple of
6721 2^k and of 'mp_bits_per_limb'. The results of interpolation will be the
6722 following negacyclic convolution of the input pieces, and the choice of
6723 N' ensures these sums aren't truncated.
6727 w[n] = / (-1) * x[i] * y[j]
6732 The points used for the evaluation are g^i for i=0 to 2^k-1 where
6733 g=2^(2N'/2^k). g is a 2^k'th root of unity mod 2^N'+1, which produces
6734 necessary cancellations at the interpolation stage, and it's also a
6735 power of 2 so the fast Fourier transforms used for the evaluation and
6736 interpolation do only shifts, adds and negations.
6738 The pointwise multiplications are done modulo 2^N'+1 and either
6739 recurse into a further FFT or use a plain multiplication (Toom-3,
6740 Karatsuba or basecase), whichever is optimal at the size N'. The
6741 interpolation is an inverse fast Fourier transform. The resulting set
6742 of sums of x[i]*y[j] are added at appropriate offsets to give the final
6745 Squaring is the same, but x is the only input so it's one transform
6746 at the evaluate stage and the pointwise multiplies are squares. The
6747 interpolation is the same.
6749 For a mod 2^N+1 product, an FFT-k is an O(N^(k/(k-1))) algorithm, the
6750 exponent representing 2^k recursed modular multiplies each 1/2^(k-1) the
6751 size of the original. Each successive k is an asymptotic improvement,
6752 but overheads mean each is only faster at bigger and bigger sizes. In
6753 the code, 'MUL_FFT_TABLE' and 'SQR_FFT_TABLE' are the thresholds where
6754 each k is used. Each new k effectively swaps some multiplying for some
6755 shifts, adds and overheads.
6757 A mod 2^N+1 product can be formed with a normal NxN->2N bit multiply
6758 plus a subtraction, so an FFT and Toom-3 etc can be compared directly.
6759 A k=4 FFT at O(N^1.333) can be expected to be the first faster than
6760 Toom-3 at O(N^1.465). In practice this is what's found, with
6761 'MUL_FFT_MODF_THRESHOLD' and 'SQR_FFT_MODF_THRESHOLD' being between 300
6762 and 1000 limbs, depending on the CPU. So far it's been found that only
6763 very large FFTs recurse into pointwise multiplies above these sizes.
6765 When an FFT is to give a full product, the change of N to 2N doesn't
6766 alter the theoretical complexity for a given k, but for the purposes of
6767 considering where an FFT might be first used it can be assumed that the
6768 FFT is recursing into a normal multiply and that on that basis it's
6769 doing 2^k recursed multiplies each 1/2^(k-2) the size of the inputs,
6770 making it O(N^(k/(k-2))). This would mean k=7 at O(N^1.4) would be the
6771 first FFT faster than Toom-3. In practice 'MUL_FFT_THRESHOLD' and
6772 'SQR_FFT_THRESHOLD' have been found to be in the k=8 range, somewhere
6773 between 3000 and 10000 limbs.
6775 The way N is split into 2^k pieces and then 2M+k+3 is rounded up to a
6776 multiple of 2^k and 'mp_bits_per_limb' means that when
6777 2^k>=mp\_bits\_per\_limb the effective N is a multiple of 2^(2k-1) bits.
6778 The +k+3 means some values of N just under such a multiple will be
6779 rounded to the next. The complexity calculations above assume that a
6780 favourable size is used, meaning one which isn't padded through
6781 rounding, and it's also assumed that the extra +k+3 bits are negligible
6782 at typical FFT sizes.
6784 The practical effect of the 2^(2k-1) constraint is to introduce a
6785 step-effect into measured speeds. For example k=8 will round N up to a
6786 multiple of 32768 bits, so for a 32-bit limb there'll be 512 limb groups
6787 of sizes for which 'mpn_mul_n' runs at the same speed. Or for k=9
6788 groups of 2048 limbs, k=10 groups of 8192 limbs, etc. In practice it's
6789 been found each k is used at quite small multiples of its size
6790 constraint and so the step effect is quite noticeable in a time versus
6793 The threshold determinations currently measure at the mid-points of
6794 size steps, but this is sub-optimal since at the start of a new step it
6795 can happen that it's better to go back to the previous k for a while.
6796 Something more sophisticated for 'MUL_FFT_TABLE' and 'SQR_FFT_TABLE'
6800 File: gmp.info, Node: Other Multiplication, Next: Unbalanced Multiplication, Prev: FFT Multiplication, Up: Multiplication Algorithms
6802 15.1.7 Other Multiplication
6803 ---------------------------
6805 The Toom algorithms described above (*note Toom 3-Way Multiplication::,
6806 *note Toom 4-Way Multiplication::) generalizes to split into an
6807 arbitrary number of pieces, as per Knuth section 4.3.3 algorithm C.
6808 This is not currently used. The notes here are merely for interest.
6810 In general a split into r+1 pieces is made, and evaluations and
6811 pointwise multiplications done at 2*r+1 points. A 4-way split does 7
6812 pointwise multiplies, 5-way does 9, etc. Asymptotically an (r+1)-way
6813 algorithm is O(N^(log(2*r+1)/log(r+1))). Only the pointwise
6814 multiplications count towards big-O complexity, but the time spent in
6815 the evaluate and interpolate stages grows with r and has a significant
6816 practical impact, with the asymptotic advantage of each r realized only
6817 at bigger and bigger sizes. The overheads grow as O(N*r), whereas in an
6818 r=2^k FFT they grow only as O(N*log(r)).
6820 Knuth algorithm C evaluates at points 0,1,2,...,2*r, but exercise 4
6821 uses -r,...,0,...,r and the latter saves some small multiplies in the
6822 evaluate stage (or rather trades them for additions), and has a further
6823 saving of nearly half the interpolate steps. The idea is to separate
6824 odd and even final coefficients and then perform algorithm C steps C7
6825 and C8 on them separately. The divisors at step C7 become j^2 and the
6826 multipliers at C8 become 2*t*j-j^2.
6828 Splitting odd and even parts through positive and negative points can
6829 be thought of as using -1 as a square root of unity. If a 4th root of
6830 unity was available then a further split and speedup would be possible,
6831 but no such root exists for plain integers. Going to complex integers
6832 with i=sqrt(-1) doesn't help, essentially because in Cartesian form it
6833 takes three real multiplies to do a complex multiply. The existence of
6834 2^k'th roots of unity in a suitable ring or field lets the fast Fourier
6835 transform keep splitting and get to O(N*log(r)).
6837 Floating point FFTs use complex numbers approximating Nth roots of
6838 unity. Some processors have special support for such FFTs. But these
6839 are not used in GMP since it's very difficult to guarantee an exact
6840 result (to some number of bits). An occasional difference of 1 in the
6841 last bit might not matter to a typical signal processing algorithm, but
6842 is of course of vital importance to GMP.
6845 File: gmp.info, Node: Unbalanced Multiplication, Prev: Other Multiplication, Up: Multiplication Algorithms
6847 15.1.8 Unbalanced Multiplication
6848 --------------------------------
6850 Multiplication of operands with different sizes, both below
6851 'MUL_TOOM22_THRESHOLD' are done with plain schoolbook multiplication
6852 (*note Basecase Multiplication::).
6854 For really large operands, we invoke FFT directly.
6856 For operands between these sizes, we use Toom inspired algorithms
6857 suggested by Alberto Zanoni and Marco Bodrato. The idea is to split the
6858 operands into polynomials of different degree. GMP currently splits the
6859 smaller operand onto 2 coefficients, i.e., a polynomial of degree 1, but
6860 the larger operand can be split into 2, 3, or 4 coefficients, i.e., a
6861 polynomial of degree 1 to 3.
6864 File: gmp.info, Node: Division Algorithms, Next: Greatest Common Divisor Algorithms, Prev: Multiplication Algorithms, Up: Algorithms
6866 15.2 Division Algorithms
6867 ========================
6871 * Single Limb Division::
6872 * Basecase Division::
6873 * Divide and Conquer Division::
6874 * Block-Wise Barrett Division::
6877 * Small Quotient Division::
6880 File: gmp.info, Node: Single Limb Division, Next: Basecase Division, Prev: Division Algorithms, Up: Division Algorithms
6882 15.2.1 Single Limb Division
6883 ---------------------------
6885 Nx1 division is implemented using repeated 2x1 divisions from high to
6886 low, either with a hardware divide instruction or a multiplication by
6887 inverse, whichever is best on a given CPU.
6889 The multiply by inverse follows "Improved division by invariant
6890 integers" by Möller and Granlund (*note References::) and is implemented
6891 as 'udiv_qrnnd_preinv' in 'gmp-impl.h'. The idea is to have a
6892 fixed-point approximation to 1/d (see 'invert_limb') and then multiply
6893 by the high limb (plus one bit) of the dividend to get a quotient q.
6894 With d normalized (high bit set), q is no more than 1 too small.
6895 Subtracting q*d from the dividend gives a remainder, and reveals whether
6896 q or q-1 is correct.
6898 The result is a division done with two multiplications and four or
6899 five arithmetic operations. On CPUs with low latency multipliers this
6900 can be much faster than a hardware divide, though the cost of
6901 calculating the inverse at the start may mean it's only better on inputs
6902 bigger than say 4 or 5 limbs.
6904 When a divisor must be normalized, either for the generic C
6905 '__udiv_qrnnd_c' or the multiply by inverse, the division performed is
6906 actually a*2^k by d*2^k where a is the dividend and k is the power
6907 necessary to have the high bit of d*2^k set. The bit shifts for the
6908 dividend are usually accomplished "on the fly" meaning by extracting the
6909 appropriate bits at each step. Done this way the quotient limbs come
6910 out aligned ready to store. When only the remainder is wanted, an
6911 alternative is to take the dividend limbs unshifted and calculate r = a
6912 mod d*2^k followed by an extra final step r*2^k mod d*2^k. This can
6913 help on CPUs with poor bit shifts or few registers.
6915 The multiply by inverse can be done two limbs at a time. The
6916 calculation is basically the same, but the inverse is two limbs and the
6917 divisor treated as if padded with a low zero limb. This means more
6918 work, since the inverse will need a 2x2 multiply, but the four 1x1s to
6919 do that are independent and can therefore be done partly or wholly in
6920 parallel. Likewise for a 2x1 calculating q*d. The net effect is to
6921 process two limbs with roughly the same two multiplies worth of latency
6922 that one limb at a time gives. This extends to 3 or 4 limbs at a time,
6923 though the extra work to apply the inverse will almost certainly soon
6924 reach the limits of multiplier throughput.
6926 A similar approach in reverse can be taken to process just half a
6927 limb at a time if the divisor is only a half limb. In this case the 1x1
6928 multiply for the inverse effectively becomes two (1/2)x1 for each limb,
6929 which can be a saving on CPUs with a fast half limb multiply, or in fact
6930 if the only multiply is a half limb, and especially if it's not
6934 File: gmp.info, Node: Basecase Division, Next: Divide and Conquer Division, Prev: Single Limb Division, Up: Division Algorithms
6936 15.2.2 Basecase Division
6937 ------------------------
6939 Basecase NxM division is like long division done by hand, but in base
6940 2^mp_bits_per_limb. See Knuth section 4.3.1 algorithm D, and
6941 'mpn/generic/sb_divrem_mn.c'.
6943 Briefly stated, while the dividend remains larger than the divisor, a
6944 high quotient limb is formed and the Nx1 product q*d subtracted at the
6945 top end of the dividend. With a normalized divisor (most significant
6946 bit set), each quotient limb can be formed with a 2x1 division and a 1x1
6947 multiplication plus some subtractions. The 2x1 division is by the high
6948 limb of the divisor and is done either with a hardware divide or a
6949 multiply by inverse (the same as in *note Single Limb Division::)
6950 whichever is faster. Such a quotient is sometimes one too big,
6951 requiring an addback of the divisor, but that happens rarely.
6953 With Q=N-M being the number of quotient limbs, this is an O(Q*M)
6954 algorithm and will run at a speed similar to a basecase QxM
6955 multiplication, differing in fact only in the extra multiply and divide
6956 for each of the Q quotient limbs.
6959 File: gmp.info, Node: Divide and Conquer Division, Next: Block-Wise Barrett Division, Prev: Basecase Division, Up: Division Algorithms
6961 15.2.3 Divide and Conquer Division
6962 ----------------------------------
6964 For divisors larger than 'DC_DIV_QR_THRESHOLD', division is done by
6965 dividing. Or to be precise by a recursive divide and conquer algorithm
6966 based on work by Moenck and Borodin, Jebelean, and Burnikel and Ziegler
6967 (*note References::).
6969 The algorithm consists essentially of recognising that a 2NxN
6970 division can be done with the basecase division algorithm (*note
6971 Basecase Division::), but using N/2 limbs as a base, not just a single
6972 limb. This way the multiplications that arise are (N/2)x(N/2) and can
6973 take advantage of Karatsuba and higher multiplication algorithms (*note
6974 Multiplication Algorithms::). The two "digits" of the quotient are
6975 formed by recursive Nx(N/2) divisions.
6977 If the (N/2)x(N/2) multiplies are done with a basecase multiplication
6978 then the work is about the same as a basecase division, but with more
6979 function call overheads and with some subtractions separated from the
6980 multiplies. These overheads mean that it's only when N/2 is above
6981 'MUL_TOOM22_THRESHOLD' that divide and conquer is of use.
6983 'DC_DIV_QR_THRESHOLD' is based on the divisor size N, so it will be
6984 somewhere above twice 'MUL_TOOM22_THRESHOLD', but how much above depends
6985 on the CPU. An optimized 'mpn_mul_basecase' can lower
6986 'DC_DIV_QR_THRESHOLD' a little by offering a ready-made advantage over
6987 repeated 'mpn_submul_1' calls.
6989 Divide and conquer is asymptotically O(M(N)*log(N)) where M(N) is the
6990 time for an NxN multiplication done with FFTs. The actual time is a sum
6991 over multiplications of the recursed sizes, as can be seen near the end
6992 of section 2.2 of Burnikel and Ziegler. For example, within the Toom-3
6993 range, divide and conquer is 2.63*M(N). With higher algorithms the M(N)
6994 term improves and the multiplier tends to log(N). In practice, at
6995 moderate to large sizes, a 2NxN division is about 2 to 4 times slower
6996 than an NxN multiplication.
6999 File: gmp.info, Node: Block-Wise Barrett Division, Next: Exact Division, Prev: Divide and Conquer Division, Up: Division Algorithms
7001 15.2.4 Block-Wise Barrett Division
7002 ----------------------------------
7004 For the largest divisions, a block-wise Barrett division algorithm is
7005 used. Here, the divisor is inverted to a precision determined by the
7006 relative size of the dividend and divisor. Blocks of quotient limbs are
7007 then generated by multiplying blocks from the dividend by the inverse.
7009 Our block-wise algorithm computes a smaller inverse than in the plain
7010 Barrett algorithm. For a 2n/n division, the inverse will be just
7014 File: gmp.info, Node: Exact Division, Next: Exact Remainder, Prev: Block-Wise Barrett Division, Up: Division Algorithms
7016 15.2.5 Exact Division
7017 ---------------------
7019 A so-called exact division is when the dividend is known to be an exact
7020 multiple of the divisor. Jebelean's exact division algorithm uses this
7021 knowledge to make some significant optimizations (*note References::).
7023 The idea can be illustrated in decimal for example with 368154
7024 divided by 543. Because the low digit of the dividend is 4, the low
7025 digit of the quotient must be 8. This is arrived at from 4*7 mod 10,
7026 using the fact 7 is the modular inverse of 3 (the low digit of the
7027 divisor), since 3*7 == 1 mod 10. So 8*543=4344 can be subtracted from
7028 the dividend leaving 363810. Notice the low digit has become zero.
7030 The procedure is repeated at the second digit, with the next quotient
7031 digit 7 (7 == 1*7 mod 10), subtracting 7*543=3801, leaving 325800. And
7032 finally at the third digit with quotient digit 6 (8*7 mod 10),
7033 subtracting 6*543=3258 leaving 0. So the quotient is 678.
7035 Notice however that the multiplies and subtractions don't need to
7036 extend past the low three digits of the dividend, since that's enough to
7037 determine the three quotient digits. For the last quotient digit no
7038 subtraction is needed at all. On a 2NxN division like this one, only
7039 about half the work of a normal basecase division is necessary.
7041 For an NxM exact division producing Q=N-M quotient limbs, the saving
7042 over a normal basecase division is in two parts. Firstly, each of the Q
7043 quotient limbs needs only one multiply, not a 2x1 divide and multiply.
7044 Secondly, the crossproducts are reduced when Q>M to Q*M-M*(M+1)/2, or
7045 when Q<=M to Q*(Q-1)/2. Notice the savings are complementary. If Q is
7046 big then many divisions are saved, or if Q is small then the
7047 crossproducts reduce to a small number.
7049 The modular inverse used is calculated efficiently by 'binvert_limb'
7050 in 'gmp-impl.h'. This does four multiplies for a 32-bit limb, or six
7051 for a 64-bit limb. 'tune/modlinv.c' has some alternate implementations
7052 that might suit processors better at bit twiddling than multiplying.
7054 The sub-quadratic exact division described by Jebelean in "Exact
7055 Division with Karatsuba Complexity" is not currently implemented. It
7056 uses a rearrangement similar to the divide and conquer for normal
7057 division (*note Divide and Conquer Division::), but operating from low
7058 to high. A further possibility not currently implemented is
7059 "Bidirectional Exact Integer Division" by Krandick and Jebelean which
7060 forms quotient limbs from both the high and low ends of the dividend,
7061 and can halve once more the number of crossproducts needed in a 2NxN
7064 A special case exact division by 3 exists in 'mpn_divexact_by3',
7065 supporting Toom-3 multiplication and 'mpq' canonicalizations. It forms
7066 quotient digits with a multiply by the modular inverse of 3 (which is
7067 '0xAA..AAB') and uses two comparisons to determine a borrow for the next
7068 limb. The multiplications don't need to be on the dependent chain, as
7069 long as the effect of the borrows is applied, which can help chips with
7070 pipelined multipliers.
projects / xonotic / xonotic.git / blob