Calculating π in Clojure (Salamin-Brent)

Took a shot at implementing a PI digit generator in Clojure using a ‘fast’ algorithm.
It seemed like a decent enough excercise to try and understand something about performance in Clojure.

MacBook Pro – Intel Core 2 Duo 2.26 GHz – 4GB RAM
Java(TM) SE Runtime Environment (build 1.6.0_15-b03-219)
Java HotSpot(TM) 64-Bit Server VM (build 14.1-b02-90, mixed mode)
Clojure 1.1.0-alpha-SNAPSHOT (Aug 20 2009) git commit f1f5ad40984d46bdc314090552b76471ee2b8d01

Clojure matches the performance of Java in this example.

The Clojure code :

(import 'java.lang.Math)
(import 'java.math.MathContext)
(import 'java.math.BigDecimal)

(defn sb-pi [places]
  "Calculates PI digits using the Salamin-Brent algorithm
   and Java's BigDecimal class."

  (let [digits (.intValue (+ 10 places)) ;; add some guard digits
        round-mode BigDecimal/ROUND_DOWN]

    (letfn [(big-sqrt[#^BigDecimal num]
             "Calculates square root using Newton's method."
             (letfn [(big-sqrt-int 
                      [#^BigDecimal num #^BigDecimal x0 #^BigDecimal x1]
                      "aux function for calculating square root"
                      (let [#^BigDecimal x0new x1
                            #^BigDecimal x1new (-> num (.divide x0new digits round-mode))
                            #^BigDecimal xsum (+ x1new x0new)
                            #^BigDecimal x1tot (-> xsum (.divide 2M digits round-mode))]
                        (if (= x0 x1)
                          (recur num x1 x1tot))))]
                num 0M (BigDecimal/valueOf
                        (Math/sqrt (. num doubleValue))))))
            (sb-pi-int [#^BigDecimal a #^BigDecimal b 
                        #^BigDecimal t #^BigDecimal x #^BigDecimal y]
             "aux function for calculating PI"
                 [#^BigDecimal y1 a
                  #^BigDecimal absum (+ a b)
                  #^BigDecimal a1 (-> absum (.divide 2M digits round-mode))
                  #^BigDecimal b1 (big-sqrt (* b y1))
                  #^BigDecimal ydiff (- y1 a1)
                  #^BigDecimal t1 (- t (* x ydiff ydiff))
                  #^BigDecimal x1 (* x 2M)]               
               (if (== a b)
                 (let [#^BigDecimal absum1 (+ a1 b1)
                       #^BigDecimal absqrd (* absum1 absum1)
                       #^BigDecimal u (* t1 4M)]
                   (-> absqrd
                       (.divide u digits round-mode)
                       (.setScale places round-mode)))
                 (recur a1 b1 t1 x1 y1))))]
      (sb-pi-int 1M (-> 1M (.divide #^BigDecimal (big-sqrt 2M) digits round-mode))
                       (/ 1M 4M) 1M nil))))

(time (println (sb-pi (Integer/parseInt (second *command-line-args*)))))

$ time clj pi.clj 1               -->       3.403 msecs
$ time clj pi.clj 10              -->       3.956 msecs
$ time clj pi.clj 100             -->      10.630 msecs
$ time clj pi.clj 1000            -->     141.937 msecs
$ time clj pi.clj 10000           -->    3316.180 msecs

The same algorithm in Java (but using iteration instead of recursion) :

import java.math.BigDecimal;
import static java.math.BigDecimal.*;

class Pi {
  private static final BigDecimal TWO = new BigDecimal(2);
  private static final BigDecimal FOUR = new BigDecimal(4);
  private static int ROUND_MODE = ROUND_DOWN;

  public static void main(String[] args) {
    long start = System.nanoTime();
    System.out.println("Elapsed time: " + 
                       ((System.nanoTime() - start) / 1E6) + " msecs");  
  // Salamin-Brent Algorithm
  public static BigDecimal pi(final int digits) {
    final int SCALE = 10 + digits;
    BigDecimal a = ONE;
    BigDecimal b = ONE.divide(sqrt(TWO, SCALE), SCALE, ROUND_MODE);
    BigDecimal t = new BigDecimal(0.25);
    BigDecimal x = ONE;
    BigDecimal y;
    while (!a.equals(b)) {
      y = a;
      a = a.add(b).divide(TWO, SCALE, ROUND_MODE);
      b = sqrt(b.multiply(y), SCALE);
      t = t.subtract(x.multiply(y.subtract(a).multiply(y.subtract(a))));
      x = x.multiply(TWO);
    return a.add(b)
      .divide(t.multiply(FOUR), SCALE, ROUND_MODE)
      .setScale(digits, ROUND_MODE);
  // square root method (Newton's)
  public static BigDecimal sqrt(BigDecimal A, final int SCALE) {
    BigDecimal x0 = new BigDecimal("0");
    BigDecimal x1 = new BigDecimal(Math.sqrt(A.doubleValue()));
    while (!x0.equals(x1)) {
      x0 = x1;
      x1 = A.divide(x0, SCALE, ROUND_MODE);
      x1 = x1.add(x0);
      x1 = x1.divide(TWO, SCALE, ROUND_MODE);
    return x1;

$ time java Pi 1         ---->         2.162 msecs
$ time java Pi 10        ---->         2.425 msecs
$ time java Pi 100       ---->         7.897 msecs
$ time java Pi 1000      ---->       150.610 msecs
$ time java Pi 10000     ---->      3009.705 msecs

Clojure : The zombie-reanimated corpse of Lisp

* updated 2009-08-21 ** updated 2009-08-20 *


Nowadays it’s frowned upon to criticize anything in society. It’s quite pervasive.
That is unless you happen to also be criticizing the same thing that a group of other people are also criticizing.
We hardly speak our minds, instead it’s more common to take a passive-aggressive and sometimes socio-pathic approach when any type of issue arises. Annoying as they can be there just aren’t enough Zed’s in the world. We’d rather be part of the herd and not admit it to ourselves. We don’t want to be shunned, we don’t want to lose friends, our homes, wives and the stuff that we accumulate.
So we go with the flow at the expense of the dream and what is better in the long term. After
all everyone’s been doing just that since the beginning of time right ? And besides
won’t we make more money if we just go along with everyone else ? Hardly, look at
the current state of the world economy. Astonishing even as we’re told to buy that it’s looking better, it has to, we don’t want to believe the alternative.


On that note, this is my satiric rant regarding what many are hyping as the future of Lisp, so be warned.
I dedicate it to the Zed’s of this world. It’s abrasiveness is meant to counter the ‘mega-hype’ that accompanies Clojure’s sales pitch. Clojure is not worthy of the Lisp name. It is the antithesis of the Lisp machine in every way. It has little regard for the intellectual rigor and cohesiveness of Lisp.

Whether it’s author and followers are even aware or willing to acknowledge this :
Clojure has inherited many of Java’s idioms and problems.
Clojure is not only built on top of the JVM. It is also built on the Java language.

For better or worse this means that many Clojure types are Java types underneath.
This is a fundamental and idiomatic clash in Clojure’s type system and this effects the way that software can be written in Clojure. Java is an imperative language while Clojure is a ‘trying much too hard to be pure’ functional language. Clojure is hyped for it’s immutable strictness of local variables to the extent that some refer to it as if it’s a ‘silver bullet’ for issues related to writing software for multi-core platforms (STM etc.).

It’s dictator doesn’t provide much of a reason for enforcing the immutability of local variables. Yes, I’ve read the long thread on CLL that ‘explains’ some of the reasons. I’m still not convinced why concurrent programming support had to be implemented at the language level. The focus on concurrent programming assumes much in the way that immutable types are forced on the programmer unless they decide to stoop to using Java’s mutable types.

Why not simply add support for concurrent programming as a library instead ?

Why not provide both immutable and mutable versions of the same data types ?

Why not let the programmer decide what paradigm works best for themselves ?

The reasoning for implementing Clojure on the JVM is understandable, there is much that the JVM has to offer : memory models, GC, overall maturity etc.
The hype about ‘getting all that for free’ is very optimistic. The JVM was designed specifically for the Java language, so draw your own conclusions.

That seems to be changing more recently
. Supporting other languages on the JVM may turn out to be ‘a good thing’ depending on what your favorite language is on that platform.

There are an abundance of Java libraries that Clojure fans will be glad to have access to now.
The hype about having access to all these libraries ‘for free’ is ironic. Many Java libraries only exist
due to limitations of the Java language itself. Also when writing real-world software, what about the effort required to align the multitude of Java libraries that assume an imperative environment with Clojure ? Can these libraries be used easily, safely and with the same performance ?
You just have to picture all the newbie Clojure programmers when they hit this wall.
No problem right ? They can just wrap all the imperative Java libraries using their newly discovered hammer : macros. Hopefully many of these old libraries are used only as a bridge until they are replaced by more shorter idiomatic functional versions.

Clojure introduces new problems(ok ‘challenges’ if you will) having discarded enough of the last half-century
of Lisp progress. Sure, Common Lisp has it’s problems but hailing Clojure as it’s future is a new low in the community. Another lisper put it so bluntly :  “Clojure is the gutted and zombie-reanimated corpse of a dream.”

If anything, Clojure is hopefully the future of the JVM. It is definitely a step up the language continuum, for those who decide to move from Java to Clojure.

The hype, guru-worship, immaturity of the language and the benevolent dictatorship of Clojure are rapidly catalyzing something in the world of commercial software development. This can be an indicator of trouble ahead. There’s good reason that Haskell’s motto has long been “Avoid success at all costs.”
Clojure is popular right now and in a few years it’ll be used in common “tech” speak by clueless managers
and inexperienced programmers who consider Clojure their first exposure to a functional “Lisp”.
Many of us will inevitably be using Clojure in our day jobs in the same way that we currently have little choice in
using whatever the required and popular lowest common denominator is at a company.
At least Clojure is a real improvement over the status quo, I can appreciate that.

Another slice of π ?

The previous Common Lisp and Haskell functions to generate the digits of PI where only accurate between 10000 and 20000 digits. The algorithm uses an approximation where we discard a certain number of ‘guard’ digits to get an accurate result. Some background regarding how the number of guard digits is calculated : There is ‘wobble’ in the number of contaminated digits depending on the number of input digits. When only the order of magnitude of rounding error is of interest, ulps(units in the last place) and ε (machine precision) may be used interchangeably, since they differ by at most a factor of β (the ‘wobble’ in the number of contaminated digits). For example, when a floating-point number is in error by n ulps, that means that the number of contaminated digits is about logβn. [1] [2]

Here are new versions of the functions using a guard digits calculation :

Common Lisp :

  "Accurately calculates PI digits using Machin's formula
   with fixed point arithmetic and guard digits."
      ((arccot-recur (xsq n xpower op)
         (let ((term (floor xpower n))
               (opfun (if op #'+ #'-)))
           (if (= term 0)
               (funcall opfun
                        (arccot-recur xsq (+ n 2) (floor xpower xsq) (not op))
       (arccot (x unity)
         (let ((xpower (floor (/ unity x))))
           (arccot-recur (* x x) 1 xpower t))))

    (if (> digits 0)
        (let* ((guard (+ 10 (ceiling (log digits 10))))
               (unity (expt 10 (+ digits guard)))
               (thispi (* 4 (- (* 4 (arccot 5 unity)) (arccot 239 unity)))))
          (floor thispi (expt 10 guard))))))

Haskell :

Accurately calculates PI digits using Machin's formula
with fixed point arithmetic and variable guards digits.

arccot :: Integer -> Integer -> Integer
arccot x unity = arccot' 0 (unity `div` x) 1 1
      arccot' sum xpower n sign
          | xpower `div` n == 0 = sum
          | otherwise           =
              arccot' (sum + sign * (xpower `div` n))
                      (xpower `div` (x * x)) (n + 2) (- sign)

machin_pi :: Integer -> Integer
machin_pi digits =
    if digits < 1 then 0 else
        pi' `div` (10 ^ guard)
              guard = 10 + (ceiling (logBase 10 (fromInteger digits)))
              unity = 10 ^ (digits + guard)
              pi' = 4 * (4 * arccot 5 unity - arccot 239 unity)

“What Every Computer Scientist Should Know About Floating-Point Arithmetic”, by David Goldberg, March 1991

1. Relative error and Ulps

2. Guard digits