Andrew Que Sites list Photos
Projects Contact
Main

December 26, 2013

Multiplication of Huge Numbers

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.
//-----------------------------------------------------------------------------
templateclass 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'.
//-----------------------------------------------------------------------------
templateclass 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.

In C multiplication overflow is ignored. If one multiplies two 32-bit numbers, and the result is larger than 32-bits the upper bits go into the bit bucket. As long as the result will fit into 32-bits, nothing is lost. If the result will be larger than 32-bits, one must use 64-bit data types. However, as of 2011 the C standard has no integer data types larger than 64-bit. So if larger numbers are desired one either has to use floating-point, or make their own multiply. At work this month I ran into a problem where I needed to multiply two 64-bit numbers that would overflow a 64-bit result. After the work day was complete, I wrote a simple function to do this, and we will examine how it works in this article.

When doing multiplication of binary numbers, the result will never be larger than twice the number of bits of the data types used. That is, two 32-bit numbers multiplied together will never have a result larger than 64-bits. The reason can be expressed mathematically.

Where n is the number of bits.

So when multiplying two 64-bit values it is known the result will always fit into a 128-bit. We can construct this result using two 64-bit values—an upper word and a lower word.

Doing the multiplication one has to think back to their grade school days and long multiplication. Let us first solve a simple problem doing long multiplication to solve 12x34.

 

1

2

 

3

4

 

4

8

3

6

0

4

0

8

The steps are to first multiply the least-significant digits of the 34 by the least-significant digit of 12 (in this case 4 and 2) and place the result in the left most column. Then the most-significant digits of 34 by the least-significant digit of 12 (3 and 2) and place that in the second column from the left. Next is the most-significant digit of 34 by the least-significant of 12 (3 and 2) which is stored in the second column from the left, and lastly the most-significant digit of 34 by the most-significant digit of 12 (3 and 1) which is stored in the third column from the left. Sum the two results with carries and we have our end result.

These steps are precisely what is needed in order to multiply two 64-bit numbers. But rather than digits, each column will represent 32-bits. If you still think of each column as a digit, the following chart is an algebraic representation of the above. To change to 32-bits just consider that rather than each column counting from 0 to 9 before carrying the digit to the next column, it counts from 0 to 232 – 1 before carrying to the next column.

 

au

al

 

bu

bl

 

bl*au

bl*al

bu*au

bu*al

0

Here au is the upper 32-bits of the first word (a), and al is the lower 32-bits. The only problem with looking at the setup this way is that it does not show the word sizes. So let us arrange things like this:

   

au

al

   

bu

bl

   

bl*al

 

bl*au

 
 

bu*al

 

bu*au

   

Now we can see that the results of the multiplications are twice the width of the two values being multiplied, and how that lines up with the other values in the finial summation. Again, each column is 32-bits and the four total columns means we have a 128-bit result. We split the input into two 32-bit numbers so that each individual multiplication never has a result larger than 64-bits. This way there is never an overflow. The trade off is that we now have to do 4 multiplications and several additions.

So here is the source code for doing the multiplication.

//******************************************************************************
// Name: mult128.h
// Uses: 128-bit multiplication using 64-bit input.
// Date: 2013/12/08
// Author: Andrew Que (http://www.DrQue.net/)
// Compiler: C99 compliant.
// Revision:
//   1.0 - 2013/12/08 - Creation.
//
//                                Public domain
//******************************************************************************
#ifndef MULT128_H
#define MULT128_H

#include <stdint.h>

// 128-bit unsigned integer.
typedef struct
{
  uint64_t upper;
  uint64_t lower;
} UINT128;

//------------------------------------------------------------------------------
// USES:
//   Multiply two 64-bit unsigned integers to form valueA 128-bit result.
// INPUT:
//   valueA - First multiplicand.
//   valueB - Second multiplicand.
// OUTPUT:
//   128-bit unsigned integer of multiply.
//------------------------------------------------------------------------------
static inline UINT128 multiply_uint128( uint64_t valueA, uint64_t valueB )
{
  int const SHIFT = 32;
  uint64_t const MASK = 0xFFFFFFFF;
  UINT128 result;
  uint64_t product[ 4 ];
  uint64_t multiply;

  //
  // Product is taken by dividing each 64-bit word into two 32-bit, and doing
  // multiplication on the parts.
  //    a b   <- Upper/lower word of 'ValueA'
  //  * c d   <- Upper/lower word of 'ValueB'
  //
  // The complete multiply is then the sum of each of the parts.  Each
  // multiplication spans two words of the result as such:
  //   [3] [2] [1] [0]  <- Resulting full product, held in 'product'.
  //           [b * d]
  //       [a * d]
  //       [b * c]
  //   [a * c]
  //

  // 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 );

  // Store result.
  result.upper = ( product[ 3 ] << SHIFT ) | product[ 2 ];
  result.lower = ( product[ 1 ] << SHIFT ) | product[ 0 ];

  return result;
}

#endif // MULT128_H

There you have it. The code is easily adapted to other integer sizes. By changing the types to uint32_t, and the mask and shift values this function could produce a 64-bit result. C already has built in support for 64-bit types, but it could be done nonetheless.

The source code is here.

1 comment has been made.

From Erica

December 26, 2013 at 4:36 AM

Short and sweet, cute and neat.

December 24, 2013

Digits in a number of n-bits

Here's a little math coding snip. Say you need to create a string of an integer number. How many digits will it have?

The primitive C integer types are 8, 16, 32, and 64 bits. The largest number that can be represented by this number of bits is 2n – 1 where n is the number of bits. The subtract one is due to the fact the count starts at zero.

How does this help determine the number of digits the string will contain? To understand that, let us first examine how many digits various numbers have.

We use Hindu-Arabic numerals to represent numbers. Most people say Arabic numerals, but it was actually ancient Indian mathematicians who developed the system. Arab traders introduced the numbering system to Europe, and hence the name credited the number system to them. I think they both deserve credit hence I refer to them as Hindu-Arabic numerals. These numerals are in base 10 as we have 10 digits, 0 to 9.

Logarithms will tell us how many characters are needed to represent any number in any number base. For example, 1000 has 4 digits. We can also use the log base 10 to calculate the same thing.

Where d is the number of digits. The log base 10 of 1000 has an even output. To get the number of digits for any number will require the floor function to be added:

For example, 12345 has 5 digits. The log of 12345 is ≈4.091. The floor function is just the integer part of the number, so 4. Add one to that and we have our 5 digits.

How does this help with our original problem? Well computers store numbers in powers of two. Most understand this as binary, but it is really powers of two. For example the binary value 1101 = 23 + 22 + 20 = 13. This binary value is 4 bits, and the largest value that can be represented using 4 bits is 24 – 1 = 15. If we know how many digits are in 24 then we know any the maximum number of digits any 4 bits value can have.

This can be used for any number of bits by substituting the exponent of 4 for n:

This can be simplified using logarithm identities:

And the log base 10 of 2 is just a constant, and happens to be ≈0.30103.

That is a fairly simple equation to calculate. Using 0.30103 as the constant produces accurate values up to around 300,000 bits—far more than is needed for any C primitive. Using 0.3 is good enough for all C primitives and works up to 255 bits.

To stick with integer types to avoid floating point values we can add in a multiply by 10.

That's great for C primitive but not too useful—there are plenty of string functions and so few types to worry about. What about very large numbers with an arbitrary number of bits? For this we can use fixed-point math.

Here c is the constant representing log10( 2 ), and f is the fixed-point shift. Any algebra student will tell you that the 2f will cancel, which is true. The reason we use this is because the fixed shift and the constant will be combine to form a new constant: C = 2fc = 2f * log10( 2 ). So the equation becomes:

Remember that a divide by any power of two is actually just a bit shift to the right, and that is very cheap to compute. Now we just need a value for C for as many bits as we would like to have accuracy.

Shift (f) Value (C) Accuracy (bits) Accuracy (bytes)
32 10,343,311,892 274,877,906,936 34,359,738,367
30 2,585,827,973 274,877,906,936 34,359,738,367
24 40,403,562 536,870,904 67,108,863
16 157,826 262,136 32,767
8 617 65,528 8,191

So there it is. Even using an 8-bit shift will allow the calculation to be accurate to number over 8,000 bytes in length—more than 19,000 digits long. Using a 16-bit shift is good for 32,000 byte number—79,000 digits in length. And 30 bits enough for 34 billion bytes—almost 662 billion digits if you ever had the desire to print that many.

   Xiphos has impressive mechanical.  A couple of weeks ago the breaks went out on his truck, and after a lot of effort he traced out some broken break lines the need to be replaced and fighting the sea of rust that is his truck body he managed to change them out.  Now that his truck works, he is working on a car that needs new gaskets.  Here he is cleaning part of the engine block on the sink, removing the old destroyed gasket.