Manipulation trés rapide de trés grandes nombres entiers + fonctions de base de l'arithmétique modulaire,nombres premiers...

Soyez le premier à donner votre avis sur cette source.

Vue 5 418 fois - Téléchargée 873 fois

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.

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

Ajouter un commentaire Commentaires
Messages postés
615
Date d'inscription
dimanche 13 août 2006
Statut
Membre
Dernière intervention
13 décembre 2018
3
Je viens de découvrir ton travail...
Super travail Kamel!
ça m'aide à comprendre l'assembleur en plus dont j'ai besoin en ce moment!
Messages postés
637
Date d'inscription
mardi 22 avril 2003
Statut
Membre
Dernière intervention
9 janvier 2017

Bonjour,

Bravo cela fonctionne très bien sur D7.

Un détail:
sur les deux boutons " a^^b " et " a^^b mod "
c'est un seul signe " ^ " qu'il faut voir .

Merci
Salutations
Messages postés
267
Date d'inscription
mardi 24 juillet 2007
Statut
Membre
Dernière intervention
7 juin 2018

Re-Bonjour Kamel78,

Je retire ce que viens d'écrire : Je n'avais pas remarqué le Sous-ZIP nommé 'HugeComputationDemo.zip' qui était inclu dans le ZIP et qui contient les nouveautés citées dans l'historique ci-dessus.

A mon avis il est inutile de surcharger le ZIP avec le sous-dossier nommé 'Demo old'.

Par ailleurs je signale : [Erreur fatale] GeneratorDemo.pas(7): Fichier non trouvé : 'PBNumEdit.dcu'
Vérification faite on peut virer PBNumEdit du uses car PBNumEdit1 n'est pas utilisé.

En plus on trouve Demo.dpr et HugeComputation.dpr qui semblent être la même chose : C'est lequel des deux qui est le meilleur ???

Pour info : le Help.doc (daté du 01/06/2012 ???) malgré cette datation bizarre n'est pas à jour.

A+.
Messages postés
267
Date d'inscription
mardi 24 juillet 2007
Statut
Membre
Dernière intervention
7 juin 2018

Bonjour Kamel78,

Je viens de télécharger la dernière version du 11 mai 2012 02:48:59 dont l'historique ci-dessus dit :
"ajout des fonctions : 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; Ajout d'une bien meilleur application de démonstration !"
... mais il y a dû y avoir une erreur de manipulation car on y retrouve seulement le sous-dossier nommé 'Demo old' et vérification faite on n'y trouve aucune des fonctions précitées.

En tous cas : Chapeau pour if IsHugPrime(M22,0.99) qui donne "M22 est Premier, Mis : 56 s 718,904 ms" : ça décoiffe !!!(Intel Core i7-2700K à 3,5 GHz.

A+.
Messages postés
324
Date d'inscription
samedi 18 décembre 2004
Statut
Membre
Dernière intervention
30 octobre 2020
2
J'ai une version de mon logiciel OPIX où le test de primalité est effectué avec la bibliothèque GMP connue pour être optimisée. Pour le 22-ième nombre premier de Mersenne : M22 = 2^9941 - 1 qui a 2993 chiffres décimaux le logiciel Demo.exe ici présent de Kamel78 et GMP dans mon logiciel mettent tous les deux 2 minutes et une douzaine de secondes pour vérifier la primalité de M22 sur mon ordinateur qui est un peu vieux. C'est un bien beau résultat.
Afficher les 26 commentaires

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.