Nonsimple Continued Fractions

If a continued fraction has numerators unequal to 1, or nonpositive denominators apart from an initial zero denomimator, then we say it is nonsimple or nonregular.

Some reals admit no discernible pattern in their simple continued fraction expansion, yet possess elegant nonsimple continued fraction expansions. For example, \(\pi = [3; 7, 15, 1, 292, 1, 1, ?!?]\), but

\[ \frac{4}{\pi} = 1 + \frac{1}{3 + \frac{4}{5 + \frac{9}{7 + ...}}} \]

Also

\[ \pi = 3 + \frac{1}{6 + \frac{9}{6 + \frac{25}{6 + ...}}} \]

Let

\[ \alpha = b_1 + \frac{a_2}{b_2 + \frac{a_3}{b_3 + ...}} \]

be a nonsimple continued fraction. We sometimes write the above more compactly as

\[ \alpha = b_1 + \frac{a_2}{b_2 +} \frac{a_3}{b_3 + } ... \]

By induction, we can show the convergents \(p_i\) are related by the recurrence \(p_{i+1} = b_{i+1} p_i + a_{i+1} p_{i-1}\) (where \(a_1 = 1\)) and similarly for \(q_i\). Also \(p_{i+1} q_i - p_i q_{i+1} = (-1)^i a_1...a_i\). Unlike simple continued fractions we have no guarantees \(p_i, q_i\) are coprime, nor that the sequence converges.

This suggests the following procedure for computing convergents for a sufficiently well-behaved nonsimple continued fraction. We start a table as before, but with an additional row for the numerator written above the denominators. For example, for \(4/\pi\) we write:

1

1

4

9

16

25

1

3

5

7

9

11

0

1

1

0

Again we work on the rows from left to right. For each entry, we compute (number to left) \(\times\) (lower heading) + (number 2 to the left) \(\times\) (upper heading):

1

1

4

9

16

25

1

3

5

7

9

11

0

1

1

4

24

204

2220

1

0

1

3

19

160

1744

We divide all the numbers involved in the last two convergents by 4 before proceeding.

1

1

4

9

16

25

1

3

5

7

9

11

0

1

1

4

24

51

555

7380

1

0

1

3

19

40

436

5796

Here’s a C program to compute the first few convergents of \(\pi\) from the above expansion. We initialize the table with \( \left(\begin{smallmatrix} 4&0\\0&1 \end{smallmatrix} \right) \) so we are computing the 1/4 of the reciprocal of the expansion.

#include <stdio.h>

int main() {
  int num, denom;     // Numerator, denominator of continued fraction.
  int pold, p, pnew;  // Numerator of convergents.
  int qold, q, qnew;  // Denominator of convergents.

  // pi = 1/4 of the reciprocal of the expansion.
  pold = 4; p = 0;
  qold = 0; q = 1;

  // Compute the first cases separately, before the pattern kicks in.
  num = 1; denom = 1;
  pnew = p * denom + pold * num;  // The recurrence relations.
  qnew = q * denom + qold * num;

  printf("p/q = %d/%d\n", pnew, qnew);
  pold = p; qold = q;
  p = pnew; q = qnew;

  num = 1; denom = 3;
  pnew = p * denom + pold * num;
  qnew = q * denom + qold * num;

  printf("p/q = %d/%d\n", pnew, qnew);
  pold = p; qold = q;
  p = pnew; q = qnew;

  // Iterate a few times only, because we'll soon outgrow sizeof(int).
  int i;
  for (i = 1; i < 5; i++) {
    denom += 2;        // The next odd number.
    num += 2 * i + 1;  // The next square
    pnew = p * denom + pold * num;
    qnew = q * denom + qold * num;
    // We should divide by the GCD. See next program.
    printf("p/q = %d/%d\n", pnew, qnew);
    pold = p; qold = q;
    p = pnew; q = qnew;
  }

  return 0;
}

Assuming it behaves well, we can convert a generalized continued fraction to a simple continued fraction on the fly. We modify the above code to use GMP, a multiprecision library, and to output a continued fraction rather than convergents:

#include <stdio.h>
#include <gmp.h>

int main() {
  mpz_t num, denom;

  mpz_t pold, p, pnew;
  mpz_t qold, q, qnew;
  mpz_t t0, t1;  // t is for "temporary variable".

  void recur() {
    // Compute the recurrence relations.
    mpz_mul(pnew, p, denom);
    mpz_addmul(pnew, pold, num);
    mpz_mul(qnew, q, denom);
    mpz_addmul(qnew, qold, num);

    // Remove any nontrivial GCD.
    mpz_gcd(t0, pnew, qnew);
    mpz_gcd(t1, p, q);
    mpz_gcd(t0, t0, t1);
    if (mpz_cmp_ui(t0, 1)) {
      mpz_divexact(pnew, pnew, t0);
      mpz_divexact(qnew, qnew, t0);
      mpz_divexact(p, p, t0);
      mpz_divexact(q, q, t0);
    }

    // Compare integer parts.
    // Use pold, qold as temporary variables.
    // They will be overwritten soon anyway.
    mpz_fdiv_qr(pold, t0, p, q);
    // We keep it simple by performing a second division and remainder.
    // Faster: multiply and add the quotient just found to see if it
    // equals the integer part of pnew/qnew.
    mpz_fdiv_qr(qold, t1, pnew, qnew);
    if (!mpz_cmp(pold, qold)) {
      // If equal then output and subtract the quotient,
      // then invert before the next iteration.
      gmp_printf(" %Zd", pold);
      mpz_set(pold, q);
      mpz_set(qold, t0);
      mpz_set(p, qnew);
      mpz_set(q, t1);
    } else {
      mpz_set(pold, p); mpz_set(p, pnew);
      mpz_set(qold, q); mpz_set(q, qnew);
    }
  }

  mpz_init(num); mpz_init(denom);
  mpz_init(pold); mpz_init(p); mpz_init(pnew);
  mpz_init(qold); mpz_init(q); mpz_init(qnew);
  mpz_init(t0); mpz_init(t1);

  mpz_set_ui(pold, 4); mpz_set_ui(p, 0);
  mpz_set_ui(qold, 0); mpz_set_ui(q, 1);

  mpz_set_ui(num, 1); mpz_set_ui(denom, 1);
  recur();

  mpz_set_ui(num, 1); mpz_set_ui(denom, 3);
  recur();

  int i;
  for (i = 1; i < 100; i++) {
    mpz_add_ui(denom, denom, 2);
    mpz_add_ui(num, num, (i << 1) + 1);
    recur();
  }
  printf("\n");

  return 0;
}

The code requires the nested function feature of gcc (see Breuel for motivation and implementation). The output:

3 7 15 1 292 1 1 1 2 1 3 1 14 2 1 1 2 2 2 2 1 84 2 1 1 15 3 13 1 4 2 6 6 99 1 2 2 6 3 5 1 1 6 8 1 7 1 2 3 7 1 2 1 1 12 1 1 1 3 1 1 8 1 1 2 1 6 1 1 5 2 2 3 1 2

Negative Convergents

Negative terms can lead to negative convergents, and we must tread carefully. Consider the graph of \(y = \frac{p x - p'}{q x - q'}\) for nonnegative \(x\). As before, we are looking at part of a hyperbola, and again the \(y\)-intercept is \(p'/q'\), but this time \(y\) approaches \(\pm \infty\) at \(x = q'/q\), the sign depending on which convergent is larger, and the direction of approach. As \(x \rightarrow \infty\), again we have \(y \rightarrow p/q\). Thus the answer lies outside \(p,q\) and \(p', q'\). We therefore avoid computing terms from pairs of convergents whose denominators have opposite signs.

Pathological Cases

Generalized continued fractions can be difficult to work with. Consider a simple continued fraction \([a_1; a_2, ...]\) where we lift the restriction \(a_2, ... \gt 0\). Then \(m, 0, n\) may be replaced by the single term \(m + n\). Thus a sequence such as \(1, 2, 3, 0, -3, -2, -1\) may be replaced by the single term \(0\). To further confuse the situation, we may replace \(-a, -b, -c, -d\) with \(-a-1,1,b-1,c,d\).

Fortunately, we can ignore most of these issues. Gosper states that among the generalized continued fractions that arise in practice, the main troublesome cases are those involving the occasional 0, and those of the form \(1, -1, 1, -1, 1, -1\). The latter may be dropped entirely.


Ben Lynn blynn@cs.stanford.edu 💡