Created

Embed URL

HTTPS clone URL

SSH clone URL

You can clone with HTTPS or SSH.

Download Gist

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

View .gitignore
1 2 3 4 5 6
CFLAGS = -W -Wall -Wextra -O3 -g -ansi
 
sqrt2 : sqrt2.c
 
run : sqrt2
./$^ | bc
View .gitignore
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
#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.