The other day I wrote about doing multiplication of two 64-bit numbers. What about doing multiplication of arbitrary long numbers? Typically when this is the case, one uses an arbitrary precision math library. There are several available and they are usually the best way to take care of the times when large numbers are needed. However, for academic reasons I thought I'd explore this topic.

Multiplying two large number means one first has to define how the numbers are stored. Naturally arrays will be used. But how are the arrays arranged? They could be made of any bit width. And there is the question of endianness—do we store the most significant word first or last? What if we desire all of this to be variable?

To accomplish this we can use C++ templates. We need templates because we do not know the data types of the arrays, but by using template we can allow this to be figured out at compile time.

// Name: multiplyArray.h

// Uses: Multiply arrays of arbitrary length and type.

// Date: 2013/12/20

// Author: Andrew Que (http://www.DrQue.net/)

// Compiler: C++98 compliant.

// Revision:

// 0.9 - 2013/12/20 - Creation.

// To be done:

// + Add support for signed numbers.

//

//

// The MIT License (MIT)

//

// Copyright (c) 2013 Andrew Que (http://www.DrQue.net/)

//

// Permission is hereby granted, free of charge, to any person obtaining a copy

// of this software and associated documentation files (the "Software"), to deal

// in the Software without restriction, including without limitation the rights

// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell

// copies of the Software, and to permit persons to whom the Software is

// furnished to do so, subject to the following conditions:

//

// The above copyright notice and this permission notice shall be included in

// all copies or substantial portions of the Software.

//

// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR

// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,

// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE

// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER

// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,

// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN

// THE SOFTWARE.

//

//******************************************************************************

#ifndef MULTIPLYARRAY_H

#define MULTIPLYARRAY_H

//-----------------------------------------------------------------------------

// USES:

// Multiply two unsigned words and return the result in an upper and lower

// word of the same width.

// INPUT:

// valueA - First multiplicand.

// valueB - Second multiplicand.

// upperWord - Location to store result.

// lowerWord - Location to store result.

// OUTPUT:

// upperWord - Upper bits of result.

// lowerWord - Lower bits of result.

//-----------------------------------------------------------------------------

template< class TYPE >

void multiplyWords

(

TYPE valueA,

TYPE valueB,

TYPE * upperWord,

TYPE * lowerWord

)

{

int const BITS_PER_BYTE = 8;

int const WORD_SIZE = sizeof( TYPE );

// The shift is half the word size in bits, used to select the upper word.

int const SHIFT = WORD_SIZE * BITS_PER_BYTE / 2;

// The mask is half the bits of the word, used to select the lower word.

TYPE const MASK = ( (TYPE)1 << SHIFT ) - 1;

TYPE multiply;

TYPE product[ 4 ];

// b * d

multiply = ( valueA & MASK ) * ( valueB & MASK );

product[ 0 ] = ( multiply & MASK );

// a * d

multiply =

( valueA >> SHIFT ) * ( valueB & MASK ) + ( multiply >> SHIFT );

product[ 1 ] = ( multiply & MASK );

product[ 2 ] = ( multiply >> SHIFT );

// b * c

multiply = product[ 1 ] + ( valueA & MASK ) * ( valueB >> SHIFT );

product[ 1 ] = ( multiply & MASK );

// a * c

multiply =

product[ 2 ] + ( valueA >> SHIFT ) * ( valueB >> SHIFT ) + ( multiply >> SHIFT );

product[ 2 ] = ( multiply & MASK );

product[ 3 ] = ( multiply >> SHIFT );

*upperWord = ( product[ 3 ] << SHIFT ) | product[ 2 ];

*lowerWord = ( product[ 1 ] << SHIFT ) | product[ 0 ];

}

//-----------------------------------------------------------------------------

// USES:

// Multiply two unsigned arrays of words.

// INPUT:

// valueA - First multiplicand.

// valueB - Second multiplicand.

// product - Storage for result. Must be allocated and large enough for

// result.

// wordsA - Number of words in 'valueA'.

// wordsB - Number of words in 'valueB'.

// isBigEndian - True for big endian (most significant word last), false for

// little endian.

// OUTPUT:

// product - The product of 'valueA' and 'valueB'.

//-----------------------------------------------------------------------------

template< class TYPE >

void multiplyArray

(

TYPE const * valueA,

TYPE const * valueB,

TYPE * product,

int wordsA,

int wordsB,

bool isBigEndian

)

{

int const PRODUCT_WORDS = wordsA + wordsB;

// Start product at zero.

for ( int index = 0; index < PRODUCT_WORDS; ++index )

product[ index ] = 0;

// All words in value A.

for ( int countA = 0; countA < wordsA; ++countA )

{

int indexA = countA;

if ( ! isBigEndian )

indexA = wordsA - countA - 1;

// All words in value B.

for ( int countB = 0; countB < wordsB; ++countB )

{

int indexB = countB;

if ( ! isBigEndian )

indexB = wordsB - countB - 1;

TYPE upperWord;

TYPE lowerWord;

multiplyWords< TYPE >

(

valueA[ indexA ],

valueB[ indexB ],

&upperWord,

&lowerWord

);

// Accumulate this value into product.

for ( int carryCount = ( countA + countB ); carryCount < PRODUCT_WORDS; ++carryCount )

{

int carryIndex = carryCount;

if ( ! isBigEndian )

carryIndex = PRODUCT_WORDS - carryCount - 1;

TYPE accumulator = lowerWord + product[ carryIndex ];

if ( accumulator < lowerWord )

upperWord += 1;

product[ carryIndex ] = accumulator;

lowerWord = upperWord;

upperWord = 0;

}

}

}

}

#endif // MULTIPLYARRAY_H

Let's examine the code. The primary function is *multiplyArray*. This template function can except any type that handles math operation and bit operations. Like all template function, it must be fully define in the header file so the compiler to create the correct code when the function is used. There are three loops nested in one an other. The first two loop through each word of their respective array. Inside the first two loops is where the multiplication of the two selected words take place. Note there is a call to a multiply template function. This template function is almost identical to the one I wrote about yesterday with the difference being it is now a template function that can handle whatever data type is specified. In the future I will write about why this is useful.

Once multiplied, the results must be accumulated into the running sum. This is done in the last loop. One item to pay attention to is the accumulator, and that test to see if the accumulator is less than the lower word. C provides no way to tell if adding two numbers has overflowed—there is no “add plus carry”. To check for overflow, the result will be less than either of the values added together. If that happens the code sets the upper word to account for the carry.

It should noted that no matter what values are being multiplied the upper word will never overflow by adding one. The upper word only starts with any overflow from the multiplication, which is always less than the total the upper word could be. After that, the upper word is either one or zero.

The design of this code is setup to be either byte endian. If the endianness of the words is setup the same as the machine endianness, one can treat the arrays as unions—an array of 32-bit integers could also be 8-bit bytes. The design is also completely portable and machine independent. Any C++98 compliant compiler will work.

For testing this unit I created some functions that turn the array into a string, and came up with some test vectors I compared against values computed with GP/PARI Calculator. Maybe one day I will make a simple library that includes the string functions and test vectors. For now my test code is not available.

I have released this source code under the MIT license. I chose this license because I feel the work I am doing is so basic it doesn't deserve to be maintained as open-source. At the same time it is complete enough to be used as a library. The MIT license lets people use the source by itself, or in a derived work, open or closed source, and keeps my name in the credits. Why is that important? While there is some pride in having my name attached to my work, it's also a matter of promotion. I am a contract software engineer, and if someone sees my work and wants to hire me to do more they will know who to contact if my name is in the license. It is a long-shot, but I would select a programmer for hire based on software I saw they had written.