3
\$\begingroup\$

I am developing a C++ big integer library for C++11, 14, 17, 20, etc., as an exercise.

The Toom-3 algorithm is resisting me a little bit. AFAICT, it does produce the right results but not as fast as I expected.

My questions are:

  • Do you see any way to improve performance?
    (Implementing Toom 2.5 was suggested by user greybeard in the now deleted comments, and I will do that for sure).
  • Do you have any other comments (style, naming of variables, etc)?

Multiplication functions

Here are the 2 functions that, I think, matter to explain my problem:

  • mult is in charge of picking what it thinks is the most suitable algorithm given the operands' size; 3 algorithms are implemented so far (if you are unfamiliar with them, you can get some introduction in the documentation I dropped on the linked repository).

Note that I have set a high threshold. At such value, the benefit of calling mult_toom3 over mult_karatsuba should be visible even if called only once, and then never again in the recursion.

  • mult_toom3 performs Toom-3 calculation, calling mult recursively.

In the following mult_toom3() implementation, I have added some time measures (when //#define MeasureTime is uncommented) to show where processing time is being spent. Almost all of it is in the recursion (and none in the calculations around).

void bigint::bigint_t::mult(
    limbs::iterator       product, limbs::iterator       product_end,
    limbs::const_iterator a_begin, limbs::const_iterator a_end,
    limbs::const_iterator b_begin, limbs::const_iterator b_end
) {
    constexpr size_t karatsubaThreshold =  128,
                     toom3Threshold     = 2000;

    ptrdiff_t a_len    = std::distance(a_begin, a_end),
              b_len    = std::distance(b_begin, b_end),
              calc_len = std::min(a_len, b_len);

    if (calc_len >= karatsubaThreshold && b_len * 3 < a_len * 2) {
    // Unbalanced operand case (a >> b) -> Split the longer one to better use algorithms
        auto a_mid = a_begin + b_len;
        auto product_mid = product + (2 * b_len);
        mult(product, product_mid, b_begin, b_end, a_begin, a_mid);
        limbs tailResult; tailResult.resize(a_len, 0);
        mult(tailResult.begin(), tailResult.end(), b_begin, b_end, a_mid, a_end);
        addInto(product + b_len, tailResult.cbegin(), tailResult.cend());
    }
    else if (calc_len >= karatsubaThreshold && a_len * 3 < b_len * 2) {
    // Unbalanced operand case (a << b) -> Split the longer one to better use algorithms
        auto b_mid = b_begin + a_len;
        auto product_mid = product + (2 * a_len);
        mult(product, product_mid, a_begin, a_end, b_begin, b_mid);
        limbs tailResult; tailResult.resize(b_len, 0);
        mult(tailResult.begin(), tailResult.end(), a_begin, a_end, b_mid, b_end);
        addInto(product + a_len, tailResult.cbegin(), tailResult.cend());
    }
    else if (calc_len >= toom3Threshold)
        mult_toom3(product, product_end, a_begin, a_end, b_begin, b_end);
    else if (calc_len >= karatsubaThreshold)
        mult_karatsuba(product, product_end, a_begin, a_end, b_begin, b_end);
    else
        mult_vanilla(product, product_end, a_begin, a_end, b_begin, b_end);
}

// See algorithm at https://github.com/AtmoFX/bigint/blob/master/documentation/multiplication.md#Toom-Cook
// Notation as per the documentation: lowercase variable = strings of limbs; uppercase variable = polynomial
void bigint::bigint_t::mult_toom3(
    limbs::iterator       product_begin, limbs::iterator       product_end,
    limbs::const_iterator a_begin      , limbs::const_iterator a_end,
    limbs::const_iterator b_begin      , limbs::const_iterator b_end
)
{
//#define MeasureTime
#ifdef MeasureTime
    using std::chrono::steady_clock;
    using std::chrono::duration_cast;
    using std::chrono::duration;
    using std::chrono::milliseconds;

    auto milestone1 = steady_clock::now();
#endif
    size_t size_a  = std::distance(a_begin, a_end),
           size_b  = std::distance(b_begin, b_end),
           l       = std::min(size_a, size_b),
           l_third = (l + 2) / 3;

    limbs::const_iterator
        a1_begin = a_begin  + l_third,
        a2_begin = a1_begin + l_third,
        b1_begin = b_begin  + l_third,
        b2_begin = b1_begin + l_third;
    limbs::iterator
        p0_begin = product_begin,
        p1_begin = p0_begin + l_third,
        p2_begin = p1_begin + l_third,
        p3_begin = p2_begin + l_third,
        p4_begin = p3_begin + l_third;

    bigint_t a0, a1, a2, b0, b1, b2;
    a0.value.resize(l_third); a1.value.resize(l_third); a2.value.resize(size_a - 2 * l_third);
    b0.value.resize(l_third); b1.value.resize(l_third); b2.value.resize(size_b - 2 * l_third);
    std::copy(a_begin, a1_begin, a0.value.begin()); std::copy(a1_begin, a2_begin, a1.value.begin()); std::copy(a2_begin, a_end, a2.value.begin());
    std::copy(b_begin, b1_begin, b0.value.begin()); std::copy(b1_begin, b2_begin, b1.value.begin()); std::copy(b2_begin, b_end, b2.value.begin());

#ifdef MeasureTime
    auto milestone2 = steady_clock::now();
#endif

    bigint_t va_1 = a2 + a0 - a1              , vb_1 = b2 + b0 - b1,
             va1  = a2 + a1 + a0              , vb1 = b2 + b1 + b0,
             va2  = (a2 << 2) + (a1 << 1) + a0, vb2 = (b2 << 2) + (b1 << 1) + b0;

#ifdef MeasureTime
    auto milestone3 = steady_clock::now();
#endif

    //P(-1), P(1) and P(2)
    bigint_t P_1 = va_1 * vb_1,
             P1  = va1  * vb1,
             P2  = va2  * vb2;


    // p0 = P(0) and p4 = P(infinity)
    bigint_t p0  = a0 * b0,
             p4  = a2 * b2;

#ifdef MeasureTime
    auto milestone4 = steady_clock::now();
#endif

    bigint_t p0_plus_p4 = p0 + p4,

             p2 = ((P_1 + P1) >> 1) - p0_plus_p4,
             p3 = (p0 - 14 * p4 + P2 - ((p2 + P1) << 1)) / 6,
             p1 = P1 - (p0_plus_p4 + p2 + p3);

    addInto(p0_begin, p0.value.cbegin(), p0.value.cend());
    addInto(p1_begin, p1.value.cbegin(), p1.value.cend());
    addInto(p2_begin, p2.value.cbegin(), p2.value.cend());
    addInto(p3_begin, p3.value.cbegin(), p3.value.cend());
    addInto(p4_begin, p4.value.cbegin(), p4.value.cend());

#ifdef MeasureTime
    auto milestone5 = steady_clock::now();
    std::cerr
        << "Memory allocation + copy  : " << duration_cast<milliseconds>(milestone2 - milestone1).count() << " ms\n"
        << "Variable calculation (+/-): " << duration_cast<milliseconds>(milestone3 - milestone2).count() << " ms\n"
        << "Recursion                 : " << duration_cast<milliseconds>(milestone4 - milestone3).count() << " ms\n"
        << "Equation resolution       : " << duration_cast<milliseconds>(milestone5 - milestone4).count() << " ms\n";
#undef MeasureTime
#endif
}

Main function for test

If you want to test for yourselves, you'll need to download the full source available on my repository and e.g. use a main() function such as the one below.

It has 2 test cases:

  • Squaring of a number.
    Operands have the same size and there is definitely a performance gain.
  • Factorial calculation.
    The factorial algorithm is designed to multiply big operands together (it does not do ((((1 * 2) * 3) * 4) ...) but it spends a non-negligible amount of time doing multiplications with at least one operand under the Toom-3 threshold.
    Although it is close to impossible to tell how much, we expect - and since the last edit, there is - some performance gain.
#include <bigint.h>
#include <chrono>
#include <iostream>

int main(int argc, char* argv[])
{
    using std::chrono::steady_clock;
    using std::chrono::duration_cast;
    using std::chrono::duration;
    using std::chrono::milliseconds;

    {
        auto p = bigint::power(453256, 25000);
        std::cerr << "Starting multiplication...\n";
        auto start = steady_clock::now();
        //s and q will use 1 round of Toom-3 multiplication each.
        auto s = p * p;
        auto q = s * s;
        auto end = steady_clock::now();

        std::cerr << "Total multiplication time: " << duration_cast<milliseconds>(end - start).count() << " ms\n";

        std::cout << p.toString() << '\n';
        //return 0;
    }

    std::cerr << "\n\n\n";

    {
        std::cerr << "Starting factorial calculation...\n";
        auto start = steady_clock::now();
        // f uses Toom-3 for its calculation.
        auto f = bigint::factorial(200000);
        auto end = steady_clock::now();

        std::cerr << "Total multiplication time: " << duration_cast<milliseconds>(end - start).count() << " ms\n";

        std::cout << f.toString() << '\n';
        return 0;
    }
}

Measures

Edit: With comments provided to my question, I have managed to create some performance gains since the original post. Below are the measurements made for the multiplication and the factorial function now.
This is most likely not fully optimized yet.

For the multiplication case (453,25625,000 squared, then squared again):

  • When Toom-3 is commented out (i.e. mult_karatsuba is called in every scenario where mult_toom3 should have been called):
    Starting multiplication...
    Total multiplication time: 141 ms
    
  • When Toom-3 is not commented:
    Total multiplication time: 88 ms
    

For the factorial (200,000!)

  • When Toom-3 is commented:
    Starting factorial calculation...
    Total multiplication time: 3811 ms
    
  • When Toom-3 is not commented:
    Starting factorial calculation...
    Total multiplication time: 3021 ms
    
New contributor
Atmo is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
\$\endgroup\$
3
  • 1
    \$\begingroup\$ Wouldn't it start calling mult_karatsuba() after just a few levels of recursion? So it wouldn't be surprising to see the speed of your implementation of Toom3 be equal to that of Karatsuba for values that are just a bit larger than 10000 limbs. \$\endgroup\$
    – G. Sliepen
    Commented Jun 8 at 8:10
  • \$\begingroup\$ @G.Sliepen Yes, mult_karatsuba is very quickly invoked but the operands of the multiplication test case were chosen to demonstrate that at these scales, the performance can visibly be improved by even 1 call to mult_toom3. I added some more details about a few other tweaks that could turn 1M! into hundreds of calls to mult_toom3, yet with no visible gain of performance (as opposed to the first test case). \$\endgroup\$
    – Atmo
    Commented Jun 8 at 10:12
  • \$\begingroup\$ The time between milestones 4 and 3 includes vanilla and Karatsuba as well as further Toom3 processing \$\endgroup\$
    – greybeard
    Commented Jun 8 at 16:34

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.