return to home

some of my

Adventures with Pi

Temp storage, good on-line pdf compressor https://www.freepdfconvert.com/compress-pdf#

by Ed Thelen

The only Vice I admit to ;-)) I compute Pi
REALLY serious Pi computations
Pi calculated to 50 trillion digits' with a pair of 32-core AMD Epyc chips, 1TB RAM, 510TB disk space
Also see Pi on a decimal character machine
and I want to remember Cornell scientific preprint server .

Table of Contents:

Introduction
A fast PC in 1992 - 66 MHz
adventures with a 32 bit ARM development board
Using a much much faster algorithm in 2019 - "Go For Broke"
a 2020 adventure by Ronald Mak
a 2020 adventure by Ed Thelen
Dream 2020 - Using a cell phone ;-) with fast 64 bit ARM processor and much more memory.
Outside Resources

Introduction

Long ago, about 1964, and far away, well - OK - Phoenix, Arizona
I was getting paid to play with computers :-)) how about that for fun ?? !!
So I got to enjoy computing PI

A fast PC in 1992 - 66 MHz

So along the way, about 1992 I got an INTEL based PC with with a double clocked 33 MHz CPU - which had an internal rate of an astounding 66 MHz !! Well, just how fast was that ?? My obvious test was to compute Pi to 1,000,000 places as I had done in 1967 with the then world's fastest computer ( the Control Data 6600 ).
Unfortunately, the DOS world had/has a severe memory addressing limitations, so a work around was to buy something called a DOS Extender from Phar Lap that permitted easier large memory addressing.
So - the tale starts - being lazy about I/O I decided to format lines and have some language ?BASIC? write the results to disk. There is an advantage in old fashioned machines, and the INTEL 80x86 family in the divide instruction - both the quotient and remainder are available from one divide :-)) In many methods of computing Pi to many significant places, having the remainder available is very helpful in arbitrary precision divide. Unfortunately, no modern language can take advantage of this quirk. You divide once to get the quotient and you have to do additional processing to get the remainder -
See Techie Notes
To make the assembly code, there was a facility in the BASIC to assemble code and interface the code with symbolic names, etc.
Unfortunately, I don't remember much about the attempt.

I do remember I was sharp on 80386 type code as I sometimes used it at work.
A major tussle was to get the address symbols in the "higher level language" useful in the assembly language program.
I printed out the result in nice format, a stack of paper a little over 1 inch thick.
and gave it to Dan Larson as a house warming present - he seemed disappointed :-((

adventures with a 32 bit ARM development board

Introduction
Overview and results
Development environment
Compute the arc tan (1/n)
Crippled Divide and work around
Notes
The total "C" code
Pi to 20,000 digits
Introduction

Long ago, about 1964, and far away, well - OK - Phoenix, Arizona
I was getting paid to play with computers :-)) how about that for fun ?? !!
So I got to enjoy computing PI
Recently there has been a flood of very low priced "evaluation" boards so that you can try out and become familiar with some vendor's computer chips. I finally decided to become modern and try one.
This Quick Start Guide with Lessons 0-10 on video seemed attractive. Lesson 0 tells how to obtain the free development software with control interface and simulator. The cheap USB powered $13 board with flashing lights caught my attention ;-)) - a 32 bit machine w evaluation kit . A TM4C123GH6PM, a 32-bit ARM Cortex-M4-based microcontroller with 256-kB Flash memory, 32 kB SRAM and 80-MHz operation.
Lesson 0 also tells how to order this board. (The control interface can also work with many other boards.)

Overview and results

So I ordered the board from Ti, and started watching and doing the lessons, using the simulator, until the hardware arrived.

When the hardware arrived, I had trouble getting the software to recognize the hardware. ( OK, I'm not real good at following directions.) Finally called up and was guided to some unobvious box to click. Magic happened - all became happy - the simulated code in the lessons worked just fine in the Ti hardware.

Lesson # 10 complete - time to compute Pi ;-))
But how to display/output the result???
The results of some experimentation yielded:

doing "printf("Hello World xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx\n");
to the console output box in the development screen is slow, about 10 characters/second.
stdio.h seemed stripped of file I/O functions such as "fopen( );" and "fclose( );"
so that path seemed out
I settled on displaying the decimal values as part of a HEX formatted display, example

Each hex digit contains the decimal value of the decimal value in that position. That way the hex memory dump (available "instantly" is readable and comparable to decimal values for Pi in books and on the Internet :-)) - No need for a printer ;-))
Unfortunately, the formatting is not ideal, nor is the formatting in the references ;-)

OK, we seem to be getting valid looking results compared to references :-))
Now, the usual thing is to
- optimize the code for faster results
- optimize memory for more possible digits in available memory
Tweaking things, and reusing one of the computing fields to display results, I was able to make this board with 32K bytes of SRAM output "decimal" 20,000 digits (see picture above for a little of the output) in under 3 minutes :-))
(The compiled software is executed out of e-prom, not part of the SRAM.)
And none of the chips on the board felt warm to the touch ;-))

Development environment

As mentioned above, I used Quick Start Guide. Which recommended a development system from IAR Systems. Use the long term but size limited option as mentioned in the Quick Start Guide.

The two main user displays are
edit
run

Compute the arctan (1/n)

We will use the 300 year old Machin formula for computing Pi.

which uses two arctangents.
A convenient formula for the arctan is the series

One of the arctangents is for the angle 1/5 radians, so the formula for that is
arctan(1/5) = 1/5 - 1/(3 x 5^3) + 1/(5 x 5^5) - ...
which can be rewritten

arctan(1/5) = 1/5 - 1/((2n+1) x 5^(2n+1) + ... - ...
n=1 n=2 n=3
Lots of divisions, with alternating signs of successive terms.
This can conveniently computed by:

zero arrays of unsigned integers called d1, d2, and acc
set i to 1
assume a decimal point to the right of the left (high order) integer in array d1
place a 1 as the left integer in array d1
unsigned divide array d1 by 5 giving array d1
unsigned divide array d1 by i giving array d2
add 2 to i
add (alternately subtract) array d2 and array acc giving array acc
unsigned divide array d1 by 5^2 or 25 giving array d1
if array d1 not equal to zero, go to 6 above
The arctan of 1/5 (in this case) is in array acc
A slightly optimized version in "C" follows:
// -------------------------- // int atan( divisor, initial_add_sub ), // does NOT clear accumulator !! // division by 239^2 works just fine void atan( int recip_radians, int sub1st_flag ); void atan( int recip_radians, int sub1st_flag ){ unsigned int loc_divisor = recip_radians*recip_radians ; int loc_flag = sub1st_flag; int all0s = 0; // initialize atan_d1 field zerof( atan_d1, f_len); // zero most working fields zerof( atan_d2, f_len); atan_d1[0] = 1; // prepare for recip radians div( atan_d1, atan_d1, recip_radians, all0s, f_len ); // d1 correct if ( (loc_flag&1) == 0) add( atan_d1, atan_acc, all0s, f_len ); // add else sub( atan_d1, atan_acc, all0s, f_len ); // subtract loc_flag ^= 1; // flip flag // fields correctly set for taylor series unsigned int divn = 3; // initial t2 divisor // divide d1 by divisor while (all0s != (f_len)){ //void div( unsigned short *fi, unsigned short *fo, // unsigned int divisor, int start, int len ); div( atan_d1, atan_d1, loc_divisor, all0s, f_len ); // divide d1 giving d1 div( atan_d1, atan_d2, divn, all0s, f_len ); // divide d1 giving d2 divn += 2; // update, prep for next round if (divn > max_ndiv) max_ndiv = divn; // remember max size if ( (loc_flag&1) == 0) add( atan_d2, atan_acc, all0s, f_len ); // add else sub( atan_d2, atan_acc, all0s, f_len ); // subtract loc_flag ^= 1; // flip flag // if ( (atan_d1[all0s] ^ (atan_d1[all0s+1] ) == 0) { if ( atan_d1[all0s] == 0) ++all0s; // start at next word next time } }

Crippled Divide and work around

"In the goode olde dayze", when computers were simple, they had two basic computing registers, usually called "A" and "B" or "Q". They were often concatenated for shifting.
A multiply, which can yield about twice as many digits or bits as the inputs usually used the concatenated "A" and "B" registers for the results.
A divide usually used the concatenated "A" and "B" registers for the dividend, and left the quotient in the "A" register and the remainder in the "B" register.
Most Pi algorithms do lots of dividing, and to do extended precision work, you need the remainder.
Then two things happened, both inconvenient for Pi finders

"Higher Level Languages" did not make the remainders available, or if they did, they caused another divide to get the remainder (again)
Computers came with "register files", 6 or 14 or 29 or ... general purpose registers -
Some of the computers placed the remainder in the next odd # register if the quotient was in an even numbered register
However, modern RISC computers just trash the remainder, if you want it, compute it.
Oddly, since most people don't want remainders, and High Level Languages don't know what to do with 'em, the results are acceptable to most people.
Work Around
Assume your divisor is positive, and no larger than than 1/2 your register lengths
Example (32 bit registers with 16 bit divisors max, and an unsigned divide instruction)
then the following "C" code will work
// -------------------------- void div( unsigned short *fi, unsigned short *fo, unsigned int divisor, int start, int len ); // "fi" is input array, "fo" is output array, may be same area, // "start" is where to start dividing, to save time because high order becomes all zeros // "len" is total length of working array void div( unsigned short *fi, unsigned short *fo, unsigned int divisor, int start, int len ){ int index = start; // get to register int endi = len; unsigned int remainder = 0; // storage for remainder unsigned int t = 0; // storage for quotient unsigned int dividend = 0; // dividend /// ( a 32 bit register, with register zeroing rules ) while (index <= endi) { dividend = fi[index]; // low 16 bits of dividend, //hi order zeroed because it is a register dividend |= remainder << (sizeof(unsigned short)*8); // form hi order of dividend (register) t = dividend/divisor; // divide to get quotient fo[index] = t; // quotient, send to array of 16 bit unsigned ints remainder = dividend - (t*divisor); // form remainder ++index; // go hi order to low order } }

Notes

unsigned int ;-)) The compiler seems to do a good job using unsigned multiplies and divides :-)) The high order bit is not treated as a special bit - so I can use all 16 bits using 239^2 = 57121 as a divisor :-))
The available SRAM in the machine is 32K bytes, starting at 0x20000000 . (The program is stored and run from e-prom, in low memory.)
- I wound up using 3 fields of 5006 16 bit words, reusing the zeroed d1 array as the display array.
- This producing 20,000 decimal characters.
- The last 8 characters matched values from the Internet which I accepted as valid.
- The time with the big divisors was 178 seconds.
Randy was wondering if he could multiply by the reciprocal of 5 instead of dividing by 5. The leading hex digits of 1/5 is 0.333333333333333333333333333333 ... in hex The leading hex digits of 1/25 is also a repeating fraction as are 1/239 and 1/(239^2) . I think a good consideration, but no cigar :-((

> How many decimal digits is that array worth?
I was figuring on outputting one decimal digit for every hex digit. Then I thought that wasteful - then thought some more and decided that I could really get 9 decimal digits (plus a little) for every 8 hex digits.
Now - how did I figure that ?
     4 hex digits (unsigned) = 65,536 decimal which is 4 decimal digits plus a little over half a digit
     8 hex digits (unsigned) = 4,294,967,296 decimal which is 9 decimal digits plus a little under half a digit
     so, being truly lazy, I emit 4 decimal digits per 4 hex digits
Experimenting with printing more digits until low order digits are incorrect indicate that it is "safe" to continue printing 14 % more digits, but 16 % more yields low order errors.
Now - about getting those digits printed, and in my usual form
     nnnnn nnnnn nnnnn nnnnn nnnnn nnnnn nnnnn nnnnn nnnnn nnnnn * 10^-eeeeee0
What about errors creeping in from the truncated least significant part?
I think that the the division process sweeps digitization errors "out the low end", but the adding and subtracting into the accumulator causes error to accumulate. Using 5006 16 bit integers, there are about 16,000 iterations of dividing by 25 and adding/subtracting the division by (2n-1) and less than 10% of that in the process of dividing by 239^2. Assuming a maximum digitation error of 1 count per addition/subtraction would yield a maximum error of 1 16 bit int or less than 4 hex digits ( about 65,000 ).
We are dropping off 1 decimal digit per 8 decimal digits "printed". As there are 5006/2 8 hex digit units, for our 20,000 decimal digit conversion, we are dropping of about 2,500 decimal digits - a much more than safe output truncation :-))

The total "C" code

/* Pi x */ // there is a limit of 32K bytes SRAM starting at 0x20000000 // this version does over 20,000 decimal digits in 170 seconds // // game plan // Pi/4 = 4arctan(1/5)-arctan(1/239) // pi = 16arctan(1/5)-4arctan(1/239) // 5^2 = 25 239^2 = 57121 2^16 = 65,536 // the compiler generates unsigned divide and multiply which work well // a) make 3 large areas, atan_d1[], atan_d2[], atan_acc[] // area atan_d1[] is reused as an output field // b) make subroutines div( *fi, *fo, divisor, start, len ), // int add( *fi, *fo, start, len ), // how test for overflow? // sub( *fi, *fo, start, len ), // how test for underflow? // atan( divisor, initial_add_sub ), // status is some error - // how to handle division by 239^2 ?? // do /239 twice ?? // zerof( *fi, len ), #define f_len 5006 // actual field length - 2, allow two guard ints // // 5006 = 4096 + 512 + 256 + 128 +16 - 2, // time was 2 min +50 sec max_divn = 33,499 // over 20,,000 decimal digits // int is 16 bits - static unsigned int max_ndiv; static int add_sub_error; static unsigned short atan_d1[f_len+2], atan_d2[f_len+2], atan_acc[f_len+2]; // guard at end // note: atan_d1 is re-used as the decimal output field // -------------------------- void zerof( unsigned short *fi, int lim ); void zerof( unsigned short *fi, int lim ){ int counter = 0; while ( counter < lim) { fi[counter] = 0x0000; ++counter; } } // -------------------------- // add fi into fo void add( unsigned short *fi, unsigned short *fo, int zeros, int len ); void add( unsigned short *fi, unsigned short *fo, int zeros, int len ){ int carry = 0; // yes, initialize carry int index = len; // start at low end int t; while ( (index >= zeros) || (carry != 0) ) { t = carry + fi[index] + fo[index]; fo[index] = t; carry = t >> (sizeof(unsigned short)*8); --index; } if ( index < -1 ) add_sub_error |= 1; // bit 1 is add overflow } // -------------------------- // sub fi from fo // error if return not 0 void sub( unsigned short *fi, unsigned short *fo, int zeros, int len ); void sub( unsigned short *fi, unsigned short *fo, int zeros, int len ){ int borrow = 0; // yes, initialize borrow int index = len; // start at low end int t; while ( index >= zeros ) { t = fo[index] - borrow - fi[index]; fo[index] = t; borrow = t >> (sizeof(unsigned short)*8) & 1; --index; } if ( index < -1 ) add_sub_error |= 2; // bit 2 is subtract underflow } // -------------------------- // void mpy( *fio, start, multiplier ), // multiply in place by multiplier // return is int part[0] void mpy( unsigned short *fio, int start, int multiplier ); void mpy( unsigned short *fio, int start, int multiplier ){ unsigned int carry = 0; // hi part of product int index = start; // start at low end unsigned int t; while ( index >= 0 ) { // t = carry + fi[index] + fo[index]; t = ( fio[index] * multiplier) + carry; fio[index] = t; carry = t >> (sizeof(unsigned short)*8); --index; } } // -------------------------- void div( unsigned short *fi, unsigned short *fo, unsigned int divisor, int start, int len ); // "fi" is input array, "fo" is output array, may be same area, // "start" is where to start dividing, to save time because high order becomes all zeros // "len" is total length of working array void div( unsigned short *fi, unsigned short *fo, unsigned int divisor, int start, int len ){ int index = start; // get to register int endi = len; unsigned int remainder = 0; // storage for remainder unsigned int t = 0; // storage for quotient unsigned int dividend = 0; // dividend /// ( a 32 bit register, with register zeroing rules ) while (index <= endi) { dividend = fi[index]; // low 16 bits of dividend, //hi order zeroed because it is a register dividend |= remainder << (sizeof(unsigned short)*8); // form hi order of dividend (register) t = dividend/divisor; // divide to get quotient fo[index] = t; // quotient, send to array of 16 bit unsigned ints remainder = dividend - (t*divisor); // form remainder ++index; // go hi order to low order } } // -------------------------- // int atan( divisor, initial_add_sub ), // does NOT clear accumulator !! // division by 239^2 works just fine void atan( int recip_radians, int sub1st_flag ); void atan( int recip_radians, int sub1st_flag ){ unsigned int loc_divisor = recip_radians*recip_radians ; int loc_flag = sub1st_flag; int all0s = 0; // initialize atan_d1 field zerof( atan_d1, f_len); // zero most working fields zerof( atan_d2, f_len); atan_d1[0] = 1; // prepare for recip radians div( atan_d1, atan_d1, recip_radians, all0s, f_len ); // d1 correct if ( (loc_flag&1) == 0) add( atan_d1, atan_acc, all0s, f_len ); // add else sub( atan_d1, atan_acc, all0s, f_len ); // subtract loc_flag ^= 1; // flip flag // fields correctly set for taylor series unsigned int divn = 3; // initial t2 divisor // divide d1 by divisor while (all0s != (f_len)){ //void div( unsigned short *fi, unsigned short *fo, // unsigned int divisor, int start, int len ); div( atan_d1, atan_d1, loc_divisor, all0s, f_len ); // divide d1 giving d1 div( atan_d1, atan_d2, divn, all0s, f_len ); // divide d1 giving d2 divn += 2; // update, prep for next round if (divn > max_ndiv) max_ndiv = divn; // remember max size if ( (loc_flag&1) == 0) add( atan_d2, atan_acc, all0s, f_len ); // add else sub( atan_d2, atan_acc, all0s, f_len ); // subtract loc_flag ^= 1; // flip flag // if ( (atan_d1[all0s] ^ (atan_d1[all0s+1] ) == 0) { if ( atan_d1[all0s] == 0) ++all0s; // start at next word next time } } // form ascii array from atan_acc, atan_acc is destroyed // conv_bin_char( int len); void conv_bin_char( int len); void conv_bin_char( int len){ int t1,t2,t3,i,wdivisor; int index = 0; zerof( atan_d1, f_len); // changed from ascii while ( index <= len) { t1 = 0; // initialize display t2 = 0; // initialize working wdivisor = 1000; i = 3; t3 = atan_acc[0]; while (wdivisor > 0 ) { t2 = t3/wdivisor; // form a "decimal" digit t1 |= t2 << (i*4); // insert next "decimal" digit t3 = t3 - (t2*wdivisor); // remainder wdivisor /= 10; --i; } atan_d1[index] = t1; // insert "decimal" into output field atan_acc[0] = 0; // clear hi order accum mpy( atan_acc, f_len - (index/2), 10000 ); // multiply, shorting field ++index; // go to next word } } //--------------------------------------------------------------- int main() { max_ndiv = 0; // clear warns and errors add_sub_error = 0; zerof( atan_acc, f_len); // clear result atan( 5, 0 ); // form arctan (1/5) mpy(atan_acc, f_len, 4 ); // follow Machin atan( 239, 1 ); // form arctan (1/239) mpy(atan_acc, f_len, 4 ); // follow Machin conv_bin_char( f_len); // "print" it out return 0; }

Pi to 20,000 digits using the Ti 32 bit ARM development board, Work done in 2012 ??

Images of the "decimalized hex digits" in a memory dump of SRAM in the Ti board are available in 11 frames as follows:
1st, 2nd, 3rd, 4th, 5th, 6th, 7th, 8th, 9th, 10th, 11th
I added a little formatting and printing code to the above "C" program which sent about 20,000 digits of Pi to the "Console" display of the user interface, at about 10 characters/second. And here they are :-))
00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 xE -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 xE -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 xE -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 xE -350 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 xE -400 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 xE -450 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 xE -500 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 xE -550 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132 xE -600 00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 xE -650 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 xE -700 42019 95611 21290 21960 86403 44181 59813 62977 47713 09960 xE -750 51870 72113 49999 99837 29780 49951 05973 17328 16096 31859 xE -800 50244 59455 34690 83026 42522 30825 33446 85035 26193 11881 xE -850 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 xE -900 59825 34904 28755 46873 11595 62863 88235 37875 93751 95778 xE -950 18577 80532 17122 68066 13001 92787 66111 95909 21642 01989 xE -1000 38095 25720 10654 85863 27886 59361 53381 82796 82303 01952 xE -1050 03530 18529 68995 77362 25994 13891 24972 17752 83479 13151 xE -1100 55748 57242 45415 06959 50829 53311 68617 27855 88907 50983 xE -1150 81754 63746 49393 19255 06040 09277 01671 13900 98488 24012 xE -1200 85836 16035 63707 66010 47101 81942 95559 61989 46767 83744 xE -1250 94482 55379 77472 68471 04047 53464 62080 46684 25906 94912 xE -1300 93313 67702 89891 52104 75216 20569 66024 05803 81501 93511 xE -1350 25338 24300 35587 64024 74964 73263 91419 92726 04269 92279 xE -1400 67823 54781 63600 93417 21641 21992 45863 15030 28618 29745 xE -1450 55706 74983 85054 94588 58692 69956 90927 21079 75093 02955 xE -1500 32116 53449 87202 75596 02364 80665 49911 98818 34797 75356 xE -1550 63698 07426 54252 78625 51818 41757 46728 90977 77279 38000 xE -1600 81647 06001 61452 49192 17321 72147 72350 14144 19735 68548 xE -1650 16136 11573 52552 13347 57418 49468 43852 33239 07394 14333 xE -1700 45477 62416 86251 89835 69485 56209 92192 22184 27255 02542 xE -1750 56887 67179 04946 01653 46680 49886 27232 79178 60857 84383 xE -1800 82796 79766 81454 10095 38837 86360 95068 00642 25125 20511 xE -1850 73929 84896 08412 84886 26945 60424 19652 85022 21066 11863 xE -1900 06744 27862 20391 94945 04712 37137 86960 95636 43719 17287 xE -1950 46776 46575 73962 41389 08658 32645 99581 33904 78027 59009 xE -2000 94657 64078 95126 94683 98352 59570 98258 22620 52248 94077 xE -2050 26719 47826 84826 01476 99090 26401 36394 43745 53050 68203 xE -2100 49625 24517 49399 65143 14298 09190 65925 09372 21696 46151 xE -2150 57098 58387 41059 78859 59772 97549 89301 61753 92846 81382 xE -2200 68683 86894 27741 55991 85592 52459 53959 43104 99725 24680 xE -2250 84598 72736 44695 84865 38367 36222 62609 91246 08051 24388 xE -2300 43904 51244 13654 97627 80797 71569 14359 97700 12961 60894 xE -2350 41694 86855 58484 06353 42207 22258 28488 64815 84560 28506 xE -2400 01684 27394 52267 46767 88952 52138 52254 99546 66727 82398 xE -2450 64565 96116 35488 62305 77456 49803 55936 34568 17432 41125 xE -2500 15076 06947 94510 96596 09402 52288 79710 89314 56691 36867 xE -2550 22874 89405 60101 50330 86179 28680 92087 47609 17824 93858 xE -2600 90097 14909 67598 52613 65549 78189 31297 84821 68299 89487 xE -2650 22658 80485 75640 14270 47755 51323 79641 45152 37462 34364 xE -2700 54285 84447 95265 86782 10511 41354 73573 95231 13427 16610 xE -2750 21359 69536 23144 29524 84937 18711 01457 65403 59027 99344 xE -2800 03742 00731 05785 39062 19838 74478 08478 48968 33214 45713 xE -2850 86875 19435 06430 21845 31910 48481 00537 06146 80674 91927 xE -2900 81911 97939 95206 14196 63428 75444 06437 45123 71819 21799 xE -2950 98391 01591 95618 14675 14269 12397 48940 90718 64942 31961 xE -3000 56794 52080 95146 55022 52316 03881 93014 20937 62137 85595 xE -3050 66389 37787 08303 90697 92077 34672 21825 62599 66150 14215 xE -3100 03068 03844 77345 49202 60541 46659 25201 49744 28507 32518 xE -3150 66600 21324 34088 19071 04863 31734 64965 14539 05796 26856 xE -3200 10055 08106 65879 69981 63574 73638 40525 71459 10289 70641 xE -3250 40110 97120 62804 39039 75951 56771 57700 42033 78699 36007 xE -3300 23055 87631 76359 42187 31251 47120 53292 81918 26186 12586 xE -3350 73215 79198 41484 88291 64470 60957 52706 95722 09175 67116 xE -3400 72291 09816 90915 28017 35067 12748 58322 28718 35209 35396 xE -3450 57251 21083 57915 13698 82091 44421 00675 10334 67110 31412 xE -3500 67111 36990 86585 16398 31501 97016 51511 68517 14376 57618 xE -3550 35155 65088 49099 89859 98238 73455 28331 63550 76479 18535 xE -3600 89322 61854 89632 13293 30898 57064 20467 52590 70915 48141 xE -3650 65498 59461 63718 02709 81994 30992 44889 57571 28289 05923 xE -3700 23326 09729 97120 84433 57326 54893 82391 19325 97463 66730 xE -3750 58360 41428 13883 03203 82490 37589 85243 74417 02913 27656 xE -3800 18093 77344 40307 07469 21120 19130 20330 38019 76211 01100 xE -3850 44929 32151 60842 44485 96376 69838 95228 68478 31235 52658 xE -3900 21314 49576 85726 24334 41893 03968 64262 43410 77322 69780 xE -3950 28073 18915 44110 10446 82325 27162 01052 65227 21116 60396 xE -4000 66557 30925 47110 55785 37634 66820 65310 98965 26918 62056 xE -4050 47693 12570 58635 66201 85581 00729 36065 98764 86117 91045 xE -4100 33488 50346 11365 76867 53249 44166 80396 26579 78771 85560 xE -4150 84552 96541 26654 08530 61434 44318 58676 97514 56614 06800 xE -4200 70023 78776 59134 40171 27494 70420 56223 05389 94561 31407 xE -4250 11270 00407 85473 32699 39081 45466 46458 80797 27082 66830 xE -4300 63432 85878 56983 05235 80893 30657 57406 79545 71637 75254 xE -4350 20211 49557 61581 40025 01262 28594 13021 64715 50979 25923 xE -4400 09907 96547 37612 55176 56751 35751 78296 66454 77917 45011 xE -4450 29961 48903 04639 94713 29621 07340 43751 89573 59614 58901 xE -4500 93897 13111 79042 97828 56475 03203 19869 15140 28708 08599 xE -4550 04801 09412 14722 13179 47647 77262 24142 54854 54033 21571 xE -4600 85306 14228 81375 85043 06332 17518 29798 66223 71721 59160 xE -4650 77166 92547 48738 98665 49494 50114 65406 28433 66393 79003 xE -4700 97692 65672 14638 53067 36096 57120 91807 63832 71664 16274 xE -4750 88880 07869 25602 90228 47210 40317 21186 08204 19000 42296 xE -4800 61711 96377 92133 75751 14959 50156 60496 31862 94726 54736 xE -4850 42523 08177 03675 15906 73502 35072 83540 56704 03867 43513 xE -4900 62222 47715 89150 49530 98444 89333 09634 08780 76932 59939 xE -4950 78054 19341 44737 74418 42631 29860 80998 88687 41326 04721 xE -5000 56951 62396 58645 73021 63159 81931 95167 35381 29741 67729 xE -5050 47867 24229 24654 36680 09806 76928 23828 06899 64004 82435 xE -5100 40370 14163 14965 89794 09243 23789 69070 69779 42236 25082 xE -5150 21688 95738 37986 23001 59377 64716 51228 93578 60158 81617 xE -5200 55782 97352 33446 04281 51262 72037 34314 65319 77774 16031 xE -5250 99066 55418 76397 92933 44195 21541 34189 94854 44734 56738 xE -5300 31624 99341 91318 14809 27777 10386 38773 43177 20754 56545 xE -5350 32207 77092 12019 05166 09628 04909 26360 19759 88281 61332 xE -5400 31666 36528 61932 66863 36062 73567 63035 44776 28035 04507 xE -5450 77235 54710 58595 48702 79081 43562 40145 17180 62464 36267 xE -5500 94561 27531 81340 78330 33625 42327 83944 97538 24372 05835 xE -5550 31147 71199 26063 81334 67768 79695 97030 98339 13077 10987 xE -5600 04085 91337 46414 42822 77263 46594 70474 58784 77872 01927 xE -5650 71528 07317 67907 70715 72134 44730 60570 07334 92436 93113 xE -5700 83504 93163 12840 42512 19256 51798 06941 13528 01314 70130 xE -5750 47816 43788 51852 90928 54520 11658 39341 96562 13491 43415 xE -5800 95625 86586 55705 52690 49652 09858 03385 07224 26482 93972 xE -5850 85847 83163 05777 75606 88876 44624 82468 57926 03953 52773 xE -5900 48030 48029 00587 60758 25104 74709 16439 61362 67604 49256 xE -5950 27420 42083 20856 61190 62545 43372 13153 59584 50687 72460 xE -6000 29016 18766 79524 06163 42522 57719 54291 62991 93064 55377 xE -6050 99140 37340 43287 52628 88963 99587 94757 29174 64263 57455 xE -6100 25407 90914 51357 11136 94109 11939 32519 10760 20825 20261 xE -6150 87985 31887 70584 29725 91677 81314 96990 09019 21169 71737 xE -6200 27847 68472 68608 49003 37702 42429 16513 00500 51683 23364 xE -6250 35038 95170 29893 92233 45172 20138 12806 96501 17844 08745 xE -6300 19601 21228 59937 16231 30171 14448 46409 03890 64495 44400 xE -6350 61986 90754 85160 26327 50529 83491 87407 86680 88183 38510 xE -6400 22833 45085 04860 82503 93021 33219 71551 84306 35455 00766 xE -6450 82829 49304 13776 55279 39751 75461 39539 84683 39363 83047 xE -6500 46119 96653 85815 38420 56853 38621 86725 23340 28308 71123 xE -6550 28278 92125 07712 62946 32295 63989 89893 58211 67456 27010 xE -6600 21835 64622 01349 67151 88190 97303 81198 00497 34072 39610 xE -6650 36854 06643 19395 09790 19069 96395 52453 00545 05806 85501 xE -6700 95673 02292 19139 33918 56803 44903 98205 95510 02263 53536 xE -6750 19204 19947 45538 59381 02343 95544 95977 83779 02374 21617 xE -6800 27111 72364 34354 39478 22181 85286 24085 14006 66044 33258 xE -6850 88569 86705 43154 70696 57474 58550 33232 33421 07301 54594 xE -6900 05165 53790 68662 73337 99585 11562 57843 22988 27372 31989 xE -6950 87571 41595 78111 96358 33005 94087 30681 21602 87649 62867 xE -7000 44604 77464 91599 50549 73742 56269 01049 03778 19868 35938 xE -7050 14657 41268 04925 64879 85561 45372 34786 73303 90468 83834 xE -7100 36346 55379 49864 19270 56387 29317 48723 32083 76011 23029 xE -7150 91136 79386 27089 43879 93620 16295 15413 37142 48928 30722 xE -7200 01269 01475 46684 76535 76164 77379 46752 00490 75715 55278 xE -7250 19653 62132 39264 06160 13635 81559 07422 02020 31872 77605 xE -7300 27721 90055 61484 25551 87925 30343 51398 44253 22341 57623 xE -7350 36106 42506 39049 75008 65627 10953 59194 65897 51413 10348 xE -7400 22769 30624 74353 63256 91607 81547 81811 52843 66795 70611 xE -7450 08615 33150 44521 27473 92454 49454 23682 88606 13408 41486 xE -7500 37767 00961 20715 12491 40430 27253 86076 48236 34143 34623 xE -7550 51897 57664 52164 13767 96903 14950 19108 57598 44239 19862 xE -7600 91642 19399 49072 36234 64684 41173 94032 65918 40443 78051 xE -7650 33389 45257 42399 50829 65912 28508 55582 15725 03107 12570 xE -7700 12668 30240 29295 25220 11872 67675 62204 15420 51618 41634 xE -7750 84756 51699 98116 14101 00299 60783 86909 29160 30288 40026 xE -7800 91041 40792 88621 50784 24516 70908 70006 99282 12066 04183 xE -7850 71806 53556 72525 32567 53286 12910 42487 76182 58297 65157 xE -7900 95984 70356 22262 93486 00341 58722 98053 49896 50226 29174 xE -7950 87882 02734 20922 22453 39856 26476 69149 05562 84250 39127 xE -8000 57710 28402 79980 66365 82548 89264 88025 45661 01729 67026 xE -8050 64076 55904 29099 45681 50652 65305 37182 94127 03369 31378 xE -8100 51786 09040 70866 71149 65583 43434 76933 85781 71138 64558 xE -8150 73678 12301 45876 87126 60348 91390 95620 09939 36103 10291 xE -8200 61615 28813 84379 09904 23174 73363 94804 57593 14931 40529 xE -8250 76347 57481 19356 70911 01377 51721 00803 15590 24853 09066 xE -8300 92037 67192 20332 29094 33467 68514 22144 77379 39375 17034 xE -8350 43661 99104 03375 11173 54719 18550 46449 02636 55128 16228 xE -8400 82446 25759 16333 03910 72253 83742 18214 08835 08657 39177 xE -8450 15096 82887 47826 56995 99574 49066 17583 44137 52239 70968 xE -8500 34080 05355 98491 75417 38188 39994 46974 86762 65516 58276 xE -8550 58483 58845 31427 75687 90029 09517 02835 29716 34456 21296 xE -8600 40435 23117 60066 51012 41200 65975 58512 76178 58382 92041 xE -8650 97484 42360 80071 93045 76189 32349 22927 96501 98751 87212 xE -8700 72675 07981 25547 09589 04556 35792 12210 33346 69749 92356 xE -8750 30254 94780 24901 14195 21238 28153 09114 07907 38602 51522 xE -8800 74299 58180 72471 62591 66854 51333 12394 80494 70791 19153 xE -8850 26734 30282 44186 04142 63639 54800 04480 02670 49624 82017 xE -8900 92896 47669 75831 83271 31425 17029 69234 88962 76684 40323 xE -8950 26092 75249 60357 99646 92565 04936 81836 09003 23809 29345 xE -9000 95889 70695 36534 94060 34021 66544 37558 90045 63288 22505 xE -9050 45255 64056 44824 65151 87547 11962 18443 96582 53375 43885 xE -9100 69094 11303 15095 26179 37800 29741 20766 51479 39425 90298 xE -9150 96959 46995 56576 12186 56196 73378 62362 56125 21632 08628 xE -9200 69222 10327 48892 18654 36480 22967 80705 76561 51446 32046 xE -9250 92790 68212 07388 37781 42335 62823 60896 32080 68222 46801 xE -9300 22482 61177 18589 63814 09183 90367 36722 20888 32151 37556 xE -9350 00372 79839 40041 52970 02878 30766 70944 47456 01345 56417 xE -9400 25437 09069 79396 12257 14298 94671 54357 84687 88614 44581 xE -9450 23145 93571 98492 25284 71605 04922 12424 70141 21478 05734 xE -9500 55105 00801 90869 96033 02763 47870 81081 75450 11930 71412 xE -9550 23390 86639 38339 52942 57869 05076 43100 63835 19834 38934 xE -9600 15961 31854 34754 64955 69781 03829 30971 64651 43840 70070 xE -9650 73604 11237 35998 43452 25161 05070 27056 23526 60127 64848 xE -9700 30840 76118 30130 52793 20542 74628 65403 60367 45328 65105 xE -9750 70658 74882 25698 15793 67897 66974 22057 50596 83440 86973 xE -9800 50201 41020 67235 85020 07245 22563 26513 41055 92401 90274 xE -9850 21624 84391 40359 98953 53945 90944 07046 91209 14093 87001 xE -9900 26456 00162 37428 80210 92764 57931 06579 22955 24988 72758 xE -9950 46101 26483 69998 92256 95968 81592 05600 10165 52563 75678 xE-10000 56672 27966 19885 78279 48488 55834 39751 87445 45512 96563 xE-10050 44348 03966 42055 79829 36804 35220 27709 84294 23253 30225 xE-10100 76341 80703 94769 94159 79159 45300 69752 14829 33665 55661 xE-10150 56787 36400 53666 56416 54732 17043 90352 13295 43529 16941 xE-10200 45990 41608 75320 18683 79370 23488 86894 79151 07163 78529 xE-10250 02345 29244 07736 59495 63051 00742 10871 42613 49745 95615 xE-10300 13849 87137 57047 10178 79573 10422 96906 66702 14498 63746 xE-10350 45952 80824 36944 57897 72330 04876 47652 41339 07592 04340 xE-10400 19634 03911 47320 23380 71509 52220 10682 56342 74716 46024 xE-10450 33544 00515 21266 93249 34196 73977 04159 56837 53555 16673 xE-10500 02739 00749 72973 63549 64533 28886 98440 61196 49616 27734 xE-10550 49518 27369 55882 20757 35517 66515 89855 19098 66653 93549 xE-10600 48106 88732 06859 90754 07923 42402 30092 59007 01731 96036 xE-10650 22547 56478 94064 75483 46647 76041 14632 33905 65134 33068 xE-10700 44953 97907 09030 23460 46147 09616 96886 88501 40834 70405 xE-10750 46074 29586 99138 29668 24681 85710 31887 90652 87036 65083 xE-10800 24319 74404 77185 56789 34823 08943 10682 87027 22809 73624 xE-10850 80939 96270 60747 26455 39925 39944 28081 13736 94338 87294 xE-10900 06307 92615 95995 46262 46297 07062 59484 55690 34711 97299 xE-10950 64090 89418 05953 43932 51236 23550 81349 49004 36427 85271 xE-11000 38315 91256 89892 95196 42728 75739 46914 27253 43669 41532 xE-11050 36100 45373 04881 98551 70659 41217 35246 25895 48730 16760 xE-11100 02988 65925 78662 85612 49665 52353 38294 28785 42534 04830 xE-11150 83307 01653 72285 63559 15253 47844 59818 31341 12900 19992 xE-11200 05981 35220 51173 36585 64078 26484 94276 44113 76393 86692 xE-11250 48031 18364 45369 85891 75442 64739 98822 84621 84490 08777 xE-11300 69776 31279 57226 72655 56259 62825 42765 31830 01340 70922 xE-11350 33436 57791 60128 09317 94017 18598 59993 38492 35495 64005 xE-11400 70995 58561 13498 02524 99066 98423 30173 50358 04408 11685 xE-11450 52653 11709 95708 99427 32870 92584 87894 43646 00504 10892 xE-11500 26691 78352 58707 85951 29834 41729 53519 53788 55345 73742 xE-11550 60859 02908 17651 55780 39059 46408 73506 12322 61120 09373 xE-11600 10804 85485 26357 22825 76820 34160 50484 66277 50450 03126 xE-11650 20080 07998 04925 48534 69414 69775 16493 27095 04934 63938 xE-11700 24322 27188 51597 40547 02148 28971 11777 92376 12257 88734 xE-11750 77188 19682 54629 81268 68581 70507 40272 55026 33290 44976 xE-11800 27789 44236 21674 11918 62694 39650 67151 57795 86756 48239 xE-11850 93917 60426 01763 38704 54990 17614 36412 04692 18237 07648 xE-11900 87834 19689 68611 81558 15873 60629 38603 81017 12158 55272 xE-11950 66830 08238 34046 56475 88040 51380 80163 36388 74216 37140 xE-12000 64354 95561 86896 41122 82140 75330 26551 00424 10489 67835 xE-12050 28588 29024 36709 04887 11819 09094 94533 14421 82876 61810 xE-12100 31007 35477 05498 15968 07720 09474 69613 43609 28614 84941 xE-12150 78501 71807 79306 81085 46900 09445 89952 79424 39813 92135 xE-12200 05586 42219 64834 91512 63901 28038 32001 09773 86806 62877 xE-12250 92397 18014 61343 24457 26400 97374 25700 73592 10031 54150 xE-12300 89367 93008 16998 05365 20276 00727 74967 45840 02836 24053 xE-12350 46037 26341 65542 59027 60183 48403 06811 38185 51059 79705 xE-12400 66400 75094 26087 88573 57960 37324 51414 67867 03688 09880 xE-12450 60971 64258 49759 51380 69309 44940 15154 22221 94329 13021 xE-12500 73912 53835 59150 31003 33032 51117 49156 96917 45027 14943 xE-12550 31515 58854 03922 16409 72291 01129 03552 18157 62823 28318 xE-12600 23425 48326 11191 28009 28252 56190 20526 30163 91147 72473 xE-12650 31485 73910 77758 74425 38761 17465 78671 16941 47764 21441 xE-12700 11126 35835 53871 36101 10232 67987 75641 02468 24032 26483 xE-12750 46417 66369 80663 78576 81349 20453 02240 81972 78564 71983 xE-12800 96308 78154 32211 66912 24641 59117 76732 25326 43356 86146 xE-12850 18654 52226 81268 87268 44596 84424 16107 85401 67681 42080 xE-12900 88502 80054 14361 31462 30821 02594 17375 62389 94207 57136 xE-12950 27516 74573 18918 94562 83525 70441 33543 75857 53426 98699 xE-13000 47254 70316 56613 99199 96826 28247 27064 13362 22178 92390 xE-13050 31760 85428 94373 39356 18891 65125 04244 04008 95271 98378 xE-13100 73864 80584 72689 54624 38823 43751 78852 01439 56005 71048 xE-13150 11949 88423 90606 13695 73423 15590 79670 34614 91434 47886 xE-13200 36041 03182 35073 65027 78590 89757 82727 31305 04889 39890 xE-13250 09923 91350 33732 50855 98265 58670 89242 61242 94736 70193 xE-13300 90772 71307 06869 17092 64625 48423 24074 85503 66080 13604 xE-13350 66895 11840 09366 86095 46325 00214 58529 30950 00090 71510 xE-13400 58236 26729 32645 37382 10493 87249 96699 33942 46855 16483 xE-13450 26113 41461 10680 26744 66373 34375 34076 42940 26682 97386 xE-13500 52209 35701 62638 46485 28514 90362 93201 99199 68828 51718 xE-13550 39536 69134 52224 44708 04592 39660 28171 56551 56566 61113 xE-13600 59823 11225 06289 05854 91450 97157 55390 02439 31535 19090 xE-13650 21071 19457 30024 38801 76615 03527 08626 02537 88179 75194 xE-13700 78061 01371 50044 89917 21002 22013 35013 10601 63915 41589 xE-13750 57803 71177 92775 22597 87428 91917 91552 24171 89585 36168 xE-13800 05947 41234 19339 84202 18745 64925 64434 62392 53195 31351 xE-13850 03311 47639 49119 95072 85843 06583 61935 36932 96992 89837 xE-13900 91494 19394 06085 72486 39688 36903 26556 43642 16644 25760 xE-13950 79147 10869 98431 57337 49648 83529 27693 28220 76294 72823 xE-14000 81537 40996 15455 98798 25989 10937 17126 21828 30258 48112 xE-14050 38901 19682 21429 45766 75807 18653 80650 64870 26133 89282 xE-14100 29949 72574 53033 28389 63818 43944 77077 94022 84359 88341 xE-14150 00358 38542 38973 54243 95647 55568 40952 24844 55413 92394 xE-14200 10001 62076 93636 84677 64130 17819 65937 99715 57468 54194 xE-14250 63348 93748 43912 97423 91433 65936 04100 35234 37770 65888 xE-14300 67781 13949 86164 78747 14079 32638 58738 62473 28896 45643 xE-14350 59877 46676 38479 46650 40741 11825 65837 88784 54858 14896 xE-14400 29612 73998 41344 27260 86061 87245 54523 60643 15371 01127 xE-14450 46809 77870 44640 94758 28034 87697 58948 32824 12392 92960 xE-14500 58294 86191 96670 91895 80898 33201 21031 84303 40128 49511 xE-14550 62035 34280 14412 76172 85830 24355 98300 32042 02451 20728 xE-14600 72535 58119 58401 49180 96925 33950 75778 40006 74655 26031 xE-14650 44616 70508 27682 77222 35341 91102 63416 31571 47406 12385 xE-14700 04258 45988 41990 76112 87258 05911 39356 89601 43166 82831 xE-14750 76323 56732 54170 73420 81733 22304 62987 99280 49085 14094 xE-14800 79036 88786 87894 93054 69557 03072 61900 95020 76433 49335 xE-14850 91060 24545 08645 36289 35456 86295 85313 15337 18386 82656 xE-14900 17862 27363 71697 57741 83023 98600 65914 81616 40494 49650 xE-14950 11732 13138 95747 06208 84748 02365 37103 11508 98427 99275 xE-15000 44268 53277 97431 13951 43574 17221 97597 99359 68525 22857 xE-15050 45263 79628 96126 91572 35798 66205 73408 37576 68738 84266 xE-15100 40599 09935 05000 81337 54324 54635 96750 48442 35284 87470 xE-15150 14435 45419 57625 84735 64216 19813 40734 68541 11766 88311 xE-15200 86544 89377 69795 66517 27966 23267 14810 33864 39137 51865 xE-15250 94673 00244 34500 54499 53997 42372 32871 24948 34706 04406 xE-15300 34716 06325 83064 98297 95510 10954 18362 35030 30945 30973 xE-15350 35834 46283 94763 04775 64501 50085 07578 94954 89313 93944 xE-15400 89921 61255 25597 70143 68589 43585 87752 63796 25597 08167 xE-15450 76438 00125 43650 23714 12783 46792 61019 95585 22471 72201 xE-15500 77723 70041 78084 19423 94872 54068 01556 03599 83905 48985 xE-15550 72354 67456 42390 58585 02167 19031 39526 29445 54391 31663 xE-15600 13453 08939 06204 67843 87785 05423 93905 24731 36201 29476 xE-15650 91874 97519 10114 72315 28932 67725 33918 14660 73000 89027 xE-15700 76896 31148 10902 20972 45207 59167 29700 78505 80717 18638 xE-15750 10549 67973 10016 78708 50694 20709 22329 08070 38326 34534 xE-15800 52038 02786 09905 56900 13413 71823 68370 99194 95164 89600 xE-15850 75504 93412 67876 43674 63849 02063 96401 97666 85592 33565 xE-15900 46391 38363 18574 56981 47196 21084 10809 61884 60545 60390 xE-15950 38455 34372 91414 46513 47494 07848 84423 77217 51543 34260 xE-16000 30669 88317 68331 00113 31086 90421 93903 10801 43784 33415 xE-16050 13709 24353 01367 76310 84913 51615 64226 98475 07430 32971 xE-16100 67469 64066 65315 27035 32546 71126 67522 46055 11995 81831 xE-16150 96376 37076 17991 91920 35795 82007 59560 53023 46267 75794 xE-16200 39363 07463 05690 10801 14942 71410 09391 36913 81072 58137 xE-16250 81357 89400 55995 00183 54251 18417 21360 55727 52210 35268 xE-16300 03735 72652 79224 17373 60575 11278 87218 19084 49006 17801 xE-16350 38897 10770 82293 10027 97665 93583 87589 09395 68814 85602 xE-16400 63224 39372 65624 72776 03789 08144 58837 85501 97028 43779 xE-16450 36240 78250 52704 87581 64703 24581 29087 83952 32453 23789 xE-16500 60298 41669 22548 96497 15606 98119 21865 84926 77040 39564 xE-16550 81278 10217 99132 17416 30581 05545 98801 30048 45629 97651 xE-16600 12124 15363 74515 00563 50701 27815 92671 42413 42103 30156 xE-16650 61653 56024 73380 78430 28655 25722 27530 49998 83701 53487 xE-16700 93008 06260 18096 23815 16136 69033 41111 38653 85109 19367 xE-16750 39383 52293 45888 32255 08870 64507 53947 39520 43968 07906 xE-16800 70868 06445 09698 65488 01682 87434 37861 26453 81583 42807 xE-16850 53061 84548 59037 98217 99459 96811 54419 74253 63443 99602 xE-16900 90251 00158 88272 16474 50068 20704 19376 15845 47123 18346 xE-16950 00726 29339 55054 82395 57137 25684 02322 68213 01247 67945 xE-17000 22644 82091 02356 47752 72308 20810 63518 89915 26928 89108 xE-17050 45557 11266 03965 03439 78962 78250 01611 01532 35160 51965 xE-17100 59042 11844 94990 77899 92007 32947 69058 68577 87872 09829 xE-17150 01352 95661 39788 84860 50978 60859 57017 73129 81553 14951 xE-17200 68146 71769 59760 99421 00361 83559 13877 78176 98458 75810 xE-17250 44662 83998 80600 61622 98486 16935 33738 65787 73598 33616 xE-17300 13384 13385 36842 11978 93890 01852 95691 96780 45544 82858 xE-17350 48370 11709 67212 53533 87586 21582 31013 31038 77668 27211 xE-17400 57269 49518 17958 97546 93992 64219 79155 23385 76623 16762 xE-17450 75475 70354 69941 48929 04130 18638 61194 39196 28388 70543 xE-17500 67774 32242 76809 13236 54494 85366 76800 00010 65262 48547 xE-17550 30558 61598 99914 01707 69838 54831 88750 14293 89089 95068 xE-17600 54530 76511 68033 37322 26517 56622 07526 95179 14422 52808 xE-17650 16517 16677 66727 93035 48515 42040 23817 46089 23283 91703 xE-17700 27542 57508 67655 11785 93950 02793 38959 20576 68278 96776 xE-17750 44531 84040 41855 40104 35134 83895 31201 32637 83692 83580 xE-17800 82719 37831 26549 61745 99705 67450 71833 20650 34556 64403 xE-17850 44904 53627 56001 12501 84335 60736 12227 65949 27839 37064 xE-17900 78426 45676 33881 88075 65612 16896 05041 61139 03906 39601 xE-17950 62022 15368 49410 92605 38768 87148 37989 55999 91120 99164 xE-18000 64644 11918 56827 70045 74243 43402 16722 76445 58933 01277 xE-18050 81586 86952 50694 99364 61017 56850 60167 14535 43158 14801 xE-18100 05458 86056 45501 33203 75864 54858 40324 02987 17093 48091 xE-18150 05562 11671 54684 84778 03944 75697 98042 63180 99175 64228 xE-18200 09873 99876 69732 37695 73701 58080 68229 04599 21236 61689 xE-18250 02596 27304 30679 31653 11494 01764 73769 38735 14093 36183 xE-18300 32161 42802 14976 33991 89835 48487 56252 98752 42387 30775 xE-18350 59555 95546 51963 94401 82184 09984 12489 82623 67377 14672 xE-18400 26061 63364 32964 06335 72810 70788 75816 40438 14850 18841 xE-18450 14318 85988 27694 49011 93212 96827 15888 41338 69434 68285 xE-18500 90066 64080 63140 77757 72570 56307 29400 49294 03024 20498 xE-18550 41656 54797 36705 48558 04458 65720 22763 78404 66823 37985 xE-18600 28271 05784 31975 35417 95011 34727 36257 74080 21347 68260 xE-18650 45022 85157 97957 97647 46702 28409 99561 60156 91089 03845 xE-18700 82450 26792 65942 05550 39587 92298 18526 48007 06837 65041 xE-18750 83656 20945 55434 61351 34152 57006 59748 81916 34135 95567 xE-18800 19649 65403 21872 71602 64859 30490 39787 48958 90661 27250 xE-18850 79482 82769 38953 52175 36218 50796 29778 51461 88432 71922 xE-18900 32238 10158 74445 05286 65238 02253 28438 91375 27384 58923 xE-18950 84422 53547 26530 98171 57844 78342 15822 32702 06902 87232 xE-19000 33005 38621 63479 88509 46954 72004 79523 11201 50432 93226 xE-19050 62827 27632 17790 88400 87861 48022 14753 76578 10581 97022 xE-19100 26309 71749 50721 27248 47947 81695 72961 42365 85957 82090 xE-19150 83073 32335 60348 46531 87302 93026 65964 50137 18375 42889 xE-19200 75579 71449 92465 40386 81799 21389 34692 44741 98509 73346 xE-19250 26793 32107 26868 70768 06263 99193 61965 04409 95421 67627 xE-19300 84091 46698 56925 71507 43157 40793 80532 39252 39477 55744 xE-19350 15918 45821 56251 81921 55233 70960 74833 29234 92103 45146 xE-19400 26437 44980 55961 03307 99414 53477 84574 69999 21285 99999 xE-19450 39961 22816 15219 31488 87693 88022 28108 30019 86016 54941 xE-19500 65426 16968 58678 83726 09587 74567 61825 07275 99295 08931 xE-19550 80521 87292 46108 67639 95891 61458 55058 39727 42098 09097 xE-19600 81729 32393 01067 66386 82404 01113 04024 70073 50857 82872 xE-19650 46271 34946 36853 18154 69690 46696 86939 25472 51941 39929 xE-19700 14652 42385 77625 50047 48529 54768 14795 46700 70503 47999 xE-19750 58886 76950 16124 97228 20403 03995 46327 88306 95976 24936 xE-19800 15101 02436 55535 22306 90612 94938 85990 15734 66102 37122 xE-19850 35478 91129 25476 96176 00504 79749 28060 72126 80392 26911 xE-19900 02777 22610 25441 49221 57650 45081 20677 17357 12027 18024 xE-19950 29681 06203 77657 88371 66909 10941 80744 87814 04907 55178 xE-20000 20385 65390 99104 77594 1413 Finding a little more memory for slightly larger arrays, and printing 9/8 x 4 decimal digits per 16 bit word in the array, 20385 65390 99104 77594 14132 15432 84406 25030 18027 57169 xE-20050 65082 09642 73484 14695 72639 78842 56008 45312 14065 93580 xE-20100 90412 71135 92004 19759 85136 25479 61606 32288 73618 13673 xE-20150 73244 50607 92441 17639 97597 46193 83584 57491 59880 97667 xE-20200 44709 30065 46342 42346 06342 37474 66608 04317 01260 05205 xE-20250 59284 93695 94143 40814 68529 81505 39471 78900 45183 57551 xE-20300 54125 22359 05906 87264 87863 57525 41911 28887 73717 66374 xE-20350 86027 66063 49603 53679 47026 92322 97186 83277 17393 23619 xE-20400 20077 74522 12624 75186 98334 95151 01986 42698 87847 17193 xE-20450 96649 76907 08252 17423 36566 27259 28440 62043 02141 13719 xE-20500 92278 52699 84698 84770 23238 23840 05565 55178 89087 66136 xE-20550 01304 77098 43861 16870 52310 55314 91625 17283 73272 86760 xE-20600 07248 17298 76375 69816 33541 50746 08838 66364 06934 70437 xE-20650 20668 86512 75688 26614 97307 88657 01568 50169 18647 48854 xE-20700 16791 54596 50723 42877 30699 85371 39043 00266 53078 39877 xE-20750 63850 32381 82155 35597 32353 06860 43010 67576 08389 08627 xE-20800 04984 18885 95138 09103 04235 95782 49514 39885 90113 18583 xE-20850 58406 67472 37029 71497 85084 14585 30857 81339 15627 07603 xE-20900 56390 76394 73114 55495 83226 69457 02494 13983 16343 32378 xE-20950 97595 56808 56836 29725 38679 13275 05554 25244 91943 58912 xE-21000 84050 45226 95381 21791 31914 51350 09938 46311 77401 79715 xE-21050 12283 78546 01160 35955 40286 44059 02496 46693 07077 69055 xE-21100 48102 88502 08085 80087 81157 73817 19174 17760 17330 73855 xE-21150 47580 06056 01433 77432 99012 72867 72530 43182 51975 79167 xE-21200 92969 96504 14607 06645 71258 88346 97979 64293 16229 65520 xE-21250 16879 73000 35646 30457 93088 40327 48077 18115 55330 90988 xE-21300 70255 05207 68046 30346 08658 16539 48769 51960 04408 48206 xE-21350 59673 79473 16808 64156 45650 53004 98816 16490 57883 11543 xE-21400 45485 05266 00698 23093 15777 65003 78070 46612 64706 02145 xE-21450 75057 93270 96204 78256 15247 14591 89652 23608 39664 56241 xE-21500 05195 51052 23572 39739 51288 18164 05978 59142 79148 16542 xE-21550 63289 20042 81609 13693 77737 22299 98332 70820 82969 95573 xE-21600 77273 75667 61552 71139 22588 05520 18988 76201 14168 00546 xE-21650 87365 58063 34716 03734 29170 39079 86396 52296 13128 01782 xE-21700 67971 72898 22936 07028 80690 87768 66059 32527 46378 40539 xE-21750 76918 48082 04102 19447 19713 86925 60841 62451 12398 06201 xE-21800 13184 54124 47820 50110 79876 07171 55683 15407 88654 39041 xE-21850 21087 30324 02010 68534 19472 30476 66672 17498 69868 54707 xE-21900 67812 05124 73679 24791 93150 85644 47753 79853 79973 22344 xE-21950 56122 78584 32968 46647 51333 65736 92387 20146 47236 79427 xE-22000 87004 25032 55589 92688 43495 92876 12400 75587 56946 41370 xE-22050 56251 40011 79713 31662 07153 71543 60068 76477 31867 55871 xE-22100 48783 98908 10742 95309 41060 59694 43158 47753 97009 43988 xE-22150 39491 44323 53668 53920 99468 79645 06653 39857 38887 86614 xE-22200 76294 43414 01049 88899 31600 51207 67810 35886 11660 20296 xE-22250 11936 39682 13496 07501 11649 83278 56353 16145 16845 76956 xE-22300 87109 00299 97698 41263 26650 23477 16728 65737 85790 85746 xE-22350 64607 72283 41540 31144 15294 18804 78254 38761 77079 04300 xE-22400 01566 98677 67957 60909 96693 60755 94965 15273 63498 11896 xE-22450 41304 33116 62774 71233 88174 06037 31743 97054 06703 10967 xE-22500 67657 48695 35878 96700 31925 86625 94105 10533 58438 46560 xE-22550 23391 79674 92678 44763 70847 49783 33655 57900 73841 91473 xE-22600 19886 27135 25954 62518 16043 42253 72996 28632 67496 82405 xE-22650 80602 96421 14638 64368 64224 72488 72834 34170 44157 34824 xE-22700 81833 30164 05669 59668 86676 95634 91416 32842 64149 74533 xE-22750 34999 94800 02669 98758 88159 35073 57815 19588 99005 39512 xE-22800

Using a much much faster algorithm in 2019 - "Go For Broke"

"Go For Broke"

In memory of the 442nd Regimental Combat Team who fought valiantly against the Nazi in WWII in-spite of the fact that over 100,000 U.S. citizens of Japanese ancestry were absolutely screwed by Roosevelt. If you are in a "detention camp" with no income, you can't pay the taxes on your farm or home or dockage fees for your fishing boat - so the government takes 'em - and you lose "everything"!
Our church in Minnesota heard of this and provided names and camp addresses of "Japanese" Americans. My sister wrote to a young lady until she was released near the end of WWII.
I heard of a few Californians who helped their absent neighbors by paying their taxes and fees while their friends were "detained" in "camps".

On Saturday afternoon, Jan 5th 2019, Ron Mak, an Adjunct Professor at San Jose State University at San Jose State University e-mailed:

On my Windows PC with an Intel Core i7 CPU, I was able to calculate a billion decimal digits of pi in only 1.24 hours using the �PiFast� software from this website: http://numbers.computation.free.fr/Constants/PiProgram/pifast.html I thought it might take all night, but it didn�t. These are amazing advancements in computer hardware and in pi calculation formulas. The early mainframes took hours (if not days) to compute just a few hundred digits of pi using classic arctangent formulas. I zipped and uploaded my output file to http://www.cs.sjsu.edu/~mak/pi/pi-1G.txt.zip Here�s the head of the file: ~/~mak/pi: /usr/bin/head -25 pi-1G.txt Program : PiFast version 4.3 (fix 1), by Xavier Gourdon Computation of 1000000000 digits of Pi Method used : Chudnovsky Size of FFT : 32768 K Physical memory used : ~ 1577920 K Disk memory used : ~ 1907.35 Meg ------------------------------------------------------------ Computation run information : Start : Sat Jan 05 02:29:42 2019 End : Sat Jan 05 03:44:17 2019 Duration : 4475.09 seconds Time spent in disk swapping : 434 s (reading 21, writing 412) ============================================================ Total computation time : 4475.09 seconds (~ 1.24 hours) ============================================================ Pi with 1000000000 digits : Pi = 3. 1415926535 8979323846 2643383279 5028841971 6939937510 : 50 5820974944 5923078164 0628620899 8628034825 3421170679 : 100 8214808651 3282306647 0938446095 5058223172 5359408128 : 150 4811174502 8410270193 8521105559 6446229489 5493038196 : 200 4428810975 6659334461 2847564823 3786783165 2712019091 : 250 4564856692 3460348610 4543266482 1339360726 0249141273 : 300 And its tail: ~/~mak/pi: /usr/bin/tail -5 pi-1G.txt 8755703450 2943416861 5993243541 9973181435 5060392734 : 999999900 6434543524 2766553567 4357021939 6394581990 5483278746 : 999999950 7139868209 3196353628 2046127557 1517139511 5275045519 : 1000000000 How do I know this is correct? Here�s a file from MIT: http://www.cs.sjsu.edu/~mak/pi/pi-billion-mit.txt.zip The tail of the MIT file (there are no line breaks in the file): ~/~mak/pi: tail -r -c -50 pi-billion-mit.txt 71398682093196353628204612755715171395115275045519 Since the last 50 digits match, I�ll assume the rest of the digits match, too. The one millionth decimal digit of pi is a 1. The one billionth decimal digit is a 9. (I have much work to do today. But computing pi is a lot more fun!) � Ron http://www.cs.sjsu.edu/~mak/

I have always used John Machin's simple arctan formula from 1706 since it is very easy to code, and with a little help, I can prove it ;-)
But lets see what the modern world of fast formulas and complex algorithms has to offer. So I down loaded the PiFast code and docs and tried it out?
Using my old Pentium machine (10 years old?) my time was over 3 times longer than Ron's,
about 4.5 hours, rather than months of computing - IMPRESSIVE
Ron also sent the following documents which help explain the speed-up algorithms.
    - BinarySplittingMethod.pdf
    - FFT-BasedMultiplication.pdf
    - ArbitraryPrecisionComputation.pdf

A postscript from Ron:

Just a thought: Since the algorithm pages data out to disk, which is very slow compared to main memory, what if you created a RAM disk and told the algorithm to use that? Then how fast could it calculate pi?
Since I have 16 GB of main memory in my Windows laptop, I tried to tell the algorithm not to use the disk at all. But then the program crashed.
Apparently, the PiFast website is old (version 4.3 of the software was released in July 2003). Back then, the record was a billion digits in 3 hours 48 minutes. Maybe Ed�s machine is of that prehistoric era?

Info from Randy Thelen:

I was hoping to find the source code to the algorithm. But, the PiFast author has not released it. :-(

a 2020 adventure by Ronald Mak

Subject: Very fast pi calculation From: Ronald Mak/SJSU Date: Thu, August 06, 2020 1:08 am To: ...
One of the skills I teach my students is how to download and install a Linux software package that is distributed in source form. They must run the configure and make scripts to build and install the package. One of the packages I use is the Multiple Precision Integers and Rationals Library (MPIR) which enables a C or C++ program to perform integer and floating-point calculations with unlimited precision. See http://mpir.org/. We used to call this "string decimal arithmetic" - very much like how the IBM 1401 does arithmetic in hardware. MPIR is highly optimized and is distributed as x86-64 assembly code.
An assignment I often give is to use MPIR to implement the Borwein pi formula.
See https://en.wikipedia.org/wiki/Borwein%27s_algorithm - scroll down to "Quartic convergence (1985)". This formula is easy to code and iteratively computes the value 1/pi, and each iteration increases the number of correct decimal digits by a factor of 4. After 10 iterations, you will have one million correct digits, which you then invert to get the first million digits of pi.
On my MacBook Pro laptop (2.4 GHz 8-core Intel Core i9) running Ubuntu 19.04, I can compute a million digits of pi using the Borwein formula and MPIR's C library in 2.364 seconds, which does not include the printing time. See the attached C program BigPi.c and its output BigPi.txt. How do I know the output is correct? I compared each and every digit to what Dr. Evil computed -
see http://3.141592653589793238462643383279502884197169399375105820974944592.com/index314159.html . Well, at least the last block of 10 digits match: 5779458151
Isn't it amazing how fast our computers have become? In 1949, it took the ENIAC 70 hours to compute 2,035 digits of pi.
But wait, it gets even better! Like most modern computers, my laptop is multicore. How can I exploit this technology? I also teach an intro to multithreaded programming. In fact, I warn my students that if all they know is single-threaded programming, they will soon become obsolete.
If you look carefully at Borwein's quartic formula, you can see that besides the main thread, you can have a second thread that computes the y values and a third thread that computes the powers of 2. Computing the powers of 2 is fast, and that thread does its job (storing the powers into an array) quickly and dies. Then in each iteration of the formula, in order to compute the next value of a, the main thread must wait for the y thread to produce the next value of y. But as the main thread is computing an a value, the y thread can concurrently compute the next y value. Therefore, computing each iteration's a value overlaps with computing the y value for the next iteration. I use semaphores and mutexes from the pthreads library to synchronize the operations of these threads.
See the attached multithreaded C program BigPiMT.c and its output BigPiMT.txt. The multithreaded C program computes one million digits of pi in 1.994 seconds. The text file BigPiTimings.txt is an easy comparison. Despite the overhead of managing the multithreading, the MT version is over 15% faster.
Enjoy! It's past time for this thread to go to sleep ...
-- Ron

Ronald Mak http://www.cs.sjsu.edu/~mak/ Department of Computer Engineering Department of Computer Science Department of Applied Data Science San Jose State University San Jose, CA USA 95192 San Jose State University Research Foundation NASA Ames Research Center

"On two occasions I have been asked, 'Pray, Mr. Babbage, if you put into the machine wrong figures, will the right answers come out?' I am not able rightly to apprehend the kind of confusion of ideas that could provoke such a question."
-- Charles Babbage, 1791-1871, World's first computer scientist
"The purpose of computing is insight, not numbers."
-- Richard Hamming, 1915-1998, Mathematician
"We choose to go to the moon in this decade and do the other things, not because they are easy, but because they are hard ..."
-- President John F. Kennedy, September 12, 1962

Attachments
BigPi.c
BigPi.txt
BigPiMT.c
BigPiMT.txt
BigPiTimings.txt

Followed a day later by

I made a simple improvement BigPiMT-2.c to my multithreaded pi program. The y thread is the bottleneck. During each iteration of the formula, it takes more time than the main thread. So to balance the workload, I moved calculating y^2 from the end of the y thread into the main thread. My timings now for a million digits of pi:

Single threaded: 2.364 seconds
Original multithreaded: 1.994 seconds
Rebalanced multithreaded: 1.899 seconds
Nearly another 5% improvement in time.
Van suggested looking into Coarray Fortran which has high-level support for parallel processing.

a 2020 adventure by Ed Thelen
Using an Apple iMac with 3GHz 6-Core Intel Core i5
- A "C" implementation of Machin's algorithm
- Timings for my Intel processor
- Much abbreviated output

Well its December 11, 2020 - Pi again - fast Intel processor, C language
My son Randy works at Apple and had something to do with the Apple "M1" processor announced November 10, 2020.
This "M1" processor is apparently lightning fast but uses much less electric power than the Intel processor in my machine. Apple pushes "performance per watt".
I thought it would be fun to compare the Intel and Apple processors running the Machin Pi algorithm written in identical C on both machines. See https://en.wikipedia.org/wiki/Machin-like_formula for discussion and proof of Machin's method.
Randy is gearing up to use the much more complicated Ramanujan formula, also here, which requires multiple length multipliers, divisors, and also extreme precision powers and square roots. I am much too lazy (and incompetent?) to get so complex :-(

A "C" implementation of Machin's algorithm

// Created by Ed Thelen on 11/13/20. // Randy reminded me of *pointers ;-) #include #include // for time data format // // need to run on a 64 bit machine so "long" is 64 bits, not 32 // // 2^16 = 65,536 // 2^32 = 4,294,967,296 (2^32)-1 = hex$ FFFFFFFF // 2^64 = 18,446,744,073,709,551,616 /* MyPi x */ // // game plan - Machin formula for Pi // Pi/4 = 4arctan(1/5)-arctan(1/239) // pi = 16arctan(1/5)-4arctan(1/239) // 5^2 = 25 239^2 = 57121 2^16 = 65,536 - divisors fit into 16 bits // a) make 3 large areas, atan_d1[], atan_d2[], atan_acc[] // area atan_d1[] is reused as an output field // b) make subroutines div( *fi, *fo, divisor, start, len ), // int add( *fi, *fo, start, len ), // sub( *fi, *fo, start, len ), // atan( divisor, initial_add_sub ), // status is some error - // zerof( *fi, len ), // // lets try to use 31 bit fields in a 64 bit word // #define f_len 54010 // actual field length - 2, 500,000 dec dig, allow two guard ints //#define f_len 108020 // actual field length - 2, 1,000,000 dec dig, allow two guard ints //#define f_len 216040 // actual field length - 2, 2,000,000 dec dig, allow two guard ints //#define f_len 240 // quick test #define DataSize 31 // size of data in each long long word (64 bits) #define DataMask 0x7fffffff // 32/4 - 8 characters static unsigned long long max_ndiv; // I forgot what this is for ????? static int add_sub_error; // quit if some form of add or subtract error static unsigned long long atan_d1[f_len+2], atan_d2[f_len+2], atan_acc[f_len+2]; // guard at end // note: atan_d1 is re-used as the decimal output field // -------------------------- void zerof( unsigned long long *fi, int lim ){ // zero designated field int counter = 0; while ( counter < lim) { fi[counter++] = 0x00000000; // 64 bits of zero //++counter; } } // -------------------------- // add fi into fo void add( unsigned long long *fi, unsigned long long *fo, int zeros, int len ){ unsigned long long carry = 0; // yes, initialize carry int index = len; // start at low end unsigned long long t; while ( (index >= 0) || (carry != 0) ) { t = carry + fi[index] + fo[index]; fo[index--] = t & DataMask; //printf(" carry = \n"); carry = t >> DataSize; //--index; } if ( index < -1 ){ add_sub_error = 1; } // bit 1 is add overflow //printf(" exit add eflag=%X \n", add_sub_error); } // -------------------------- // sub fi from fo // error if return not 0 void sub( unsigned long long *fi, unsigned long long *fo, int zeros, int len ){ unsigned long long borrow = 0; // yes, initialize borrow int index = len; // start at low end unsigned long long t; //printf(" entering sub, index = %d \n", index); while ( index >= 0 ) { //printf(" click %d \n", index); t = fo[index] - borrow - fi[index]; fo[index--] = t & DataMask; borrow = t >> DataSize & 1; //--index; } if ( index < -1 ) add_sub_error |= 2; // bit 2 is subtract underflow //printf(" exit sub eflag=%d \n", add_sub_error); } // -------------------------- // multiply in place by multiplier, return is int part[0] void mpy( unsigned long long *fio, int start, unsigned long long multiplier ){ unsigned long long carry = 0; // hi part of product int index = start; // start at low end unsigned long long t; while ( index >= 0 ) { // t = carry + fi[index] + fo[index]; t = ( fio[index] * multiplier) + carry; fio[index--] = t & DataMask; carry = t >> DataSize; //--index; } } // -------------------------- // divide field by u long long void div( unsigned long long *fi, unsigned long long *fo, unsigned long long divisor, int start, int len ); // "fi" is input array, "fo" is output array, may be same area, // "start" is where to start dividing, to save time because high order becomes all zeros // "len" is total length of working array void div( unsigned long long *fi, unsigned long long *fo, unsigned long long divisor, int start, int len ){ int index = start; // get to register int endi = len; unsigned long long remainder = 0; // storage for remainder unsigned long long t = 0; // storage for quotient unsigned long long dividend = 0; // dividend /// ( a 64 bit register, with register zeroing rules ) //printf(" entering div \n"); while (index <= endi) { dividend = fi[index]; // low 31 bits of dividend, //hi order zeroed because it is a register dividend |= remainder << DataSize; // form hi order of dividend (register) t = dividend/divisor; // divide to get quotient fo[index++] = t & DataMask; // quotient, send to array of 16 bit unsigned ints remainder = dividend - (t*divisor); // form remainder //printf("%llu ",fo[index]); //++index; // go hi order to low order } //printf(" exit div\n"); } // -------------------------- // int myatan( divisor, initial_add_sub ), // does NOT clear accumulator !! // division by 239^2 works just fine void myatan( int recip_radians, int sub1st_flag ){ unsigned long long loc_divisor = recip_radians*recip_radians ; unsigned long long loc_flag = sub1st_flag; unsigned long long all0s = 0; // initialize atan_d1 field zerof( atan_d1, f_len); // zero most working fields zerof( atan_d2, f_len); atan_d1[0] = 1; // prepare for recip radians div( atan_d1, atan_d1, recip_radians, all0s, f_len ); // d1 correct if ( (loc_flag&1) == 0) add( atan_d1, atan_acc, all0s, f_len ); // add else sub( atan_d1, atan_acc, all0s, f_len ); // subtract loc_flag ^= 1; // flip flag // fields correctly set for taylor series unsigned int divn = 3; // initial t2 divisor // divide d1 by divisor while (all0s != (f_len)){ //void div( unsigned long *fi, unsigned long *fo, // unsigned int divisor, int start, int len ); div( atan_d1, atan_d1, loc_divisor, all0s, f_len ); // divide d1 giving d1 div( atan_d1, atan_d2, divn, all0s, f_len ); // divide d1 giving d2 divn += 2; // update, prep for next round if (divn > max_ndiv) max_ndiv = divn; // remember max size //printf("I'm here 1 \n"); if ( (loc_flag & 1) == 0) add( atan_d2, atan_acc, all0s, f_len ); // add else sub( atan_d2, atan_acc, all0s, f_len ); // subtract loc_flag ^= 1; // flip flag if ( atan_d1[all0s] == 0) ++all0s; // start at next word next time } } int main(int argc, const char * argv[]) { // insert code here... printf("Hello, World! - Lets Try MyPi again\n"); int i; i = sizeof(unsigned long long); printf(" machine size = %d 8 bit char\n", i ); if (i != 8) { printf(" machine not 64 bit type, FATAL ERROR \n"); return (-1); } printf(" machine is 64 bit type, \n"); printf(" Starting, field length = %d \n", f_len); // time variables time_t start, stop; start = time(NULL); stop = start; max_ndiv = 0; // clear warns and errors add_sub_error = 0; zerof( atan_acc, f_len); // clear result myatan( 5, 0 ); // form arctan (1/5) mpy(atan_acc, f_len, 4 ); // follow Machin start = time(NULL); printf(" end 4*arctan(1/5), elapsed time = %ld \n", start-stop ); stop = start; myatan( 239, 1 ); // form arctan (1/239) mpy(atan_acc, f_len, 4 ); // follow Machin start = time(NULL); printf(" end 4*arctan(239) + 16*arctan(1/5), elapsed time = %ld \n", start-stop ); stop = start; if ( add_sub_error != 0){ printf(" ERROR - add_sub_error = %d \n",add_sub_error); return -1; } unsigned long long wdigits; // start print out int line; int u1; int exp; exp = 0; printf("starting output printing \n"); int num_line = (f_len/5.4)-1 ; printf(" number of lines = %d\n", num_line); printf(" "); printf(" %llu E 0\n", atan_acc[0]); for (line=0; line < num_line; ++line){ for (u1=0; u1 < 10; ++u1){ atan_acc[0] = 0; mpy(atan_acc, f_len-(line*4), 100000 ); wdigits = atan_acc[0]; // to be printted if (wdigits < 10000) printf("0"); // now pad leading zeros if (wdigits < 1000 ) printf("0"); if (wdigits < 100 ) printf("0"); if (wdigits < 10 ) printf("0"); printf("%llu ", atan_acc[0]); } exp -= 50; printf(" E%6d\n", exp); } printf("OK, we are done! ;-) \n"); start = time(NULL); printf(" end output, elapsed time = %ld \n", start-stop ); return 0; }

Timings for my Intel processor

Field Length
64 bit units # Digits Out
Atan 1/5
(sec) Atan 1/289
(sec) Tot. Comp.
(sec) Output
(sec) Total
(sec) Total
(min)
54,010 500,000 303 88 391 9 400 6.67
108,020 1,000,000 1,208 352 1,560 39 1,599 26.65
216,040 2,000,000 4,837 1,439 6,276 155 6,431 107.18

Here is a report on optimization of Machin Pi by unrolling loops. - I chose to repeat stuff in loop by repeating code 5 times before doing the while type branch - add, subtract, and multiply went easily, divide would hang the program. I screwed around for too many hours trying to unroll the divide, and quit. As sour grapes, the divide is bulkier and would not gain as high a percent improvement. OK - timing - after all that fuss there was about a 3+ percent improvement in speed of the arctangents. Maybe the branch prediction logic works too well to make the unrolling optimization worth the effort .
Experiment using Stormer's related formula
Pi/4 = 44arctan(1/57) + 7arctan(1/239) - 12arctan(682) + 24arctan(1/12943)

Field Length
64 bit units # Digits Out
Tot. Comp.
(sec) Total
(sec) Total
(min)
54,010 500,000 326 336 5.6
Stormer's Total Compute time is 16% faster -

2,000,000 Digits of Pi - Much abbreviated

ed@Eds-iMac Mac-C % clang MyPi.c -o MiPy.o ed@Eds-iMac Mac-C % ./MiPy.o Hello, World! - Lets Try MyPi again machine size = 8 8 bit char machine is 64 bit type, Starting, field length = 216040 end 4*arctan(1/5), elapsed time = 4837 end 4*arctan(239) + 16*arctan(1/5), elapsed time = 1439 starting output printing number of lines = 40006 3 E 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 E -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 E -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 E -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 E -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 E -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 E -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 E -350 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 E -400 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 E -450 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 E -500 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 E -550 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132 E -600 00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 E -650 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 E -700 42019 95611 21290 21960 86403 44181 59813 62977 47713 09960 E -750 51870 72113 49999 99837 29780 49951 05973 17328 16096 31859 E -800 50244 59455 34690 83026 42522 30825 33446 85035 26193 11881 E -850 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 E -900 59825 34904 28755 46873 11595 62863 88235 37875 93751 95778 E -950 18577 80532 17122 68066 13001 92787 66111 95909 21642 01989 E -1000 38095 25720 10654 85863 27886 59361 53381 82796 82303 01952 E -1050 03530 18529 68995 77362 25994 13891 24972 17752 83479 13151 E -1100 and on and on ... 21624 84391 40359 98953 53945 90944 07046 91209 14093 87001 E -9900 26456 00162 37428 80210 92764 57931 06579 22955 24988 72758 E -9950 46101 26483 69998 92256 95968 81592 05600 10165 52563 75678 E-10000 56672 27966 19885 78279 48488 55834 39751 87445 45512 96563 E-10050 44348 03966 42055 79829 36804 35220 27709 84294 23253 30225 E-10100 76341 80703 94769 94159 79159 45300 69752 14829 33665 55661 E-10150 and more and more ... 30600 54697 16387 94317 11968 73484 68873 81866 56751 27929 E-99900 85750 16363 41131 46275 30499 01913 56468 23804 32997 06957 E-99950 70150 78933 77286 58035 71279 09137 67420 80565 54936 24646 E-100000 41260 02437 96845 43777 33902 64725 12819 41632 00768 48736 E-100050 25176 40659 67540 69362 17588 79307 85591 64787 77274 73927 E-100100 20029 10342 94956 24476 61308 20072 92507 34529 17076 42266 E-100150 will this ever stop ??? 99745 65987 36982 87428 26826 74844 21213 71546 82327 51128 E-499900 53416 52370 31653 07043 25828 33721 12371 37609 69593 99375 E-499950 49536 22322 22197 46596 19332 52907 40424 87602 51381 95242 E-500000 69739 10175 63719 75343 00447 96178 25043 11533 15067 58256 E-500050 27353 43476 25253 91425 15275 70478 78437 67885 24187 19663 E-500100 Good Grief Charlie Brown !!! 64874 79017 72741 25678 19055 55621 80504 87674 69911 40839 E-999850 97791 93765 42320 62337 47173 24703 36976 33579 25891 51526 E-999900 03156 14033 32127 28491 94418 43715 06965 52087 54245 05989 E-999950 56787 96130 33116 46283 99634 64604 22090 10610 57794 58151 E-1000000 30927 56283 20845 31584 65200 10277 97235 61292 30126 05863 E-1000050 53601 16492 09902 58745 55214 03969 79115 34022 41589 81324 E-1000100 50572 29232 33213 10453 92469 98366 11275 09161 34187 79722 E-1000150 OK, you win, I'm exhausted ! 53675 96380 19091 94175 86559 31287 39602 79125 10596 54044 E-1999900 96212 15177 02095 78971 06655 25923 69719 33822 82267 49132 E-1999950 29071 74473 58925 65046 16637 35632 36871 06519 14572 97909 E-2000000 61217 31251 17978 46672 66882 73303 91053 48653 07389 01742 E-2000050 82734 37618 99789 18962 35310 76701 58514 24532 85425 10122 E-2000100 76696 94946 58077 27822 87268 83325 84497 96975 22136 06776 E-2000150 OK, we are done! ;-) end output, elapsed time = 155

Dream 2020 - Using a cell phone ;-) with fast 64 bit ARM processor and much more memory.

My sons and I tease each other "a lot".
Mentioning that Apple cellphones use 64 bit ARM processors, son Randy sent the following:

Looking at this web page: https://www.johndcook.com/blog/2011/03/14/algorithm-record-pi-calculation/ I found the Ramanujan algorithm for computing Pi: Here's a program without multi-precision arithmetic: A = 6 - 4 * math.sqrt(2) Y = math.sqrt(2) - 1 def iter(a, y, n): y_sub = math.pow(1 - math.pow(y, 4), 1/4) y_num = 1 - y_sub y_den = 1 + y_sub y_next = y_num / y_den a_part1 = math.pow(1 + y_next, 4) * a a_part2 = math.pow(2, 2 * n + 3) * y_next a_part3 = 1 + y_next + math.pow(y_next, 2) a_next = a_part1 - a_part2 * a_part3 return a_next, y_next for n in range(0, 5): print(f"n = {n}, pi = {1/A}") A, Y = iter(A, Y, n) print(f"pi = {1/A}") And its output: n = 0, pi = 2.9142135623730985 n = 1, pi = 3.1415926462135473 n = 2, pi = 3.1415926535897927 n = 3, pi = 3.1415926535897927 n = 4, pi = 3.1415926535897927 pi = 3.1415926535897927

It would be easy to imply that Randy was restless with the ancient Machin algorithm I was using.
And he seems ready to take on extended precision multiply and divide using extended length multipliers and divisors ??
(Using Machin I had avoided this added complexity.)
Randy is my youngest, age mid-50s, apparently in his prime. I'm 89, and NOT in my prime :-((
So, what to do - look on the Internet of course ;-))
- Starting FORTH, on-line
- - Randy is good at FORTH, which is light weight way to get something going on new computers ;-)
- "The Art of Computer Programming", Vol 2, Seminumerical Algorithms, Knuth, 3rd edition (I have it)
- - In extended precision work, divide is the most difficult - just ask any 4th grader
- - For your edification (I hope there is such a word) I offer pages Knuth-272-3 and Knuth-274-5
- - - (To avoid conflict, Knuth defined his own machine and mnemonics.)
- the Chudnovsky_algorithm based on Ramanujan
- Multiple Precision Arithmetic On Graphics Processing Units by Niall Emmart
- - Randy says using GPUs is another kettle of sharp toothed fish !!
A tale - late 1970s
I was working for Measurex, as was a mathematician named Ed Barlow. One day he visited my office (we didn't have cubicles - ala Dilbert) about something. I complained that I used Machin's formula to compute Pi but couldn't prove it. Off the top of his head, Ed wrote a proof on my chalk board!!
This derivation is similar to my memory of what was left on the board, for several weeks ;-)
Sunday, Nov 1, 2020, afternoon
Randy started coding from Knuth. I shouted encouraging words via the Internet. Then we noticed that Knuth and others talk of extended precision *integer* divide, and we think we want to deal with decimal numbers including the decimal fraction parts. Do we have to deal with scaling??

From Randy, evening Nov 2, 2020
Ah ha! It's alive!!
The digits computed by the app are listed below, but here's a chart that shows the relationship between the Number of "digits" in the Fixed Point number vs. the number of Pi digits computed (and are accurate):

The phrase "calc_pi(x)" indicates the number of iterations of the core algorithm that I'm running. It turns out 6 is the maximum I can handle reliably: the term 2^(2*n+3) where n == the iteration count being processed; creates a "fixed point" number that reaches my integer upper bounds at n = 6: 2x6+3 = 15.
N = 40; calc_pi(20) 00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 6446- ----- ----- ----- xE -200 N = 80; calc_pi(20) 00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 xE -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 xE -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 xE -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 xE -350 78925 90360 01133 05305 48820 46--- ----- ----- ----- ----- xE -400 N = 120; calc_pi(4) 00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 xE -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 xE -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 xE -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 xE -350 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 xE -400 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 xE -450 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 xE -500 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 xE -550 60943 70277 05392 1717- ----- ----- ----- ----- ----- ----- xE -600 N = 160; calc_pi(4) 00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 xE -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 xE -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 xE -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 xE -350 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 xE -400 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 xE -450 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 xE -500 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 xE -550 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132 xE -600 00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 xE -650 14684 40901 22495 34301 46549 58537 10507 92279 689-- ----- xE -700 N = 200; calc_pi(5) 00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 xE -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 xE -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 xE -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 xE -350 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 xE -400 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 xE -450 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 xE -500 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 xE -550 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132 xE -600 00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 xE -650 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 xE -700 42019 95611 21290 21960 86403 44181 59813 62977 47713 09960 xE -750 51870 72113 49999 99837 29780 49951 05973 17328 16096 31859 xE -800 50244 59455 34690 83026 42522 30825 33446 85035 26193 11881 xE -850 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 xE -900 59825 34904 28755 46873 11595 62863 88235 37875 93751 95778 xE -950 1857- ----- ----- ----- ----- ----- ----- ----- ----- ----- xE -1000 Just to show you I'm not making this up, here are the last few lines from the N = 200 output without the hyphens: 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 xE -900 59825 34904 28755 46873 11595 62863 88235 37875 93751 95778 xE -950 18575 29773 50061 87223 87560 61244 28383 73698 61150 64969 xE -1000 57764 42430 73331 97085 12995 73206 29583 66305 46466 77572 xE -1050 80590 51107 53656 85989 88633 76933 67506 82472 65886 57965 xE -1100 73962 04075 30754 68915 58874 48156 21447 12363 16938 18476 xE -1150 20704 43841 47694 60868 31500 14925 67759 18110 92766 69837 xE -1200 48644 39354 40404 38202 70356 12211 49138 19588 65922 90197 xE -1250 79741 98156 55687 70390 94441 83079 84413 99622 33603 01562 xE -1300 The data on the line labeled "xE -1000" shows the data going bad on the 5th character of the line. :-) And now, Pi to 1000 places: N = 210; calc_pi(5) 00003 xE 0 14159 26535 89793 23846 26433 83279 50288 41971 69399 37510 xE -50 58209 74944 59230 78164 06286 20899 86280 34825 34211 70679 xE -100 82148 08651 32823 06647 09384 46095 50582 23172 53594 08128 xE -150 48111 74502 84102 70193 85211 05559 64462 29489 54930 38196 xE -200 44288 10975 66593 34461 28475 64823 37867 83165 27120 19091 xE -250 45648 56692 34603 48610 45432 66482 13393 60726 02491 41273 xE -300 72458 70066 06315 58817 48815 20920 96282 92540 91715 36436 xE -350 78925 90360 01133 05305 48820 46652 13841 46951 94151 16094 xE -400 33057 27036 57595 91953 09218 61173 81932 61179 31051 18548 xE -450 07446 23799 62749 56735 18857 52724 89122 79381 83011 94912 xE -500 98336 73362 44065 66430 86021 39494 63952 24737 19070 21798 xE -550 60943 70277 05392 17176 29317 67523 84674 81846 76694 05132 xE -600 00056 81271 45263 56082 77857 71342 75778 96091 73637 17872 xE -650 14684 40901 22495 34301 46549 58537 10507 92279 68925 89235 xE -700 42019 95611 21290 21960 86403 44181 59813 62977 47713 09960 xE -750 51870 72113 49999 99837 29780 49951 05973 17328 16096 31859 xE -800 50244 59455 34690 83026 42522 30825 33446 85035 26193 11881 xE -850 71010 00313 78387 52886 58753 32083 81420 61717 76691 47303 xE -900 59825 34904 28755 46873 11595 62863 88235 37875 93751 95778 xE -950 18577 80532 17122 68066 13001 92787 66111 95909 21642 01989 xE -1000 38--- 34294 21374 96925 31925 09955 18466 70195 33374 98796 xE -1050 - Randy
I have plans to extend the fixed point data representation to include arbitrarily large integer portions; but there's also another way to solve this, but it's a bit delicate: the 2^(2*n+3) is then multiplied by Yn+1; and that means it can be distributed as part 2^(x) and part a bit shift of Yn+1. So, there are options that are available.
The punch line is that it takes 210 sixteen bit digits to compute Pi to 1002 places!

Outside Resources - a note, apparently both names "Multiple Precision" and "Arbitrary Precision" are acceptable.

Multiple precision divide - Knuth to the rescue

Knuth V2 pages 272&3
Knuth V2 pages 274&5

"Simple Algorithm for Arbitrary-Precision Integer Division" on Youtube
Host also suggests going to http://justinparrtech.com/JustinParr-Tech/an-algorithm-for-arbitrary-precision-integer-division/
or go to http://justinparrtech.com and search for "Integer Division". (I could not find a search box or button.)
Also https://www.youtube.com/watch?v=eCaXlAaN2uE, 11. Integer Arithmetic, Karatsuba Multiplication, skip the first 30 minutes - Intel based