(a * b) mod c with int64

[circular]

Try with this unit I've just made :

[KpjComp]

@circular,

Nice unit!!,  I like the use of operator overloading there.

Implementing two's compliment with int128's, would be fun too.  ,

Using the same technique, a version that's N'bit might come in handy.  It could then maybe be used for cross platform X Bit cryptography routines.  Here the unsigned/signed could maybe be part of the structure too.

[circular]

Thanks Kpj.

It is easy to add a sign, for example by creating a SInt128 (signed int 128) record that contains the UInt128 value and a "negative" boolean. Then, it is rather easy to do operators for this new structure.

There is already two's complement. For example if you write a - b, b is transformed this way. But well, there is no operator for just -b. It can be easily added by doing the same as in a-b, which is apply "not" and then inc.

[BeniBela]

Try with this unit I've just made :

Great, it works.

(even better than the asm version, since that seems to cause SIGFPE for a*b results larger than 2^100 )

But is it really necessary to do that much shifting in the division function?

[KpjComp]

But is it really necessary to do that much shifting in the division function?

Well @circular has kind of emulated what a CPU does, when you call IMUL and friends, internally the CPU is doing all this Bit Shifting for you, and explains why Multiply and Division often take a fair few CPU cycles.  The maximum number of shifts is then of course proportional to the number bits, but at least is not proportional to the potential size of number you can make.

[circular]

Shifting is a fast operation. There is not that much shift in fact. Just one way and the other.

It may be optimized by computing the upper bit position, but I'm satisfied that it works and gives the right result. And I don't have so much time to spend on the optimization.

But as I already said, it should be already quite fast. But you have not answered the question of KpjComp :

Is speed important?

Yes KpjComp, it's just like multiply instruction.

[BeniBela]

 But you have not answered the question of KpjComp :

Is speed important?

Honestly, I don't need it at all anymore.

I just wanted to solve a project euler problem in which you have to test 50 million numbers for primality,  (which requires around 100 of these modulo multiplication for each Rabin Miller test). But now I have solved it with the asm version...

But perhaps it is important for some later euler problems

[circular]

lol ok
Well don't count on me to optimize it then. 

[circular]

Hello, I read this topic and download int128.pas example! It's very good implementation!
Please help me, how can i use biggest number, such as >1024 bits?
Thanks.

@UserOfLazarus;
I would say that the first step would be to remove the shortcuts dw0, dw1 etc. and use only array notation. Then to transform with for loops and using a constant for the number of dwords in the number. For 128 bits, it is 128/32 = 4.

After this transformation is complete, you could change the constant. For example, for 1024 bits, that would be 1024/32=32.

[guest60499]

It seems like someone has created a unit for libgmp bindings: http://wiki.freepascal.org/gmp.

[useroflazarus]

modified version of int128 by circular 

notes:
1) 0...N bits 
2) added two functions of modular arithmetic (PowMod and ModInverse)
3) example of generate probably prime numbers (use Fermat's little theorem)

[circular]

Good attempt.

I guess that the multiply function is not working. I would suggest to do something like:

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

1. operator * (Value1: TLongW;  Value2: TLongW): TLongW;
2. var
3.   buf : qword;
4.   carry : longword;
5.   i, j : Integer;
6.   ctemp : TLongW;
7. begin
8.  
9.  result := 0;
10.  ctemp := 0;
11.  
12.  buf := 0;
13.  
14.  ctemp := 0;
15.  for i:=0 to size1-1 do begin
16.    if Value2.dw[j] <> 0 then
17.    for j:=0 to size2-1 do
18.      if Value1.dw[i]<>0 then begin
19.      buf := qword(Value1.dw[i]) * qword(Value2.dw[j]);
20.      carry := buf shr 32;
21.  
22.      ctemp.dw[i+j] := buf and $ffffffff;
23.      if i+j+1 <= MAXSIZE then ctemp.dw[i+j+1] := carry;
24.  
25.      result += ctemp;
26.  
27.      ctemp.dw[i+j] := 0;
28.      if i+j+1 <= MAXSIZE then ctemp.dw[i+j+1] := 0;
29.    end;
30.  end;
31. end;

Also, when checking the shift parameter of shr/shl, using "shift.dw[0] >= (MAXSIZE+1)*32" will not tell you if there is something in dw[1], etc.

[useroflazarus]

Also, when checking the shift parameter of shr/shl, using "shift.dw[0] >= (MAXSIZE+1)*32" will not tell you if there is something in dw[1], etc.

for example 1024 (or 2048 bits number) shift left is not big (dw[0] maximum is 2^32 = 4294967295) and 1024 much less than 4294967295.
4294967295 bits is very big for operands, so don't have to worry, or am I wrong?

[circular]

You don't have to worry of course.

Yet for now you would have:

1 shl 4294967296 = 1
1 shr 4294967296 = 1

Every detail count in my opinion. However, the problem I highlight about the multiplication is more important.

[useroflazarus]

slightly corrected multiplication.

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

1. type
2.   tq = packed record
3.     case boolean of
4.     true:  (qw : qword);
5.     false: (lw, hw : longword);
6.   end;
7.  
8. operator * (Value1: TLongW;  Value2: TLongW): TLongW;
9. var
10.   buf: tq;
11.   i, j : integer;
12.   size1, size2 : longword;
13. begin
14.  
15.  result := 0;
16.  size1 := GetSize(Value1);
17.  size2 := GetSize(Value2);
18.  
19.  for i:=0 to size1-1 do begin
20.    buf.qw := 0;
21.    for j:=0 to size2-1 do begin
22.       buf.qw := qword(Value1.dw[i]) * qword(Value2.dw[j]) + qword(buf.hw) + qword(result.dw[i+j]);
23.       result.dw[i+j] := buf.lw;
24.    end;
25.    result.dw[i+j+1] := buf.hw;
26.  end;
27. end;
28.  
29.  

i suppose, that better variant of multiplication is https://en.wikipedia.org/wiki/Karatsuba_algorithm