README

   1 = frac =
   2
   3 Frac is a C library for computing with http://crypto.stanford.edu/pbc/notes/contfrac/[continued fractions].
   4
   5 == Installation ==
   6
   7 The GMP library is required. Type:
   8
   9   $ make
  10
  11 to compile.
  12
  13 Then try "pi 2000" to generate 2000 digits of pi using a continued
  14 fraction. It takes longer than I had hoped, but it's faster than
  15
  16   $ echo "scale = 2000; 4*a(1)" | bc -l
  17
  18 Also try "hakmem 1000" to compute 1000 digits of the example continued fraction
  19 constant from HAKMEM, which compares favourably with:
  20
  21   $ echo "scale=1000
  22           pi=4*a(1)
  23           x=2*sqrt(5)
  24           (sqrt(3/pi^2+e(1)))/((e(x)-1)/(e(x)+1)-s(69))" | bc -l
  25
  26 Unlike frac, bc can only output the final answer in one go after all
  27 calculations are complete.
  28
  29 There's no API documentation; see the source.
  30
  31 == Continued fractions ==
  32
  33 The binary representation of 1/3 is infinite. Floating-point numbers, our
  34 standard tool for working with reals, cannot even handle numbers a young child
  35 could understand.
  36
  37 In contrast, every rational can be uniquely and exactly represented by a finite
  38 continued fraction (provided we forbid those with last term 1). Continued
  39 fractions are immune from oddities such as 1/3 = 0.333... and 1 = 0.999...,
  40 living in a realm free from the tyranny of binary, or any other base.
  41
  42 Every irrational number can also be uniquely expressed as an infinite continued
  43 fraction. Many frequently encountered algebraic and transcendental numbers have
  44 easily computable continued fraction expansions which can be truncated to yield
  45 approximations with built-in error bounds, unlike floating-point
  46 approximations, which must be accompanied at all times with their error bounds.
  47
  48 Furthermore, suppose we use floating-point to compute some answer to 100
  49 decimal places. What if we now want 200 decimal places? Even if we possess the
  50 first 100 digits and carefully preserve the program state, we must start from
  51 scratch and redo all our arithmetic with 200 digit precision. If we had used
  52 continued fractions, we could arbitrarily increase the precision of our answer
  53 without redoing any work, and without any tedious error analysis.
  54
  55 But we must pay a price for these advantages. We have a representation from
  56 which we can instantly infer error bounds. A representation for which we may
  57 arbitrarily increase precision without repeating any calculations. A
  58 representation which in some sense is the closest we can hope to get to an
  59 exact number using rationals. Is it any wonder therefore that binary operations
  60 on continued fractions, whose output must also uphold these high standards, are
  61 clumsy and convoluted?
  62
  63 These drawbacks are almost certainly why these fascinating creatures remain
  64 obscure. Yet surely there are some problems that scream for continued
  65 fractions. How about a screen saver that computes more and more digits of pi,
  66 each time picking up where it left off?  Or imagine a real-time application
  67 where all continued fractions simply output the last computed convergent on a
  68 timer interrupt. Under heavy system load, approximations are coarser but the
  69 show goes on.
  70
  71 === Disclaimer ===
  72
  73 Lest my salesmanship backfire, let me temper my anti-floating-point rhetoric.
  74 Increasing the precision of floating-point operations without redoing work is
  75 in fact possible for many common cases. For example, Newton's method is
  76 self-correcting, meaning we may use low precision at first, and increase it for
  77 future iterations. Even pi enjoys such methods: there exists a formula
  78 revealing any hexadecimal digit of pi without computing any previous digits,
  79 though it requires about the same amount of work.
  80
  81 Moreover, several floating-point algorithms converge quadratically or faster,
  82 thus good implementations will asymptotically outperform continued fractions
  83 as these often converge linearly.
  84
  85 Nonetheless, for lower accuracies, the smaller overhead may give continued
  86 fractions the edge in certain problems such as finding the square root of 2.
  87 Additionally, precision decisions are automatic, for example, one simply needs
  88 enough bits to hold the last computed convergents. Built-in error analysis
  89 simplifies and hence accelerates development.
  90
  91 == Bugs and other issues ==
  92
  93 I intend to devote time to projects that operate on real numbers rather than
  94 study real numbers for its own sake. But perhaps then I'll find a perfect
  95 application for continued fractions, which will push me to fix deficiencies in
  96 this code:
  97
  98 - I want to replace the condition variable with a semaphore to silence Helgrind
  99   (whose documentation states that condition variables almost unconditionally
 100   and invariably raise false alarms).
 101
 102 - Computing the continued fraction expansion of pi from the Chudnovsky
 103   brothers' Ramanujan formula would be much faster. More constants could
 104   benefit from using efficiently computable sequences of narrower intervals
 105   for their continued fraction expansions.
 106
 107 - The API could use cosmetic surgery.
 108
 109 - The Makefile is fragile and annoying to maintain and use.
 110
 111 - The testing infrastructure needs upgrading.
 112
 113 - Much more, e.g. see TODOs sprinkled throughout the code.
 114
 115 == Design ==
 116
 117 I was inspired by Doug McIlroy's paper,
 118 ``http://swtch.com/~rsc/thread/squint.pdf['Squinting at Power Series']'', a
 119 captivating introduction to the communicating sequential processes (CSP) model
 120 of concurrent programming: some processes generate coefficients of a power
 121 series, while other processes read coefficients on input channels and write
 122 their output on another channel as soon as possible. Each output coefficient is
 123 written as soon as it is known, before reading more input. Like Unix pipes,
 124 these processes are connected together to compute sums, products, and so on.
 125
 126 For example, a process that repeatedly spits out 1 on its output channel
 127 represents the power series 1 + 'x' + 'x'^2^ + ... This is similar to
 128 running
 129
 130  $ yes 1
 131
 132 This process could be fed into another process that produces
 133 the derivative of its input. A Unix equivalent might be:
 134
 135  $ yes 1 | awk \'{ print (NR-1)*$0 }\'
 136
 137 Continued fractions likewise naturally fit the CSP model, as Gosper suggets.
 138 The simplest are processes that output terms according to some
 139 pattern while maintaining barely any internal state. Other processes can then
 140 consume these terms to compute convergents, decimal expansions, sums, products,
 141 square roots and so on.
 142
 143 Thus we spawn at least one thread per continued fraction or operation. These
 144 threads communicate via crude demand channels, sending terms of continued
 145 fractions back and forth.
 146
 147 For example, consider the
 148 code for generating 'e' = 2.71828...
 149
 150 ..............................................................................
 151 // e = [2; 1, 2, 1, 1, 4, 1, ...]
 152 static void *e_expansion(cf_t cf) {
 153   int even = 2;
 154   cf_put_int(cf, even);
 155
 156   while(cf_wait(cf)) {
 157     cf_put_int(cf, 1);
 158     cf_put_int(cf, even);
 159     even += 2;
 160     cf_put_int(cf, 1);
 161   }
 162
 163   return NULL;
 164 }
 165 ..............................................................................
 166
 167 The function +cf_put+ places a term on the output channel, and +cf_wait+ is our
 168 implementation of demand channels, a sort of cooperative multitasking. Without
 169 it, not only would we have to destroy and clean up after the thread ourselves,
 170 but more seriously, the thread might consume vast amounts of resources
 171 computing unwanted terms. The +cf_wait+ functions instructs the thread to stay
 172 idle. Our threads call this function often, and if it returns zero, our threads
 173 clean themselves up and exit.
 174
 175 Threads are the future, if not already the present. Multicore systems are
 176 already commonplace, and as time passes, the number of cores per system will
 177 steadily march upward. Happily, this suits continued fractions.