README

   1 = Continued Fraction Library =
   2 Ben Lynn, September 2008
   3
   4 Frac is a C library for computing with continued fractions.
   5
   6 == Continued Fractions ==
   7
   8 1/3 has no finite floating point representation.  How can a standard tool for
   9 working with real numbers be unable to handle a number a young child could
  10 understand?
  11
  12 Enter continued fractions. Finite continued fractions uniquely (if you forbid
  13 the last term from being 1) and exactly represent rationals.  Contrast this
  14 with 1/3 = 0.333... and 1 = 0.999...  They live in a realm free from the
  15 tyranny of base 2. Commonly encountered irrational numbers, both algebraic and
  16 transcendental, are readily represented with infinite continued fractions, and
  17 when we truncate them we have an approximation along with an error bound. In
  18 contrast, floating point numbers must be accompanied at all times with their
  19 error bounds.
  20
  21 Repeated digits in decimal representation occur for fractions, but why waste
  22 such an interesting phenomenon on something as mundane as a rational? With
  23 continued fractions, repeated terms occur only for square roots. Other infinite
  24 but patterned sequences represent other celebrated numbers, such as e and pi.
  25
  26 Let's say we use floating-point to compute some answer to 100 decimal places.
  27 What if we now want 200 decimal places? Even if we possess the first 100
  28 digits, and carefully preserved the program state, we must start from scratch
  29 and redo all our arithmetic with 200 digit precision.
  30
  31 If only we had used continued fractions, for then we could arbitrarily
  32 increase the precision of our answer without redoing any work. Furthermore,
  33 the error analysis becomes trivial.
  34
  35 We must pay a price for the advantages of continued fractions.  We
  36 have a representation from which we can instantly infer error bounds. A
  37 representation for which we may arbitrarily increase precision without redoing
  38 any work. A representation which in some sense is the closest we can hope to
  39 get to an exact number. Is it any wonder therefore that binary operations,
  40 whose output must conform to this wondrous representation, are clumsy and
  41 convoluted, and that deriving continued fraction expansions for many
  42 reals is time-consuming?
  43
  44 These drawbacks are almost certainly why these fascinating creatures remain
  45 obscure.  Yet surely there are some problems that scream for continued
  46 fractions?  How about a screen saver that computes more and more digits of pi,
  47 each time picking up where it left off?  Or imagine a real-time application
  48 where all continued fractions simply output the last computed convergent on a
  49 timer interrupt.  Under heavy system load, approximations are coarser but the
  50 show goes on.
  51
  52 == Installation ==
  53
  54 I use the GMP library for multi-precision arithmetic. Type
  55
  56   $ make
  57
  58 to compile. Try "pi 2000" to generate 2000 digits of pi using a continued
  59 fraction. It takes longer than I had hoped, but it's faster than
  60
  61   $ echo "scale = 2000; 4*a(1)" | bc -l
  62
  63 There's no API documentation; see the source.
  64
  65 == Demand channels ==
  66
  67 Threads are the future, if not already the present. Once, single core CPUs grew
  68 faster at an astounding pace. But limits are being approached, and engineers
  69 are turning to multicore designs to surmount technical obstacles. Multicore
  70 systems are already commonplace, and as time passes, the number of cores per
  71 system will steadily march upward. Happily, this suits continued fractions.
  72
  73 Inspired by Doug McIlroy's "Squinting at Power Series" (search for it), we
  74 spawn at least one thread per continued fraction, and also per operation on
  75 continued fractions. The threads communicate via crude demand channels, each
  76 providing a few output terms upon each request.
  77
  78 This scheme allows us to quickly write code for generating continued
  79 fraction expansions or operating on them. For example, consider the code for
  80 generating 'e' = 2.71828...
  81
  82 ..............................................................................
  83 // e = [2; 1, 2, 1, 1, 4, 1, ...]
  84 static void *e_expansion(cf_t cf) {
  85   int even = 2;
  86   cf_put_int(cf, even);
  87
  88   while(cf_wait(cf)) {
  89     cf_put_int(cf, 1);
  90     cf_put_int(cf, even);
  91     even += 2;
  92     cf_put_int(cf, 1);
  93   }
  94
  95   return NULL;
  96 }
  97 ..............................................................................
  98
  99 The function +cf_put+ places a term on the output channel, and +cf_wait+ is our
 100 implementation of demand channels, a sort of cooperative multitasking.  Without
 101 it, not only would we have to destroy and clean up after the thread ourselves,
 102 but more seriously, the thread might consume vast amounts of resources on
 103 unwanted terms.  The +cf_wait+ functions instructs the thread to stay idle.
 104 Our threads call this function often, and if it returns zero, our threads
 105 clean themselves up and exit.
 106
 107 == Bugs and other issues ==
 108
 109 I intend to devote time to projects that operate on real numbers rather than
 110 study real numbers for its own sake. But perhaps then I'll find a perfect
 111 application for continued fractions, which will push me to fix deficiencies
 112 in this code.
 113
 114 - The API could use cosmetic surgery.
 115
 116 - I want to replace the condition variable with a semaphore to silence Helgrind
 117   (whose documentation states that condition variables almost unconditionally
 118   and invariably raise false alarms).
 119
 120 - The Makefile is fragile and annoying to maintain.
 121
 122 - The testing infrastructure needs upgrading.
 123
 124 - I'd like a split function as described by McIlroy, though for continued
 125   fractions a "tee"-like function is more apt: rather than spawning a thread
 126   for each term held back, we spawn and watch over two threads on invocation to
 127   service the two demand channels, storing all backlogged terms in one place.
 128
 129 === Increasing floating-point precision ===
 130
 131 A disclaimer to some of my anti-floating-point rhetoric: increasing the
 132 precision of floating-point operations without redoing work is in fact possible
 133 for many common cases. For example, Newton's method is self-correcting, meaning
 134 we may use low precision at first, and increase it for future iterations.
 135 Even pi enjoys such methods: there is a formula computing any hexadecimal
 136 digit of pi without computing any previous digits, though it requires about the
 137 same amount of work.
 138
 139 However, I'd prefer continued fractions to Newton's method when writing code
 140 to compute digits of, say, sqrt(2). Precision decisions are made automatically:
 141 one simply needs enough bits to hold the last computed convergents.