Arbitrary precision arithmetic and pi

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); 
}