Copper Thoughts

A web λog by Isaac Hodes.

Fun With Set Theory and Lisp

A few weeks ago a provoking question was posted on StackOverflow asking if, for example, Common-Lisp's numberp could be written using only the primitives that John McCarthy presented in his now famous essay introducing LISP.

My answer was yes, but it wasn't the kind of yes you might expect. In his essay, McCarthy did not even mention numbers. There are no primitives to deal with them; only primitives to construct and deconstruct lists (i.e car, cdr, cons), as well as a few to enable the creation and naming of functions (i.e. lambda, label). But there are no number primitives: no direct interaction with the ALU, no way to directly access and manipulate integers and floats. So we have figure something else out.

It's All Lists From Here

So how can you represent numbers?

I'm going to eschew using the troublesome cons cell in my explanations. Instead, we'll use the more common list notation that's both easier to explain and easier to write. Of course, I'll be using the cons function (how could I not?) and you should at least have an idea of how that works.

It turns out that there are a number of ways to represent the natural numbers. Paul Graham postulates one way of representing them in his essay The Roots of Lisp. I'll state here, slightly tidied up, his suggestion:

  • Let the empty list be zero.
  • Let 1 = (cons 0 ())
  • Let 2 = (cons 1 ())
  • etc.

In this way we can represent all the natural numbers. It'd be trivial to write a function number? that would test to see if if its argument had the form of a number. There are a number (hah!) of things we could do with this representation, but I think we can do better.

Edit: Phil Hagelberg reminds me of another, very important, representation of the natural numbers I'd be remiss to leave out: the Church Encoding. Church Encoding is way of representing the natural numbers using only lambda calculus: there are no lists involved. Instead, a number n is representing by a higher order lambda abstraction (a function), and returns the result of that function being composed with itself n times. There's actually a fair amount that can be done with this representation, but it's outside the scope of this essay. For the curious, the afore-linked Wikipedia article does a fine job of covering the basics.

The More Lists, the Merrier

But really, that's not enough lists, is it? (((((((((()))))))))) is not nearly an unwieldy enough definition for 10. We need more lists.

Let's take a step back and look at set theory; a branch of math and CS that's informed a lot of underpinning of CS and programming in general. Set theory has already provided us with a handy way of defining numbers based on lists. Well not lists, sets, but sets are basically lists. In fact, contrary to popular/programming belief, sets can hold the same thing multiple times, but in set theory that means nothing. We'll say the same thing about lists: "when you say my list a contains b and b, I just hear that b is in a".

So how do we represent numbers using only sets? Well, we start with the empty set {} being zero. But, and here's where things get weird, we now define a function S called the "successor function". The successor of n is the union of n and the set of n, where union is a function of n and m that returns the set which contains all the elements in n and all the elements in m. S(n) = {n and all the elements in n}. To make this clearer, let's look at a few examples.

S(0) = union(0 {0}) = {0} ;; we'll call this 1
S(1) = union(1 {1}) = union({0} {{0}}) 
       = {0 {0}} ;; and this 2
S(2) = {0 1 2}  ;; 3
S(3) = {0 1 2 3} ;; 4

At the very least, this is an interesting definition of the natural numbers. The predecessors of each number are contained within it. And the set of all the natural numbers has a form similar to an individual natural number. Neat.

Let's see how we can represent this with lists. Our goal is a function, S, that returns the successor of a given list.

(define S
  (lambda (n)
    (union n (cons n ()))))

Perfect, though we don't exactly have a union function yet:

(define union 
  (lambda (n m)
    (cond ((null? n) m)
           (else (cons (car n)
                       (union (cdr n) m))))))

Now we do. You can test it out in your own REPL, or online with TryScheme. We're going to move on now.

So What?

Thats a valid question. At this point, it seems as though we've merely complicated Mr. Graham's elegant and simple proposition. It turns out that we've actually made things easier, in the long run.

Before we even get to addition and the like, let's look at one example of where we've made things easier and more elegant. Let's look at an ordering on the natural numbers; <. With this definition, we don't need anything else to define it: n < m if (and only if) n is in m. That's a straight forward definition, and one we can check rather quickly with a Lisp function you might have on hand anyway (The Little Schemer calls it member?).

But first we need a better eq?. The current eq? operates only on atoms; we need an equal? that works on our list representation of numbers. To keep things simple, we're going to assume our lists are ordered. That will be the case if we build our numbers with just S.

(define equal?
  (lambda (n m)
    (cond ((and (null? n) (null? m)) #t)
          ((or (null? n) (null? m)) #f)
          (else (and (equal? (car n) (car m))
                     (equal? (cdr n) (cdr m)))))))

(define member?
  (lambda (a l)
    (cond ((null? l) #f)
          (equal? a (car l)) #t)
          (else (member? a (cdr l))))))

(define <
    (lambda (n m)
       (member? n m)))

This works because all our natural numbers are either the empty list or a list of natural numbers, thanks to our nicely defined S function.

It turns out that this successor function is also nice for defining and working with ideas like ordinals and cardinals, and just generally makes out theorems easier to prove and reason about. But, without writing a book on it, suffice it to say there's a certain elegance about it. I'll leave it at that.

In Addition…

But back to our problem: numbers on a simple Lisp machine. We've now got a suitable representation of the natural numbers (though I've neglected to introduce the axioms required to fully justify this; let's just use common sense and carry on), and we have a way of comparing two numbers.

We also have a way of adding one to a number in S. How about adding two natural numbers together?

(define zero?
  (lambda (n)
    (equal? () n)))

(define A 
  (lambda (n m)
    (cond ((zero? m) n)
          (else (S (A n (P m)))))))

That is a nice recursive definition, but we have a mysterious new P function now, as well. P stands for "predecessor", and it gives us the natural number than m is a successor of.

Before defining P, let me present the set theoretical definition of addition:

  • A(n, 0) = n
  • A(n, S(m)) = S(A(n, m))

This is also recursive, but it has the unfair advantage of what the Clojurist in me calls the "de-structuring" of its operands. A can check to see if its second operand is a successor, and if it is will take the successor of the result of another A's value.

With Lisp, we can check to see is m is zero, and if it isn't, we know that it is a successor, we just don't know what it is a successor of. That is where P comes in.

(define R
  (lambda (fn acc l)
     (cond ((null? l) acc)
           (else (R fn (fn acc (car l)) (cdr l))))))

(define P
  (lambda (n)
    (cond ((zero? n) #f) ;; really, undefined
          (else (R union () n)))))

And out of nowhere comes a higher-order function! While set theory may have some expressive advantages over Lisp, we've got some special sauce, too. In this case (R union () n) is equivalent to the set theoretical "union over a set" that the Zermelo-Fraenkel Union Axiom gives provides. Most people call R "reduce".

You might notice that we get duplicates in our lists, but as I said before, we don't pay those any attention. To keep things clean and working well, we'll implicitly apply a helper function, distinct, to the result of every operation. distinct is left as an exercise to the reader; but this is a very important exercise: some of our function will not work if distinct isn't applied to the result of each function.

With a proper P defined, A now works flawlessly. Now let's define M, multiplication and B, subtraction:

(define M
  (lambda (n m)
    (cond ((zero? m) ())
          (else (A n (M n (P m)))))))

(define B                                                         
  (lambda (n m)
    (cond ((zero? m) n)
          (else (P (B n (P m)))))))

It gets easier and easier as we build off of our previous functions, and combine them to yield more complex, more useful abstractions. Finally, let's write our number predicate;

(define length
  (lambda (n)
    (cond ((null? n) ())
          (else (S (length (cdr n)))))))

(define number?
  (lambda (n)
    (equal? n (length n))))

And there you have it; we've got the basics of our natural numbers, represented using McCarthy's basic LISP.

But, Wait! There's More

Ah, we've left off the integers. And floats!

Fortunately, this isn't a big deal. With our newly defined natural numbers, we can represent integers as a tuple of natural numbers, where the first natural number in the tuple represents the "negative" portion of the integer, and the second natural number represents the "positive" portion of the integer.

In this manner we can represent all the integers. Of course, we'll have to create new functions for adding, subtracting and the like, but that becomes a matter of building new, higher level abstractions with our existing functions.

For example, we can define adding two integers like this:

(define first
  (lambda (n)
    (car n)))

(define second 
  (lambda (n)
    (car (cdr n))))

(define makeint
  (lambda (n p)
    (cons n (cons p ()))))

(define A2
  (lambda (x y)
     (makeint (A (first x) (first y)) 
              (A (second x) (second y)))))

Something important to keep in mind, however, is that different tuples can represent the same integer. For instance, the integer (2, 5) is the same as the integer (0, 3); they're both what we would have called "3" back in the natural numbers section. In fact, it makes sense to call both of them "3" here. But the second example here is preferable, I'd say. This is something that is dealt with in set theory by equivalence relations; we can do the same thing in Lisp. Below I'll define a function which returns the "simplest" version of a given integer; one where either the negative portion is zero, the positive portion is zero, or they are both zero.

(define simplint
  (lambda (z)
    (cond ((or (zero? (first z)) (zero? (second z))) z)
          (else (makeint (P (first z)) 
                         (P (second z)))))))

Subtraction, multiplication and other operations on integers begin to fall into place with these definitions. Things don't get more complicated as you move further away from your base representation (unless you try to formally prove all of my assertions): if anything, they get easier. This is a consequence of appropriate abstraction.

So far, we have arbitrarily large integers available to us (limited in practice only by memory and processing power). But we're still missing a type of number most programmers expect to have available to them: the float. Well, we can do that too.

First, let's make sure we know what a float can and cannot represent. Floats can't represent all real numbers: that is to say, every float is representable by a ratio of two numbers. In fact, traditional floats can't even represent all rationals: for instance, 1/3 is 0.333… repeating forever. We would run out of bits to represent it on a real computer. But we can represent all numbers a float can, and more, with ratios. And, as I'm sure has now become clear, we definitely have ratios in the form of a tuple of integers, the first representing the numerator and the second the denominator. These are the rational numbers.

Again, an equivalence relation on the rationals is called for, as are new operations. This time, though, they are left as an exercise for the reader.

Final Words

Though a lot more can be said about the intersection of set theory and Lisp, I'll stop while this post is still of Internet-length. If there's enough interest (mostly on my part), perhaps this will be extended into a series of essay on the subject, mirroring the progression of a elementary set theory book.

This essay is a product of my investigation into the possibility of a truly programmable and reprogrammable Lisp environment, one in which such a definition of the natural, integer and rational numbers would not only be feasible and efficient, but idiomatic. I'm a long way from explaining how such a system would work, but I thought I should record at least some of the byproducts of my thought process so that someone might obtain some value from my mental wanderings.


I should note that while I used the operators and and or, neither were in McCarthy's original LISP. Consider (and a b) a stand-in for (cond (a b) (else #f)) and (or a b) for (cond (a #t) (else b)), and all is well.

If there are errors in my explanations or code samples, please let me know. I certainly simplified much of the content so that it could fit in about 2000 words, but it should be correct.

The Golden Ratio is Evil

What is the Golden Ratio?

The Golden Ratio is a number that's been fascinating artists, mathematicians and philosophers for millennia. The fascination started with Mr. Phidias , circa 400 BCE, the sculptor involved in making the statues of the Parthenon. He considered the ratio a rule by which to craft the perfect human form, and it's after him that we know the Golden Ratio as φ, the Greek letter phi.

A few mathematicians and philosophers later, Fibonacci discovered that the ratio, from term to term, of his eponymous sequence approaches φ. Pretty cool, that.

Used in architecture, painting, drawing, and many works of adulatory prose, φ has garnered a lot of attention over the ages. But there's a darker side to φ

Calculating the Golden Ratio

As already mentioned, the Golden Ratio can be calculated using the Fibonacci sequence. In order to do that, though, we need the Fibonacci sequence. I've decided to make a lazy version of the sequence, so we can think of it as an infinite sequence. I split my fibo into two parts: fibo-from, which gives you the entire sequence after the point you give it, and fibo, which gives you the entire sequence.

(defn fibo-from [n-2 n-1]
  (lazy-seq
   (cons (+ n-2 n-1)
         (fibo-from n-1 (+ n-2 n-1)))))

(defn fibo []
  (concat [0 1]
          (fibo-from 0 1)))

Now we can continue on to deriving φ: Let's look at the first 20 ratios:

(map (fn [[x y]] (/ y x)) 
     (drop 1 (partition 2 1 (take 22 (fibo)))))

If we map double (the function) across the ratios we get in return, we can cleary see that we're getting closer and closer to the Golden Ratio. Of course, we'll never actually reach it. On the subject of ancient Greece, you could say φ tantalizes us.

Another way of calculating our fickle ratio is with a simple continued fraction. I've generalized the simple-continued-fraction function from my previous post on finding e, but the idea is the same:

(defn generalized-continued-fraction
  "Returns the value of the GCF of the sequences of numerators
   and denominators. 'denomseq must have one more element than
   'numseq."
  [numseq denomseq]
  (let [sn (reverse numseq) sd (reverse denomseq) ]
    (loop [hn (first sn) sn (rest sn) hd (first sd) sd (rest sd)]
      (cond (empty? sd) hd
            :else (recur (first sn) (rest sn)
                         (+ (first sd) (/ hn hd)) (rest sd))))))

(defn simple-continued-fraction
  "Returns the final value of the SCF of 'sequence."
  [sequence]
  (generalized-continued-fraction
    (take (dec (count sequence)) (repeat 1)) sequence))

The sequence that represents the Golden Ratio is rather simple and beautiful, and also surprising if you haven't seen it before. It is 1, 1, 1, 1, . . . and so forth. Therefore using simple-continued-fraction to find it is trivial:

(double (simple-continued-fraction (take 1000 (repeat 1))))
=> 1.6180339…

Evil Lurks

But in this seemingly beautiful number, a dark secret hides. But first we need to understand what makes a number evil. It's a simple idea, really: if the sum of the first n digits of the fractal part of a decimal number is 666, the number is evil.

But in order to get precise results, we need to get a precise representation of the Ratio. Here, we get to harness the power of Java! In this case, we're going to use the BigDecimal class.

(let [n (simple-continued-fraction (take 1000 (repeat 1)))]
  (.divide (BigDecimal. (numerator n))
           (BigDecimal. (denominator n))
           200 ;; our precision
           BigDecimal/ROUND_HALF_UP)) 
=> 1.6180339887498948482045868343656381177203091798057628621
354486227052604628189024497072072041893911374847540880753868
917521266338622235369317931800607667263544333890865959395829
0563832266131992829026788M

Which looks about right.

Now we need to see about getting the digits from the fractal portion of the number. To save time, I'm going to cheat a little, and transform the result into a string and obtain the digits from that. Assume that I've stored the result of the previous calculation in the var gold.

(map #(Integer. (str %))
     (drop 2 (str gold)))
=> (6 1 8 0 3 3 9 8 8 7 4 . . . ) ;; and so on

Now let's see how many numbers we have to take to get the number of the beast. Assume that our sequence of fractal digits found is stored in the var goldseq.

(defn sum [& xs] (reduce + xs))
(loop [maybebeast (sum (take 1 goldseq)) n 1]
  (cond (= maybebeast 666) n
        :else (recur (apply sum (take (inc n) goldseq))
                     (inc n))))
=> 146

And there you have it, the sum of the first 146 fractal digits of φ is 666. Verified: φ is evil. Of course, if it had not been evil, or it had taken the sum of more than 200 fractal digits to equal 666, this wouldn't have worked. It would have been easy to test though: continue testing the sum of the first n digits of the number until the sum exceeds 666. (Which could fail if we continue to get zeroes. Per the halting problem, we might never know if a number is evil, depending on its fractal component. Unless we could prove something else about it. But that's another story…)

But it's not just the fact that φ reminds us of one of Hades' cruelest trials (the trial of Tantalus), or that the first 143 digits sum to 666. No, there is yet more evil to be found in φ. If you take the ratios of the sections of a pentagram, you'll again end up with the Golden Ratio. Evil everywhere!


In all seriousness, φ is a lot of fun to play with, and that's about it. If you were serious about calculating φ properly, you'd find the positive root the quadratic equation which defines it. That equation is:

…and solving for the positive root yields the Ratio.

Regardless, I find sequences and continued fractions a lot more interesting and conducive to play in Lisp with.

Let me leave you with this piece of advice: beware the evil phi!

The Clojurist's Guide to Java

Judging from the buzz on Twitter, the mailing list, Hacker News and other places where this kind of buzz occurs, Clojure's popularity is rapidly growing. The thing is, according to Chas Emerick's State of Clojure published earlier this summer, most of the people coming to Clojure don't have much experience with Java. This would be fine but, sooner or later, you'll think you need Java. How often you actually do to use Java is debatable; most of the time, something basic you find yourself needing may very well be in clojure.contrib or on Github, but chances are you'll be recommended a Java library or two to accomplish a task you may not know how to in Clojure. And that's okay.

Using Java from Clojure is actually pretty painless, contrary to what people who haven't tried doing it say, but first you need to have at least a basic grasp on how to use Java. If you've had any experience with another object oriented language, then getting started with Java is just a matter of pulling up your sleeves and searching through the JavaDocs for what you need. Of course, you may need a refresher too. And using Java from Clojure can be a little strange at first.

This guide doesn't aim to comprehensively cover all the ways to use Java from Clojure, nor does it attempt to explain the intricacies of Java the language or the JVM. Rather, it should serve as a "Getting Started with Java from Clojure" guide that will hopefully enable you to more easily navigate the Java documentation and use Java in your Clojure projects when the need arises. One of the nice things about Clojure is that you don't really need to know Java to use it from Clojure.

So without further ado:

A Clojurist's Guide to Java

Classes

A class is a bundle of methods (functions which act on the class) that can serve as a data type. For example, to create a new class of the type Double: (Double. 1.2) which initializes the class Double (the period after the class name is syntactic sugar for calling the class constructor methods, which initialize the class with the values you provide) with the value 1.2. (Note that we don't have to import or qualify the Double class as it's a part of java.lang, which is included by default in your namespace.)

Now, look at the Double's constructor documentation in the Java 6 API:

Double

public Double(double value)

Constructs a newly allocated Double object that represents 
the primitive double argument.

Parameters:
value - the value to be represented by the Double.

So you can see what happened there. You "built" a new Double with value 1.2, which is a double (the primitive). A little confusing there, but really a Double is a class that represents a double and allows you to do things relating to doubles with its methods.

Static Methods

For instance, to parse a Double value out of a string, we can use the static method (meaning we don't need a particular instance of Double, we can just call it like we called sqrt from the Math class) parseDouble(String s):

(Double/parseDouble "1.2") => 1.2

Nothing too tricky there. Notice how the argument of the method is a String, so we pass a string as an argument to the method.

Nonstatic Methods

Say we want to use a Java class that we initialized to something. Not too difficult:

(-> (String. "Hey there") ;; make a new String object
    (.toUpperCase)) ;; pass it to .toUpperCase (look up -> to see what it does)
                    ;; toUpperCase is a non-static method
=> "HEY THERE"

So now we've used a method which is not static, and which requires a real, live String object to deal with. Let's look at how the docs say it works:

toUpperCase

public String toUpperCase()

Converts all of the characters in this String to upper case 
using the rules of the default locale. This method is
equivalent to toUpperCase(Locale.getDefault()).

Returns:
the String, converted to uppercase.

So here we have a method which returns a string (as shown by the "String" after the public in the definition, and takes no parameters. But wait! It does take a parameter. In Python, it'd be the implicit parameter self: this is called this in Java. We could also use the method like this: (.toUpper (String. "Hey there")) and get the same result.

More on Methods

Since you deal with mutable data and classes in Java, you need to be able to apply functions to classes (instances of classes, really) and not expect a return value.

For instance, say we're dealing with a JFrame from the javax.swing library. We might need to do a number of things to it, not with it (you generally operate with values, not on them in functional languages). We can, like this:

(doto (JFrame. "My Frame!");; clever name
      (.setContentPane ... here we'd add a JPanel or something to the JFrame)
      (.pack) ;; this simply arranges the stuff in the frame
      (.setVisibleTrue)) ;; this makes the Frame visible

doto just passes its first argument to all the other functions you supply it, and passes it as the first argument to them. So here we're just doing a lot of things to the JFrame that don't return anything in particular. All these methods are listed as methods of the JFrame in the documentation (or its superclasses don't worry about those yet).

Alternatively, we could do something like this:

(def frame (JFrame. "I'm a Frame!"))
(do ;; do just evaluates everything in its body, in order
      ;; and returns the value of the final expression
    (.pack frame) ;; = frame.pack(); in Java
    (. setVisible frame true)0 ;; = frame.setVisible(true);

Notes

A few useful functions you should be aware of when playing with Java in Clojure:

  • (class thething) ;; returns the class name of thething; useful to see what you're actually operating on.
  • -> and ->> ;; threading macros; check out their documentation, along with doto
  • (instance? AClass thething) ;; true if thething is an instance of AClass

Wrapping up

You should now be prepared to explore the JavaDocs yourself. Here you'll find everything that is available to you in a standard Java 1.6 install. There will be new concepts, but a quick Google search should answer most of your questions, and you can always come back here with specific ones.

Be sure to look into the other important Clojure functions like proxy and reify as well as extend-type and its friends. I don't often use them, but when I need to, they can be invaluable. I still am understanding them myself, in fact. There's a ton out there, but it's mostly a problem of volume rather than complexity. It's not a bad problem to have.

Additional reading:

New Blog

I've broken some of the web.

Let me apologize now: I'm sorry. Now let me explain.

The original Copper Thoughts was a way of staking out a claim on the Internet. I wanted a place to publish stuff, so I threw up a Wordpress installation and let it go at that, telling myself that one day I would take the time to create a "proper" custom blogging engine. A couple months ago, I finally did. It was in Clojure, and it was a neat idea, but it wasn't perfect, so I scrapped it and continued to use Wordpress. A couple weeks ago, I switched over the post structure of the blog to something more intelligible, but broke links in the process. It wasn't a big deal, as most of those posts got little to no traffic, so it didn't matter much. The links will be safe from this point on, but I've only keep the past couple posts around. The rest have been lain to rest.

For now, until I've created the perfect baked blogging engine, I will be posting raw HTML of each post when I'm ready to publish them, and will update the archive and index pages by hand. This should at least provide some measure of impetus for me to quickly complete Clogs (the baked blogging engine).

But know this, Gentle Reader: I won't break any more of the web.

E and Clojure

A few articles made the rounds a number of weeks ago, all talking about computing Euler's Number, e, using Clojure. It wasn't until I saw the idea on Programming Praxis, though, that I decided to just do it.

It's honestly a trivial problem; in fact, here's a little solution right up front:

(defn to-one []
  (loop [times 0 sum 0]
    (cond
       (> sum 1) times
       :else (recur (inc times) (+ sum (rand))))))
<p>(defn to-one-avg [n]
  (/ (reduce + (repeatedly n to-one))
      n 1.0))</p>
user> (time (to-one-avg 100000000)_
"Elapsed time: 59539.291 msecs"
2.7182 5238

Error as compared to Math/E: 0.000067438. Not great.

That's a solution based off of the cool fact that, on average, of the numbers of random numbers between 0 and 1 is, well, 2.7182… e. (See this blog post for a little explanation.) The problem is, it takes a lot of time to get an only moderately accurate answer. Even doing it in C results in less than stellar results, and still sluggish performance:

See it in a Gist, as it takes up too much room.

ihodes@Machine:~/Desktop
>> time ./efinder 100000000
Our estimate for e is 2.7182752500.
Total time: 0m12.386s

That's not so hot. There's got to be a better way.


An alternative for finding e can be represented nicely as the sum of a sequence:

sumfore

And it can be repreented just as nicely in Clojure! My sigma and factorial (fact) functions are not standard in Clojure, but they're both quite straightforward. Not only is it a pretty solution, but it's quite fast, and is just as accurate as our built-in Math/E.

(defn sigma
  "Sum of 'f across the range 's to 'e, inclusive."
  [s e f]
  (reduce + (map f (range s (inc e)))))<p>(defn fact
  "Returns the factorial of 'n."
  [n]
  (reduce * (range 1 (inc n))))</p>
user> (time (double (sigma 0 100 #(/ 1 (fact %)))))
"Elapsed time: 7.998 msecs"
2.718281828459045

Error as compared to Math/E: 4.440892098e-16. Not bad.

But we've barely scratched the surface of how much of e is possible to calculate, and how quickly.


Yet another method for finding e involves continued fractions. Lucky for us, the continued fractions involved are simple. No, really, that's what they're called. The one we're after looks like this:

continuedefrac

One thing to notice is that we can represent that continued fraction as a sequence. The the numerators are always 1, so we just need a list of the integer part of the denominators. In this case, [1 0 1 1 2 1 1 4 1 1 6 1 1 8 1 1 … ]. That's fun enough to build in Clojure:

(defn efracs
  [n]
  (cons 1 (interleave
                (filter even? (range n))
                (repeat 1) (repeat 1))))

The trick becomes using this sequence build our continued fraction and calculate e. What first comes to mind is writing a recursive function that would destroy our stack. But instead of building it from the top down, let's think backwards. Let's instead construct our continued fraction from the bottom up, and end up at our final answer with a well-behaved stack.

(defn simple-continued-fraction
  "Returns the final value of the SCF of 'coll."
  [coll]
  (let [coll (reverse coll)]
    (loop [coll (rest coll) h (first coll)]
      (cond (= 1 (count coll)) (+ (first coll) (/ h))
            :else (recur (rest coll) (+ (first coll) (/ h)))))))
(time (simple-continued-fraction (efracs 30)))
"Elapsed time: 0.986 msecs"
2.718281828459045

Scrumtralescent.


Now, I hate to break it to you, but these are not the methods used to calculate e to 1e18 digits. You can't even store that many digits in RAM; instead, you need to calculate e digit by digit, and write each digit to a file on your hard drive.

There's a particular program called PiFast33 that was written by Xavier Gordon around the year 2000. That program, or variants of it, is what is being used to calculate e with such great precision. According to this site, which appears to be connected to Gordon in some fashion, PiFast33 uses the alternating version of the series described above to calculate e. The error shrinks rapidly, allowing the digits of e to be verified and written to the storage file quickly and efficiently.

It looks like an interesting program to try to write; maybe I'll get the chance to hack it up and post about it soon.