Compute 100,000 digits of sqrt(2) in pure ANSI C

.gitignore

sqrt2
*.txt

Makefile

CFLAGS = -W -Wall -Wextra -O3 -g -ansi

sqrt2 : sqrt2.c

run : sqrt2
	./$^ | bc

sqrt2.c

#include <stdio.h>
#include <stdlib.h>
#include <inttypes.h>

typedef struct {
    size_t size;
    uint16_t *v;
} bigint;

bigint *bigint_alloc(size_t size, uint16_t init) {
    bigint *n = malloc(sizeof(bigint));
    n->size = size;
    n->v = calloc(size, sizeof(uint16_t));
    n->v[0] = init;
    return n;
}

void bigint_free(bigint *n)
{
    free(n->v);
    free(n);
}

void bigint_trim(bigint *n)
{
    while (n->v[n->size - 1] == 0 && n->size > 0) n->size--;
    n->v = realloc(n->v, n->size * sizeof(uint16_t));
}

void bigint_print(FILE *out, bigint *n)
{
    int i;
    for (i = n->size - 1; i >= 0; i--) {
        fprintf(out, (unsigned) i == n->size - 1 ? "%X" : "%04X", n->v[i]);
    }
    if (n->size == 0) puts("0");
}

bigint *bigint_add(bigint *a, bigint *b)
{
    size_t i, max = a->size > b->size ? a->size : b->size;
    bigint *n = bigint_alloc(max + 1, 0);
    for (i = 0; i < max; i++) {
        if (i > a->size)
            n->v[i] = b->v[i];
        else if (i > b->size)
            n->v[i] = a->v[i];
        else {
            uint32_t m = a->v[i] + b->v[i];
            n->v[i] += m & 0xffff;
            n->v[i + 1] = m >> 16;
        }
    }
    bigint_trim(n);
    return n;
}

/* Add a short to the bigint at index p, taking care of overflow. */
void bigint_adds(bigint *n, size_t p, uint16_t s)
{
    uint32_t m = s;
    while (m > 0) {
        m += n->v[p];
        n->v[p++] = m & 0xffff;
        m >>= 16;
    };
}

bigint *bigint_multiply(bigint *a, bigint *b)
{
    bigint *n = bigint_alloc(a->size + b->size + 1, 0);
    size_t ai, bi;
    for (bi = 0; bi < b->size; bi++) {
        for (ai = 0; ai < a->size; ai++) {
            uint32_t m = a->v[ai] * b->v[bi];
            bigint_adds(n, bi + ai, m & 0xffff);
            bigint_adds(n, bi + ai + 1, m >> 16);
        }
    }
    bigint_trim(n);
    return n;
}

bigint *two, *n, *d;

/* Bring n/d closer to sqrt(2): a(n+1) = a(n)/2 + 1/a(n) */
void iterate()
{
    bigint *nd = bigint_multiply(n, d);
    bigint *nd2 = bigint_multiply(nd, two);
    bigint *nn = bigint_multiply(n, n);
    bigint *dd = bigint_multiply(d, d);
    bigint *dd2 = bigint_multiply(dd, two);
    bigint *sum = bigint_add(nn, dd2);
    bigint_free(n);
    bigint_free(d);
    n = sum;
    d = nd2;
    bigint_free(nd);
    bigint_free(nn);
    bigint_free(dd);
    bigint_free(dd2);
}

int main()
{
    /* Initialize. */
    two = bigint_alloc(1, 2);
    n = bigint_alloc(1, 1);
    d = bigint_alloc(1, 1);

    /* Iterate to cover 100,000 digits. */
    int i;
    for (i = 0; i < 19; i++) iterate();

    printf("scale=100000\nibase=16\n");
    bigint_print(stdout, n);
    printf("/");
    bigint_print(stdout, d);
    printf("\n");
    return 0;
}