{
Delphi 4 內建 Int64 是最新功能，可是更大的位數呢？前一陣子因為研究密碼
學的緣故必須使用到任意位數的整數運算，因此就花了點時間研究了一下。
說起來這個大數 Library 有點半弔子，因為那時寫得實在很興奮，以為可以自己
可以寫一套非常完整的 Public-Key Security System ，包括 Encryption/
Decryption/Sign/Verify signature 等等的功能，後來見識到 PGP 的威力之後
就發現其實這不過雕蟲小技。可是對一般人來說，這個 Library 算是相當足夠了。
提供了：加、減、乘、除的基本運算、十進位字串及大數間的轉換、十六進位字串
及大數間的轉換、算出a^b mod n、比較大小、shift left/right，等等功能。
注意事項：
1.大數使用 Dynamic Array ，因此變數拷貝的時候僅僅拷貝參考而已，不是本體，
特別注意（請參照 Delphi 4使用手冊）
2.此大數 Library 不提供負數、小數，請自己再外包一層 :Q，因為密碼學似乎不
需要用到負數以及小數。
}
unit BigNum;
interface
uses
SysUtils, Classes, Windows;
type
{ TDigits: a bignum is no more than a dynamic array containing variable
number of cardinals(unsigned 32-bit integers) }
TDigits = array of Cardinal;
EBigError = class(Exception)
end;
var
BigZero: TDigits;
BigOne: TDigits;
BigTwo: TDigits;
procedure _BigNormalize(var a: TDigits);
function BigAdd(a, b: TDigits): TDigits;
function _BigSubAt(var a: TDigits; const b: TDigits; ka, kb: Integer): Integer;
function BigSub(a, b: TDigits): TDigits;
function BigMul(a, b: TDigits): TDigits; overload;
function BigMul(a: TDigits; b: Cardinal): TDigits; overload;
function TrimLeftChar(const S: string; C: Char): string;
function BigToStrHex(a: TDigits): string;
function BigConvertHex(S: string): TDigits; overload;
function BigConvert(S: string): TDigits; overload;
function BigCopy(a: TDigits; paddings: Integer = 0): TDigits;
function BigCompareAt(a, b: TDigits; ka, kb: Integer): Integer;
function BigCompare(a, b: TDigits): Integer; overload;
procedure _BigDivBy2(var a: TDigits);
function BigDivBy2(a: TDigits; Count: Integer = 1): TDigits; overload;
function BigDiv(u, b: TDigits; var Remainder: TDigits): TDigits; overload;
function BigIsOdd(a: TDigits): Boolean;
function BigDiv(a, b: TDigits): TDigits; overload;
function BigMod(a, b: TDigits): TDigits; overload;
function BigExpMod(x, y, n: TDigits): TDigits;
function BigGCD(x, y: TDigits): TDigits;
function BigLCM(x, y: TDigits): TDigits;
function BigMax(a: array of TDigits): TDigits;
function BigMin(a: array of TDigits): TDigits;
procedure BigSwap(var a, b: TDigits);
function BigPower(x: TDigits; y: Integer): TDigits;
function BigToStr(a: TDigits): string;
implementation
procedure BigError(Msg: string);
begin
raise EBigError.Create(Msg);
end;
function Max(a, b: Integer): Integer;
begin
if a > b then
Result := a
else
Result := b;
end;
function NormalLength(a: TDigits): Integer;
var
I: Integer;
begin
I := High(a);
while (I >= 0) and (a[I] = 0) do
Dec(I);
Result := I + 1;
end;
{ _BigNormalize: Cuts off the trailing zero digits. }
procedure _BigNormalize(var a: TDigits);
var
I: Integer;
begin
I := High(a);
while (I >= 0) and (a[I] = 0) do
Dec(I);
SetLength(a, I+1);
end;
function BigNormalize(a: TDigits): TDigits;
begin
Result := BigCopy(a);
_BigNormalize(a);
end;
{ BigAdd: Performs addition between 'a' and 'b' and returns the result bignum }
function BigAdd(a, b: TDigits): TDigits;
var
I: Integer;
Sum: Int64;
Carry: Integer;
begin
SetLength(Result, Max(Length(a), Length(b)));
Carry := 0;
for I:=0 to High(Result) do
begin
Sum := Carry;
if I < Length(a) then
Inc(Sum, Int64(a[I]));
if I < Length(b) then
Inc(Sum, Int64(b[I]));
Result[I] := Int64Rec(Sum).Lo;
Carry := Int64Rec(Sum).Hi;
end;
_BigNormalize(Result);
end;
{ _BigSubAt: Performs subtraction in 'a' and 'b' with 'a' and 'b' firstly
shifted by 'ka' and 'kb'.
Returns:
0 if a = b
1 if a > b
-1 if a < b.
If 'b' is less than
'a', then the result should be represented in 2^32's complement
}
function _BigSubAt(var a: TDigits; const b: TDigits; ka, kb: Integer): Integer;
var
I: Integer;
Sum: Int64;
Carry: Integer;
EqualFlag: Boolean;
begin
EqualFlag := True;
Carry := 0;
for I:=0 to High(a) do
begin
Sum := - Carry;
if (I-ka <= High(a)) and (I-ka >= 0) then
Inc(Sum, Int64(a[I-ka]));
if (I-kb <= High(b)) and (I-kb >= 0) then
Dec(Sum, Int64(b[I-kb]));
if Sum <> 0 then
EqualFlag := False;
if Sum < 0 then
begin
a[I] := Sum + (Int64(1) shl 32);
Carry := 1;
end
else
begin
Carry := 0;
a[I] := Sum;
end;
end;
Result := 1;
if Carry <> 0 then
Result := -1;
if EqualFlag then
Result := 0;
end;
{ BigSub: Subtracts 'b' from 'a' and returns the result. If 'b' is less than
'a', then the result should be represented in 2^32's complement }
function BigSub(a, b: TDigits): TDigits;
begin
Result := BigCopy(a);
_BigSubAt(Result, b, 0, 0);
_BigNormalize(Result);
end;
{ BigMul: Performs multiplication between bignums 'a' and 'b' and returns the
result bignum}
function BigMul(a, b: TDigits): TDigits; overload;
var
I, J: Integer;
Temp: array of Int64;
Sum: Int64;
Carry: Cardinal;
begin
SetLength(Temp, Length(a) + Length(b) + 1);
for I:=0 to High(Temp) do
Temp[I] := 0;
for I:=0 to High(a) do
for J:=0 to High(b) do
Inc(Temp[I+J], Int64(a[I]) * Int64(b[J]));
SetLength(Result, Length(Temp));
Carry := 0;
for I:=0 to High(Temp) do
begin
Sum := Temp[I] + Carry;
Result[I] := Int64Rec(Sum).Lo;
Carry := Int64Rec(Sum).Hi;
end;
_BigNormalize(Result);
end;
function BigMul(a: TDigits; b: Cardinal): TDigits; overload;
var
k: TDigits;
begin
SetLength(k, 1);
k[0] := b;
Result := BigMul(a, k);
end;
const
HexStr = '0123456789ABCDEF';
function TrimLeftChar(const S: string; C: Char): string;
var
I, L: Integer;
begin
L := Length(S);
I := 1;
while (I <= L) and (S[I] <= C) do Inc(I);
Result := Copy(S, I, Maxint);
end;
{ BigToStrHex: Converts a bignum into a hexidecimal string }
function BigToStrHex(a: TDigits): string;
var
I: Integer;
function NibbleToStr(Nib: Cardinal): string;
var
I: Integer;
K: Cardinal;
begin
SetLength(Result, 8);
for I:=8 downto 1 do
begin
K := Nib and $F;
Result[I] := HexStr[K+1];
Nib := Nib shr 4;
end;
end;
begin
if Length(a) = 0 then
begin
Result := '0';
Exit;
end;
Result := TrimLeftChar(NibbleToStr(a[High(a)]), '0');
for I:=High(a)-1 downto 0 do
Result := Result + NibbleToStr(a[I]);
if Result = '' then
Result := '0';
end;
{ BigToStrHex: Converts a bignum into a hexidecimal string }
function BigDebug(a: TDigits): string;
var
I: Integer;
function NibbleToStr(Nib: Cardinal): string;
var
I: Integer;
K: Cardinal;
begin
SetLength(Result, 8);
for I:=8 downto 1 do
begin
K := Nib and $F;
Result[I] := HexStr[K+1];
Nib := Nib shr 4;
end;
end;
begin
Result := '';
for I:=High(a) downto 0 do
Result := Result + NibbleToStr(a[I]) + ' ';
end;
{ BigConvertHex: Converts a hexidecimal string into a bignum }
function BigConvertHex(S: string): TDigits; overload;
var
I, J: Integer;
function StrToNibble(P: string): Cardinal;
var
I: Integer;
K: Cardinal;
Base: Cardinal;
begin
Base := 1;
Result := 0;
for I:=8 downto 1 do
begin
K := Pos(P[I], HexStr);
if K = 0 then
BigError('BigConvertHex: Invalid character: ' + P[I]);
Result := Result + (K-1) * Base;
Base := Base shl 4;
end;
end;
begin
SetLength(Result, (Length(S) + 7) div 8);
S := UpperCase(Trim(S));
I := Length(S);
J := 0;
while True do
begin
if I >= 8 then
Result[J] := StrToNibble(Copy(S, I-8+1, 8))
else
begin
Result[J] := StrToNibble(StringOfChar('0', 8-I) + Copy(S, 1, I));
Exit;
end;
Dec(I, 8);
Inc(J);
end;
end;
{ BigConvert: Converts a decimal string into a bignum }
function BigConvert(S: string): TDigits; overload;
var
bcd: array of Byte;
I: Integer;
L: Integer;
Base: Cardinal;
K: Integer;
begin
S := Trim(S);
L := Length(S);
SetLength(bcd, L+1);
SetLength(Result, 0);
SetLength(Result, L div 4 + 1);
for I:=1 to L do
bcd[L-I] := Ord(S[I]) - Ord('0');
Base := 1;
K := 0;
while True do
begin
if bcd[0] and 1 <> 0 then
Result[K] := Result[K] or Base;
for I:=0 to L-1 do
begin
bcd[I] := bcd[I] shr 1;
if bcd[I+1] and 1 <> 0 then
Inc(bcd[I], 5);
end;
Base := Base shl 1;
if Base = 0 then
begin
Base := 1;
Inc(K);
if K > High(Result) then
Break;
end;
end;
_BigNormalize(Result);
end;
function BigCopy(a: TDigits; paddings: Integer = 0): TDigits;
var
I: Integer;
begin
SetLength(Result, Length(a) + paddings);
for I:=0 to High(a) do
Result[I] := a[I];
end;
{ BigCompare: Compares the values in 'a' and 'b' with their digits shifted
by 'ka' and 'kb' first. Assuming 'A' equals to 'a' shifting 'ka' and 'B'
equals to 'b' shifting 'kb' then this function
returns 1: when A > B
-1: when A < B
0: when A = B }
function BigCompareAt(a, b: TDigits; ka, kb: Integer): Integer;
var
I, L: Integer;
da, db: Cardinal;
begin
L := Max(High(a)+ka, High(b)+kb);
for I:=L downto 0 do
begin
if (I-ka <= High(a)) and (I-ka >= 0) then
da := a[I-ka]
else
da := 0;
if (I-kb <= High(b)) and (I-kb >= 0) then
db := b[I-kb]
else
db := 0;
if da > db then
begin
Result := 1;
Exit;
end else if da < db then
begin
Result := -1;
Exit;
end;
end;
Result := 0;
end;
{ BigCompare: Compares a bignum 'a' with another bignum 'b'.
Returns 1: when a > b
-1: when a < b
0: when a = b }
function BigCompare(a, b: TDigits): Integer; overload;
begin
Result := BigCompareAt(a, b, 0, 0);
end;
{ BigCompare: Compares a bignum 'a' with a cardinal 'b'.
Returns 1: when a > b
-1: when a < b
0: when a = b }
function BigCompare(a: TDigits; b: Cardinal): Integer; overload;
var
I: Integer;
begin
if Length(a) = 0 then
begin
if b = 0 then
Result := 0
else
Result := -1;
Exit;
end;
Result := 1;
for I:=1 to High(a) do
if a[I] <> 0 then
Exit;
if a[0] > b then
Result := 1
else if a[0] < b then
Result := -1
else
Result := 0;
end;
{ _BigDivBy2: Divides a bignum by 2. Using 'shr' and is fast. }
procedure _BigDivBy2(var a: TDigits);
var
I: Integer;
begin
for I:=0 to High(a)-1 do
a[I] := (a[I] shr 1) or ((a[I+1] and 1) shl 31);
I := High(a);
a[I] := a[I] shr 1;
end;
function BigDivBy2(a: TDigits; Count: Integer = 1): TDigits; overload;
var
I: Integer;
begin
Result := BigCopy(a);
for I:=1 to Count do
_BigDivBy2(Result);
_BigNormalize(Result);
end;
function BigDiv1(k, b: TDigits; var m: TDigits): TDigits; overload;
var
D: Integer;
QLower, QUpper, Q: Int64;
P: TDigits;
a: TDigits;
Sum: Int64;
Leader: Cardinal;
begin
if BigCompare(b, 0) = 0 then
BigError('BigDiv: Division by zero');
{ Before performing operation on 'k', I make a copy of it and store it in 'a',
this is because dynamic arrays are shared among variables refer to it while
we don't want the caller to expect his variables is modified unsolicited }
a := BigCopy(k);
_BigNormalize(a);
_BigNormalize(b);
{ 'D' is the alignment position on 'a' in respect of 'b'. }
D := High(a) - High(b) + 1;
if D < 0 then
D := 0;
{ Set the length of the 'Result' bignum. The maximum length of the 'Result'
bignum should not exceed High(a) - High(b) + 1. }
SetLength(Result, D);
{ 'Sum' is the the sum of the leading two nibbles in 'a' that is validated
each time when 'a' is modified. }
Sum := 0;
{ 'Leader' is the leading nibble in 'b' and is used to determine the upper
bound of the quotient. }
Leader := b[High(b)];
{ Now begins the main loop processing each digit in 'a' }
while True do
begin
Dec(D);
{ If D < 0, the division must have been over and the remaining digits in
'a' should be the modulo }
if D < 0 then
begin
m := a;
_BigNormalize(m);
_BigNormalize(Result);
Exit;
end;
{ Shifts the digits in 'Sum' }
Int64Rec(Sum).Hi := Int64Rec(Sum).Lo;
Int64Rec(Sum).Lo := a[D+High(b)];
{ 'QUpper' and 'QLower' determines the upperbound and lowerbound of the
real quotient 'Q'. QUpper assumes the remaining bignum as 'Sum'
"appended" by a series zeros. QLower assumes the remaining bignum as
'Sum+1' "appended" by a series zeros. Note that QLower does not include
'Sum div (Leader + 1)'. Because the real quotient is discrete, we simply
assign the lowerbound of the quotient to 'Sum div (Leader + 1) + 1'.
This also helps when we're performing binary search in the following
lines. }
QUpper := Sum div Leader;
QLower := Sum div (Leader + 1) + 1;
{ By obtaining 'QUpper' and 'QLower', we are then able to use binary
search to determine the exact position of 'Q'. }
while True do
begin
Q := (QUpper + QLower) shr 1;
{ If QUpper and QLower differ only by one, we should stop immediately.
This is because the exact value of Q may not be obtained due to the
remaining digits in 'a'. Since Q is always in the valid range of the
real quotient, there should be no reason that Q is other than the real
quotient when the error gap is only 1. }
if QUpper - QLower <= 1 then
begin
Q := QUpper;
repeat
P := BigMul(b, Q);
if BigCompareAt(a, P, 0, D) >= 0 then
Break;
Dec(Q);
until False;
{ Store the calculated real quotient into Result }
Result[D] := Q;
{ Since the real quotient is obtained in Q, we now could subtract
P (P = b * Q) from a and proceed with next round. }
_BigSubAt(a, P, 0, D);
Break;
end;
P := BigMul(b, Q);
{ Binary search scheme }
case BigCompareAt(a, P, 0, D) of
0:
begin
QUpper := Q;
QLower := Q;
end;
1: QLower := Q;
-1: QUpper := Q;
end;
end;
{ Assigns the new value of 'Sum'. }
Sum := a[D+High(b)];
end;
end;
function GetNibble(a: TDigits; I: Integer): Cardinal;
begin
if (I >= 0) and (I < Length(a)) then
Result := a[I]
else
Result := 0;
end;
procedure _BigShr(var a: TDigits; Count: Integer); forward;
procedure _BigShl(var a: TDigits; Count: Integer); forward;
procedure _BigShr(var a: TDigits; Count: Integer);
var
I: Integer;
Seg: Integer;
M, N: Integer;
begin
if Count < 0 then
begin
_BigShl(a, -Count);
Exit;
end;
Seg := Count div 32;
N := Count mod 32;
M := 32 - N;
if N = 0 then
for I:=0 to High(a) do
a[I] := GetNibble(a, I+Seg)
else
for I:=0 to High(a) do
a[I] := (GetNibble(a, I+Seg) shr N) or (GetNibble(a, I+Seg+1) shl M);
_BigNormalize(a);
end;
procedure _BigShl(var a: TDigits; Count: Integer);
var
I: Integer;
Seg: Integer;
M, N: Integer;
begin
if Count < 0 then
begin
_BigShr(a, -Count);
Exit;
end;
Seg := Count div 32;
N := Count mod 32;
M := 32 - N;
SetLength(a, Length(a) + Seg + 1);
if N = 0 then
for I:=High(a) downto 0 do
a[I] := GetNibble(a, I-Seg) shl N
else
for I:=High(a) downto 0 do
a[I] := (GetNibble(a, I-Seg) shl N) or (GetNibble(a, I-Seg-1) shr M);
_BigNormalize(a);
end;
procedure _BigAlign(var a: TDigits; var Count: Integer);
var
I: Integer;
N, M: Integer;
Leader, Mask: Cardinal;
begin
_BigNormalize(a);
N := 0;
Leader := a[High(a)];
Mask := Cardinal(1) shl 31;
while Leader and Mask = 0 do
begin
Inc(N);
Mask := Mask shr 1;
end;
M := 32 - N;
for I:=High(a) downto 1 do
a[I] := (a[I] shl N) or (a[I-1] shr M);
a[0] := a[0] shl N;
Count := N;
end;
function BigDiv(u, b: TDigits; var Remainder: TDigits): TDigits; overload;
var
I: Integer;
K: Integer;
Mask: Cardinal;
Count: Integer;
D: Integer;
a: TDigits;
begin
if BigCompare(b, 0) = 0 then
BigError('BigDiv: Division by zero');
a := BigCopy(u, Length(b));
Mask := Cardinal(1) shl 31;
SetLength(Result, 0);
SetLength(Result, High(a) - High(b));
K := High(Result);
Count := 32*Length(u);
D := Length(u);
for I:=0 to Count-1 do
begin
_BigShl(a, 1);
if BigCompareAt(a, b, 0, D) >= 0 then
begin
_BigSubAt(a, b, 0, D);
Result[K] := Result[K] or Mask;
end;
Mask := Mask shr 1;
if Mask = 0 then
begin
Mask := Cardinal(1) shl 31;
Dec(K);
end;
end;
_BigShr(a, Count);
Remainder := a;
_BigNormalize(Result);
end;
function BigIsDivisible(a: TDigits; b: Cardinal): Boolean;
var
Sum: Int64;
I: Integer;
begin
Sum := 0;
for I:=High(a) downto 0 do
begin
Int64Rec(Sum).Hi := Int64Rec(Sum).Lo;
Int64Rec(Sum).Lo := a[I];
Sum := Sum mod b;
end;
Result := Sum = 0;
end;
function BigDiv(a: TDigits; b: Cardinal; var Remainder: Cardinal): TDigits; overload
var
Sum: Int64;
I: Integer;
begin
Sum := 0;
SetLength(Result, Length(a));
for I:=High(a) downto 0 do
begin
Int64Rec(Sum).Hi := Int64Rec(Sum).Lo;
Int64Rec(Sum).Lo := a[I];
Result[I] := Sum div b;
Sum := Sum mod b;
end;
Remainder := Sum;
end;
function BigMod(a: TDigits; b: Cardinal): Cardinal; overload;
begin
BigDiv(a, b, Result);
end;
procedure BigInc(var a: TDigits; b: Cardinal);
var
I: Integer;
begin
I := 0;
while True do
begin
if High(a) < I then
begin
SetLength(a, I+1);
a[I] := 0;
end;
if a[I] = $FFFFFFFF then
a[I] := 0
else
begin
Inc(a[I]);
Exit;
end;
end;
end;
procedure BigDec(var a: TDigits);
var
I: Integer;
begin
I := 0;
while True do
begin
if High(a) < I then
BigError('BigDec: Decrement underflow');
if a[I] = 0 then
a[I] := $FFFFFFFF
else
begin
Dec(a[I]);
_BigNormalize(a);
Exit;
end;
end;
end;
function BigIsOdd(a: TDigits): Boolean;
begin
if Length(a) = 0 then
Result := False
else
Result := a[0] and 1 <> 0;
end;
function BigIsEven(a: TDigits): Boolean;
begin
if Length(a) = 0 then
Result := True
else
Result := a[0] and 1 = 0;
end;
function BigDiv(a, b: TDigits): TDigits; overload;
var
m: TDigits;
begin
Result := BigDiv(a, b, m);
end;
function BigMod(a, b: TDigits): TDigits; overload
begin
BigDiv(a, b, Result);
end;
function BigExpMod(x, y, n: TDigits): TDigits;
var
s, t, u: TDigits;
begin
s := BigOne;
t := x;
u := BigCopy(y);
while BigCompare(u, 0) > 0 do
begin
if BigIsOdd(u) then
s := BigMod(BigMul(s, t), n);
_BigDivBy2(u);
t := BigMod(BigMul(t, t), n);
end;
Result := s;
end;
function BigGCD(x, y: TDigits): TDigits; overload;
begin
if (BigCompare(x, 0) = 0) and (BigCompare(y, 0) = 0) then
BigError('BigGCD: Cannot compute the greatest common divisor between two zero bignum.');
while True do
begin
if BigCompare(x, 0) = 0 then
begin
Result := y;
Exit;
end;
if BigCompare(y, 0) = 0 then
begin
Result := x;
Exit;
end;
if BigCompare(x, y) > 0 then
x := BigMod(x, y)
else
y := BigMod(y, x);
end;
end;
function BigGCD(x: array of TDigits): TDigits; overload;
var
g: TDigits;
I: Integer;
begin
if Length(x) = 0 then
BigError('BigGCD: Cannot compute the greatest common divisor when no arguments are given.');
g := x[0];
for I:=1 to High(x) do
begin
g := BigGCD(g, x[I]);
if BigCompare(g, 1) = 0 then
Break;
end;
Result := g;
end;
function BigLCM(x, y: TDigits): TDigits;
begin
Result := BigDiv(BigMul(x, y), BigGCD(x, y));
end;
{ BigMax: Computes the maximum value of a given list. Note that this function
DOES NOT returns copied value of the maximum. }
function BigMax(a: array of TDigits): TDigits;
var
I: Integer;
begin
if Length(a) = 0 then
BigError('BigMax: Cannot compute maximum value given no argument at all.');
Result := a[0];
for I:=1 to High(a) do
if BigCompare(Result, a[I]) > 0 then
Result := a[I];
end;
{ BigMin: Computes the minimum value of a given list. Note that this function
DOES NOT returns copied value of the minimum. }
function BigMin(a: array of TDigits): TDigits;
var
I: Integer;
begin
if Length(a) = 0 then
BigError('BigMin: Cannot compute minimum value given no argument at all.');
Result := a[0];
for I:=1 to High(a) do
if BigCompare(Result, a[I]) < 0 then
Result := a[I];
end;
{ BigSwap: Exchanges values between 'a' and 'b'. }
procedure BigSwap(var a, b: TDigits);
var
t: TDigits;
begin
t := a;
a := b;
b := t;
end;
function BigPower(x: TDigits; y: Integer): TDigits;
var
I: Integer;
begin
Result := BigOne;
for I:=1 to y do
Result := BigMul(Result, x);
end;
{ BigToStr: Converts a bignum into a decimal string }
function BigToStr(a: TDigits): string;
var
bcd: array of Byte;
K: Integer;
Base: Cardinal;
I: Integer;
Carry: Cardinal;
L: Integer;
begin
if Length(a) = 0 then
begin
Result := '0';
Exit;
end;
SetLength(bcd, Length(a) * 10);
Base := Cardinal(1) shl 31;
K := High(a);
while True do
begin
if a[K] and Base = 0 then
Carry := 0
else
Carry := 1;
for I:=0 to High(bcd) do
begin
bcd[I] := (bcd[I] shl 1) or Carry;
if bcd[I] < 10 then
Carry := 0
else
begin
Carry := 1;
Dec(bcd[I], 10);
end;
end;
Base := Base shr 1;
if Base = 0 then
begin
Base := Cardinal(1) shl 31;
Dec(K);
if K < 0 then
Break;
end;
end;
I := High(bcd);
while (I >= 0) and (bcd[I] = 0) do
Dec(I);
L := I+1;
SetLength(Result, L);
for I:=0 to L-1 do
Result[L-I] := Char(bcd[I] + Ord('0'));
end;
{ ExtBinEuclid: Performs extended Euclidean algorithm (from page 247) }
procedure _ExtBinEuclid(var u, v, u1, u2, u3: TDigits);
var
k: Integer;
t1, t2, t3: TDigits;
begin
if BigCompare(u, v) < 0 then
BigSwap(u, v);
k := 0;
while BigIsEven(u) and BigIsEven(v) do
begin
_BigShr(u, 1);
_BigShr(v, 1);
Inc(K);
end;
u1 := BigCopy(BigOne);
u2 := BigCopy(BigZero);
u3 := BigCopy(u);
t1 := BigCopy(v);
t2 := BigSub(u, BigOne);
t3 := BigCopy(v);
repeat
repeat
if BigIsEven(u3) then
begin
if BigIsOdd(u1) or BigIsOdd(u2) then
begin
u1 := BigAdd(u1, v);
u2 := BigAdd(u2, u);
end;
_BigShr(u1, 1);
_BigShr(u2, 1);
_BigShr(u3, 1);
end;
if BigIsEven(t3) or (BigCompare(u3, t3) < 0) then
begin
BigSwap(u1, t1);
BigSwap(u2, t2);
BigSwap(u3, t3);
end;
until BigIsEven(u3);
while (BigCompare(u1, t1) < 0) or (BigCompare(u2, t2) < 0) do
begin
u1 := BigAdd(u1, v);
u2 := BigAdd(u2, u);
end;
_BigSubAt(u1, t1, 0, 0);
_BigSubAt(u2, t2, 0, 0);
_BigSubAt(u3, t3, 0, 0);
until BigCompare(t3, 0) > 0;
_BigShl(u1, k);
_BigShl(u2, k);
_BigShl(u3, k);
end;
function BigExtEuclid(u1, v1: TDigits): TDigits;
var
a, b, gcd, u, v: TDigits;
begin
u := BigCopy(u1);
v := BigCopy(v1);
_ExtBinEuclid(u, v, a, b, gcd);
if BigCompare(gcd, 1) = 0 then
Result := BigSub(u, b)
else
Result := BigZero;
end;
function BigRandom(Nibbles: Integer): TDigits; overload;
var
I: Integer;
begin
SetLength(Result, Nibbles);
for I:=0 to Nibbles-1 do
Result[I] := Cardinal(Random(-1)+1);
while Result[High(Result)] = 0 do
Result[High(Result)] := Random(0);
end;
function BigRandom(p: TDigits): TDigits; overload;
var
L: Integer;
begin
L := NormalLength(p);
if L = 0 then
begin
SetLength(Result, 0);
Exit;
end;
Result := BigRandom(L);
Result[L-1] := Cardinal(Random(p[L-1]));
end;
function RabbinMillerPrimeTest(p: TDigits): Boolean;
var
p_1: TDigits;
a, z, m: TDigits;
b, j: Cardinal;
label
Step4;
begin
Result := True;
b := 0;
p_1 := BigSub(p, BigOne);
m := BigCopy(p_1);
while BigIsEven(m) do
begin
Inc(b);
_BigShr(m, 1);
end;
Write('@');
a := BigRandom(1);
j := 0;
z := BigExpMod(a, m, p);
Write('!');
if (BigCompare(z, 1) = 0) or (BigCompare(z, p_1) = 0) then
begin
Result := True;
Exit;
end;
Step4:
Write('#');
if (j > 0) and (BigCompare(z, 1) = 0) then
begin
Result := False;
Exit;
end;
Inc(j);
if (j < b) and (BigCompare(z, p_1) <> 0) then
begin
z := BigExpMod(z, BigTwo, p);
goto Step4;
end;
if BigCompare(z, p_1) = 0 then
begin
Result := True;
Exit;
end;
if (j = b) and (BigCompare(z, p_1) <> 0) then
Result := False;
end;
var
Primes : array[1..1026] of Integer =
(3, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127,
131, 137, 139, 149, 151, 157, 163, 167, 173, 179, 181, 191, 193, 197, 199, 211, 223, 227, 229, 233, 239, 241, 251,
257, 263, 269, 271, 277, 281, 283, 293, 307, 311, 313, 317, 331, 337, 347, 349, 353, 359, 367, 373, 379, 383, 389,
397, 401, 409, 419, 421, 431, 433, 439, 443, 449, 457, 461, 463, 467, 479, 487, 491, 499, 503, 509, 521, 523, 541,
547, 557, 563, 569, 571, 577, 587, 593, 599, 601, 607, 613, 617, 619, 631, 641, 643, 647, 653, 659, 661, 673, 677,
683, 691, 701, 709, 719, 727, 733, 739, 743, 751, 757, 761, 769, 773, 787, 797, 809, 811, 821, 823, 827, 829, 839,
853, 857, 859, 863, 877, 881, 883, 887, 907, 911, 919, 929, 937, 941, 947, 953, 967, 971, 977, 983, 991, 997, 1009,
1013, 1019, 1021, 1031, 1033, 1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097, 1103, 1109, 1117, 1123,
1129, 1151, 1153, 1163, 1171, 1181, 1187, 1193, 1201, 1213, 1217, 1223, 1229, 1231, 1237, 1249, 1259, 1277, 1279,
1283, 1289, 1291, 1297, 1301, 1303, 1307, 1319, 1321, 1327, 1361, 1367, 1373, 1381, 1399, 1409, 1423, 1427, 1429,
1433, 1439, 1447, 1451, 1453, 1459, 1471, 1481, 1483, 1487, 1489, 1493, 1499, 1511, 1523, 1531, 1543, 1549, 1553,
1559, 1567, 1571, 1579, 1583, 1597, 1601, 1607, 1609, 1613, 1619, 1621, 1627, 1637, 1657, 1663, 1667, 1669, 1693,
1697, 1699, 1709, 1721, 1723, 1733, 1741, 1747, 1753, 1759, 1777, 1783, 1787, 1789, 1801, 1811, 1823, 1831, 1847,
1861, 1867, 1871, 1873, 1877, 1879, 1889, 1901, 1907, 1913, 1931, 1933, 1949, 1951, 1973, 1979, 1987, 1993, 1997,
1999, 2003, 2011, 2017, 2027, 2029, 2039, 2053, 2063, 2069, 2081, 2083, 2087, 2089, 2099, 2111, 2113, 2129, 2131,
2137, 2141, 2143, 2153, 2161, 2179, 2203, 2207, 2213, 2221, 2237, 2239, 2243, 2251, 2267, 2269, 2273, 2281, 2287,
2293, 2297, 2309, 2311, 2333, 2339, 2341, 2347, 2351, 2357, 2371, 2377, 2381, 2383, 2389, 2393, 2399, 2411, 2417,
2423, 2437, 2441, 2447, 2459, 2467, 2473, 2477, 2503, 2521, 2531, 2539, 2543, 2549, 2551, 2557, 2579, 2591, 2593,
2609, 2617, 2621, 2633, 2647, 2657, 2659, 2663, 2671, 2677, 2683, 2687, 2689, 2693, 2699, 2707, 2711, 2713, 2719,
2729, 2731, 2741, 2749, 2753, 2767, 2777, 2789, 2791, 2797, 2801, 2803, 2819, 2833, 2837, 2843, 2851, 2857, 2861,
2879, 2887, 2897, 2903, 2909, 2917, 2927, 2939, 2953, 2957, 2963, 2969, 2971, 2999, 3001, 3011, 3019, 3023, 3037,
3041, 3049, 3061, 3067, 3079, 3083, 3089, 3109, 3119, 3121, 3137, 3163, 3167, 3169, 3181, 3187, 3191, 3203, 3209,
3217, 3221, 3229, 3251, 3253, 3257, 3259, 3271, 3299, 3301, 3307, 3313, 3319, 3323, 3329, 3331, 3343, 3347, 3359,
3361, 3371, 3373, 3389, 3391, 3407, 3413, 3433, 3449, 3457, 3461, 3463, 3467, 3469, 3491, 3499, 3511, 3517, 3527,
3529, 3533, 3539, 3541, 3547, 3557, 3559, 3571, 3581, 3583, 3593, 3607, 3613, 3617, 3623, 3631, 3637, 3643, 3659,
3671, 3673, 3677, 3691, 3697, 3701, 3709, 3719, 3727, 3733, 3739, 3761, 3767, 3769, 3779, 3793, 3797, 3803, 3821,
3823, 3833, 3847, 3851, 3853, 3863, 3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931, 3943, 3947, 3967,
3989, 4001, 4003, 4007, 4013, 4019, 4021, 4027, 4049, 4051, 4057, 4073, 4079, 4091, 4093, 4099, 4111, 4127, 4129,
4133, 4139, 4153, 4157, 4159, 4177, 4201, 4211, 4217, 4219, 4229, 4231, 4241, 4243, 4253, 4259, 4261, 4271, 4273,
4283, 4289, 4297, 4327, 4337, 4339, 4349, 4357, 4363, 4373, 4391, 4397, 4409, 4421, 4423, 4441, 4447, 4451, 4457,
4463, 4481, 4483, 4493, 4507, 4513, 4517, 4519, 4523, 4547, 4549, 4561, 4567, 4583, 4591, 4597, 4603, 4621, 4637,
4639, 4643, 4649, 4651, 4657, 4663, 4673, 4679, 4691, 4703, 4721, 4723, 4729, 4733, 4751, 4759, 4783, 4787, 4789,
4793, 4799, 4801, 4813, 4817, 4831, 4861, 4871, 4877, 4889, 4903, 4909, 4919, 4931, 4933, 4937, 4943, 4951, 4957,
4967, 4969, 4973, 4987, 4993, 4999, 5003, 5009, 5011, 5021, 5023, 5039, 5051, 5059, 5077, 5081, 5087, 5099, 5101,
5107, 5113, 5119, 5147, 5153, 5167, 5171, 5179, 5189, 5197, 5209, 5227, 5231, 5233, 5237, 5261, 5273, 5279, 5281,
5297, 5303, 5309, 5323, 5333, 5347, 5351, 5381, 5387, 5393, 5399, 5407, 5413, 5417, 5419, 5431, 5437, 5441, 5443,
5449, 5471, 5477, 5479, 5483, 5501, 5503, 5507, 5519, 5521, 5527, 5531, 5557, 5563, 5569, 5573, 5581, 5591, 5623,
5639, 5641, 5647, 5651, 5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717, 5737, 5741, 5743, 5749, 5779,
5783, 5791, 5801, 5807, 5813, 5821, 5827, 5839, 5843, 5849, 5851, 5857, 5861, 5867, 5869, 5879, 5881, 5897, 5903,
5923, 5927, 5939, 5953, 5981, 5987, 6007, 6011, 6029, 6037, 6043, 6047, 6053, 6067, 6073, 6079, 6089, 6091, 6101,
6113, 6121, 6131, 6133, 6143, 6151, 6163, 6173, 6197, 6199, 6203, 6211, 6217, 6221, 6229, 6247, 6257, 6263, 6269,
6271, 6277, 6287, 6299, 6301, 6311, 6317, 6323, 6329, 6337, 6343, 6353, 6359, 6361, 6367, 6373, 6379, 6389, 6397,
6421, 6427, 6449, 6451, 6469, 6473, 6481, 6491, 6521, 6529, 6547, 6551, 6553, 6563, 6569, 6571, 6577, 6581, 6599,
6607, 6619, 6637, 6653, 6659, 6661, 6673, 6679, 6689, 6691, 6701, 6703, 6709, 6719, 6733, 6737, 6761, 6763, 6779,
6781, 6791, 6793, 6803, 6823, 6827, 6829, 6833, 6841, 6857, 6863, 6869, 6871, 6883, 6899, 6907, 6911, 6917, 6947,
6949, 6959, 6961, 6967, 6971, 6977, 6983, 6991, 6997, 7001, 7013, 7019, 7027, 7039, 7043, 7057, 7069, 7079, 7103,
7109, 7121, 7127, 7129, 7151, 7159, 7177, 7187, 7193, 7207, 7211, 7213, 7219, 7229, 7237, 7243, 7247, 7253, 7283,
7297, 7307, 7309, 7321, 7331, 7333, 7349, 7351, 7369, 7393, 7411, 7417, 7433, 7451, 7457, 7459, 7477, 7481, 7487,
7489, 7499, 7507, 7517, 7523, 7529, 7537, 7541, 7547, 7549, 7559, 7561, 7573, 7577, 7583, 7589, 7591, 7603, 7607,
7621, 7639, 7643, 7649, 7669, 7673, 7681, 7687, 7691, 7699, 7703, 7717, 7723, 7727, 7741, 7753, 7757, 7759, 7789,
7793, 7817, 7823, 7829, 7841, 7853, 7867, 7873, 7877, 7879, 7883, 7901, 7907, 7919, 7927, 7933, 7937, 7949, 7951,
7963, 7993, 8009, 8011, 8017, 8039, 8053, 8059, 8069, 8081, 8087, 8089, 8093, 8101, 8111, 8117, 8123, 8147, 8161,
8167, 8171, 8179, 8191);
function BigPrimeTest(p: TDigits; nTest: Integer = 5): Boolean;
var
I: Integer;
begin
Result := False;
for I:=Low(Primes) to High(Primes) do
if BigIsDivisible(p, Primes[I]) then
Exit;
Write('-');
for I:=1 to nTest do
if not RabbinMillerPrimeTest(p) then
Exit;
Result := True;
end;
function BigGeneratePrime(Nibbles: Integer): TDigits;
var
Round: Integer;
begin
Round := 0;
while True do
begin
Result := BigRandom(Nibbles);
Result[Nibbles-1] := Result[Nibbles-1] or Cardinal(1 shl 31);
Result[0] := Result[0] or 1;
Write('*');
if BigPrimeTest(Result, 5) then
begin
WriteLn('(r=', round, ')');
Exit;
end;
Write('.');
Inc(Round);
end;
end;
procedure TestGCD(a, b: string);
begin
WriteLn('gcd(',a , ',', b,')=',BigToStr(BigGCD(BigConvert(a), BigConvert(b))));
end;
procedure TestDiv(sa, sb: string); overload;
var
a, b, r, m: TDigits;
begin
a := BigConvert(sa);
b := BigConvert(sb);
r := BigDiv(a, b, m);
WriteLn('BigDiv(',sa , ',', sb,')=', BigToStr(r), ' ... ', BigToStr(m));
end;
procedure TestDiv(sa: string; b: Cardinal); overload;
var
a, r: TDigits;
m: Cardinal;
begin
a := BigConvert(sa);
r := BigDiv(a, b, m);
WriteLn('BigDiv(',sa , ',', b,')=', BigToStr(r), ' ... ', m);
end;
procedure TestRandom;
var
s: string;
begin
while True do
begin
WriteLn('BigRandom=' + BigToStr(BigRandom(16)));
// ReadLn(s);
if s <> '' then
Exit;
end;
end;
procedure TestPrime;
var
I: Integer;
begin
for I:=1 to 5 do
begin
WriteLn('Prime=' + BigToStr(BigGeneratePrime(4)));
end;
end;
procedure TestShr;
var
a: TDigits;
begin
a := BigConvertHex('1');
while BigCompare(a, 0) > 0 do
begin
_BigShl(a, 1);
WriteLn(BigToStr(a));
ReadLn;
end;
end;
procedure TestExpMod(x, y, n: string);
begin
WriteLn(x, '^', y, ' mod ', n, ' = ', BigToStr(BigExpMod(BigConvert(x), BigConvert(y), BigConvert(n))));
end;
procedure Main;
begin
TestDiv('1', '1');
TestDiv('0', '1');
TestDiv('99999', '1');
TestDiv('50000000000000000000000003', 5);
TestDiv('10000000000000000', '9');
TestDiv('10000000000000000', '4294967296');
TestDiv('4312092353817979174', '444877445811470254');
TestDiv('4312092353817979174', '444877445811470254');
TestDiv('10000000000000000000000000000000000000', '7');
TestExpMod('7', '100', '1234567890');
TestExpMod('9', '100', '1234567890');
TestGCD('44455555874312443214','4432174432919929448122');
TestGCD('1234', '5678');
TestGCD('195265', '266198');
WriteLn('gcd(12,6) = ', BigToStr(BigGCD(BigConvert('12'),BigConvert('6'))));
// TestRandom;
// Randomize;
TestPrime;
end;
procedure Test;
var
a, b, c: Int64;
procedure Unavoid(var p);
begin
end;
begin
a := 10;
Unavoid(a);
b := 5;
Unavoid(b);
c := a div b;
Unavoid(c);
end;
initialization
Test;
BigDebug(BigZero);
BigOne := BigConvert('1');
BigTwo := BigConvert('2');
Main;
end.
--
Chaos is the best description of the constant state of human society.
Therefore, dynamic balance is required for us to to survive when vital fault
occurrs. So in the society, we don't chase for peace and order, but existence
and survival, instead.
--
※ 發信站: 批踢踢實業坊(ptt.twbbs.org)
◆ From: h82.s80.ts30.hi