New Big Integer library is almost finished

[Thaddy]

I tried to do that, but I could not find the message numbers anywhere!  Where did you find them?

You can get the error numbers by compiling with -vewhnlq
This shows error numbers (q) for errors, warnings, hints and notes + line numbers.
You can also grep errore as mentioned. I should have a little tool for that somewhere which uses regular expressions.

[srvaldez]

...
Beside that the following arc tan series evaluation might be faster - at least it has a substantially lower Lehmer value

44*arctan( 1/57 ) + 7*arctan( 1/239 ) - 12*arctan( 1/682 ) + 24*arctan( 1/12943 )

Maybe you want to give it a try.

Cheers,
MathMan

I just tried it and the times are MUCH longer, the first few times I killed the program because I thought that the arctan function was in an infinite loop

@ad1mt, please try the following version, try each arctan formulas, why is the second formula about 22 times slower ?
it's this term that takes a very long time ArcTanInv(T,96,12943)

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. function ArcTanInv(const t: Multi_Int_XV; m,a: UInt64): Multi_Int_XV;
15. var
16.     r    : Multi_Int_XV;
17.     s    : Multi_Int_XV;
18.     b, k : UInt64;
19. begin
20.     s := t*(m*a); // s := 10**d * m * a
21.     b := a*a + 1;
22.     a := b + b;
23.     k := 0;
24.     r := 0;
25.     repeat
26.         s := s div b;
27.         if s=0 then
28.             Break;
29.         r := r + s;
30.         Inc(k,2);
31.         s := s * k;
32.         Inc(b,a);
33.     until FALSE;
34.     result:=r;
35. end;
36.  
37. const prec = 10000; // digits of precision
38.  
39. var
40.  
41.     big_int_size,
42.     start_time,
43.     end_time        :int32;
44.     delta           :double;
45.  
46.     t, x        :Multi_Int_XV;
47.     n:UInt32;
48.  
49. begin
50.     big_int_size:= 2 + (round(2 * prec * 3.321928094887362) div 64);
51.     Multi_Int_Initialisation(big_int_size);
52.  
53.     n := prec; // prec decimals
54.  
55.     start_time:= GetTickCount64;
56.     // t := 10**(n+4) (+4 digits for rounding errors)
57.     t:=10;
58.     t := t**(n+4);
59.     // Pi = 48 ArcTan(1/18) + 32 ArcTan(1/57) - 20 ArcTan(1/239)
60.    
61.     x := ArcTanInv(t,48,18) + ArcTanInv(t,32,57) - ArcTanInv(t,20,239);
62.     //x := ArcTanInv(T,176,57) + ArcTanInv(T,28,239) - ArcTanInv(T,48,682) + ArcTanInv(T,96,12943);
63.     // suppress the 4 added digits
64.     x := x div 100000;
65.     end_time:= GetTickCount64;
66.  
67.     writeln;
68.     writeln('Trunc(Pi * 10**' + IntToStr(n) + ')');
69.     writeln(x.ToStr);
70.     delta:= (end_time - start_time) / 1000;
71.     writeln('time elapsed is ', delta:1:6, ' seconds');
72.  
73. end.
74.  

[srvaldez]

I think that the problem is that there's no arithmetic support for Int64, Int64 gets promoted to a Multi_Int_XV and then the operation is performed in Multi_Int_XV precision

[MathMan]

I think that the problem is that there's no arithmetic support for Int64, Int64 gets promoted to a Multi_Int_XV and then the operation is performed in Multi_Int_XV precision

Ah, strange. I was aware of the issue that a limiting factor is the size of the half-limb (UInt32 on a 64 bit system). Therefore I picked a series which could be calculated with 32 bit divisions - the maximum divisor should be 12943^2 = 167521249 < 2^32-1. Maybe I should have looked more closely into your series expansion for arctan ...

Cheers,
MathMan

[ad1mt]

I think that the problem is that there's no arithmetic support for Int64, Int64 gets promoted to a Multi_Int_XV and then the operation is performed in Multi_Int_XV precision

That's more or less what is going on.
The divisor (b) always gets promoted to Multi_Int_XV, however...
When the value is < (2^32), it fits into a single 32-bit word in the Multi_Int_XV array. Then it gets processed by a "fast track" loop in the division function.
But when the value is >= (2^32), it requires two or more 32-bit words in the Multi_Int_XV array. Then it gets processed by the much more complex & slower "long division" loop.
I've never timed them before, and it has surprised me just how much slower the "long division" loop is compared to "fast track" loop. I'm going to do some more investigation to see if there are some efficiency gains that I might have missed, or bugs I might have introduced when I "re-engineered" the Knuth algorithm. I don't think there can be any logic bugs, but there might be some "efficiency bugs" (if that makes sense).

[MathMan]

code borrowed from https://www.ellipsa.eu/public/fnx/fnx.html

I think FNX is very good. If it worked in 32bit environments, I would not have created my library.

Hi ad1mt,

May I ask what exactly was not working with FNX in a 32bit environment? Was there a compiler issue, were calculations not working correct?

Cheers,
MathMan

[ad1mt]

May I ask what exactly was not working with FNX in a 32bit environment? Was there a compiler issue, were calculations not working correct?

It says in the source and/or documentation that it's "untested" in 32bit environment.
However, I tried compiling it on the 32bit compiler, in case I got lucky. It issued several warnings about not being tested for 32bits, then the compilation failed with "fnx.intutils.pas(2134,5) Error: Internal error 200706094".
I did have a quick look at the source code, with a view to making it work, but the code was so complex, I was way out of my league.

[TRon]

Error: Internal error 200706094".

Internal (compiler) errors are always reason to report in the bug-tracker.

compiler/constexp.pas

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

1. ...
2. else if c.uvalue>qword(high(int64)) then
3.     internalerror(200706094)
4.   else
5.  ...
6.  

[Thaddy]

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

1. operator := (const c:Tconstexprint):int64;
2.  
3. begin
4.   if c.overflow then
5.     internalerror(200706093)
6.   else if c.signed then
7.     result:=c.svalue
8.   else if c.uvalue>qword(high(int64)) then
9.     internalerror(200706094)
10.   else
11.     result:=int64(c.uvalue);
12. end;

the error disappears when you do this, but I am afraid I do not know why, yet:

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

1. operator := (const c:Tconstexprint):int64;
2.  
3. begin
4.   if c.overflow then
5.     internalerror(200706093)
6.   else if c.signed then
7.     result:=c.svalue
8.   else if c.uvalue>=qword(high(int64)) then
9.     internalerror(200706094)
10.   else
11.     result:=int64(c.uvalue);
12. end;

I do not even know if that is legal. Just tried it.
But it should cover just the - in hex notation - left most difference between $F and $E. I think that is the intention.
High(int64) starts with $E.
It is also the only place in the compiler code where that specific error is mentioned.
Maybe just:

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

1.   else c.uvalue>high(int64) then // left most bit is set to one, which makes it $F

since c.uvalue is unsigned.

[MathMan]

The offending line in "fnx.intutils.pas" is

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

1.     Dec(A,UInt64(Low(Int64)));
2.  

'A' is defined as UInt64 in the function header

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

1.   function U64Cbrt(A: UInt64): UInt64;
2.  

Further on Low(Int64) = $80000000000000000 = -2^63. The compiler seems to stumble on the cast to UInt64.

When I compiled this on a 64 bit system (x86-64, FPC 3.2.2) no compiler error was thrown.

ad1mt - which version of FPC did you use when you tried to compile FNX?

Cheers,
MathMan

[srvaldez]

I tested Multi_int, Mparth and fnx with the Pi Chudnovsky program, and while fnx was the fastest, when I tried to go above 50000 digits it would crash with an overflow error, with Mparith I tested up to 1,000,000 digits in just over 17 seconds
 

[ad1mt]

ad1mt - which version of FPC did you use when you tried to compile FNX?

FPC v3.2.2 in Lazarus v2.2.6

[ad1mt]

I tested Multi_int, Mparth and fnx with the Pi Chudnovsky program, and while fnx was the fastest, when I tried to go above 50000 digits it would crash with an overflow error, with Mparith I tested up to 1,000,000 digits in just over 17 seconds

Maybe FNX is not as good as I thought   

[srvaldez]

actually, you can change the maximum size of the FNX BigInteger, in fnx.bigint.pas change the MAX_SIZE and MAX_BIT_SIZE to the following

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

1.     MAX_SIZE = Int32(1) shl 18;
2.     MAX_BIT_SIZE = Int32(MAX_SIZE shl 6);
3.  

now you can compile and run the Pi Chudnovsky program to 1,000,000 digits, it's about twice as fast as Mparith
but it would seem that some kind of check should be made in the FNX BigInteger to prevent a crash
<edit>
Mparith is different in that the size of the BigInteger grows as needed and is not of fixed size, that can be an advantage as you don't have to worry about overflow

[ad1mt]

I've made version 4.70 available.
It has some obscure but major bug fixes in the division functines.

I discovered the bugs while investigating the speed difference between the "fast track" and "long track" paths.
I still don't understand reason for the big speed difference. I will keep looking.