Huge integer class using base 2^32 (was 256) follow up

Question

Thanks again to everyone for their very useful suggestions concerning a huge integer class I posted here several days ago. I have made many of the suggested changes, perhaps the biggest being in the implementation. In particular, I am now using an array of uint32_t type to store digits with each member representing a base \$2^{32}\$ digit. This replaces the uint8_t used in the previous version. The radix complement is now base - \$2^{32}\$ complement. This change of type has improved speed performance dramatically. For example, in the previous version calculating \$1000!\$ took about 43 seconds on my laptop. It is now done in much less than 1 second! I have not made the class a template, as suggested, but that is quite easy to do. I have not made the thousands separator adhere to local practice. I acknowledge the importance of this and it should be done, but, right now, I just need something that makes sense to a numerate human.

Other changes are:

1. The class is in its own namespace.
2. The non-class inline utility functions are now in an anonymous namespace to restrict access to the translation unit (file).
3. A conversion constructor now creates/converts from long long int rather than from long int, to make the code more friendly to other operating systems.
4. Cleaned up the code for the constructor from a C string, which now validates the supplied string. (I did not change type to std::string because 99% of the time the string supplied will be a literal.)
5. Fixed undefined (okay, plain buggy) behaviour in isZero() and toRawString(). Many thanks to those who pointed this out.
6. Added static member functions getMinimum() and getMaximum() to return the minimum and maximum, respectively, of the integers that are representable in the base \$2^{32}\$ scheme.
7. Added a member function numDecimalDigits() which returns the number of decimal digits in the decimal representation of *this.
8. Added explicit casts in the inner loops of the arithmetic operators +, -, etc., to ensure that the operands are promoted to uint64_t before carrying out addition, multiplication, etc., so that carries propagate correctly.
9. Updated the driver program, which now runs much faster and calculates many, many more digits.
10. Other minor code aesthetic changes.

The default implementation has numDigits hard coded to 300, which gives 2890 decimal digits. You can go considerably higher than this if your hardware can handle it. However, conversion to long double will silently break (if you indeed need it) outputting NaNs, unless you take steps to put the CPU into a state that faults on these errors.

Thanks again to all who contributed to these changes. This project has produced a working class that I think is quite usable, and I am very encouraged by that.

The new code follows. I'd appreciate your comments and suggestions for improvement. Thank you for your time.

HugeInt.h

/*
 * HugeInt.h
 * 
 * Definition of the huge integer class
 * 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 <iostream>

namespace iota {

class HugeInt {
public:
    HugeInt() = default;
    HugeInt(const 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&); // TODO:

    // increment and decrement operators
    HugeInt& operator+=(const HugeInt&);
    HugeInt& operator-=(const HugeInt&);
    HugeInt& operator*=(const HugeInt&);
   // HugeInt& operator/=(const HugeInt&); TODO:
    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 int      numDigits{300};       // max. no. radix 2^32 digits
    static const uint64_t twoPow32{4294967296}; // 2^32, for convenience
    uint32_t              digits[numDigits]{0}; // radix 2^32 digits; default 0

    // private utility functions
    bool     isZero() const;
    bool     isNegative() const;
    HugeInt& radixComplementSelf();  
    HugeInt  radixComplement() const;
    HugeInt  shortDivide(uint32_t) const;
    uint32_t shortModulo(uint32_t) const;
    HugeInt  shortMultiply(uint32_t) const;
    HugeInt& shiftLeftDigits(int);
};

} /* namespace iota */

#endif /* HUGEINT_H */

HugeInt.cpp

/*
 * HugeInt.cpp
 *
 * Implementation of the HugeInt class. See comments in HugeInt.h for
 * details of representation, etc.
 *
 * February, 2020
 *
 * RADIX 2^32 VERSION
 * 
 */

#include <cstdlib>   // for abs(), labs(), etc.
#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstring>
#include <stdexcept>
#include <cmath>
#include "HugeInt.h"

/*
 * Non-member utility functions (in anonymous namespace -- file scope only).
 * 
 */

namespace { /* anonymous namespace */

/*
 * Simple function to check for non-digit characters in a C string.
 * 
 * Returns true if string contains all digit characters; otherwise
 * false.
 * 
 */

inline bool validate_digits(const char *const str) {
    bool retval = true;

    for (size_t i = 0; i < std::strlen(str); ++i) {
        if (std::isdigit(static_cast<unsigned char>(str[i])) == 0) {
            retval = false;
            break;
        } 
    }

    return retval;
}

/**
 * get_carry32
 *
 * Return the high 32 bits of a 64-bit uint64_t.
 * Return this 32-bit value as a uint32_t.
 * 
 * @param value
 * @return 
 */

inline uint32_t get_carry32(uint64_t value) {
    return static_cast<uint32_t>(value >> 32 & 0xffffffff);
}

/**
 * get_digit32
 * 
 * Return the low 32 bits of the two-byte word stored as an int.
 * Return this 32-bit value as a uint32_t.
 * 
 * @param value
 * @return 
 */

inline uint32_t get_digit32(uint64_t value) {
    return static_cast<uint32_t>(value & 0xffffffff);
}

} /* anonymous namespace */



/*
 * Member functions in namespace iota
 *
 */

namespace iota {

/**
 * Constructor (conversion constructor)
 *
 * Construct a HugeInt from a long long int.
 *
 */ 

HugeInt::HugeInt(const long long int x) {
    if (x == 0LL) {
        return;
    }

    long long int xp{std::llabs(x)};
    int i{0};

    // Successively determine units, 2^32's, (2^32)^2's, (2^32)^3's, etc.
    // storing them in digits[0], digits[1], digits[2], ...,
    // respectively. That is units = digits[0], 2^32's = digits[1], etc.
    while (xp > 0LL) {
        digits[i++] = xp % twoPow32;
        xp /= twoPow32;
    }

    if (x < 0LL) {
        radixComplementSelf();
    }
}

/**
 * Constructor (conversion constructor)
 *
 * Construct a HugeInt from a null-terminated C string representing the
 * base 10 representation of the number. The string is assumed to have 
 * the form "[+/-]31415926", including an optional '+' or '-' sign. 
 *
 * WARNING: No spaces are allowed in the decimal string containing numerals.
 * 
 * 
 * @param str
 */

HugeInt::HugeInt(const char *const str) {
    const int len{static_cast<int>(std::strlen(str))};

    if (len == 0) {
        throw std::invalid_argument{"empty decimal string in constructor."};
    }

    // Check for explicit positive and negative signs and adjust accordingly.
    // If negative, we flag the case and perform a radix complement at the end.
    bool flagNegative{false};
    int  numDecimalDigits{len};
    int  offset{0};

    if (str[0] == '+') {
        --numDecimalDigits;
        ++offset;
    } 

    if (str[0] == '-') {
        flagNegative = true;
        --numDecimalDigits;
        ++offset;
    }

    // validate the string of numerals
    if (!validate_digits(str + offset)) {
        throw std::invalid_argument{"string contains non-digit in constructor."};
    }

    // Loop (backwards) through each decimal digit, digit[i], in the string, 
    // adding its numerical contribution, digit[i]*10^i, to theNumber. Here i 
    // runs upwards from zero, starting at the right-most digit of the string 
    // of decimal digits.
    uint32_t digitValue{0};
    HugeInt  theNumber{0LL};
    HugeInt  powerOfTen{1LL}; // initially 10^0 = 1

    for (int i = 0; i < numDecimalDigits; ++i) {
        digitValue = static_cast<uint32_t>(str[len - 1 - i]) - '0';
        theNumber += powerOfTen.shortMultiply(digitValue);
        powerOfTen = powerOfTen.shortMultiply(10);
    }

    if (flagNegative) {
        theNumber.radixComplementSelf();
    }

    for (int i = 0; i < numDigits; ++i) {
        digits[i] = theNumber.digits[i];
    }
}

/**
 * Copy constructor (could be defaulted)
 * 
 * @param rhs
 */

HugeInt::HugeInt(const HugeInt& rhs) {
    // TODO: perhaps call copy assignment?
    for (int i = 0; i < numDigits; ++i)
        digits[i] = rhs.digits[i];
}

/**
 * Assignment operator
 * 
 * @param rhs
 * @return 
 */

const HugeInt& HugeInt::operator=(const HugeInt& rhs) {
    if (&rhs != this) {
        for (int i = 0; i < numDigits; ++i) {
            digits[i] = rhs.digits[i]; 
        }
    }

    return *this;
}

/**
 * Unary minus operator
 * 
 * @return 
 */

HugeInt HugeInt::operator-() const {
    return radixComplement();
}

/**
 * operator long double() 
 *
 * Use with static_cast<long double>(hugeint) to convert hugeint to its
 * approximate (long double) floating point value.
 * 
 */
HugeInt::operator long double() const {
    long double sign{1.0L};
    HugeInt     copy{*this};

    if (copy.isNegative()) {
        copy.radixComplementSelf();
        sign = -1.0L;
    }

    long double retval{0.0L};
    long double pwrOfBase{1.0L}; // Base = 2^32; (2^32)^0 initially

    for (int i = 0; i < numDigits; ++i) {
        retval += copy.digits[i] * pwrOfBase;
        pwrOfBase *= twoPow32;
    }

    return retval*sign;
}

/**
 * Operator +=
 *
 * NOTE: With the conversion constructors provided, also
 *       provides operator+=(long int) and
 *                operator+=(const char *const)
 * 
 * @param increment
 * @return 
 */

HugeInt& HugeInt::operator+=(const HugeInt& increment) {
    *this = *this + increment;
    return *this;
}

/**
 * Operator -=
 * 
 * NOTE: With the conversion constructors provided, also
 *       provides operator-=(long int) and
 *                operator-=(const char *const)
 * 
 * 
 * @param decrement
 * @return 
 */

HugeInt& HugeInt::operator-=(const HugeInt& decrement) {
    *this = *this - decrement;
    return *this;
}

/**
 * Operator *=
 * 
 * NOTE: With the conversion constructors provided, also
 *       provides operator*=(long int) and
 *                operator*=(const char *const)
 * 
 * @param multiplier
 * @return 
 */

HugeInt& HugeInt::operator*=(const HugeInt& multiplier) {
    *this = *this * multiplier;
    return *this;
}

/**
 * Operator ++ (prefix)
 * 
 * @return 
 */

HugeInt& HugeInt::operator++() {
    *this = *this + 1LL;
    return *this;
}

/**
 * Operator ++ (postfix)
 * 
 * @param 
 * @return 
 */

HugeInt HugeInt::operator++(int) {
   HugeInt retval{*this};
   ++(*this);

   return retval;
}

/**
 * Operator -- (prefix)
 * 
 * @return 
 */

HugeInt& HugeInt::operator--() {
   *this = *this - 1LL;
   return *this;
}

/**
 * Operator -- (postfix)
 * 
 * @param 
 * @return 
 */

HugeInt HugeInt::operator--(int) {
   HugeInt retval{*this};
   --(*this);

   return retval;
}


////////////////////////////////////////////////////////////////////////////
// Input/Output                                                           //
////////////////////////////////////////////////////////////////////////////

/**
 * toRawString()
 * 
 * Format a HugeInt as string in raw internal format, i.e., as a sequence 
 * of base-2^32 digits (each in decimal form, 0 <= digit <= 2^32 - 1).
 *  
 * @return 
 */

std::string HugeInt::toRawString() const {
    int istart{numDigits - 1};

    for ( ; istart >= 0; --istart) {
        if (digits[istart] != 0) {
            break;
        }
    }

    std::ostringstream oss;

    if (istart == -1) // the number is zero
    {
        oss << digits[0];
    } else {
        for (int i = istart; i >= 0; --i) {
            oss << std::setw(10) << std::setfill('0') << digits[i] << " ";
        }
    }

    return oss.str();
}

/**
 * toDecimalString()
 * 
 * Format HugeInt as a string of decimal digits. The length of the decimal 
 * string is estimated (roughly) by solving for x:
 *
 *     (2^32)^N = 10^x    ==>    x = N log_10(2^32) = N * 9.63296 (approx)
 *
 * where N is the number of base 2^32 digits. A safety margin of 5 is added
 * for good measure.
 * 
 * @return 
 */

std::string HugeInt::toDecimalString() const {
    std::ostringstream oss;

    // Special case HugeInt == 0 is easy
    if (isZero()) {
        oss << "0";
        return oss.str();
    }

    // set copy to the absolute value of *this
    // for use in shortDivide and shortModulo
    HugeInt tmp;

    if (isNegative()) {
        oss << "-";
        tmp = this->radixComplement();
    } else {
        tmp = *this;
    }

    // determine the decimal digits of the absolute value 
    int       i{0};
    const int numDecimal{static_cast<int>(numDigits * 9.63296) + 5};
    uint32_t  decimalDigits[numDecimal]{0};

    while (!tmp.isZero()) {
        decimalDigits[i++] = tmp.shortModulo(10);
        tmp = tmp.shortDivide(10);
    }

    // output the decimal digits
    for (int j = i - 1; j >= 0; --j) {
        if (j < i - 1) {
            if ((j + 1) % 3 == 0) {
                oss << ','; // thousands separator
            }
        }

        oss << decimalDigits[j];
    }

    return oss.str();
}

/////////////////////////////////////////////////////////////////////////////
// Useful informational member functions                                   //
/////////////////////////////////////////////////////////////////////////////

/**
 * getMinimum()
 * 
 * Return the minimum representable value for a HugeInt. Static member
 * function.
 * 
 * @return 
 */

HugeInt HugeInt::getMinimum() {
    HugeInt retval;

    retval.digits[numDigits - 1] = 2147483648;

    return retval;
}

/**
 * getMaximum()
 * 
 * Return the maximum representable value for a HugeInt. Static member 
 * function.
 * 
 * @return 
 */

HugeInt HugeInt::getMaximum() {
    HugeInt retval;

    retval.digits[numDigits - 1] = 2147483648;

    --retval;

    return retval;
}

/**
 * numDecimalDigits()
 * 
 * Return the number of decimal digits this HugeInt has. 
 * 
 * We use a simple algorithm using base-10 logarithms. Consider, e.g., 457, 
 * which we can write as 4.57 * 10^2. Taking base-10 logs: 
 *         
 *          log10(4.57 * 10^2) = log10(4.57) + 2.
 * 
 * Since 0 < log10(4.57) < log10(10) = 1, we need to round up (always) to get
 * the extra digit, corresponding to the fractional part in the eq. above. 
 * Hence the use of ceil below. Values of x in the range -10 < x < 10 are dealt
 * with as a special case.
 * 
 * @return 
 */

int HugeInt::numDecimalDigits() const {

    if (-10 < *this && *this < 10) {
        return 1;
    }
    else {
        long double approx = static_cast<long double>(*this);
        return static_cast<int>(std::ceil(std::log10(std::fabs(approx))));
    }
}

////////////////////////////////////////////////////////////////////////////
// friend functions                                                       //
////////////////////////////////////////////////////////////////////////////

/**
 * friend binary operator +
 *
 * Add two HugeInts a and b and return c = a + b.
 *
 * Note: since we provide a conversion constructor for long long int's, this 
 *       function, in effect, also provides the following functionality by 
 *       implicit conversion of long long int's to HugeInt
 *
 *       c = a + <some long long int>    e.g.  c = a + 2412356LL
 *       c = <some long long int> + a    e.g.  c = 2412356LL + a
 * 
 * @param a
 * @param b
 * @return 
 */

HugeInt operator+(const HugeInt& a, const HugeInt& b) {
    uint32_t carry{0};
    uint64_t partial{0};
    HugeInt sum;

    for (int i = 0; i < HugeInt::numDigits; ++i) {
        partial =   static_cast<uint64_t>(a.digits[i]) 
                  + static_cast<uint64_t>(b.digits[i]) 
                  + static_cast<uint64_t>(carry);
        carry = get_carry32(partial);
        sum.digits[i] = get_digit32(partial);
    }

    return sum;
}

/**
 * friend binary operator-
 *
 * Subtract HugeInt a from HugeInt a and return the value c = a - b.
 *
 * Note: since we provide a conversion constructor for long long int's, this 
 *       function, in effect, also provides the following functionality by 
 *       implicit conversion of long long int's to HugeInt:
 *
 *       c = a - <some long long int>    e.g.  c = a - 2412356LL
 *       c = <some long long int> - a    e.g.  c = 2412356LL - a
 * 
 * @param a
 * @param b
 * @return 
 */

HugeInt operator-(const HugeInt& a, const HugeInt& b) {
    return a + (-b);
}

/**
 * friend binary operator *
 *
 * Multiply two HugeInt numbers. Uses standard long multipication algorithm
 * adapted to base 2^32. See comments on implicit conversion before 
 * HugeInt operator+(const HugeInt&, const HugeInt& ) above.
 * 
 * @param a
 * @param b
 * @return 
 */

HugeInt operator*(const HugeInt& a, const HugeInt& b) {
    HugeInt product{0LL};
    HugeInt partial;

    for (int i = 0; i < HugeInt::numDigits; ++i) {
        partial = a.shortMultiply(b.digits[i]);
        product += partial.shiftLeftDigits(i);
    }

    return product;
}

////////////////////////////////////////////////////////////////////////////
// Relational operators (friends)                                         //
////////////////////////////////////////////////////////////////////////////

/**
 * Operator ==
 * 
 * @param lhs
 * @param rhs
 * @return 
 */

bool operator==(const HugeInt& lhs, const HugeInt& rhs) {
   HugeInt diff{rhs - lhs};

   return diff.isZero();
}

/**
 * Operator !=
 * 
 * @param lhs
 * @param rhs
 * @return 
 */

bool operator!=(const HugeInt& lhs, const HugeInt& rhs) {
   return !(rhs == lhs);
}

/**
 * Operator <
 * 
 * @param lhs
 * @param rhs
 * @return 
 */

bool operator<(const HugeInt& lhs, const HugeInt& rhs) {
   HugeInt diff{lhs - rhs};

   return diff.isNegative();
}

/**
 * Operator >
 * 
 * @param lhs
 * @param rhs
 * @return 
 */

bool operator>(const HugeInt& lhs, const HugeInt& rhs) {
    return rhs < lhs;
}

/**
 * Operator <=
 * 
 * @param lhs
 * @param rhs
 * @return 
 */

bool operator<=(const HugeInt& lhs, const HugeInt& rhs) {
    return !(lhs > rhs);
}

/**
 * Operator >=
 * 
 * @param lhs
 * @param rhs
 * @return 
 */

bool operator>=(const HugeInt& lhs, const HugeInt& rhs) {
    return !(lhs < rhs);
}


////////////////////////////////////////////////////////////////////////////
// Private utility functions                                              //
////////////////////////////////////////////////////////////////////////////

/**
 * isZero()
 * 
 * Return true if the HugeInt is zero, otherwise false.
 * 
 * @return 
 */

bool HugeInt::isZero() const {
    int i{numDigits - 1};

    for ( ; i >= 0; --i) {
        if (digits[i] != 0) {
            break;
        }
    }

    return i == -1;
}

/**
 * isNegative()
 * 
 * Return true if a number x is negative (x < 0). If x >=0, then
 * return false.
 * 
 * NOTE: In the radix-2^32 complement convention, negative numbers, x, are 
 *       represented by the range of values: (2^32)^N/2 <= x <=(2^32)^N - 1.
 *       Since (2^32)^N/2 = (2^32/2)*(2^32)^(N-1) = 2147483648*(2^32)^(N-1), 
 *       we need only check whether the (N - 1)'th base 2^32 digit is at 
 *       least 2147483648. 
 * 
 * @return 
 */

bool HugeInt::isNegative() const {
    return digits[numDigits - 1] >= 2147483648;
}

/**
 * shortDivide:
 * 
 * Return the result of a base 2^32 short division by divisor, where 
 * 0 < divisor <= 2^32 - 1, using the usual primary school algorithm 
 * adapted to radix 2^32.
 *
 * WARNING: assumes both HugeInt and the divisor are POSITIVE.
 * 
 * @param divisor
 * @return 
 */

HugeInt HugeInt::shortDivide(uint32_t divisor) const {
    uint64_t j;
    uint64_t remainder{0};
    HugeInt quotient;

    for (int i = numDigits - 1; i >= 0; --i) {
        j = twoPow32 * remainder + static_cast<uint64_t>(digits[i]);
        quotient.digits[i] = static_cast<uint32_t>(j / divisor);
        remainder = j % divisor;
    }

    return quotient;
}

/**
 * shortModulo
 * 
 * Return the remainder of a base 2^32 short division by divisor, where 
 * 0 < divisor <= 2^32 - 1.
 *
 * WARNING: assumes both HugeInt and the divisor are POSITIVE.
 * 
 * @param divisor
 * @return 
 */

uint32_t HugeInt::shortModulo(uint32_t divisor) const {
    uint64_t j;
    uint64_t remainder{0};

    for (int i = numDigits - 1; i >= 0; --i) {
        j = twoPow32 * remainder + static_cast<uint64_t>(digits[i]);
        remainder = j % divisor;
    }

    return static_cast<uint32_t>(remainder);
}

/**
 * shortMultiply
 * 
 * Return the result of a base 2^32 short multiplication by multiplier, where
 * 0 <= multiplier <= 2^32 - 1.
 *
 * WARNING: assumes both HugeInt and multiplier are POSITIVE.
 * 
 * @param multiplier
 * @return 
 */

HugeInt HugeInt::shortMultiply(uint32_t multiplier) const {
    uint32_t carry{0};
    uint64_t tmp;
    HugeInt product;

    for (int i = 0; i < numDigits; ++i) {
        tmp = static_cast<uint64_t>(digits[i]) * multiplier + carry;
        carry = get_carry32(tmp);
        product.digits[i] = get_digit32(tmp);
    }

    return product;
}

/**
 * shiftLeftDigits
 *
 * Shift this HugeInt's radix-2^32 digits left by num places, filling
 * with zeroes from the right.
 * 
 * @param num
 * @return 
 */

HugeInt& HugeInt::shiftLeftDigits(int num) {
    if (num == 0) {
        return *this;
    }

    for (int i = numDigits - num - 1; i >= 0; --i) {
        digits[i + num] = digits[i];
    }

    for (int i = 0; i < num; ++i) {
        digits[i] = 0;
    }

    return *this;
}

/**
 * radixComplementSelf()
 *
 * Perform a radix complement on the object in place (changes object).
 * 
 * @return 
 */

HugeInt& HugeInt::radixComplementSelf() {
    if (!isZero()) {
        uint64_t sum{0};
        uint32_t carry{1};

        for (int i = 0; i < numDigits; ++i) {
            sum =   static_cast<uint64_t>(twoPow32 - 1) 
                  - static_cast<uint64_t>(digits[i]) 
                  + static_cast<uint64_t>(carry);
            carry = get_carry32(sum);
            digits[i] = get_digit32(sum);
        }
    }

    return *this;
}

/**
 * radixComplement()
 * 
 * Return the radix-2^32 (base-2^32) complement of HugeInt.
 * 
 * @return 
 */

HugeInt HugeInt::radixComplement() const {
    HugeInt result{*this};

    return result.radixComplementSelf();
}

/**
 * operator<<
 * 
 * Overloaded stream insertion for HugeInt.
 * 
 * @param output
 * @param x
 * @return 
 */

std::ostream& operator<<(std::ostream& output, const HugeInt& x) {
    output << x.toDecimalString();

    return output;
}

/**
 * operator >>
 * 
 * Overloaded stream extraction for HugeInt.
 * 
 * @param input
 * @param x
 * @return 
 */

std::istream& operator>>(std::istream& input, HugeInt& x) {
    std::string str;

    input >> str;
    x = HugeInt(str.c_str());

    return input;
}

} /* namespace iota */

Sample driver code:

/*
 * Simple driver to test a few features of the HugeInt class.
 * 
 * Improved version.
 * 
 */

#include <iostream>
#include <iomanip>
#include <limits>
#include "HugeInt.h"

iota::HugeInt read_bounded_hugeint(const iota::HugeInt&, const iota::HugeInt&);
iota::HugeInt factorial_recursive(const iota::HugeInt&);
iota::HugeInt factorial_iterative(const iota::HugeInt&);
iota::HugeInt fibonacci_recursive(const iota::HugeInt&);
iota::HugeInt fibonacci_iterative(const iota::HugeInt&);
void preamble();

// limits to avoid overflow
const iota::HugeInt FACTORIAL_LIMIT{1100LL};
const iota::HugeInt FIBONACCI_LIMIT{13000LL};

int main() {

    preamble(); // blah

    iota::HugeInt nfac = read_bounded_hugeint(0LL, FACTORIAL_LIMIT);

    iota::HugeInt factorial = factorial_iterative(nfac);
    long double factorial_dec = static_cast<long double>(factorial);

    std::cout << "\nThe value of " << nfac << "! is:\n";
    std::cout << factorial << '\n';
    std::cout << "\nThis value has " << factorial.numDecimalDigits() 
              << " decimal digits.\n";
    std::cout.precision(std::numeric_limits<long double>::digits10);
    std::cout << "\nIts decimal approximation is: " << factorial_dec << "\n\n";


    iota::HugeInt nfib = read_bounded_hugeint(0LL, FIBONACCI_LIMIT);

    iota::HugeInt fibonacci = fibonacci_iterative(nfib);
    long double fibonacci_dec = static_cast<long double>(fibonacci);

    std::cout << "\nThe " << nfib << "th Fibonacci number is:\n";
    std::cout << fibonacci << '\n';
    std::cout << "\nThis value has " << fibonacci.numDecimalDigits() 
              << " decimal digits.\n";
    std::cout << "\nIts decimal approximation is: " << fibonacci_dec << '\n';

    std::cout << "\nComparing these two values we observe that ";
    if (factorial == fibonacci) {
        std::cout << nfac << "! == Fibonacci_{" << nfib << "}\n";
    }

    if (factorial < fibonacci) {
        std::cout << nfac << "! < Fibonacci_{" << nfib << "}\n";
    }

    if (factorial > fibonacci) {
        std::cout << nfac << "! > Fibonacci_{" << nfib << "}\n";
    }

    iota::HugeInt sum = factorial + fibonacci;
    iota::HugeInt diff = factorial - fibonacci;

    std::cout << "\nTheir sum (factorial + fibonacci) is:\n";
    std::cout << sum << '\n';
    std::cout << "\n\twhich is approximately " << static_cast<long double>(sum);
    std::cout << '\n';

    std::cout << "\nTheir difference (factorial - fibonacci) is:\n";
    std::cout << diff << '\n';
    std::cout << "\n\twhich is approximately " << static_cast<long double>(diff);
    std::cout << '\n';

    iota::HugeInt x{"-80538738812075974"};
    iota::HugeInt y{"80435758145817515"};
    iota::HugeInt z{"12602123297335631"};

    iota::HugeInt k = x*x*x + y*y*y + z*z*z;

    std::cout << "\nDid you know that, with:\n";
    std::cout << "\tx = " << x << '\n';
    std::cout << "\ty = " << y << '\n';
    std::cout << "\tz = " << z << '\n';
    std::cout << "\nx^3 + y^3 + z^3 = " << k << '\n';

    return 0;
}


/**
 * read_bounded_hugeint
 * 
 * Read and return a value in the range min <= value <= max.
 * Dies after 5 retries.
 * 
 * @param min
 * @param max
 * @return 
 */
iota::HugeInt read_bounded_hugeint(const iota::HugeInt& min, 
                                   const iota::HugeInt& max) {
    iota::HugeInt value;
    bool fail;
    int retries = 0;

    do {
        try {
            std::cout << "Enter an integer (" << min << " - " << max << "): "; 
            std::cin >> value;

            if (value < min || value > max) {
                fail = true;
                ++retries;
            }
            else {
                fail = false;
            }
        }
        catch (std::invalid_argument& error) {
            std::cout << "You entered an invalid HugeInt value.";
            std::cout << " Please use, e.g., [+/-]1234567876376763.\n";
            //std::cout << "Exception: " << error.what() << '\n';
            fail = true;
            ++retries;
        }  
    } while (fail && retries < 5);

    if (retries == 5) {
        std::cerr << "Giving up...\n";
        exit(EXIT_FAILURE);
    }

    return value;
}

/**
 * factorial_recursive:
 * 
 * Recursive factorial function using HugeInt. Not too slow.
 * 
 * @param n
 * @return 
 */

iota::HugeInt factorial_recursive(const iota::HugeInt& n) {
    const iota::HugeInt one{1LL};

    if (n <= one) {
        return one;
    } else {
        return n * factorial_recursive(n - one);
    }
}

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;
}

/**
 * fibonacci_recursive:
 * 
 * Recursively calculate the n'th Fibonacci number, where n>=0.
 * 
 * WARNING: S l o w . . .
 * 
 * @param n
 * @return 
 */
iota::HugeInt fibonacci_recursive(const iota::HugeInt& n) {
    const iota::HugeInt zero;
    const iota::HugeInt one{1LL};

    if ((n == zero) || (n == one)) {
        return n;
    } 
    else {
        return fibonacci_recursive(n - 1LL) + fibonacci_recursive(n - 2LL);
    }  
}

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;
}

void preamble() {
    long double min = static_cast<long double>(iota::HugeInt::getMinimum());
    long double max = static_cast<long double>(iota::HugeInt::getMaximum());

    std::cout.precision(std::numeric_limits<long double>::digits10);
    std::cout <<"**************************************************************"
              <<"*************\n\n";
    std::cout << "The range of integers, x, that can be represented in " 
              << "the default HugeInt\nconfiguration is, approximately\n"
              << "      " << min << " <= x <= " << max << '\n';

    std::cout << "\nThe precise values of the upper and lower limits can be "
              << "found using\nHugeInt::getMinimum()/HugeInt::getMaximum().\n";

    std::cout << "\nThe maximum number of decimal digits of an integer "
              << "representable with\na HugeInt is: " 
              << iota::HugeInt::getMaximum().numDecimalDigits()
              << "\n\n";   
    std::cout <<"**************************************************************"
              <<"*************\n\n";
}

Answer 1

I recommend including your internal headers first, before any standard headers. This helps expose any accidental dependencies that make it hard to use your types in another program.

So in HugeInt.cpp, I'd write

#include "HugeInt.h"

#include <cstdlib>   // for abs(), labs(), etc.
#include <iostream>
#include <iomanip>
#include <sstream>
#include <cstring>
#include <stdexcept>
#include <cmath>

---

Don't include <iostream> in the header; that's overkill. We can include <iosfwd> instead, which is much lighter, and just gives us the forward declarations we need to specify std::istream& and std::ostream& in our signatures, without bringing in the full template definitions that are in <iostream>. That means that translation units not doing I/O don't carry the overhead.

---

Spelling (throughout): uint32_t and uint64_t are in the std namespace. Your compiler is allowed to also declare them in the global namespace, but is not required to, so you have a portability bug.

---

Perhaps we should provide a divmod() function that gives quotient and remainder in one operation, then use that to implement division and modulo? The "short" version of this would certainly be useful, and save duplication when executing toDecimalString().

---

We use a string-stream to implement toDecimalString(), and then use that to implement operator<<(). I think we should do that the other way around: use operator<<() to implement toDecimalString(). Then streamed output wouldn't need to create a temporary string. Consider using my stream-saver if you want to preserve its manipulator state.

---

In toDecimalString(), we have this loop that produces one digit at a time:

while (!tmp.isZero()) {
    decimalDigits[i++] = tmp.shortModulo(10);
    tmp = tmp.shortDivide(10);
}

We can reduce the number of divisions ninefold, since shortModulo() accepts std::uint32_t. We can divide by 1 billion rather than ten, if we're careful about how we print the first "block":

// determine the decimal digits of the absolute value
constexpr int numDecimal{static_cast<int>(numDigits * 1.07032) + 1};
std::array<std::uint32_t, numDecimal> decimal;

int i{0};
while (!tmp.isZero()) {
    decimal[i++] = tmp.shortModulo(1'000'000'000);
    tmp = tmp.shortDivide(1'000'000'000);
}

// output the decimal digits, in threes
oss << std::setw(0) << std::setfill('0');
--i;
// first digits
auto const d0 = decimal[i];
if (d0 > 1'000'000) {
    oss << (d0 / 1'000'000) << ',' << std::setw(3);
}
if (d0 > 1'000) {
    oss << (d0 / 1'000 % 1'000) << ',' << std::setw(3);
}
oss << (d0 % 1'000);
// subsequent digits
while (i--) {
    auto const d = decimal[i];
    oss << ',' << std::setw(3) << (d / 1'000'000)
        << ',' << std::setw(3) << (d / 1'000 % 1'000)
        << ',' << std::setw(3) << (d % 1'000);
}

Answer 2

API

• In HugeInt(long long int), there's no need for the const. In the function declaration it is ignored, at least by the GNU C++ compiler.
• Instead of isZero, I prefer to write == 0. I'd rather implement two == operators that compare HugeInt, long long int and vice versa.
• Same for isNegative, I'd write that as < 0.
• I'd also write huge_int % 100 instead of calling huge_int.shortModulo(100), as the former is much shorter.
• Same for the other mathematical operators.
• Instead of spelling out every digit of 4294967296, why don't you just write uint64_t(1) << 32? That would make the comment redundant.
• I'd rename twoPow32 to base or radix.
• The name radixComplementSelf sounds like a bloated description for negate. It should probably be renamed to the latter.

Implementation

• validate_digits should rather be called is_all_digits since it doesn't validate anything, it just tests without rejecting.
• In validate_digits, the retval variable can be removed if you just return from within the loop.
• Still in validate_digits, what answer do you expect for the empty string, true or false?
• In get_carry32 the & 0xffffffff is unnecessary since uint32_t is an exact-width integer type.
• Same for get_digit32.
• In get_digit32, the comment "two-byte word" probably refers to the old version that had uint8_t as its underlying type.
• Is there a particular reason for using x == 0LL instead of the simpler x == 0? For integer literals I find the latter easier to read, but of course your tools for static analysis may disagree.
• std::llabs(x) looks a lot like it would invoke undefined behavior for LLONG_MIN.
• In the const char * constructor, len should be of type std::size_t instead of int.
• Still in the const char * constructor, the value of uint32_t digitValue{0} is unused. This variable should be moved inside the for loop, to make its scope as small as possible.
• Still in the const char * constructor, I think it's more efficient to process the digits from left to right, using Horner's method. This could save you both of the temporary HugeInt variables. It certainly feels strange to have two variables of the same type that you are currently constructing, as if there were some hidden recursion.
• In the copy constructor, I wouldn't call the parameter rhs since there is no corresponding lhs variable anywhere nearby.
• In the operator=, I would leave out the &rhs == this test. It's unlikely that a variable gets assigned to itself.
• In operator long double, you should add the missing spaces in the retval*sign expression. The rest looks simple and fine.
• Instead of defining operator+= in terms of operator+, the more efficient variant is the other way round: first implement operator+= as a basic operation and then define operator+ using +=. It typically saves a few variable assignments and memory allocations.
• I wonder whether you should really implement operator++(int) and operator--(int). I don't think they are needed often. I'd wait until I really need them, just out of curiosity.
• In operator-- you should rather call operator-= instead of operator-, since that will be more efficient after you rewrote the code.
• Instead of the currently unused HugeInt::toRawString I'd rather provide HugeInt::toHexString since that can be implemented easier and doesn't need any padding. Did you find any use for the current toRawString?
• In toDecimalString, in the case of zero there is no need to allocate an std::ostringstream. You can just write: if (isZero()) return "0";.
• Still in toDecimalString, the variable name tmp is always bad. Temporary, sure, but temporary what?
• In getMinimum and getMaximum, you should rather write uint32_t(1) << 31 instead of spelling out the digits. It's simpler to read.
• The comment above operator!= contains exactly zero useful characters. It should be removed, and the other comments as well.
• In isZero, it's a waste of time that you have to iterate through all the 300 digits just to check if any of them is nonzero. It would be far more efficient to have a size_t len field in the HugeInt. Then you would only need to check whether len <= 1 && digits[0] == 0.
• In shortModulo, in addition to the "WARNING: assumes" comment, you should also add assertions to the code: assert(divisor > 0) and assert(*this >= 0).

Test

In addition to the main function with an interactive test session, you should have lots of automatic unit tests that cover all normal and edge cases. Having these tests and running them regularly provides a safety net for all kinds of refactorings.

Performance

Calculating the factorial of 1000 and printing it should be possible in less than a millisecond. For example, the following quickly-written Kotlin code runs in 5 seconds on my machine, which means 0.5 milliseconds per iteration:

class BigIntTest {
    @Test
    fun name() {
        fun fac(n: Int): BigInteger {
            var fac = BigInteger.valueOf(1)
            for (i in 1..n) {
                fac *= BigInteger.valueOf(i.toLong())
            }
            return fac
        }

        for (i in 1..10000)
            fac(1000).toString()
    }
}

A similar Go program is even faster, it needs only 0.137 milliseconds per iteration:

package main

import (
    "math/big"
    "testing"
)

func fac(n int) *big.Int {
    fac := big.NewInt(1)
    for i := 1; i <= n; i++ {
        fac.Mul(fac, big.NewInt(int64(i)))
    }
    return fac
}

var result string

func Benchmark(b *testing.B) {
    for i := 0; i < b.N; i++ {
        result = fac(1000).String()
    }
}

Your code should become similarly fast.

Answer 3

The multiplication approach could be improved, even without non-portable optimization (significant gains are possible from that too, for example using _mulx_u64). Sadly it would make the code for multiplication less pretty, it looks nice now and with the following approach it wouldn't look as nice.

A not so nice thing that is happening here is that the order of the computation forces the creation of big partial products which are added one by one, each taking an explict shift and a full BigInt addition. For the factorial benchmark this is not a concern, as one operand is always tiny. For a benchmark that stresses "balanced" multiplication, you could for example extend the Fibonacci tests to exponentiating a 2x2 matrix.

An alternative arrangement is to generate the small partial products (one limb of the multiplier times one limb of the multiplicand) in the order of their weight, so a group of partial products (and a carry) of equal weight can be summed immediately without lots of temporary storage and then result in a limb of the final result. Handling carries is a bit tricky.

For clarity, here is a diagram of that ordering:

(source: cryptographic hardware and embedded systems)