Multiply and divide of variable bit floating point numbers

Question

I am developing a BigNumber library which allows for creation of large numbers (Integers and Floating). You can find the public repository here

I've implemented add and subtract for both BigFloating and BigInteger, however multiply and divide only for BigInteger

The bits for a floating point number are stored in a std::vector and are in the format:

sign (1 bit)|binary integer(variable bits)|binary fraction(variable bits)

Such that the number 5.5 would have the bits 0 101 1

What are some algorithms for multiplying and dividing numbers with this format?

i.e.

(5.5)   101 1
* (5.5) 101 1
-------------
= (30.25) 11110 01

or

(5.5)   101 1
/ (5.5) 101 1
-------------
= (1)     1

The functions to be implemented are found in BigCommon.cpp and are:

std::tuple<std::vector<bool>, size_t, size_t> BigCommon::multiplyBits(const std::vector<bool> &lhsBits, const std::vector<bool> &rhsBits, const size_t &integerBitsSize, const size_t &mantissaBitsSize)

and

std::tuple<std::vector<bool>, size_t, size_t> BigCommon::divideBits(const std::vector<bool> &lhsBits, const std::vector<bool> &rhsBits, const size_t &integerBitsSize, const size_t &mantissaBitsSize)

Update

I have implemented the multiplyBits algorithm like so:

std::tuple<std::vector<bool>, size_t, size_t> BigCommon::multiplyBits(const std::vector<bool> &lhsBits, const std::vector<bool> &rhsBits, const size_t &integerBitsSize, const size_t &mantissaBitsSize)
{
 std::vector<bool> result;
 result.insert(result.begin(), lhsBits[0] ^ rhsBits[0]);
 size_t newIntegerBitsSize = 0;
 size_t newMantissaBitsSize = mantissaBitsSize + mantissaBitsSize;
 std::vector<bool> lhsBinary(lhsBits.begin() + 1, lhsBits.end());
 std::vector<bool> rhsBinary(rhsBits.begin() + 1, rhsBits.end());
 std::vector<bool> multResult = multiplyBinaryVectors(lhsBinary, rhsBinary);
 newIntegerBitsSize = multResult.size() - newMantissaBitsSize;
 result.insert(result.begin() + 1, multResult.begin(), multResult.end());
 return {result, newIntegerBitsSize, newMantissaBitsSize};
};

Now just for divide!

Update 2

I have successfully implemented division using the following algorithm:

Code redacted in favour of answer.

Update 3

After some testing I found the division algorithm doesn't work with some types of numbers, here are some test cases: 5 / 0.27 10 / 100. Definitely to do with divideBinaryVectors

Answer 1

What you describe is a binary fixed point representation with variable(unbounded) size both above and below the fixed binary point (the split between integer and fractional bits.

For this kind of representation, multiplication is pretty easy -- just concatenate the integer and fractional bits of each number, multiply, and split into integer/fraction again. If your input number have j and k fractional bits respectively, the result will have j+k fractional bits.

Divide is harder as you need to decide how far to continue the divide (how many additional fractional bits) to compute when the result is not exatly representable in binary. Fundamentally, it is just long division -- concatenate the integer/fraction bits, divide (a process that may not terminate) and then split the result. So for example, 10/5.5 gives you

           1.11010 00101 11010 ....
       /----------------------
101.1 / 1010.00000 00000 00000
         101 1
         100 10000 00000 00000
          10 11
           1 11000 00000 00000
           1 011
              1100 00000 00000
              1011
                 1 00000 00000
                   1011
                    1010 00000
                     101 1
                     100 10000
                      10 11
                       1 11000
                       1 011
                          1100
                          1011
                             1

with that repeating 10 bit value in the result. You need to decide when you want to cut it off and round the fraction.

Answer 2

My non-reputable answer to What are some algorithms for multiplying and dividing big numbers using a std::vector<bool> for sign and significand and a size_t for integer bits?

Best advice for evaluation faster than abysmally slow:
Convert to a useful format, calculate, and convert back.

For calculating "in" this format:
Multiplication is relatively immune to differences in number representations.
Algorithms to know are long multiplication and Toom, namely Toom-2/Karatsuba.

Non-performing/pre-checking division (lumped with restoring division more often than not) isn't much more effort than non-restoring due to comparison allowing early out: for every quotient bit an average about three bits need to be inspected.

Beyond missing specification, I find two problems with the implementation attempt in BigCommon.cpp as of 2023/05/04:

• the way to terminate calculation based on parameter mantissaBitsSize is more than dubious
• assuming remainder >= divisor does the right thing, use subtract.
  Just ex-oring remainder and divisor bits is meaningless regardless of direction of iteration.

Answer 3

The cpp_bin_float class from the Boost.Multiprecision library implements arbitrarily large floating point arithmetic with binary representation as you seek. It also supports special edge cases for infinities and NaNs.

The division implementation is here. It's a bit complex, but you can run a simple example code under debugger to understand the logic. Basically, as was already said here, for division you just assume the numbers as integers and then calculate the location of the point separating the fractional part. But first you have to make the decision on the target precision you want to obtain. boost's implementation also uses arbitrary long integer arithmetic behind the scenes.

A note about your implementation:
It's strongly recomended to avoid usage of std::vector<bool>. A better alternative might be boost::dynamic_bitset.

Answer 4

To perform division on fixed-point binary numbers we can use the following algorithm, a modified version of long-division that correctly updates integerBitsSize and fractionBitsSize

std::tuple<std::vector<bool>, size_t, size_t> BigCommon::divideBits(const std::vector<bool> &lhsBits, const std::vector<bool> &rhsBits, const size_t &integerBitsSize, const size_t &fractionBitsSize)
{
 if (isZero(rhsBits))
 {
  throw std::domain_error("Division by zero");
 }

 bool resultSign = lhsBits.front() ^ rhsBits.front();

 std::vector<bool> lhs(lhsBits.begin() + 1, lhsBits.end());
 std::vector<bool> rhs(rhsBits.begin() + 1, rhsBits.end());

 // Normalize dividend and divisor to remove leading zeros
 normalizeBits(lhs);
 normalizeBits(rhs);

 std::vector<bool> dividend = lhs;
 std::vector<bool> divisor = rhs;

 std::vector<bool> result;
 std::vector<bool> tempDividend;

 size_t resultIntegerBitsSize = 0;
 size_t resultFractionBitsSize = 0;

 size_t dividendSize = dividend.size();
 size_t dividendIndex = 0;

 while (true)
 {
  tempDividend.insert(tempDividend.end(), dividendIndex < dividendSize ? dividend[dividendIndex] : false);
  dividendIndex++;
  if (binaryVectorLessThanOrEqualsTo(divisor, tempDividend))
  {
   result.insert(result.end(), true);
   tempDividend = subtractBinaryVectors(tempDividend, divisor);
  }
  else
  {
   result.insert(result.end(), false);
  }
  if (dividendIndex <= dividendSize)
  {
   resultIntegerBitsSize++;
  }
  else
  {
   resultFractionBitsSize++;
   if (resultFractionBitsSize >= 128)
   {
    break;
   }
  }
 }

 while (resultIntegerBitsSize > 0 && !result.front())
 {
  result.erase(result.begin());
  resultIntegerBitsSize--;
 }

 while (resultFractionBitsSize > 0 && !result.back())
 {
  result.erase(result.end() - 1);
  resultFractionBitsSize--;
 }

 result.insert(result.begin(), resultSign);
 return {result, resultIntegerBitsSize, resultFractionBitsSize};
};

It has been tested on the failing test cases 5 / 0.27 and 10 / 100 and is returning correct results.