New Big Integer library is almost finished

[ad1mt]

I've just made version 4.79 available.
It has obscure, but important, bug fixes related to the problems discussed here... https://forum.lazarus.freepascal.org/index.php/topic,69092.0.html

I'm confident that the latest fixes will enable substantial improvements to the efficiency of the code, in the near future.

[ad1mt]

I've just made version 4.82 available.
This has significant improvements throughout the code.

The improvements are due to the use of management operators for extended/advanced records. This has significantly increased the robustness, and has eliminated several problems/bugs which had proved very difficult to fix. It has also enabled faster execution times.

As always, bug reports & suggestions are always welcome.

[srvaldez]

thanks ad1mt 
I noticed that the thread title changed a bit 

[srvaldez]

ad1mt
my previous tests run much faster except for Pi - Chudnovsky https://forum.lazarus.freepascal.org/index.php/topic,68290.msg534208.html#msg534208
no matter how high I set the precision it crashes, it could very well be a bug in my code, don't know at the moment

[ad1mt]

I've just made version 4.84 available.
This has much faster algorithms for converting from decimal and heximal strings.
Plus many bug fixes.

[ad1mt]

my previous tests run much faster except for Pi - Chudnovsky https://forum.lazarus.freepascal.org/index.php/topic,68290.msg534208.html#msg534208
no matter how high I set the precision it crashes, it could very well be a bug in my code, don't know at the moment

@srvaldez
Thanks for your report.
Please try the version 4.84 that I've just uploaded.

Also please adjust the size limit you have set. Previously, the size you specified was a "base size", and the library set a limit of twice the base size. Now, the size you specify is the limit, which means you have to double the sizes you specify in your code.

This is a side-effect of the faster algorithms. The latest versions of the library have much better control/checking of the size of Multi_Int_XV types. This enabled me to have a base size of just 2, with the big int variables enlarging or reducing in size depending on the size of the value they contained. This has resulted in significant speed increases, because big int variables that only contain a small value now no longer have large base size containers. However, the library no longer has a concept of a base size any more, and the size you specify is the limit.

I have included both of your Pi calculation programs in my test suite, and they work ok for me.

[srvaldez]

hello ad1mt 
please give an example of setting the limit, suppose that I want to set the limit to 1024 decimal digits

[ad1mt]

hello ad1mt 
please give an example of setting the limit, suppose that I want to set the limit to 1024 decimal digits

The limit is specified in half-words. So in a 64bit world, a half-word is 32bits, and in a 32bit world a half-word is 16bits.
The maximum value is calculated as = ( (2^W) ^ N) - 1; where W=half-word size in bits, and N=number of half-words.
So, in a 64bit world, a size/limit of 10 has a max value of

Code: Text  [Select][+][-]

1. ((2^32)^10)-1
2. = 2135987035920910082395021706169552114602704522356652769947041607822219725780640550022962086936575

You need to think in powers of 2.
To answer your specific question, a size/limit of 107 would give a maximum value slightly greater than 1030 decimal digits (but less than the full range of 1031 digits). A size/limit of 106 would give a maximum value slightly greater than 1021 decimal digits (but less than 1022 digits).

I'll have a think, and see if I can come up with a way of making it simpler.

[srvaldez]

thanks ad1mt
so for 64-bit limit = ceil(digits/(32*log10(2))) and for 32-bit limit = ceil(digits/(16*log10(2)))

[srvaldez]

there's something wrong, the following should compute Pi using the Chudnovsky algorithm, only the first 14 digits are correct
the exact same program works with Mparith using mp_int

Code: Pascal  [Select][+][-]

1. {$mode objfpc}
2. {$MODESWITCH NESTEDCOMMENTS+}
3. {$warn 6058 off}
4. {$warn 5025 off}
5.  
6. program speed_test;
7. uses    sysutils
8. ,       strutils
9. ,       strings
10. ,       math
11. ,       Multi_Int
12. ;
13.  
14. const prec = 50; // digits of precision
15.  
16. var
17.     c3_over_24, temp : Multi_Int_XV;
18.     big_int_size,
19.     start_time,
20.     end_time        :int32;
21.     delta           :double;
22.  
23.     t, pi        :Multi_Int_XV;
24.     n:UInt32;
25.    
26. procedure bs(a : uint32; b : uint32; var pab : Multi_Int_XV; var qab : Multi_Int_XV; var tab_ : Multi_Int_XV);
27.     (*
28.     computes the terms for binary splitting the chudnovsky infinite series
29.  
30.     a(a) = +/- (13591409 + 545140134*a)
31.     p(a) = (6*a-5)*(2*a-1)*(6*a-1)
32.     b(a) = 1
33.     q(a) = a*a*a*c3_over_24
34.  
35.     returns p(a,b), q(a,b) and t(a,b) *)
36.    
37.     var
38.         m, t : uint32;
39.         pam, qam, tam, pmb, qmb, tmb : Multi_Int_XV;
40.        
41.     begin
42.     if (b - a) = 1 then
43.     begin
44.         // directly compute p(a,a+1), q(a,a+1) and t(a,a+1)
45.         if a = 0 then
46.         begin
47.             pab := 1;
48.             qab := 1;
49.         end
50.         else
51.         begin
52.             t:=6*a;
53.             pab := t-1;
54.             pab := pab*(a+a-1);
55.             pab := pab*(t-5);
56.             temp := a;
57.             qab := temp*temp;
58.             qab := qab*temp;
59.             qab := qab*c3_over_24;
60.         end;
61.         tab_ := 545140134;
62.         tab_ := tab_*a;
63.         tab_ := tab_+13591409;
64.         tab_ := tab_*pab; // a(a) * p(a)
65.         if odd(a) then tab_ := tab_*(-1); // = -tab
66.     end
67.     else
68.     begin
69.         // recursively compute p(a,b), q(a,b) and t(a,b)
70.         // m is the midpoint of a and b
71.         m := (a + b) div 2;
72.         // recursively calculate p(a,m), q(a,m) and t(a,m)
73.         bs(a, m, pam, qam, tam);
74.         // recursively calculate p(m,b), q(m,b) and t(m,b)
75.         bs(m, b, pmb, qmb, tmb);
76.         // now combine
77.         pab := pam*pmb;
78.         qab := qam*qmb;
79.         temp := qmb*tam;
80.         tab_ := pam*tmb;
81.         tab_ := tab_+temp;
82.     end;
83. end;
84.  
85. procedure pi_chudnovsky(var pi : Multi_Int_XV; digits : uint32);
86.         (*
87.         compute int(pi * 10**digits)
88.  
89.         this is done using chudnovsky's series with binary splitting
90.         *)
91.     var
92.         one_squared, p, q, sqrtc, t : Multi_Int_XV;
93.         n : uint32;
94.         digits_per_term : double;
95.  
96.     //     how many terms to compute
97.     begin
98.         //digits_per_term := c3_over_24;
99.         digits_per_term := 14.18164746272548; //Ln(digits_per_term/6/2/6)*0.4342944819032518;
100.         n := round(digits/digits_per_term + 1);
101.     //     calclate p(0,n) and q(0,n)
102.         bs(0, n, p, q, t);
103.         one_squared := 10;
104.         one_squared := one_squared**(2*digits);
105.         sqrtc := one_squared*10005;
106.         sqrtc := SqRoot(sqrtc);
107.         pi := sqrtc*426880;
108.         pi := pi*q;
109.         pi := pi div t;
110. End;
111.  
112. begin
113.     big_int_size:= 3 *ceil( prec / (32*0.3010299956639812));
114.     Multi_Int_Set_XV_Limit(big_int_size); // don't know whether it's needed or not
115.     Multi_Int_Initialisation(big_int_size);
116.    
117.    
118.     start_time:= GetTickCount64;
119.    
120.     c3_over_24 := 640320;
121.     temp := c3_over_24*c3_over_24;
122.     c3_over_24 := temp*c3_over_24;
123.     c3_over_24 := c3_over_24 div 24;
124.        
125.     pi_chudnovsky(pi, prec);
126.    
127.     end_time:= GetTickCount64;
128.  
129.     writeln;
130.     writeln('Trunc(Pi * 10**' + IntToStr(prec) + ')');
131.     writeln(pi.ToStr);
132.     delta:= (end_time - start_time) / 1000;
133.     writeln('time elapsed is ', delta:1:6, ' seconds');
134. end.
135.  

[srvaldez]

it works with Multi_Int_4.71 with a slight modification due to ambiguous overloaded function in an Uint32 expression

Code: Pascal  [Select][+][-]

1. {$mode objfpc}
2. {$MODESWITCH NESTEDCOMMENTS+}
3. {$warn 6058 off}
4. {$warn 5025 off}
5.  
6. program speed_test;
7. uses    sysutils
8. ,       strutils
9. ,       strings
10. ,       math
11. ,       Multi_Int
12. ;
13.  
14. const prec = 50; //50000; // digits of precision
15.  
16. var
17.     c3_over_24, temp : Multi_Int_XV;
18.     big_int_size,
19.     start_time,
20.     end_time        :int32;
21.     delta           :double;
22.  
23.     t, pi        :Multi_Int_XV;
24.     n:UInt32;
25.    
26. procedure bs(a : uint32; b : uint32; var pab : Multi_Int_XV; var qab : Multi_Int_XV; var tab_ : Multi_Int_XV);
27.     (*
28.     computes the terms for binary splitting the chudnovsky infinite series
29.  
30.     a(a) = +/- (13591409 + 545140134*a)
31.     p(a) = (6*a-5)*(2*a-1)*(6*a-1)
32.     b(a) = 1
33.     q(a) = a*a*a*c3_over_24
34.  
35.     returns p(a,b), q(a,b) and t(a,b) *)
36.    
37.     var
38.         m, t, tmp : uint32;
39.         pam, qam, tam, pmb, qmb, tmb : Multi_Int_XV;
40.        
41.     begin
42.     if (b - a) = 1 then
43.     begin
44.         // directly compute p(a,a+1), q(a,a+1) and t(a,a+1)
45.         if a = 0 then
46.         begin
47.             pab := 1;
48.             qab := 1;
49.         end
50.         else
51.         begin
52.             t:=6*a;
53.             pab := t-1;
54.             tmp:=a+a-1;
55.             pab := pab*tmp;
56.             tmp:=t-5;
57.             pab := pab*tmp;
58.             temp := a;
59.             qab := temp*temp;
60.             qab := qab*temp;
61.             qab := qab*c3_over_24;
62.         end;
63.         tab_ := 545140134;
64.         tab_ := tab_*a;
65.         tab_ := tab_+13591409;
66.         tab_ := tab_*pab; // a(a) * p(a)
67.         if odd(a) then tab_ := tab_*(-1); // = -tab
68.     end
69.     else
70.     begin
71.         // recursively compute p(a,b), q(a,b) and t(a,b)
72.         // m is the midpoint of a and b
73.         m := (a + b) div 2;
74.         // recursively calculate p(a,m), q(a,m) and t(a,m)
75.         bs(a, m, pam, qam, tam);
76.         // recursively calculate p(m,b), q(m,b) and t(m,b)
77.         bs(m, b, pmb, qmb, tmb);
78.         // now combine
79.         pab := pam*pmb;
80.         qab := qam*qmb;
81.         temp := qmb*tam;
82.         tab_ := pam*tmb;
83.         tab_ := tab_+temp;
84.     end;
85. end;
86.  
87. procedure pi_chudnovsky(var pi : Multi_Int_XV; digits : uint32);
88.         (*
89.         compute int(pi * 10**digits)
90.  
91.         this is done using chudnovsky's series with binary splitting
92.         *)
93.     var
94.         one_squared, p, q, sqrtc, t : Multi_Int_XV;
95.         n : uint32;
96.         digits_per_term : double;
97.  
98.     //     how many terms to compute
99.     begin
100.         digits_per_term := c3_over_24;
101.         digits_per_term := Ln(digits_per_term/6/2/6)*0.4342944819032518;
102.         n := round(digits/digits_per_term + 1);
103.     //     calclate p(0,n) and q(0,n)
104.         bs(0, n, p, q, t);
105.         one_squared := 10;
106.         one_squared := one_squared**(2*digits);
107.         sqrtc := one_squared*10005;
108.         sqrtc := SqRoot(sqrtc);
109.         pi := sqrtc*426880;
110.         pi := pi*q;
111.         pi := pi div t;
112. End;
113.  
114. begin
115.     big_int_size:= 4 + (round(3 * prec * 3.321928094887362) div 64);
116.     Multi_Int_Initialisation(big_int_size);
117.  
118.     start_time:= GetTickCount64;
119.    
120.     c3_over_24 := 640320;
121.     temp := c3_over_24*c3_over_24;
122.     c3_over_24 := temp*c3_over_24;
123.     c3_over_24 := c3_over_24 div 24;
124.    
125.     pi_chudnovsky(pi, prec);
126.    
127.     end_time:= GetTickCount64;
128.  
129.     writeln;
130.     writeln('Trunc(Pi * 10**' + IntToStr(prec) + ')');
131.     writeln(pi.ToStr);
132.     delta:= (end_time - start_time) / 1000;
133.     writeln('time elapsed is ', delta:1:6, ' seconds');
134.  
135. end.
136.  

[srvaldez]

actually it doesn't work with version 4.71 either, it I set the precision to 50000 then Pi is accurate to bout 2100 digits, if I set the precision to 2150 then only about 154 digits are accurate

[ad1mt]

actually it doesn't work with version 4.71 either, it I set the precision to 50000 then Pi is accurate to bout 2100 digits, if I set the precision to 2150 then only about 154 digits are accurate

@srvaldez
Many thanks for your bug reports. It never ceases to amaze me that I can build & use a big test suite, but then someone else writes a relatively small program that breaks my library. You really are helping a lot... please keep going!

I'm working through past versions of the library, to see which, if any, work. So far, I have to go back to version 4.70 before your latest Pi program works (I'm calling this Pi version 3). Version 4.71 is partially broken, and version 4.72 gives the same incorrect output as all later versions.

Now I have both a working & non-working version of the library to compare, which will greatly help the debugging process, because I can compare exactly what code changes were made between those versions. I can also step through the code in each version at the point where the results diverge, to see exactly what the cause is.

Some programmers dislike debugging, but I find it very satisfying!

[srvaldez]

😁👍

[ad1mt]

the following should compute Pi using the Chudnovsky algorithm, only the first 14 digits are correct

@srvaldez
That was easier than I expected... its a good thing that I have strict version control!

I've just made version 4.85 available. Your latest Pi program now gives the same results as your previous Pi programs.

Looking forward to the next bug!