README

   1 = Continued Fraction Library =
   2
   3 Ben Lynn, September 2008
   4
   5 == Introduction ==
   6
   7 1/3 has no finite floating point representation.  How can a standard tool for
   8 working with real numbers be unable to handle a number a young child could
   9 understand?
  10
  11 Enter continued fractions. Finite continued fractions uniquely (if you forbid
  12 the last term from being 1) and exactly represent rationals.  Contrast this
  13 with 1/3 = 0.333... and 1 = 0.999...  They live in a realm free from the
  14 tyranny of base 2. Commonly encountered irrational numbers, both algebraic and
  15 transcendental, are readily represented with infinite continued fractions, and
  16 when we truncate them we have an approximation along with an error bound. In
  17 contrast, floating point numbers must be accompanied at all times with their
  18 error bounds.
  19
  20 Repeated digits in decimal representation occur for fractions, but why waste
  21 such an interesting phenomenon on something as mundane as a rational? With
  22 continued fractions, repeated terms occur only for square roots. Other infinite
  23 but patterned sequences represent other celebrated numbers, such as e and pi.
  24
  25 As with floating point numbers, we truncate infinite continued fractions to compute on them. But unlike floating point numbers, when more precision is desired
  26 we can compute more terms without redoing any work. Imagine a real-time
  27 application where on a timer interrupt, all continued fractions simply output
  28 the last computed convergent. Under heavy system load, approximations are
  29 coarser but the show goes on.
  30
  31 == Installation ==
  32
  33 I use the GMP library for multi-precision arithmetic. Type
  34
  35  $ make
  36
  37 to compile. The "pi" binary generates 5000 digits of pi using a continued
  38 fraction. It takes longer than I had hoped.
  39
  40 There's no API documentation yet; see the source.
  41
  42 == Demand channels ==
  43
  44 Threads are the future, if not already the present. Once, single core CPUs grew
  45 faster at an astounding pace. But limits are being approached, and engineers
  46 are turning to multicore designs to surmount technical obstacles. Multicore
  47 systems are already commonplace, and as time passes, the number of cores per
  48 system will steadily march upward. Happily, this suits continued fractions.
  49
  50 Inspired by Doug McIlroy's "Squinting at Power Series" (search for it), we
  51 spawn at least one thread per continued fraction, and also per operation on
  52 continued fractions. The threads communicate via crude demand channels, each
  53 providing a few output terms upon each request.
  54
  55 This scheme allows us to quickly write code for generating continued
  56 fraction expansions or operating on them. For example, consider the code for
  57 generating 'e' = 2.71828...
  58
  59 ..............................................................................
  60 // e = [2; 1, 2, 1, 1, 4, 1, ...]
  61 static void *e_expansion(cf_t cf) {
  62   mpz_t even, one;
  63   mpz_init(even); mpz_init(one);
  64   mpz_set_ui(even, 2); mpz_set_ui(one, 1);
  65
  66   cf_put(cf, even);
  67
  68   while(cf_wait(cf)) {
  69     cf_put(cf, one);
  70     cf_put(cf, even);
  71     mpz_add_ui(even, even, 2);
  72     cf_put(cf, one);
  73   }
  74
  75   mpz_clear(one);
  76   mpz_clear(even);
  77   return NULL;
  78 }
  79 ..............................................................................
  80
  81 Most of the verbosity is due to GMP. I plan to support plain ints as well by
  82 automatically upgrading them, so the above would be:
  83
  84 ..............................................................................
  85 // e = [2; 1, 2, 1, 1, 4, 1, ...]
  86 static void *e_expansion(cf_t cf) {
  87   int even = 2;
  88   cf_put_int(cf, even);
  89
  90   while(cf_wait(cf)) {
  91     cf_put_int(cf, 1);
  92     cf_put_int(cf, even);
  93     even += 2;
  94     cf_put_int(cf, 1);
  95   }
  96
  97   return NULL;
  98 }
  99 ..............................................................................
 100
 101 So we only need to know about two functions when programming for this library:
 102 +cf_put+ places a term on the output channel, and +cf_wait+ is our
 103 implementation of demand channels, a sort of cooperative multitasking.  Without
 104 it, not only would we have to destroy and clean up after the thread ourselves,
 105 but more seriously, the thread might consume vast amounts of resources on
 106 unwanted terms.  The +cf_wait+ functions instructs the thread to stay idle.
 107 Our threads call this function often, and if it it returns zero, our threads
 108 clean themselves up and exit.
 109
 110 == Bugs and other issues ==
 111
 112 I intend to devote time to projects that operate on real numbers rather than
 113 study real numbers for its own sake. But perhaps then I'll find a perfect
 114 application for continued fractions, which will push me to fix deficiencies
 115 in this code.
 116
 117 - The API could use serious cosmetic surgery.
 118
 119 - I want to replace the condition variable with a semaphore to silence Helgrind
 120   (whose documentation states that condition variables almost unconditionally
 121   and invariably raise false alarms).
 122
 123 - The Makefile is fragile and annoying to maintain.
 124
 125 - The testing infrastructure needs upgrading.
 126
 127 - I'd like a split function as described by McIlroy, though for continued
 128   fractions a "tee"-like function is more apt: rather than spawning a thread
 129   for each term held back, we spawn and watch over two threads on invocation to
 130   service the two demand channels, storing all backlogged terms in one thread.