#include "gmp.h" #include "gmpxx.h" mpf_set_default_prec(4096); mpf_class num; mpf_sqrt_ui(num.get_mpf_t(), 3); gmp_printf("%.*Ff\n", 1000, num.get_mpf_t());
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)
Formula Machina: 2.01c
Formula Bailey - Borwein - Plaffa: 1.42c
Bellard Formula: 1.82c
Formula Chudnovsky: 0.22c
Formula Chudnovsky (without factorial): 0.082c
#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: %.*Ff\n", 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: %fs\n", 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)); }
Source: https://habr.com/ru/post/309674/
All Articles