Arbitrary-precision rational arithmetic in C++

Introduction

This is a follow-up to a post from several years ago describing a C++ implementation of arbitrary-precision unsigned integer arithmetic. This weekend I extended this to also support arbitrary-precision signed integers and rational numbers. Although this started as an educational tool, it now feels a bit more complete, and actually usable for the combinatorics and probability applications of the sort that are frequently discussed here.

I tried to stick to the original objectives of relatively simple and hopefully readable code, with stand-alone, header-only implementation, as freely available in the public domain as legally possible.

The code is available here, as well as on GitHub, in three header files:

  • #include "math_Unsigned.h" defines a math::Unsigned type representing the natural numbers with all of the sensible arithmetic, bitwise, and relational operators, essentially everything except bitwise one’s complement… although more on this shortly.
  • #include "math_Integer.h" defines an Integer type with a sign-and-magnitude implementation in terms of Unsigned, with all corresponding operators, including bitwise operators having two’s complement semantics assuming “infinite sign extension.”
  • #include "math_Rational.h" defines a Rational type implemented in terms of Integer numerator and denominator.

This was a fun exercise; there were interesting challenges in developing each of the three classes. As discussed previously, the unsigned type handles the actual arbitrary-precision representation (implemented as a vector<uint32_t> of digits in base 2^{32}), where division is by far the most complex operation to implement efficiently.

The implementation of the signed integer type is relatively straightforward… except for the bitwise operators. Assuming a sign-and-magnitude representation (using an Unsigned under the hood), it is an interesting exercise to work out how to implement bitwise and, or, xor, and not, so that they have two’s complement semantics even for negative operands. In the process, I had to add an “AND NOT” operator to the original underlying unsigned type (there is actually a built-in operator &^ for this in Go).

With this machinery in place, the rational type is the simplest to implement. The only wrinkle here is that a few additional constructors are needed, since user-defined conversions from the more primitive integral types (e.g., Rational from Integer, Integer from int32_t, etc.) are only implicitly applied “one level deep.”

Example application: Are seven riffle shuffles enough?

To test and demonstrate use of these classes, consider riffle shuffling a standard poker deck of 52 playing cards. How many shuffles are sufficient to “fully randomize” the deck? A popular rule of thumb, attributed to Bayer and Diaconis, is that seven riffle shuffles are recommended. (See a longer list of references here, along with some simpler counting arguments that at least six shuffles are certainly necessary.)

This recommendation is based on analysis of the Gilbert-Shannon-Reeds model of a single riffle shuffle, and of the total variation distance between probability distributions Q^m and U, where Q^m is the distribution of arrangements of the deck after m GSR riffle shuffles, and U is the desired uniform distribution where every arrangement is equally likely. We can compute this total variation distance exactly as a function of the number m of shuffles, as demonstrated in the following example code:

#include "math_Rational.h"
#include <iostream>
using namespace math;

Integer factorial(int n)
{
    Integer f = 1;
    for (int k = 1; k <= n; ++k)
    {
        f *= k;
    }
    return f;
}

Integer binomial(int n, int k)
{
    if (0 <= k && k <= n)
    {
        return factorial(n) / factorial(k) / factorial(n - k);
    }
    else
    {
        return 0;
    }
}

Integer power(int base, int exp)
{
    Integer n = 1;
    for (int k = 0; k < exp; ++k)
    {
        n *= base;
    }
    return n;
}

Integer eulerian(int n, int k)
{
    Integer r = 0;
    for (int j = 0; j < k + 2; ++j)
    {
        r += (power(-1, j) * binomial(n + 1, j) * power(k + 1 - j, n));
    }
    return r;
}

Rational total_variation_distance(int cards, int shuffles)
{
    Rational q = 0;
    for (int r = 1; r <= cards; ++r)
    {
        Rational a = Rational(
            binomial((1 << shuffles) + cards - r, cards),
            Integer(1) << (cards * shuffles)) - Rational(1, factorial(cards));
        q += (eulerian(cards, r - 1) * (a < 0 ? -a : a));
    }
    return q / 2;
}

int main()
{
    int cards = 52;
    for (int shuffles = 0; shuffles <= 15; ++shuffles)
    {
        std::cout << shuffles << " " <<
            total_variation_distance(cards, shuffles).to_double() << std::endl;
    }
}

The following figure shows the results. Total variation distance ranges from a maximum of one (between discrete distributions with disjoint support) to a minimum of zero, in this case corresponding to an exactly uniform distribution of arrangements of the deck.

Total variation distance vs. number of GSR riffle shuffles of a standard 52-card deck.

We can see the sharp threshold behavior, where total variation distance transitions from near one to near zero over just a few shuffles, first dropping below 1/2 at seven shuffles.

This entry was posted in Uncategorized. Bookmark the permalink.

3 Responses to Arbitrary-precision rational arithmetic in C++

  1. Pingback: Arbitrary-precision integer arithmetic in C++ | Possibly Wrong

  2. Nice work. You’ve shared enough code that you just need to make a GitHub or GitLab account already. 🙂 Sharing zip files of code on a generic file host is sooo 1990’s.

    What made you pick the Unlicense? (I’m wondering if I had something to do with it, being that it’s my favorite.)

    • Thanks! The file archive is mostly due to inertia, and a side effect of a mathematician trying to masquerade as a software developer :). I’ve had a GitHub account for a while, just haven’t taken the time to move everything, and wanted to keep it all in the same place.

      But I have created a repo on GitHub for this, and updated the post to link to it. The license is an interesting issue– I’ve changed my thinking on this over the years, and at times been admittedly rather cavalier in my attempt at dedicating to the public domain. The Unlicense seems to be the current best such “serious” attempt.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.