### Solution from PH

This is a simple speed up using a larger digital base (it should be nearly 4 times faster than the original solution.) So this is to further the argument that simple human analysis still rules the day. I did not try the aggressive reciprocal substitution method posed by NB, since that requires going to non-portable code.

The original code performed everything in base 10… I just cranked the base to as large as could be done while making sure the intermediate calculations all fit into a 32 bit number, and still being divisible by 10 (so that I didn’t have to add in a base conversion step.) So the code follows the exact same algorithm, but uses base 1000 numbers.

This was a demonstration of how human analysis can beat the hell out of any compiler analysis (there is no compiler in the world, real or imagined that could possibly do the transformation of the code that I did). Its original construction was not general enough, and I had to work through several unexpected bugs before getting it to work correctly.

OTOH, this code example is not necessarily representative of anything other than itself. I recall that the whole thing came down to a one-line bottleneck. While this kind of thing does happen in real world code, usually not to this extreme.

#include <stdlib.h> #include <stdio.h> #define BASE (10000) void ComputePi (int numdigits, int * pi) { int alength = 40 * numdigits / 3; int * a = (int*) malloc (alength * sizeof(int)); int piLength = 0; int nines = 0; int predigit = 0; int i, j, p; for(i = 0; i < alength; ++i) a[i] = 2; for (j = 0; j < numdigits; ++j) { int q = 0; p = 2 * alength – 1; for (i = alength-1; i >= 1; i–) { int x = BASE*a[i] + q*(i+1); a[i] = x % p; q = x / p; p -= 2; } i = a[0]; a[0] = q % BASE; q = (q / BASE) + i; if (q == (BASE-1)) { ++nines; } else { int k; if (q == BASE) { pi[piLength] = predigit + 1; for (k = 1; k <= nines; ++k) pi[piLength + k] = 0; predigit = 0; } else { pi[piLength] = predigit; for (k = 1; k <= nines; ++k) pi[piLength + k] = BASE-1; predigit = q; } piLength += nines + 1; nines = 0; } } free (a); pi[piLength] = predigit; } int main (int argc, char * argv[]) { int numdigits; int * pi; FILE* fp = fopen (argv[2], “w”); int i, c, t; if (argc <= 1) { fprintf (stderr, “usage: pi #DIGITS [FILE]”); return 1; } t = atoi (argv[1]); numdigits = (t / 4) + 2; pi = (int *) malloc ((numdigits + 1)*sizeof (int)); ComputePi (numdigits, pi); if (argc > 2) { fp = fopen (argv[2], “w”); } else { /* Dunno if this is portable, but it’ll work in Windoze/DOS/UNIX */ fp = fdopen (1, “w”); } if (fp == NULL) { fprintf (stderr, “Cannot open %s\n”, argv[2]); return 2; } c = 0; fprintf (fp, “%d”, pi[1]); c++; for (i=2; c + 4 <= t; i++) { fprintf (fp, “%04d”, pi[i]); c+=4; } t -= c; if (t > 0) { char endbuff[8]; sprintf (endbuff, “%04d”, pi[i]); endbuff[t] = ‘\0’; fprintf (fp, “%s”, endbuff); } fputc (‘\n’, fp); fclose (fp); free (pi); return 0; }

Pages: « Prev 1 2 3 4 5 6 7 8 9 Next »

Discuss (16 comments)