A _pythonic_ pidigts calculator for The Computer Lanaguage Benchmarks Game, using builtin arbitrary precision integers

pidgits.py

# A _pythonic_ pidigts calculator for The Computer Lanaguage Benchmarks Game
# Unlike the other python versions of this program, this uses python's builtin integers.
# I'd argue this is the 'pythonic' way to do this, since it takes advantage of builtin types.
# Derived from the C version and the Chappel version, but with more descriptive names.
from itertools import count, zip_longest, islice, repeat
import sys

class Term(object):
	__slots__ = ("accumulator", "numerator", "denominator") # got2gofast!
	def __init__(self, accumulator, numerator, denominator):
		self.accumulator = accumulator
		self.numerator = numerator
		self.denominator = denominator

	def next_term(self, k):
		k2 = k * 2 + 1
		self.accumulator += self.numerator * 2
		self.accumulator *= k2
		self.denominator *= k2
		self.numerator *= k
		return self

	def extract_digit(self, nth):
		digit = ((self.numerator * nth) + self.accumulator) // self.denominator
		assert digit >= 0 and digit <= 9, \
			"Invalid digit {} extracted at index {} from {!r}".format(digit, nth, self)
		return digit

	def eliminate_digit(self, digit):
		self.accumulator -= self.denominator * digit
		self.accumulator *= 10
		self.numerator *= 10

	def __repr__(self):
		return "Term(accumulator={}, numerator={}, denominator={})"\
			.format(self.accumulator, self.numerator, self.denominator)

def generatePiDigits():
	"""An infinite generator for the digits of PI"""
	term = Term(0, 1, 1)
	term_generator = map(term.next_term, count(1))
	while True:
	    next(term_generator)
	    if term.numerator > term.accumulator:
	        continue
	    digit = term.extract_digit(3)
	    if digit != term.extract_digit(4):
	        continue
	    # At this point 'digit' is our result
	    # Before we yield it, eliminate our digit
	    term.eliminate_digit(digit)
	    yield digit

#Precompute the string-representation of the digits
_DIGIT_STRINGS=tuple(str(x) for x in range(10))
digitToString=getattr(_DIGIT_STRINGS, '__getitem__')

def printDigits(digits, amount, output=sys.stdout):
    digits=map(digitToString, digits)
    remaining_digits=amount%10
    for count in range(10, amount, 10):
        output.write(''.join(islice(digits, 10)))
        output.write('\t: ')
        output.write(str(count))
        output.write('\n')
    if remaining_digits:
        output.write(''.join(islice(digits, remaining_digits)).ljust(10))
        output.write('\t: ')
        output.write(str(amount))
        output.write('\n')

amount = int(sys.argv[1]) if len(sys.argv) > 1 else 27

printDigits(generatePiDigits(), amount)

pidigits.c

#include <stdio.h>
#include <stdint.h>
#include <stdbool.h>
#include <errno.h>
#include <gmp.h>
#include <stdlib.h>
#include <assert.h>

typedef struct Term Term;

Term *create_term(
    uint64_t accumulator,
    uint64_t numerator,
    uint64_t denominator,
    uint64_t k
);

Term *create_initial_term() {
    // Term(accumulator=0, numerator=1, denominator=1, k=0)
    return create_term(0, 1, 1, 0);
}

unsigned char next_pi_digit(Term *term);

void print_pi_digits(const unsigned long long amount) {
    Term *term = create_initial_term();
    unsigned long long remaining = amount;
    // While we have at least ten digits, we don't need to pad anything
    for (; remaining >= 10; remaining -= 10) {
        for (int i = 0; i < 10; i++) {
            putchar('0' + next_pi_digit(term));
        }
        printf("\t:%llu\n", amount - remaining + 10);
    }
    // Output remaining, and pad the ending if we have anything left
    if (remaining > 0) {
        assert(remaining < 10);
        for (int i = 0; i < 10; i++) {
            if (i < remaining) {
                putchar('0' + next_pi_digit(term));
            } else {
                putchar(' ');
            }
        }
        printf("\t:%llu\n", amount);
    }
}

int main(int argc, char **argv) {
    unsigned long long amount = 27; // Default to printing the first 27 digits
    if (argc > 1) {
        // Retarted C way to safely parse numbers (ignore this)
        errno = 0;
        amount = strtoull(argv[1], NULL, 10);
        switch (errno) {
            case 0:
                // Successfully parsed it
                break;
            case ERANGE:
                printf("Can't handle more than 2^64 digits of pi: %s", argv[1]);
                return 1;
            default:
                printf("Invalid number of digits: %s", argv[1]);
                return 1;
        }
        if (amount == 0) {
            printf("Can't print zero digits of PI!");
            return 1;
        }
    }
    print_pi_digits(amount);
    return 0; 
}


//
// Math Implementation, ordered by increasing voodo-level
//

struct Term {
    mpz_t accumulator;
    mpz_t numerator;
    mpz_t denominator;
    uint64_t k;
};

void next_term(Term* term);

unsigned char extract_digit(Term *term, uint64_t digit_index);

void eliminate_digit(Term *term, unsigned char digit);

inline unsigned char next_pi_digit(Term *term) {
    unsigned char digit;
    do {
        next_term(term);
        // if (term.numerator > term.accumulator) continue;
        if (mpz_cmp(term->numerator, term->accumulator) > 0) continue;
        digit = extract_digit(term, 3);
    } while (digit != extract_digit(term, 4));
    eliminate_digit(term, digit);
    return digit;
}

inline void next_term(Term *term) {
    uint64_t k = ++term->k;
    uint64_t k2 = k * 2 + 1;
    /*
     * Effective Pesudocode:
     * long k2 = k * 2 = 1;
     * term.accumulator += term.numerator * 2
     * term.accumulator *= k2
     * term.denominator *= k2
     * term.numerator *= k
     */
    // Optimization TODO: investigate bitshifting term.numerator instead of multiplying
    mpz_addmul_ui(term->accumulator, term->numerator, 2);
    mpz_mul_ui(term->accumulator, term->accumulator, k2);
    mpz_mul_ui(term->denominator, term->denominator, k2);
    mpz_mul_ui(term->numerator, term->numerator, k);
}


inline void eliminate_digit(Term *term, unsigned char digit) {
    /*
     * Effective Pesudocode
     * self.accumulator -= self.denominator * digit
	 * self.accumulator *= 10
	 * self.numerator *= 10
	 */
    assert(digit <= 9);
    mpz_submul_ui(term->accumulator, term->denominator, digit);
    mpz_mul_ui(term->accumulator, term->accumulator, 10);
    mpz_mul_ui(term->numerator, term->numerator, 10);
}

// NOTE: These are assumed to be initialized anytime a calculation is called
mpz_t tmp1, tmp2;
bool initialized = false;

inline unsigned char extract_digit(Term *term, uint64_t digit_index) {
    /*
     * Effective Pesudocode:
     * u8 digit = ((term.numerator * digit_index) + term.accumulator) / term.denominator
     */
    // Optimization: joggling between tmp1 and tmp2, so GMP won't have to use temp buffers
    // We assume these've been initailzied by the 'create_initial_term' function'
    assert(initialized);
    assert(digit_index < 10);
    mpz_mul_ui(tmp1, term->numerator, digit_index);
    mpz_add(tmp2, tmp1, term->accumulator);
    mpz_tdiv_q(tmp1, tmp2, term->denominator);
    assert(mpz_get_ui(tmp1) <= 9);
    return (unsigned char) mpz_get_ui(tmp1);
}

// Initialization code (boring)

void initialize() {
    assert(!initialized);
    mpz_init(tmp1);
    mpz_init(tmp2);
    initialized = true;
}

inline Term *create_term(
    uint64_t accumulator,
    uint64_t numerator,
    uint64_t denominator,
    uint64_t k
) {
    Term *term = calloc(1, sizeof(Term));
    mpz_init_set_ui(term->accumulator, accumulator);
    mpz_init_set_ui(term->numerator, numerator);
    mpz_init_set_ui(term->denominator, denominator);
    term->k = k;
    // Lazy-initialie our temporary variables
    if (!initialized) initialize();
    return term;
}