Sunday, June 23, 2013

Calculating Fibonacci Numbers

Yeah, I know, it's been done, it's boring, nothing new to see here. I understand, this is probably all known, but I just wanted to show this.

As I'm sure you are all aware, the Fibonacci sequence can be defined recursively via the equation:

\[ F(n) = F(n-1) + F(n-2)~~~~~~~F(n) = n~~\forall~~n \in \{0, 1\} \]

Now, this is fine and good, but if you translate this formula into a recursive algorithm in the most naive way possible, it is also extremely computationally expensive. The fact that each number requires the knowledge of the two previous numbers means that a naive algorithm to compute this function runs in exponential time. The Lispy way (IMHO) to deal with this exponential explosion is to use memoization as a dynamic programming method. This will cache the previous results and turn this function into a linear time algorithm. Another common way to compute this is to run through the entire sequence until you get the value you are looking for. This also runs in linear time but (usually) with a smaller constant factor and has constant space requirements. Unfortunately, it requires a substantial change from the intuitive description of the problem, making the code much harder to understand, which is why I don't like it as much as the memoized recursion. However, in the name of clarity when exploring the complexity, I will be using the "sum up from the base case method".

But, anyway, my point was that you can easily write a linear time Fibonacci number calculator. That is what we would say under normal circumstances, but this isn't really true, or rather, it is only true if you assume constant time multiplication. For large \(n\), it becomes clear that this algorithm is not, actually, linear. This is because it takes \(O(N)\) time to add two \(N\) digit numbers, and because at least some of the numbers that we wish to add are of order \(F(n)\) (which have \(O(\log F(n)) \sim O(n)\) digits), we see now that the actual time is more like \(O(n\log F(n)) \sim O(n^2)\). This might seem a bit surprising, but look and see for yourself.

A quadratic fit to the time scaling of the Fibonacci sequence generator function

But, if we are willing to throw clarity out the window all together, we can also use a different technique to compute the number. We can use the closed form solution to the Fibonacci sequence, Binet's formula (all of this is from Wikipedia):

\[ F(n) = \frac {\psi^n - \phi^n}{\psi - \phi} \]

…where…

\[ \psi = \frac{1+\sqrt 5}{2} \\ \phi = \frac{1-\sqrt 5}{2} \]

And, because \(\psi = -\phi^{-1}\), we can simply say:

\[ F(n) = \frac {\phi^n - (- \phi)^{-n}}{\sqrt 5} \]

Now, I hear you thinking, these aren't going to work if you actually want the exact number, because these have nasty things like \(\sqrt 5\) and and \(\phi^{-n}\) in them. And you would be correct that this is not something that you can simply type into any programming language and expect a valid number back. But the fact of the matter is, these formulas are actually correct, it is the language and the language's insistence that computation on general real numbers be performed using approximate floating point arithmetic that causes such an issue in the first place. If, however, we were to turn to the venerable bc, an arbitrary precision calculator that is already on your computer if you have a UNIX like operating system, then we can try it out for various \(n\) by setting the precision to a "high enough" value (the scale variable sets how many decimal digits of precision you have).

scale=100

phi = (1+sqrt(5))/2

define fib (n) {
    (phi^n - (-phi)^(-n))/sqrt(5) }

fib(5)
4.999999999999999999999999999999999999999999999999999999999999999999\
9999999999999999999999999999999988

fib(10)
54.99999999999999999999999999999999999999999999999999999999999999999\
99999999999999999999999999999999728

fib(100)
354224848179261915074.9999999999999999999999999999999999999999999999\
999999999999999999999999999999981555491928230697552907

fib(200)
280571172992510140037611932413038677189524.9999999999999999999999999\
99999999999999999999999999999997069407044903016136967720432076973572\
3021244

fib(400)
17602368064501396646822694539241125077038438330449219188672599289657\
5345044216019674.999999999999996317359669914941458065726733504366353\
6782163342769596482054432795143577733011851910134

fib(500)
13942322456169788013972438287040728395007025658769730726410896294832\
5571622863290691557658876222517646901.483330817850463056986072289669\
50027087266539483120658105031310816022996182603489404860217143670694\
09

We see that as the numbers are computed in the most significant digits, the error accumulates in the least significant. Once the error meets the important numbers, which has happened at the 500th Fibonacci number, we have reached the point where we actually don't know if the last digits are correct or not. At that point we need to increase the scale and try again.

Also, if you try this at home, you will find that the computation takes longer for larger values of \(n\), despite this being a closed form expression. This isn't really surprising, the closed for expression may seem magical, but it still needs to be implemented with algorithms that multiply, divide, exponentiate, and "square root" very large or very precise numbers. More on this later.

Incidentally, the CLISP implementation of Common Lisp also has arbitrary precision arithmetic and can be used to perform the same calculation much faster (also, you can also do this in other imps with some success using rational numbers if you take the time to write a sqrt function that will return a rational approximation of the square root to within the needed accuracy). If you try this at home using CLISP, note that the scale value in bc designates decimal digits while (long-float-digits) in CLISP designates binary digits (i.e. the CLISP one needs to be larger for the same precision. For example, CLISP needs a value of 3322 to give around 1000 decimal digits of precision).

(setf (long-float-digits) 350)

(let ((phi (/ (+ 1 (sqrt 5L0)) 2L0)))
  (defun fib (n)
    (/ (- (expt phi n) (expt (- phi) (- n)))
       (sqrt 5L0))))

(fib 100)
3.542248481792619150750000000000000000000000000000000000000000000000000\
000000000000000000000000000000000005L20

(fib 200)
2.805711729925101400376119324130386771895250000000000000000000000000000\
0000000000000000000000000000000000084L41

(fib 400)
1.760236806450139664682269453924112507703843833044921918867259928965753\
4504421601967500000000000000000000109L83

(fib 500)
1.394232245616978801397243828704072839500702565876973072641089629483255\
7162286329069155765887622252129412605L104

This works okay, but we have to constantly check to make sure that we haven't run out of precision in our number. We can fix this by computing the number of significant digits needed to accurately compute \(F(n)\) by first computing an approximate \(F(n)\) with floating point numbers and using that approximate value as a guess of how many \(b\)-ary digits you need in have by taking the \(\log_b\).

(defun fib-approx (n)
  (let ((phi (/ (+ 1 (sqrt 5d0)) 2d0)))
    (/ (- (expt (float phi 0d0) n) (expt (- (float phi 0d0)) (- n)))
       (sqrt 5d0))))

(defun fib (n)
  (setf (long-float-digits) (max 64 (round (* 1.5 (log (fib-approx n) 2)))))
  (let ((phi (/ (+ 1 (sqrt 5L0)) 2L0)))
    (nth-value 0 (round (/ (- (expt phi n) (expt (- phi) (- n)))
                           (sqrt 5L0))))))

(fib 5)
5

(fib 100)
354224848179261916032

(fib 1000)
43466557686937456435688527675040625802564660517371780402481729089536\
55541794905189040387984007925516929592259308032263477520968962323987\
33224711616429964409065331879382989696499285160037044761377951668492\
28875

However, this only really works up to the point when the double floats start to overflow (Fibonacci numbers less than around 300 digits). After that our fast approximation no longer works as the Fibonacci numbers are too big to fit into any machine sized number. Is there a way to just get an approximate ballpark figure for the size of a Fibonacci number without actually calculating it? Yes. All we need to do is look at the what is important in the equation when \(n\) is large. If \(n\) is large, \((-\phi)^{-n}\) is extremely small since \(\phi^{-1}\) is smaller than 1. This means that we can approximate the large \(n\) Fibonacci numbers as:

\[ F(n\gg 1) = \frac{\phi^n}{\sqrt 5} \]

Now, we know that the number of \(b\)-ary digits in a number is roughly equal to the \(\log_{b}\) of the number. This means that the number of digits in a large Fibonacci number is:

\[ \log F(n \gg 1) = \log \phi^n - \log \sqrt 5 = n\log\phi - \log\sqrt 5 \]

Adding this to our Fibonacci number calculator:

(defun log-fib-approx (n &optional (base 2))
  (let ((phi (/ (+ 1 (sqrt 5d0)) 2d0)))
    (- (* n (log phi base)) (log (sqrt 5) base))))

(defun fib (n)
  (setf (long-float-digits) (max 64 (round (* 1.5 (log-fib-approx n 2)))))
  (let ((phi (/ (+ 1 (sqrt 5L0)) 2L0)))
    (nth-value 0 (round (/ (- (expt phi n) (expt (- phi) (- n)))
                           (sqrt 5L0))))))

(fib 5)
5

(fib 100)
354224848179261915075

(fib 1000)
43466557686937456435688527675040625802564660517371780402481729089536555417\
94905189040387984007925516929592259308032263477520968962323987332247116164\
2996440906533187938298969649928516003704476137795166849228875

(fib 10000)
33644764876431.......84058580903081179711592936853976188104504236323110912

At this point, we might ask how all this compares to calculating the number via the recurrence equation? On the one hand, using intelligent algorithms we have found an \(O(n^2)\) time algorithm based on the recurrence relation. On the other, it isn't totally clear what the scaling behavior of the proposed closed form algorithm is. In the Binet's formula based method, we have (roughly in order carried out in the calculation):

  1. An arbitrary precision square root
  2. An fixed to arbitrary precision addition (adding one to a value)
  3. An arbitrary precision to fixed division (divide by two)
  4. Two arbitrary precision exponentiations
  5. An arbitrary precision subtraction
  6. An arbitrary precision division

Looking at this list, the obvious prime suspects for the bottleneck are the square root and the exponentiation, but let's look at each. Wikipedia tells me that multiplication and division have the same complexity lower bounded with \(O(N\log N)\) (probably) and that addition and subtraction run in \(O(N)\) time.

Using simple divide and conquer exponentiation gives an \(O(N\log N \log n)\) time scaling for exponent \(n\). The only thing that is left is the scaling of the square root. Looking through the GMP docs tells me that this can be done in \(O(N\log N)\) time (again assuming the Schönhage-Strassen conjecture for multiplication). This means that the most expensive part of the calculation is the exponentiation, at \(O(\log F(n)\log \log F(n) \log n)\) time. What does that mean for our algorithm?

The Fibonacci sequence grows in an exponential manner. This is clear from recurrence, \(F(n+1) = F(n) + F(n-1) \sim 2F(n)\). Roughly speaking, increasing \(n\) by a constant factor has the effect of multiplying \(F\) by a constant factor, which is basically the definition of "exponential". But, in our complexity, We are taking the \(\log\) of that number, which "undoes" the exponential function. More concretely, we actually calculated a large \(n\) approximation of \(\log F(n)\) above, which turned out to scale linearly with \(n\). All of this boils down to the relationship I mentioned before: \(O(\log F(n)) \sim O(n)\). Thus, we expect this calculation to be \(O(n (\log n)^2)\), whereas the complexity for the recurrence based results took \(O(n^2)\). This result demonstrates that Binet's formula is, in fact, an asymptotic winner in this computational contest. When you plot the timing data, it seems very linear (I certainly cannot seem to resolve the \((\log n)^2\) portion of the scaling behavior).

If you are actually interested in calculating Fibonacci numbers very quickly and are using CLISP, then it is probably a safe bet to use the Binet's formula method rather than the more intuitive algorithm even for small numbers. On CLISP, the differences in run times are around two orders of magnitude faster than the intuitive method (although it seems to be much less predictable in terms of run time in any particular run, my prime suspect on this is the GC). Using SBCL (a much more sophisticated compiler), the intuitive algorithm runs significantly faster than on CLISP, but still is no match for Binet's method (on CLISP).

Thursday, June 6, 2013

Some Customizations to W3m-Mode

For many tasks, I find that it is easier to use a text mode browser in Emacs than it is to use an external browser. This is clearly the case when you would like to copy some text from a web page into an Emacs buffer, like when you are grabbing code examples from a manual or tutorial. One of the main uses I have of this is to quickly browse the Common Lisp Hyperspec. I have found that W3m and W3m-Mode is quite good for these tasks. W3m is a text mode browser that exists independently from Emacs. W3m-Mode is a mode that connects to W3m, sends requests from Emacs to the external process, and receives the output and displays it in an Emacs buffer. Edit: If you want to set Emacs up to use this browser, or any browser, customize the browse-url-browser-function variable.

Of course W3m won't work with a lot of "modern" websites, anything that uses Flash (of course), or anything that utilizes Javascript extensively. Even with all of those barriers, it works surprisingly well because much of the information that is out there on the Internet (and most of the information that I consume using this browser) is actually just text. Further, thanks to HTML5 and the return of the separation of content (HTML) and design (CSS) many websites can be rendered usefully in text.

However, as much as I have learned to like W3m-Mode, there are some annoying aspects of the mode that are the default. Namely, I would like to be able to move through the buffer with the standard directional keys, as well as switch tabs in the way that most other browsers do via "C-<tab>", and move forward and backward in history via "M-<right>" and "M-<left>", respectively. With a little help from the Emacs folks on /r/emacs, here are some bindings that will set this up for you.

;; W3m for browsing
(defun w3m-mode-options ()
  (local-set-key (kbd "<left>") 'backward-char)
  (local-set-key (kbd "<right>") 'forward-char)
  (local-set-key (kbd "<up>") 'previous-line)
  (local-set-key (kbd "<down>") 'next-line)
  (local-set-key (kbd "M-<left>") 'w3m-view-previous-page)
  (local-set-key (kbd "M-<right>") 'w3m-view-next-page)
  (local-set-key (kbd "C-<tab>") 'w3m-next-buffer)
  (local-set-key (kbd "q") 'bury-buffer))
(add-hook 'w3m-mode-hook 'w3m-mode-options)

That last local-set-key is to rebind the "q" key. By default, "q" tells W3m-Mode to kill all W3m sessions and close the window. This isn't the right thing to do for a couple of reasons. First, if I want to kill all W3m sessions, I will do so by killing the buffer (which will kill the sessions). I don't need a separate key for that. Second, closing windows without the user expressly asking for a window to be closed is not good Emacs interface style. The correct thing to do is to go away (whatever that means) and leave the window to be filled with another buffer. The only exception to the rule is popup windows, windows that are meant to be extremely short lived (less than a minute) and which created the window in which the reside. W3m is most definitely not a popup window. Looking at this, I felt that the best solution is to map "q" to bury-buffer and leave it at that.

There are a few other bindings that are different from the standard key bindings that a browser uses, including making a new tab and kill that tab. These exist and are reasonable (they fit in with the standard Emacs user interface), just check the help for W3m-Mode via "C-h m".

Saturday, June 1, 2013

A Multiple Cursor Trick and an Improvement

I have been baking, in the back of my mind, a way to make multiple cursors more powerful. What I wanted was a way to move through a buffer, and mark certain places where I want to place my multiple cursors, and then make my changes. I think this idea is basically baked enough to reveal that… it is not actually necessary to modify the multiple-cursors library to achieve this goal.

If you want to do this, you can simply insert a sequence of characters that is unique at each point that you wish to place a cursor (something that you would never use, like some crazy Unicode character, say "ಠ", you just have to find an easy way to insert it into the buffer). Then, select that tagging sequence and use mc/mark-all-like-this and edit away. This actually works in a pinch. For instance, you can do things like this:

However, I guess that this is slightly more hackish than I like, so I came up with a new, better method. What I decided would be pretty awesome is to have multiple-cursors be able to mark spots off the top of the mark ring (either via popping or just walking the mark ring). I defined a function mc/mark-pop that will just this. I use the suggested multiple cursors and expand region key bindings and bound mc/mark-pop to "C-S-p" (which makes sense on a Dvorak keyboard layout, if you use Qwerty, you might use "C-?"). This means that you can do some pretty awesome stuff like this:

I have submitted a pull request to Magnars which has been accepted. This is still a bit inconvenient to use, but it has promise for being an excellent building block for future multiple cursor editing tricks.