New Big Integer library is almost finished

[ad1mt]

I've just made version 4.74 available.
This has bug fixes, and faster subtraction algorithm.

[srvaldez]

just for fun https://www.flyingcoloursmaths.co.uk/ask-uncle-colin-is-the-fibonacci-series-witchcraft/

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 fibmax=256;
15.  
16. function numberOfDigits(const n : longint) : longint;
17. var d : double;
18. begin
19.     if n = 1 then numberOfDigits:= 1
20.     else
21.     begin
22.         d := 0.2089876402499788*n-0.3494850021680092;
23.         numberOfDigits:=ceil(d);
24.     end;
25. end;
26.  
27. var
28.     fib, one, nines : Multi_Int_XV;
29.     big_int_size,
30.     start_time,
31.     end_time, lngth :longint;
32.     delta           :double;
33.     i, j            :longint;
34.     s: ansistring;
35.    
36.  
37.  
38. begin
39.     lngth := numberOfDigits(fibmax);
40.     big_int_size:= 4 + (round(2 * (fibmax*lngth+lngth) * 3.321928094887362) div 64);
41.     Multi_Int_Initialisation(big_int_size);
42.  
43.     start_time:= GetTickCount64;
44.     one:=10;
45.     one:=one**(fibmax*lngth+lngth);
46.     s:=DupeString('9', lngth-1)+'8'+DupeString('9', lngth);
47.     nines:=s;
48.     fib:=one div nines;
49.     s:=fib.ToStr;
50.     s:=DupeString('0', lngth-1)+s;
51.     j:=1;
52.     for i:=1 to (fibmax div 2) do
53.     begin
54.         writeln(MidStr(s, j, lngth), '   ', MidStr(s, j+lngth, lngth));
55.         j+=(2*lngth);
56.     end;
57.     end_time:= GetTickCount64;
58.  
59.     writeln;
60.  
61.     delta:= (end_time - start_time) / 1000;
62.     writeln('time elapsed is ', delta:1:6, ' seconds');
63.  
64. end.
65.  

[ad1mt]

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

1. program app;
2.  
3. uses Multi_Int;
4.  
5. procedure test1;
6. var
7.   i: Int32;
8.   m: Multi_Int_X3;
9. begin
10.   writeln('# Int32 (signed)');
11.   // assign highest possible signed int
12.   i := high(integer);    
13.   writeln('i    = ', i:11);
14.   // and overflow it
15.   i += 1;        
16.   writeln('i +1 = ', i:11);
17.  
18.   writeln;
19.                
20.   writeln('# Multi_Int_X3 (?)');
21.   // assign highest possible signed int
22.   // m := high(integer);  
23. // Please do this instead
24.   m:= Multi_Int_X3_MAXINT;
25.   writeln('m    = ', m.ToStr:11);
26.   // and overflow it
27.   m += 1;    
28.   writeln('m +1 = ', m.ToStr:11);
29.  
30.   writeln;
31.  
32.   writeln('value correct? ', i = int64(m), ', i = ', i, ', m = ', int64(m));
33.   writeln('bits the same? ', binstr(i, 64) = binstr(int64(m), 64));
34.   writeln('i bits = ', binstr(i, 64));
35.   writeln('m bits = ', binstr(int64(m), 64));
36.   writeln('m bits = ', m.ToBin:64);
37. end;
38.  
39. procedure test2;
40. var
41.   i: Int32;
42. begin
43.   writeln('# Int32');
44.   i := -123;  
45.   writeln('value = ', i);
46.   writeln('bits                   = ', binstr(i, 33), ' value = ', Int64(i):11);
47.   i := -123;
48.   i := i xor (1 shl 31);
49.   writeln('bits xored without ABS = ', binstr(i, 33), ' value = ', Int64(i):11);
50.   i := -123;
51.   i := abs(i);
52.   i := i xor (1 shl 31);
53.   writeln('bits xored    with ABS = ', binstr(i, 33), ' value = ', Int64(i):11);
54. end;
55.  
56. procedure test3;
57. var
58.   i: Multi_Int_X3;
59. begin
60.   writeln('# Multi_Int_X3');
61.   i := -123;
62.   writeln('value = ', i.ToStr);
63.   writeln('bits                   = ', i.ToBin:33, ' value = ', Int64(i):11); // using Int64 because of overflow
64.   i := -123;
65.   i := i xor (1 shl 31);
66.   writeln('bits xored without ABS = ', i.ToBin:33, ' value = ', Int64(i):11);
67.   i := -123;
68.   i := abs(i);
69.   i := i xor (1 shl 31);
70.   writeln('bits xored    with ABS = ', i.ToBin:33, ' value = ', Int64(i):11);
71. end;
72.  
73. procedure test4;
74. var
75.   i: Multi_Int_X3;
76. begin
77.   writeln('# Multi_Int_X3');
78.   writeln('# cast to Int64 before converting to binary');
79.   i := -123;
80.   writeln('value = ', i.ToStr);
81.   writeln('bits                   = ', binstr(int64(i), 33):33, ' value = ', Int64(i):11); // using Int64 because of overflow
82.   i := -123;
83.   i := i xor (1 shl 31);
84.   writeln('bits xored without ABS = ', binstr(int64(i), 33):33, ' value = ', Int64(i):11);
85.   i := -123;
86.   i := abs(i);
87.   i := i xor (1 shl 31);
88.   writeln('bits xored    with ABS = ', binstr(int64(i), 33):33, ' value = ', Int64(i):11);
89. end;
90.  
91. begin
92.   writeln('****** test 1 ******');
93.   test1;
94.   writeln;
95.   writeln;
96.  
97.   writeln('****** test 2 ******');
98.   test2;
99.   writeln;
100.   writeln;
101.  
102.   writeln('****** test 3 ******');
103.   test3;
104.   writeln;
105.   writeln;
106.  
107.   writeln('****** test 4 ******');
108.   test4;
109.   writeln;
110.   writeln;
111.  
112.   readln;
113. end.

@Fibonacci thanks for your bug report. It has highlighted some bugs, and I've just made version 4.75 available to fix them.

There are a couple of issues with your program.

My library will not work as expected with the built-in high function; I am looking for a solution but it might not be possible. Please use the Multi_Int_X3_MAXINT constant for this... see the alternative code above. Once that has been changed, the M variable then overflows as expected. The built-in high function returns the highest index of an array type, but I'm not sure what it is doing with my big integer types.

One of the changes I've made in the latest version is that Multi_Int_RAISE_EXCEPTIONS_ENABLED is now set = FALSE by default. This means that if you want exceptions to be raised by my library, you must set this variable = TRUE in your code.
EDIT: I've since decided this was a mistake see my later post.

The other issue is that you are trying to compare bit paterns of integers of differing sizes and different internal representations. For example, the bit pattern of  -1 stored in a 32bit or 64bit integer will be completely different to the bit pattern of -1 stored on a big integer. This is because the internal representation is not the same. In a 32bit or 64bit integer, the representation is two's complement. Whereas my representation the number is stored as an absolute value with a separate sign flag.

I am mistaken, please let me know.

[ad1mt]

just for fun https://www.flyingcoloursmaths.co.uk/ask-uncle-colin-is-the-fibonacci-series-witchcraft/

@srvaldez... you do not say whether your program worked!

[srvaldez]

yes it works in 64-bit, didn't check 32-bit until just now, there's an access violation

[srvaldez]

it works in 32-bit if you double the precision

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

1. big_int_size:= 4 + (round(4 * (fibmax*lngth+lngth) * 3.321928094887362) div 64);

[ad1mt]

I've just version 4.76 available. This fixes the bugs revealed by the latest bug report.

One of the changes I've made in the latest version is that Multi_Int_RAISE_EXCEPTIONS_ENABLED is now set = FALSE by default. This means that if you want exceptions to be raised by my library, you must set this variable = TRUE in your code.

I decided this was a mistake, and Multi_Int_RAISE_EXCEPTIONS_ENABLED is now set = TRUE by default.

[ad1mt]

it works in 32-bit if you double the precision

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

1. big_int_size:= 4 + (round(4 * (fibmax*lngth+lngth) * 3.321928094887362) div 64);

There was nothing wrong with your program, and it should have worked in the 32bit world. It does now.

[ad1mt]

I've just made version 4.77 available.
This has bug fixes and some code tidying-up.

[ad1mt]

ArcTanInv(T,96,12943)

@srvaldez or @Mathman
I'm curious about the maths behind the ArcTanInv function. Its obviously using integer maths operations to do decimal maths calculations, but I don't understand how.
Please could you explain, in very simple terms, how its working, or point me to a source that would explain it?
Thanks.
It is rather amusing that I'm trying to write a big integer maths library, and yet my maths is not that good!   
The answer is that although I'm not that good it, I have always been interested in maths.

[MathMan]

@ad1mt

It is quite simple (I hope). If you calculate ArcTanInv in decimal notation you'd get 0.xxxxxx..xxx with k decimal digits behind the decimal point.

Now scale up the calculation by multiplication with 10^k and you get an integer with k decimal digits in front of the decimal point. Effectively both represent the same thing.

That is why you always find an initial term where 10^k is calculated as a new representative for '1' and all the subsequent divisions / multiplications work on this.

Hope this explains it somehow?

Cheers,
MathMan

[ad1mt]

I've just made version 4.78 available.
The faster subtract algorithm is now used in the add function. It also has bug fixes.
One day, I hope to have fixed all the bugs.

[ad1mt]

Now scale up the calculation by multiplication with 10^k and you get an integer with k decimal digits in front of the decimal point. Effectively both represent the same thing.

@Mathman
Thanks for your reply...
So do you decide in advance how many decimal places you want to calculate to, then scale all values up by that factor?
So for example, to calculate to 6 decimal places, you multiply all values by 1 million, then at the end, divide all values by 1 million?

[MathMan]

Now scale up the calculation by multiplication with 10^k and you get an integer with k decimal digits in front of the decimal point. Effectively both represent the same thing.

@Mathman
Thanks for your reply...
So do you decide in advance how many decimal places you want to calculate to, then scale all values up by that factor?
So for example, to calculate to 6 decimal places, you multiply all values by 1 million, then at the end, divide all values by 1 million?

I'll answer in light of the capabilities of your library, if that is ok. Your statement is nearly correct with the exception of the final division. Instead you print out the integer result and state something like "This is [whatever] * 10^6". Take a look at the examples of srvaldez for pi calculation - you'll find exactly that.

[ad1mt]

then at the end, divide all values by 1 million?

What I meant to say and should have said is... display the number as though it had been divided by 1 million.