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 New version 2.4 released 27 April 2013: see New in Latest Version 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.

Just released! bdcalc
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.
UpdatedUpdated version 2.0 of bdcalc released 19 October 2014.

Contents

Recommended reading

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.

New in Latest Version

Version 2.4 (Released 27 April 2013)

Changes in Version 2.3

Version 2.3 (Released 11 November 2011)

Changes in Version 2.2

Version 2.2 (Released 31 July 2008)

Changes in Version 2.1

Version 2.1 (Released 23 August 2006)

Changes in Version 2.0

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.

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-13 by D.I. Management Services Pty Limited <www.di-mgt.com.au>, and is used with permission.

If you use it in connection with a web page, please include a link to our web site:

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

There is no need to write and ask permission, but we'd be delighted to hear from developers who use the code in their applications. For the full wording of the copyright notice and licence conditions see bigdigitsCopyright.txt.

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.

Using the BigDigits Software

A simple example that multiplies the number 0xdeadbeefface by 2

#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);
	bdPrintHex("0x", w, "\n");
	/* Free all objects we made */
	bdFree(&u);
	bdFree(&v);
	bdFree(&w);
	return 0;
}

This should display

0x1bd5b7ddff59c
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. Internal memory allocation is handled automatically. 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 below. 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

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.

Please read this small print before using the software

David Ireland and DI Management Services Pty Limited make no representations concerning either the merchantability of this software or the suitability of this software for any particular purpose. It is provided "as is" without express or implied warranty of any kind. Our liability will be limited exclusively to the refund of the money you paid us for the software, namely nothing. By using the software you expressly agree to such a waiver. If you do not agree to the terms, do not use the software. For the full details see bigdigitsCopyright.txt.

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.
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.

If you define the NO_ALLOCS preprocessor macro then the "mp" functions will 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 digits 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

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 digit size you want to handle (the 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_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
A simple example in C++. See also Compiling with C++.

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,.... There are a few "experimental" functions included (mpIsNegative, mpChs and mpAbs) intended for signed integers, but so far we have not developed that any further.

In particular, 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
Thanks to Taylor Francis for reminding us about this issue.

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.

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 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:

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: 23 October 2013

1. bdShortMult function in bigd.c

Before line 609 insert the following

602 int bdShortMult(T w, T u, bdigit_t d)
603     /* Compute w = u * d */
604 {
605    DIGIT_T overflow;
606    size_t dig_size = u->ndigits;
607
608    assert(w && u);
    /***************************************************/
    // [2013-06-05]: catch empty u; make sure w is empty
    if (dig_size == 0 || d == 0)
    {   /* u == 0 */
        bd_resize(w, 0);
        return 0;
    }
    /***************************************************/
609    bd_resize(w, dig_size+1);
610
611    overflow = mpShortMult(w->digits, u->digits, d, dig_size);

If u is zero and has never been used but w has been used before and is not empty, the result in w may not be zero. Changes in detail. To reproduce:

BIGD u, w;
u = bdNew();
w = bdNew();
bdSetShort(w, 666);	   // w = 666
bdSetZero(u);              // u = 0
bdShortMult(w, u, 10);     // w = 0 * 10
bdPrintDecimal("w=", w, " (expecting 0)\n");
bdFree(&u);
bdFree(&w);
w=666 (expecting 0)

Contact

Any comments, feedback, questions: please send us a message.

This page last updated 19 October 2014