Вычисление 1000000 знаков числа Пи. На iPhone

import sys
import math
from decimal import *
from math import factorial
from time import time

# http://thelivingpearl.com/2013/05/28/computing-pi-with-python/
# http://www.numberworld.org/misc_runs/pi-5t/details.html

# Bellard's Formula PI

def bellard(n):
    pi = Decimal(0)
    k = 0
    while k < n:
      pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) - Decimal(1)/(4*k+3))
      k += 1
    pi = pi * 1/(2**6)
    return pi

# Bailey-Borwein-Plouffe formula

def bbp(n):
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6)))
        k += 1
    return pi

# http://www.craig-wood.com/nick/articles/pi-chudnovsky/

def chudnovsky(n):
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
        k += 1
    pi = pi * Decimal(10005).sqrt()/4270934400
    pi = pi**(-1)
    return pi

def chudnovsky2(n):
    pi = Decimal(13591409)
    ak = Decimal(1)
    k = 1
    while k < n:
        ak *= -Decimal((6*k-5)*(2*k-1)*(6*k-1))/Decimal(k*k*k*26680*640320*640320)

        val = ak*(13591409 + 545140134*k)
        
        d = Decimal((6*k-5)*(2*k-1)*(6*k-1))/Decimal(k*k*k*26680*640320*640320)
        
        pi += val
        k += 1
    pi = pi * Decimal(10005).sqrt()/4270934400
    pi = pi**(-1)
    return pi

# arctan
# http://www.craig-wood.com/nick/articles/pi-machin/

def arctan(x):
    """
    Calculate arctan(1/x)

    arctan(1/x) = 1/x - 1/3*x**3 + 1/5*x**5 - ... (x >= 1)

    This calculates it in fixed point, using the value for one passed in
    """
    x2 = x*x
    x = Decimal(x)

    total = Decimal(0)
    sign = 1
    for i in range(1, 512, 2):
        #total += sign / (i * x)
        total += sign / (i * x ** i)
        sign = -sign
        #x *= x2
        #print(total)
    return total

def pi_machin(n):
    return 4*(4*arctan(5) - arctan(239))

def pi_gauss(one):
    return 4*(12*arctan(18) + 8*arctan(57) - 5*arctan(239))

if __name__ == "__main__":
  
  N = 1000
  
  getcontext().prec = N
  
  print ""
  print "Atan"
  start = time()
  my_pi = pi_machin(N)
  print "Pi = {}".format(str(my_pi))
  print("Time", time()-start)
  
  
  print "BBP"
  start = time()
  my_pi = bbp(N)
  print "Pi = {}".format(str(my_pi))
  print("Time", time()-start)
  
  print "Bellard"
  start = time()
  my_pi = bellard(N)
  print "Pi = {}".format(str(my_pi))
  print("Time", time()-start)
  
  print ""
  print "Chudnovsky"
  start = time()
  my_pi = chudnovsky(N/14)
  print "Pi = {}".format(str(my_pi))
  print("Time", time()-start)
  
  print ""
  print "Chudnovsky2"
  start = time()
  my_pi = chudnovsky2(N/14)
  print "Pi = {}".format(str(my_pi))
  print("Time", time()-start)

#include "gmp.h"
#include "gmpxx.h"

- (void) calcPi {
  /*
   # Python source
   pi = Decimal(13591409)
   ak = Decimal(1)
   k = 1
   while k < n:
      ak *= -Decimal((6*k-5)*(2*k-1)*(6*k-1))/Decimal(k*k*k*26680*640320*640320)
   
      val = ak*(13591409 + 545140134*k)
      pi += val
      k += 1
   #print "Iteration: {} of {}".format(k, n)
   pi = pi * Decimal(10005).sqrt()/4270934400
   pi = pi**(-1)
   return pi
   */
  
  unsigned long int bits = [self getBitsSize];
  mpf_set_default_prec(bits);
  
  NSDate *date1 = [NSDate date];
  
  mpf_class pi = 13591409;
  mpf_class ak = 1;
  mpf_class v1, v2, v3, tmp1, d = 0, d1 = 0;
  
  unsigned long int N = (unsigned long int)self.numDigits/14;
  for(unsigned long int k=1; k<N; k++) {
    // (6*k-5)*(2*k-1)*(6*k-1)
    v1 = 6*k-5;
    mpf_mul_ui(v2.get_mpf_t(), v1.get_mpf_t(), 2*k-1);
    mpf_mul_ui(v1.get_mpf_t(), v2.get_mpf_t(), 6*k-1); // v1 = (6*k-5)*(2*k-1)*(6*k-1)

    // k*k*k*26680*640320*640320
    v2 = k;
    mpf_mul_ui(v3.get_mpf_t(), v2.get_mpf_t(), k);
    mpf_mul_ui(v2.get_mpf_t(), v3.get_mpf_t(), k);
    mpf_mul_ui(tmp1.get_mpf_t(), v2.get_mpf_t(), 26680); // tmp <= k*26680
    mpf_mul_ui(v2.get_mpf_t(), tmp1.get_mpf_t(), 640320);
    mpf_mul_ui(tmp1.get_mpf_t(), v2.get_mpf_t(), 640320);
    
    // v2 <= v1/tmp = (6*k-5)*(2*k-1)*(6*k-1)/(k*k*k*26680*640320*640320)
    mpf_div(v2.get_mpf_t(), v1.get_mpf_t(), tmp1.get_mpf_t());
    mpf_neg(v1.get_mpf_t(), v2.get_mpf_t()); // v1 = -v1
    
    mpf_mul(tmp1.get_mpf_t(), ak.get_mpf_t(), v1.get_mpf_t()); // tmp <= ak*v1
    mpf_set(ak.get_mpf_t(), tmp1.get_mpf_t()); // ak = ak*v1
    
    v1 = 545140134;
    mpf_mul_ui(v2.get_mpf_t(), v1.get_mpf_t(), k); // v2 = 545140134*k
    mpf_add_ui(v1.get_mpf_t(), v2.get_mpf_t(), 13591409); // v1 = 545140134*k + 13591409
    mpf_mul(tmp1.get_mpf_t(), v1.get_mpf_t(), ak.get_mpf_t()); // tmp1 = v1*ak
        
    mpf_add(v1.get_mpf_t(), tmp1.get_mpf_t(), pi.get_mpf_t()); // v1 = tmp1 + pi
    mpf_set(pi.get_mpf_t(), v1.get_mpf_t());  // pi = v1
    
    if (k % 5 == 0) {
      // Release CPU when app is inactive, otherwise app will be killed by iOS
      if ([[UIApplication sharedApplication] applicationState] != UIApplicationStateActive)
        [NSThread sleepForTimeInterval:3];
    }
  }
  
  mpf_sqrt_ui(d1.get_mpf_t(), 10005); // d1 = sqrt(10005)
  mpf_mul(d.get_mpf_t(), d1.get_mpf_t(), pi.get_mpf_t()); // d = d1*pi
  mpf_div_ui(d1.get_mpf_t(), d.get_mpf_t(), 4270934400); // d1 = d/4270934400
  
  mpf_ui_div(d.get_mpf_t(), 1, d1.get_mpf_t()); // d = 1/d1

  double diff = [[NSDate date] timeIntervalSinceDate:date1];
  //gmp_printf("pi: %.*Ffn", self.numDigits, d.get_mpf_t()); // Not works for big numbers, only for debugging
  [self valueSave:d.get_mpf_t() time:diff];
  
  NSLog(@"Done-pi, time = %.2f", diff);
  
  self.totalTime = diff;
}

- (void)valueSave:(mpf_t)val time:(double)time{
  NSString *docsDirectory = [NSSearchPathForDirectoriesInDomains(NSDocumentDirectory, NSUserDomainMask, YES) objectAtIndex:0];
  NSString *name = [NSString stringWithFormat:@"value-%llu.txt", self.numDigits];
  NSString *path = [docsDirectory stringByAppendingPathComponent:name];
  
  FILE *fp = fopen ([path UTF8String], "w+");
  char buf[255];
  sprintf(buf, "Time: %fsn", time);
  fwrite(buf , sizeof(char), strlen(buf), fp);
  mpf_out_str(fp, 10, 0, val); // 10-основание системы счисления
  fclose(fp);
  
  NSLog(@"Saved in: %@", path);
}

- (unsigned long int)getBitsSize {
  // Calculate the number of bits, required to store numDigits (empiric way)
  return (unsigned long int)(2.5*self.numDigits/log(2));
}