Re: Calculating digits of π (Salamin-Brent)

[jng27]

This updated version is 2x as fast as the previous version :

(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 (+ 10 places) ;; add some guard digits
        round-mode BigDecimal/ROUND_DOWN]

    (defn big-sqrt[#^BigDecimal num]
      "Calculates square root using Newton's method."

      (defn big-sqrt-int
        [#^BigDecimal num
         #^BigDecimal x0
         #^BigDecimal x1]
        (let [x0new x1
              x1new (. num divide x0new digits round-mode)
              x1tot (. (+ x1new x0new)
                       divide 2M digits round-mode)]
          (if (= x0 x1)
            x1tot
            (recur num x1 x1tot))))

      (big-sqrt-int
       num 0M (. BigDecimal valueOf (Math/sqrt (. num doubleValue)))))

    (defn sb-pi-int
      [#^BigDecimal a #^BigDecimal b #^BigDecimal
       t #^BigDecimal x #^BigDecimal y]

      (let
          [#^BigDecimal y1 a
           #^BigDecimal a1 (. (+ a b) divide 2M digits round-mode)
           #^BigDecimal b1 (big-sqrt (* b y1))
           #^BigDecimal t1 (- t
                              (* x
                                 (* (- y1 a1)
                                    (- y1 a1))))
           #^BigDecimal x1 (* x 2M)]

        (if (. a equals b)
          (. (. (* (+ a1 b1) (+ a1 b1))
             divide (* t1 4M) digits round-mode)
             setScale places round-mode)
          (recur a1 b1 t1 x1 y1))))

      (sb-pi-int 1M
                 (. 1M divide (big-sqrt 2M) digits round-mode)
                 (/ 1M 4M)
                 1M
                 nil)))

      ;(big-sqrt (new BigDecimal 2))))

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


$ time clj pi.clj 1               --> 0.977s
$ time clj pi.clj 10              --> 0.975s
$ time clj pi.clj 100             --> 0.979s
$ time clj pi.clj 1000            --> 1.089s
$ time clj pi.clj 10000           --> 11.759s




The same algorithm in Java is faster by a rough factor of 3x :



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) {
    System.out.println(pi(Integer.parseInt(args[0])));
  }

  // 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)
      .multiply(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         ---->      0.367s
$ time java Pi 10        ---->      0.366s
$ time java Pi 100       ---->      0.370s
$ time java Pi 1000      ---->      0.518s
$ time java Pi 10000     ---->      3.425s

--~--~---------~--~----~------------~-------~--~----~
You received this message because you are subscribed to the Google
Groups "Clojure" group.
To post to this group, send email to clojure@googlegroups.com
Note that posts from new members are moderated - please be patient with your 
first post.
To unsubscribe from this group, send email to
clojure+unsubscr...@googlegroups.com
For more options, visit this group at
http://groups.google.com/group/clojure?hl=en
-~----------~----~----~----~------~----~------~--~---