Baking-Pi Challenge - Understanding & Improving

Question

I spent some time yesterday writing the solution for this challenge published on Reddit, and was able to get through it without cheating, but I was left with a couple of questions. Reference material here.

This is my code.

(ns baking-pi.core
  (:import java.math.MathContext))

(defn modpow [n e m]
  (.modPow (biginteger n) (biginteger e) (biginteger m)))

(defn div [top bot]
  (with-precision 34 :rounding HALF_EVEN 
    (/ (bigdec top) (bigdec bot))))

(defn pow [n e]
  (.pow (bigdec n) (bigdec e) MathContext/DECIMAL128))

(defn round
  ([n] (.round (bigdec n) MathContext/DECIMAL128))
  ([n & args] (->> [n args] (flatten) (map round))))

(defn left [n d]
  (letfn [(calc [k] (let [bot (+' (*' 8 k) d)
                          top (modpow 16 (-' n k) bot)]
                      (div top bot)))]
    (->> (inc' n)
         (range 0)
         (map calc)
         (reduce +'))))

(defn right [n d]
  (letfn [(calc [[sum'' sum' k]]
                (let [sum' (if (nil? sum') 0M sum')
                      top (pow 16 (-' n k))
                      bot (+' (*' k 8) d)
                      delta (div top bot)]
                  [sum' (+' sum' delta) (inc' k)]))
          (pred [[sum'' sum' k]]
                (cond (or (nil? sum'') (nil? sum')) true
                      (apply == (round sum'' sum')) false
                      :else true))]
    (->> [nil nil (inc' n)]
         (iterate calc)
         (drop-while pred)
         (first)
         (second))))

(defn bbp [n]
  (letfn [(part [m d] (*' m (+' (left n d) (right n d))))]
    (let [sum (-' (part 4 1) (part 2 4) (part 1 5) (part 1 6))]
      (-> sum
          (-' (long sum))
          (*' 16)
          (mod 16)
          (Long/toHexString)))))

I have 2 questions.

1. The wiki makes the following statement. Since my calculation is accurate up to 34 digits after the decimal, how can I leverage it to produce more hexadecimal digits of PI per bbp call?

   in theory, the next few digits up to the accuracy of the calculations used would also be accurate

2. My algorithm relied on BigInteger's modPow for modular exponentiation (based on the following quote), and BigDecimals everywhere else. It is also slow. Bearing in mind that I don't want to lose meaningful accuracy per question #1, what is the best way to speed this program up and make it valid clojurescript as well as clojure?

   To calculate 16 n − k mod (8k + 1) quickly and efficiently, use the modular exponentiation algorithm.

EDIT: Changed from 3 questions to 2. Managed to answer first question on my own.

Matutinal

Answer 1

1. if you want more bits computed per bpp call

   then you have to change your equation from 1/(16^k) base to bigger one. You can do it by summing 2 iterations (k and k+1) so you have something like

   (...)/16^k  + (...)/16^(k+1)
   (...)/256^k
   

   but in this case you need more precise int operations. It is usually faster to use the less precise iterations

2. if you look at the basic equation then you see you do not need bigint for computation at all

   that is why this iterations are used but the output number is bigint of course. So you do not need to compute modular arithmetics on bigint.

   I do not know how optimized are the one you used ... but here are mine:

   • Modular arithmetics and finite field NTT optimizations

3. if you want just speed and not infinite precision then use other PSLQ equations

   My understanding of PSLQ is that it is algorithm to find relation between real number and integer iterations.

here is my favorite up to 800 digits of Pi algorithm and here is extracted code from it in case the link broke down:

//The following 160 character C program, written by Dik T. Winter at CWI, computes pi  to 800 decimal digits. 
int a=10000,b,c=2800,d,e,f[2801],g;main(){for(;b-c;)f[b++]=a/5;
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);}

Goda

Summing up two operations results in a common denominator of 16^(k+1)=16*16^k, not 32^k=(2^k)*(16^k). – Stephanotis

Power binds stronger than multiplication, 16*16^k is 16*(16^k) and not (16*16)^k. Of course, if one sums up the terms from k or k+1 to 2k, then the common denominator is 16^(2k)=256^k. – Stephanotis

16^(2*k) = 256^k and 16^(2*k+1) = 16*(256^k) though. So leveraging the odd/evenness of consecutive terms, you can end up with an expression with a 256^n denominator. – Housemaid

AFAIK, the C code you've attached is not a PSLQ algorithm at all. The parent article was talking about different methods, and the short C program was attached at the end simply because it's a cool program, not because it's PSLQ. I can be sure since I've just reimplementated it. The code used the old-fashioned arctan() series in 32-bit fixed-point arithmetic and calculates 4 decimal digits per iteration, see crypto.stanford.edu/pbc/notes/pi/code.html. – Alpine