Using the Multiple Precision Library
Dealing with infinite-precision integers
John Rogers
JR is a programmer in the Seattle area. He can be contacted on CompuServe at 72634,2402.
ANSI C is a fine language--as long as you don't have to count too high. Although most implementations of C boast 32-bit integer arithmetic, many applications are beginning to require higher precision. (Microsoft, for instance, makes over a billion dollars a year, and that dollar amount barely fits in 32 bits. What if Bill wants to count in pennies to keep the auditors happy?) However, even the IEEE double-precision floating-point format only gives 15 digits of precision, as I understand it.
What's needed is a way to deal with multiple-precision integers, independent of the machine's word size. The multiple-precision (MP) integer library available with UNIX V7, UNIX SVR4, 4.3BSD, and other versions of UNIX provide "infinite-precision" signed integer operations for C programs. I've also ported the GMP library to Windows NT. In this article, I'll describe how to use the scantily documented MP routines, along with providing sample code, portability information, some MP "helpers," and a few other hints.
Declaring Variables
To declare variables for use with the MP library, you generally declare them as pointers to the MINT (multiple-precision integer) type. The MINT type is often defined in <mp.h> as a structure. Even if the fields in the structure are documented in your version of MP, I recommend avoiding direct references to them in your programs. The field names and usage differ greatly between MP implementations.
You should also be aware that the MINT type name (all uppercase letters) is only available as mint (all lowercase letters) in a few versions. I recommend you define NEED_MINT when using those versions, and write code like that shown in Example 1; after that, you can use the MINT type. I do this in the header file for my MP helpers (see MPHelp.h in Listing One). In this article, I will only use the MINT type name.
Table 1 lists the MP functions supported by various versions of MP, along with entries for the type name(s) supported by each, and whether or not the header file <mp.h> is available. Limitations of the various versions are also noted. The details in this table should help you write portable code using MP.
One rule is important when writing code with MINT variables: Make sure every variable is initialized by calling itom or xtom, both of which return pointers to MINT variables. Throughout this article, I use phrases like "refers to a MINT variable" when I really mean "is passed a pointer to a MINT variable." In short, none of the MP functions are passed MINT variables directly.
The itom (integer-to-MINT) routine is available in every version of MP. It is by far the most common way to initialize MINT variables, as itom returns a pointer to a MINT. Beware that itom really only accepts a value of type short in most versions (despite what some man pages imply), and at least one version will not accept the "most negative short integer" (SHRT_MIN, I assume). Being allowed to initialize only with a short does not limit the final size of the variable.
The xtom (hex-string-to-MINT) routine is not widely available, and it is the only way to initialize a MINT variable besides itom. Therefore, xtom returns a pointer to a MINT variable. The only argument to xtom is a pointer to a string containing hex digits. For example, to declare and initialize some variables, use code like that in Example 2. Once you have initialized a MINT variable, you can update its value repeatedly with a variety of MP functions. The MP library will automatically expand and shrink the memory allocated for that MINT variable. In some versions, you can explicitly give back the memory associated with a variable by using the mfree function. Note that mfree effectively deinitializes the MINT variable.
You should be aware of one other aspect related to memory allocation: MP routines do not inform your code in any way if they run out of memory. There's no return code, no global-error variable, no callback, no signal, no exception. The routines die in whatever manner they choose, probably writing a message to stderr in English, and possibly leaving a core file if they are running on a UNIX system. MP routines are fine for many applications, but think twice before using them in something like an air-traffic control system.
Passing out-of-range arguments to the MP routines will generally have similar results as running out of memory. So, write robust code, and test it well. Beware of passing references to the same MINT as input and output for the same function call. (Using the same MINT for two different inputs seems to work fine, however.)
The MP compare (mcmp) routine and the standard-C string-compare (strcmp) routine are very similar. Both return type int when comparing two values. Both return 0, when the values are equal; less than 0, when the first value is smaller; and greater than 0, when the first value is larger. See MPHelp.c (Listing Two) for samples of using the mcmp routine.
Add, Subtract, Multiply, Divide
Every version of MP has support for the basic math operations add, subtract, multiply, and divide. The MINT add routine (madd) takes its first two arguments as inputs, adds them, and stores the result in the third argument. The madd routine, like almost every other MP routine, has no return value and assumes that all of its arguments (even the output arguments) have already been initialized. The MINT subtract routine, msub, is similar to madd in Example 3. It subtracts the value of its second argument from its first argument and sets its third argument accordingly. The MINT multiply routine (mult) simply multiplies its first argument by its second argument and updates the third argument with the product. All three arguments must refer to MINT variables.
The MP library has two divide routines that differ only in name and the data types of two arguments. For the MINT divide routine, mdiv, all arguments refer to MINT variables. For the short divide routine, sdiv, the data type of the divisor is short. The remainder in sdiv is represented as a pointer to a short. Other than those two exceptions in sdiv, all other arguments to the divide routines must refer to MINT variables. For both routines, the first two arguments are inputs. The first argument is the dividend and the second is the divisor. The divide routines also have two outputs each. The third argument will be set to the quotient, and the fourth will be set to the remainder. Be careful when passing negative numbers to the divide routines, because various versions of mdiv and sdiv set the sign for the remainder using different rules. As with itom, be aware that at least one version of sdiv will not accept the "most negative short integer" (SHRT_MIN, I assume) as the divisor. Listing Two (MPHelp.c) gives examples of calling both MP divide routines, mdiv and sdiv.
I/O
The MP versions of which I am aware include eight different MINT I/O routines. However, only two of those routines are available in all versions, and one of those two is plagued by problems. The two common routines are MINT output (mout) and MINT input (min). The mout and min routines do decimal I/O using stdout. The min routine in some MP versions may return an int value of EOF (end of file), but this return value is often not documented or missing altogether. For portability, I recommend avoiding the min return value and calling the stdio feof function instead. The feof routine and the EOF #define are declared in <stdio.h>. The name min is doubly cursed; in some older versions of Microsoft C an incompatible min macro is defined in <stdlib.h>.
If you have access to them, I recommend the Berkeley-style m_in and m_out routines, which support any file (not just stdin/stdout) and more number bases (2 through 10 in 4.3BSD).
Table 2 gives more-detailed information on all eight of the MP I/O functions. It documents the return values, argument lists, number bases, and files supported by each function. IsPrime.c (Listing Three) provides a sample call to mout, while PowTable.c (Listing Four) contains calls to m_out.
Square Roots
All MP versions have a MINT square-root routine, msqrt. The only input is the first argument, which must refer to a MINT having a value greater than or equal to 0. The outputs are the other two arguments. The second argument must point to a MINT set with the square root, and the third argument must refer to a MINT set with the remainder. MPHelp.c (Listing Two) illustrates a call to msqrt.
GCD
Although the function name may differ, every version of MP has a greatest common divisor (GCD) routine, usually named gcd. In one version it is mgcd, probably to avoid conflicts with an incompatible math-library routine.
The gcd and mgcd routines each have two inputs and one output. All arguments refer to MINT variables. The first and second arguments are inputs; the third will be set to the greatest common divisor of the other two. I recommend that you avoid calling either routine with inputs that have 0 or negative values.
Powers
Every version of MP has an rpow (regular power) routine. The first input (referring to a MINT) is the base, the second is the exponent. The data type of the exponent is often undocumented. I recommend only passing a short for this argument. The third argument must refer to a MINT that will be set to the first argument raised to the second argument. Try to avoid 0 and negative powers.
All versions of MP have another power function, usually called pow. At least one MP version calls it mpow (also because of a math-library-function incompatibility). The arguments for pow and mpow are identical; all refer to MINT values. The first three arguments are all inputs. The first argument is the base, the second is the exponent, and the third is a value the result of which is to be computed modulo. The fourth argument is the output; it will be set with the result. Unfortunately, the "modulo" argument of this function makes it more or less useless.
The code sample for the power routines produces a table of the powers of two, in various number bases, up to any arbitrary number. (This is part of the reason for my interest in the MP routines. I'm writing a reference book and I want to include a table of powers of two, up to at least 64 bits, but my C compiler quits at 32.) See PowTable.c (Listing Four) for an example of calling rpow.
Miscellaneous Routines
One miscellaneous routine I've found useful is move, which copies one MINT value to another. The first argument is the input; the second is the output. This routine is often used because it is not wise to pass the same MINT value to an MP function as both an input and an output. Instead you introduce temporary variables and call the move routine here and there.
The Berkeley version of MP also provides a function named invert which gets three arguments, all of which must refer to MINT variables. The first two arguments are inputs; the third is the output. The 4.3BSD man page says "invert computes c such that a*c mod b=1, for a and b relatively prime." I don't know of any use for this function yet; I'm only including it for completeness.
The last miscellaneous MP function is the MINT-to-hex-string (mtox) routine. The only argument to mtox is an input, referring to the MINT variable to be converted. The return value from mtox is a pointer to a newly allocated string. The area allocated for the string may be released by calling the usual function (free, declared in <stdlib.h>).
Conclusions
The MP library includes some powerful functions that give you the means to expand numbers in programs in a fairly portable way. You now have enough information to make effective use of them in your programs. Good luck!
Table 1: MP functions versus version of MP.
UNIX UNIX Research UNIX 4.3 GMP
V7 Tenth Edition SVR4 BSD 1.3.2
header file <mp.h> no yes yes yes yes
type MINT no no yes yes yes
type mint no yes no no no
fmin no yes no yes1 no
fmout no yes no yes no
gcd yes2 no8 yes yes yes
invert no no no yes no
itom yes5 yes9 yes yes5 yes
madd yes yes yes yes yes
mcmp no yes7 yes7 yes7 yes
mdiv yes yes yes yes yes17,19
mfree no yes yes no yes
mgcd no15 yes no15 no15 no15
min yes1 yes yes1 yes1 yes1,20
m_in no no no yes1 no
mout yes yes yes yes yes19
m_out no no no yes14 no
move no yes no yes yes
mpow no13 yes no13 no13 no13
msqrt yes yes yes yes yes
msub yes yes yes yes yes
mtox no no yes16 no yes18
mult yes yes yes yes yes
omin no no no yes1 no
omout no no no yes no
pow yes no10 yes yes yes
rpow yes3 yes6 yes11 yes12 yes11
sdiv yes4 yes9 yes yes yes17
xtom no no yes16 no yes
1 fmin, min, m_in, omin: Return type and/or value not documented.
2 gcd: Prototype is given in man page, but not documented.
3 rpow: Second argument (exponent) is documented as "MINT *";
perhaps this is wrong, as all other versions seem to use short or
int.
4 sdiv: Second argument (divisor) type not documented; probably
short as in other versions.
5 itom: Argument type is not documented; probably short.
6 rpow: Second argument (exponent) type not given. Perhaps short or
int was meant.
7 mcmp: Return type not documented; probably int.
8 gcd: Probably conflicts with a libm/math.h routine. In this
version, use mgcd instead.
9 itom, sdiv: Don't call with SHRT_MIN (the most negative short
value).
10 pow: Probably conflicts with pow in libm/math.h; use mpow
instead.
11 rpow: Second argument (power) documented as short.
12 rpow: Second argument (power) documented as int; use short for
portability.
13 mpow: Use pow instead.
14 m_out: Only supports bases 2--10.
15 mgcd: Use gcd instead.
16 mtox, xtom: It is undocumented whether hex uses upper- and/or
lowercase letters.
17 mdiv, sdiv: Remainder is same sign as the dividend.
18 mtox: Generates lowercase hex letters; negative numbers have
leading minus sign.
19 mdiv, mout: Buggy on systems with 16-bit int.
20 min: Ignores leading spaces and tabs.
Example 1: Define NEED_MINT when using the MINT type name.
#include <mp.h> #ifdef NEED_MINT typedef mint MINT; #endif
Example 2: Declaring and initializing variables.
MINT * One; MINT * Three; MINT * Sum; ... One = itom(1); Three = itom(3); Sum = itom(7); /* dummy value */
Table 2: MP I/O routines.
Input Routines Base File Prototype 10 any int? fmin( MINT * Number, FILE * File ) 10 stdin int? min( MINT * Number ) 2--10 any int? m_in( MINT * Number, int Base, FILE * File ) 8 stdin int? omin( MINT * Number ) Output Routines Base File Prototype 10 any void fmout( MINT * Number, FILE * File ) 10 stdout void mout( MINT * Number ) 2--10 any void m_out( MINT * Number, int Base, FILE * File ) 8 stdout void omout( MINT * Number )
Example 3: The madd function.
madd( One, /* first input */ Three, /* second input */ Sum); /* result - set */
Listing One
/* Copyright (c) 1994 by John Rogers. All rights reserved. * FUNCTION - MPHelp.h declares routines which are "helpers" to users of the * multiple precision (MP) library. The routines declared here are: * miseven(number), misodd(number), misprime(number), miszero(number) * All of these functions return "boolean" values in "int" types, where 0 * means false and 1 means true. * AUTHOR - JR (John Rogers), CompuServe: 72634,2402 * Internet: [email protected] */ #ifndef MPHELP_H #define MPHELP_H #include <mp.h> /* MINT or mint typedef. */ /* Define NEED_MINT in makefile if local mp.h does not provide the typedef. * In that case, we define the MINT type here. */ #ifdef NEED_MINT typedef mint MINT; #endif /* ROUTINES, in alphabetical order: */ int /* Note: return values are 0=false and 1=true. */ miseven( const MINT * Number); int /* Note: return values are 0=false and 1=true. */ misodd( const MINT * Number); int /* Note: return values are 0=false and 1=true. */ misprime( const MINT * Number); int /* Note: return values are 0=false and 1=true. */ miszero( const MINT * Number); #endif /* MPHELP_H */
Listing Two
/* Copyright (c) 1994 by John Rogers. All rights reserved. * FUNCTION - MPHelp.c contains the MP "helpers": miseven(number), * misodd(number), misprime(number), miszero(number). All of these functions * return "boolean" values in "int" types, where 0 is false and 1 means true. * AUTHOR - JR (John Rogers), CompuServe: 72634,2402 * Internet: [email protected] */ #include <assert.h> /* assert(). */ #include <mp.h> /* MINT typedef. */ #include "mphelp.h" /* My prototypes. */ /* ROUTINES, in alphabetical order: */ int /* Note: return values are 0=false and 1=true. */ miseven( const MINT * Number) { MINT * Quotient; short Remainder; int ReturnValue; /* 0=false, 1=true. */ /* Initialize MINT variable (any value will do), so MP routines * allocate space. */ Quotient = itom(7); /* Dummy value. */ /* Divide Number by two and look at remainder. */ sdiv( Number, /* dividend - input */ 2, /* divisor - input */ Quotient, /* quotient - output */ & Remainder); /* remainder - output */ /* We'll return "true" if-and-only-if Remainder is zero. */ ReturnValue = ( Remainder == 0 ); mfree( Quotient ); /* Free space used by temp. */ return (ReturnValue); } /* miseven */ int /* Note: return values are 0=false and 1=true. */ misodd( const MINT * Number) { return ( !miseven( Number ) ); } int /* Note: return values are 0=false and 1=true. */ misprime( const MINT * Candidate) { int CompareResult; MINT * Constant_Two = itom(2); MINT * Divisor = NULL; MINT * MaxDivisor = NULL; MINT * Quotient = NULL; MINT * Remainder = NULL; int ReturnValue; /* 0=false, 1=true. */ /* Check for easy cases: * -infinity <= x <= 1 not prime * x = 2 prime * x > 2, x is even not prime */ CompareResult = mcmp(Candidate, Constant_Two); if (CompareResult < 0) { /* Anything less than 2 isn't prime. */ ReturnValue = 0; /* false */ goto Cleanup; } else if (CompareResult == 0) { /* Exactly two, yes that is a prime. */ ReturnValue = 1; /* true */ goto Cleanup; } assert( CompareResult > 0 ); if (miseven(Candidate)) { ReturnValue = 0; /* false */ goto Cleanup; } /* Well, all that's left is the hard stuff. Try all of the odd divisors * from 3 up to the square root of the candidate. */ assert( misodd( Candidate ) ); Divisor = itom( 3 ); MaxDivisor = itom( 1 ); Quotient = itom( 1 ); /* don't care value */ Remainder = itom( 1 ); /* don't care value */ msqrt( Candidate, /* input value */ MaxDivisor, /* square root */ Remainder ); /* remainder */ for ( ; ; ) { /* loop forever */ /* Try dividing this one. */ mdiv( Candidate, /* dividend */ Divisor, Quotient, Remainder); /* Does this divisor divide evenly? */ if ( miszero( Remainder ) ) { /* If we were dividing by Candidate, then Candidate must be prime. */ if (mcmp( Candidate, Divisor ) == 0) { ReturnValue = 1; /* true */ goto Cleanup; } else { /* Otherwise, if this divisor divides evenly, it factors the * candidate, which therefore cannot be prime. */ ReturnValue = 0; /* false */ goto Cleanup; } } /* Have we gone as far as we can? If so, this must be prime! */ if (mcmp( Divisor, MaxDivisor ) >= 0) { ReturnValue = 1; /* true */ goto Cleanup; } /* Bump to next odd divisor. We shouldn't use MP routines to * update Divisor in place, so use a temporary for result of madd(). */ #define RandomTemp Quotient madd( Divisor, /* a value */ Constant_Two, /* 2nd value */ RandomTemp ); /* the sum */ move( RandomTemp, /* source */ Divisor ); /* dest */ assert( misodd( Divisor ) ); } /* loop forever */ Cleanup: /* Free memory used by temp vars. */ if (Constant_Two != NULL) { mfree( Constant_Two ); } if (Divisor != NULL) { mfree( Divisor ); } if (MaxDivisor != NULL) { mfree( MaxDivisor ); } if (Quotient != NULL) { mfree( Quotient ); } if (Remainder != NULL) { mfree( Remainder ); } return (ReturnValue); } /* misprime */ int /* Note: return values are 0=false and 1=true. */ miszero( const MINT * Number) { MINT * Constant_Zero = itom(0); int ReturnValue; /* 0=false, 1=true. */ /* We'll use the standard mcmp (MP compare) function for this. * mcmp(a,b) returns <0 if a<b * =0 if a=b * >0 if a>b */ if ( mcmp( Number, Constant_Zero ) == 0 ) { ReturnValue = 1; /* true */ } else { ReturnValue = 0; /* false */ } mfree( Constant_Zero ); return (ReturnValue); } /* miszero */
Listing Three
/* Copyright (c) 1994 by JR (John Rogers). All rights reserved. * FUNCTION - This program takes a given number and * determines whether or not it is prime. * SYNTAX - IsPrime -n number * AUTHOR - JR (John Rogers) CompuServe: 72634,2402 * Internet: [email protected] */ #include <assert.h> /* assert(). */ #include <mp.h> /* MINT typedef, itom(). */ #include "mphelp.h" /* misprime(). */ #include "sample.h" /* Die(), StringToM(). */ #include <stdio.h> /* printf(). */ #include <stdlib.h> /* EXIT_SUCCESS. */ #include <unistd.h> /* getopt(), optind, optarg. */ #define USAGE \ "This program determines if a number is prime.\n" \ "Usage: IsPrime -n number\n" \ "Or: IsPrime -?\n\n" \ "where:\n" \ " -n number gives number to check\n" \ " -? displays this message\n\n" \ "Author: JR (John Rogers).\n" int main( int argc, char * argv[]) { const char * MyOpts = "n:N:"; MINT * Number = NULL; int ThisOpt; int ValueIsPrime; /* 0=false, 1=true */ /* Do initial setup and argument handling. */ while ((ThisOpt=getopt(argc,argv,MyOpts)) != EOF) { switch (ThisOpt) { case 'n': case 'N': Number = StringToM( optarg ); assert(Number != NULL); break; case '?': Die( USAGE ); /*NOTREACHED*/ default: Die( "bad return value from getopt"); /*NOTREACHED*/ } } /* Handle missing number. */ if (Number == NULL) { Die( USAGE ); /*NOTREACHED*/ } /* Call an MP helper to do the hard work. */ ValueIsPrime = misprime( Number ); /* Tell user what we found out. */ mout( Number ); if (ValueIsPrime) { printf("NUMBER IS PRIME!\n"); } else { printf("NUMBER IS NOT PRIME!\n"); } /* Free temp storage. */ if (Number != NULL) { mfree( Number ); } /* All done; set exit code. */ if (ValueIsPrime) { return(EXIT_SUCCESS); } else { return(EXIT_FAILURE); } } /* main */
Listing Four
/* Copyright (c) 1994 by JR (John Rogers). All rights reserved. * FUNCTION - This program generates a table of powers * of 2, using the MP library routines. * SYNTAX - PowTable -n number * AUTHOR - JR (John Rogers) CompuServe: 72634,2402 * Internet: [email protected] */ #include <assert.h> /* assert(). */ #include <mp.h> /* MINT, itom, rpow, mfree. */ #include "sample.h" /* Die(), StringToShort(). */ #include <stdio.h> /* stdout. */ #include <stdlib.h> /* EXIT_SUCCESS. */ #include <unistd.h> /* getopt(), optind, optarg. */ /* Only define this if m_out is supported at all. */ /* #define MP_SUPPORTS_M_OUT */ /* Only define this if m_out supports bases > 10. */ /* #define MP_OUTPUT_SUPPORTS_LARGE_BASES */ #define MY_USAGE \ "Usage: PowTable -n number\n" \ "Author: JR (John Rogers).\n\n" static void DisplayOneTableEntry( short PowerToCompute ) { MINT * Two = itom(2); MINT * PowerOf2 = itom(42); /* Dummy value. */ /* Compute regular power (2 ** N). */ rpow( Two, /* number to raise */ PowerToCompute, /* exponent (short) */ PowerOf2); /* result */ /* Write the power of two for this one. */ printf( "\n\n2 ** %d is:\n", (int) PowerToCompute ); #ifdef MP_SUPPORTS_M_OUT /* Write 2**N in binary, octal, decimal, hex. */ m_out( PowerOf2, 2, stdout ); m_out( PowerOf2, 8, stdout ); m_out( PowerOf2, 10, stdout ); #ifdef MP_OUTPUT_SUPPORTS_LARGE_BASES m_out( PowerOf2, 16, stdout ); #endif #else /* Output in decimal (only base supported). */ mout( PowerOf2 ); #endif mfree( Two ); /* Free temps */ mfree( PowerOf2 ); } /* DisplayOneTableEntry */ static void GeneratePowerTable( short MaxPower ) /* >= 1 */ { short CurrentPower = 1; for (;;) { DisplayOneTableEntry( CurrentPower ); ++CurrentPower; /* Have we gone far enough? */ if (CurrentPower > MaxPower) { break; /* done; go cleanup and return */ } } } /* GeneratePowerTable */ int main( int argc, char * argv[]) { short MaxPower = 0; const char * MyOpts = "n:N:"; int ThisOpt; /* Do initial setup and argument handling. */ while ((ThisOpt=getopt(argc,argv,MyOpts)) != EOF) { switch (ThisOpt) { case 'n': case 'N': MaxPower = StringToShort( optarg ); assert(MaxPower != 0); break; case '?': Die( MY_USAGE ); /*NOTREACHED*/ default: Die( "bad return value from getopt"); /*NOTREACHED*/ } } /* Handle missing argument. */ if (MaxPower == 0) { Die( MY_USAGE ); /*NOTREACHED*/ } /* Generate the table. */ GeneratePowerTable( MaxPower ); return(EXIT_SUCCESS); } /* main */
Copyright © 1995, Dr. Dobb's Journal
