| /** |
| * Licensed to the Apache Software Foundation (ASF) under one |
| * or more contributor license agreements. See the NOTICE file |
| * distributed with this work for additional information |
| * regarding copyright ownership. The ASF licenses this file |
| * to you under the Apache License, Version 2.0 (the |
| * "License"); you may not use this file except in compliance |
| * with the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.hadoop.examples.pi.math; |
| |
| import java.io.IOException; |
| import java.util.ArrayList; |
| import java.util.Arrays; |
| import java.util.Collections; |
| import java.util.Iterator; |
| import java.util.List; |
| import java.util.Map; |
| import java.util.TreeMap; |
| |
| import org.apache.hadoop.examples.pi.Container; |
| import org.apache.hadoop.examples.pi.Util; |
| |
| /** |
| * Bellard's BBP-type Pi formula |
| * 1/2^6 \sum_{n=0}^\infty (-1)^n/2^{10n} |
| * (-2^5/(4n+1) -1/(4n+3) +2^8/(10n+1) -2^6/(10n+3) -2^2/(10n+5) |
| * -2^2/(10n+7) +1/(10n+9)) |
| * |
| * References: |
| * |
| * [1] David H. Bailey, Peter B. Borwein and Simon Plouffe. On the Rapid |
| * Computation of Various Polylogarithmic Constants. |
| * Math. Comp., 66:903-913, 1996. |
| * |
| * [2] Fabrice Bellard. A new formula to compute the n'th binary digit of pi, |
| * 1997. Available at http://fabrice.bellard.free.fr/pi . |
| */ |
| public final class Bellard { |
| /** Parameters for the sums */ |
| public enum Parameter { |
| // \sum_{k=0}^\infty (-1)^{k+1}( 2^{d-10k-1}/(4k+1) + 2^{d-10k-6}/(4k+3) ) |
| P8_1(false, 1, 8, -1), |
| P8_3(false, 3, 8, -6), |
| P8_5(P8_1), |
| P8_7(P8_3), |
| |
| /* |
| * 2^d\sum_{k=0}^\infty (-1)^k( 2^{ 2-10k} / (10k + 1) |
| * -2^{ -10k} / (10k + 3) |
| * -2^{-4-10k} / (10k + 5) |
| * -2^{-4-10k} / (10k + 7) |
| * +2^{-6-10k} / (10k + 9) ) |
| */ |
| P20_21(true , 1, 20, 2), |
| P20_3(false, 3, 20, 0), |
| P20_5(false, 5, 20, -4), |
| P20_7(false, 7, 20, -4), |
| P20_9(true , 9, 20, -6), |
| P20_11(P20_21), |
| P20_13(P20_3), |
| P20_15(P20_5), |
| P20_17(P20_7), |
| P20_19(P20_9); |
| |
| final boolean isplus; |
| final long j; |
| final int deltaN; |
| final int deltaE; |
| final int offsetE; |
| |
| private Parameter(boolean isplus, long j, int deltaN, int offsetE) { |
| this.isplus = isplus; |
| this.j = j; |
| this.deltaN = deltaN; |
| this.deltaE = -20; |
| this.offsetE = offsetE; |
| } |
| |
| private Parameter(Parameter p) { |
| this.isplus = !p.isplus; |
| this.j = p.j + (p.deltaN >> 1); |
| this.deltaN = p.deltaN; |
| this.deltaE = p.deltaE; |
| this.offsetE = p.offsetE + (p.deltaE >> 1); |
| } |
| |
| /** Get the Parameter represented by the String */ |
| public static Parameter get(String s) { |
| s = s.trim(); |
| if (s.charAt(0) == 'P') |
| s = s.substring(1); |
| final String[] parts = s.split("\\D+"); |
| if (parts.length >= 2) { |
| final String name = "P" + parts[0] + "_" + parts[1]; |
| for(Parameter p : values()) |
| if (p.name().equals(name)) |
| return p; |
| } |
| throw new IllegalArgumentException("s=" + s |
| + ", parts=" + Arrays.asList(parts)); |
| } |
| } |
| |
| /** The sums in the Bellard's formula */ |
| public static class Sum implements Container<Summation>, Iterable<Summation> { |
| private static final long ACCURACY_BIT = 50; |
| |
| private final Parameter parameter; |
| private final Summation sigma; |
| private final Summation[] parts; |
| private final Tail tail; |
| |
| /** Constructor */ |
| private <T extends Container<Summation>> Sum(long b, Parameter p, int nParts, List<T> existing) { |
| if (b < 0) |
| throw new IllegalArgumentException("b = " + b + " < 0"); |
| if (nParts < 1) |
| throw new IllegalArgumentException("nParts = " + nParts + " < 1"); |
| final long i = p.j == 1 && p.offsetE >= 0? 1 : 0; |
| final long e = b + i*p.deltaE + p.offsetE; |
| final long n = i*p.deltaN + p.j; |
| |
| this.parameter = p; |
| this.sigma = new Summation(n, p.deltaN, e, p.deltaE, 0); |
| this.parts = partition(sigma, nParts, existing); |
| this.tail = new Tail(n, e); |
| } |
| |
| private static <T extends Container<Summation>> Summation[] partition( |
| Summation sigma, int nParts, List<T> existing) { |
| final List<Summation> parts = new ArrayList<Summation>(); |
| if (existing == null || existing.isEmpty()) |
| parts.addAll(Arrays.asList(sigma.partition(nParts))); |
| else { |
| final long stepsPerPart = sigma.getSteps()/nParts; |
| final List<Summation> remaining = sigma.remainingTerms(existing); |
| |
| for(Summation s : remaining) { |
| final int n = (int)((s.getSteps() - 1)/stepsPerPart) + 1; |
| parts.addAll(Arrays.asList(s.partition(n))); |
| } |
| |
| for(Container<Summation> c : existing) |
| parts.add(c.getElement()); |
| Collections.sort(parts); |
| } |
| return parts.toArray(new Summation[parts.size()]); |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public String toString() { |
| int n = 0; |
| for(Summation s : parts) |
| if (s.getValue() == null) |
| n++; |
| return getClass().getSimpleName() + "{" + parameter + ": " + sigma |
| + ", remaining=" + n + "}"; |
| } |
| |
| /** Set the value of sigma */ |
| public void setValue(Summation s) { |
| if (s.getValue() == null) |
| throw new IllegalArgumentException("s.getValue()" |
| + "\n sigma=" + sigma |
| + "\n s =" + s); |
| if (!s.contains(sigma) || !sigma.contains(s)) |
| throw new IllegalArgumentException("!s.contains(sigma) || !sigma.contains(s)" |
| + "\n sigma=" + sigma |
| + "\n s =" + s); |
| sigma.setValue(s.getValue()); |
| } |
| |
| /** get the value of sigma */ |
| public double getValue() { |
| if (sigma.getValue() == null) { |
| double d = 0; |
| for(int i = 0; i < parts.length; i++) |
| d = Modular.addMod(d, parts[i].compute()); |
| sigma.setValue(d); |
| } |
| |
| final double s = Modular.addMod(sigma.getValue(), tail.compute()); |
| return parameter.isplus? s: -s; |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public Summation getElement() { |
| if (sigma.getValue() == null) { |
| int i = 0; |
| double d = 0; |
| for(; i < parts.length && parts[i].getValue() != null; i++) |
| d = Modular.addMod(d, parts[i].getValue()); |
| if (i == parts.length) |
| sigma.setValue(d); |
| } |
| return sigma; |
| } |
| |
| /** The sum tail */ |
| private class Tail { |
| private long n; |
| private long e; |
| |
| private Tail(long n, long e) { |
| this.n = n; |
| this.e = e; |
| } |
| |
| private double compute() { |
| if (e > 0) { |
| final long edelta = -sigma.E.delta; |
| long q = e / edelta; |
| long r = e % edelta; |
| if (r == 0) { |
| e = 0; |
| n += q * sigma.N.delta; |
| } else { |
| e = edelta - r; |
| n += (q + 1)*sigma.N.delta; |
| } |
| } else if (e < 0) |
| e = -e; |
| |
| double s = 0; |
| for(;; e -= sigma.E.delta) { |
| if (e > ACCURACY_BIT || (1L << (ACCURACY_BIT - e)) < n) |
| return s; |
| |
| s += 1.0 / (n << e); |
| if (s >= 1) s--; |
| n += sigma.N.delta; |
| } |
| } |
| } |
| |
| /** {@inheritDoc} */ |
| @Override |
| public Iterator<Summation> iterator() { |
| return new Iterator<Summation>() { |
| private int i = 0; |
| |
| /** {@inheritDoc} */ |
| @Override |
| public boolean hasNext() {return i < parts.length;} |
| /** {@inheritDoc} */ |
| @Override |
| public Summation next() {return parts[i++];} |
| /** Unsupported */ |
| @Override |
| public void remove() {throw new UnsupportedOperationException();} |
| }; |
| } |
| } |
| |
| /** Get the sums for the Bellard formula. */ |
| public static <T extends Container<Summation>> Map<Parameter, Sum> getSums( |
| long b, int partsPerSum, Map<Parameter, List<T>> existing) { |
| final Map<Parameter, Sum> sums = new TreeMap<Parameter, Sum>(); |
| for(Parameter p : Parameter.values()) { |
| final Sum s = new Sum(b, p, partsPerSum, existing.get(p)); |
| Util.out.println("put " + s); |
| sums.put(p, s); |
| } |
| return sums; |
| } |
| |
| /** Compute bits of Pi from the results. */ |
| public static <T extends Container<Summation>> double computePi( |
| final long b, Map<Parameter, T> results) { |
| if (results.size() != Parameter.values().length) |
| throw new IllegalArgumentException("m.size() != Parameter.values().length" |
| + ", m.size()=" + results.size() |
| + "\n m=" + results); |
| |
| double pi = 0; |
| for(Parameter p : Parameter.values()) { |
| final Summation sigma = results.get(p).getElement(); |
| final Sum s = new Sum(b, p, 1, null); |
| s.setValue(sigma); |
| pi = Modular.addMod(pi, s.getValue()); |
| } |
| return pi; |
| } |
| |
| /** Compute bits of Pi in the local machine. */ |
| public static double computePi(final long b) { |
| double pi = 0; |
| for(Parameter p : Parameter.values()) |
| pi = Modular.addMod(pi, new Sum(b, p, 1, null).getValue()); |
| return pi; |
| } |
| |
| /** Estimate the number of terms. */ |
| public static long bit2terms(long b) { |
| return 7*(b/10); |
| } |
| |
| private static void computePi(Util.Timer t, long b) { |
| t.tick(Util.pi2string(computePi(b), bit2terms(b))); |
| } |
| |
| /** main */ |
| public static void main(String[] args) throws IOException { |
| final Util.Timer t = new Util.Timer(false); |
| |
| computePi(t, 0); |
| computePi(t, 1); |
| computePi(t, 2); |
| computePi(t, 3); |
| computePi(t, 4); |
| |
| Util.printBitSkipped(1008); |
| computePi(t, 1008); |
| computePi(t, 1012); |
| |
| long b = 10; |
| for(int i = 0; i < 7; i++) { |
| Util.printBitSkipped(b); |
| computePi(t, b - 4); |
| computePi(t, b); |
| computePi(t, b + 4); |
| b *= 10; |
| } |
| } |
| } |