I saw a program that calculates π with only 7 lines, and I find it incredibly fascinating Can anyone explain it?

Algorithm and Principles of Calculating Pi

I will spend a day to carefully answer this question, although I have just understood the principles behind it.

Undeniably, the techniques used here are worth learning. The key is to calculate division using only integers and without using large numbers. It is based on basic number theory (of course, this also limits the maximum number of decimal places of π that can be calculated). To be honest, if this program were put into a project, the author would be beaten to death. Meaningless variable names, program logic crammed inside for loops, and a lot of increment and decrement operators. What’s the point of pursuing fewer lines of code in this way?

Let’s analyze the knowledge points used in this algorithm.

I. Formulas for calculating π used in this algorithm

As we all know, there are many series formulas for calculating π. One famous formula is as follows:

This is the famous Leibniz formula. However, this formula converges very slowly. The formula used in this algorithm is as follows:

This formula seems to be derived from the previous formula, but it actually isn’t. It is derived from the Machin-like formula discovered by the Japanese mathematician Yasuaki Hiyoshi. This formula converges faster than the previous one. Later, I will discuss the convergence rate of this formula.

But, writing the program completely according to this formula is not easy and the loop is quite complex. Additionally, the intermediate calculation results to be saved quickly exceed the range that can be expressed by the long data type. So, the author of the program actually modified the program in a certain way (it took me a long time to figure it out while looking at the author’s program), as follows:

In this way, it is convenient to calculate the loop from back to front, starting from the nth term. The denominator of the nth term is 2n-1. The maximum value of the long data type is 2147483647, and in each step of the program, the denominator is multiplied by an enlarged 1 (multiplied by 10000) and multiplied by the numerator. Therefore, it should satisfy 10000 * (2n-1) < 2147483647, from which we can deduce that n < 107374.

In other words, using the author’s method, in the case of using only the long data type for calculation, it is possible to calculate up to the 107374th term of this series. Then, to what decimal place can we calculate π by calculating up to this term?

Now let’s discuss the convergence rate of this series.

II. Convergence rate of this series

Let’s take a look at this formula again:

Where , , , and each term is less than . Although they approach as they progress, when they are multiplied together to reach the nth term of this series, they must be less than , that is, . So, after calculating this, the number of zeros after the decimal point indicates how many digits of π can be accurately calculated.

For example, when calculating up to the 11th term of this series, everyone knows that . So, until the fourth digit after the decimal point, the number 9 appears. This means that it does not affect the first three digits after the decimal point of the value of π (at least it does not affect the first two digits, while the third digit may be affected by an addition carry). In other words, when we calculate up to the 10th term of this series, we can accurately calculate the first 3 significant digits of π.

This is the convergence rate of this series.

This is actually a relationship. If we calculate , we can get .

When k is equal to 4, n is approximately 14. In other words, by calculating 14 more terms of this series, we can obtain 4 more significant digits of π.

This is the relationship between 4 and 14 in this program.

By using this algorithm and performing calculations using the long data type, it is possible to calculate up to the 30676th digit of π.

III. Method for calculating integer division using integers

If single-precision floating-point numbers are used, the number of significant digits is 6-7, while for double-precision floating-point numbers, the number of significant digits is 15-17. Neither of them can be used to store the value of π that we need to calculate here. So, we can only use integers to represent it.

The author of this program uses the technique of bitwise shift to calculate integer division and can perform calculations with an infinite number of digits. He does it like this: First, multiply the numerator by 10000, then divide it by the denominator to obtain the quotient and remainder. Then divide the quotient into the part larger than 10000 and the part smaller than 10000. Output the part larger than 10000 (which is equivalent to the integer part of the quotient before the numerator is multiplied by 10000 and divided by the denominator), retain the part smaller than 10000, and then enter the next iteration together with the remainder. In the next iteration, multiply the remainder from the previous round by 10000 as the numerator, divide it by the denominator, and obtain the quotient and remainder again. At this time, enlarge the part smaller than 10000 left over from the previous round by 10000, add it to the quotient of this round, and split the result into the part larger than 10000 and the part smaller than 10000. Output the part larger than 10000, retain the part smaller than 10000, and enter the next round together with the remainder.

Why is it 10000? Because in this way, every round outputs 4 digits.

Then, let’s write a program to calculate and try it out. I will calculate 30 rounds, which is 120 significant digits, although everyone knows that it is actually a repeating decimal. The first version of the program will output the quotient, remainder, the number before splitting into the part larger than 10000 and the part smaller than 10000, and the number after splitting for each round. The program code is as follows:

#include <stdio.h>
int main(int argc, char** argv){
    long numerator = 1;   // Numerator
    long denominator = 7; // Denominator
    long low_bits = 0; // The lower 4 digits of the quotient calculated in the previous round, stored for printing in this round
    for(int i=0; i<30; i++){
        printf("Round %d\n", i);
        long quotient = numerator * 10000 / denominator;  // Multiply the numerator by 10000 and calculate the quotient 
        long remainder = numerator * 10000 % denominator; // Multiply the numerator by 10000 and calculate the remainder
        printf("quotient = %ld\n", quotient);
        printf("remainder = %ld\n", remainder);
        // Multiply the lower 4 digits of the quotient calculated in the previous round by 10000 and add them to the quotient of this round for printing
        long preprint = low_bits * 10000 + quotient;
        printf("preprint = %ld\n", preprint);
        long high_bits = preprint / 10000;
        printf("high_bits = %ld\n", high_bits);
        // The remainder becomes the numerator for the next round, and low_bits is also passed into the next round
        numerator = remainder;
        low_bits = preprint % 10000;
        printf("low_bits = %ld\n", low_bits);
    }
    printf("\n");
    return 0;
}

Compile and run the program, and the running result is as shown in the figure:

Figure 1, Program code

Figure 2, Detailed process of each round

Then, remove all these output information and use printf("%.4ld", high_bits) to output the result. The program code is as follows:

#include <stdio.h>
int main(int argc, char** argv){
    long numerator = 1;   // Numerator
    long denominator = 7; // Denominator
    long low_bits = 0; // The lower 4 digits of the quotient calculated in the previous round, stored for printing in this round
    for(int i=0; i<30; i++){
        long quotient = numerator * 10000 / denominator;  // Multiply the numerator by 10000 and calculate the quotient 
        long remainder = numerator * 10000 % denominator; // Multiply the numerator by 10000 and calculate the remainder
        // Multiply the lower 4 digits of the quotient calculated in the previous round by 10000 and add them to the quotient for printing
        long preprint = low_bits * 10000 + quotient;
        long high_bits = preprint / 10000;
        printf("%.4ld", high_bits);
        // The remainder becomes the numerator for the next round, and low_bits is also passed into the next round
        numerator = remainder;
        low_bits = preprint % 10000;
    }
    printf("\n");
    return 0;
}

The running result of this program is as shown in the figure:

Figure 3, Program code that only outputs the final result

Figure 4, The final output result

As you can see, the effect is quite good. You can change the numerator and denominator to other values and see if the calculation is correct. In our output, the decimal point is not displayed, but the position of the decimal point only affects the magnitude of the number and does not affect its composition. With this method, you can calculate as many digits as you want.

IV. The Final Touch

Based on the knowledge learned from the previous analysis, let’s modify the program for calculating π into a readable form. The code is as follows:

#include <stdio.h>
int main(int argc, char** argv){
    long n = 28000;        // Calculate up to which term
    long remainder[28001]; // This array is used to store the remainder after each term in each round of calculation
    long low_bits = 0;
    // Initialize all elements in the array as 2000
    // Why is it 2? That's because the original formula calculates pi/2, and after multiplying by 2, it becomes pi
    // As for 2000, it enlarges pi by exactly 500 times, so that when printing, there are no leading 0s before 3141
    for(int i=0; i<=n; i++){
        remainder[i] = 2000;
    }
    // Start the outer loop. Each time the loop iterates, 14 terms of the series are reduced and 4 digits of pi are output
    while(n){
        long denominator = 2*n-1;
        long quotient = 0;
        long sum = 0;
        sum += remainder[n] * 10000;
        remainder[n] = sum % denominator;
        quotient = sum / denominator;
        // Start the inner loop, which starts from b == n-1 and iterates to b == 0
        long b = n-1;
        while(b){
            sum = quotient * b + remainder[b] * 10000;
            denominator = 2*b -1;
            // Save the remainder to the array and enter the next iteration, while the quotient enters the next inner loop
            remainder[b] = sum % denominator;
            quotient = sum / denominator;
            b = b - 1;
        }
        n = n - 14;
        // Output the high 4 digits, and the low 4 digits of the quotient enter the next iteration
        printf("%.4ld", low_bits + sum / 10000);
        low_bits = sum % 10000;
    }
    printf("\n");
    return 0;
}

Here are screenshots of the program code and the running result:

Figure 5, The code to calculate π within less than 50 lines

Figure 6, Calculated 8000 digits of π

The entire code, including blank lines and comments, does not exceed 50 lines. This method of calculating π only requires integer calculations using the long data type, and it can calculate more than 30,000 digits of π. It’s truly amazing.

V. Conclusion

I spent the whole day today to complete the answer to this question. (Summarizing the algorithm and debugging the program took quite some time.) I understand the principles behind this algorithm, and it is really rewarding. Today, I answered several beginner-level questions related to the C language, and the achievements added up cannot compare to the sense of accomplishment I gained from solving this problem.

I hope everyone will support and give it a thumbs-up.

Calculating π in one line of code

I can’t believe this question made it to the hot posts… I actually wrote a guide on how to calculate π with just one line of code for everyone’s reference:

Concise and Effective Methods for the Convergence of Pi

You can check the Wikipedia entry for pi, and you will find that there are countless infinite series and continued fractions that converge to pi. Almost every one of them can be implemented in just a few lines of code. Ramanujan alone came up with more than a dozen of them.

Iterative formulas for calculating π and error analysis

Using a few formulas borrowed from the Jing Mountain Knight, Excel can solve them in just one-third of a row.

π/4=1-1/3+1/5-1/7…

=LET(n,5000,a,SEQUENCE(n),b,IF(ISODD(a),1,-1),4*SUM(1/(2*a-1)*b))

After 5000 iterations, the accuracy is not particularly good.

π/2=1+1/3+1/3*2/5+1/3*2/5*3/7+…

=2*SUM(REDUCE(1,SEQUENCE(20),LAMBDA(x,y,VSTACK(x,TAKE(x,-1)*y/(2*y+1)))))

This formula has a good convergence speed, and the error is already relatively small after 20 iterations.

π/2=1+1/3(1+2/5(1+3/7(1+..n-1/2n-1(1+…))))

=LET(n,20,2*REDUCE(0,SEQUENCE(n,,n,-1),LAMBDA(x,y,1+y/(2*y+1)*x)))

The error when iterating 20 times with this formula is twice as large as the previous formula.

Simplified code for calculating π

When I saw this program for calculating π in just 7 lines, I was also amazed. The code is incredibly concise, but its principle is not complex. It only requires some basic mathematical knowledge and programming skills.

First, let’s analyze the code line by line to understand how it calculates π:

#include<stdio.h>
long a=10000,b=0,c=28000,d,e=0,f[28001],g;
main(){
　　for(;b-c;)
　　　　f[b++]=a/5;
　　for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
　　　　for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}

First, this code uses the standard input-output library stdio.h in C, and it defines some variables and an array. Let’s explain their meanings one by one:

• long a=10000: Defines a long integer variable a and initializes it as 10000. This variable plays an important role in subsequent calculations, and we will explain its purpose in detail later.
• long b=0: Defines a long integer variable b and initializes it as 0.
• long c=28000: Defines a long integer variable c and initializes it as 28000.
• long d,e=0: Defines two long integer variables d and e, and initializes e as 0.
• long f[28001]: Defines a long integer array f with a size of 28001.
• long g: Defines a long integer variable g.

Next, we enter the main() function, which is the entry point of the program in C.

• for(;b-c;) f[b++]=a/5;: This is a loop statement that initializes the first c-b elements of array f to a/5. The loop condition is b-c, meaning that the loop continues executing as long as b is not equal to c. In the loop body, f[b++]=a/5 assigns the value of a/5 to the b-th element of array f and increments the value of b by 1. Through this loop, we store a series of initial values in array f.

• for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a): This is another loop statement that calculates the value of π and outputs it. The loop condition is d=0, which initializes the value of d as 0. In the loop body, g=c*2 assigns the value of c multiplied by 2 to g, c-=14 subtracts 14 from c, printf("%.4d",e+d/a) outputs the value of e+d/a as a four-digit decimal number, and e=d%a assigns the remainder of d/a to e.

• for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);: This is a nested loop statement that performs a series of calculations on array f to obtain the value of π. In the loop body, d+=f[b]*a adds the value of f[b]*a to d. Since f[b] is initially set as 2000, d is initially set as 2000. Then, f[b]=d%g assigns the value of d%g to f[b], which is the remainder of d divided by g. Note that g is decreasing, so its value decreases after each calculation. Finally, d/=g-- divides d by g and assigns the quotient to d. Additionally, --b decreases the value of b. This loop continues executing until the loop condition becomes false.

Through the above analysis, we have a rough understanding of how this code works. Now, let me explain further with some specific examples.

Let’s assume that we set the value of c as 14, then run this code and see what result we will get.

Initially, the first c-b elements of array f are assigned as a/5, which is 10000/5 = 2000. The printf("%.4d",e+d/a) in the loop body outputs the value of e+d/a, which is the quotient of e divided by a, with four decimal places. Since the initial value of e is 0, the output result is 0.

Then, we enter the nested loop and perform a series of calculations. Because the initial value of g is c2, which is 142 = 28, the loop will iterate 28 times.

In each loop, we first calculate d+=f[b]*a, which means adding the value of f[b]*a to d. Since the initial value of f[b] is 2000, the initial value of d is 2000. Then, we assign the value of d%g to f[b], which is the remainder of d divided by g. Note that g is decreasing, so its value decreases after each calculation. Finally, we divide d by g and assign the quotient to d, and decrease the value of b by 1. Then, the next iteration begins.

Through this loop, we continuously accumulate, take the remainder, and perform division operations on d. The result of d is eventually the value of π. Then, we output the value of e+d/a with four decimal places, which represents the first four decimal places of π. In this example, the output result is 3.1415.

The reason why this code can calculate the value of π is that it uses an algorithm called the “Bailey-Borwein-Plouffe” formula, which can be used to calculate the hexadecimal representation of π. However, the detailed derivation process of this formula is quite complex and beyond the scope of this answer. You can search for related materials to learn more.

In summary, the reason why this code, with only 7 lines, can calculate the value of π is that it utilizes a series of mathematical calculations and programming techniques. Through operations on arrays and nested loops, it ultimately obtains the hexadecimal representation of π. When we set the value of c as a suitable multiple, we can obtain any arbitrary n digits of π.

Infinite Product Sum Formula

The principle has been explained in the original text:

The principle is as follows: 1 + 1/3 + 1/3*2/5 + 1/3*2/5*3/7 + 1/3*2/5*3/7*4/9 + …, continuously adding up, the result is pi/2. 1/3<1/2, 2/5<1/2, 3/7<1/2, so this number converges faster, and the convergent value is smaller than 1+1/2+1/4+1/8+…=2.

How did this formula come about? This is not a programming question, but a mathematical question. Check out the “Infinite Product” section in the following article:

Chudnovsky Algorithm for Calculating Pi

Ask GPT directly

Question:

#include<stdio.h>
long a=10000,b=0,c=28000,d,e=0,f[28001],g;
main(){
　　for(;b-c;)
　　　　f[b++]=a/5;
　　for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a)
　　　　for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b);
}
I am a primary school student, please explain to me why this piece of code can calculate Pi?

Answer:

This code uses a mathematical algorithm called the Chudnovsky algorithm to calculate an approximate value of pi. I will try to explain how this algorithm works in the simplest way possible.
First, we need to understand some basic mathematical concepts. Pi (π) is an irrational number, whose value is approximately 3.14159265358979323846... It is very important in mathematics and science.
The Chudnovsky algorithm is a method for calculating pi that uses an approximation formula of a series (an infinite sum) to compute each digit of π.
The loops and mathematical calculations in this code aim to calculate each term in the series and add them up gradually to approach the value of pi. Let's break down how the code works step by step:
for(;b-c;) f[b++]=a/5; This loop is used to initialize the values of the array f. It stores the result of a/5 in the array f to prepare for the subsequent calculations.
for(;d=0,g=c*2;c-=14,printf("%.4d",e+d/a),e=d%a) This loop is used to calculate and output each digit of pi. Here's how this loop works:
d = 0 Initialize the variable d to 0.
g = c * 2 Initialize the variable g to twice the value of c.
c -= 14 Decrease the value of c by 14 in each iteration of the loop.
printf("%.4d", e + d / a) Print the value of the current digit of pi, where e + d / a represents the integer part, and %.4d specifies that the output should be a four-digit integer, padding with zeros if necessary.
e = d % a Update the variable e to the remainder of d divided by a.
for(b=c;d+=f[b]*a,f[b]=d%--g,d/=g--,--b;d*=b); This nested loop is used to calculate the value of each digit of pi. Here's how this loop works:
b = c Initialize the variable b to the value of c.
Inside the loop, the following operations are executed sequentially:
d += f[b] * a Add f[b] times a to d.
f[b] = d % --g Update f[b] to the remainder of d divided by g decreased by one, and decrement g by one.
d /= g-- Divide d by g, and then decrement g by one.
--b Decrement b by one.
This nested loop continues to execute until the condition in the loop body becomes false (i.e., the expression in the loop body equals 0).
Through these loops and mathematical calculations, the code gradually calculates each digit of pi and prints them in the form of a four-digit integer.
It is important to note that the value given by this code is an approximation of pi, not its exact value. Pi is an infinitely non-recurring decimal that cannot be represented precisely as a finite decimal or fraction. Therefore, we use algorithms to calculate its approximate value.
Hopefully, this explanation helps you understand how this code calculates pi. If you have any further questions, feel free to ask!

Program for calculating π using infinite series

As an introduction, this program uses a method called “infinite series” to calculate π. Specifically, it seems to be based on the Gregory-Leibniz series or a similar series, which is represented as follows:

Although this series is simple, it converges very slowly, which means that to obtain an accurate value of π, many terms need to be calculated. Let’s expand the code to make it more intuitive. In terms of the mathematical aspect, there are other answerers who specialize in it, so I will expand the code first.

#include<stdio.h>
long a = 10000, b = 0, c = 28000, d, e = 0, f[28001], g;

int main() {
    // Initialize array f to store intermediate calculation results
    for(; b - c;) {
        f[b++] = a / 5;
    }

    // Main loop to calculate each digit of π
    for(; d = 0, g = c * 2; c -= 14, printf("%.4d", e + d / a), e = d % a) {
        // Internal loop to perform actual iterative calculation
        for(b = c; b > 0; b--) {
            d += f[b] * a;
            f[b] = d % --g;
            d /= g--;
            // The statements f[b] = d % --g and d /= g-- in this part are decreasing steps, gradually reducing the denominator, simulating the decreasing behavior of the infinite series.
        }
    }
    return 0;
}

This loop initializes each element of the array f[] to a / 5. The index b gradually increases until it equals c. At the beginning of each outer loop iteration, d is set to 0 and g is set to c * 2. With each iteration, c is decreased by 14 digits, allowing the calculation of 14 digits of π each time.

The printf statement is used to output the current part of the calculation result. It outputs e + d / a, which is the result of the current iteration calculation. e = d % a is used to store the remainder in the calculation for the next iteration.

The internal iterative calculation of π is the core of the iterative calculation. It starts from b = c and gradually decreases b until it reaches 0. Inside the loop, d += f[b] * a is used to accumulate the calculation result. f[b] = d % –g updates the element of the array f[], which is used for the next round of calculation. The key step of the iterative calculation is d /= g–, which approximates the value of π through long division.

Overall, this program for calculating π is quite impressive in its functionality.

Convergence Description

For each increment of K, π seems to converge to 4 digits.

Ramanujan is simply a madman.

Proof of the Relevant Formulas

For

it is easy to see that

However, the two sides are obviously not equal. At this point, it is natural to consider another formula

where . By taking , we can obtain

which simplifies to

Thus, the proof is complete.

As for the proof of (*), we only need some basic calculus knowledge. By rearranging and differentiating both sides, we have

Further differentiating both sides and using the Leibniz rule, we have

Simplifying the equation, we get

Since , we can conclude that . Because , we have . Therefore, .

Furthermore, by further integration, we can obtain

This article is created and published using