Skip to content

@skeeto /.gitignore
Created

Embed URL

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
Compute 100,000 digits of sqrt(2) in pure ANSI C
sqrt2
*.txt
CFLAGS = -W -Wall -Wextra -O3 -g -ansi
sqrt2 : sqrt2.c
run : sqrt2
./$^ | bc
#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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.