var
DMulS:
function(R,A: Pointer; AC: Cardinal; B: Pointer; BC: Cardinal): Integer = _DMulS;
function DMulK(R,A: PDigits; AC: Integer; B: PDigits; BC: Integer; T: PDigits): Integer;
register;
// Karatsuba Mul, don't call directly use DMulN()
// R[] = A[] * B[], A are AC digits and B are BC digits long
// T temporary workspace must be 2(AC + BC) digits
// R must be (AC + BC) digits
// returns R.Count
// we use goto's because faster
var
I,J,K,A0,A1,R0,R1: Integer;
label
Null, Finish;
begin
NCalcCheck;
// first stage, remove leading zero's and recalc sizes,
I := AC + BC;
if A[AC -1] = 0
then AC := DNorm(0, A, AC -1);
if AC > 0
then
begin
if B[BC -1] = 0
then BC := DNorm(0, B, BC -1);
if BC = 0
then AC := 0;
end else BC := 0;
Result := AC + BC;
if Result < I
then
begin
while I > Result
do
begin
Dec(I);
R[I] := 0;
end;
// DFill(0, @R[Result], I - Result);
if Result = 0
then Exit;
end;
// swap such that A >= B
if AC < BC
then
begin
I := AC; AC := BC; BC := I;
Pointer(J) := A; A := B; B := Pointer(J);
end;
// smaller Mul possible ?
if BC < Mul_KBE
then
begin
Result := DMulS(R, A, AC, B, BC);
Exit;
end;
// second stage, do Karatsuba dependend on the sizes of A and B
if AC = BC
then
begin
A0 := AC
shr 1;
J := 0;
if not Odd(AC)
then
begin
I := DCmpN(A, @A[A0], A0);
if I = 1
then DSubN(R, A, @A[A0], A0)
else
if I = -1
then
begin
DSubN(R, @A[A0], A, A0);
J := 1;
end else DFill(0, R, A0);
I := DCmpN(B, @B[A0], A0);
if I = 1
then DSubN(@R[A0], B, @B[A0], A0)
else
if I = -1
then
begin
DSubN(@R[A0], @B[A0], B, A0);
J := J
xor 1;
end else DFill(0, @R[A0], A0);
DMulK(T, R, A0, @R[A0], A0, @T[AC]);
DMulK(R, A, A0, B, A0, @T[AC]);
DMulK(@R[AC], @A[A0], A0, @B[A0], A0, @T[AC]);
if J <> 0
then J := DAddN(T, T, R, AC)
else J := -DSubN(T, R, T, AC);
Inc(J, DAddN(T, T, @R[AC], AC));
Inc(J, DAddN(@R[A0], @R[A0], T, AC));
if J <> 0
then DIncD(J, @R[AC + A0], A0);
end else
begin
A1 := AC - A0;
K := A[A0];
if K <> 0
then Dec(K, DSubN(R, A, @A[A1], A0))
else
begin
I := DCmpN(A, @A[A1], A0);
if I = 1
then DSubN(R, A, @A[A1], A0)
else
if I = -1
then
begin
DSubN(R, @A[A1], A, A0);
J := 1;
end else DFill(0, R, A0);
end;
R[A0] := K;
K := B[A0];
if K <> 0
then Dec(K, DSubN(@R[A1], B, @B[A1], A0))
else
begin
I := DCmpN(B, @B[A1], A0);
if I = 1
then DSubN(@R[A1], B, @B[A1], A0)
else
if I = -1
then
begin
DSubN(@R[A1], @B[A1], B, A0);
J := J
xor 1;
end else DFill(0, @R[A1], A0);
end;
R[AC] := K;
R0 := AC +1;
DMulK(T, R, A1, @R[A1], A1, @T[R0]);
DMulK(R, A, A1, B, A1, @T[R0]);
DMulK(@R[R0], @A[A1], A0, @B[A1], A0, @T[R0]);
if J <> 0
then DAddN(T, R, T, R0)
else DSubN(T, R, T, R0);
R1 := AC -1;
if DAddN(T, T, @R[R0], R1) <> 0
then
begin
I := T[R1] + 1;
T[R1] := I;
if I = 0
then Inc(T[AC]);
end;
if DAddN(@R[A1], @R[A1], T, R0) <> 0
then
DIncD(1, @R[R0 + A1], A0);
end;
goto Finish;
end;
A0 := (AC +1)
shr 1;
if A0 < BC
then
begin
A1 := AC - A0;
R0 := A0 + A0;
R1 := A1 + BC - A0;
if A0 > A1
then
begin // A1 - A0
if A[A1] = 0
then
begin
I := DSubC(R, @A[A0], A, A1);
if I = 1
then R[A1] := TDigit(-1)
else R[A1] := 0;
end else
begin
R[A1] := TDigit( -Integer(A[A1]) - DSubN(R, @A[A0], A, A1) );
I := 1;
end;
end else
I := DSubC(R, @A[A0], A, A0);
if I = 0
then goto Null;
J := BC - A0;
// B1
if A0 > J
then // B0 - B1
begin
DMoveZero(@B[A0], T, J, A0);
J := DSubC(@R[A0], B, T, A0);
end else
J := DSubC(@R[A0], B, @B[A0], A0);
if J = 0
then goto Null;
DMulK(T, R, A0, @R[A0], A0, @T[R0]);
if I = -1
then
begin
I := -J;
if I = -1
then I := -DSubN(@T[A0], @T[A0], R, A0)
else I := 0;
end else
begin
I := J;
if DSubN(@T[A0], @T[A0], @R[A0], A0) = 0
then Inc(I);
if I = 1
then
begin
DSubN(@T[A0], @T[A0], R, A0);
I := 0;
end;
end;
J := DMulK(R, A, A0, B, A0, @T[R0]);
if (J > 0)
and (DAddN(T, T, R, J) <> 0)
then
begin
K := R0 - J;
if K > 0
then I := I + DIncD(1, @T[J], K)
else Inc(I);
end;
J := DMulK(@R[R0], @A[A0], A1, @B[A0], BC - A0, @T[R0]);
if (J > 0)
and (DAddN(T, T, @R[R0], J) <> 0)
then
begin
K := R0 - J;
if K > 0
then I := I + DIncD(1, @T[J], K)
else Inc(I);
end;
I := I + DAddN(@R[A0], @R[A0], T, R0);
if I > 0
then
DIncD(I, @R[R0 + A0], A0);
goto Finish;
end;
// A is twice of B
if A0 = BC
then
begin
Dec(A0);
A1 := AC - A0;
R0 := A0 + A0;
R1 := A0 + A1;
// !!!!
// A1 - A0
R[A0] := A[R0];
I := A1 - A0;
if I = 2
then
R[A0 + 1] := A[R0 + 1];
if DSubN(R, @A[A0], A, A0) <> 0
then
DDecD(1, @R[A0], I);
// B0 - B1
DMove(B, @R[A1], A0);
TDigit(I) := R[A1];
TDigit(K) := TDigit(I) - B[A0];
R[A1] := TDigit(K);
if TDigit(K) > TDigit(I)
then I := -DDecD(1, @R[A1 +1], A0 -1)
else I := 0;
DMulK(T, R, A1, @R[A1], A0, @T[AC]);
if I = -1
then
DSubN(@T[A0], @T[A0], R, A1);
J := DMulK(R, A, A0, B, A0, @T[AC]);
if (J > 0)
and (DAddN(T, T, R, J) <> 0)
then
begin
K := R1 - J;
if K > 0
then I := I + DIncD(1, @T[J], K)
else Inc(I);
end;
J := DMulK(@R[R0], @A[A0], A1, @B[A0], BC - A0, @T[AC]);
if (J > 0)
and (DAddN(T, T, @R[R0], J) <> 0)
then
begin
K := R1 - J;
if K > 0
then I := I + DIncD(1, @T[J], K)
else Inc(I);
end;
I := I + DAddN(@R[A0], @R[A0], T, R1);
if I > 0
then
begin
A0 := A0 + R1;
DIncD(I, @R[A0], AC + BC - A0);
end;
goto Finish;
end;
// 1.) large digits MUL, remarks see below, A is more as twice long as B
if not Odd(Cardinal(AC)
div Cardinal(BC))
then
begin
DMulK(R, A, BC, B, BC, T);
DMove(@R[BC], T, BC);
I := BC;
end else I := 0;
Dec(AC, BC + BC);
repeat
if I > AC
then Break;
DMulK(@R[I], @A[I], BC, B, BC, @T[I]);
Inc(I, BC);
if I > AC
then Break;
DMulK(@T[I - BC], @A[I], BC, B, BC, @T[I + BC]);
Inc(I, BC);
until False;
AC := AC + BC + BC - I;
if AC > 0
then DMulK(@R[I], B, BC, @A[I], AC, @T[I + BC])
else Dec(I, BC);
if DAddN(@R[BC], @R[BC], T, I) <> 0
then
if AC > 0
then
DIncD(1, @R[BC + I], AC);
goto Finish;
Null:
// A0 - A1 = 0 or B1 - B0 = 0,
// special case, we can leave some calculations
DMulK(R, A, A0, B, A0, T);
DMulK(@R[R0], @A[A0], A1, @B[A0], BC - A0, T);
I := DAddN(T, R, @R[R0], R1);
// T01 := [R01] + R23
if DAddN(@R[R1 + A0], @R[R1 + A0], @R[R1], R0 - R1) <> 0
then // [R12] := [R12] + [R01]
I := I + DIncD(1, @R[R0 + A0], R1 - A0);
I := I + DAddN(@R[A0], @R[A0], T, R1);
// [R12] := [R12] + T01
if I > 0
then
DIncD(I, @R[R1 + A0], A0);
Finish:
// correct result, speedup outer Karatsuba
// if R[Result -1] = 0 then Dec(Result);
Result := DNorm(0, R, Result);
end;
{ Remark on point 1.) Karatsuba/TOOM large Digit Mul,
here are A.Count > 2 * B.Count,
Instead to split symmetric and call DMulK() recursive again we use a
large digit MUL of B.Count Digitsize. The Results of each large Digit MUL
are stored interleaved into R and T
A0A1 A2A3 A4A5 A6A7 * B0B1
C0C1 C2C3 <- A0A1 * B0B1
D0D1 D2D3 <- A2A3 * B0B1
E0E1 E2E3 <- A4A5 * B0B1
F0F1 F2F3 <- A6A7 * B0B1
stored interleaved into R and T and added with only ONE final addition
R C0C1 D0D1 D2D3 F0F1 F2F3
+ T C2C3 E0E1 E2E3
= R R0R1 R2R3 R4R5 R6R7 R8R9
In above case we avoid many useloss moves/coping and reduce it to maximal
one move of C2C3 from R to T.
Same method is applied for asymmetric Toom-Cook Multiplication.
}
var
DSqrS:
function(R,A: Pointer; C: Cardinal): Integer = _DSqrS;
function DSqrK(R,A: PDigits; C: Integer; T: PDigits): Integer;
// karatsuba Square, don't call directly use DSqrN()
// R[] = A[]², A are C digits long, R are 2C digits long,
// T temporary must be 2C digits long
var
I,D,J: Integer;
begin
NCalcCheck;
I := C + C;
if A[C -1] = 0
then
begin
C := DNorm(0, A, C -1);
D := C + C;
DFill(0, @R[D], I - D);
Result := D;
if D = 0
then Exit;
end else Result := I;
if C < Sqr_KBE
then
begin
Result := DSqrS(R, A, C);
Exit;
end;
D := C
shr 1;
if Odd(C)
then
begin
I := D +1;
TDigit(J) := A[D];
// code with implicit comparsion in subtraction and postprocessed correction if overflow is detected
if DSubN(R, A, @A[I], D) <> 0
then
if J = 0
then DCplN(R, R, D)
// we must inverse our overflow
else Dec(J);
// code with preprocessing comparsion to detect overflows,
// suppose <50% of our computation produce overflows then above code avoids >50% comparsions
{ if J = 0 then
begin
case DCmpN(A, @A[I], D) of
1: DSubN(R, A, @A[I], D);
-1: DSubN(R, @A[I], A, D);
0: DFill(0, R, D);
end;
end else Dec(J, DSubN(R, A, @A[I], D));}
R[D] := TDigit(J);
Inc(C);
DSqrK(T, R, I, @T[C]);
DSqrK(R, A, I, @T[C]);
DSqrK(@R[C], @A[I], D, @T[C]);
DSubN(T, R, T, C);
if DAddN(T, T, @R[C], C -2) <> 0
then
begin
TDigit(J) := T[C - 2] + 1;
T[C -2] := TDigit(J);
if J = 0
then Inc(T[C -1]);
end;
if DAddN(@R[I], @R[I], T, C) <> 0
then
DIncD(1, @R[C + I], D);
end else
begin
// code with implicit comparsion in subtraction
I := DSubC(R, @A[D], A, D);
if I <> 0
then
begin
if I > 0
then DCplN(R, R, D);
DSqrK(T, R, D, @T[C]);
end;
// code with comparsion
{ I := DCmpN(A, @A[D], D);
if I <> 0 then
begin
if I > 0 then DSubR(A, @A[D], D, R)
else DSubR(@A[D], A, D, R);
DSqrK(T, R, D, @T[C]);
end;}
DSqrK(@R[C], @A[D], D, @T[C]);
J := DSqrK(R, A, D, @T[C]);
if I <> 0
then
begin
// I := DSubR(R, T, C, T) + DAddN(T, @R[C], C) + DAddN(@R[D], T, C);
I := DSubN(T, T, @R[C], C);
if (J > 0)
and (DSubN(T, T, R, J) <> 0)
then
Inc(I, DDecD(1, @T[J], C - J));
Dec(I, DSubN(@R[D], @R[D], T, C));
end else
if J > 0
then I := DAddN(T, R, @R[C], C) + DAddN(@R[D], @R[D], T, C)
else I := DAddN(@R[D], @R[D], @R[C], C);
if I > 0
then DIncD(I, @R[D + C], D);
end;
// correct result, speedup outer Karatsuba
if R[Result -1] = 0
then Dec(Result);
end;