Record que arrasa el anterior sobre cálculo de cifras del número pi

Hace unos mese publiqué un post sobre cómo un Fabrice Bellard con un PC "normalito" había batido el record de computación de cifras del número pi, elevando dicho record hasta 2,7 billones de cifras.
Hoy, con grata sorpresa leo en microsiervos que Nicholas Sze acaba de fulminar ese record y otro que se produjo hace sólo unas semanas, calculando hasta 2.000 billones de cifras.
Personalmente esto no tiene tanto mérito como la hazaña de Bellard, ya que Sze se ha valido de 1000 ordenadores de la granja de Yahoo durante 23 días, para alcanzar dicho record.
Mientras que Bellard utilizó el método de Chudnovsky para aproximar pi, Sze ha utilizado otro método adaptado para utilizar la computación paralela basada en el framework MapReduce (creado por Google).

Personalmente, hace meses cuando leí la noticia de Bellard intenté hacerme mi propia implementación en Java del algoritmo para verificar en mis propias carnes su funcionamiento, y la verdad es que conseguí aproximar 100.000 cifras en un tiempo relativamente pequeño.

Probé varias aproximaciones, una de ellas basada en la clase BigSquareRoot que encontré aquí. La implementación del código en Java, por si os interesa. La clase CalculaPi:


import java.math.BigDecimal;
import java.math.BigInteger;
import java.math.RoundingMode;


public class CalculaPi {


private static final BigDecimal DA12 = BigDecimal.valueOf(13591409l * 12l);


private static BigDecimal DC32;

private static final BigDecimal DDOCE = BigDecimal.valueOf(12);


private static final BigInteger IA = BigInteger.valueOf(13591409l);
private static final BigInteger IB = BigInteger.valueOf(545140134l);
private static final BigInteger IC3P24 = BigInteger.valueOf(10939058860032000l);
private static final BigInteger IDOS = BigInteger.valueOf(2);
private static final BigInteger IMENOSUNO = BigInteger.valueOf(-1);
private static final BigInteger IUNO = BigInteger.valueOf(1);
private static final BigInteger ICERO = BigInteger.valueOf(0);
private static final BigInteger ISEIS = BigInteger.valueOf(6);
private static final BigInteger IMENOSCINCO = BigInteger.valueOf(-5);


private static final int ESCALA = 100000;

public static void main(String[] args) {
BigInteger n;
BigSquareRoot app = new BigSquareRoot ();
n = BigInteger.valueOf(640320l);
app.setScale (ESCALA);
DC32 = app.get(n).multiply(BigDecimal.valueOf(640320l));

System.out.println("Raiz: " + DC32);
for(int i=999; i<1003; i++) {
BigDecimal tmp = paraN(BigInteger.valueOf(i));
System.out.println("Para " + i + " es: " + tmp.toPlainString());
}
}

public static BigDecimal paraN(BigInteger n) {
PQT temporal = calcular(ICERO, n);
return DC32.multiply(temporal.q).divide(DDOCE.multiply(temporal.t).
add(DA12.multiply(temporal.q)), ESCALA, RoundingMode.CEILING);
}


public static PQT calcular(BigInteger n1, BigInteger n2) {
if(n1.add(IUNO).equals(n2)) {
return new PQT( p(n2), q(n2), a(n2).multiply(p(n2)));
}
BigInteger m = n1.add(n2).divide(IDOS);
PQT uno = calcular(n1, m);
PQT dos = calcular(m, n2);
return new PQT(uno.p.multiply(dos.p), 
  uno.q.multiply(dos.q), 
  uno.t.multiply(dos.q).add(uno.p.multiply(dos.t)));
}

public static BigInteger a(BigInteger n) {
BigInteger ret = IA.add(IB.multiply(n));
if(n.remainder(IDOS).equals(IUNO)) { // Optimizar con AND, ver el último bit y se sabe si divisible por 2
ret = ret.negate();
}
return ret;
}

public static BigInteger q(BigInteger n) {
return n.multiply(n).multiply(n).multiply(IC3P24);
}

public static BigInteger p(BigInteger n) {
return (IDOS.multiply(n).add(IMENOSUNO)).multiply(n.multiply(ISEIS).add(IMENOSCINCO)).multiply(n.multiply(ISEIS).add(IMENOSUNO));
}
}


La clase PQT:


import java.math.BigDecimal;
import java.math.BigInteger;


public class PQT {
public BigDecimal p,q,t;
public PQT(BigDecimal p, BigDecimal q, BigDecimal t) {
super();
this.p = p;
this.q = q;
this.t = t;
}
public PQT(BigInteger p, BigInteger q, BigInteger t) {
super();
this.p = new BigDecimal(p);
this.q = new BigDecimal(q);
this.t = new BigDecimal(t);
}
}


Enlace a la noticia en la BBC.