Manipulation trés rapide de trés grandes nombres entiers + fonctions de base de l'arithmétique modulaire,nombres premiers...
Description
MANIPULATION TRÈS RAPIDE DE TRÈS GRANDS NOMBRES ENTIERS + FONCTIONS DE BASE DE L'ARITHMÉTIQUE MODULAIRE, NOMBRES PREMIERS ...
Description
Dans cette librairie, un type Hugint est défini pour la manipulation de nombres entiers trés grands, avec les fonctions correspondantes (addition, soustraction, multiplication et division). En plus les fonctions mathématiques de base sont implémentées : arithmétique modulaire, génération de nombres premiers trés grands, exponentionation modulaire (Montgomery), résolution d'équations diophantiennes et de congruence, factorisation, ... Le code est très optimisé, la division utilise l'algorithme le plus rapide de Knuth !
Développé par Faraoun Kamel Mohamed. UDL-SBA
Pour toute remarque ou suggestion : Kamel_mh@yahoo.fr
Conclusion
L'application ci-jointe est une démonstration utilisant cette bibliothèque.
L'application en jointe est une démonstration utilisant cette bibliothèque. J'ai attacher un fichier Zip Avec le source qui contien deux démonstrations, s'il n'est pas téléchargable alors mail me at kamel_mh@yahoo.fr.
Merci.
Description
Dans cette librairie, un type Hugint est défini pour la manipulation de nombres entiers trés grands, avec les fonctions correspondantes (addition, soustraction, multiplication et division). En plus les fonctions mathématiques de base sont implémentées : arithmétique modulaire, génération de nombres premiers trés grands, exponentionation modulaire (Montgomery), résolution d'équations diophantiennes et de congruence, factorisation, ... Le code est très optimisé, la division utilise l'algorithme le plus rapide de Knuth !
Développé par Faraoun Kamel Mohamed. UDL-SBA
Pour toute remarque ou suggestion : Kamel_mh@yahoo.fr
Conclusion
L'application ci-jointe est une démonstration utilisant cette bibliothèque.
Source / Exemple :
unit HugIntegers;
{********************************************************************************************}
{ Développer par Faraoun Kamel Mohamed }
{ Université Djilali Liabes -Sidi Bel Abbes - Algérie }
{ kamel_mh@yahoo.fr }
{********************************************************************************************}
interface
uses windows,Sysutils,StdCtrls,forms,Dialogs;
Const MaxLen=4096; {must be in the form 2^^n}
MaxPrimSteps=2629;
Type
{ Define a signed Hug integer structure : -2^^(Maxlen*32)...+2^^(Maxlen*32)}
Hugint=record case integer of
0:(i32:packed array [-1..MaxLen] of Dword);
1:(i16:packed array [-2..MaxLen*2] of word);
2:(i8:packed array [-4..MaxLen*4] of Byte);
end;
HugintArray=Array of Hugint;
PHugintArray=^HugintArray;
MontgomeryStruct=Record
modu,r,v,invr:Hugint;
nbits:Dword;
end;
SafeStruct=Record
Q:Hugint;
k:Dword;
end;
PSafeStruct=^SafeStruct;
Procedure AddHug(const t1,t2,Result:Hugint);
Procedure SubHug(const t1,t2,Result:Hugint);
Procedure MulHug(const t1,t2,Result:Hugint);
Procedure DivideModHug(a,b:Hugint;var Q,res:hugint);
Procedure ModuloHug(a,b:Hugint;var res:hugint);
procedure ShlHug(a:Hugint;count:word);
procedure ShrHug(const a:Hugint;const count:word);
Procedure SqrHug(const a,result:Hugint);
Procedure PowerHug(a:hugint;var Result:hugint ;p:longword);
Procedure ModPowerHug(a,b,modulo:hugint;var Result:hugint);
Function CompareHug(const t1,t2:Hugint):Smallint;
function IsNull(a:Hugint):boolean;
Function IsNeg(a:Hugint):boolean;
Function StrToHug(s:string):Hugint;
Function HugToStr(a:Hugint):String;
Function HugToHex(a:Hugint):String;
Function NumOfBits(a:Hugint):longint;
Function NumOfDecimalDigits(a:Hugint):Word;
Procedure HCopy(a,b:Hugint);
Function HexToHug(s:string):Hugint;
Function Hugtobin(a :Hugint):string;
Function BintoHug (s:string):Hugint;
Function OctaltoHug(S:string):Hugint;
Function HugtoOctal(a:Hugint):string;
Procedure IncHug(Var a:Hugint;p:Dword);
Procedure DecHug(Var a:Hugint;p:Dword);
Procedure AndHug(a,b:Hugint;Var Result:Hugint);
Procedure OrHug(a,b:Hugint;Var Result:Hugint);
Procedure XorHug(a,b:Hugint;Var Result:Hugint);
Procedure NotHug(a:Hugint;Var Result:Hugint);
Function AreHugsRelativlyPrimes(a,b:Hugint):Boolean;
Procedure ExtendedEuclid(a,b:HugInt;Var s,t:HugInt);
Function ResolveDiopantiene(a,b,c:Hugint;var x,y:Hugint):Boolean;
Function ResolveLinearCongruenceEquation(a,b,c:Hugint;Var Sol:PHugintArray;All:boolean):integer;
Procedure InvModHug(const a,modulo:Hugint;Var Result:Hugint);
Procedure GetRandomHugOnBits(Var Result:Hugint;Numbits:Integer);
Procedure GetRandomHugLowerThan(var Result:Hugint;Limit:Hugint);
Procedure InitMontgomeryStruct(Modulo:Hugint);
Procedure MgModPowerHug(a,b,n:Hugint;var Result:Hugint);
Function IsHugPrime(Number:Hugint;Prob:Extended=0.99):Boolean;
Procedure GetHugPrimeOnBits(Numbit:integer;Var Result:Hugint;SkipStep:integer=MaxPrimSteps);
Procedure GetHug4Kp3PrimeOnBits(Numbit:integer;Var Result:Hugint;SkipStep:integer=MaxPrimSteps);
Procedure GetStrongHugPrimeOnBits(Numbit:integer;Var Result:Hugint);
Procedure SqrtHug(Number:Hugint;Var Result:Hugint);
Procedure ModSquareHug(a,modulo:Hugint;Var Result:Hugint);
Procedure ResolveMultiLinearCongSystems(Resid,Modu:HugintArray;WithoutVerif:Boolean;Var Result:HugintArray);
Procedure ResolveQuadraticCongruenceEquation(a,b,c:Hugint;ModuloFactors:HugintArray;var n:Hugint;Var Result:HugintArray);
Procedure MonteCarloFactorizeHug(a:Hugint;var factors:HugintArray;Memo:Tmemo=nil);
Procedure PollardFactorizeHug(a:hugint;var Factors:HugintArray;Memo:Tmemo=nil);
Function BinaryGCDHug(a,b:HugInt):HugInt;
Function EuclidGCDHug(a,b:HugInt):HugInt;
Procedure GetPrime1(Numbit:integer;init:Hugint;Var Result:Hugint;var nbmont,nbitr:integer);
Procedure GetPrime2(Numbit:integer;init:Hugint;Var Result:Hugint;var nbmont,nbitr:integer);
Procedure GetPrime3(Numbit:integer;init:Hugint;Var Result:Hugint;var nbmont,nbitr:integer);
Procedure GetSafePrimeOnBits(Numbit:integer;Var Result:Hugint);
Procedure GetMulSafePrimeOnBits(Numbit:integer;Var Res:Hugint;facts:PSafeStruct=nil);
Procedure GetSafeWithGenerator(Numbit:integer;Var Result,Generator:Hugint);
Var Zero,One,Two,Three,Four,ProdFirstPrimes:Hugint;
MgStruct:MontgomeryStruct;
implementation
uses Unit1,FirstPrimesTable,PrimesDemo;
Var ind:Integer;
{***************************************************************************************}
{ Addition of two signed Hug integeres }
{***************************************************************************************}
Procedure AddHug(const t1,t2,result:Hugint);
var carr,indic,zer:Dword;
sizeA,sizeB:Word;
asm
push ebx
push esi
push edi
mov zer,0
mov bx,Word ptr[edx]
mov sizeB,bx
mov bx,Word ptr [eax] { Compare the size }
mov sizeA,bx
shl sizeA,2
shl sizeB,2
cmp bx,Word ptr [edx] { of the two numbers }
jg @@1 { and set the size }
mov bx,Word ptr [edx] { of the result to }
@@1: mov Word ptr[ecx],bx { the maximum one. }
cmp bx,0
je @@fin2
mov esi,ecx { Save the result adress (ecx).}
mov carr,0
mov indic,0
mov bx,word ptr [eax+2]{ }
cmp bx,Word ptr[edx+2] { Compare the signe of the two numbers}
je @@same { jump if same sign. }
{--------------------------------------------------------}
cmp bx,0 { put the negatif }
je @@t2 { number adr in edx }
xchg eax,edx { }
push ecx
mov cx,sizeA
mov bx,sizeB
mov sizeA,bx
mov sizeB,cx
pop ecx
@@t2: mov esi,ecx { Save ecx. }
inc carr
dec indic
{--------------------------------------------------------}
@@same: mov word ptr[ecx+2],bx {set the sign of the result (will be modified if diffrent signes))}
xor ecx,ecx { Here start the addition with carry. }
mov cx,word ptr [esi]
xor ebx,ebx
@@lp: add ebx,4
mov edi,0
cmp bx,sizeB
jg @@001
mov edi,Dword ptr [edx+ebx]
@@001: xor edi,indic {perform the 2 complement }
add edi,carr
mov carr,0
adc carr,0
cmp bx,sizeA
jg @@002
add edi,Dword ptr[eax+ebx]
adc carr,0
@@002: mov Dword ptr[esi+ebx],edi { esi contain the initial ecx value (@t3)}
or zer,edi
loop @@lp
mov ecx,carr
cmp indic,0
jne @@cmpRes {Jump if the two numbers are with diffrents signs}
add word ptr[esi],cx
mov Dword ptr[esi+ebx+4],ecx
or zer,ecx
jmp @@fin
{--------------------------------------------------------}
@@cmpRes:dec ecx
mov word ptr[esi+2],cx
cmp ecx,0
je @@fin
xor ecx,ecx
mov cx,word ptr[esi]
inc carr
xor ebx,ebx
@@lp3: add ebx,4
mov edi,Dword ptr[esi+ebx]
not edi
add edi,carr
mov Dword ptr[esi+ebx],edi
mov carr,0
adc carr,0
loop @@lp3
{--------------------------------------------------------}
@@fin: cmp zer,0
jne @@lll
mov word ptr[esi],0
jmp @@fin2
@@lll: xor ecx,ecx
mov cx,word ptr [esi]
shl ecx,2
add ecx,esi
@loop: cmp dword ptr[ecx],0
jne @@fin2
dec Word ptr[esi]
sub ecx,4
jmp @loop
@@fin2: mov ecx,esi
pop edi
pop esi
pop ebx
end;
{***************************************************************************************}
{ Substraction of two signed Hug integeres }
{***************************************************************************************}
Procedure SubHug(const t1,t2,Result:Hugint);
asm
push esi
mov esi,edx
cmp eax,edx
je @@zero {compare if t1 and t2 are same }
not word ptr[edx+2] {invesre signe of t2}
call addHug {and perform an addition }
cmp esi,ecx {in the case of t2=Result don't inverse the result sign}
je @@ret
not word ptr[edx+2]
@@ret: pop esi
Ret
{--------------------------------------------------------}
@@zero: push esi { if t1 and t2 are same then return 0}
mov esi,ecx
mov ecx,DWord ptr[ecx]
and ecx,$FFFFFFFF
jz @@fin
@@lp: mov Dword ptr[esi+ecx*4],0
loop @@lp
mov word ptr[esi],0
{--------------------------------------------------------}
@@fin: mov ecx,esi
pop esi
end;
{***************************************************************************************}
{ Faster Substraction of two signed Hug integeres }
{ Work only if t1,t2>=0, and t1>t2 }
{***************************************************************************************}
Procedure FastSub(const t1,t2:Hugint);
var borr:Dword;
asm
push ebx
push esi
mov borr,0
xor esi,esi
xor ecx,ecx
mov cx,word ptr[edx]
inc cx
@lp: add esi,4
mov ebx,Dword ptr[eax+esi]
sub ebx,borr
mov borr,0
sbb ebx,Dword ptr[edx+esi]
adc borr,0
Mov Dword ptr[eax+esi],ebx
loop @lp
xor ecx,ecx
Mov CX,word ptr[eax]
shl ecx,2
@ll: cmp dWord ptr[eax+ecx],0
jne @fin
dec Word ptr[eax]
sub ecx,4
jz @fin
jmp @ll
@fin: pop esi
pop ebx
end;
{***************************************************************************************}
{ Multiplication of two signed Hug integeres based on standard grade school algorithm }
{***************************************************************************************}
Procedure MulHug(const t1,t2,Result:Hugint);
{t1, t2 and Result Must be diffrent }
var carr:Dword;
sizet1,sizet2:Dword;
asm
push ebx
push esi
push edi
mov bx,word ptr[eax+2]
xor bx,word ptr[edx+2]
mov word ptr[ecx+2],bx
mov esi,1
xor ebx,ebx
mov bx,word ptr[eax]
mov sizet1,ebx
mov bx,word ptr[edx]
mov sizet2,ebx
cmp ebx,0
jz @@isnull
cmp sizet1,0
jnz @@loop1
@@isnull:mov Dword ptr[ecx],0 {one of the operand is null}
jmp @@end4
@@loop1:cmp esi,sizet1
jg @@end1
mov edi,1
mov carr,0
@@loop2:cmp edi,sizet2
jg @@end2
mov ebx,esi
add ebx,edi
dec ebx
cmp esi,1
jg @@next
mov Dword ptr[ecx+ebx*4],0
@@next: push eax
push edx
mov eax,Dword ptr[eax+esi*4]
mul Dword ptr [edx+edi*4]
add eax,carr
adc edx,0
add eax,Dword ptr[ecx+ebx*4]
adc edx,0
mov carr,edx
mov Dword ptr[ecx+ebx*4],eax
pop edx
pop eax
inc edi
jmp @@loop2
@@end2: mov ebx,carr
add edi,esi
dec edi
mov Dword ptr[ecx+edi*4],ebx
inc esi
jmp @@loop1
@@end1: mov bx,word ptr [eax]
add bx,word ptr [edx]
@@end3: mov word ptr[ecx],bx
{--------------------------------------------------------}
mov esi,ecx
xor ecx,ecx
mov cx,word ptr [esi]
shl ecx,2
add ecx,esi
@loop: cmp dword ptr[ecx],0
jne @@end4
cmp ecx,esi
je @@end4
dec Word ptr[esi]
sub ecx,4
jmp @loop
{--------------------------------------------------------}
@@end4: pop edi
pop esi
pop ebx
end;
{***************************************************************************************}
{ Compute the Square of a hug integer }
{***************************************************************************************}
Procedure SqrHug(const a,result:Hugint);
var carr:Dword;
sizet1,sizet2:Dword;
asm
push ebx
push esi
push edi
mov word ptr[edx+2],0
mov ecx,edx
mov esi,1
xor ebx,ebx
mov bx,word ptr[eax]
mov sizet1,ebx
cmp ebx,0
jnz @@loop1
jmp @@end3
@@loop1:cmp esi,sizet1
jg @@end1
mov edi,1
mov carr,0
@@loop2:cmp edi,sizet1
jg @@end2
mov ebx,esi
add ebx,edi
dec ebx
cmp esi,1
jg @@next
mov Dword ptr[ecx+ebx*4],0
@@next: push eax
push edx
mov edx,eax
mov eax,Dword ptr[eax+esi*4]
mul Dword ptr [edx+edi*4]
add eax,carr
adc edx,0
add eax,Dword ptr[ecx+ebx*4]
adc edx,0
mov carr,edx
mov Dword ptr[ecx+ebx*4],eax
pop edx
pop eax
inc edi
jmp @@loop2
@@end2: mov ebx,carr
add edi,esi
dec edi
mov Dword ptr[ecx+edi*4],ebx
inc esi
jmp @@loop1
@@end1: mov bx,word ptr [eax]
shl bx,1
@@end3: mov word ptr[ecx],bx
@@end4: pop edi
pop esi
pop ebx
end;
{***************************************************************************************}
{ Fast division function based on Knuth division algorithm }
{***************************************************************************************}
Procedure DivideModHug(a,b:Hugint;var Q,res:hugint);
Const basediv2=32768;
base=65536;
base_1=65535;
basearray:array [false..true] of dword=(base_1,base);
var d,inv_d:Byte;
bn,bn_1,ri,ri_1,ri_2:Word;
i,j,hb,ha,hr,ptr:integer;
rv,rhat,qhat,right,left,borrow,carry:Dword;
sigA,sigB:Word;
begin
if isnull(b) then begin
Raise ERangeError.Create( 'Illegal Division by 0.');
exit;
end;
{Determination of Q and Res signes}
q.i16[-1]:=a.i16[-1] and not(b.i16[-1]);
Res.i16[-1]:=a.i16[-1] or b.i16[-1];
{a and b comparison }
sigA:=a.i16[-1];
sigB:=b.i16[-1]; {save the signes of a and b}
a.i16[-1]:=0;
b.i16[-1]:=0; {take the absolute value of a and b}
i:=CompareHug(a,b);
a.i16[-1]:=sigA;
b.i16[-1]:=sigB; {restore the signes of a and b}
if i<=0 then begin
if i=0 then begin
res.i16[-2]:=0;
q.i16[-2]:=1;
q.i32[0]:=1;
exit;
end
else begin
for i:=-1 to a.i16[-2]-1 do res.i32[i]:=a.i32[i];
q.i16[-2]:=0;
q.i32[0]:=0;
exit;
end;
end;
ha:=a.i16[-2]shl 1-1;
if a.i16[ha]=0 then dec(ha);
for i:=-1 to a.i16[-2]-1 do res.i32[i]:=a.i32[i];
Res.i32[res.i16[-2]]:=0;
{ We can make a short division }
if (b.i16[-2]=1)and(b.i16[1]=0) then begin
rv:=0;
bn:=b.i16[0];
for i:=ha downto 0 do begin
rhat:=(rv shl 16)+res.i16[i];
Q.i16[i]:=rhat div bn;
rv:=rhat mod bn;
end;
Q.i16[ha+1]:=0;
Q.i16[-2]:=res.i16[-2];
while (Q.i32[Q.i16[-2]-1]=0)and(Q.i16[-2]>0) do dec(Q.i16[-2]);
if rv=0 then Res.i16[-2]:=0
else begin
Res.i16[-2]:=1;
Res.i32[0]:=rv;
end;
exit;
end;
{ Structures initialization }
hb:=b.i16[-2]shl 1-1;
if b.i16[hb]=0 then dec(hb);
ptr:=ha-hb;
hr:=ha+1;
Q.i16[-2]:=ptr shr 1+3;
for i:=0 to Q.i16[-2] do Q.i32[i]:=0;
{Camputation of the scaling factor d >[Base/Bn] }
bn:=b.i16[hb];
bn_1:=b.i16[hb-1];
d:=0;
while (bn<basediv2) do begin
bn:=bn shl 1;
inc(d);
end;
inv_d:=16-d;
{Initialization of the divisor most signifients digits Bn and Bn-1 }
if d>0 then begin
bn:=bn or(b.i16[hb-1] shr (inv_d));
if hb>2 then bn_1:=(bn_1 shl d)or(b.i16[hb-2] shr(inv_d))
else bn_1:=(bn_1 shl d);
end;
{ The main Loop}
for i:=ptr downto 0 do begin
{computation of the divident's most significant digits ri, ri_1 and ri_2 }
ri:=res.i16[hr];
ri_1:=res.i16[hr-1];
ri_2:=res.i16[hr-2];
if d>0 then begin
ri:=(ri shl d)or(ri_1 shr(inv_d));
ri_1:=(ri_1 shl d)or(ri_2 shr (inv_d));
if (hr-3)<0 then ri_2:=ri_2 shl d
else ri_2:=(ri_2 shl d)or(res.i16[hr-3] shr(inv_d));
end;
{Estimation of the quotion digit qhat }
if ri<>bn then begin
rhat:=(ri shl 16)or ri_1;
qhat:=rhat div bn;
rhat:=rhat mod bn;
right:=(rhat shl 16) or ri_2;
left:=bn_1*qhat;
if left>right then begin
dec(qhat);
if (rhat+bn<base) and(left-bn_1>right+(bn shl 16)) then dec(qhat);
end;
end
else begin
qhat:=base_1;
rhat:=ri_1+bn;
right:=(rhat shl 16)or ri_2;
if rhat<base then begin
left:=bn_1*qhat;
if left>right then begin
dec(qhat);
if (rhat+bn<base) and(left-bn_1>right+(bn shl 16)) then dec(qhat);
end;
end;
end;
{multiplication and substration (aiai-1......ai-n):=(aiai-1......ai-n)-qhat*b }
borrow:=base;
carry:=0;
for j:=i to i+hb do begin
carry:=b.i16[j-i]*qhat+(carry shr 16);
borrow:=res.i16[j]+basearray[borrow>=base]-(carry and base_1);
res.i16[j]:=borrow;
end;
borrow:=res.i16[i+hb+1]+basearray[borrow>=base]-(carry shr 16);
res.i16[i+hb+1]:=borrow;
{the quotion digit is stored }
Q.i16[i]:=qhat;
if (borrow<base) then begin
carry:=0;
for j:=i to i+hb do begin
carry:=res.i16[j]+b.i16[j-i]+(carry shr 16);
res.i16[j]:=carry;
end;
res.i16[i+hb+1]:=res.i16[j+1]+(carry shr 16);
dec(Q.i16[i]);
end;
hr:=hr-1;
end;
Res.i16[-2]:=b.i16[-2];
if hb mod 2=0 then Res.i16[hb+1]:=0;
while (res.i32[res.i16[-2]-1]=0)and(Res.i16[-2]>0) do dec(Res.i16[-2]);
while (Q.i32[Q.i16[-2]-1]=0)and(Q.i16[-2]>0) do dec(Q.i16[-2]);
end;
{***************************************************************************************}
{ Fast modulo function based on Knuth division algorithm }
{***************************************************************************************}
Procedure ModuloHug(a,b:Hugint;var res:hugint);
Const basediv2=32768;
base=65536;
base_1=65535;
basearray:array [false..true] of dword=(base_1,base);
var d,inv_d:Byte;
bn,bn_1,ri,ri_1,ri_2:Word;
i,j,hb,ha,hr,ptr:integer;
rv,rhat,qhat,right,left,borrow,carry:Dword;
sigA,sigB:Word;
begin
if isnull(b) then begin
Raise ERangeError.Create( 'Illegal Division by 0.');
exit;
end;
{Determination of Res signes}
Res.i16[-1]:=a.i16[-1] or b.i16[-1];
{a and b comparison }
sigA:=a.i16[-1];
sigB:=b.i16[-1]; {save the signes of a and b}
a.i16[-1]:=0;
b.i16[-1]:=0; {take the absolute value of a and b}
i:=CompareHug(a,b);
a.i16[-1]:=sigA;
b.i16[-1]:=sigB; {restore the signes of a and b}
if i<=0 then begin
if i=0 then begin
res.i16[-2]:=0;
exit;
end
else begin
for i:=-1 to a.i16[-2]-1 do res.i32[i]:=a.i32[i];
exit;
end;
end;
ha:=a.i16[-2]shl 1-1;
if a.i16[ha]=0 then dec(ha);
for i:=-1 to a.i16[-2]-1 do res.i32[i]:=a.i32[i];
Res.i32[res.i16[-2]]:=0;
{ We can make a short division }
if (b.i16[-2]=1)and(b.i16[1]=0) then begin
rv:=0;
bn:=b.i16[0];
for i:=ha downto 0 do begin
rhat:=(rv shl 16)+res.i16[i];
rv:=rhat mod bn;
end;
if rv=0 then Res.i16[-2]:=0
else begin
Res.i16[-2]:=1;
Res.i32[0]:=rv;
end;
exit;
end;
{ Structures initialization }
hb:=b.i16[-2]shl 1-1;
if b.i16[hb]=0 then dec(hb);
ptr:=ha-hb;
hr:=ha+1;
{Camputation of the scaling factor d >[Base/Bn] }
bn:=b.i16[hb];
bn_1:=b.i16[hb-1];
d:=0;
while (bn<basediv2) do begin
bn:=bn shl 1;
inc(d);
end;
inv_d:=16-d;
{Initialization of the divisor most signifients digits Bn and Bn-1 }
if d>0 then begin
bn:=bn or(b.i16[hb-1] shr (inv_d));
if hb>2 then bn_1:=(bn_1 shl d)or(b.i16[hb-2] shr(inv_d))
else bn_1:=(bn_1 shl d);
end;
{ The main Loop}
for i:=ptr downto 0 do begin
{computation of the divident's most significant digits ri, ri_1 and ri_2 }
ri:=res.i16[hr];
ri_1:=res.i16[hr-1];
ri_2:=res.i16[hr-2];
if d>0 then begin
ri:=(ri shl d)or(ri_1 shr(inv_d));
ri_1:=(ri_1 shl d)or(ri_2 shr (inv_d));
if (hr-3)<0 then ri_2:=ri_2 shl d
else ri_2:=(ri_2 shl d)or(res.i16[hr-3] shr(inv_d));
end;
{Estimation of the quotion digit qhat }
if ri<>bn then begin
rhat:=(ri shl 16)or ri_1;
qhat:=rhat div bn;
rhat:=rhat mod bn;
right:=(rhat shl 16) or ri_2;
left:=bn_1*qhat;
if left>right then begin
dec(qhat);
if (rhat+bn<base) and(left-bn_1>right+(bn shl 16)) then dec(qhat);
end;
end
else begin
qhat:=base_1;
rhat:=ri_1+bn;
right:=(rhat shl 16)or ri_2;
if rhat<base then begin
left:=bn_1*qhat;
if left>right then begin
dec(qhat);
if (rhat+bn<base) and(left-bn_1>right+(bn shl 16)) then dec(qhat);
end;
end;
end;
{multiplication and substration (aiai-1......ai-n):=(aiai-1......ai-n)-qhat*b }
borrow:=base;
carry:=0;
for j:=i to i+hb do begin
carry:=b.i16[j-i]*qhat+(carry shr 16);
borrow:=res.i16[j]+basearray[borrow>=base]-(carry and base_1);
res.i16[j]:=borrow;
end;
borrow:=res.i16[i+hb+1]+basearray[borrow>=base]-(carry shr 16);
res.i16[i+hb+1]:=borrow;
if (borrow<base) then begin
carry:=0;
for j:=i to i+hb do begin
carry:=res.i16[j]+b.i16[j-i]+(carry shr 16);
res.i16[j]:=carry;
end;
res.i16[i+hb+1]:=res.i16[j+1]+(carry shr 16);
end;
hr:=hr-1;
end;
Res.i16[-2]:=b.i16[-2];
if hb mod 2=0 then Res.i16[hb+1]:=0;
while (res.i32[res.i16[-2]-1]=0)and(Res.i16[-2]>0) do dec(Res.i16[-2]);
end;
{***************************************************************************************}
{ Comparison of two signed Hug integeres }
{***************************************************************************************}
Function CompareHug(const t1,t2:Hugint):Smallint;
{Result = 0:equal; 1:t1>t2; -1:t1<t2}
asm
push ebx
push edi
xor ebx,ebx
xor cx,cx
mov bx,word ptr[eax+2] { Compare the sign }
cmp bx,Word ptr[edx+2] { jump if the same }
je @@samesig
or bx,1 {0 or 1=1; FFFF or 1 =FFFF=-1; }
mov cx,bx
jmp @@end
@@samesig:xor cx,cx
mov bx,word ptr[eax] {Compare the size}
cmp bx,word ptr[edx]
je @@samsize
sub bx,word ptr[edx]
sets cl
dec cx
not cx
or cx,1
jmp @@end
@@samsize:cmp bx,0
je @@end
@@samsize1:mov edi,Dword ptr[eax+ebx*4]
cmp edi,Dword ptr[edx+ebx*4]
jne @@result
dec bx
jz @@end
jmp @@samsize1
@@Result: setbe cl
dec cx
not cx
or cx,1
@@end: mov result,cx
pop edi
pop ebx
end;
{***************************************************************************************}
{ Right shift of a Hug integer }
{***************************************************************************************}
procedure ShrHug(const a:Hugint;const count:word);
var c:word;
asm
push esi
push edi
cmp Word ptr[eax],0
je @@fin2
{--------------------------------------------------------}
xor ecx,ecx
mov cx,Word ptr[eax]
inc cx
shl cx,2
mov Dword ptr[eax+ecx],0
{--------------------------------------------------------}
xor ecx,ecx
mov cx,count
shr cx,5
mov edi,eax
jz @small
add edi,4
mov esi,edi
//--------------- move bytes
cmp cx,Word ptr[eax]
jl @@lower
mov cx,Word ptr[eax]
@@lower:sub Word ptr[eax],cx
shl cx,2
add esi,ecx
mov cx,Word ptr[eax]
rep movsd
//----------------{rz with zeros}
mov cx,Word ptr[eax]
inc cx
shl cx,2
mov edi,eax
add edi,ecx
mov cx,count
shr cx,5
push eax
xor eax,eax
rep stosd
pop eax
mov edi,eax
//---------------
@small:mov cx,count
and cx,$1f
jz @@fin2
xor edx,edx
mov dx,Word ptr[eax]
shl edx,2
add edx,eax
@@loop:cmp eax,edx
jg @@fin
add eax,4
mov esi,Dword ptr[eax+4]
shrd Dword ptr[eax],esi,cl
mov esi,Dword ptr[eax]
jmp @@loop
@@fin: cmp Dword ptr[eax-4],0
jne @@fin2
dec word ptr[edi]
@@fin2:pop edi
pop esi
end;
{***************************************************************************************}
{ Left shift of a Hug integer }
{***************************************************************************************}
procedure ShlHug(a:Hugint;count:word);
var c:word;
asm
push esi
push ebx
push edi
mov edi,eax
{--------------------------------------------------------}
xor ecx,ecx
mov cx,Word ptr[eax]
inc cx
shl cx,2
mov Dword ptr[eax+ecx],0
{--------------------------------------------------------}
xor ecx,ecx
mov cx,count
shr cx,5
jz @@small
add edi,4
mov esi,edi
//--------------- move bytes
mov bx,cx
mov cx,Word ptr[eax]
shl cx,2
add esi,ecx
shr cx,2
add cx,bx
mov Word ptr[eax],cx
shl cx,2
add edi,ecx
shr cx,2
std
rep movsd
cld
//----------------{rz with zeros}
mov cx,count
shr cx,5
push eax
xor eax,eax
rep stosd
pop eax
mov edi,eax
//---------------
@@small:mov cx,count
and cx,$1f
jz @@fin2
xor edx,edx
mov dx,Word ptr[eax]
and edx,maxlen-1
shl edx,2
add eax,edx
@@loop: cmp eax,edi
je @@fin
mov esi,Dword ptr[eax]
shld Dword ptr[eax+4],esi,cl
sub eax,4
jmp @@loop
@@fin: shl Dword ptr[edi+4],cl
cmp Dword ptr[eax+edx+4],0
je @@fin2
inc Word ptr[edi]
@@fin2: pop edi
pop ebx
pop esi
end;
{***************************************************************************************}
{ Test if a Hug integer is Null }
{***************************************************************************************}
function IsNull(a:Hugint):boolean;
asm
mov ax,word ptr[eax]
or ax,ax
setz al
end;
{***************************************************************************************}
{ Test if a Hug integer is Negative }
{***************************************************************************************}
Function IsNeg(a:Hugint):boolean;
asm
mov ax,word ptr[eax+2]
end;
{***************************************************************************************}
{ Increment a Hug Integer }
{***************************************************************************************}
Procedure IncHug(Var a:Hugint;p:Dword);
asm
cmp word ptr[eax+2],0
je @@start
not word ptr[eax+2]
call DecHug
not word ptr[eax+2]
Ret
@@start:push ebx
xor ebx,ebx
mov ecx,eax
add ecx,4
mov bx, word ptr[eax]
cmp bx,0
je @@endlp
@@non0: add Dword ptr[ecx],edx
jnc @@end
mov edx,1
inc ebx
shl ebx,2
add ebx,eax
@@loop: add ecx,4
cmp ecx,ebx
je @@endlp
add Dword ptr[ecx],1
jnc @@end
jmp @@loop
@@endlp:mov Dword ptr[ecx],edx
inc word ptr[eax]
@@end: pop ebx
end;
{***************************************************************************************}
{ Decrement a Hug Integer }
{***************************************************************************************}
Procedure DecHug(Var a:Hugint;p:Dword);
asm
cmp word ptr[eax+2],0
je @@Start
not word ptr[eax+2]
call IncHug
not word ptr[eax+2]
Ret
@@start:push ebx
push esi
xor ebx,ebx
mov ecx,eax
add ecx,4
mov bx, word ptr[eax]
cmp bx,0
je @@endlp
@@non0: mov esi,Dword ptr[ecx]
sub esi,edx
sub Dword ptr[ecx],edx
jnb @@end
mov edx,1
inc ebx
shl ebx,2
add ebx,eax
@@loop: add ecx,4
cmp ecx,ebx
je @@endlp
sub Dword ptr[ecx],1
jnb @@end
jmp @@loop
@@endlp:not Dword ptr[ecx-4]
mov word ptr[eax+2],$ffff
@@end: cmp Word ptr[eax],1
jne @@fin1
cmp DWord ptr[eax+4],0
jne @@fin1
mov Word ptr[eax],0
@@fin1: cmp Word ptr[eax],2
jne @@fin
cmp DWord ptr[eax+8],0
jne @@fin
mov Word ptr[eax],1
@@fin: pop esi
pop ebx
end;
{***************************************************************************************}
{ Perform a logical "And" between two hug integers }
{***************************************************************************************}
Procedure AndHug(a,b:Hugint;Var Result:Hugint);
asm
push ebx
push esi
xor ebx,ebx
mov bx,word ptr[eax] {eax must represent alwayse the lowest hugint size}
cmp bx,word ptr[edx]
jl @@start
xchg eax,edx
@@start:mov ebx,1
@@loop1:mov esi,Dword ptr[eax+ebx*4]
and esi,Dword ptr[edx+ebx*4]
mov Dword ptr[ecx+ebx*4],esi
inc bx
cmp bx,Word ptr[eax]
jle @@loop1
@@loop2:cmp bx,Word ptr[edx]
jg @@end
mov Dword ptr[ecx+ebx*4],0
inc bx
jmp @@loop2
@@end: mov bx,word ptr[edx]
mov word ptr[ecx],bx
pop esi
pop ebx
end;
{***************************************************************************************}
{ Perform a logical "Or" between two hug integers }
{***************************************************************************************}
Procedure OrHug(a,b:Hugint;Var Result:Hugint);
asm
push ebx
push esi
xor ebx,ebx
mov bx,word ptr[eax] {eax must represent alwayse the lowest hugint size}
cmp bx,word ptr[edx]
jl @@start
xchg eax,edx
@@start:mov ebx,1
@@loop1:mov esi,Dword ptr[eax+ebx*4]
or esi,Dword ptr[edx+ebx*4]
mov Dword ptr[ecx+ebx*4],esi
inc bx
cmp bx,Word ptr[eax]
jle @@loop1
@@loop2:cmp bx,Word ptr[edx]
jg @@end
mov esi,Dword ptr[edx+ebx*4]
mov Dword ptr[ecx+ebx*4],esi
inc bx
jmp @@loop2
@@end: mov bx,word ptr[edx]
mov word ptr[ecx],bx
pop esi
pop ebx
end;
{***************************************************************************************}
{ Perform a logical "Xor" between two hug integers }
{***************************************************************************************}
Procedure XorHug(a,b:Hugint;Var Result:Hugint);
asm
push ebx
push esi
xor ebx,ebx
mov bx,word ptr[eax] {eax must represent alwayse the lowest hugint size}
cmp bx,word ptr[edx]
jl @@start
xchg eax,edx
@@start:mov ebx,1
@@loop1:mov esi,Dword ptr[eax+ebx*4]
xor esi,Dword ptr[edx+ebx*4]
mov Dword ptr[ecx+ebx*4],esi
inc bx
cmp bx,Word ptr[eax]
jle @@loop1
@@loop2:cmp bx,Word ptr[edx]
jg @@end
mov esi,Dword ptr[edx+ebx*4]
mov Dword ptr[ecx+ebx*4],esi
inc bx
jmp @@loop2
@@end: mov bx,word ptr[edx]
mov word ptr[ecx],bx
pop esi
pop ebx
end;
{***************************************************************************************}
{ Perform a logical "Not" between two hug integers }
{***************************************************************************************}
Procedure NotHug(a:Hugint;Var Result:Hugint);
asm
push ebx
push esi
xor ebx,ebx
@@loop: cmp bx,word ptr[eax]
jnl @@end
inc ebx
mov esi,Dword ptr[eax+ebx*4]
not esi
mov Dword ptr[edx+ebx*4],esi
jmp @@loop
@@end: mov bx,word ptr[eax]
mov word ptr[edx],bx
pop esi
pop ebx
end;
{***************************************************************************************}
{ Affect a Hug integer to another }
{***************************************************************************************}
Procedure HCopy(a,b:Hugint);
asm
push esi
push edi
mov esi,eax
mov edi,edx
xor ecx,ecx
mov cx,word ptr[eax]
inc cx
// shl cx,2
rep movsd
pop edi
pop esi
end;
{***************************************************************************************}
{ Compute the Power of a Hug integer to a longword value }
{***************************************************************************************}
Procedure PowerHug(a:hugint;var Result:hugint ;p:longword);
var i:byte;
tmp:hugint;
begin
Result.i32[-1]:=1;
Result.i32[0]:=1;
for i:=31 downto 0 do begin
Mulhug(result,result,tmp);
Hcopy(tmp,result);
if (p and(1 shl i)<>0)then begin
MulHug(Result,a,tmp);
Hcopy(tmp,Result);
end;
end;
end;
{***************************************************************************************}
{ Compute the Moular Power of a Hug integer to a Hug integer }
{***************************************************************************************}
Procedure ModPowerHug(a,b,modulo:hugint;var Result:hugint);
var i:Word;
tmp:hugint;
nbits:Longint;
begin
Result.i32[-1]:=1;
Result.i32[0]:=1;
nbits:=numofbits(b);
for i:=nbits-1 downto 0 do begin
Mulhug(result,result,tmp);
ModuloHug(tmp,modulo,Result);
if (b.i32[i shr 5] and (1 shl (i and 31))<>0) then begin
Mulhug(result,a,tmp);
ModuloHug(tmp,modulo,Result);
end;
end;
end;
{***************************************************************************************}
{ Return the number of bits of a Hug integer }
{***************************************************************************************}
Function NumOfBits(a:Hugint):longint;
var l:longint;
begin
if a.i16[-2]=0 then Result:=0
else begin
l:=(a.i16[-2]-1)shl 5;
l:=l+Trunc(Ln(a.i32[a.i16[-2]-1])/ln(2)+1);
Result:=l;
end;
end;
{***************************************************************************************}
{ Return the number of decimal digits of a Hug integer }
{***************************************************************************************}
Function NumOfDecimalDigits(a:Hugint):Word;
begin
Result:=Length(Hugtostr(a));
end;
{***************************************************************************************}
{ Determine if two Hug integers are relativly prime }
{***************************************************************************************}
Function AreHugsRelativlyPrimes(a,b:Hugint):Boolean;
begin
Result:=(CompareHug(BinaryGCDHug(a,b),One)=0);
end;
{***************************************************************************************}
{ Applay the exetended Euclid Algorithm to find s and t such That as+bt=gcd(a,b) }
{***************************************************************************************}
Procedure ExtendedEuclid(a,b:HugInt;Var s,t:HugInt);
var s0,s1,t0,t1,r,q,aa,bb,tmp:Hugint;
z:string;
begin
Hcopy(a,aa);
Hcopy(b,bb);
Hcopy(One,s0);
s1.i32[-1]:=0;
t0.i32[-1]:=0;
Hcopy(One,t1);
DivideModHug(aa,bb,q,r);
While not IsNull(r) do begin
MulHug(q,s1,tmp);
SubHug(s0,tmp,s);
MulHug(q,t1,tmp);
Subhug(t0,tmp,t);
Hcopy(s1,s0);
Hcopy(s,s1);
Hcopy(t1,t0);
Hcopy(t,t1);
Hcopy(bb,aa);
Hcopy(r,bb);
DivideModHug(aa,bb,q,r);
end;
end;
{***************************************************************************************}
{ Resolve the linear Diopantiene Equation ax+by=c }
{***************************************************************************************}
Function ResolveDiopantiene(a,b,c:Hugint;var x,y:Hugint):Boolean;
var g,r,q,tmp:Hugint;
begin
g:=BinaryGCDHug(a,b);
DivideModHug(c,g,q,r);
if IsNull(r) then begin
ExtendedEuclid(a,b,x,y);
Mulhug(x,q,tmp);
Hcopy(tmp,x);
MulHug(y,q,tmp);
Hcopy(tmp,y);
Result:=True;
end
else begin
Result:=False;
Raise ERangeError.Create( 'No integer solution for this equation ')
end;
end;
{***************************************************************************************}
{ Resolve the linear congurence ax # b[c] and return the number of solutions as result }
{***************************************************************************************}
Function ResolveLinearCongruenceEquation(a,b,c:Hugint;Var Sol:PHugintArray;All:boolean):integer;
var g,q,qq,qqq,r,x,y,tmp,aa,cc,s0,s1:Hugint;
i:integer;
NumSol:Dword;
begin
g:=BinaryGCDHug(a,c);
DivideModHug(b,g,q,r);
Hcopy(a,aa);
if IsNull(r) then begin
ModuloHug(aa,c,aa);
ModuloHug(b,c,b);
{********** ResolveDiopantiene(aa,c,b,x,y) ****}
g:=BinaryGCDHug(aa,c);
DivideModHug(b,g,qq,r);
if IsNull(r) then begin
{*********** ExtendedEuclid(aa,c,x,y) ****************}
Hcopy(c,cc);
Hcopy(One,s0);
s1.i32[-1]:=0;
DivideModHug(aa,cc,qqq,r);
While not IsNull(r) do begin
MulHug(qqq,s1,tmp);
SubHug(s0,tmp,x);
Hcopy(s1,s0);
Hcopy(x,s1);
Hcopy(cc,aa);
Hcopy(r,cc);
DivideModHug(aa,cc,qqq,r);
end;
{********************************************}
Mulhug(x,qq,tmp);
Hcopy(tmp,x);
end
else begin
Result:=-1;
Raise ERangeError.Create( 'No solution for this congruence equation ')
end;
{********************************************}
Numsol:=g.i16[0];
Setlength(Sol^,Numsol);
Hcopy(x,Sol^[0]);
if IsNeg(Sol^[0]) then AddHug(x,c,Sol^[0]);
ModuloHug(Sol^[0],c,Sol^[0]);
if All then begin
Tmp.i32[-1]:=1;
for i:=1 to NumSol-1 do begin
Tmp.i32[0]:=i;
MulHug(c,tmp,aa);
DivideModHug(aa,g,q,r);
AddHug(Sol^[0],q,Sol^[i]);
if IsNeg(Sol^[i]) then AddHug(Sol^[i],c,Sol^[i]);
ModuloHug(Sol^[i],c,Sol^[i]);
end;
end;
Result:=Numsol;
end
else begin
Result:=-1;
Raise ERangeError.Create( 'No solution for this congruence equation ')
end;
end;
{***************************************************************************************}
{ Compute the modular inverse of a Hug integer }
{***************************************************************************************}
Procedure InvModHug(const a,modulo:Hugint;Var Result:Hugint);
var g:Hugint;
Sol:PHugintArray;
begin
new(Sol);
if ResolveLinearCongruenceEquation(a,one,modulo,sol,false)=-1 then Raise ERangeError.Create( 'No inverse for the number with respect to the specified modulo ')
else Hcopy(sol^[0],Result);
end;
{***************************************************************************************}
{ Generate a random binary Hug integer on Numbit Bits }
{***************************************************************************************}
Procedure GetRandomHugOnBits(Var Result:Hugint;Numbits:Integer);
var b,i:integer;
x:Dword;
begin;
b:=Numbits shr 5;
Result.i32[-1]:=b;
For i:=0 to b-1 do Result.i32[i]:=Random($ffffffff);
if Numbits and 31=0 then Result.i32[b-1]:=Result.i32[b-1] or $80000000
else begin
inc(Result.i32[-1]);
x:=(1 shl (Numbits and 31)-1);
Result.i32[b]:=Random(x)or (1 shl ((Numbits and 31)-1));
end;
end;
{***************************************************************************************}
{ Generate a random binary block Lower than Limit }
{***************************************************************************************}
Procedure GetRandomHugLowerThan(var Result:Hugint;Limit:Hugint);
var b,i:integer;
begin
Result.i32[-1]:=Limit.i32[-1];
For i:=0 to Result.i32[-1]-1 do Result.i32[i]:=Random($FFFFFFFF);
ModuloHug(Result,Limit,Result);
while (Result.i32[Result.i16[-2]-1]=0)and(Result.i16[-2]>0) do dec(Result.i16[-2]);
end;
{***************************************************************************************}
{ Initialization of a montgomery strcutur with repect to a modulo (must be odd) }
{***************************************************************************************}
Procedure InitMontgomeryStruct(Modulo:Hugint);
var tmp,Head:DWord;
i:integer;
begin
if (Modulo.i32[0] and 1=0) then begin
Raise ERangeError.Create( 'Montgomery exponentiation can operate on odd modulo only ');
exit;
end;
with MgStruct do begin
r.i32[-1]:=Modulo.i32[-1];
for i:=0 to r.i16[-2]-1 do r.i32[i]:=0;
Head:=Modulo.i32[Modulo.i16[-2]-1];
if Head>$80000000 then begin
nbits:=r.i16[-2]shl 5;
inc(r.i16[-2]);
r.i32[r.i16[-2]-1]:=1;
end
else begin
tmp:=1;
nbits:=0;
while (tmp<=Head) do begin
tmp:=tmp shl 1;
inc(nbits);
end;
r.i32[r.i16[-2]-1]:=tmp;
nbits:=nbits+(r.i16[-2]-1) shl 5;
end;
ResolveDiopantiene(r,modulo,one,Invr,v);
if IsNeg(Invr) then AddHug(Invr,Modulo,Invr);
v.i16[-1]:=not(v.i16[-1]);
if IsNeg(v) then AddHug(v,r,v);
Hcopy(Modulo,Modu);
end;
end;
{***************************************************************************************}
{ Compute a^^b (mod n) using montgomery algorithm (n must be odd number) }
{ the montgomery approache is faster for very big numbers }
{***************************************************************************************}
Procedure MgModPowerHug(a,b,n:Hugint;var Result:Hugint);
var s,r0,m,t,fia:Hugint;
i:integer;
nbits:Longint;
begin
Hcopy(a,fia);
ShlHug(fia,MgStruct.nbits);
ModuloHug(fia,n,fia);
Hcopy(MgStruct.r,r0);
nbits:=NumOfBits(b);
For i:=nbits-1 downto 0 do begin
MulHug(r0,r0,s);
MulHug(s,MgStruct.v,t);
if t.i16[-2]>=MgStruct.r.i16[-2] then begin
t.i16[-2]:=MgStruct.r.i16[-2];
t.i32[t.i16[-2]-1]:=t.i32[t.i16[-2]-1] and (MgStruct.r.i32[MgStruct.r.i16[-2]-1]-1);
end;
while (t.i32[t.i16[-2]-1]=0)and(t.i16[-2]>0)do dec(t.i16[-2]);
MulHug(t,n,m);
AddHug(s,m,r0);
ShrHug(r0,MgStruct.nbits);
if CompareHug(r0,n)<>-1 then SubHug(r0,n,r0);
if (b.i32[i shr 5])and(1 shl (i and 31))<>0 then begin
MulHug(r0,fia,s);
MulHug(s,MgStruct.v,t);
if t.i16[-2]>=MgStruct.r.i16[-2] then begin
t.i16[-2]:=MgStruct.r.i16[-2];
t.i32[t.i16[-2]-1]:=t.i32[t.i16[-2]-1] and (MgStruct.r.i32[MgStruct.r.i16[-2]-1]-1);
end;
while (t.i32[t.i16[-2]-1]=0)and(t.i16[-2]>0)do dec(t.i16[-2]);
MulHug(t,n,m);
AddHug(s,m,r0);
ShrHug(r0,MgStruct.nbits);
if CompareHug(r0,n)<>-1 then SubHug(r0,n,r0);
end;
end;
MulHug(r0,MgStruct.invr,Result);
ModuloHug(Result,n,Result);
end;
{***************************************************************************************}
{ Test If a Hug Integer is a prime number Using the Rabin-Miller Test }
{***************************************************************************************}
Function IsHugPrime(Number:Hugint;Prob:Extended=0.99):Boolean;
Var SminusOne,Base,x,Q,minusOne:Hugint;
nbItr,Rest:Byte;
i:integer;
nbits:Longint;
Prim:Boolean;
begin
Randomize;
if (Number.i16[-2]=1)and(Number.i32[0]=2) then begin
{The number is 2}
Result:=True;
exit;
end;
if (Number.i32[0] and 1=0) then begin
{The number is Even}
Result:=False;
exit;
end;
SubHug(Number,One,minusOne);
ModPowerHug(Two,MinusOne,Number,x);
if CompareHug(x,One)<>0 then begin
Result:=false;
exit;
end;
Hcopy(MinusOne,SMinusOne);
ShrHug(SminusOne,1);
nbItr:=Round((-Ln(1-Prob)/(2*ln(2))));
i:=0;
Prim:=True;
nbits:=Numofbits(Number);
InitMontgomeryStruct(Number);
while (i<=nbItr-1)and(Prim) do begin
GetRandomHugOnBits(Base,12-1);
MgModPowerHug(Base,minusOne,Number,x);
if (x.i16[-2]<>1)or(x.i32[0]<>1) then
Prim:=false
else begin
Hcopy(SMinusOne,Q);
Rest:=0;
While (Rest=0)and Prim do begin
MgModPowerHug(Base,Q,Number,x);
if CompareHug(x,minusOne)=0 then Rest:=1
else begin
if (x.i16[-2]<>1)or(x.i32[0]<>1) then Prim:=False
else begin
Rest:=Q.i32[0] and 1;
ShrHug(Q,1);
end;
end;
end;
end;
i:=i+1;
end;
Result:=Prim;
end;
{***************************************************************************************}
{ Generate a Prime Number using Rabin-Miller Test }
{***************************************************************************************}
Procedure GetHugPrimeOnBits(Numbit:integer;Var Result:Hugint;SkipStep:integer=MaxPrimSteps);
var i,j,k:integer;
b,bb:boolean;
Begin
Randomize;
GetRandomHugOnBits(Result,Numbit);
Result.i32[0]:=Result.i32[0] or 1;
//i:=1;
Repeat
IncHug(Result,2);
bb:=(CompareHug(BinaryGCDHug(ProdFirstPrimes,Result),One)=0);
b:=bb and IsHugPrime(Result,0.99);
// i:=i+1;
// Form11.Edit1.Text:=inttostr(i);
// Application.ProcessMessages;
Until(b);
end;
Procedure GetPrime1(Numbit:integer;init:Hugint;Var Result:Hugint;var nbmont,nbitr:integer);
var i,j,k:integer;
b,bb:boolean;
Begin
Randomize;
Hcopy(init,Result);
Result.i32[0]:=Result.i32[0] or 1;
nbitr:=1;
nbmont:=0;
Repeat
IncHug(Result,2);
bb:=(CompareHug(BinaryGCDHug(ProdFirstPrimes,Result),One)=0);
b:=bb and IsHugPrime(Result,0.99);
if bb then inc(nbmont);
nbitr:=nbitr+1;
Until(b);
end;
Procedure GetPrime2(Numbit:integer;init:Hugint;Var Result:Hugint;var nbmont,nbitr:integer);
var i,j,k,step,Nstep:integer;
b,bb:boolean;
First,Small,x:Hugint;
R:Array[0..2629] of Dword;
Begin
Small.i32[-1]:=1;
Randomize;
Hcopy(init,Result);
Result.i32[0]:=Result.i32[0] or 1;
repeat
i:=1;
b:=true;
while (i<=2629)and b do begin
Small.i32[0]:=firstPrimes[i];
ModuloHug(Result,Small,x);
R[i]:=x.I32[0]*x.I16[-2];
b:=R[i]<>0;
i:=i+1;
end;
if not b then IncHug(Result,2);
until(b);
b:=IsHugPrime(Result,0.99);
NStep:=2;
nbitr:=1;
nbmont:=1;
While not b do begin
j:=1;
b:=true;
Step:=Nstep;
Nstep:=2;
IncHug(Result,Step);
while (j<=2629)do begin
R[j]:=(R[j]+step) mod FirstPrimes[j];
if (R[j]+Nstep)Mod FirstPrimes[j]=0 then Nstep:=Nstep+2;
b:=b and (R[j]<>0);
j:=j+1;
end;
if b then inc(nbmont);
b:=b and IsHugPrime(Result,0.99);
nbitr:=nbitr+1;
end;
end;
Procedure GetPrime3(Numbit:integer;init:Hugint;Var Result:Hugint;var nbmont,nbitr:integer);
var i,j,k,step,Nstep,Pas:integer;
z:Dword;
b,bb:boolean;
First,Small,x:Hugint;
R:Array[0..2629] of Dword;
Begin
Small.i32[-1]:=1;
Randomize;
Hcopy(init,Result);
Result.i32[0]:=Result.i32[0] or 1;
repeat
i:=1;
b:=true;
while (i<=2629)and b do begin
Small.i32[0]:=firstPrimes[i];
ModuloHug(Result,Small,x);
R[i]:=x.I32[0]*x.I16[-2];
b:=R[i]<>0;
i:=i+1;
end;
if not b then IncHug(Result,2);
until(b);
b:=IsHugPrime(Result,0.99);
NStep:=2;
Pas:=0;
nbitr:=1;
nbmont:=1;
While not b do begin
j:=1;
b:=true;
Pas:=Pas+NStep;
IncHug(Result,NStep);
Nstep:=2;
while (j<=2629) and b do begin
b:= b and ((R[j]+Pas) mod FirstPrimes[j]<>0);
if (R[j]+pas+Nstep)Mod FirstPrimes[j]=0 then Nstep:=Nstep+2;
j:=j+1;
end;
if b then inc(nbmont);
b:=b and IsHugPrime(Result,0.99);
nbitr:=nbitr+1;
end;
end;
{***************************************************************************************}
{ Generate a Prime Number on the form 4k+3 to be used for Blum-Glodwasser Cipher }
{***************************************************************************************}
Procedure GetHug4Kp3PrimeOnBits(Numbit:integer;Var Result:Hugint;SkipStep:integer=MaxPrimSteps);
var i,j,k:integer;
Two,Four,Small,Res:Hugint;
b:boolean;
Begin
Randomize;
Small.i16[-2]:=1;
Small.i16[-1]:=0;
GetRandomHugOnBits(Result,Numbit);
Result.i32[0]:=Result.i32[0] or 1;
if (Result.i32[0] and 3<>3) then Result.i32[0]:=Result.i32[0]+2;
k:=SkipStep;
Repeat
IncHug(Result,4);
j:=1;
b:=True;
While(j<=k)and b do begin
Small.i32[0]:=FirstPrimes[j];
ModuloHug(Result,Small,Res);
b:=b and not Isnull(Res);
j:=j+1;
end;
b:=b and IsHugPrime(Result,0.99);
Until(b);
end;
{***************************************************************************************}
{ Generate a Strong Prime Number using Gordon Algorithm }
{***************************************************************************************}
Procedure GetStrongHugPrimeOnBits(Numbit:integer;Var Result:Hugint);
var s,t,r,z,ps,p:Hugint;
tmp,i:Hugint;
bb:boolean;
Begin
GetRandomHugOnBits(s,Numbit div 2-8);
GetRandomHugOnBits(t,Numbit div 2-8);
Hcopy(one,i);
Repeat
Mulhug(t,i,tmp);
ShlHug(tmp,1);
AddHug(tmp,One,r);
IncHug(i,1);
until (IsHugPrime(r,0.99));
Subhug(r,two,tmp);
ModPowerHug(s,tmp,r,z);
MulHug(z,s,tmp);
ShlHug(tmp,1);
SubHug(tmp,One,ps);
Hcopy(One,i);
Repeat
Mulhug(r,s,tmp);
shlhug(tmp,1);
Mulhug(tmp,i,p);
addhug(p,ps,p);
shlhug(i,1);
until(NumOfBits(p)>Numbit);
bb:=IsHugPrime(p,0.999) ;
while not bb do begin
IncHug(i,1);
Mulhug(r,s,tmp);
shlhug(tmp,1);
mulhug(tmp,i,p);
addhug(ps,p,p);
bb:=IsHugPrime(p,0.99) ;
end;
Result:=p;
end;
{***************************************************************************************}
{ Generate a Safe Prime Number (2p+1) }
{***************************************************************************************}
Procedure GetSafePrimeOnBits(Numbit:integer;Var Result:Hugint);
var i,j,k,step,Nstep:integer;
b,p1,p2:boolean;
Q,Q_1,Small,x:Hugint;
R:Array[0..2629] of Dword;
Begin
Small.i32[-1]:=1;
GetRandomHugOnBits(Q,Numbit-1);
Q.i32[0]:=Q.i32[0] or 1;
k:=0;
Repeat
repeat
i:=1;
b:=(Q.i32[0] and 3=1);
while (i<=2629)and b do begin
Small.i32[0]:=firstPrimes[i];
ModuloHug(Q,Small,x);
R[i]:=x.I32[0]*x.I16[-2];
b:=(R[i]<>0)and(R[i]<>(firstPrimes[i]-1)div 2);
i:=i+1;
end;
if not b then IncHug(Q,2);
until(b);
b:=IsHugPrime(Q,0.99);
NStep:=2;
While not b do begin
j:=1;
b:=(Q.i32[0] and 3=1);
Step:=Nstep;
Nstep:=2;
IncHug(Q,Step);
while (j<=2629)do begin
R[j]:=(R[j]+step) mod FirstPrimes[j];
if (R[j]+Nstep)Mod FirstPrimes[j]=0 then Nstep:=Nstep+2;
b:=b and (R[j]<>0) and (R[j]<>(FirstPrimes[j]-1) div 2);
j:=j+1;
end;
b:=b and IsHugPrime(Q,0.99);
end;
Hcopy(Q,Result);
Hcopy(Q,Q_1);
ShlHug(Result,1);
IncHug(Result,1);
DecHug(Q_1,1);
ShrHug(Q_1,1);
inchug(Q,2);
j:=1;
b:=true;
while (j<=2629)do begin
R[j]:=(2*R[j]+1) mod FirstPrimes[j];
b:=b and (R[j]<>0);
j:=j+1;
end;
p1:=b and IsHugPrime(Result,0.99);
j:=1;
b:=true;
while (j<=2629)do begin
R[j]:=(2*R[j]-1) mod FirstPrimes[j];
b:=b and (R[j]<>0);
j:=j+1;
end;
p2:=p1 or (b and IsHugPrime(Q_1,0.99));
k:=k+1;
// Form1.Edit1.Text:=inttostr(k);
// Application.ProcessMessages;
Until(p1 or p2);
if p2 and not(p1)then Hcopy(Q_1,Result);
end;
{***************************************************************************************}
{ Generate a multiple Safe Prime Number (2*k*p+1) }
{***************************************************************************************}
Procedure GetMulSafePrimeOnBits(Numbit:integer;Var Res:Hugint;facts:PSafeStruct=nil);
var i,j,k,step,Nstep:Word;
fac:Dword;
b,p1,p2:boolean;
Q,Q2,Small,x:Hugint;
R,RQ:Array[0..2629] of Dword;
begin
Small.i32[-1]:=1;
GetRandomHugOnBits(Q,Numbit-16);
Q.i32[0]:=Q.i32[0] or 1;
k:=0;
Repeat
repeat
i:=1;
b:=True;
while (i<=2629)and b do begin
Small.i32[0]:=firstPrimes[i];
ModuloHug(Q,Small,x);
R[i]:=x.I32[0]*x.I16[-2];
b:=(R[i]<>0);
i:=i+1;
end;
if not b then IncHug(Q,2);
until(b);
b:=IsHugPrime(Q,0.99);
NStep:=2;
While not b do begin
j:=1;
b:=True;
Step:=Nstep;
Nstep:=2;
IncHug(Q,Step);
while (j<=2629)do begin
R[j]:=(R[j]+step) mod FirstPrimes[j];
if (R[j]+Nstep)Mod FirstPrimes[j]=0 then Nstep:=Nstep+2;
b:=b and (R[j]<>0);
j:=j+1;
end;
b:=b and IsHugPrime(Q,0.99);
end;
Hcopy(Q,Res);
ShlHug(Res,1);
Hcopy(Res,Q2);
ShlHug(Res,15);
incHug(Res,1);
fac:=65536;
j:=1;
b:=true;
while (j<=2629)do begin
Small.i32[0]:=firstPrimes[j];
ModuloHug(Res,Small,x);
R[j]:=x.I32[0]*x.I16[-2];
ModuloHug(Q2,Small,x);
RQ[j]:=x.I32[0]*x.I16[-2];
b:=b and (R[j]<>0);
j:=j+1;
end;
p1:=b and IsHugPrime(Res,0.99);
While(Not p1)and(fac<4294967296) do begin
AddHug(Res,Q2,Res);
fac:=fac+1;
j:=1;
b:=true;
while (j<=2629)do begin
R[j]:=(R[j]+RQ[j]) mod FirstPrimes[j];
b:=b and (R[j]<>0);
j:=j+1;
end;
P1:=b and IsHugPrime(Res);
//Form1.Edit2.Text:=inttostr(fac);
//Application.ProcessMessages;
end;
inchug(Q,2);
Until(p1);
if Facts<>Nil then begin
Facts^.k:=fac;
ShrHug(Q2,1);
Hcopy(Q2,Facts^.Q);
end;
end;
{***************************************************************************************}
{ Generate a multiple Safe Prime Number (2*k*p+1) With the corresponding Generator }
{***************************************************************************************}
Procedure GetSafeWithGenerator(Numbit:integer;Var Result,Generator:Hugint);
var k:DWORD;
Factors:Array of Dword;
sq,i:Dword;
Fstruct:PSafeStruct;
Q,res,pow,Small:Hugint;
isgen:boolean;
begin
New(Fstruct);
GetMulSafePrimeOnBits(Numbit,Result,Fstruct);
k:=Fstruct^.k;
sq:=Round(sqrt(k))+1;
Setlength(Factors,1);
Factors[0]:=2;
For i:=3 to sq do begin
if k mod i=0 then begin
Setlength(Factors,Length(Factors)+1);
Factors[Length(Factors)-1]:=i;
end;
end;
InitMontgomeryStruct(Result);
Repeat
GetRandomHugLowerThan(Generator,Result);
DecHug(Generator,1);
DivideModHug(Result,Fstruct^.Q,Q,Res);
MgModPowerHug(Generator,Q,Result,pow);
isgen:=not (CompareHug(Pow,One)=0);
i:=0;
Small.i32[-1]:=1;
while (isgen)and(i<=Length(Factors)-1) do begin
Small.i32[0]:=factors[i];
DivideModHug(Result,Small,Q,Res);
MgModPowerHug(Generator,Q,Result,pow);
isgen:=not (CompareHug(Pow,One)=0);
i:=i+1;
end;
until isgen;
end;
{***************************************************************************************}
{ Find the Integer Square root of a positive Hug Integer }
{***************************************************************************************}
Procedure SqrtHug(Number:Hugint;Var Result:Hugint);
var a,a_plus1,b,tmp:Hugint;
nbit:Longint;
stop:boolean;
begin
if IsNull(Number) then begin
Result.i32[-1]:=0;
exit;
end
else if IsNeg(Number) then begin
Raise ERangeError.Create( 'Can''t compute the square root of a negative number');
exit;
end;
nbit:=NumOfBits(Number) shr 1;
a.i32[-1]:=nbit shr 5+1;
a.i32[a.i32[-1]-1]:=1 shl (nbit and 31);
Stop:=false;
while not Stop do begin
DivideModHug(Number,a,b,tmp);
AddHug(a,b,b);
ShrHug(b,1);
SubHug(b,a,tmp);
tmp.i16[-1]:=0;
if CompareHug(Tmp,Two)=-1 then begin
MulHug(a,a,tmp);
if CompareHug(tmp,Number)<1 then begin
AddHug(a,One,a_plus1);
MulHug(a_plus1,a_plus1,tmp);
if CompareHug(tmp,Number)=1
then begin
Hcopy(a,Result);
Stop:=true;
end;
end;
end;
Hcopy(b,a);
end;
end;
{***************************************************************************************}
{ Find the Square of a hug integer "a" modulo "modulo" (Result^2==a[Modulo]) }
{***************************************************************************************}
Procedure ModSquareHug(a,modulo:Hugint;Var Result:Hugint);
var e1,e2,g,MinusOne,MinusOneDiv2,tmp,tmp2,tmp3:Hugint;
begin
{First we test if Modu is a Prime numbre}
if not IsHugPrime(modulo) then begin
Raise ERangeError.Create( 'The modulo must be a prime number ');
exit;
end;
{Then we test if a is a quadratic Residu modulo modu using Euler's Criterion}
SubHug(Modulo,One,tmp);
ShrHug(tmp,1);
ModPowerHug(a,tmp,modulo,g);
if (CompareHug(a,Modulo)<>0)and(CompareHug(g,One)<>0) then begin
Raise ERangeError.Create( Hugtostr(a)+' is not a Quadratic residu Modulo '+Hugtostr(modulo));
exit;
end;
if (Modulo.i32[0] and 3)=3 then begin
{The modulo is in the form 4k+3 so direct computation...}
AddHug(Modulo,One,Tmp);
ShrHug(Tmp,2);
ModPowerHug(a,tmp,modulo,Result);
end
else begin
{The Modulo is in the form 4k+1 so we use the Tonelli-Shanks Algorithm}
SubHug(modulo,One,MinusOne);
Hcopy(MinusOne,MinusOneDiv2);
ShrHug(Minusonediv2,1);
Hcopy(Two,g);
ModPowerHug(g,Minusonediv2,Modulo,tmp);
while CompareHug(tmp,MinusOne)<>0 do begin
incHug(g,1);
ModPowerHug(g,Minusonediv2,Modulo,tmp);
end;
Hcopy(MinusOneDiv2,e1);
e2.i32[-1]:=0;
While(e1.i32[0] and 1=0) do begin
ShrHug(e1,1);
ModPowerHug(a,e1,Modulo,tmp);
ModPowerHug(g,e2,Modulo,tmp2);
MulHug(tmp,tmp2,tmp3);
ModuloHug(tmp3,Modulo,tmp);
if CompareHug(tmp,MinusOne)=0 then begin
Addhug(e2,MinusOneDiv2,e2);
ShrHug(e2,1);
end;
end;
incHug(e1,1);
ShrHug(e1,1);
ModPowerHug(a,e1,Modulo,tmp2);
ShrHug(e2,1);
ModPowerHug(g,e2,Modulo,tmp3);
MulHug(tmp2,tmp3,tmp);
ModuloHug(tmp,Modulo,Result);
end;
end;
{***************************************************************************************}
{ Convert a Hug integer to a decimal string }
{***************************************************************************************}
Function HugToStr(a:Hugint):String;
var s,sig:string;
i:integer;
rem,bilion:Hugint;
SaveSig:Word;
begin
bilion.i32[-1]:=1;
bilion.i32[0]:=1000000000;
result:='';
SaveSig:=a.i16[-1];
if SaveSig<>0 then sig:='-' else sig:='';
a.i16[-1]:=0;
repeat
dividemodhug(a,bilion,a,rem);
if rem.i16[-2]<>0 then s:=inttostr(rem.i32[0])
else s:='0';
while length(s)<9 do s:='0'+s;
Result:=s+Result;
until(a.i32[-1]and $ffff=1)or(a.i32[-1]and $ffff=0)and(a.i32[0]<=999999999);
Result:=inttostr(a.i32[0])+result;
while (result[1]='0')and(length(Result)>1) do result:=copy(result,2,length(result));
result:=sig+result;
a.i16[-1]:=SaveSig;
end;
{***************************************************************************************}
{ Convert a Hug integer to a Hexadecimal string }
{***************************************************************************************}
Function HugToHex(a:Hugint):String;
var i:integer;
s:string;
begin
Result:='';
for i:=0 to a.i16[-2]-1 do begin
s:=inttohex(a.i32[i],8);
Result:=s+Result;
end;
if result<>'' then While (Result[1]='0')and (length(Result)>1) do
Result:=Copy(Result,2,length(Result))
else Result:='0';
end;
{***************************************************************************************}
{ Convert a Hexadecimal string to a Hug integer }
{***************************************************************************************}
Function HexToHug(s:string):Hugint;
var i,len:word;
tmp:Hugint;
begin
Hcopy(Zero,Result);
tmp.i16[-1]:=0;
len:=length(s);
for i:=len downto 1 do begin
tmp.i16[-2]:=1;
if s[i] in ['0'..'9'] then tmp.i32[0]:=Ord(s[i])-ord('0')
else if s[i]in ['A','B','C','D','E','F'] then tmp.i32[0]:=Ord(s[i])-ord('A')+10
else if s[i]in ['a','b','c','d','e','f'] then tmp.i32[0]:=Ord(s[i])-ord('a')+10
else tmp.i32[0]:=0;
ShlHug(tmp,(len-i)*4);
AddHug(Result,tmp,Result);
end;
end;
{***************************************************************************************}
{ Convert a binary string to a Hug integer }
{***************************************************************************************}
Function Hugtobin(a :Hugint):string;
var i,j:integer;
begin
Result:='';
for i:=0 to a.i16[-2] do begin
for j:=0 to 31 do if (a.i32[i] and (1 shl j))=0 then result:='0'+result
else result:='1'+result;
end;
if result<>'' then While (Result[1]='0')and (length(Result)>1) do Result:=Copy(Result,2,length(Result))
else Result:='0';
end;
{***************************************************************************************}
{ Convert a Hug integer to a binary string }
{***************************************************************************************}
Function BintoHug (s:string):Hugint;
var i:integer;
begin
Hcopy(Zero,Result);
for i:=0 to length(s) do begin
shlhug(Result,1);
if s[i]='1' then orhug(Result,One,Result);
end;
end;
{***************************************************************************************}
{ Convert a Hug integer to an Octal string }
{***************************************************************************************}
Function HugtoOctal(a:Hugint):string;
var i,j:integer;
w:word;
begin
Result:='';
While not(IsNull(a)) do begin
w:=a.i32[0] and 7;
Result:=inttostr(w)+Result;
ShrHug(a,3);
end;
if Result<>'' then While (Result[1]='0')and (length(Result)>1) do Result:=Copy(Result,2,length(Result))
else result:='0';
end;
{***************************************************************************************}
{ Convert an octal string to a Hug integer }
{***************************************************************************************}
Function OctaltoHug(S:string):Hugint;
var i,len:integer;
tmp:Hugint;
begin
Hcopy(Zero,Result);
tmp.i16[-1]:=0;
len:=length(s);
for i:=len downto 1 do begin
tmp.i16[-2]:=1;
if s[i] in ['0'..'7'] then tmp.i32[0]:=Ord(s[i])-ord('0')
else tmp.i32[0]:=0;
ShlHug(tmp,(len-i)*3);
AddHug(Result,tmp,Result);
end;
end;
{***************************************************************************************}
{ Convert decimal string to a Hug integer }
{***************************************************************************************}
Function StrToHug(s:string):Hugint;
var str:string;
i:integer;
signe:Dword;
million,tmp,t,z,h:Hugint;
begin
signe:=0;
if s[1]='-' then begin
s:=copy(S,2,length(s));
signe:=$ffff0000;
end;
Million.i32[-1]:=1;
t.i32[-1]:=1;
tmp.i32[-1]:=1;
t.i32[0]:=1;
Million.i32[0]:=1000000;
if length(s)<=7 then begin
Result.i32[0]:=strtoint(s);
if Result.i32[0]<>0 then Result.i32[-1]:=1
else Result.i32[-1]:=0;
end
else begin
str:=s;
Result.i32[-1]:=1;
Result.i32[0]:=0;
repeat
tmp.i32[0]:=strtoint(copy(str,length(str)-5,length(str)));
str:=copy(str,1,length(str)-6);
MulHug(tmp,t,z);
Addhug(Result,z,Result);
mulhug(t,million,h);
Hcopy(h,t);
until(str='');
end;
Result.i32[-1]:=Result.i32[-1] or signe;
end;
{***************************************************************************************}
{ Compute the GCD of two Hug integers using classical Euclid Algo }
{***************************************************************************************}
Function EuclidGCDHug(a,b:HugInt):HugInt;
var Qtmp:Hugint;
Swap:array[0..1]of HugInt;
i,j,k :Byte;
id:integer;
begin
Result.i32[-1]:=0;
if (IsNull(a))or (IsNull(b)) then exit
else begin
Hcopy(a,Swap[0]);
Hcopy(b,Swap[1]);
i:=0;
j:=1;
Repeat
DivideModhug(Swap[i],Swap[j],Qtmp,Swap[i]);
if not IsNull(Swap[i]) then begin
k:=i;
i:=j;
j:=k;
end;
until IsNull(Swap[i]);
for id:=-1 to Swap[j].i16[-2] do Result.i32[id]:=Swap[j].i32[id];
end;
end;
{***************************************************************************************}
{ Compute the GCD of two Hug integers using Binary Algo }
{***************************************************************************************}
Function BinaryGCDHug(a,b:HugInt):HugInt;
var u,v,g,t,z:Hugint;
n1,n2:Longint;
begin
Result.i32[-1]:=0;
if (IsNull(a))or (IsNull(b)) then exit
else begin
Hcopy(a,u);
Hcopy(b,v);
Hcopy(One,g);
While(u.i32[0] and 1=0)and(v.i32[0] and 1=0) do begin
ShrHug(u,1);
ShrHug(v,1);
ShlHug(g,1);
end;
While not(isNull(v)) do begin
if (u.i32[0] and 1=0) then ShrHug(u,1)
else if (v.i32[0] and 1=0) then ShrHug(v,1)
else begin
SubHug(u,v,t);
t.i16[-1]:=0;
ShrHug(t,1);
if CompareHug(u,v)<=0 then Hcopy(t,v) else Hcopy(t,u);
end;
end;
MulHug(u,g,Result);
end;
end;
{***************************************************************************************}
{ Resolve a multi-modulo linear congruences system using the "Chinees Remainder Theorem"}
{ The Result is given by the value:Result[0] and the modulo:Result[1] }
{***************************************************************************************}
Procedure ResolveMultiLinearCongSystems(Resid,Modu:HugintArray;WithoutVerif:Boolean;Var Result:HugintArray);
var i,j:integer;
prim:boolean;
m,tmp,tmp1:Hugint;
begin
if (Length(Resid)<>Length(Modu))or(Length(Resid)=0)or(length(Modu)=0) then
Raise ERangeError.Create( 'Invalid linear congruences system parametres ');
if Not WithoutVerif then begin
i:=0;
j:=0;
prim:=true;
while (i<=length(Modu)-1)and prim do begin
while (j<=Length(Modu)-1) and prim do begin
if i<>j then
Prim:=Prim and AreHugsRelativlyPrimes(Modu[i],Modu[j]);
j:=j+1;
end;
i:=i+1;
end;
if not prim then Raise ERangeError.Create( 'Modulus are not pirwise relativly primes');
end;
Setlength(Result,2);
Result[1].i32[-1]:=1;
Result[1].i32[0]:=1;
For i:=0 to Length(Modu)-1 do begin
MulHug(Result[1],Modu[i],tmp);
Hcopy(tmp,Result[1]);
end;
Result[0].i32[0]:=0;
for i:=0 to length(Modu)-1 do begin
DivideModHug(Result[1],Modu[i],m,tmp);
InvModHug(m,Modu[i],tmp);
MulHug(tmp,m,tmp1);
MulHug(tmp1,Resid[i],tmp);
AddHug(Result[0],tmp,Result[0]);
end;
ModuloHug(Result[0],Result[1],Result[0]);
end;
{***************************************************************************************}
{ Resolve a quadratic Congruence of the Form ax²+bx+c # 0[Mod p1*p2*...*pn] }
{ The Result is given by the modulo:n and the solution values: Result[] }
{***************************************************************************************}
Procedure ResolveQuadraticCongruenceEquation(a,b,c:Hugint;ModuloFactors:HugintArray;var n:Hugint;Var Result:HugintArray);
var i,j:Integer;
Prime:Boolean;
IndivSolutions,TmpArray,TmpRes:HugintArray;
aa,bb,inva,inv2,tmp:Hugint;
NbSol:Word;
begin
if Length(ModuloFactors)=0 then Raise ERangeError.Create( 'Invalid Quadratic congruence parametres');
i:=0;
j:=0;
Prime:=True;
{p1,p2,....,pn must be distincts primes }
While (i<=Length(ModuloFactors)-1) and (Prime) do begin
While (j<=Length(ModuloFactors)-1) and Prime do begin
if i<>j then Prime:=Prime and AreHugsRelativlyPrimes(ModuloFactors[i],ModuloFactors[j]);
j:=j+1;
end;
i:=i+1;
end;
if not Prime then Raise ERangeError.Create( 'Factors are not pirewise relativly primes...');
n.i32[-1]:=1;
n.i32[0]:=1;
For i:=0 to length(moduloFactors)-1 do begin
MulHug(n,ModuloFactors[i],tmp);
Hcopy(tmp,n);
end;
InvModHug(a,n,Inva);
InvModHug(Two,n,Inv2);
MulHug(inv2,b,tmp);
MulHug(tmp,tmp,aa);
MulHug(aa,inva,tmp);
SubHug(tmp,c,tmp);
MulHug(inva,tmp,aa);
ModuloHug(aa,n,aa);
MulHug(inv2,inva,tmp);
MulHug(tmp,b,bb);
ModuloHug(bb,n,bb);
SetLength(IndivSolutions,0);
For i:=0 to Length(ModuloFactors)-1 do begin
SetLength(IndivSolutions,Length(IndivSolutions)+1);
ModSquareHug(aa,ModuloFactors[i],IndivSolutions[High(IndivSolutions)]);
end;
NbSol:=Round(Exp(Length(IndivSolutions)*Ln(2)));
Setlength(Tmparray,Length(IndivSolutions));
Setlength(Result,NbSol);
Setlength(TmpRes,2);
For i:=0 to NbSol-1 do begin
For j:=0 to Length(TmpArray)-1 do begin
Hcopy(IndivSolutions[j],TmpArray[j]);
if (i and (1 shl j))<>0 then
SubHug(ModuloFactors[j],TmpArray[j],TmpArray[j]);
end;
ResolveMultiLinearCongSystems(TmpArray,ModuloFactors,True,TmpRes);
Hcopy(TmpRes[0],Result[i]);
SubHug(Result[i],bb,Result[i]);
while IsNeg(Result[i]) do AddHug(Result[i],n,Result[i]);
end;
end;
{***************************************************************************************}
{ Factorize a Hug integer to the product of it's prime factores }
{ using the probabilistic Monte-Carlo Algorithm }
{***************************************************************************************}
Procedure MonteCarloFactorizeHug(a:Hugint;var factors:HugintArray;Memo:Tmemo=nil);
var mi,m2i,tmp,g,Q:Hugint;
i:Dword;
begin
if IsHugPrime(a) then begin
Setlength(Factors,Length(Factors)+1);
Hcopy(a,Factors[Length(Factors)-1]);
end
else begin
GetRandomHugOnBits(mi,NumOfBits(a));
Hcopy(mi,m2i);
i:=0;
Repeat
MulHug(mi,mi,tmp);
IncHug(tmp,1);
ModuloHug(tmp,a,mi);
MulHug(m2i,m2i,tmp);
IncHug(tmp,1);
ModuloHug(tmp,a,m2i);
MulHug(m2i,m2i,tmp);
IncHug(tmp,1);
ModuloHug(tmp,a,m2i);
SubHug(m2i,mi,tmp);
tmp.i16[-1]:=0;
i:=i+1;
if (IsNull(tmp)) then begin
GetRandomHugLowerThan(mi,a);
Hcopy(mi,m2i);
Hcopy(one,g);
end
else g:=BinaryGCDHug(a,tmp);
if (i>100000000) and(CompareHug(g,one)=0) then begin
GetRandomHugLowerThan(mi,a);
Hcopy(mi,m2i);
i:=0;
end;
until(CompareHug(g,One)<>0)and(CompareHug(g,a)<>0);
DivideModHug(a,g,Q,tmp);
if not IsHugPrime(g) then MonteCarloFactorizeHug(g,factors)
else begin
Setlength(Factors,Length(Factors)+1);
Hcopy(g,Factors[Length(Factors)-1]);
if Memo<>nil then begin
Memo.Lines.Add(Hugtostr(g)+' is a prime factor');
Memo.Lines.Add('Now factorizing :'+hugtostr(Q));
end;
end;
if not IsHugPrime(Q) then MonteCarloFactorizeHug(Q,factors)
else begin
Setlength(Factors,Length(Factors)+1);
Hcopy(g,Factors[Length(Factors)-1]);
if Memo<>nil then Memo.Lines.Add(Hugtostr(Q)+' is a prime factor');
end;
end;
end;
{***************************************************************************************}
{ Factorize a Hug integer to the product of it's prime factores }
{ using the Pollard P-1 Algorithm }
{***************************************************************************************}
Procedure PollardFactorizeHug(a:hugint;var Factors:HugintArray;Memo:Tmemo=nil);
var r,tmp,i,g:Hugint;
begin
if IsHugPrime(a) then begin
Setlength(Factors,Length(Factors)+1);
Hcopy(a,Factors[Length(Factors)-1]);
if Memo<>nil then Memo.Lines.Add(Hugtostr(a)+' is a prime factor');
end
else begin
repeat
GetRandomHugOnBits(r,NumOfBits(a));
until (CompareHug(BinaryGCDHug(a,r),One)=0);
Hcopy(One,i);
repeat
ModPowerHug(r,i,a,tmp);
IncHug(i,1);
Hcopy(tmp,r);
DecHug(tmp,1);
if IsNull(tmp) then begin
repeat
GetRandomHugOnBits(r,NumOfBits(a));
until (CompareHug(BinaryGCDHug(a,r),One)=0);
Hcopy(One,i);
HCopy(One,g);
end
else g:=BinaryGCDHug(tmp,a);
until(CompareHug(g,One)<>0);
DivideModHug(a,g,r,tmp);
if not IsHugPrime(g) then PollardFactorizeHug(g,Factors)
else begin
Setlength(Factors,Length(Factors)+1);
Hcopy(g,Factors[Length(Factors)-1]);
if Memo<>nil then begin
Memo.Lines.Add(Hugtostr(g)+' is a prime factor');
Memo.Lines.Add('Now factorizing :'+hugtostr(r));
end;
end;
if not IsHugPrime(r) then PollardFactorizeHug(r,Factors)
else begin
Setlength(Factors,Length(Factors)+1);
Hcopy(g,Factors[Length(Factors)-1]);
if Memo<>nil then Memo.Lines.Add(Hugtostr(r)+' is a prime factor');
end;
end;
end;
Begin
Zero:=Strtohug('0');
One:=Strtohug('1');
Two:=Strtohug('2');
Three:=Strtohug('3');
Four:=Strtohug('4');
ProdFirstPrimes.i32[-1]:=1059;
For ind:=0 to 1058 do ProdFirstPrimes.i32[ind]:=_ProdFirstPrimes[ind];
end.
Conclusion :
L'application en jointe est une démonstration utilisant cette bibliothèque. J'ai attacher un fichier Zip Avec le source qui contien deux démonstrations, s'il n'est pas téléchargable alors mail me at kamel_mh@yahoo.fr.
Merci.
Codes Sources
A voir également
- Manipulation trés rapide de trés grandes nombres entiers + fonctions de base de l'arithmétique modulaire,nombres premiers...
- Générateur de nombre ✓ - Forum PHP
- Base de données javascript - Conseils pratiques -Javascript
- Nombres en árabe ✓ - Forum Delphi / Pascal
- Traduction d'un nombre reel en chiffres en arabe ✓ - Forum Delphi / Pascal
- Nombre de semaine par an - Forum C / C++ / C++.NET
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.