exercise.c

   1 // Compute 37 digits of sqrt(2).
   2
   3 #include <stdio.h>
   4 #include <stdlib.h>
   5
   6 int main(int argc, char **argv) {
   7   // 4 byte precision gives 18 digits,
   8   // 8 byte precision gives 37 digits.
   9   int n = 37;
  10   if (argc > 1) {
  11     n = atoi(argv[1]);
  12     if (n <= 0) n = 37;
  13   }
  14   unsigned long long p0, p1, q0, q1, r0, r1;
  15
  16   // sqrt(2) = [1; 2, 2, ...]
  17   // Handle first term here.
  18   p0 = p1 = q1 = 1;
  19   q0 = 0;
  20
  21   while (n > 0) {
  22     // Continued fraction recurrence relations.
  23     r0 = p0 + (p1 << 1);
  24     p0 = p1;
  25     p1 = r0;
  26
  27     r0 = q0 + (q1 << 1);
  28     q0 = q1;
  29     q1 = r0;
  30
  31     // digit <- floor(p0 / q0)
  32     // r0 = q0 * (digit + 1)
  33     // r1 = q1 * (digit + 1)
  34     int digit = 0;
  35     r0 = q0;
  36     r1 = q1;
  37     while (digit < 9 && r0 <= p0) {
  38       r0 += q0;
  39       r1 += q1;
  40       digit++;
  41     }
  42
  43     if (r1 > p1) {
  44       r1 -= q1;
  45       if (r1 <= p1) {
  46         n--;
  47         putchar(digit + '0');
  48         p1 -= r1;
  49         r0 -= q0;
  50         p0 -= r0;
  51
  52         // p0 <- 10 * p0, p1 <- 10 * p1
  53         p0 = ((p0 << 2) + p0) << 1;
  54         p1 = ((p1 << 2) + p1) << 1;
  55       }
  56     }
  57   }
  58   putchar('\n');
  59   return 0;
  60 }