This program implements arbitrary precision (aka bignum) arithmetic and implements a subset of the dc language in order to run a program that calculates pi digit by digit using a spigot algorithm. It is in principle capable of calculating an infinite number of digits -- it will however slow down very quickly and eventually run out of memory. (And the bignum code may have subtle bugs I'm not aware of that kick in with too large numbers...) The pi program was written on another occasion so the main challenge was implementing bignum arithmetic. Addition, subtraction and, to some extent, multiplication were relatively straightforward, however finding and implementing a suitable algorithm for division took some time. Multiplication and division use very crude algorithms with lots of potential for optimization.
#include <stdio.h> #include <stdlib.h> #include <string.h> #include <stdint.h> #include <assert.h> char piprog[] = "1sq180sr60st2si[3li*1+d1+*3*suli27*12-lq*5lr*+lt5*/d48+Psy10lqlid2*1-***" "10lulqli5*2-*lr+lylt*-**srsqlult*stli1+silmx]smlmx"; typedef struct mp mp; struct mp { char neg; int s; int ld; int ref; uint32_t *w; }; mp **stack, **stp; int stsize; mp *reg[256]; void mpfree(mp *r) { free(r->w); free(r); } void mpgrow(mp *r, int s) { if(r->s > s) return; r->s = (s + 16) & -16; r->w = (uint32_t*) realloc(r->w, r->s * 4); } int mpcmp(mp *a, mp *b) { int i; if(a->ld != b->ld) return a->ld - b->ld; for(i = a->ld; i >= 0; i--){ if(a->w[i] > b->w[i]) return 1; if(a->w[i] < b->w[i]) return -1; } return 0; } mp * mpint(int n) { mp *m; if(n < 0){ m = mpint(-n); m->neg = 1; return m; } m = (mp*) malloc(sizeof(mp)); m->ref = 0; m->neg = 0; m->s = 16; m->w = (uint32_t*) malloc(64); m->ld = 0; m->w[0] = n; return m; } mp * mpdup(mp *a) { mp *r; r = (mp*) malloc(sizeof(mp)); r->ref = 0; r->neg = a->neg; r->s = a->s; r->w = (uint32_t*) malloc(a->s * 4); r->ld = a->ld; memcpy(r->w, a->w, (a->ld + 1) * 4); return r; } void mpcpy(mp *r, mp *a) { if(r == a) return; mpgrow(r, a->ld); memcpy(r->w, a->w, (a->ld + 1) * 4); r->neg = a->neg; r->ld = a->ld; } void mpsub(mp *r, mp *a, mp *b); void mpadd(mp *r, mp *a, mp *b) { int c, nc, i, l, rld; uint32_t aw, bw; if(a->neg && !b->neg){ a->neg = 0; mpsub(r, b, a); if(r != a) a->neg = 1; return; } if(b->neg && !a->neg){ b->neg = 0; mpsub(r, a, b); if(r != b) b->neg = 1; return; } r->neg = a->neg; c = 0; if(a->ld > b->ld) l = a->ld + 1; else l = b->ld + 1; mpgrow(r, l); rld = 0; for(i = 0; i <= l; i++){ aw = i <= a->ld ? a->w[i] : 0; bw = i <= b->ld ? b->w[i] : 0; if(c) nc = aw >= ~bw; else nc = aw > ~bw; r->w[i] = aw + bw + c; if(r->w[i] != 0) rld = i; c = nc; } r->ld = rld; } void mpsub(mp *r, mp *a, mp *b) { int i, c, nc, l, rld; uint32_t aw, bw; if(a->neg && !b->neg){ b->neg = 1; mpadd(r, a, b); if(r != b) b->neg = 0; return; } if(!a->neg && b->neg){ b->neg = 0; mpadd(r, a, b); if(r != b) b->neg = 1; return; } if(a->neg && b->neg){ a->neg = b->neg = 0; mpsub(r, b, a); if(r != a) a->neg = 1; if(r != b) b->neg = 1; return; } if(mpcmp(a, b) < 0){ mpsub(r, b, a); r->neg = 1; return; } r->neg = 0; c = 0; mpgrow(r, a->ld); l = a->ld; rld = 0; for(i = 0; i <= l; i++){ aw = a->w[i]; bw = i <= b->ld ? b->w[i] : 0; if(c) nc = aw <= bw; else nc = aw < bw; r->w[i] = aw - bw - c; if(r->w[i] != 0) rld = i; c = nc; } r->ld = rld; } void digadd(uint32_t *r, uint32_t a) { int c; c = *r > ~a; *r++ += a; while(c){ c = ~*r == 0; (*r++)++; } } void mpmul(mp *r, mp *a, mp *b) { int i, j, k; uint64_t x; uint32_t aw, bw; mpgrow(r, a->ld + b->ld + 1); r->neg = a->neg ^ b->neg; memset(r->w, 0, r->s * 4); for(i = 0; i <= a->ld; i++){ aw = a->w[i]; for(j = 0; j <= b->ld; j++){ bw = b->w[j]; x = (uint64_t) aw * bw; digadd(r->w+i+j, x); digadd(r->w+i+j+1, x>>32); } } for(i = a->ld + b->ld + 1; i >= 0; i--) if(r->w[i] != 0) break; r->ld = i; } void mpsh(mp *r, mp *a, int n) { int i, k; if(n == 0){ mpcpy(r, a); return; } if(n > 0){ k = n >> 5; mpgrow(r, a->ld + k + 1); n &= 31; r->w[a->ld + k + 1] = 0; for(i = a->ld; i >= 0; i--){ if(n != 0) r->w[i+k+1] |= a->w[i] >> 32 - n; r->w[i+k] = a->w[i] << n; } for(i = 0; i < k; i++) r->w[i] = 0; r->ld = a->ld + k + (r->w[a->ld+k+1] != 0); }else{ n = -n; k = n >> 5; if(k > a->ld){ r->ld = 0; r->w[0] = 0; } n &= 31; for(i = 0; i <= a->ld; i++){ if(i >= k+1 && n != 0) r->w[i-k-1] |= a->w[i] << 32 - n; if(i >= k) r->w[i-k] = a->w[i] >> n; } r->ld = a->ld - k - (r->w[a->ld-k] == 0); } } int mpbitcnt(mp *a) { int r; uint32_t w; static char ld[16] = {0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3}; r = a->ld << 5; w = a->w[a->ld]; if(w >= 0x10000){ w >>= 16; r += 16; } if(w >= 0x100){ w >>= 8; r += 8; } if(w >= 0x10){ w >>= 4; r += 4; } return r + ld[w]; } void mpdiv(mp *q, mp *a, mp *b) { mp *tq, *ta, *tb; int n, an, bn; if(a->neg || b->neg){ an = a->neg; bn = b->neg; a->neg = b->neg = 0; mpdiv(q, a, b); q->neg = an ^ bn; a->neg = an; b->neg = bn; return; } if(b->ld == 0 && b->w[0] == 0){ printf("division by zero\n"); exit(1); } q->ld = 0; q->w[0] = 0; if(a->ld < b->ld) return; ta = mpdup(a); tb = mpdup(b); n = mpbitcnt(a) - mpbitcnt(b); tq = mpint(1); mpsh(tq, tq, n); mpsh(tb, tb, n); while(n-- >= 0){ if(mpcmp(ta, tb) >= 0){ mpsub(ta, ta, tb); mpadd(q, tq, q); } mpsh(tb, tb, -1); mpsh(tq, tq, -1); } mpfree(ta); mpfree(tb); mpfree(tq); } void mpprint(mp *a) { int i; if(a->neg) printf("-"); printf("0x%x", a->w[a->ld]); for(i = a->ld - 1; i >= 0; i--) printf("%.8x", a->w[i]); printf(" "); } void push(mp *m) { mp **n; if(stp >= stack + stsize){ n = (mp**) realloc(stack, (stsize + 16) * sizeof(mp*)); stsize += 16; stp = n + (stp - stack); stack = n; } *stp++ = m; m->ref++; } mp * pop(void) { assert(stp > stack); return *--stp; } void decref(mp *m) { if(m != NULL && --m->ref == 0) mpfree(m); } void eval(char *s) { mp *m, *a, *b; char c, *p; for(; *s != 0; s++) switch(*s){ case '0': case '1': case '2': case '3': case '4': case '5': case '6': case '7': case '8': case '9': push(mpint(strtol(s, &s, 0))); s--; break; case 's': decref(reg[c = *++s]); reg[c] = pop(); break; case 'l': push(reg[*++s]); break; case '[': s++; p = strchr(s, ']'); assert(p != NULL); m = mpint(0); mpgrow(m, p - s); memcpy(m->w, s, p - s); m->w[p - s] = 0; push(m); s = p; break; case 'x': m = pop(); eval((char*)(m->w)); decref(m); break; case '+': a = pop(); b = pop(); m = mpint(0); mpadd(m, b, a); push(m); decref(a); decref(b); break; case '-': a = pop(); b = pop(); m = mpint(0); mpsub(m, b, a); push(m); decref(a); decref(b); break; case '*': a = pop(); b = pop(); m = mpint(0); mpmul(m, a, b); push(m); decref(a); decref(b); break; case '/': a = pop(); b = pop(); m = mpint(0); mpdiv(m, b, a); push(m); decref(a); decref(b); break; case 'd': assert(stp > stack); push(stp[-1]); break; case 'P': m = pop(); putchar(m->w[0]); decref(m); break; default: printf("unknown char %c\n", *s); } } int main() { setbuf(stdout, NULL); eval(piprog); }