Calculate factorial of big number

factorial.cc

// Compile:
//
// - Linux: g++ -O2 -pthread factorial.cc
// - Mac: clang++ -O2 -std=c++11 factorial.cc

#include <vector>
#include <string>
#include <algorithm>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <thread>

template <typename NumType, NumType DigitSize, NumType DigitWidth>
class BigInt
{
public:
    typedef BigInt<NumType, DigitSize, DigitWidth> SelfType;


    BigInt() = default;
    BigInt(const SelfType& x) = default;

    explicit BigInt(const std::string& x)
    {
        size_t i = x.size();
        for(; i >= DigitWidth ; i -= DigitWidth)
        {
            Data.push_back(0);
            std::istringstream Converter(x.substr(i - DigitWidth, DigitWidth));
            Converter >> Data[Data.size() - 1];
        }

        if(i > 0)
        {
            Data.push_back(0);
            std::istringstream Converter(x.substr(0, i));
            Converter >> Data[Data.size() - 1];
        }
    }

    ~BigInt() = default;

    // static const NumType DigitSize = 100000000;
    // static const NumType DigitWidth = 8;

    void reset()
    {
        Data.clear();
        Data.push_back(0);
    }

    SelfType& operator+=(const SelfType& rhs)
    {
        addWithShift(rhs, 0);
        return *this;
    }

    SelfType& operator*=(const SelfType& rhs)
    {
        Data.reserve(Data.size() + rhs.Data.size());
        SelfType Orig(*this);
        Orig.Data.reserve(Data.size() + 1);
        reset();

        for(size_t i = 0; i < rhs.Data.size(); ++i)
        {
            SelfType Iter(Orig);
            Iter *= rhs.Data[i];
            addWithShift(Iter, i);
        }
        return *this;
    }

    static SelfType factorial(NumType n)
    {
        SelfType One(1);
        SelfType Result(1);
        SelfType x(1);
        Result.Data.reserve(1000000);
        for(NumType i = 0; i < n; ++i)
        {
            if(i % 1000 == 0)
            {
                std::cerr << ".";
                std::cerr.flush();
            }

            Result *= x;
            x += One;
        }

        return Result;
    }

    static SelfType factorialParallel(NumType n, size_t thread_count)
    {
        std::vector<std::thread> Threads;
        std::vector<SelfType> Results;
        Results.resize(thread_count);

        for(size_t i = 0; i < thread_count; i++)
        {
            Threads.push_back(std::thread(factorialThread, n, i, thread_count,
                                          &(Results[i])));
        }

        for(size_t i = 0; i < thread_count; i++)
        {
            Threads[i].join();
        }

        for(size_t i = 1; i < thread_count; i++)
        {
            std::cerr << "Multiplying " << i << "th result..." << std::endl;
            Results[0] *= Results[i];
        }

        return Results[0];
    }

private:
    explicit BigInt(NumType x)
    {
        if(x == 0)
        {
            Data.push_back(0);
        }
        else
        {
            while(x > 0)
            {
                Data.push_back(x % DigitSize);
                x /= DigitSize;
            }
        }
    }

    std::vector<NumType> Data;

    static void normalizeDigit(NumType& digit, NumType& extra)
    {
        if(digit >= DigitSize)
        {
            extra = digit / DigitSize;
            digit = digit % DigitSize;
        }
        else
        {
            extra = 0;
        }
    }

    NumType operator[] (const size_t i) const
    {
        if(i >= Data.size())
        {
            return 0;
        }
        else
        {
            return Data[i];
        }
    }

    // Lhs += rhs * DigitSize^shift
    SelfType& addWithShift(const SelfType& rhs, size_t shift)
    {
        NumType Extra = 0;
        Data.reserve(std::max(Data.size(), rhs.Data.size() + shift) + 1);
        Data.resize(std::max(Data.size(), rhs.Data.size() + shift));
        size_t i = shift;

        for(; i < std::min(Data.size(), rhs.Data.size() + shift); ++i)
        {
            Data[i] = Data[i] + rhs.Data[i - shift] + Extra;
            normalizeDigit(Data[i], Extra);
        }

        if(rhs.Data.size() + shift <= Data.size())
        {
            while(Extra > 0)
            {
                if(i >= Data.size())
                {
                    Data.push_back(0);
                }

                Data[i] = Data[i] + Extra;
                normalizeDigit(Data[i], Extra);
                ++i;
            }
        }
        else
        {
            for(; i < rhs.Data.size() + shift; ++i)
            {
                Data.push_back(0);
                Data[i] = rhs.Data[i - shift] + Extra;
                normalizeDigit(Data[i], Extra);
            }

            while(Extra > 0)
            {
                if(i >= Data.size())
                {
                    Data.push_back(0);
                }

                Data[i] = Extra;
                normalizeDigit(Data[i], Extra);
                ++i;
            }
        }

        return *this;
    }

    SelfType& operator*=(NumType rhs)
    {
        if(rhs == 1)
        {
            return *this;
        }

        if(rhs == 0)
        {
            Data.clear();
            Data.push_back(0);
            return *this;
        }

        size_t i = 0;
        NumType Extra = 0;

        for(; i < Data.size(); ++i)
        {
            Data[i] = Data[i] * rhs + Extra;
            normalizeDigit(Data[i], Extra);
        }

        if(Extra > 0)
        {
            Data.push_back(Extra);
        }

        return *this;
    }

    static void factorialThreadReversed(
        NumType n, size_t thread_idx, size_t total_threads,
        SelfType* result)
    {
        SelfType Delta(total_threads);
        (*result) = SelfType(1);
        // SelfType x(thread_idx + 1);
        result -> Data.reserve(1000000);

        for(NumType i = n; i --> 0 ;)
        {
            if(i % total_threads != thread_idx)
            {
                continue;
            }

            SelfType x(i + 1);
            if(i % 1000 == 0)
            {
                std::cerr << ".";
                std::cerr.flush();
            }

            (*result) *= x;
        }
    }

    static void factorialThread(
        NumType n, size_t thread_idx, size_t total_threads,
        SelfType* result)
    {
        SelfType Delta(total_threads);
        (*result) = SelfType(1);
        SelfType x(thread_idx + 1);
        result -> Data.reserve(1000000);

        for(NumType i = thread_idx; i < n; i += total_threads)
        {
            if(i % total_threads != thread_idx)
            {
                continue;
            }

            if(i % 1000 == 0)
            {
                std::cerr << ".";
                std::cerr.flush();
            }

            (*result) *= x;
            x += Delta;
        }
    }

    template <typename NumType_, NumType_ DigitSize_, NumType_ DigitWidth_>
    friend std::ostream& operator <<(
        std::ostream& out, const BigInt<NumType_, DigitSize_, DigitWidth_>& v);

    friend void test();
};

template<typename NumType, NumType DigitSize, NumType DigitWidth>
std::ostream& operator<<(std::ostream& os, const BigInt<NumType, DigitSize, DigitWidth>& x)
{
    // std::ios Format(nullptr);
    // Format.copyfmt(os);

    size_t i = x.Data.size() - 1;
    os << x.Data[i];

    for(; i --> 0 ;)
    {
        os << std::setfill('0')
           << std::setw(DigitWidth)
           << x.Data[i];
    }

    // os.copyfmt(Format);
    return os;
}

void test()
{
    typedef BigInt<unsigned int, 1000, 3> TheInt;
    auto x = TheInt("9999999");
    unsigned int y = 999;
    auto z = TheInt("999");
    std::cout << x << ", " << y << ", " << z << std::endl;

    x *= y;
    std::cout << x << std::endl;
    x += z;
    std::cout << x << std::endl;
    x.addWithShift(z, 0);
    std::cout << x << std::endl;
    x.addWithShift(z, 1);
    std::cout << x << std::endl;
    x.addWithShift(z, 2);
    std::cout << x << std::endl;
    x.addWithShift(z, 10);
    std::cout << x << std::endl;
    x *= z;
    std::cout << x << std::endl;

    auto a = TheInt("9999999");
    auto b = TheInt("99999999999999999999999999999999999999999999");
    a *= b;
    std::cout << a << std::endl;
}

int main()
{
    typedef BigInt<uint64_t, 100000000, 8> TheInt;
    std::cout << TheInt::factorialParallel(200000, 2) << std::endl;
    return 0;
}