DI Management Home > BigDigits

# BigDigits multiple-precision arithmetic source code

BigDigits is a free library of multiple-precision arithmetic routines written in ANSI C to carry out large natural number calculations as required in cryptography calculations. This code has been built using the algorithms in Knuth Vol 2 and Menezes as the primary references. New version 2.6 released 31 March 2016: see New in Latest Version below. 2020-04-08: If using on a 64-bit non-Windows platform, check the Errata below.

The BigDigits library includes the classical multiple-precision arithmetic algorithms from Knuth: add, subtract, multiply and divide. It also includes modular multiplication, exponentiation and inversion; bit manipulation; number theory functions such as the greatest common divisor and Jacobi symbol; 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.

Now available: bdcalc is a command-line calculator for large natural numbers. It has all the functionality of the BigDigits library available from the command-line, or from a script file. More details and downloads are available on the bdcalc page.
New version 2.2 of bdcalc released 17 June 2016.

## Contents

Affiliate disclosure: we get a small commission for purchases made through the above links

## Interfaces to the library

The library has two interfaces: a complete set of functions called using the BIGD handle (the "bd" library) which handles memory allocation automatically, and the set of core functions (the "mp" library) which requires a bit more care but is faster and has an option (NO_ALLOCS) to avoid all memory allocation calls completely.

We recommend you use the bd library, which has 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, and another option to make use of native 64-bit integers.

Check for any errata.

## Documentation

This documentation has been produced using the excellent Doxygen tool.

Version 2.6 (Released 31 March 2016)
• Added new functions to perform modular arithmetic operations, particularly useful for doing computations in a prime field for elliptic curves:
• mpModSquare, bdModSquare to compute a = x^2 (mod m)
• mpModSqrt, bdModSqrt to compute a square root modulo a prime p
• mpModHalve, bdModHalve to compute a = x / 2 (mod p)
• mpModAdd, bdModAdd and mpModSub, bdModSub to compute w = u ± v (mod m) quickly for u and v in the range [0, m-1].
• The function mpModSpecial to compute u = v (mod m) in the special case where v is known to be in the range [0, km] for k a small integer. This is faster than mpModulo.
• Added functions mpToShort, bdToShort to convert a multiprecision integer to a single digit.
• Added functions mpShortIsEqual, bdShortIsEqual to check if a multiprecision integer is equal to a single digit.
• Added the explicit constant-time comparison functions.
• mpEqual_ct(), bdIsEqual_ct()
• mpIsZero_ct(), bdIsZero_ct()
• mpCompare_ct(), bdCompare_ct()
• Reversed the change in v2.5 to make the comparison functions constant-time by default. The original comparison functions are now not constant-time, as they were before v2.5. The user must explicitly select the constant-time variant (ending _ct). To summarize: the functions mpEqual, bdIsEqual, mpIsZero, bdIsZero, mpCompare and bdCompare are not constant-time in v2.6.
Comment: We realised after we made the change in v2.5 that the majority of comparison operations in our own code were simple checks for special cases and were not comparing secret values, so they did not need the more expensive constant-time operations. There were only a few occasions where we did comparisons on secret values and these were easily found and explicitly changed to do a constant-time comparison. So, sorry, we shouldn't have changed the default behaviour in v2.5. It's back now to how it was before.
• Added the functions mpCompileTime, bdCompileTime to show the time and date of compilation.
• Added the function mpPrintBits to print the value in 0/1 bit format.
• Added the option to specify the global preprocessor definition, say, MAX_FIXED_BIT_LENGTH=640 when using the NO_ALLOCS option for the "mp" functions. This allows you to change the maximum fixed length of arrays without editing the core source code.
• Increased the size of the fixed-length buffers in the mpPrint*() functions to avoid buffer overflows
• Fixed an issue with the mpPrint*() functions when displaying zero values.
• Added the bdNewVars and bdFreeVars functions to allow "bulk" creation and destruction of "bd" variables in a single line of code.

## Changes in earlier versions

Version 2.5 (Released 22 October 2015)
• Fixed bdShortMult to catch empty u parameter.
• Fixed bdPrintBits to print "0" if bit length is zero.
• Changed comparison functions to be constant-time by default. Retained the old, quicker functions as _q variants. (Changed back in version 2.6).
• Added constant-time modular exponentiation functions bdModExp_ct and mpModExp_ct.
• Added mpPrintDecimalSigned function to display negative numbers ("mp" only).
Version 2.4 (Released 27 April 2013)
• Added markup for Doxygen documentation to produce a manual in CHM and HTML formats.
• Updated the "pretty-good" internal RNG functions in bigdRand and bigdigitsRand to improve the seeding of the ANSI X9.31 PRNG algorithm.
• Fixed the NO_ALLOCS option so it does not use malloc with the mpConv functions (thanks to Kevin Kramb for pointing this out).
• Added the new functions bdRandomOctets, bdRandomNumber and bdQuickRandBits, which do various contortions you might want to do with random numbers.
• Updated the test source code modules.
Version 2.3 (Released 11 November 2011)
• Added new functions to compute the integer cube root: bdCubeRoot and mpCubeRoot.
• Updated and improved the integer square root functions: bdSqrt and mpSqrt.
• Added the more useful print functions with an optional prefix and suffix:
• bdPrintHex prints a bigdigit number in hex format.
• bdPrintDecimal prints a bigdigit number in decimal format.
and the equivalent mp functions mpPrintHex and mpPrintDecimal.
• Added the new function bdPower to compute the n-th power of a number y = gn. Use this with caution as you can quickly run out of memory. It's meant for small exponents like n=3.
• Improved the efficiency of the greatest common divisor functions, bdGcd and mpGcd.
• Fixed the bug in mpChs (thanks to Valery Blazhnov).
• Fixed the bug in mpAlloc which would fail if you called bdSetDigit(b, 0) (thanks to "Radistao").
Version 2.2 (Released 31 July 2008)
• Added sliding-window exponentiation to speed up the ModExp functions.
• bdJacobi computes the Jacobi (Legendre) symbol.
• bdGetBit returns the value of a bit at a given position in a number.
• bdSetBit sets the bit at a given position in a number.
• bdNotBits computes bitwise a = NOT a.
• Added the new functions mpIsNegative, mpChs and mpAbs in anticipation of adding full signed integer functionality in the next a later version (but see A Caution About Negative Numbers).
• Added the NO_ALLOCS option to compile the "mp" library without using any memory allocation.
• Added the USE_64WITH32 option to use the 64-bit integer types (long long) if available on a 32-bit machine.
• Improved the zeroisation and destroy macros.
• Moved the type declarations for the exact-width types to a separate include file.
• Re-organised and re-named the ancillary source and include files (yet again!). See How to Compile.
• Deprecated the bd*Ex functions in favour of re-named "safe" bd*_s functions.
Version 2.1 (Released 23 August 2006)
• Added support for 64-bit compilers.
• Added fudge so opaque pointers will work with C++ compilers.
• bdSqrt computes the `integer' square root, s = floor(sqrt(x)).
• bdModPowerOf2 computes a = a (mod 2L).
• 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, 2n]) 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.
Version 2.0 of BigDigits includes a complete set of "bd" wrapper functions around the original BigDigits "mp" functions. In version 1, we showed a few examples of how memory management could be handled. In version 2.0, all memory management is handled automatically behind opaque pointers. You can now use the bd functions to do all your multiple-precision computations without declaring the size of any array. Of course, you can still use the sub-set of mp functions if you wish. We've added memory allocation where necessary for temporary variables, so there is now no requirement to hard-code any maximum array size as there was in version 1.

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.

As of version 2.5 (October 2015) we have re-issued this source code under the Mozilla Public License, v. 2.0.

We have no objection whatsoever to developers using this code in commercial software or as part of an open source project provided the MPL licence remains unchanged in our source code modules.

[For reference, here is the old copyright and licence notice for version 2.4 (April 2013) and earlier.]

<a href="https://www.di-mgt.com.au/bigdigits.html">BigDigits multiple-precision arithmetic source code</a>

## Using the BigDigits Software

A simple example that multiplies the 48-bit number 0xdeadbeefface (244837814106830) by 2 (we keep this example small so we can reproduce it easily)

#include "bigd.h"
int main(void)
{
BIGD u, v, w;
/* Create new BIGD objects */
u = bdNew();
v = bdNew();
w = bdNew();

/* Compute 2 * 0xdeadbeefface */
bdSetShort(u, 2);
bdMultiply(w, u, v);

/* Display the result in hex and decimal */
bdPrintHex("0x", w, "\n");
bdPrintDecimal("", w, "\n");

/* Free all objects we made */
bdFree(&u);
bdFree(&v);
bdFree(&w);
return 0;
}

This should display

0x1bd5b7ddff59c
489675628213660
which you can verify using the Hex mode of the Windows Calculator in Scientific mode.

All multiple-precision variables need to be declared with a

BIGD v;
statement and created with
v = bdNew();
When finished, free the resources with
bdFree(&v);

Between calling bdNew() and bdFree(), you can use the variable to carry out any multiple-precision operation using any function in the bigd.h interface.

Note the bdFree() function is passed the address of the variable, which is also set to NULL at the same time to avoid accidental reuse.

### A more compact version

Use the new bdNewVars and bdFreeVars functions added in version 2.6 to create and free BIGD objects in bulk (just don't forget the NULL at the end of each list like, er, we did).

#include "bigd.h"
int main(void)
{
BIGD u, v, w;
/* Create new BIGD objects */
bdNewVars(&u, &v, &w, NULL);

/* Compute 2 * 0xdeadbeefface */
bdSetShort(u, 2);
bdMultiply(w, u, v);

/* Display the result in hex and decimal */
bdPrintHex("0x", w, "\n");
bdPrintDecimal("", w, "\n");

/* Free all objects we made */
bdFreeVars(&u, &v, &w, NULL);
return 0;
}

### Remarks

In the bd functions, all internal memory allocation is handled automatically. The mp functions require the user to allocate fixed-length arrays of sufficient length: this is faster but less safe.

This code is all written in ANSI C (except the optional ASM and Win32-specific bits). C++ programmers should be able to incorporate this code into their programs with ease - see Compiling with C++. There is a complete list of functions in the documentation. For more details on compiling, see How to Compile.

CAUTION: there are a few "gotchas" - some variables get trashed before use so be careful using the same variable twice in a function call, see Overlapping Variables and other "Gotchas".

The download contains the full source code files for the multiple-precision library functions, test programs, a Makefile for Linux/cygwin gcc users, and a list of the md5sums for all the files in the library. Check the Errata below for changes.

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.
bigdtypes.h
A common set of macros for exact-width types required by both bigdigits.h and bigd.h. Does not need to be explicity included.
bigdRand.h
A separate interface for the internal random number functions that rely on spRandom. Include this if your code makes use of the bdRandomBits or bdRandDigit functions.
bigdRand.c
Separate source code for the internal RNG functions bdRandomBits or bdRandDigit.
bigdigitsRand.c
Source code for the internal RNG spBetterRand and mpRandomBits functions. 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.
Do not mix these source files with those from older versions of BigDigits.

## How to Compile

### Compiling for the "bd" library

To compile a program that uses the bigd functions without using the internal RNG functions, create a source file, say, bdcode.c:
/* bdcode.c */
#include "bigd.h"
int main()
{
BIGD b;
b = bdNew();
/* ... */
bdFree(&b);
return 0;
}
and compile your source file bdcode.c along with the following files
bigd.c
bigd.h
bigdigits.c
bigdigits.h
bigdtypes.h

If your program requires the internal RNG functions (bdRandDigit or bdRandomBits) then you will need to add

bigdRand.c
bigdRand.h
bigdigitsRand.c
bigdigitsRand.h

### Creating a static library of BIGD functions

We prefer to create a static library of BIGD functions to use in our applications. The only interface then required is the bigd.h include file and, if required, bigdRand.h.

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

bigd.c
bigdigits.c
bigdRand.c
bigdigitsRand.c

To use the faster ASM alternative, define the pre-processor constant USE_SPASM or, alternatively, define USE_64WITH32 to use the native 64-bit integers.

To use the resulting library, include the lines

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

in your source code and link to the bigd.lib library. All the bd functions can then be called from your program.

### Using the "mp" library

To create a program that just uses the mp functions, make a source file mpcode.c:
/* 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;
}
and compile mpcode.c along with the following files
bigdigits.c
bigdigits.h
bigdtypes.h

If you use the "mp" RNG functions spBetterRand or mpRandomBits, then you need to add

bigdigitsRand.c
bigdigitsRand.h

The "mp" functions require you to allocate fixed arrays of digits instead of relying on the automatic allocation routines that the "bd" library uses. Note that some "mp" functions still use `malloc()` to create temporary arrays.

To avoid any use of memory allocation functions, define the NO_ALLOCS preprocessor macro. The "mp" functions will then use fixed arrays and not use any memory allocation functions at all. You can set the value of MAX_FIXED_DIGITS in the bigdigits.h file. It is currently set to cope with integers of up to 8192 bits. If you don't define NO_ALLOCS, then this maximum is ignored and temporary arrays will be allocated as required.

### Preprocessor definitions

In both the "bd" and "mp" libraries you can define one of the preprocessor macros

• USE_SPASM to use the faster ASM routines (on Intel processors where the __asm option is available with your compiler); or
• USE_64WITH32 to use the 64-bit integers if available (e.g. long long)

In the "mp" library only, you can define NO_ALLOCS to avoid using the memory allocation functions completely. If this is set, then you can set the MAX_FIXED_DIGITS definition for the maximum integer size you want to handle (the current default is 8192 bits).

You can detect which of these preprocessor definitions has been applied by using the bdVersion or mpVersion functions (they return the same results).

2300 = none used
2301 = USE_SPASM
2302 = USE_64WITH32
2305 = NO_ALLOCS
2306 = NO_ALLOCS+USE_SPASM
2307 = NO_ALLOCS+USE_64WITH32

The USE_SPASM option takes precedence over USE_64WITH32. The option NO_ALLOCS cannot be used with the "bd" library, only with "mp".

### Summary of files required for compiling

The internal RNG functions are bdRandomBits or bdRandDigit in the "bd" library; and mpRandomBits or spBetterRand in the "mp" library.

LibraryNo internal RNG functionsUses an internal RNG function
"mp" only
bigdigits.c
bigdigits.h
bigdtypes.h
bigdigits.c
bigdigitRand.c
bigdigits.h
bigdtypes.h
bigdigitRand.h
Full "bd"
bigd.c
bigdigits.c
bigd.h
bigdigits.h
bigdtypes.h
bigd.c
bigdigits.c
bigdRand.c
bigdigitRand.c
bigd.h
bigdigits.h
bigdtypes.h
bigdRand.h
bigdigitRand.h

### Print format specifiers

We use BIGD-compatible versions of the C99 format specifiers, so you can replace print expressions of the form

bdigit_t b;
/* ... */
printf("Value=%lx\n", b);
with the more generic
printf("Value=%" PRIxBIGD "\n"), b);
where PRIxBIGD is our own version of the C99 PRIx32.

### Compiling with C++

This code is written in ANSI C which is (almost, but not quite) a sub-set of C++. The opaque pointer techniques in bigd.h and bigd.c based on [HANSON97] rely on the separate namespaces for structures and type names, which fall over when compiled with C++. There is a `fudge' in the code thanks to Ian Tree which fixes this. See the example in t_bdCPP.cpp. Do not try to change the extensions of the core source code files like bigd.c to .cpp. This will not work.

### Compiling in Linux and other Unix environments

See the Makefile in the download for Linux/Cygwin environments using gcc. These are not necessarily examples of the most efficient Makefiles possible.

### Compiling on a Mac

2012-01-23: Arno Hautala reports success compiling and running the included test suite on Mac OS X.

I did have to make one addition to line 36 of bigdtypes.h to add " || defined(__APPLE__)" so that stdint.h would be properly included.

Update: This change has been incorporated in version 2.4.

### Compiling on a 64-bit machine

The code has been changed to use the C99 <stdint.h> types uint32_t and uint16_t instead of "unsigned long" and "unsigned short". These are implemented via some convoluted preprocessor instructions that seem to work on most test systems we've tried. You may need to tweak these. Thanks to Terry R. Friedrichsen and Allen Zhao for their help on this.

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.

### Using with an 8-bit microcontroller

Alfredo Ortega has modified the BigDigits code to run on an 8-bit microprocessor. See Implementation of the RSA algorithm in a 8-bit microcontroller.

Alternative: try compiling the "mp" library with the NO_ALLOCS preprocessor macro defined. This removes all calls to memory allocation functions and just uses fixed arrays. You can define the maximum length of array you want to use.

### Adding PKCS#1 and elliptic curve cryptography on binary and prime fields

2009-12-10: Dang Nguyen Duc has extended an earlier version of the BigDigits library to a more complete library for implementing public key cryptography. In particular, adding PKCS#1, elliptic curve on binary and prime fields. He has posted his code to libmcrypto: Math Library for Crypto Students.

“ Thank you very much for writing very clean and easy-to-follow code. That was much more useful for my study than using GNU MP. Along the way, I have added a few functions to your bigdigits library... I still attempt to maintain the cleanness and easy-to-follow property of your original bigdigit library without optimization or inline assembly. A study tool is my top priority. ”

## Compatibility

The source code and tests have been compiled and verified on Windows 7/XP/Win98 operating systems using MSVC++ versions 5, 6, 7, 8 (Express 2005) and 9 (MSVC++ 2008); Borland C++ 5.5.1 for Win32; Cygwin using gcc version 3.3.3 on W2K; and Linux (Fedora Core 1) using gcc version 3.3.2 (Red Hat Linux 3.3.2-1).

Allen Zhao reports success compiling version 2.1 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. Arno Hautala reports success (23 January 2012) compiling on Mac OS X.

## Test Code

Test source file names start with "t_". They test various aspects of the code with perhaps somewhat boring tests.
t_bdSimple.c
A simple example that computes 2 * 0xdeadbeefface (244837814106830) and displays the result.
t_bdTest.c
A selection of tests that tests most functions 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. (We use 1025 bits to test the generator for odd numbers.)
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. Again, this just tests the functions. You could re-write this to use a fixed RSA key pair.
t_bdRsaCrack.c
Code to "crack" RSA if the user has done encryption in a certain (poor) way.
t_bdRsaFactorN.c
How to factor the RSA modulus given the secret exponent d.
t_bdRSA_blinded.c
Carry out RSA calculations using constant-time exponentiation and blinding to provide resistance against simple power attacks (SPA) and differential power analysis (DPA).
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_bdRandomOctets.c
Examples using the bdRandomOctets and bdRandomNumber functions.
t_bdQuickRandBits.c
Examples comparing the bdQuickRandBits and bdRandomBits functions.
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_mpJacobi.c
Carries out some tests of the mpJacobi function using an example taken from ANSI X9.31 Appendix D.5.2.
t_spExtras.c
Carries out some tests on the functions in spExtras.c.
t_bdCPP.cpp

## Overlapping Variables and other "Gotchas"

The mp functions are rather raw and user-unfriendly, but faster. They assume that the user has allocated memory of the correct length and understands the quirks of the various functions. For example, mpMultiply requires the two multiplicands x and y to be of the same ndigits length but the product p must be twice that length. The user is meant to know and understand this or face the consequences of a GPF error when there is overflow.

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 bdAdd_s() and bdMultiply_s() 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.

## A Caution About Negative Numbers

The BigDigits library is designed to work with the natural numbers ℕ; that is, the non-negative integers 0,1,2,....

Negative numbers will "work" with the fixed-array "mp" functions if you are happy using twos-complement arithmetic. There are a few "experimental" functions included (mpIsNegative(), mpChs() and mpAbs()) intended for signed integers, but so far we have not developed that very far. We added the function mpPrintDecimalSigned() in v2.5, but only for the fixed-length "mp" functions.

But be careful with the memory-managed "bd" functions. If you go below zero with the "bd" functions the results are undefined. You should either make sure you always do cut-off subtraction that yields a non-negative result (i.e. to compute a - b, first check that ba; if not return zero or flag an error) or check the "carry/overflow" result from the bdSubtract or bdDecrement operation and flag an error if this is not zero.

We give an example of the problem in the test code t_bdTests.c
/* Cope, sort of, with `negative' numbers */
printf("Go ``negative''...\n");
bdSetZero(u);
overflow = bdShortSub(w, u, 2);
pr_msg("w=0-2=", w);
printf("overflow=%" PRIuBIGD "\n", overflow);
overflow = bdIncrement(w);
pr_msg("w+1=", w);
printf("overflow=%" PRIuBIGD "\n", overflow);
overflow = bdIncrement(w);
pr_msg("(Note error!) w+1=", w);
printf("overflow=%" PRIuBIGD "\n", overflow);
which gives the output
Go ``negative''...
w=0-2=fffffffe
overflow=1
w+1=ffffffff
overflow=0
(Note error!) w+1=00000001 00000000
overflow=1

How you could do cut-off subtraction so the result is never negative. Check the return value to see if you underflowed.

int cutoff_subtraction(BIGD w, BIGD u, BIGD v)
{
int carry = 0;
if (bdCompare(u, v) < 0) {
bdSetZero(w);
carry = 1;
}
else {
bdSubtract(w, u, v);
}
return carry;
}

## Random Number Functions

There are three occasions when you might need to use a random number generator in connection with the BIGDIGITs package.

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)) or better, srand(time(NULL) ^ clock()) 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 method 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, both in bigdRand.h. 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

By David Ireland
I know there are several packages out there that do this already. I wanted to write my own set partly to prove to myself I could do it, and partly so we could have a set we could use in our own commercial applications without restriction. You are welcome to use it in your applications, including commercial ones, provided you adhere to the copyright conditions.

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 - look!). 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 code can still be compiled using only a "pure" ANSI C compiler. We've added a few Win32-specific options to handle error messages and to get an accurate 64-bit time, but these are ignored if the WIN32 compiler flag is not defined.

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. However, adding sliding-window exponentiation in version 2.2 made some significant improvements. Adding an mpSquare function 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. In version 2.2 you can alternatively use the USE_64WITH32 preprocessor macro to make use of 64-bit types like long long if available.

We used Knuth Vol 2 [KNUTH] and Menezes [MENEZES] as the primary references plus the other documents noted in the references. All functions in BigDigits were derived from first principles using the algorithms described in the primary references. We are advised that this code complies with the definition of an original work for the purposes of copyright.

Given that this code has been in the public domain for years without any significant bugs being reported, we are very reluctant to make big changes to the core "mp" code. Some points of interest on Version 2:-
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.
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=232 is u = (7fffffff 80000001 00000000 00000000)32, v = (80000000 80000002 00000005)32. See the solution to 4.3.1 Example 22 in Knuth, p625. Fortunately the values given by DEK for b=216 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 15 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

Thank you in particular to the following people for their suggestions and advice that have contributed to the development of this code:
• Jesse Chisholm
• Sandy Dunlop
• Terry R. Friedrichsen
• Peter Hebert
• James Ku
• Ian Tree
• Rickard Westman
• Allen Zhao
• Valery Blazhnov
• Kevin Kramb

## References

Sources of the algorithms used in BigDigits.

## Bibliography

Other sources we used in devising tests and other techniques used in the code.

## Errata

Last updated: 8 April 2020

### 1. Typo in bigdigits.c for uint64_t

bigdigits.c line 193: change HAVE_C99INCLUDES to HAVE_C99_INCLUDES.
190 /* Make sure we have a uint64_t available */
191 #if defined (_WIN32) || defined(WIN32)
192 typedef unsigned __int64 uint64_t;
193 #elif !defined(HAVE_C99INCLUDES) && !defined(HAVE_SYS_TYPES)
194 typedef unsigned long long int uint64_t;
195 #endif
In fact, the lines 190 to 195 are redundant and can be removed.

### 2. Error in t_bdTest.c

Add a NULL to the end of the lists in lines 167 and 1052.

Before:

167   bdNewVars(&u, &v, &w, &p, &q, &r, &x, &y);   // WRONG!
...
1052  bdFreeVars(&u, &v, &w, &p, &q, &r, &x, &y);  // WRONG!
After:
167   bdNewVars(&u, &v, &w, &p, &q, &r, &x, &y, NULL);
...
1052  bdFreeVars(&u, &v, &w, &p, &q, &r, &x, &y, NULL);