src/examples/org/apache/hadoop/examples/pi/math/Bellard.java

/**
 * 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;
    }
  }
}