HugeInt division using base 2^32 implementation

Question

I previously submitted my code for a HugeInt class here. That version lacked a long division implementation.

I have now implemented (long) division using Donald Knuth's algorithm D, following the code that appeared in Hacker's Delight fairly closely. Since the algorithm relies on some a priori non-obvious facts and its implementation details are sometimes subtle (signed/unsigned implicit conversions, for example), I have annotated the code in unsigned_divide() quite heavily. (I would appreciate opinion on whether I have done too much, perhaps interrupting "algorithmic flow".) I have checked the new function unsigned_divide() quite extensively, and it seems to work correctly, but there are a lot of integers!

I would appreciate a review of the code (given below) with a view to improving its performance and to see if there is perhaps a special/edge case I have overlooked (apart from division by 0, which is purposefully not checked). I include below simple driver code for testing division specifically.

The full code can be found at https://github.com/Richard-Mace/huge-integer-class and contains a few minor improvements since the last iteration (based on many good suggestions I received here), and some pruning of redundant functions.

Thank you in advance for your time and assistance.

Code from HugeInt.cpp: (header file below)

/**
 * friend binary operator /
 * 
 * Return the quotient of two HugeInt numbers. Uses utility function 
 * unsigned_divide, which employs Donald Knuth's long division algorithm. See 
 * comments on implicit conversion before 
 * HugeInt operator+(const HugeInt&, const HugeInt&) above, which are 
 * applicable here also.
 * 
 * @param a
 * @param b
 * @return 
 */
HugeInt operator/(const HugeInt& a, const HugeInt& b) {    
    if (a < 0) {
        if (b < 0) {
            return unsigned_divide(-a, -b, nullptr);
        }
        else {
            return -unsigned_divide(-a, b, nullptr);
        }
    }
    else {
        if (b < 0) {
            return -unsigned_divide(a, -b, nullptr);
        }
        else {
            return unsigned_divide(a, b, nullptr);
        }
    }
}

/**
 * friend binary operator %
 * 
 * Return the remainder from the division of two HugeInt numbers. Uses utility 
 * function unsigned_divide. Adheres to the C/C++ convention that the sign of 
 * the remainder is the same as the sign of the dividend. See comments on 
 * implicit conversion before HugeInt operator+(const HugeInt&, const HugeInt&) 
 * above, which are applicable here also.
 * 
 * @param a
 * @param b
 * @return 
 */

HugeInt operator%(const HugeInt& a, const HugeInt& b) {
    HugeInt remainder;

    if (a < 0) {
        if (b < 0) {
            unsigned_divide(-a, -b, &remainder);
            return -remainder;
        }
        else {
            unsigned_divide(-a, b, &remainder);
            return -remainder;
        }
    }
    else {
        if (b < 0) {
            unsigned_divide(a, -b, &remainder);
            return remainder;
        }
        else {
            unsigned_divide(a, b, &remainder);
            return remainder;
        }
    }
}

/**
 * unsigned_divide: (private utility function)
 * 
 * Unsigned division of a by b giving quotient q = [a/b] and remainder r, such 
 * that 
 *                   a = q * b + r,    where 0 <= r < b.
 * 
 * Dividend a is assumed non-negative (a >= 0) and divisor b is positive 
 * definite (b > 0). If the number of base-2^32 digits in b is 1, then short 
 * division is used. Otherwise Donald Knuth's Algorithm D is used. 
 * 
 * The implementation of Knuth's algorithm here is very similar to that which 
 * appears in the book Hacker's Delight by Henry S. Warren, but borrows some 
 * ideas from janmr's blog entry (to which credit is duly given):
 *  
 *   https://janmr.com/blog/2014/04/basic-multiple-precision-long-division/
 * 
 * If remainder is not a nullptr, then the remainder r is returned in space 
 * allocated by the caller.
 * 
 * WARNING: no checks on the validity of a and b are made for performance 
 * reasons.
 * 
 * @param a
 * @param b
 * @return 
 */

HugeInt unsigned_divide(const HugeInt& a, const HugeInt& b, 
                        HugeInt* const remainder)
{   
    HugeInt dividend{a};
    HugeInt divisor{b};

    // Determine the number of base-2^32 digits in dividend and divisor.
    int n{HugeInt::numDigits_};
    for ( ; n > 0 && divisor.digits_[n - 1] == 0; --n);

    int m{HugeInt::numDigits_};
    for ( ; m > 0 && dividend.digits_[m - 1] == 0; --m);

    // Technically, m can equal 0 here, if 'a' (the dividend) = 0. This is no 
    // problem as it will be caught and handled by CASE 1 below.

    // CASE 1: m < n => quotient = 0; remainder = dividend.
    HugeInt quotient;

    if (m < n) {
        if (remainder != nullptr) {
            if (*remainder != 0) {
                *remainder == 0LL;
            }

            for (int i = 0; i < m; ++i) {
                remainder->digits_[i] = dividend.digits_[i];
            }
        }

        return quotient;
    }

    // CASE 2: Divisor has only one base-2^32 digit (n = 1). Do a short 
    //         division and return.
    if (n < 2) {
        std::uint64_t partial{0};

        for (int i = m - 1 ; i >= 0; --i) {
            partial = HugeInt::base_ * partial 
                    + static_cast<std::uint64_t>(dividend.digits_[i]);
            quotient.digits_[i] = 
                    static_cast<std::uint32_t>(partial / divisor.digits_[0]);
            partial %= divisor.digits_[0];
        } 

        if (remainder != nullptr) {
            if (*remainder != 0) {
                *remainder == 0LL;
            }

            remainder->digits_[0] = partial;
        }

        return quotient;
    }

    // CASE 3: m >= n and the number of digits, n, in the divisor is >= 2. 
    // Proceed with long division using Donald Knuth's Algorithm D.
    //
    // Determine power-of-two normalisation factor, d = 2^shifts, necessary for
    // d * divisor.digits[n-1] >= base_ / 2.
    int shifts{0};
    std::uint32_t vn{divisor.digits_[n - 1]};

    while (vn < (HugeInt::base_ >> 1)) {
        vn <<= 1;
        ++shifts;
    }

    // Scale the divisor and dividend by factor d, using shifts for efficiency. 
    // This scaling does not affect the quotient, but it ensures that
    // q_k <= qhat <= q_k + 2 (see later).
    for (int i = n - 1; i > 0; --i) {
        divisor.digits_[i] = (divisor.digits_[i] << shifts) | 
          (static_cast<std::uint64_t>(divisor.digits_[i - 1]) >> (32 - shifts));
    }
    divisor.digits_[0] = divisor.digits_[0] << shifts;

    // Prepend a (m+1)'th zero-value digit to the dividend, then shift.
    dividend.digits_[m] = 
        static_cast<std::uint64_t>(dividend.digits_[m - 1]) >> (32 - shifts);
    for (int i = m - 1; i > 0; --i) {
        dividend.digits_[i] = (dividend.digits_[i] << shifts) |
         (static_cast<std::uint64_t>(dividend.digits_[i - 1]) >> (32 - shifts));
    }
    dividend.digits_[0] = dividend.digits_[0] << shifts;

    // Do the long division using the primary school algorithm, estimating
    // partial quotients with a two most significant digit approximation for 
    // the dividend and a single most significant digit approximation for the 
    // divisor.
    for (int k = m - n; k >= 0; --k) {
        std::uint64_t rhat = dividend.digits_[k + n] * HugeInt::base_ 
            + static_cast<std::uint64_t>(dividend.digits_[k + n - 1]);

        std::uint64_t qhat = rhat / divisor.digits_[n - 1];

        rhat %= divisor.digits_[n - 1];

        // Digit q_k estimated by qhat must satisfy 0 <= q_k <= base_ - 1. 
        // If too large, decrement and adjust remainder rhat accordingly.
        if (qhat == HugeInt::base_) {
            qhat -= 1;
            rhat += divisor.digits_[n - 1];
        }

        // Compare with a "second order" approximation to the partial quotient. 
        // If this comparison indicates that qhat overestimates, decrement,
        // adjust remainder rhat and repeat.
        while (rhat < HugeInt::base_ && (qhat * divisor.digits_[n - 2]
                > HugeInt::base_ * rhat + dividend.digits_[k + n - 2])) {
            qhat -= 1;
            rhat += divisor.digits_[n - 1];
        }

        // We have an estimate qhat for the true digit q_k that satisfies
        // q_k <= qhat <= q_k + 1. Calculate the corresponding remainder 
        // (a_{k+n} ... a_{k}) - qhat * (b_{n-1}...b_{0}) for this partial
        // quotient, storing the result in digits a_{k+n}... a_{k} of the 
        // dividend. Care is taken with the carries. The overwritten digits 
        // accrue, and eventually become, the complete remainder.

        std::int64_t carry{0};     // signed; carry > 0, borrow < 0
        std::int64_t widedigit;    // signed
        for (int i = 0; i < n; ++i) {
            std::uint64_t product = static_cast<std::uint32_t>(qhat) 
                            * static_cast<std::uint64_t>(divisor.digits_[i]);

            widedigit = (dividend.digits_[k + i] + carry) 
                        - (product & 0xffffffffLL);

            dividend.digits_[k + i] = widedigit; // assigns 2^32-complement
                                                 // if widedigit < 0

            carry = (widedigit >> 32) - (product >> 32);
        }

        widedigit = dividend.digits_[k + n] + carry;
        dividend.digits_[k + n] = widedigit;           // 2^32-complement if
                                                       // widedigit < 0

        // Accept and store the tentative quotient digit.
        quotient.digits_[k] = qhat;

        // However, since q_k <= qhat <= q_k + 1, either we have the correct
        // digit, or we need to decrement. To resolve this, check if there was 
        // a borrow on determining the final k + n digit of the remainder. If 
        // no, we have q_k = qhat and we are done. Otherwise, qhat = q_k + 1, 
        // and we need to decrement and add the divisor to digits k + n ... k 
        // of the dividend (now the remainder).
        if (widedigit < 0) {
            quotient.digits_[k] -= 1;
            widedigit = 0;
            for (int i = 0; i < n; ++i) {
                widedigit += static_cast<std::uint64_t>(dividend.digits_[k + i])
                           + divisor.digits_[i];
                dividend.digits_[k + i] = widedigit;
                widedigit >>= 32;
            }

            dividend.digits_[k + n] += carry;
        }
    } /* end main loop over k */

    // We are done. Return the remainder?
    if (remainder != nullptr) {
        if (*remainder != 0) {
            *remainder == 0LL;
        }

        // Denormalise dividend, which now contains the full remainder 
        // (stored in n - 1 digits). 
        for (int i = 0; i < n - 1; ++i) {
            remainder->digits_[i] = (dividend.digits_[i] >> shifts) |
                    (static_cast<std::uint64_t>(dividend.digits_[i + 1]) 
                        << (32 - shifts));
        }

        remainder->digits_[n - 1] = dividend.digits_[n - 1] >> shifts;
    }

    return quotient; 
}

HugeInt.h:

/*
 * HugeInt.h
 * 
 * Definition of the huge integer class
 * Richard Mace, February, 2020
 * 
 * RADIX 2^32 VERSION
 *
 * Huge integers are represented as N-digit arrays of uint32_t types, where
 * each uint32_t value represents a base-2^32 digit. By default N = 300, which 
 * corresponds to a maximum of 2890 decimal digits. Each uint32_t contains 
 * a single base-2^32 digit in the range 0 <= digit <= 2^32 - 1. If `index' 
 * represents the index of the array of uint32_t digits[N], 
 * i.e., 0 <= index <= N - 1, and 'value' represents the power of 2^32 
 * corresponding to the radix 2^32 digit at 'index', then we have the following 
 * correspondence:
 *
 * index  |...... |     4    |     3    |     2    |     1    |     0    |
 * -----------------------------------------------------------------------
 * value  |...... | (2^32)^4 | (2^32)^3 | (2^32)^2 | (2^32)^1 | (2^32)^0 |
 *
 * The physical layout of the uint32_t array in memory is:
 *
 * uint32_t digits[N] = {digits[0], digits[1], digits[2], digits[3], ... }
 *
 * which means that the units (2^32)^0 appear first in memory, while the power 
 * (2^32)^(N-1) appears last. This LITTLE ENDIAN storage represents the 
 * number in memory in the REVERSE order of the way we write decimal numbers, 
 * but is convenient.
 *
 * Negative integers are represented by their radix complement. With the 
 * base 2^32 implementation here, we represent negative integers by their base 
 * 2^32 complement. With this convention the range of 
 * non-negative integers is:
 *                      0 <= x <= (2^32)^N/2 - 1
 * The range of base 2^32 integers CORRESPONDING to negative values in the
 * base 2^32 complement scheme is:
 *                      (2^32)^N/2 <= x <= (2^32)^N - 1 
 * So -1 corresponds to (2^32)^N - 1, -2 corresponds to (2^32)^N - 2, and so on.
 * 
 * The complete range of integers represented by a HugeInt using radix 
 * complement is:
 * 
 *                     -(2^32)^N/2 <= x <= (2^32)^N/2 - 1
 */


#ifndef HUGEINT_H
#define HUGEINT_H

#include <string>
#include <iosfwd>

namespace iota {

class HugeInt {
public:
    HugeInt() = default;
    HugeInt(long long int);  // conversion constructor from long long int
    explicit HugeInt(const char* const); // conversion constructor from C string
    HugeInt(const HugeInt&);    // copy/conversion constructor

    // assignment operator
    const HugeInt& operator=(const HugeInt&);

    // unary minus operator
    HugeInt operator-() const;

    // conversion to long double
    explicit operator long double() const;

    // basic arithmetic
    friend HugeInt operator+(const HugeInt&, const HugeInt&);
    friend HugeInt operator-(const HugeInt&, const HugeInt&);
    friend HugeInt operator*(const HugeInt&, const HugeInt&);
    friend HugeInt operator/(const HugeInt&, const HugeInt&);
    friend HugeInt operator%(const HugeInt&, const HugeInt&);


    // increment and decrement operators
    HugeInt& operator+=(const HugeInt&);
    HugeInt& operator-=(const HugeInt&);
    HugeInt& operator*=(const HugeInt&);
    HugeInt& operator/=(const HugeInt&);
    HugeInt& operator%=(const HugeInt&);
    HugeInt& operator++();     // prefix
    HugeInt  operator++(int);  // postfix
    HugeInt& operator--();     // prefix
    HugeInt  operator--(int);  // postfix

    // relational operators
    friend bool operator==(const HugeInt&, const HugeInt&);
    friend bool operator!=(const HugeInt&, const HugeInt&);
    friend bool operator<(const HugeInt&, const HugeInt&);
    friend bool operator>(const HugeInt&, const HugeInt&);
    friend bool operator<=(const HugeInt&, const HugeInt&);
    friend bool operator>=(const HugeInt&, const HugeInt&);

    // input/output 
    std::string toRawString() const;
    std::string toDecimalString() const;
    friend std::ostream& operator<<(std::ostream&, const HugeInt&);
    friend std::istream& operator>>(std::istream&, HugeInt&);

    // informational
    int numDecimalDigits() const;
    static HugeInt getMinimum();
    static HugeInt getMaximum();

private:
    static const std::size_t   numDigits_{300};   // max. no. base 2^32 digits
    static const std::uint64_t base_{1ULL << 32}; // 2^32, for convenience
    std::uint32_t              digits_[numDigits_]{0}; // base 2^32 digits

    // private utility functions
    bool          isZero() const;
    bool          isNegative() const;
    HugeInt&      radixComplement();  
    HugeInt       shortMultiply(std::uint32_t) const;
    HugeInt       shortDivide(std::uint32_t, std::uint32_t* const) const;
    friend HugeInt unsigned_divide(const HugeInt&, const HugeInt&, 
                                   HugeInt* const);
    HugeInt&      shiftLeftDigits(int);
};

} /* namespace iota */

#endif /* HUGEINT_H */

Some rough test code:

/*
 * Some simple tests of the HugeInt division algorithm.
 * 
 * February 2020.
 */

#include "HugeInt.h"
#include <iostream>
#include <string>

// Calculate n!
//
iota::HugeInt factorial_iterative(const iota::HugeInt& n) {
    iota::HugeInt result{1LL};

    if (n == 0LL) {
        return result;
    }

    for (iota::HugeInt i = n; i >= 1; --i) {
        result *= i;
    }

    return result;
}
// Calculate the n'th Fibonacci number
//
iota::HugeInt fibonacci_iterative(const iota::HugeInt& n) {
    const iota::HugeInt zero;
    const iota::HugeInt one{1LL};

    if ((n == zero) || (n == one)) {
        return n;
    }

    iota::HugeInt retval;
    iota::HugeInt fib_nm1 = one;
    iota::HugeInt fib_nm2 = zero;

    for (iota::HugeInt i = 2; i <= n; ++i) {
        retval = fib_nm1 + fib_nm2;
        fib_nm2 = fib_nm1;
        fib_nm1 = retval;
    }

    return retval;
}

//
// Return the nth Fermat number. (n >= 1)
//
// Fermat(n) = 2^(2n) + 1.
//
//

iota::HugeInt fermat_number(long int n) {
    long long int exponent = (1LL << n);

    iota::HugeInt retval{2LL};
    for (long long int i = 1; i < exponent; ++i) {
        retval *= 2LL;
    }

    return retval + 1LL;
}


// Check the results of a division by computing
//
//      check = quotient * divisor + remainder
//
// and verifying that check == dividend.
//

void check_division(std::string case_name, 
                    iota::HugeInt dividend, 
                    iota::HugeInt divisor)
{
    iota::HugeInt quotient = dividend / divisor;
    iota::HugeInt remainder = dividend % divisor;

    iota::HugeInt check = quotient * divisor + remainder;

    std::cout << "TESTING " << case_name << " -----------------------------------------------------------------------\n\n";
    std::cout << '\t' << dividend << " / " << divisor << '\n';
    std::cout << "\tquotient = " << quotient << '\n';
    std::cout << "\tremainder = " << remainder << '\n';

    if (check != dividend) {
        std::cout << "\nERROR: " << case_name << " failure!\n\n";
    }
    else {
        std::cout << "\nPASS: " << case_name << " succeeded!\n\n";
    }
}


int main() {
    // zero divisor //////////////////////////////////////////////////////////// 
    iota::HugeInt dividend;
    iota::HugeInt divisor{"123456789098765432101234567890"};

    check_division("zero dividend (1)", dividend,  divisor);
    check_division("zero dividend (2)", dividend, -divisor);

    ///short division //////////////////////////////////////////////////////////
    dividend = static_cast<iota::HugeInt>("31415926123456789098765432101234567890987654321");
    divisor = (1ULL << 32) - 1LL;

    check_division("short division (1)",  dividend,  divisor);
    check_division("short division (2)",  dividend, -divisor);
    check_division("short division (3)", -dividend,  divisor);
    check_division("short division (4)", -dividend, -divisor);

    // long division ///////////////////////////////////////////////////////////
    dividend = static_cast<iota::HugeInt>("314159261234567890987654321012345678909876543210987657878726353557575751098");
    divisor = static_cast<iota::HugeInt>("9086565656538783989846661928745638291028749009092778710909291");

    check_division("long division (1)", dividend,  divisor);
    check_division("long division (2)", dividend, -divisor);
    check_division("long division (3)", -dividend,  divisor);
    check_division("long division (4)", -dividend, -divisor);

    // some fun (use a wide terminal) //////////////////////////////////////////
    dividend = factorial_iterative(1000LL);
    divisor  = fibonacci_iterative(10000LL);

    check_division("fun 1 division (1)", dividend,  divisor);
    check_division("fun 1 division (2)", dividend, -divisor);
    check_division("fun 1 division (3)", -dividend,  divisor);
    check_division("fun 1 division (4)", -dividend, -divisor);

    // some more fun
    dividend = factorial_iterative(1100LL);
    divisor  = fibonacci_iterative(13000LL);

    check_division("fun 2 division (1)", dividend,  divisor);
    check_division("fun 2 division (2)", dividend, -divisor);
    check_division("fun 2 division (3)", -dividend,  divisor);
    check_division("fun 2 division (4)", -dividend, -divisor);

    // fermat test: This choice of dividend and divisor turns out to be a (very)
    // rare case where we have to correct for an over-borrow in unsigned_divide 
    // by decrementing qhat and giving divisor to the remainder.

    dividend = fermat_number(9LL);
    divisor = static_cast<iota::HugeInt>("7455602825647884208337395736200454000000170665201");
    check_division("special case", dividend, divisor);

    return 0;
}

Answer 1

About the comments

Some of the comments in your code are very nice. For example, when explaining the different cases in unsigned_divide(). However, avoid comments that are just literal repetitions of the code just below. There are some borderline cases of that in your code, such as:

// Accept and store the tentative quotient digit.

And:

// Prepend a (m+1)'th zero-value digit to the dividend, then shift.

Also, if a block of code can be replaced by a call to a function with a descriptive name, that can then avoid needing a comment.

Add a description to @param and @return values

It takes little effort to write an actual description for the parameters and return value in the Doxygen comments, and it will make the resulting documentation better. For example, in the documentation for operator/, write:

 * @param a Left hand side of the division.
 * @param b Right hand side of the division.
 * @return  Returns a / b.

Also remember to document all parameters, for example I se no @param remainder in the documentation of unsigned_divide(). You can turn on warnings in Doxygen, so it will tell you if you forgot to document all parameters of each function for example.

Consider splitting off parts of unsigned_divide() into separate functions

This function is quite long, and could perhaps be split into separate functions. For example, the three cases could be separated, and in the third case, the actual long division algorithm could be put into its own function as well.

Consider adding a user-defined literal operator

You could add a user-defined literal operator that just wraps the conversion constructor, so it becomes easy to write a literal HugeInt in code that uses this class. For example:

HugeInt operator "" _huge(unsigned long long value) {
    return HugeInt(value);
}

HugeInt operator "" _huge(const char *text) {
    return HugeInt(text);
}

And then you could write:

HugeInt the_answer = 6_huge * "7"_huge;

Consider using std::vector<uint32_t> to store the digits

One issue I see with this class, which will impact performance as well, is that you allocate a static array of 1200 bytes to hold the digits. This wastes a lot of memory for smaller numbers, and will obviously fail when the numbers grow large enough to exceed this size.

There are several reasons why this can have negative effect on performance, even though it seems at first that a static array has none of the overhead of a std::vector. However, if you have many HugeInts (maybe you have an array of them), you are using a lot of memory that has to be paged in. Also, stack space could be limited, and when doing calculations it's quite normal to have several variables on the stack, so just having 4 HugeInts would already consume 2 pages. This becomes a real problem when you have some recursive algorithm: the total stack usage is then the deepest recursion level times the memory used by a single iteration.

The large size can also negatively impact memory layout, making it harder for the CPU's caches and prefetcher. Consider a function that starts with:

void somefunc(...) {
    int x;
    HugeInt y;
    int z;
    ...

There is likely a huge gap between z and x in the above code fragment.

I think it's likely in real use cases that a significant fraction of HugeInts will actually only hold small values. It might be worthwhile to investigate if some kind of vector container with small size optimization could be used. This would keep small HugeInts fast, and for the large ones, the time taken by operations on them will dwarf the performance loss of the indirection of a dynamic allocation of the digits.