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