BigDigits multipleprecision arithmetic source code
BigDigits is a free library of multipleprecision 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.
The BigDigits library includes the classical multipleprecision 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 RabinMiller Probabilistic Primality Test procedure from FIPS186 and ANSI X9.42 to show that a large integer is probably prime; and other utilities to manipulate and handle multipleprecision numbers.
New version 2.2 of bdcalc released 17 June 2016.
Contents
Recommended reading
 Art of Computer Programming, Volume 2: Seminumerical Algorithms by Donald E. Knuth
 Handbook of Applied Cryptography by Alfred J. Menezes, Paul C. van Oorschot and Scott A. Vanstone
 A Course in Computational Algebraic Number Theory by Henri Cohen
 C Interfaces and Implementations by David R. Hanson
 Elliptic Curves in Cryptography by Ian F. Blake, Gadiel Seroussi and Nigel P. Smart
Join Amazon Student FREE TwoDay Shipping for College Students
Affiliate disclosure: we get a small commission for purchases made through the above links
 Interfaces to the library
 Documentation
 New in the latest version
 Changes in earlier versions
 Copyright and commercial use
 Using the BigDigits Software
 Download
 Files in the Download
 How to Compile
 Compiling for the "bd" library
 Creating a static library of BIGD functions
 Using the "mp" functions
 Preprocessor definitions
 Summary of files required for compiling
 Print format specifiers
 Compiling with C++
 Compiling in Linux and other Unix environments
 Compiling on a Mac
 Compiling on a 64bit machine
 Using with an 8bit microcontroller
 Compatibility
 Test code
 Overlapping Variables and other "Gotchas"
 A Caution About Negative Numbers
 Random Number Functions
 History
 Programmer's Comments
 Acknowledgements
 References
 Bibliography
 Errata
 Contact
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 singledigit multiply and divide operations, and another option to make use of native 64bit integers.
Check for any errata.
Documentation
 Documention in CHM format (zipped, 117 kB)
 Documention in HTML format (follow the "Files" link):
 bigd.h  the "bd" functions
 bigdigits.h  the "mp" functions
 bigdRand.h  "bd" random number functions
 bigdigitsRand.h  "mp" random number functions
This documentation has been produced using the excellent Doxygen tool.
New in the latest version
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 computea = x^2 (mod m)
mpModSqrt
,bdModSqrt
to compute a square root modulo a prime pmpModHalve
,bdModHalve
to computea = x / 2 (mod p)

mpModAdd
,bdModAdd
andmpModSub
,bdModSub
to computew = u ± v (mod m)
quickly for u and v in the range[0, m1]
.  The function
mpModSpecial
to computeu = 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 thanmpModulo
.
 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 constanttime 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 constanttime by default.
The original comparison functions are now not constanttime, as they were before v2.5. The user must explicitly select
the constanttime variant (ending _ct).
To summarize: the functions
mpEqual
,bdIsEqual
,mpIsZero
,bdIsZero
,mpCompare
andbdCompare
are not constanttime 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 constanttime 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 constanttime 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 theNO_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 fixedlength buffers in the
mpPrint*()
functions to avoid buffer overflows  Fixed an issue with the
mpPrint*()
functions when displaying zero values.  Added the
bdNewVars
andbdFreeVars
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) Changed license conditions to Mozilla Public License, v. 2.0.
 Fixed
bdShortMult
to catch empty u parameter.  Fixed
bdPrintBits
to print "0" if bit length is zero.  Changed comparison functions to be constanttime by default. Retained the old, quicker functions as
_q
variants. (Changed back in version 2.6).  Added constanttime modular exponentiation functions
bdModExp_ct
andmpModExp_ct
.  Added
mpPrintDecimalSigned
function to display negative numbers ("mp" only).
 Added markup for Doxygen documentation to produce a manual in CHM and HTML formats.
 Updated the "prettygood" 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 thempConv
functions (thanks to Kevin Kramb for pointing this out).  Added the new functions
bdRandomOctets
,bdRandomNumber
andbdQuickRandBits
, which do various contortions you might want to do with random numbers.  Updated the test source code modules.
 Added new functions to compute the integer cube root:
bdCubeRoot
andmpCubeRoot
.  Updated and improved the integer square root functions:
bdSqrt
andmpSqrt
.  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.
mpPrintHex
andmpPrintDecimal
.  Added the new function
bdPower
to compute the nth power of a numbery = g^{n}
. 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
andmpGcd
.  Fixed the bug in
mpChs
(thanks to Valery Blazhnov).  Fixed the bug in
mpAlloc
which would fail if you calledbdSetDigit(b, 0)
(thanks to "Radistao").
 Added slidingwindow exponentiation to speed up the ModExp functions.
 Added new bd 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
andmpAbs
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 64bit integer types (long long
) if available on a 32bit machine.  Improved the zeroisation and destroy macros.
 Moved the type declarations for the exactwidth types to a separate include file.
 Reorganised and renamed the ancillary source and include files (yet again!). See How to Compile.
 Deprecated the
bd*Ex
functions in favour of renamed "safe"bd*_s
functions.
 Added support for 64bit 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(01).bdRandomBits
generates a random BIGD number up to n bits long (i.e. r = [0, 2^{n}]) using the internal RNG.
 Removed the shift limit on
bdShiftLeft
andbdShiftRight
. This was previously limited to a maximum of 32 bits. It will now shift by any value  well, any value up toSIZE_MAX
.  Reorganised 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.422003 Annex F.1.1 Range of bases in MillerRabin test.  Made minor improvements to the
spBetterRand
function.  Fixed bugs in
bdConvertHex
,bdGetBit
andbdSetBit
.  Made various minor changes in the code to avoid warning messages when compiling with the
Wall pedantic
options.
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
As of version 2.5 (October 2015) we have reissued 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.]
If you use it, we'd really appreciate a link back to this page:
<a href="http://www.dimgt.com.au/bigdigits.html">BigDigits multipleprecision arithmetic source code</a>
Please forward any comments or bug reports to <contact>. The latest version of the source code can be downloaded from www.dimgt.com.au/bigdigits.html or here.
Using the BigDigits Software
A simple example that multiplies the 48bit 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); bdConvFromHex(v, "deadbeefface"); 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 489675628213660which you can verify using the Hex mode of the Windows Calculator in Scientific mode.
All multipleprecision 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 multipleprecision 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); bdConvFromHex(v, "deadbeefface"); 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 fixedlength arrays of sufficient length: this is faster but less safe.
This code is all written in ANSI C (except the optional ASM and Win32specific 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".
Download
 Latest version 2.6: BigDigits2.6.1.zip (117 kB), released under MPLv2 on 31 March 2016.
 Manual in HTML Help (CHM) format (zipped, 135 kB).
The download contains the full source code files for the multipleprecision 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.
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 thebdigit_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 exactwidth types required by both
bigdigits.h
andbigd.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
orbdRandDigit
functions.  bigdRand.c
 Separate source code for the internal RNG functions
bdRandomBits
orbdRandDigit
.  bigdigitsRand.c
 Source code for the internal RNG
spBetterRand
andmpRandomBits
functions. Called by the functions inbigdRand.c()
. We keep this separate so you can substitute your own "best" random function to taste, should you wish.  spExtra.c/.h
 Singledigit versions of some of the multipleprecision 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
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 preprocessor constant USE_SPASM or, alternatively, define USE_64WITH32 to use the native 64bit 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 filempcode.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 64bit 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.
Library  No internal RNG functions  Uses 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 BIGDcompatible 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 subset 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
20120123: 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 64bit 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 64bit 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 8bit microcontroller
Alfredo Ortega has modified the BigDigits code to run on an 8bit microprocessor. See Implementation of the RSA algorithm in a 8bit 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
20091210: 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 easytofollow 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 easytofollow 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.21).Allen Zhao reports success compiling version 2.1 in both 32 and 64bit modes on Solaris 8, HPUX 11i/PARISC; 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 32bit/64bit 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 1025bit 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 PKCS1v1.5 encoding and a freshlygenerated 1024bit RSA key. Again, this just tests the functions. You could rewrite 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 constanttime 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 1862 "Digital Signature Standard (DSS)", 27 January 2000.
 t_bdRDSA.c
 Reproduces some test vectors of the rDSA algorithm from ANSI X9.311998 "Digital Signatures Using Reversible Public Key Cryptography for the Financial Services Industry (rDSA)", September 9, 1998. Appendix D.1 Odd e = 3 with 1024bit n and D.3 Odd e =3 with 2048bit n.
 t_bdRandomOctets.c

Examples using the
bdRandomOctets
andbdRandomNumber
functions.  t_bdQuickRandBits.c

Examples comparing the
bdQuickRandBits
andbdRandomBits
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 508bit 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
 A simple example in C++. See also Compiling with C++.
Overlapping Variables and other "Gotchas"
The mp functions are rather raw and userunfriendly, 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 nonnegative integers 0,1,2,...
.
Negative numbers will "work" with the fixedarray "mp" functions if you are happy using twoscomplement 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 fixedlength "mp" functions.
But be careful with the memorymanaged "bd" functions. If you go below zero with the "bd" functions the results are undefined.
You should
either make sure you always do cutoff subtraction that yields a nonnegative result
(i.e. to compute a  b, first check that b ≤ a; 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.
t_bdTests.c
/* Cope, sort of, with `negative' numbers */
printf("Go ``negative''...\n");
bdSetZero(u);
overflow = bdShortSub(w, u, 2);
pr_msg("w=02=", 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=02=fffffffe overflow=1 w+1=ffffffff overflow=0 (Note error!) w+1=00000001 00000000 overflow=1Thanks to Taylor Francis for reminding us about this issue.
How you could do cutoff 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.
 To generate a few pseudorandom numbers in a limited range to carry out a few tests.
 To generate a large set of random digits to use in a large number of tests.
 To generate a cryptographicallysecure random number in a security application.
 A few pseudorandom numbers for simple tests

In the first case, you can use the standard ANSI
rand()
function from thestdlib
library and seed it withsrand(time(NULL))
or better,srand(time(NULL) ^ clock())
You need to be careful if you want a number greater thanRAND_MAX
, which can be as small as 0x7fff on some systems. We have used this method in thebdSetRandTest
function and thespSimpleRand(DIGIT_T lower, DIGIT_T upper)
function which usesrand()
to generate a random number bytebybyte 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 thespBetterRand()
function, which is based on notquitecryptographicallysecure techniques. This is called from the BIGD library via the synonymbdRandDigit()
and is used bybdRandomBits
, both inbigdRand.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 nonrandom behaviour. The output should pass all FIPS1402 tests for randomness. The function does, however, rely on static variables and so should be used with care in a multithreaded environment. We have kept the code for this function in a separate module so it can easily be replaced if you wish.  Cryptographicallysecure RNGs

In the third case, you would need to call a cryptographicallysecure function. In the
bdRandomSeeded()
andbdGeneratePrime()
functions  which are intended to be used in secure cryptography packages  you pass a reference to a separate RNG callback function which can be as cryptographicallysecure as you wish. Using a separate function enables you to handle multithreaded 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 multipleprecision package, not a cryptographic module). The function has a required formatBD_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 IrelandI 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:
 Use only "pure" ANSI C
 Keep it independent of machine endianness
 Make sure it's reasonably understandable as standalone code
 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.
However, adding slidingwindow 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 frequentlycalled singledigit 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
64bit types like long long
if available.
Programmer's Comments
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 likempModuloTemp
to reuse the same memory for temporary variables for the memoryintensivempModExp
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=2^{32} 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=2^{16} readily convert upwards for 32bits.  spDivide
 We are pretty certain that the "Add Back" step will never happen in our singleprecision
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 32bit 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 32bits). 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
 "Radistao"
 Kevin Kramb
References
Sources of the algorithms used in BigDigits. [KNUTH] Donald E. Knuth, Art of Computer Programming, Volume 2: Seminumerical Algorithms, 3rd ed, AddisonWesley, 1998.
 [MENEZES] Alfred J. Menezes, Paul C. van Oorschot, Scott A Vanstone, Handbook of Applied Cryptography (Discrete Mathematics and Its Applications), CRC Press, 1997, <http://www.cacr.math.uwaterloo.ca/hac/>.
 [COHEN] Henri Cohen, A Course in Computational Algebraic Number Theory, Springer, 1996.
 [SCHNEIER] Bruce Schneier, Applied Cryptography: Protocols, Algorithms, and Source Code in C, 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).
 [FIPS1862] Digital Signature Standard (DSS), FIPS PUB 1862, U.S. Department of Commerce/National Institute of Standards and Technology, 2000 (latest PDF).
 [ANSIX942] ANSI X9.422003 Public Key Cryptography for the Financial Services Industry: Agreement of Symmetric Keys Using Discrete Logarithm Cryptography, American National Standards Institute, 2003.
 [CORON] JeanSebastian Coron, Resistance Against Differential Power Analysis for Elliptic Curve Cryptosystems, Cryptographic Hardware and Embedded Systems, Lecture Notes in Computer Science Volume 1717, 1999, pp 292302. (PDF)
 [BLAKE] Blake, Ian F., Seroussi, G., N.P. Smart, Elliptic Curves in Cryptography, Cambridge University Press, 1999.
Bibliography
Other sources we used in devising tests and other techniques used in the code. [ANSIX931] ANSI X9.311998 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. This is virtually identical to the later edition of Cryptography Engineering
 [HANSON97] David R. Hanson, C Interfaces and Implementations: Techniques for Creating Reusable Software, AddisonWesley, 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 languages  C, 2nd ed., ISO/IEC, 1999.
 [PKCS1] PKCS#1, RSA Cryptography Standard, RSA Laboratories, Version 2.1, June 2002 (republished as RFC3447).
 [PKCSEX] Burton S Kaliski, Jr, Some Examples of the PKCS Standards, RSA Laboratories Technical Note, Nov 1993.
Errata
Last updated: 11 September 20161. Error in t_bdTest.c
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);
Contact
Any comments, feedback, questions: please send us a message.
This page last updated 11 September 2016