New Big Integer library is finished

[ad1mt]

I have been working on a new Big Integer library, which I now consider to be finished.
It is designed to be easy to use, with minimal changes required to existing code.
It is written entirely in Pascal, with no assembly or C code, so should be portable to any platform that Free Pascal supports. It has been beta tested on 32bit Solaris Sun, 32bit Raspberry Pi running Linux, and 64bit Intel Windows. It has also been used in another project of mine over several months.
It is reasonably fast, but there are other big integer libraries that are significantly faster.
Available on GITHub  https://github.com/ad1mt/Multi-Word-Int
Or from my web page  https://mark-taylor.me.uk/index.php?page=Software
I have made the code public domain; no copyright, no license, no warranty. Do whatever you like with this code.
When I first started work on this library, I expected it to be fairly easy/quick, but it took me almost 2 years to finish it; so I hope some people will find it useful and save them a lot of work.
Although I consider the code to be finished, bug reports and suggestions are still welcome.

[srvaldez]

hello ad1mt 
there's a problem, possibly with conversion from longint to Multi_Int_XV, see the comments in the code

Code: [Select]

{$MODESWITCH NESTEDCOMMENTS+}

program speed_test;
uses sysutils
, strutils
, strings
, math
, Multi_Int
;

var

big_int_size,
start_time,
end_time :int32;
delta :double;

MV_i,
MV_j,
MV_k :Multi_Int_XV;
n, mm:longint;

begin
big_int_size:= 1876;
Multi_Int_Initialisation(big_int_size);
mm:=99999999;

start_time:= GetTickCount64;

MV_i:=1;
for n:=2 to 10000 do
begin
MV_i:=MV_i*n;
end;
end_time:= GetTickCount64;
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');
MV_j := '99999999'; // mm or 99999999 crashes the div loop
start_time:= GetTickCount64;
for n:=1 to 10 do
begin
MV_i:=MV_i div MV_j; //mm ; // with mm it crashes
end;
end_time:= GetTickCount64;
writeln(MV_i.ToStr);
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');

end.


[abouchez]

Nice code. And thanks for sharing!
I am only concerned about the fact that I did not find a lot of unit tests, which are imho mandatory for such a general purpose library.

I like the fact that you use what I call "half-word" types, i.e. 16-bit on 32-bit CPU, and 32-bit on 64-bit CPU.
This is what I did for my RSA calculation unit, and even in "pure pascal", the FPC compiler gives good performance.
https://github.com/synopse/mORMot2/blob/master/src/crypt/mormot.crypt.rsa.pas#L57

[ad1mt]

I am only concerned about the fact that I did not find a lot of unit tests, which are imho mandatory for such a general purpose library.

I have a test program, but I did not make it available, because I did not think anyone would be interested.
unit_test_v12_X.pas is now available on Github in the Demo folder.
N.B. because of the bug that has been reported nearby, the unit_test program will need to be updated soon.

[ad1mt]

hello ad1mt 
there's a problem, possibly with conversion from longint to Multi_Int_XV, see the comments in the code

Thank you for the bug report.
O.M.G! 
I will post a fix ASAP.

[ad1mt]

I like the fact that you use what I call "half-word" types, i.e. 16-bit on 32-bit CPU, and 32-bit on 64-bit CPU.

I did not really like having to do this, but I could not think of a better/faster method.
Another method would be to use the full size of 32/64 bit words, but I would then have to detect overflow, and then when overflow was detected, redo the calculation again in a safe way.
I assumed that the second method would be slower, although I did not test this, and maybe I should have. However, even if the second method was faster, it would have made the code significantly more complex, and it was complex enough already.

[ad1mt]

hello ad1mt 
there's a problem, possibly with conversion from longint to Multi_Int_XV, see the comments in the code

Thank you for the bug report.
O.M.G! 
I will post a fix ASAP.

Version 4.39 is now on Github and my webpage.
This fixes the bug reported by srvaldez.
Many thanks for your bug report.
Note that the unit test program I made available earlier, does not yet have a regression test for this particular bug; an updated unit test program will follow.

[srvaldez]

thanks ad1mt 

[ad1mt]

thanks ad1mt 

Please let me know if you find any other problems or if you have any suggestions.

[srvaldez]

hello ad1mt, silly tests often reveal bugs
the following silly test sets MV_i to factorial 10000, then in the following loop it divides MV_i by the sqrt(sqrt(sqrt(sqrt(sqrt(factorial 10000))))) 32 times, the result should be 1 but it's zero after the 17th divide

Code: [Select]

{$MODESWITCH NESTEDCOMMENTS+}

program speed_test;
uses sysutils
, strutils
, strings
, math
, Multi_Int
;

var

big_int_size,
start_time,
end_time :int32;
delta :double;

MV_i,
MV_k :Multi_Int_XV;
n:longint;
s:ansistring;

begin
big_int_size:= 1876;
Multi_Int_Initialisation(big_int_size);

s:='22800589663853877505317123043453362340893519704710196572342276129978';
s+='5396584763579883599204575708254300465107640442948149273572670233883';
s+='7284359538516468653306625429160207551503354344345787111445922584002';
s+='8181165619959806798093081539066769221754339703446744180324866502719';
s+='0583529084870239014806639408757137138171908408383470514328939737730';
s+='7786595490401129699001354712004975431883261644867724606315936597415';
s+='4172816638070902792312674816613501437397527507108784136297892120724';
s+='2258609187180759380286550432287007389509166022499552376394375731290';
s+='3695257424344652576833807767548298955014117685289979591337921315574';
s+='4143851718821477323546913115004666541473134847568172342426632844406';
s+='6932325217707899195155340149974932811967498981615287349348710134095';
s+='3136158397898370157895394723520831996509093434513075798856068558092';
s+='5418302202627073384336628756554287325388628135750737038028444664630';
s+='9660903249431243633032650772924592366776659937552528684543292162042';
s+='0537820044936663363199059424287956068270893800711294561314155959675';
s+='9088339153782291900491510646687410015510285225141061882566364889308';
s+='439408424494044324417996846046572131526273';

MV_k:=s; //sqrt(sqrt(sqrt(sqrt(sqrt(factorial 10000)))))

start_time:= GetTickCount64;

MV_i:=1;
for n:=2 to 10000 do
begin
MV_i:=MV_i*n;
end;
end_time:= GetTickCount64;
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');

start_time:= GetTickCount64;
for n:=1 to 16 do
begin
MV_i:=MV_i div MV_k;
end;
end_time:= GetTickCount64;

writeln;
write('should be 1 but it''s ');
writeln(MV_i.ToStr);
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');

end.


[ad1mt]

I've loaded an updated version 4.39 onto Github and my webpage (at 18:40 UK time).
This fixes another bug I just found while testing.
Note that the unit test program I made available earlier, was also updated at the same time.
If you downloaded version 4.39 prior to 18:40 UK time, please download again.
N.B. this does not fix the 2nd bug reported by srvaldez...

the following silly test sets MV_i to factorial 10000, then in the following loop it divides MV_i by the sqrt(sqrt(sqrt(sqrt(sqrt(factorial 10000))))) 32 times, the result should be 1 but it's zero after the 17th divide

Thanks srvaldez for your help testing; I think you are the first.
I will fix this bug ASAP...

[ad1mt]

I like the fact that you use what I call "half-word" types, i.e. 16-bit on 32-bit CPU, and 32-bit on 64-bit CPU.
This is what I did for my RSA calculation unit, and even in "pure pascal", the FPC compiler gives good performance.
https://github.com/synopse/mORMot2/blob/master/src/crypt/mormot.crypt.rsa.pas#L57

I've downloaded the source from Github, but I can't find your code for big integer type definitions and calculations.
I would like to study your division algorithm, and maybe steal some ideas... 
Thanks.

[BeniBela]

I like the fact that you use what I call "half-word" types, i.e. 16-bit on 32-bit CPU, and 32-bit on 64-bit CPU.

I did not really like having to do this, but I could not think of a better/faster method.
Another method would be to use the full size of 32/64 bit words, but I would then have to detect overflow, and then when overflow was detected, redo the calculation again in a safe way.
I assumed that the second method would be slower, although I did not test this, and maybe I should have. However, even if the second method was faster, it would have made the code significantly more complex, and it was complex enough already.

Well, it would be easy in assembly

mul already gives a 128-bit result for a 64-bit multiplication

[ad1mt]

Well, it would be easy in assembly
mul already gives a 128-bit result for a 64-bit multiplication

Yes, but...
1. I wanted maximum portability across all CPU's/word-sizes.
2. A long time ago, I wrote a large program in assembly, and I would never want to program in assembly again!

[srvaldez]

hello ad1mt
there's a FreeBasic implementation of BigInt including division https://github.com/stephanbrunker/big_integer/blob/master/modules/arithfunctions.bi