📜 ⬆️ ⬇️

Calculating 1,000,000 characters of Pi. On iphone

Mathematics is a very interesting and entertaining science, in which the calculation of the number Pi takes a special place. The history of its computation takes more than 2 thousand years, and the accuracy of computations varies from 256/81 in ancient Egypt and 339/108 in the Vedas, to Jamshid al-Kashi, who calculated 16 characters in the 15th century. What is at least the story of William Shanks, who spent 20 years calculating 700 characters of Pi, but then it turned out that in the second part of the calculations he was wrong ... But the text is not about that, but about algorithms. It became interesting, is it possible to calculate Pi on the iPhone? And if so, with what accuracy?


You can say right away - a million is not the limit. It is possible and more. Details and implementation under the cut.

Implementation of calculations


Even an animal of the Erinaceidae family (i.e. a hedgehog) makes it clear that if we want to calculate Pi to a million signs - such as float is not enough for us. And even double too (where was the sarcasm tag?). So you need to look for a library to work with long numbers. On iOS, Objective C is used (and Swift too, but in this case we don’t need it), which is backward compatible with C, which makes the task somewhat easier - there are quite a few different C libraries. A small search led to the GMP library, which on the one hand, does what we need, on the other hand, is very easy to use.

For example, to display 1000 sqrt (3) characters, the following code is sufficient:
')
#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()); 

A complete list of functions for working with float can be found here .

Algorithm


With the library figured out. The next question is which algorithm to use. The question is not completely idle, because There are many ways to calculate Pi. Formula Bellard:


Formula Bailey - Borway - Plaffa:


Formula Machin:


Formula Chudnovsky (do not ask me how he derived it - I do not know):


The last formula looks like the most “terrible”, but at the same time, the fastest.

Code


To test the algorithms, they were all assembled into one file on Python (it also supports “long numbers” in the Decimal type).

Source code under the spoiler
 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) 


Script startup results:
Formula Machina: 2.01c
Formula Bailey - Borwein - Plaffa: 1.42c
Bellard Formula: 1.82c
Formula Chudnovsky: 0.22c
Formula Chudnovsky (without factorial): 0.082c

As you can see, even with the use of factorial "in the forehead" (which is quite long), Chudnovsky's formula is 5-10 times faster than the previous ones. Converting it to a form without factorial (using the previous value) accelerates several times more. In general, the winner is obvious, but there is one "but" - the amount of RAM. If you “go on the record,” and for example consider a billion characters of Pi, then the question of storing data in RAM becomes critical. The iPhone has only 1GB of memory, and 512MB is available for the user program. If you request approximately 600MB, iOS simply closes the application. When counting a million characters, this is not so critical (the program takes up no more than 30Mb in memory), but if you take more, it will become critical. And here the formula of Chudnovsky will be in a noticeable loss due to the complexity of the operations performed, simpler algorithms may well be useful here.

Speaking of a billion signs. As shown by the launch in the simulator, when trying to create a number with a billion characters, the program tries to allocate more than 5GB of memory, and at the same time, of course, it drops. By estimation, the maximum number, which in principle can be calculated on the iPhone, is about 200 million. marks. On Android, you can get more (who would doubt :). Of course, the calculation will take more than one day, but that is another question.

Below is the source code in Objective-C. Since The work with the GMP library is somewhat cumbersome, the code is not very beautiful, but you can understand the principle (and you can write in any language like in Fortran :).

Code under the spoiler
 #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)); } 


For those who want to experiment on their own, here is the source code of the project for Xcode. Of course, the above algorithm is not perfect, for example, it is quite possible to parallelize at least 2 streams. Or even try using the Accelerate Framework .

results


When you run on the iPhone got out some funny problems. The first is that the program did not want to work in the background, the process of calculations just stopped. On the test iPhone 5S, there was iOS8, I don’t know if the multitasking principle has changed in iOS9, but for 8k, the solution turned out to be simple - start GPS and subscribe to events from the location service. At the same time, iOS “thinks” that we need to constantly receive user coordinates, and allows the program to work even in the minimized state.

The second problem came out a minute later after the first. As it turned out, iOS is quite “smart”, and it considers a suspicious process, eating in the background of 100% processor resources. Apparently, iOS believes that the program has “hung”, and “kills” it in about a minute of such work. I had to insert sleep into the loop so that the “duty cycle” of the CPU load was about 50/50. After that, the program could count in the background, although of course the speed was twice as low. As a result, the iPhone was simply connected to the charger and left for the night, the program did not close.

And finally, the promised result: on the iPhone 5S, a million characters of Pi were calculated to be 49,296 seconds, or about 14 hours. By the way, on the simulator with Core i7 the calculation time is 4200 seconds, i.e. about 11 times faster. It is interesting to estimate how much will be calculated 100 million characters. I leave it as a homework to those who read to here. It would also be interesting to see results for newer (than 5S) iPhone models. The program during the calculations shows the estimated execution time, the result is saved in the iTunes file sharing folder.

Edit
In the comments prompted to compare the program with implementations on Google Play. As it turned out, it uses the algorithm Arithmetic-Geometric Mean (AGM), which is indeed much faster. With the latest version, the calculation of a million characters on the iPhone 5S took 68 seconds .
A description of the algorithm and implementation in Python is here: http://www.johndcook.com/blog/2012/06/03/calculating-pi-with-agm-and-mpmath/ . C / C ++ code for GMP is available here: https://rosettacode.org/wiki/Arithmetic-geometric_mean/Calculate_Pi .

Those who want to test their memory result. Here on the site is a text file containing a billion characters of Pi, which itself was used to verify the correctness of the program results.

Why is this necessary?


As for the world record (which is 5 trillion characters for the desktop, though the desktop was not ordinary ), it is unlikely to be beaten on the smartphone. Yes, and like this table of records for smartphones like no one else did. However, the matter is certainly not in the record - as you can see, such calculations are a rather interesting engineering task, where the questions of optimization and memory, and speed of calculations, and multithreading remain open. In general, there is a place to “wiggle”, which is just the most interesting in programming.

Source: https://habr.com/ru/post/309674/


All Articles