Anyone interested in testing a new Big Integer library?

[PascalDragon]

I am now wondering if I need to code a 16bit option!   

Well, FPC does support compiling for i8086, AVR and Z80.

[ad1mt]

I've been trying to take Henry Warren's division algorithm in C, and code it in Pascal.

Here's a web page that discusses the C code, and variations of it...
  https://skanthak.homepage.t-online.de/division.html

I have made much progress, but I'm still having problems. My head is starting to spin, because my Pascal code is behaving differently to how I imagine the C code should work, in ways I don't understand. In some parts of the code, I cannot make sense of what the C code is trying to do, so have no idea what is going wrong.

Is there anyone out there who is interested in debugging extremely complex algorithms, who would be interested helping me crack this problem?

If I can get it working, I'm going to use it in my Pascal big integer library, instead of my own mediocre division algorithm I'm currently using. Once the library/unit is finished & tested, I intend to contribute it to the Free Pascal community.
Thanks.

PS...
1. there is a variation of the algorithm in the Free Pascal RTL, but it is too specialized for me to use easily.
2. I advise caution if taking the C code directly from the above web site, it is different to the code in the "Hacker's Delight" book I am working from.
3. You could argue that if I cannot solve/crack this problem myself, then I am being too ambitious in trying to create a big integer library! And I would struggle to argue against that!   
4. Oh, how I hate C!  This problem exemplifies everything I hate about that abominable language! It creates the illusion of a safe, high-level language, but in reality it is just bit of superficial gloss on top of an assembly language. Almost every security bug in Unix/Linux over the last 30 years can be directly attributed to bad design of the C language and its run-time libraries.

[cpicanco]

I've been trying to take Henry Warren's division algorithm in C, and code it in Pascal.

Here's a web page that discusses the C code, and variations of it...
  https://skanthak.homepage.t-online.de/division.html

I have made much progress, but I'm still having problems. My head is starting to spin, because my Pascal code is behaving differently to how I imagine the C code should work, in ways I don't understand. In some parts of the code, I cannot make sense of what the C code is trying to do, so have no idea what is going wrong.

Is there anyone out there who is interested in debugging extremely complex algorithms, who would be interested helping me crack this problem?

If I can get it working, I'm going to use it in my Pascal big integer library, instead of my own mediocre division algorithm I'm currently using. Once the library/unit is finished & tested, I intend to contribute it to the Free Pascal community.
Thanks.

PS...
1. there is a variation of the algorithm in the Free Pascal RTL, but it is too specialized for me to use easily.
2. I advise caution if taking the C code directly from the above web site, it is different to the code in the "Hacker's Delight" book I am working from.
3. You could argue that if I cannot solve/crack this problem myself, then I am being too ambitious in trying to create a big integer library! And I would struggle to argue against that!   
4. Oh, how I hate C!  This problem exemplifies everything I hate about that abominable language! It creates the illusion of a safe, high-level language, but in reality it is just bit of superficial gloss on top of an assembly language. Almost every security bug in Unix/Linux over the last 30 years can be directly attributed to bad design of the C language and its run-time libraries.

Hi ad1mt, have you tried ChatGPT (or GitHub Copilot)? It is really good in translating, describing and summarization of well defined contexts. Of course, you must have the expertise (which is not my case in this particular problem right now) to guide and evaluate ChatGPT's answers (or Copilot completions). Take a look at this example in attachment. I used the link you provided as a context and ask him to translate the ANSI C optimized algorithm.

[MathMan]

I've been trying to take Henry Warren's division algorithm in C, and code it in Pascal.

Here's a web page that discusses the C code, and variations of it...
  https://skanthak.homepage.t-online.de/division.html

I have made much progress, but I'm still having problems. My head is starting to spin, because my Pascal code is behaving differently to how I imagine the C code should work, in ways I don't understand. In some parts of the code, I cannot make sense of what the C code is trying to do, so have no idea what is going wrong.

Is there anyone out there who is interested in debugging extremely complex algorithms, who would be interested helping me crack this problem?

If I can get it working, I'm going to use it in my Pascal big integer library, instead of my own mediocre division algorithm I'm currently using. Once the library/unit is finished & tested, I intend to contribute it to the Free Pascal community.
Thanks.

PS...
1. there is a variation of the algorithm in the Free Pascal RTL, but it is too specialized for me to use easily.
2. I advise caution if taking the C code directly from the above web site, it is different to the code in the "Hacker's Delight" book I am working from.
3. You could argue that if I cannot solve/crack this problem myself, then I am being too ambitious in trying to create a big integer library! And I would struggle to argue against that!   
...

Some feedback on your trial from my side

* I wouldn't dig into that too heavily. This is "only" a fixed size division of a double 'limb' by a single 'limb' value. Inside the algorithms they still use div & mod, which is not optimal <= 'limb' here stands for your elementary arithmetic unit - in your case either a 32 or 64 bit unsigned integer. The term was coined by Knuth and you may have stumbled across it already.

* implementations of long division that I know usually replace the actual division by a multiplication with an inverse. It actually isn't that complicated, once you understood the Algebra behind it. Depending on the size of the divisor (1 or >1) there are two types of inverses used and there is some literature available how to calculate these.

* I could provide a reference implementation of long division in Pascal. But this is tailored to a completely different environment - so you still would have a substantial amount of integration ahead of you. And I'm not going to do that for you 

* or I could try to guide you through the logic, leaving the implementation to you. Fair warning - you'd still have to read and understand some papers for this.

Cheers,
MathMan

[ad1mt]

Hi ad1mt, have you tried ChatGPT (or GitHub Copilot)? It is really good in translating, describing and summarization of well defined contexts. Of course, you must have the expertise (which is not my case in this particular problem right now) to guide and evaluate ChatGPT's answers (or Copilot completions). Take a look at this example in attachment. I used the link you provided as a context and ask him to translate the ANSI C optimized algorithm.

Thanks for the suggestion... I've heard of people using ChatGPT to write code, but I've never tried it myself.
Did you actually run the example code you sent me?
I fed it two "problematic" numbers I'm using as a test case while trying to debug my code...

u1:= 9223372036854775808; u2:= 0; v1:= 2305843009213693948;

the pascal div operator gives the correct result...
u1 div v1 = 4 remainder 16

The A.I. generated DivLLU function gives...
Quotient: 3
Remainder: 12

The code appears to be an A.I. "halucination"...

[ad1mt]

Thanks for your reply.

* I wouldn't dig into that too heavily. This is "only" a fixed size division of a double 'limb' by a single 'limb' value. Inside the algorithms they still use div & mod, which is not optimal

Are you referring to the Warren (Knuth) "algorithm D" here?  The reason I was using this algorithm, was that it appeared to be a reasonable balance between complexity/size of code and speed of execution, and therefore better than my own home-made division algorithm.
An algorithm might be "the best", but if I'm unable to get the code working, my library will not get finished.

* implementations of long division that I know usually replace the actual division by a multiplication with an inverse. It actually isn't that complicated, once you understood the Algebra behind it. Depending on the size of the divisor (1 or >1) there are two types of inverses used and there is some literature available how to calculate these.

Yes, I've read about this idea while researching. However, I can't visualise how to make inverses work with integers. Is this solved by "scaling" everything so that all the inverses become whole numbers?
Another thing I can't understand is, how do you deal with inverses that are an infinite series?  Truncating or rounding would defeat the purpose of the integer arithmetic. Do you solve this by some kind of "lazy evaluation" that cancels out all the scaling factors in the final calculations?

* I could provide a reference implementation of long division in Pascal. But this is tailored to a completely different environment - so you still would have a substantial amount of integration ahead of you. And I'm not going to do that for you 

Thanks for your offer, I would be very grateful if you could do that.
I have my doubts about being able to understand the code... I've already looked at the division functions inside the MPArith and FNX libraries, and the code seems incredibly complex.
However, please do send the reference code, and I will take a look.  Which algorithm does your reference code use?  And how many lines of code are there?

[MathMan]

Thanks for your reply.

* I wouldn't dig into that too heavily. This is "only" a fixed size division of a double 'limb' by a single 'limb' value. Inside the algorithms they still use div & mod, which is not optimal

Are you referring to the Warren (Knuth) "algorithm D" here?  The reason I was using this algorithm, was that it appeared to be a reasonable balance between complexity/size of code and speed of execution, and therefore better than my own home-made division algorithm.
An algorithm might be "the best", but if I'm unable to get the code working, my library will not get finished.

Fair enough - if you tackle this as a learning experience, then you should stick to what you understand and can explain.

* implementations of long division that I know usually replace the actual division by a multiplication with an inverse. It actually isn't that complicated, once you understood the Algebra behind it. Depending on the size of the divisor (1 or >1) there are two types of inverses used and there is some literature available how to calculate these.

Yes, I've read about this idea while researching. However, I can't visualise how to make inverses work with integers. Is this solved by "scaling" everything so that all the inverses become whole numbers?
Another thing I can't understand is, how do you deal with inverses that are an infinite series?  Truncating or rounding would defeat the purpose of the integer arithmetic. Do you solve this by some kind of "lazy evaluation" that cancels out all the scaling factors in the final calculations?

Regarding you last point, yes and no  You calculate a single 'limb' sized inverse and with that you estimate the next quotient 'limb' via multiplication. It can be shown that this estimated quotient 'limb' is either exact or one to large. So you have to check that subtracting the Dvsr times the estimated quotient limb from the Dvd generates a carry - in which case the estimation was one to large and you have to readjust the remaining Dvd by adding back the Dvsr and decrement the quotient estimation. This way all stays in integers, no infinite series or late evaluations are required - just some standard post-processing.

* I could provide a reference implementation of long division in Pascal. But this is tailored to a completely different environment - so you still would have a substantial amount of integration ahead of you. And I'm not going to do that for you 

Thanks for your offer, I would be very grateful if you could do that.
I have my doubts about being able to understand the code... I've already looked at the division functions inside the MPArith and FNX libraries, and the code seems incredibly complex.
However, please do send the reference code, and I will take a look.  Which algorithm does your reference code use?  And how many lines of code are there?

The algorithm used is just an optimised implementation of Knuth "Algorithm D". I don't have the exact LoC figure handy, as there are several helper functions required beside the main routine.- shot from the hip ~500. I'll prepare a package - will take a bit - and post it here.

[cpicanco]

Hi ad1mt, have you tried ChatGPT (or GitHub Copilot)? It is really good in translating, describing and summarization of well defined contexts. Of course, you must have the expertise (which is not my case in this particular problem right now) to guide and evaluate ChatGPT's answers (or Copilot completions). Take a look at this example in attachment. I used the link you provided as a context and ask him to translate the ANSI C optimized algorithm.

Thanks for the suggestion... I've heard of people using ChatGPT to write code, but I've never tried it myself.
Did you actually run the example code you sent me?
I fed it two "problematic" numbers I'm using as a test case while trying to debug my code...

u1:= 9223372036854775808; u2:= 0; v1:= 2305843009213693948;

the pascal div operator gives the correct result...
u1 div v1 = 4 remainder 16

The A.I. generated DivLLU function gives...
Quotient: 3
Remainder: 12
f
The code appears to be an A.I. "halucination"...

Yes, I did some tests just to compile it. As I said, you should use it for translations (granular questions as much as possible), for templating, and auto-completion of a well defined (algoritic) line of thought. However, as I said, I do not have the expertise to validate the specific case here. So, I am claiming that, with your expertise and vocabulary, AI will be much powerfull. With knownledge, you can elaborate better and granular questions and proper evaluate the results. It can save you a lot of time saving.

[cpicanco]

Hi ad1mt, have you tried ChatGPT (or GitHub Copilot)? It is really good in translating, describing and summarization of well defined contexts. Of course, you must have the expertise (which is not my case in this particular problem right now) to guide and evaluate ChatGPT's answers (or Copilot completions). Take a look at this example in attachment. I used the link you provided as a context and ask him to translate the ANSI C optimized algorithm.

Thanks for the suggestion... I've heard of people using ChatGPT to write code, but I've never tried it myself.
Did you actually run the example code you sent me?
I fed it two "problematic" numbers I'm using as a test case while trying to debug my code...

u1:= 9223372036854775808; u2:= 0; v1:= 2305843009213693948;

the pascal div operator gives the correct result...
u1 div v1 = 4 remainder 16

The A.I. generated DivLLU function gives...
Quotient: 3
Remainder: 12
f
The code appears to be an A.I. "halucination"...

Yes, I did some tests just to compile it. As I said, you should use it for translations (granular questions as much as possible), for templating, and auto-completion of a well defined (algorithmic?) and consistent line of thought. However, as I said, I do not have the expertise to validate the specific case here. So, I am claiming that, with your expertise and vocabulary, AI will be much powerfull. With knownledge, you can elaborate better and granular questions and proper evaluate the results. It can save you a lot of time.

[ad1mt]

Fair enough - if you tackle this as a learning experience, then you should stick to what you understand and can explain.

I've found that if I type in code and it doesn't work immediately, I then have to understand how it works in order to figure out why its not working.  Trying and failing several times to get the Warren/Knuth code to work has been a real eye-opener for me. This is the first time I've ever worked with a piece of code that I was unable to understand or get working; and its not very many lines of code.

The algorithm used is just an optimised implementation of Knuth "Algorithm D". I don't have the exact LoC figure handy, as there are several helper functions required beside the main routine.- shot from the hip ~500. I'll prepare a package - will take a bit - and post it here.

Thanks for your help.
I'm curious to know where your code originates from...or did you design it yourself?

[cpicanco]

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

1. // Copyleft © 2011-2023, Stefan Kanthak <stefan.kanthak@nexgo.de>
2.  
3. // Divide a 128-bit integer dividend, supplied as a pair of 64-bit
4. // integers in u0 and u1, by a 64-bit integer divisor, supplied in v;
5. // return the 64-bit quotient and optionally the 64-bit remainder in *r.
6. unit division;
7.  
8. interface
9.  
10. function DivLLU(u0, u1, v: QWord; var r: QWord): QWord;
11.  
12. implementation
13.  
14. function DivLLU(u0, u1, v: QWord; var r: QWord): QWord;
15. var
16.   qhat, rhat, uhat: QWord;
17.   q0, q1: Cardinal;
18.   s: Byte;
19. begin
20.   if u1 >= v then
21.   begin
22.     if @r <> nil then
23.       r := High(QWord);
24.     Exit(High(QWord));
25.   end;
26.  
27.   s := 63-BsrQWord(v);
28.   if s <> 0 then
29.   begin
30.     v := v shl s;
31.     u1 := u1 shl s;
32.     u1 := u1 or (u0 shr (64 - s));
33.     u0 := u0 shl s;
34.   end;
35.  
36.   qhat := u1 div (v shr 32);
37.   rhat := u1 mod (v shr 32);
38.  
39.   while ((qhat shr 32) <> 0) or
40.         ((qhat and High(Cardinal)) * (v and High(Cardinal)) >
41.          ((rhat shl 32) or (u0 shr 32))) do
42.   begin
43.     Dec(qhat);
44.     Inc(rhat, v shr 32);
45.     if (rhat shr 32) <> 0 then Break;
46.   end;
47.  
48.   q1 := qhat and High(Cardinal);
49.  
50.   uhat := ((u1 shl 32) or (u0 shr 32)) - q1 * v;
51.  
52.   qhat := uhat div (v shr 32);
53.   rhat := uhat mod (v shr 32);
54.  
55.   while ((qhat shr 32) <> 0) or
56.         ((qhat and High(Cardinal)) * (v and High(Cardinal)) >
57.          ((rhat shl 32) or (u0 and High(Cardinal)))) do
58.   begin
59.     Dec(qhat);
60.     Inc(rhat, v shr 32);
61.     if (rhat shr 32) <> 0 then Break;
62.   end;
63.  
64.   q0 := qhat and High(Cardinal);
65.  
66.   if @r <> nil then
67.     r := (((uhat shl 32) or (u0 and High(Cardinal))) - q0 * v) shr s;
68.  
69.   Result := (QWord(q1) shl 32) or q0;
70. end;
71.  
72. end.
73.  
74.  

So, now it works with basic stuff thanks to ChatGPT 4 (through Microsoft Edge's Copilot).

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

1. program algorithmd;
2.  
3. uses division;
4.  
5. procedure RunTest(ATestName: string; const u0, u1, v: QWord);
6. var
7.   Quotient, ExpectedQuotient, ExpectedRemainder: QWord;
8.   Remainder : QWord = 0;
9. begin
10.   writeln(ATestName);
11.   writeln('u0: ', u0);
12.   writeln('u1: ', u1);
13.   writeln('v: ', v);
14.  
15.   Quotient := DivLLU(u0, u1, v, Remainder);
16.  
17.   if v = 0 then begin
18.     ExpectedQuotient := 0;
19.     ExpectedRemainder := 0;
20.   end else begin
21.     ExpectedQuotient := u0 div v;
22.     ExpectedRemainder := u0 mod v;
23.   end;
24.  
25.   writeln('Quotient: ', quotient, ' Expected: ', ExpectedQuotient);
26.   writeln('Remainder: ', remainder, ' Expected: ', ExpectedRemainder);
27.  
28.   if (Quotient = ExpectedQuotient) and (Remainder = ExpectedRemainder) then
29.     writeln('Test Passed')
30.   else
31.     writeln('Test Failed');
32.  
33.   writeln;
34. end;
35.  
36.  
37.  
38. begin
39.   RunTest('// Hard case',
40.     9223372036854775808, 0, 2305843009213693948);
41.   RunTest('// Quotient should be 5, Remainder should be 0',
42.     10, 0, 2);
43.   RunTest('// Quotient should be 33, Remainder should be 1',
44.     100, 0, 3);
45.   RunTest('// Quotient and Remainder should be calculated',
46.     1234567890123456789, 0, 9876543210);
47.   RunTest('// Quotient should be 1, Remainder should be 0',
48.     High(QWord), 0, High(QWord));
49.   RunTest('// Quotient should be High(QWord), Remainder should be 0',
50.     High(QWord), 0, 1);
51.  
52.   WriteLn('Test Basic 128-bit Division');
53.   RunTest('// Quotient should be 9223372036854775808 ??, Remainder should be 0',
54.     0, 1, 2);
55.  
56.   RunTest('// Should return valid quotient and remainder',
57.     High(QWord), High(QWord), High(QWord));
58.  
59.   RunTest('// Should return valid quotient and remainder',
60.     Random($FFFFFFFFFFFFFFFF),
61.     Random($FFFFFFFFFFFFFFFF),
62.     Random($FFFFFFFFFFFFFFFF));
63.  
64.   ReadLn; // Wait for user input before closing
65. end.
66.  

[MathMan]

The algorithm used is just an optimised implementation of Knuth "Algorithm D". I don't have the exact LoC figure handy, as there are several helper functions required beside the main routine.- shot from the hip ~500. I'll prepare a package - will take a bit - and post it here.

Thanks for your help.
I'm curious to know where your code originates from...or did you design it yourself?

Here we go, see attached. I think it will be easiest to start with the function 'lDivMod' in 'LongDivMod.pas' and drill your way down to understand how things work.

'LongDivMod.pas' contains everything that is required for the long division with remainder. I heavily commented everything to help understanding. I was correct with my estimation of 500 LoC in total. I added a unit test, based on 'fpunittest', to verify the implementation - maybe you can use / extend that for your own library functions. I also added a benchmark utility that lets you bench each function used in 'LongDivMod.pas', which may help you in detecting regressions while you try to integrate into your own lib. Both, the unit test and the benchmark, are command line utilities.

The implementation I did myself. But of course I researched into algorithms and looked left and right how others approached the same topic.

Cheers,
MathMan

[Thaddy]

the pascal div operator gives the correct result...
u1 div v1 = 4 remainder 16

The A.I. generated DivLLU function gives...
Quotient: 3
Remainder: 12

The code appears to be an A.I. "halucination"...

No, because you did not specify Integer diivide to openAI chatGPT. I got a almost correct answer,
In general you need to provide as much context as you can and then the answer is more likely to be close to correct.
There is a clear bias towards C like notation.

[ad1mt]

'LongDivMod.pas' contains everything that is required for the long division with remainder.

The implementation I did myself. But of course I researched into algorithms and looked left and right how others approached the same topic.

Cheers,
MathMan

Many thanks for this... I will take a look shortly.
I'm taking a short break from division, to switch to a different problem, which you might also be able to help me with...
How to reliably convert from real/float types to large integer. I've discovered to my dismay that this is not trivial. I will shortly ask for help in a different thread.

[MathMan]

How to reliably convert from real/float types to large integer. I've discovered to my dismay that this is not trivial. I've asked for help in a different thread.

It would be helpful if you provide a link to that thread, cause I'm unable to locate it.