Interpreting the legendary super-cool C program that calculates π

When I went to college, I circulated such a super-good C program. I could use only three lines of code to calculate the number of π to 800 digits after the decimal point. I joked that it was written by aliens, it's really amazing. At that time, everyone did not understand this article. I read an article yesterday to explain this code. Today I tried it for a long time and finally understood it. So I wrote it down and communicated with you.
This section of C code is as follows:
# Include "stdio. H"
Long A = 10000, B, c = 2800, d, e, f [2801], G;
Void main (){
For (; B-c;) f [B ++] = A/5;
For (; D = 0, G = C * 2; C-= 14, printf ("%. 4D", e + D/A), E = D %)
For (B = C; D + = f [B] * A, F [B] = D % -- g, D/= g --, -- B; D * = B );
}
I added the article explaining this code to the end. The author is Wang Cong. He looks very good at reading his blog. My main analysis is from his article. Thank you! The doubts of so many years have finally been solved.
From my perspective, I am not used to turning the for statement into a while for analysis. I think it looks quite good now. Of course, I still need to adjust some of the statements. I think we can fully understand the following issues: algorithms (two points), errors (two points), and others (four points ).
I. Algorithms
1. π Formula
π/2 = 1 + 1! /3 !! + 2! /5 !! + 3! /7 !! +... + K! /(2 * k + 1 )!! +...
This formula does not seem to have been learned in high numbers. Note: 5! = 1*2*3*4*5 !! = 1*3*5
2. Program Implementation of the formula (storing the remainder in an array)
Expand and adjust the above formula
π/2 = 1 + 1! /3 !! + 2! /5 !! + 3! /7 !! +... + K! /(2 * k + 1 )!!
= 1 + 1/3 + (1*2)/(3*5) + (1*2*3)/(3*5*7) +... + (1*2 *... * k)/(3*5 *... * (2 k + 1 ))
= 1 + 1/3 + (1/3) * (2/5) + (1/3) * (2/5) * (3/7) +... + (1/3) * (2/5 )*... * (K/(2 k + 1 ))
= (1/3) * (2/5 )*... * (K/(2 k + 1) +... + (1/3) * (2/5) * (3/7) + (1/3) * (2/5) + 1/3 + 1
= (K/(2 k + 1 ))*... * (2/5) * (1/3) +... + (3/7) * (2/5) * (1/3) + (2/5) * (1/3) + 1/3 + 1
= (K/(2 k + 1) + 1) * (k-1)/(2 (k-1) + 1) + 1 )*...) * 3/7 + 1) x 2/5 + 1) * 1/3 + 1)/1
= (1/(2 k + 1) * k + 1)/(2 (k-1) + 1) * (k-1) + 1 )/...) /7*3 + 1)/5*2 + 1)/3*1 + 1)/1
Now, we have to do division, multiplication, adding 1, division, multiplication, and adding 1 until division is divided by 1. Then we can complete it in a loop, and use a counter I, 1 divided by 2I + 1, multiply by I, add 1, and loop. However, if we want to output 800 decimal places, we can only use the method of division of large numbers: divide them into N rows for execution, which saves the remainder of each division in the formula, output after the commodity accumulation; the next step is to increase the accuracy of the stored remainder and then divide it, save the remainder, and output after the commodity accumulation, and so on. According to the mathematical operation, when k = 2800, it can be precise to 800 digits after the decimal point (actually accurate to 844 digits after the decimal point), then we set an array of 2800, initialize all to 1 (that is, set the remainder to 1, that is, 1 in the formula), and then continuously perform division, multiplication, and addition of remainder in the loop, and save the remainder in the process, it is the same as manual division.
One thing to note is to improve the accuracy. For example, F [2800], the initial value is 1, divided by 5601, multiplied by 2800, And the integer part of the quotient is added to 1 in the next operation, save the remainder, multiply it by 10 next, divide it by 5601, and multiply it by 2800. The integer part of the quotient obtained is the last digit of the result of the previous operation. If we want to output four digits at a time, we need to multiply the number by 10000. The obtained quotient's integer part will not exceed four digits, that is, the last four digits of the result of the previous operation. Note that the result of the above formula is π/2, then we initialize the array to 2. Because we need to output 4 bits at a time, it is increased by 10000 times.
It can be implemented by the following code:
# Define M 2800
# Include <stdio. h>
Long B, d, f [M + 1], I;
Void main (){
For (; B <= m;) f [B ++] = 1000*2;
For (; I <800/4; I ++, printf ("%. 4D", F [0] + D/10000), f [0] = D % 10000)
For (D = 0, B = m; D + = f [B] * 10000, F [B] = D % (2 * B + 1 ), d/= (2 * B + 1), B; B --)
D * = B;
}
If it is written like this, you can't understand it again:
# Define M 2800
# Include <stdio. h>
Void main ()
{
Long Yushu [M + 1], K = 0, I = 0, j = 0, sum = 0, Pi = 0;
For (k = 0; k <= m; k ++)
Yushu [k] = 1000*2;
For (I = 1; I <= 800/4; I ++)
{
For (j = m, sum = 0; j> 0; j --)
{
Sum + = Yushu [J] * 10000;
Yushu [J] = sum % (2 * j + 1 );
Sum/= (2 * j + 1 );
Sum * = J;
}
Sum + = Yushu [J] * 10000;
Yushu [J] = sum % (2 * j + 1 );
Printf ("%. 4D", PI + sum/10000 );
Pi = sum % 10000;
}
}
The first piece of code is slightly deformed, is to replace B with the B-1, you can get
# Define M 2800
# Include <stdio. h>
Long B, d, e, f [M + 1], I;
Void main (){
For (; B-m;) f [++ B] = 1000*2;
For (; I <800/4; I ++, printf ("%. 4D", e + D/10000), E = D % 10000)
For (D = 0, B = m; D + = f [B] * 10000, F [B] = D % (2 * b-1 ), d/= (2 * b-1), B> 1; B --)
D * = (b-1 );
}
Then we can continue to get the original code of the legendary paragraph. This process uses some blind-eye obfuscation methods, especially intentionally setting two errors without affecting the accuracy of the calculation results.
Ii. Errors
1. array Initialization
For (; B-c;) f [B ++] = A/5;
This line of code sets the array f [0] to f [2799] to 2000, F [2800] is not initialized, or 0, in actual calculation, this item does not play a role, because the result of this item is too small. However, in terms of computing logic, it is very confusing. In fact, F [0] has never been used, therefore, the initialization should be to initialize f [1] to f [2800] to 2000, so that the precision will only change after the decimal point, because the result is accurate to 848 digits, therefore, this accuracy is meaningless, but the algorithm is correct. That is
For (; B-c;) f [++ B] = A/5;
2. Division calculation process (the start point of Division calculation)
If you are careful, you can find that the initialization of the third for loop is a little different. The original code uses B = C, and we use B = m, that is, B = 2800, the difference is that every division is done according to the formula, and the original code misses one item each time. Because the results of each calculation are saved in the remainder array, therefore, this omission has a very small impact on the results, so it is negligible. However, such obfuscation is too powerful, and the original algorithms are almost totally different. Therefore, this error is critical to understanding the code.
Iii. Others
1. output result (printf)
In the code, D is the accumulation of operators in each operation, that is, the result to be output. According to the formula and our algorithm, D will not exceed 8 bits, in addition, the Left 4 bits and the last remainder e won't exceed 8 bits, so we use printf ("%. 4D ", e + D/10000) outputs the left four digits and stores the right four digits of D in E = D % 10000. Printf ("%. 4D ",...) Output a four-digit integer with 0 at the leftmost.
2. Two multiplier values (1000 and 10000)
The origins of the two numbers have been explained in the program implementation of the formula.
3. Functions of C
In the second for loop, C-= 14 is used to control the number of loops. The actual impact is that the initialization part of the third for loop is B = C. The above analysis shows the processing (error) precision is almost unaffected. In fact, C-= 14 has no special meaning, but because 2800/(800/4) = 14, we can use counters instead, just like in our code.
4. Role of G
G is the 2 k + 1 in the formula, which can be replaced by 2 * B + 1 in our program.
In our program, 2800 is the K in the formula. It is used to adjust the accuracy of the calculation result. If you can debug the program for a large number of times, I will use 6 for debugging. Counter I is also used. I <800/4 means to display 800 bits. In fact, we show 799 bits after the decimal point, and 3 before the decimal point is 800 bits in total.
The Analysis of F [2800] in the following article is actually not correct. The description of C-= 14 is as follows: "The role of C ", in the second for loop, the author did not mention the error of B = C, and the expansion of the formula is not yet in place, which is not consistent with the algorithm implementation in the program.
Interpreting Pi's weird code
Cong Wang
25th November, 2005
Institute of Post and Telecommunications, si'an, p.r. China
Network Engineering DEP.
Introduction
There is a weird PI program circulating on the Internet. Although there are only three rows, the PI value can be obtained to reach a total of 800 digits before the decimal point. This program is as follows:
/* A year of obfuscated C Contest :*/
# Include <stdio. h>
Long A = 10000, B, c = 2800, d, e, f [2801], G;
Main (){
For (; B-c;) f [B ++] = A/5;
For (; D = 0, G = C * 2; C-= 14, printf ("%. 4D", e + D/A), E = D %)
For (B = C; D + = f [B] * A, F [B] = D % -- g, D/= g --, -- B; D * = B );
}
/* (This program can calculate the PI value before the decimal point, a total of 800 digits)
(This program is recorded in sci. Math FAQ. The original author is unknown )*/
At first glance, this program is quite scary. Don't be flustered. Here we will show you how it works and some tips for writing weird C Programs. Pai_^
Expand and simplify
We know that in C, the for loop and the while loop can replace each other.
For (statement1; statement2; statement3 ){
Statements;
}
The preceding for statement can be replaced by the following while statement:
Statement1;
While (statement2 ){
Statements;
Statement3;
}
In addition, to write weird C Programs, the comma operator is undoubtedly a good assistant. Its role is:
Calculate the values of each expression from left to right, and return the value of the rightmost expression.
Embedding it in a for loop is one of the common tips for writing weird code. Therefore, the above program can be expanded:
# Include <stdio. h>/* 1 */
/* 2 */
Long A = 10000, B, c = 2800, d, e, f [2801], G;/* 3 */
Main () {/* 4 */
While (B-c! = 0) {/* 5 */
F [B] = A/5;/* 6 */
B ++;/* 7 */
}/* 8 */
D = 0;/* 9 */
G = C * 2;/* 10 */
While (G! = 0) {/* 11 */
B = C;/* 12 */
D + = f [B] * A;/* 13 */
F [B] = D % -- g;/* 14 */
D = D/g --;/* 15 */
-- B;/* 16 */
While (B! = 0) {/* 17 */
D = D * B + F [B] * A;/* 18 */
F [B] = D % -- g;/* 19 */
D = D/g --;/* 20 */
-- B;/* 21 */
}/* 22 */
C-= 14;/* 23 */
Printf ("%. 4D", e + D/A);/* 24 */
E = D % A;/* 25 */
D = 0;/* 26 */
G = C * 2;/* 27 */
}/* 28 */
}/* 29 */
Is it nice now?
Further simplification
You should be aware that the value of a is always 10000, so we can replace a with 10000. Then, observe g carefully. In the outer loop, when it is used for division or remainder, it is always equal to 2 * C-1, and B is always initialized to C. In an inner layer loop, B decreases by 1 each time, and g decreases by 2 each time. You will think of it now? Replace g with 2 * B-1! Try it on your behalf. That's right! In addition, we can also make a little simplification. The D = 0 of the 26th rows is redundant. We will merge it into the 13th rows, and the 13th rows can be rewritten: D = f [B] * ;. So the program can be changed:
# Include <stdio. h>
Long B, c = 2800, d, e, f [2801];
Main (){
While (B-c! = 0 ){
F [B] = 2000;
B ++;
}
While (C! = 0 ){
B = C;
D = f [B] x 10000;
F [B] = D % (B * 2-1 );
D = D/(B * 2-1 );
-- B;
While (B! = 0 ){
D = D * B + F [B] * 10000;
F [B] = D % (B * 2-1 );
D = D/(B * 2-1 );
-- B;
}
C-= 14;
Printf ("%. 4D", e + D/10000 );
E = D %10000;
}
}
Two variables are missing!
In-depth analysis
Now, we will go to the substantive analysis. Certain mathematical knowledge is very necessary. First, you must know that the following formula can be used to calculate pi:
PI/2 = 1 + 1! /3 !! + 2! /5 !! + 3! /7 !! +... + K! /(2 * k + 1 )!! +...
Pi has enough precision as long as there are enough items. As for why, we leave it to mathematicians.
Anyone who has written a high-precision division program knows how to use an integer array for division, instead of decimal places. In fact, it is very easy to think about how you perform division by hand. Divide the number by the divisor and place the obtained number into the result. Multiply the remainder by 10 and proceed to the next division until the remainder is zero or the required number of digits has been reached.
The original program uses so much mathematical knowledge that it is complicated and difficult to understand because it cleverly places algorithms in a loop. Let's start to analyze the program. First, let's start with the mathematical formula. We seek Pi, and the formula gives PI/2. Therefore, we multiply both sides of the formula by 2 at the same time:
Pi = 2*1 + 2*1! /3 !! + 2*2! /5 !! + 2*3! /7 !! +... + 2 * k! /(2 * k + 1 )!! +...
Next, we rewrite it into another form and expand it:
Pi = 2*1 + 2*1! /3 !! + 2*2! /5 !! + 2*3! /7 !! +... + 2 * n! /(2 * n + 1 )!!
= 2*(n-1)/(2 * N-1) * (n-2)/(2 * n-3) * (n-3)/(2 * N-5 )*... * 3/7x2/5*1/3
+ 2*(n-2)/(2 * n-3) * (n-3)/(2 * N-5) *... * 3/7*2/5*1/3
+ 2*(n-3)/(2 * N-5) *... * 3/7*2/5*1/3
+ 2*3/7*2/5*1/3
+ 2*2/5*1/3
+ 2*1/3
+ 2*1
Look at the formula, we can see that B Corresponds to the formula of N, 2 * B-1 corresponds to 2 * n-1. B starts from 2800, that is, n = 2800. (As to why n = 2800, the first 800 digits of Pi are not covered here .) Check the corresponding part of the program:
D = D * B + F [B] * 10000;
F [B] = D % (B * 2-1 );
D = D/(B * 2-1 );
D is used to store the integer part of the division result. It is accumulated, so the last d is the integer part we want. F [B] is used to store the remainder of the computation to B.
You may not understand it yet. First, why is the array containing 2801 elements? This is simple, because the author wants to use F [1] ~ F [2800], and the array subscript of C language starts from 0, and f [0] cannot be used. Second, why should we Initialize all array elements except f [2800] to 2000? What is the role of 10000? This is interesting. Because from printf ("%. 4D ", e + D/10000); we can see that D/10000 is the number before the first digit of D, and E = D % 4th ;, E is the last four digits of D. Furthermore, E and D are in a different loop. So the printed result is exactly the corresponding 4 bits of the PI we want! At the beginning, the reason why we initialize the array element to 2000 is that we can enlarge the PI by 1000 times to ensure that the integer part has four digits, and that 2 is the 2 multiplied by both sides of our formula! So it's 2000! Note that the remainder must be multiplied by 10000 instead of 10! F [2800] is 0 because the first multiplication is n-1, that is, 2799, rather than 2800! Every 4 bits are calculated, C minus 14, which ensures that a total of 4*2800/14 = 800 bits can be printed. However, this will not affect the accuracy of the results! I am not very familiar with mathematics and cannot give a definite answer. (If anyone knows, please email me to the xiyou.wangcong@gmail.com to tell me .)
By mistake, I saw an "accurate" adapted based on the above program on the Internet (this accurate refers to the absence of any element in the f [] array .) Print the 2800-bit program as follows:
Long B, c = 2800, d, e, f [2801], G;
Int main (INT argc, char * argv [])
{
For (B = 0; B <C; B ++)
F [B] = 2;
E = 0;
While (C> 0)
{
D = 0;
For (B = C; B> 0; B --)
{
D * = B;
D + = f [B] * 10;
F [B] = D % (B * 2-1 );
D/= (B * 2-1 );
}
C-= 1;
Printf ("% d", (e + D/10) % 10 );
E = D % 10;
}
Return 0;
}
Try compressing the above program into three lines.
Conclusion
End the full text with one sentence in knuth Turing's speech:
We have seen that computer programming is an art, because it applies accumulated knowledge to the world, because it requires skill and ingenuity, and especially because it produces objects of beauty. A programmer who subconsciously views himself as an artist will enjoy what he does and will do it better.
I would like to share with you! Pai_^
(EOF)