Baking-Pi Challenge - Understanding & Improving

2019-04-12 15:38发布

问题:

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.

回答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 favourite PSLQ up to 800 digits of Pi

    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);}