Recent

Author Topic: New Big Integer library is almost finished  (Read 18378 times)

ad1mt

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #15 on: September 30, 2024, 10:05:21 am »
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

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #16 on: September 30, 2024, 10:23:55 am »
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

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #17 on: September 30, 2024, 10:25:30 am »
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

  • Full Member
  • ***
  • Posts: 117
Re: New Big Integer library is finished
« Reply #18 on: September 30, 2024, 12:18:11 pm »
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.
« Last Edit: September 30, 2024, 12:41:41 pm by srvaldez »

ad1mt

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #19 on: October 01, 2024, 09:09:35 am »
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

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #20 on: October 03, 2024, 06:33:24 pm »
It has become clear that my Big Integer library is not finished...  :(

ad1mt

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #21 on: October 08, 2024, 10:12:45 am »
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

  • Sr. Member
  • ****
  • Posts: 405
Re: New Big Integer library is finished
« Reply #22 on: October 08, 2024, 11:35:52 am »

...

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

  • Full Member
  • ***
  • Posts: 117
Re: New Big Integer library is finished
« Reply #23 on: October 08, 2024, 02:38:59 pm »
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

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #24 on: October 08, 2024, 03:22:43 pm »
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!  :)
« Last Edit: October 08, 2024, 04:40:15 pm by ad1mt »

ad1mt

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #25 on: October 08, 2024, 03:23:18 pm »
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

  • Hero Member
  • *****
  • Posts: 5468
    • Bart en Mariska's Webstek
Re: New Big Integer library is finished
« Reply #26 on: October 08, 2024, 04:20:22 pm »
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

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #27 on: October 08, 2024, 04:26:29 pm »
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.
« Last Edit: October 08, 2024, 04:37:37 pm by ad1mt »

ad1mt

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #28 on: October 08, 2024, 04:29:06 pm »
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

  • Sr. Member
  • ****
  • Posts: 327
    • Mark Taylor's Home Page
Re: New Big Integer library is finished
« Reply #29 on: October 08, 2024, 04:44:56 pm »
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.

 

TinyPortal © 2005-2018