BigDigits multiple-precision arithmetic source code

This library includes the classical multiple-precision arithmetic algorithms from Knuth: add, subtract, multiply and divide. It also includes modular multiplication, exponentiation and inversion; number theory functions such as greatest common divisor; the Rabin-Miller Probabilistic Primality Test procedure from FIPS-186 and ANSI X9.42 to show that a large integer is probably prime; and other utilities to manipulate and handle multiple-precision numbers.

The new version 2 includes a complete set of wrapper functions called using the BIGD handle (the "bd" functions). These are all the functions you should need. Memory allocation is handled automatically, except the initial creation and final release of resources. There is an option to use Intel Assembler to speed up the critical single-digit multiply and divide operations.

Contents

• New in Latest Version
• Changes in Version 2.0
• Copyright and commercial use
• Using the BigDigits Software
• Download
• Files in the Download
• How to Compile
  • Summary of files required for "bd" functions
  • Creating a static library of BIGD functions
  • Using the "mp" functions
  • Print format specifiers
  • Compiling with C++
  • Compiling in Linux and other environments
  • Compiling on a 64-bit machine
  • Using with an 8-bit microcontroller
• Compatibility
• Dependencies
• Test functions
• Overlapping Variables and other "Gotchas"
• Random Number Functions
• History
• Programmer's Comments
• Acknowledgements
• References
• Bibliography
• Function list
• Errata
• Contact
• Register

New in Latest Version

• Added support for 64-bit compilers.
• Added fudge so opaque pointers will work with C++ compilers.
• Added new functions:
  • bdSqrt computes the `integer' square root, s = floor(sqrt(x)).
  • bdModPowerOf2 computes a = a (mod 2^L).
  • bdXorBits computes bitwise a = b XOR c.
  • bdOrBits computes bitwise a = b OR c.
  • bdAndBits computes bitwise a = b AND c.
  • bdVersion returns version number as an integer = major*1000+minor*100+release*10+uses_asm(0|1).
  • bdRandomBits generates a random BIGD number up to n bits long (i.e. r = [0, 2ⁿ]) using the internal RNG.
• Removed the shift limit on bdShiftLeft and bdShiftRight. This was previously limited to a maximum of 32 bits. It will now shift by any value - well, any value up to SIZE_MAX.
• Re-organised the source and include files for the internal random number generator functions to avoid having to include the spRandom modules when they are not needed. See How to Compile.
• Fixed up inconsistencies in the documentation for the bdConv* functions.
• Made minor change to mpIsPrime to comply with ANSI X.42-2003 Annex F.1.1 Range of bases in Miller-Rabin test.
• Made minor improvements to the spBetterRand function.
• Fixed bugs in bdConvertHex, bdGetBit and bdSetBit.
• Made various minor changes in the code to avoid warning messages when compiling with the -Wall -pedantic options.

Changes in Version 2.0

Some speed enhancements have been introduced, including an option to use assembler calls for the critical multiply and divide operations, which makes an appreciable difference in speed on Intel processors. For more details see History. We've also added a more robust random number generator for doing tests.

Copyright and Commercial Use

Our intention is to help developers write useful applications for end users. If you can use this code in your applications to make a buck or a quid out of doing that, good on you. We do, however, retain all rights in the source code, including the rights to distribute or post the source code.

You are welcome to use compiled versions of this code as part of your own executable files and to distribute unlimited copies of such executable files for any purposes including commercial ones. All we ask is that you keep the copyright notices intact in the source code and ensure that the following characters remain in any object or executable files you distribute and are displayed clearly in any accompanying documentation:

Contains multiple-precision arithmetic code originally written by David Ireland, copyright (c) 2001-6 by D.I. Management Services Pty Limited <www.di-mgt.com.au>, and is used with permission.

<A href="http://www.di-mgt.com.au/crypto.html">Cryptography Software Code</a>

Please forward any comments or bug reports to <contact>. The latest version of the source code can be downloaded from www.di-mgt.com.au/bigdigits.html or here. To receive notification of updates (and nothing else!) please register your copy.

Using the BigDigits Software

#include "bigd.h"
int main()
{
	BIGD u, v, w;
	/* Create new BIGD objects */
	u = bdNew();
	v = bdNew();
	w = bdNew();
	bdSetShort(u, 2);
	bdConvFromHex(v, "deadbeefface");
	bdMultiply(w, u, v);
	bdPrint(w, BD_PRINT_NL);
	/* Free all objects we made */
	bdFree(&u);
	bdFree(&v);
	bdFree(&w);
	return 0;
}

0001bd5b 7ddff59c

All multiple-precision variables need to be declared with a

   BIGD v;

   v = bdNew();

   bdFree(&v);

Download

Files in the Download

bigd.h

Include file for BIGD functions, i.e. the "bd" functions with full memory management. The only external interface you should need. Exports the BIGD handle and the bdigit_t typedef for a single digit.

bigd.c

Source code for bd functions. This is a wrapper around all the mp functions in bigdigits.c. It completely hides all references to the mp functions and the DIGIT_T typedef.

bigdigits.h

Include file for the mp functions. Exports the DIGIT_T typedef.

bigdigits.c

Source code for the mp functions.

spBigdigits.c

The core ANSI spMultiply and spDivide functions called by just about every function in this package. These functions can be replaced by the faster assembler versions in spASM.c.

spASM.c

Assembly versions of the core spMultiply and spDivide functions. If your system supports calls to Intel MASM functions, using this instead of spBigdigits will result in significant speed improvements. Requires the preprocessor flag USE_SPASM to be set.

bigdRand.h

A separate interface for the internal random number functions that rely on spRandom. Include this if your code makes use of bdRandomBits or bdRandDigit.

bigdRand.c

Separate source code for the internal RNG functions bdRandomBits or bdRandDigit.

spRandom.c

Source code for the internal RNG spBetterRand function. Called by the functions in bigdRand.c(). We keep this separate so you can substitute your own "best" random function to taste, should you wish.

spExtra.c/.h

Single-digit versions of some of the multiple-precision functions. Not required for either the BIGD or BIGDIGITS libraries but are included for interest. Useful as an introduction to understand how the mp functions work.

t_*.c

Source code for various test functions. See Test Functions.

How to Compile

(NOTE: This has changed in Version 2.1.)

/* bdcode.c */
#include "bigd.h"
int main()
{
	BIGD b;
	b = bdNew();
	/* ... */
	bdFree(&b);
	return 0;
}

bigd.c (+.h)
bigdigits.c (+.h)
spBigdigits.c

Alternatively, if your system supports Intel assembly-language statements and the __asm keyword, then replace the file spBigdigits.c with spASM.c and define the pre-processor constant USE_SPASM. This works in MSVC and results in quite significant speed improvements.

Summary of files required for "bd" functions

The internal RNG functions are bdRandomBits or bdRandDigit.

bigd.h
bigdigits.h

bigd.h
bigdigits.h
bigdRand.h
spRandom.h

bigd.c
bigdigits.c
spBigdigits.c OR spASM.c
yourSourceFile.c

bigd.c
bigdigits.c
spBigdigits.c OR spASM.c
yourSourceFile.c
bigdRand.c
spRandom.c

Creating a static library of BIGD functions

In MSVC, compile the following files to create the static library bigd.lib:

bigd.c
bigdigits.c
bigdRand.c
spRandom.c
spBigdigits.c

bigd.c
bigdigits.c
bigdRand.c
spRandom.c
spASM.c

#include "bigd.h
#include "bigdRand.h"

Using the "mp" functions

/* mpcode.c */
#include "bigdigits.h"
int main()
{
	DIGIT_T u[10], v[10], w[20];
	/* set u and v ... */
	mpMultiply(w, u, v, 10);
	/* ... */
	return 0;
}

bigdigits.c
spBigdigits.c
spRandom.c (only if calling spBetterRand)

Print format specifiers

bdigit_t b;
/* ... */
printf("Value=%lx\n", b);

printf("Value=%" PRIxBIGD "\n"), b);

Compiling with C++

"Render unto C++ the things which are C++'s, and unto ANSI-C the things that are ANSI-C's"

Compiling in Linux and other environments

Compiling on a 64-bit machine

We'd be interested to hear from anyone who has changed the code to use 64-bit digits directly. That is, changed all the uint32_t's to uint64_t, and uint16_t's to uint32_t; changed BITS_PER_DIGIT to 64, MAX_DIGIT to UINT64_MAX, PRIu32 to PRIu64, HIBITMASK to 0x8000 0000 0000 0000 UL, and so forth. One day when we can afford to buy a 64-bit machine, we'll test this ourselves.

Using with an 8-bit microcontroller

Compatibility

Allen Zhao reports success compiling in both 32- and 64-bit modes on Solaris 8, HPUX 11i/PA-RISC; AIX 4.3.3 (32bit) and AIX 5.1(32/64); IRIX 6.5 (32bit) and IRIX64 6.5 (32/64); RH Linux 7.2 / RHEL 3 / RHEL 4 with both gcc and Intel icc (v9.1) for x86/x86_64/ia64 (any combination of 32-bit/64-bit compilation). Thanks Allen. Gustavo del Dago reports success (21 Nov 2006) compiling on the OS/400 platform.

Dependencies

# bigd.c [1]
|
+--# bigdigits.c [2]
|  |
|  +-- spBigdigits.c
|      OR
|      spASM.c [3]
|
+--# bigdRand.c [4]
   |
   +-- spRandom.c [5]

[1] bigd library

exposes the "bd" functions, the opaque BIGD handle, and bdigit_t typedef. Includes automatic memory allocation. A wrapper around the bigdigits library that completely hides all the "mp" functions and the DIGIT_T typedef.

[2] bigdigits library

exposes the fixed-length "mp" functions and the DIGIT_T typedef. Can be used independently of the bigd library. Requires the programmer to declare the lengths of all variables and ensure they are all the correct size.

[3] spASM functions

Either (a) compile with the ANSI C-compatible source file spBigdigits.c or (b) compile with the Intel Assembler-compatible source file spASM.c and define the global preprocessor constant USE_SPASM.

[4] bigdRand functions

An interface to the optional internal RNG functions bdRandomBits and bdRandDigit. Compile the source file if required.

[5] spRandom functions

Provides the core "internal" random number function spBetterRand which uses a not-quite-cryptographically-secure variant of the ANSI X9.31 algorithm. Compile the source file if required.

Test Functions

t_bdSimple.c

A simple example that computes 2 * 0xdeadbeefface and displays the result.

t_bdTest.c

A selection of tests that (should) test every function in the bd library at least once. Uses assert checking to catch errors.

t_bdRSA.c

Generates a new 1025-bit RSA key pair and tests both the RSA encryption and signing operations on random data blocks. It carries out decryption using both the inversion and the CRT methods and compares the time taken. The example callback RNG function shown is a bit of hack using the insecure stdlib rand() function.

t_bdRSA1.c

A variant on t_bdRSA.c. This example takes a short message from the command line and encrypts it in the RSA encryption block using PKCS1-v1.5 encoding and a freshly-generated 1024-bit RSA key. The callback RNG function uses the new bdRandomBits function.

t_bdDSA.c

Reproduces the test vector of the Digital Signature Algorithm (DSA) from Appendix 5 FIPS PUB 186-2 "Digital Signature Standard (DSS)", 27 January 2000.

t_bdRDSA.c

Reproduces some test vectors of the rDSA algorithm from ANSI X9.31-1998 "Digital Signatures Using Reversible Public Key Cryptography for the Financial Services Industry (rDSA)", September 9, 1998. Appendix D.1 Odd e = 3 with 1024-bit n and D.3 Odd e =3 with 2048-bit n.

t_mpTest.c

Carries out some tests on the mp library functions.

t_mpRSA508.c

Uses the mp library functions to reproduce the example of RSA encryption and signature using 508-bit test vectors from "Some Examples of the PKCS Standards", An RSA Laboratories Technical Note, Burton S. Kaliski Jr., 1 November 1993.

t_spExtras.c

Carries out some tests on the functions in spExtras.c.

t_bdCPP.cpp

A simple example in C++. See also Compiling with C++.

Overlapping Variables and other "Gotchas"

For speed and efficiency reasons, some functions require that certain source parameters do not overlap with the destination parameter. The user is also meant to remember this and know the safe combinations. There are assert checks included in the mp source code to trap such errors.

The bd functions are designed to handle the memory allocation issues behind the scenes but the overlap issues are still present for some functions for speed reasons. For example, bdAdd(x, x, y) is OK, but bdAdd(y, x, y) is not, and bdMultiply(y, x, y) is OK, but bdMultiply(x, x, y) is not (we took inspiration in consistency from certain C library and Windows API functions here). If in doubt, use separate variables or use the bdAddEx() and bdMultiplyEx() functions which are safe but slower (or read the instructions more carefully!). Think of this as similar to the difference between the memcpy and memmove functions in the standard C string library.

Random Number Functions

1. To generate a few pseudo-random numbers in a limited range to carry out a few tests.
2. To generate a large set of random digits to use in a large number of tests.
3. To generate a cryptographically-secure random number in a security application.

A few pseudo-random numbers for simple tests

In the first case, you can use the standard ANSI rand() function from the stdlib library and seed it with srand(time(NULL)). You need to be careful if you want a number greater than RAND_MAX, which can be as small as 0x7fff on some systems. We have used this in the bdSetRandTest function and the spSimpleRand(DIGIT_T lower, DIGIT_T upper) function which uses rand() to generate a random number byte-by-byte in the range between `lower' and `upper'. This is adequate for the odd few random numbers in tests that don't need to be secure.

To fill up many arrays of digits with random values for tests

In the second case, using rand() is not suitable because it starts repeating itself too soon. For this purpose we have added the spBetterRand() function, which is based on not-quite-cryptographically-secure techniques. This is called from the BIGD library via the synonym bdRandDigit() and is used by bdRandomBits. This function is perfectly adequate for generating large numbers of random digits to use in tests with an infinitesimally small chance of repetition or non-random behaviour. The output should pass all FIPS-140-2 tests for randomness. The function does, however, rely on static variables and so should be used with care in a multi-threaded environment. We have kept the code for this function in a separate module so it can easily be replaced if you wish.

Cryptographically-secure RNGs

In the third case, you would need to call a cryptographically-secure function. In the bdRandomSeeded() and bdGeneratePrime() functions - which are intended to be used in secure cryptography packages - you pass a reference to a separate RNG call-back function which can be as cryptographically-secure as you wish. Using a separate function enables you to handle multi-threaded issues and use the secure RNG of your choice. We have included a couple of example functions (which aren't secure themselves) in the tests for generating RSA keys (remember we're demonstrating a multiple-precision package, not a cryptographic module). The function has a required format BD_RANDFUNC with parameters to pass the output buffer, required size and seed, which your own function would need to match.

You are welcome to use the included random number functions at your own pleasure and risk. We don't want to enter protracted discussions about how "random" these functions are, but we will respond to obvious errors. Please carry out your own tests to verify the functions provided or substitute your own ones that you trust.

History

I wrote the original BigDigits Version 1 back in 2001. I did most of the research and drafting while sitting overnight with my son who was 8 years old and in hospital at the time (he's fine now). The original design goals were:-

1. Use only "pure" ANSI C
2. Keep it independent of machine endianness
3. Make sure it's reasonably understandable as stand-alone code
4. Avoid having to declare temporary storage where possible

The original requirement to avoid allocating memory was not practical, though. Even worse, having a fixed maximum size constant hidden away in the include file was a cause of problems for some users. So we've added dynamic memory allocation for those "mp" functions that needed temporary variables of varying length, and then for the "bd" functions we added full internal memory handling. We've tested these functions in a variety of situations for memory leakage. Provided you use the bdFree() before you exit, there should be no memory leaks.

We experimented with using Montgomery Exponentiation for an improved ModExp function, but we ended up up with little or no improvement and some very hairy code, so we dropped it. We did add an mpSquare function which gives about a 40% speed improvement against using multiplication, but again at the cost of some complicated code. Using this improves the ModExp function overall by about 10% on the average.

The one experiment that did produce a significant increase in speed was to use Assembler for the frequently-called single-digit multiply and divide functions. This resulted in up to a 40% speed improvement in tests. The original ANSI C version is still the default.

Programmer's Comments

malloc

Using malloc to allocate memory for temporary variables makes for an easier life, but causes a performance hit for those functions that call them repeatedly. Note how we have used our own internal functions like mpModuloTemp to re-use the same memory for temporary variables for the memory-intensive mpModExp function.

Overlapping variables

We agonised for a long time whether to have bdFunction as the "safe" version with bdFunctionEx as the faster version for experts, or vice versa. In the end, we realised that so much of our own code that uses earlier versions of the library would have to be changed that we decided to use the opposite and possibly more confusing convention. (Postscript: With hindsight, perhaps we should have used bdFunction_s.)

Assertion checking

The "bd" functions do assertion checks that a NULL BIGD parameter has not been passed. The "mp" functions do assertion checks to ensure that overlap problems will not occur. We recommend that you leave assertion checking turned on, even in your production code.

mpDivide

Writing the Division algorithm (Algorithm D, Knuth, p272) in elegant C code was tricky. This is probably the most tested part of the library. A suitable example to test the "Add Back" at step D6 of the algorithm for b=2³² is u = (7fffffff 80000001 00000000 00000000)₃₂, v = (80000000 80000002 00000005)₃₂. See the solution to 4.3.1 Example 22 in Knuth, p625. Fortunately the values given by DEK for b=2¹⁶ readily convert upwards for 32-bits.

spDivide

We are pretty certain that the "Add Back" step will never happen in our single-precision ANSI C spDivide function but we are not brave enough to remove the code that handles it. We think this requires at least 3 digits in the divisor, v, and we only have two. At worst there are 3 times 5 unnecessary machine instructions; at best it catches something that 'theoretically' shouldn't happen. Comments - and proofs - would be welcome.

ASM

Using assembler to handle the spDivide and spMultiply functions on a 32-bit machine produced a significant improvement in performance (as you would expect when you get the processor to use its own instruction instead of your own arcane version in ANSI C). You still need to be careful when doing unsigned division in assembler as the DIV operation will throw a #DE exception if there is overflow (i.e. when the resulting quotient is larger than 32-bits). See how we handle it in the code in spASM.c. This exception crops up very rarely in normal operations - typical tests could go for hours without running into it. Why doesn't Intel just set an overflow flag instead of throwing an exception?

Acknowledgements

• Jesse Chisholm
• Sandy Dunlop
• Terry R. Friedrichsen
• Peter Hebert
• James Ku
• Ian Tree
• Rickard Westman
• Allen Zhao

References

• [KNUTH] Donald E. Knuth, The Art of Computer Programming, Vol 2 Seminumerical Algorithms, 3rd ed, Addison-Wesley, 1998.
• [MENEZES] Alfred J. Menezes, Paul C. van Oorschot, Scott A Vanstone, Handbook of Applied Cryptography, CRC Press, 1997, <http://www.cacr.math.uwaterloo.ca/hac/>.
• [SCHNEIER] Bruce Schneier, Applied Cryptography, 2nd ed, Wiley, 1996.
• [CLARK] Andrew Clark, Implementation Issues in Cryptography, <http://sky.fit/qut.edu.au/~aclark/itn536/>, (accessed March 2001).
• [GOURDON] Xavier Gourdon and Pascal Sebah, Arbitrary precision computation, <http://numbers.computation.free.fr/Constants/constants.html>, (accessed March 2001).
• [FIPS-186-2] Digital Signature Standard (DSS), FIPS PUB 186-2, U.S. Department of Commerce/National Institute of Standards and Technology, 2000.
• [ANSIX942] ANSI X9.42-2003 Public Key Cryptography for the Financial Services Industry: Agreement of Symmetric Keys Using Discrete Logarithm Cryptography, American National Standards Institute, 2003.

Bibliography

• [ANSIX931] ANSI X9.31-1998 Digital Signatures using Reversible Public Key Cryptography for the Financial Services Industry (rDSA), American National Standards Institute, 1998.
• [FERGUSON03] Niels Ferguson and Bruce Schneier, Practical Cryptography, John Wiley, 2003.
• [HANSON97] David R. Hanson, C Interfaces and Implementations: Techniques for Creating Reusable Software, Addison-Wesley, 1997, <http://www.cs.princeton.edu/software/cii/>.
• [INTEL2] The Intel Architecture Software Developer's Manual, Volume 2: Instruction Set Reference, Intel Corporation, 1999.
• [ISO99] ISO/IEC 9899:1999 Programming lanuages - C, 2nd ed., ISO/IEC, 1999.
• [PKCS1] PKCS #1, RSA Cryptography Standard, RSA Laboratories, Version 2.1, June 2002.
• [PKCS-EX] Burton S Kaliski, Jr, Some Examples of the PKCS Standards, RSA Laboratories Technical Note, Nov 1993.

Function List

Functions in bigd.h:

bdAdd
bdAddEx
bdAndBits
bdBitLength
bdCompare
bdConvFromDecimal
bdConvFromHex
bdConvFromOctets
bdConvToDecimal
bdConvToHex
bdConvToOctets
bdDecrement
bdDivide
bdDivideEx
bdFree
bdGcd
bdGeneratePrime
bdGetBit
bdIncrement
bdIsEqual
bdIsEven
bdIsOdd
bdIsPrime
bdIsZero
bdModExp
bdModInv
bdModMult
bdModPowerOf2
bdModulo
bdModuloEx
bdMultiply
bdMultiplyEx
bdNew
bdOrBits
bdPrint
bdRabinMiller
bdRandomSeeded
bdSetBit
bdSetEqual
bdSetRandTest
bdSetShort
bdSetZero
bdShiftLeft
bdShiftRight
bdShortAdd
bdShortCmp
bdShortDiv
bdShortMod
bdShortMult
bdShortSub
bdSizeof
bdSqrt
bdSquare
bdSquareEx
bdSubtract
bdSubtractEx
bdVersion
bdXorBits

Functions in bigdRand.h:

bdRandDigit
bdRandomBits

Errata

There are no errata for version 2.1.

Contact

To comment on this Contact DI Management. You can register your copy of BigDigits to receive notification of updates.    [Top]