1 /**
2  * Licensed to the Apache Software Foundation (ASF) under one
3  * or more contributor license agreements.  See the NOTICE file
4  * distributed with this work for additional information
5  * regarding copyright ownership.  The ASF licenses this file
6  * to you under the Apache License, Version 2.0 (the
7  * "License"); you may not use this file except in compliance
8  * with the License.  You may obtain a copy of the License at
9  *
10  *     http://www.apache.org/licenses/LICENSE-2.0
11  *
12  * Unless required by applicable law or agreed to in writing, software
13  * distributed under the License is distributed on an "AS IS" BASIS,
14  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
15  * See the License for the specific language governing permissions and
16  * limitations under the License.
17  */
18 package org.apache.hadoop.examples.pi.math;
19 
20 import java.io.IOException;
21 import java.util.ArrayList;
22 import java.util.Arrays;
23 import java.util.Collections;
24 import java.util.Iterator;
25 import java.util.List;
26 import java.util.Map;
27 import java.util.TreeMap;
28 import java.util.NoSuchElementException;
29 
30 import org.apache.hadoop.examples.pi.Container;
31 import org.apache.hadoop.examples.pi.Util;
32 
33 /**
34  * Bellard's BBP-type Pi formula
35  * 1/2^6 \sum_{n=0}^\infty (-1)^n/2^{10n}
36  * (-2^5/(4n+1) -1/(4n+3) +2^8/(10n+1) -2^6/(10n+3) -2^2/(10n+5)
37  *  -2^2/(10n+7) +1/(10n+9))
38  *
39  * References:
40  *
41  * [1] David H. Bailey, Peter B. Borwein and Simon Plouffe.  On the Rapid
42  *     Computation of Various Polylogarithmic Constants.
43  *     Math. Comp., 66:903-913, 1996.
44  *
45  * [2] Fabrice Bellard.  A new formula to compute the n'th binary digit of pi,
46  *     1997.  Available at http://fabrice.bellard.free.fr/pi .
47  */
48 public final class Bellard {
49   /** Parameters for the sums */
50   public enum Parameter {
51     // \sum_{k=0}^\infty (-1)^{k+1}( 2^{d-10k-1}/(4k+1) + 2^{d-10k-6}/(4k+3) )
52     P8_1(false, 1, 8, -1),
53     P8_3(false, 3, 8, -6),
54     P8_5(P8_1),
55     P8_7(P8_3),
56 
57     /*
58      *   2^d\sum_{k=0}^\infty (-1)^k( 2^{ 2-10k} / (10k + 1)
59      *                               -2^{  -10k} / (10k + 3)
60      *                               -2^{-4-10k} / (10k + 5)
61      *                               -2^{-4-10k} / (10k + 7)
62      *                               +2^{-6-10k} / (10k + 9) )
63      */
64     P20_21(true , 1, 20,  2),
65     P20_3(false, 3, 20,  0),
66     P20_5(false, 5, 20, -4),
67     P20_7(false, 7, 20, -4),
68     P20_9(true , 9, 20, -6),
69     P20_11(P20_21),
70     P20_13(P20_3),
71     P20_15(P20_5),
72     P20_17(P20_7),
73     P20_19(P20_9);
74 
75     final boolean isplus;
76     final long j;
77     final int deltaN;
78     final int deltaE;
79     final int offsetE;
80 
Parameter(boolean isplus, long j, int deltaN, int offsetE)81     private Parameter(boolean isplus, long j, int deltaN, int offsetE) {
82       this.isplus = isplus;
83       this.j = j;
84       this.deltaN = deltaN;
85       this.deltaE = -20;
86       this.offsetE = offsetE;
87     }
88 
Parameter(Parameter p)89     private Parameter(Parameter p) {
90       this.isplus = !p.isplus;
91       this.j = p.j + (p.deltaN >> 1);
92       this.deltaN = p.deltaN;
93       this.deltaE = p.deltaE;
94       this.offsetE = p.offsetE + (p.deltaE >> 1);
95     }
96 
97     /** Get the Parameter represented by the String */
get(String s)98     public static Parameter get(String s) {
99       s = s.trim();
100       if (s.charAt(0) == 'P')
101         s = s.substring(1);
102       final String[] parts = s.split("\\D+");
103       if (parts.length >= 2) {
104         final String name = "P" + parts[0] + "_" + parts[1];
105         for(Parameter p : values())
106           if (p.name().equals(name))
107             return p;
108       }
109       throw new IllegalArgumentException("s=" + s
110           + ", parts=" + Arrays.asList(parts));
111     }
112   }
113 
114   /** The sums in the Bellard's formula */
115   public static class Sum implements Container<Summation>, Iterable<Summation> {
116     private static final long ACCURACY_BIT = 50;
117 
118     private final Parameter parameter;
119     private final Summation sigma;
120     private final Summation[] parts;
121     private final Tail tail;
122 
123     /** Constructor */
Sum(long b, Parameter p, int nParts, List<T> existing)124     private <T extends Container<Summation>> Sum(long b, Parameter p, int nParts, List<T> existing) {
125       if (b < 0)
126         throw new IllegalArgumentException("b = " + b + " < 0");
127       if (nParts < 1)
128         throw new IllegalArgumentException("nParts = " + nParts + " < 1");
129       final long i = p.j == 1 && p.offsetE >= 0? 1 : 0;
130       final long e = b + i*p.deltaE + p.offsetE;
131       final long n = i*p.deltaN + p.j;
132 
133       this.parameter = p;
134       this.sigma = new Summation(n, p.deltaN, e, p.deltaE, 0);
135       this.parts = partition(sigma, nParts, existing);
136       this.tail = new Tail(n, e);
137     }
138 
partition( Summation sigma, int nParts, List<T> existing)139     private static <T extends Container<Summation>> Summation[] partition(
140         Summation sigma, int nParts, List<T> existing) {
141       final List<Summation> parts = new ArrayList<Summation>();
142       if (existing == null || existing.isEmpty())
143         parts.addAll(Arrays.asList(sigma.partition(nParts)));
144       else {
145         final long stepsPerPart = sigma.getSteps()/nParts;
146         final List<Summation> remaining = sigma.remainingTerms(existing);
147 
148         for(Summation s : remaining) {
149           final int n = (int)((s.getSteps() - 1)/stepsPerPart) + 1;
150           parts.addAll(Arrays.asList(s.partition(n)));
151         }
152 
153         for(Container<Summation> c : existing)
154           parts.add(c.getElement());
155         Collections.sort(parts);
156       }
157       return parts.toArray(new Summation[parts.size()]);
158     }
159 
160     /** {@inheritDoc} */
161     @Override
toString()162     public String toString() {
163       int n = 0;
164       for(Summation s : parts)
165         if (s.getValue() == null)
166           n++;
167       return getClass().getSimpleName() + "{" + parameter + ": " + sigma
168           + ", remaining=" + n + "}";
169     }
170 
171     /** Set the value of sigma */
setValue(Summation s)172     public void setValue(Summation s) {
173       if (s.getValue() == null)
174         throw new IllegalArgumentException("s.getValue()"
175             + "\n  sigma=" + sigma
176             + "\n  s    =" + s);
177       if (!s.contains(sigma) || !sigma.contains(s))
178         throw new IllegalArgumentException("!s.contains(sigma) || !sigma.contains(s)"
179             + "\n  sigma=" + sigma
180             + "\n  s    =" + s);
181       sigma.setValue(s.getValue());
182     }
183 
184     /** get the value of sigma */
getValue()185     public double getValue() {
186       if (sigma.getValue() == null) {
187         double d = 0;
188         for(int i = 0; i < parts.length; i++)
189           d = Modular.addMod(d, parts[i].compute());
190         sigma.setValue(d);
191       }
192 
193       final double s = Modular.addMod(sigma.getValue(), tail.compute());
194       return parameter.isplus? s: -s;
195     }
196 
197     /** {@inheritDoc} */
198     @Override
getElement()199     public Summation getElement() {
200       if (sigma.getValue() == null) {
201         int i = 0;
202         double d = 0;
203         for(; i < parts.length && parts[i].getValue() != null; i++)
204           d = Modular.addMod(d, parts[i].getValue());
205         if (i == parts.length)
206           sigma.setValue(d);
207       }
208       return sigma;
209     }
210 
211     /** The sum tail */
212     private class Tail {
213       private long n;
214       private long e;
215 
Tail(long n, long e)216       private Tail(long n, long e) {
217         this.n = n;
218         this.e = e;
219       }
220 
compute()221       private double compute() {
222         if (e > 0) {
223           final long edelta = -sigma.E.delta;
224           long q = e / edelta;
225           long r = e % edelta;
226           if (r == 0) {
227             e = 0;
228             n += q * sigma.N.delta;
229           } else {
230             e = edelta - r;
231             n += (q + 1)*sigma.N.delta;
232           }
233         } else if (e < 0)
234           e = -e;
235 
236         double s = 0;
237         for(;; e -= sigma.E.delta) {
238           if (e > ACCURACY_BIT || (1L << (ACCURACY_BIT - e)) < n)
239             return s;
240 
241           s += 1.0 / (n << e);
242           if (s >= 1) s--;
243           n += sigma.N.delta;
244         }
245       }
246     }
247 
248     /** {@inheritDoc} */
249     @Override
iterator()250     public Iterator<Summation> iterator() {
251       return new Iterator<Summation>() {
252         private int i = 0;
253 
254         /** {@inheritDoc} */
255         @Override
256         public boolean hasNext() {return i < parts.length;}
257         /** {@inheritDoc} */
258         @Override
259         public Summation next() throws NoSuchElementException {
260           if (hasNext()) {
261             return parts[i++];
262           } else {
263             throw new NoSuchElementException("Sum's iterator does not have next!");
264           }
265         }
266         /** Unsupported */
267         @Override
268         public void remove() {throw new UnsupportedOperationException();}
269       };
270     }
271   }
272 
273   /** Get the sums for the Bellard formula. */
getSums( long b, int partsPerSum, Map<Parameter, List<T>> existing)274   public static <T extends Container<Summation>> Map<Parameter, Sum> getSums(
275       long b, int partsPerSum, Map<Parameter, List<T>> existing) {
276     final Map<Parameter, Sum> sums = new TreeMap<Parameter, Sum>();
277     for(Parameter p : Parameter.values()) {
278       final Sum s = new Sum(b, p, partsPerSum, existing.get(p));
279       Util.out.println("put " + s);
280       sums.put(p, s);
281     }
282     return sums;
283   }
284 
285   /** Compute bits of Pi from the results. */
computePi( final long b, Map<Parameter, T> results)286   public static <T extends Container<Summation>> double computePi(
287       final long b, Map<Parameter, T> results) {
288     if (results.size() != Parameter.values().length)
289       throw new IllegalArgumentException("m.size() != Parameter.values().length"
290           + ", m.size()=" + results.size()
291           + "\n  m=" + results);
292 
293     double pi = 0;
294     for(Parameter p : Parameter.values()) {
295       final Summation sigma = results.get(p).getElement();
296       final Sum s = new Sum(b, p, 1, null);
297       s.setValue(sigma);
298       pi = Modular.addMod(pi, s.getValue());
299     }
300     return pi;
301   }
302 
303   /** Compute bits of Pi in the local machine. */
computePi(final long b)304   public static double computePi(final long b) {
305     double pi = 0;
306     for(Parameter p : Parameter.values())
307       pi = Modular.addMod(pi, new Sum(b, p, 1, null).getValue());
308     return pi;
309   }
310 
311   /** Estimate the number of terms. */
bit2terms(long b)312   public static long bit2terms(long b) {
313     return 7*(b/10);
314   }
315 
computePi(Util.Timer t, long b)316   private static void computePi(Util.Timer t, long b) {
317     t.tick(Util.pi2string(computePi(b), bit2terms(b)));
318   }
319 
320   /** main */
main(String[] args)321   public static void main(String[] args) throws IOException {
322     final Util.Timer t = new Util.Timer(false);
323 
324     computePi(t, 0);
325     computePi(t, 1);
326     computePi(t, 2);
327     computePi(t, 3);
328     computePi(t, 4);
329 
330     Util.printBitSkipped(1008);
331     computePi(t, 1008);
332     computePi(t, 1012);
333 
334     long b = 10;
335     for(int i = 0; i < 7; i++) {
336       Util.printBitSkipped(b);
337       computePi(t, b - 4);
338       computePi(t, b);
339       computePi(t, b + 4);
340       b *= 10;
341     }
342   }
343 }
344