/*
  Author:  Pate Williams (c) 1997

  "Algorithm 1.3.13 (Lehmer). Given a real number
  x by two rational numbers a / b and a' / b'
  such that a / b <= x <= a' / b', this algorithm
  computes the continued fraction expansion of x
  and stops exactly when it is not possible to
  determine the next partial quotient from the
  given approximates a / b and a' / b', and it
  gives lower and upper bounds for this next
  partial quotient". - Henri Cohen - See "A Course
  in Computational Algebraic Number Theory" by
  Henri Cohen page 22.
*/

#include <limits.h>
#include <stdio.h>

void Lehmer(long a, long b, long ap, long bp)
{
  long ai, q, qp, r, rp;

  L2:
  q = a / b;
  r = a - q * b;
  rp = ap - bp * q;
  if (rp < 0 || rp >= bp) {
    qp = ap / bp;
    goto L4;
  }
  ai = q;
  printf("%ld\n", ai);
  a = b;
  b = r;
  ap = bp;
  bp = rp;
  if (b != 0 && bp != 0) goto L2;
  if (b == 0 && bp == 0) goto L4;
  if (b == 0) {
    q = LONG_MAX;
    qp = ap / bp;
  }
  if (bp == 0) {
    q = a / b;
    qp = LONG_MAX;
  }
  L4:
  if (q > qp)
    printf("%ld <= a[i] <= %ld\n", qp, q);
  else
    printf("%ld <= a[i] <= %ld\n", q, qp);
}

int main(void)
{
  long a = 355, b = 113, ap = 314159265l, bp = 100000000l;

  /* approximate pi */
  Lehmer(a, b, ap, bp);
  return 0;
}