Continued Fractions

Exercise: Write a C program to compute 37 digits of the square root of 2, without arbitrary precision arithmetic nor multiplication.

Sound tough? Not if you know continued fractions. One solution lies below.

Some of these pages summarize Bill Gosper’s famous unpublished paper in the field, though I work left to right instead of right to left. Most of the remaining pages were cobbled together from ancient scribbles in one of my notebooks, where I neglected to record the sources.

Solution to exercise: I’m assuming long long gives at least 8-byte integers on the platform.

Multiplication by a constant is acceptable, but for extra credit my program instead shifts and adds.

The inner loop iterates at most 9 times, so contribute only a constant factor to the time complexity. I have not analyzed my program fully, but I believe if it were rewritten to use arbitrary precision addition, then computing n digits takes O(n^2) time, asymptotically slower than using Newton’s method and a clever multiplication routine. This is inevitable: continued fractions converge linearly for square roots while Newton’s method converges quadratically.

// Compute 37 digits of sqrt(2).

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

int main(int argc, char **argv) {
  // 4 byte precision gives 18 digits,
  // 8 byte precision gives 37 digits.
  int n = 37;
  if (argc > 1) {
    n = atoi(argv[1]);
    if (n <= 0) n = 37;
  }
  unsigned long long p0, p1, q0, q1, r0, r1;

  // sqrt(2) = [1; 2, 2, ...]
  // Handle first term here.
  p0 = p1 = q1 = 1;
  q0 = 0;

  while (n > 0) {
    // Continued fraction recurrence relations.
    r0 = p0 + (p1 << 1);
    p0 = p1;
    p1 = r0;

    r0 = q0 + (q1 << 1);
    q0 = q1;
    q1 = r0;

    // digit <- floor(p0 / q0)
    // r0 = q0 * (digit + 1)
    // r1 = q1 * (digit + 1)
    int digit = 0;
    r0 = q0;
    r1 = q1;
    while (digit < 9 && r0 <= p0) {
      r0 += q0;
      r1 += q1;
      digit++;
    }

    if (r1 > p1) {
      r1 -= q1;
      if (r1 <= p1) {
        n--;
        putchar(digit + '0');
        p1 -= r1;
        r0 -= q0;
        p0 -= r0;

        // p0 <- 10 * p0, p1 <- 10 * p1
        p0 = ((p0 << 2) + p0) << 1;
        p1 = ((p1 << 2) + p1) << 1;
      }
    }
  }
  putchar('\n');
  return 0;
}

Ben Lynn blynn@cs.stanford.edu 💡

Continued Fractions

Contents