New Big Integer library is almost finished

[ad1mt]

hello ad1mt, silly tests often reveal bugs

There are no such things as silly tests, just silly bugs!
I've uploaded version 4.40, which fixes more bugs.
Many thanks for your feedback.

[ad1mt]

the following 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

There is a bug in your code. Factorial 10000 is not a perfect square, so every time round the division loop, a remainder is lost, and the result becomes zero more quickly than you would expect.
If you calculate MV_i:= MV_k ** 32; instead of factorial 10000, then the division loop works as expected.

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

1. program bug_report_2;
2. {$MODESWITCH NESTEDCOMMENTS+}
3. uses    sysutils
4. ,               strutils
5. ,               strings
6. ,               math
7. ,               Multi_Int
8. ;
9. var
10.         big_int_size,
11.         start_time,
12.         end_time                :int32;
13.         delta                   :double;
14.         MV_i,
15.         MV_f,
16.         MV_k,
17.         MV_e            :Multi_Int_XV;
18.         n                       :longint;
19.         s                       :ansistring;
20.         rem                     :boolean;
21. begin
22. big_int_size:= 4000;
23. Multi_Int_Initialisation(big_int_size);
24.  
25. s:='22800589663853877505317123043453362340893519704710196572342276129978';
26. s+='5396584763579883599204575708254300465107640442948149273572670233883';
27. s+='7284359538516468653306625429160207551503354344345787111445922584002';
28. s+='8181165619959806798093081539066769221754339703446744180324866502719';
29. s+='0583529084870239014806639408757137138171908408383470514328939737730';
30. s+='7786595490401129699001354712004975431883261644867724606315936597415';
31. s+='4172816638070902792312674816613501437397527507108784136297892120724';
32. s+='2258609187180759380286550432287007389509166022499552376394375731290';
33. s+='3695257424344652576833807767548298955014117685289979591337921315574';
34. s+='4143851718821477323546913115004666541473134847568172342426632844406';
35. s+='6932325217707899195155340149974932811967498981615287349348710134095';
36. s+='3136158397898370157895394723520831996509093434513075798856068558092';
37. s+='5418302202627073384336628756554287325388628135750737038028444664630';
38. s+='9660903249431243633032650772924592366776659937552528684543292162042';
39. s+='0537820044936663363199059424287956068270893800711294561314155959675';
40. s+='9088339153782291900491510646687410015510285225141061882566364889308';
41. s+='439408424494044324417996846046572131526273';
42.  
43. MV_k:=s;
44.  
45. MV_i:= MV_k ** 32;   // factorial 10000 is a different value to (MV_k ** 32)
46.  
47. writeln('MV_i = ',MV_i.ToStr);
48.  
49. n:= 32;
50. while (n > 0) do
51.         begin
52.         MV_i:=MV_i div MV_k;
53.         Dec(n);
54.         end;
55.  
56. writeln(n,' MV_i = ',MV_i.tostr);
57. end.
58.  

Many thanks for your feedback, please keep on testing my library, if you have the time.

[ad1mt]

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

Thanks... I will have a look at it.

[srvaldez]

ad1mt
I know that it's not a perfect square and I disagree that it's a bug in my code, the division should round toward zero, try the following using the Pari/gp calculator http://pari.math.u-bordeaux.fr/download.html

k=sqrtint(sqrtint(sqrtint(sqrtint(sqrtint(10000!)))));
f=10000!;
for(i=1,32,f=f\k);
f

sqrtint = integer square root
f \ k = integer division

or using Mparith https://forum.lazarus.freepascal.org/index.php/topic,68447.msg528998.html#msg528998

Code: [Select]

program mp_test;

uses
mp_types, mp_base, mp_numth, mp_calc, sysutils, windows, math;

var
n,
start_time,
end_time :int32;
delta :double;

i, k:mp_int;
s:ansistring;

begin
mp_init(k);

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';
mp_read_decimal_astr(k, s);

start_time:= GetTickCount64;
mp_init(i);
mp_fact(10000, i);
end_time:= GetTickCount64;
writeln('Factorial(10000) = ', mp_adecimal(i));
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');

start_time:= GetTickCount64;
for n:=1 to 32 do
begin
mp_div(i, k, i);
end;
end_time:= GetTickCount64;
writeln('Factorial(1) = ', mp_adecimal(i));
delta:= (end_time - start_time) / 1000;
writeln('time elapsed is ', delta:1:6, ' seconds');

mp_clear(k);
mp_clear(i);
end.


[ad1mt]

I disagree that it's a bug in my code.

Sorry... I was mistaken.
With the latest bug fixes, my code produces the same results as mp_arith.
Although significantly more slowly   

[ad1mt]

It has become clear that my Big Integer library is not finished... 

[ad1mt]

sqrt(sqrt(sqrt(sqrt(sqrt(factorial 10000)))))

Code: [Select]

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)))))


Can I ask where you got that value for sqrt(sqrt(sqrt(sqrt(sqrt(factorial 10000)))))?
Thanks.

[MathMan]

...

Can I ask where you got that value for sqrt(sqrt(sqrt(sqrt(sqrt(factorial 10000)))))?
Thanks.

I would assume that srvaldez got it from Mparith or Pari - at least it is correct.

Cheers,
MathMan

[srvaldez]

I think that MathMan is right, probably pari, don't remember now 

ad1mt, now that you got the lib working you can fine tune it for performance 

[ad1mt]

I think that MathMan is right, probably pari, don't remember now 

ad1mt, now that you got the lib working you can fine tune it for performance 

Unfortunately, there might still be bugs... my result for the square root of !10000 is different to the mp_arith value; so that needs investigating. It might be to do with whether the result is rounded up or rounded down; I am rounding down with a remainder.

I'm also looking at speeding up the library. In the unreleased version I'm working on, I've implemented the Babylonian algorithm for the square root function (replacing my own home-grown version), and its now 20x faster.

I've looked at other division algorithms, but I either can't understand them, or they seem incredibly complex (e.g mp_arith) and would therefore be difficult to re-code. I won't attempt to do a "blind" line-by-line copy of code again without understanding it. I tried to do that with the Warren-Knuth C language division algorithm, and when it didn't work, I could not figure out why it was failing, or how to fix it. I was only able to get it working by "re-inventing" the algorithm using pencil & paper; and then I was able to re-code it in Pascal.

My library was supposed to be basic and easy-to-use, with maximum portability. It is not intended to be the fastest... I'm happy to leave that award to other libraries.

I don't want to spend a lot more time on it. I've got other things I'd lke to work on in my remaining years.
Feel free to help if you like! 

[ad1mt]

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

Thanks for that tip... I'll have a look at it.

[Bart]

First of all: kudo's to you.

Wow, 16712 lines of code in a single unit!
That's hard to study.
Maybe it would be an idea to split the code uning include files, eg. 1 include for defining operators, 1 for conversions, 1 for arithmetic (maybe even spilt that: multiplication related (power, factorial, square), division related (divide, square root other roots), adding, subtracting, bit operations).
I also see lots of repetitive code (checking for "Defined" or e.g. overflow flag), that could be factored out.

Bart

[ad1mt]

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

Correct me if I'm mistaken... it looks like its using the same Warren-Knuth division algorithm I'm using.

[ad1mt]

you can fine tune it for performance 

Unfortunately, tuning the code will not make a significant difference... maybe around 10% max.
To get a substantial speedup (e.g. 5x or more faster), different algorithms are needed.

[ad1mt]

I also see lots of repetitive code (checking for "Defined" or e.g. overflow flag), that could be factored out.

By "factored out", do you mean made into a source code macro?
Or do you mean made into a separate procedure/function?
Separate procedures/functions would slow down the code.