nint_1.pas

cipher 5.1。一个几乎包含了所有常见的加密算法的控件

{Copyright:      Hagen Reddmann  HaReddmann at T-Online dot de
 Author:         Hagen Reddmann
 Remarks:        this Copyright must be included
 known Problems: none
 Version:        5.1,  Part II from Delphi Encryption Compendium
                 Delphi 5

 * THIS SOFTWARE IS PROVIDED BY THE AUTHORS ''AS IS'' AND ANY EXPRESS
 * OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
 * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE
 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
 * WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
 * OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE,
 * EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

}
unit NInt_1;

interface

uses SysUtils, NInts, Console, NMath, ConsoleForm;

type
  TCFEFunc = function(Index: Integer): Integer; register;
  TIIntegerPIMethod = (piFastChudnovsky, piIterChudnovsky, piAGM, piFastMachin, piIterMachin);

procedure NCFE(var P,Q: IInteger; CFEFunc: TCFEFunc; Last: Integer);

function  CFE_Euler(Index: Integer): Integer;
function  CFE_GoldenRatio(Index: Integer): Integer;
function  CFE_Tan1(Index: Integer): Integer;

procedure NLn2(var R: IInteger);
procedure NLn10(var R: IInteger);
procedure NArcTan(var R: IInteger; const U,V: IInteger); overload;
procedure NArcTan(var R: IInteger; V: Integer); overload;
procedure NArcTanh(var R: IInteger; const V: IInteger);
procedure NSin(var R: IInteger; const U,V: IInteger);
procedure NSinh(var R: IInteger; const U,V: IInteger);
procedure NCos(var R: IInteger; const U,V: IInteger);
procedure NCosh(var R: IInteger; const U,V: IInteger);
procedure NTan(var R: IInteger; const U,V: IInteger);
procedure NTanh(var R: IInteger; const U,V: IInteger);
procedure NExp(var A: IInteger; U: Integer = 1; V: Integer = 1); overload;
procedure NExp(var A: IInteger; const U,V: IInteger); overload;
function  NPi(var A: IInteger; Decimals: Cardinal; Method: TIIntegerPIMethod = piFastChudnovsky): Cardinal;

procedure NFactorial1(var A: IInteger; N: Cardinal);

implementation

procedure NCFE(var P,Q: IInteger; CFEFunc: TCFEFunc; Last: Integer);
// continued fraction expansion
var
  I: Integer;
  T: IInteger;
begin
  NSet(P, 1);
  NSet(Q, 0);
  for I := Last downto 0 do
  begin
    NMul(T, P, CFEFunc(I));
    NAdd(Q, T);
    NSwp(Q, P);
  end;
end;

function CFE_Euler(Index: Integer): Integer;
// can be checked via NExp(.., 1, 1)
begin
  if Index = 0 then Result := 2 else
    if Index mod 3 <> 2 then Result := 1
      else Result := ((Index + 1) div 3) shl 1;
end;

function CFE_GoldenRatio(Index: Integer): Integer;
begin
  Result := 1;
end;

function CFE_Tan1(Index: Integer): Integer;
begin
  if Odd(Index) then Result := Index
    else Result := 1;
end;
{
procedure NLn2(var R: IInteger);
// R = R * Ln(2), = 2 * arctanh(1/3)
var
  S,T: IInteger;
  K: Cardinal;
begin
  K := 1;
  NShl(S, R, 33);
  NDiv(S, 3);
  NSet(R, 0);
  while NSgn(S) <> 0 do
  begin
    NDiv(T, S, K);
    NAdd(R, T);
    NDiv(S, 9);
    Inc(K, 2);
  end;
  NShr(R, 32);
end; }

procedure NLn2(var R: IInteger);
// R = R * Ln(2)
//   = R * 144 * arctanh(1/251) + 54 * arctanh(1/449) - 38 * arctanh(1/4801) + 62 * arctanh(1/8749)
var
  A,B: IInteger;
begin
  NShl(R, 32); // 1 digit rounding
  NMul(A, R, 144); NArcTanh(A, NInt( 251));
  NMul(B, R,  54); NArcTanh(B, NInt( 449)); NAdd(A, B);
  NMul(B, R,  38); NArcTanh(B, NInt(4801)); NSub(A, B);
  NMul(B, R,  62); NArcTanh(B, NInt(8749)); NAdd(A, B);
  NShr(R, A, 32);
end;

procedure NLn10(var R: IInteger);
// R = R * Ln(10)
//   = R * 478 * arctanh(1/251) + 180 * arctanh(1/449) - 126 * arctanh(1/4801) + 206 * arctanh(1/8749)
var
  A,B: IInteger;
begin
  NShl(R, 32); // 1 digit rounding
  NMul(A, R, 478); NArcTanh(A, NInt( 251));
  NMul(B, R, 180); NArcTanh(B, NInt( 449)); NAdd(A, B);
  NMul(B, R, 126); NArcTanh(B, NInt(4801)); NSub(A, B);
  NMul(B, R, 206); NArcTanh(B, NInt(8749)); NAdd(A, B);
  NShr(R, A, 32);
end;

procedure NArcTanh(var R: IInteger; const V: IInteger);
// R = R * arctanh(1 / V)
// binary splitting adaption of J鰎g Arndt's formulas

// method 1: arctanh(1 / V) = sum[n=0..inf] 1 / (2n +1) * 1 / V^(2n +1)
// method 2: arctanh(1 / V) = sum[n=0..inf] (-4)^n*n!^2 / (2n +1)! * V / (V^2 -1)^(n +1)
var
  S: IInteger;

  procedure DoAtanh1(N: Integer; var D: TIIntegerSplitData);
  // 1. method, S = V^2
  begin
    NSet(D.B, N + N +1);
    if N > 0 then D.Q := S       // we assign here instead of use of NSet(Q, S)
      else D.Q := V;             // this save memory and speedup by "copy on demand"
  end;

  procedure DoAtanh2(N: Integer; var D: TIIntegerSplitData);
  // 2. method, S = V^2 -1
  begin
    if N > 0 then
    begin
      NSet(D.P, -(N + N));
      NMul(D.Q, S, N + N +1);
    end else
    begin
      D.P := V;
      D.Q := S;
    end;
  end;

var
  L: Integer;
  P,Q: IInteger;
begin
  if NSgn(R) <> 0 then
  begin
    NSqr(S, V);
    NDec(S, 1); // comment out for method 1
    Assert(NSgn(S) > 0, 'NArcTanh(), require V > 1');
    L := Round(NLn(R) / NLn(S)) +1;
    NBinarySplitting(P, Q, L, @DoAtanh2);
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

(*
procedure NArcTan(var R: IInteger; N: Integer);
// R = R * ArcTan(1 / N)
// iterative version
var
  K,V: Integer;
  S,T: IInteger;
begin
  K := 0;
  V := N * N;
  NDiv(S, R, N);
  NSet(R, 0);
  while NSgn(S) <> 0 do
  begin
    NDiv(T, S, 2 * K + 1);
    NNeg(T, Odd(K));
    NAdd(R, T);
    NDiv(S, V);
    Inc(K);
  end;
end; *)

procedure NArcTan(var R: IInteger; const U,V: IInteger); overload;
// R = R * arctan(U / V)
// binary splitting adaption of J鰎g Arndt's formulas

// method 1: arctan(U / V) = sum[n=0..inf] (-1)^n / (2n +1) * U / V^(2n +1)
// method 2: arctan(1 / V) = sum[n=0..inf] 4^n * n!^2 / (2n +1)! * V / (V^2 +1)^(n +1)

var
  sV,sU: IInteger;

  procedure DoAtan1(N: Integer; var D: TIIntegerSplitData);
  // 1. method, S = V^2
  begin
    NSet(D.B, N + N +1);
    NNeg(D.B, Odd(N));
    if N > 0 then
    begin
      D.P := sU;
      D.Q := sV;
    end else
    begin
      D.P := U;
      D.Q := V;
    end;
  end;

  procedure DoAtan2(N: Integer; var D: TIIntegerSplitData);
  // 2. method, S = V^2 +1
  begin
    if N > 0 then
    begin
      NSet(D.P, N + N);
      NMul(D.Q, sV, N + N +1);
    end else
    begin
      D.P := V;
      D.Q := sV;
    end;
  end;

resourcestring
  sNArcTan = 'NArcTan(), require U <> 0 and |U| < |V|';
var
  P,Q: IInteger;
  L: Extended;
begin
  if (NSgn(U) = 0) or (NCmp(U, V, True) >= 0) then NRaise(@sNArcTan);
  if NSgn(R) <> 0 then
  begin
    NSqr(sV, V);
    if Abs(NSgn(U, True)) <> 1 then
    begin
      NSqr(sU, U);
      L := NSize(R) / Abs(NLn(sU) - NLn(sV)) +1;
      NBinarySplitting(P, Q, Round(L), @DoAtan1);
    end else
    begin
      NInc(sV);
      L := NSize(R) / NLn(sV) +1;
      NBinarySplitting(P, Q, Round(L), @DoAtan2);
    end;
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

procedure NArcTan(var R: IInteger; V: Integer); overload;
// R = R * arctan(1 / V)
// binary splitting adaption of J鰎g Arndt's formulas

// method: arctan(1 / V) = sum[n=0..inf] 4^n * n!^2 / (2n +1)! * V / (V^2 +1)^(n +1)
var
  S: Int64;

  procedure DoAtan(N: Integer; var D: TIIntegerSplitData);
  begin
    if N > 0 then
    begin
      Inc(N, N);
      NSet(D.P, N);
      Inc(N, 1);
      NSet(D.Q, S);
      NMul(D.Q, N);
    end else
    begin
      NSet(D.P, V);
      NSet(D.Q, S);
    end;
  end;

resourcestring
  sNArcTan = 'NArcTan(), require V <> 0';
var
  P,Q: IInteger;
  L: Extended;
begin
  if V = 0 then NRaise(@sNArcTan);
  if NSgn(R) <> 0 then
  begin
    S := Int64(V) * Int64(V) + 1; // avoid D4 Bug for S := V * V +1 !!
    Assert(S > 0);
    L := S;
    L := NLn(R) / Ln(L) +1;
    NBinarySplitting(P, Q, Round(L), @DoAtan);
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

procedure NSin(var R: IInteger; const U,V: IInteger);
// R = R * sin(U / V)
var
  sU,sV: IInteger;

  procedure DoSIN(N: Integer; var D: TIIntegerSplitData); register;
  begin
    if N > 0 then
    begin
      D.P := sU;                         // -U^2
      NSet(D.Q, Int64(N * (N + N +1)));  // 2V^2 * (2n^2 +n)             n(2n +1)
      NMul(D.Q, sV);
    end else
    begin
      NSet(D.P, U);                    // (-1)^n * x^(2n + 1) / (2n + 1)!
      NSet(D.Q, V);
    end;
  end;

resourcestring
  sNSin = 'NSin(), requiere V <> 0 and U <> 0';
var
  P,Q: IInteger;
  C,D,L: Extended;
begin
  if (NSgn(V) = 0) or (NSgn(U) = 0) then NRaise(@sNSin);
  if NSgn(R) <> 0 then
  begin
    NSqr(sU, U);
    NSqr(sV, V);
    NMul(sV, 2);
  // compute series length
    D := NLn(R);
    C := NLn(sU) - NLn(sV);
    L := 1;
    while D > 0 do
    begin
      D := D + C - Ln(L * (L + L +1));
      L := L + 1;
    end;
    NNeg(sU);
    NBinarySplitting(P, Q, Round(L), @DoSIN);  // P / Q = Sin(U / V)
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

procedure NSinh(var R: IInteger; const U,V: IInteger);
// R = R * sinh(U / V)
var
  sU,sV: IInteger;

  procedure DoSINH(N: Integer; var D: TIIntegerSplitData); register;
  begin
    if N > 0 then
    begin
      D.P := sU;                         // -U^2
      NSet(D.Q, Int64(N * (N + N +1)));  // 2V^2 * (2n^2 +n)
      NMul(D.Q, sV);
    end else
    begin
      NSet(D.P, U);                      //  U
      NSet(D.Q, V);                      //  V
    end;
  end;

resourcestring
  sNSinh = 'NSinh(), requiere V <> 0 and U <> 0';
var
  P,Q: IInteger;
  C,D,L: Extended;
begin
  if (NSgn(V) = 0) or (NSgn(U) = 0) then NRaise(@sNSinh);
  if NSgn(R) <> 0 then
  begin
    NSqr(sU, U);
    NSqr(sV, V);
    NMul(sV, 2);
  // compute series length
    D := NLn(R);
    C := NLn(sU) - NLn(sV);
    L := 1;
    while D > 0 do
    begin
      D := D + C - Ln(L * (L + L +1));
      L := L + 1;
    end;
    NBinarySplitting(P, Q, Round(L), @DoSINH);  // P / Q = Sin(U / V)
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

procedure NCos(var R: IInteger; const U,V: IInteger);
// R = R * cos(U / V)
var
  sU,sV: IInteger;

  procedure DoCOS(N: Integer; var D: TIIntegerSplitData); register;
  begin
    if N > 0 then
    begin
      D.P := sU;                         // -U^2
      NSet(D.Q, Int64(N * (N + N - 1))); //  2V^2 * (2n^2 -n)
      NMul(D.Q, sV);
    end;
  end;

resourcestring
  sNCos = 'NCos(), requiere V <> 0 and U <> 0';
var
  P,Q: IInteger;
  C,D,L: Extended;
begin
  if (NSgn(V) = 0) or (NSgn(U) = 0) then NRaise(@sNCos);
  if NSgn(R) <> 0 then
  begin
    NSqr(sU, U);
    NSqr(sV, V);
    NMul(sV, 2);
  // compute series length
    D := NLn(R);
    C := NLn(sU) - NLn(sV);
    L := 1;
    while D > 0 do
    begin
      D := D + C - Ln(L * (L + L +1));
      L := L + 1;
    end;
    NNeg(sU);
    NBinarySplitting(P, Q, Round(L), @DoCOS);  // P / Q = Sin(U / V)
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

procedure NCosh(var R: IInteger; const U,V: IInteger);
// R = R * cosh(U / V)
var
  sU,sV: IInteger;

  procedure DoCOSH(N: Integer; var D: TIIntegerSplitData); register;
  begin
    if N > 0 then
    begin
      D.P := sU;                         // -U^2
      NSet(D.Q, Int64(N * (N + N - 1))); //  2V^2 * (2n^2 -n)
      NMul(D.Q, sV);
    end;
  end;

resourcestring
  sNCosh = 'NCosh(), requiere V <> 0 and U <> 0';
var
  P,Q: IInteger;
  C,D,L: Extended;
begin
  if (NSgn(V) = 0) or (NSgn(U) = 0) then NRaise(@sNCosh);
  if NSgn(R) <> 0 then
  begin
    NSqr(sU, U);
    NSqr(sV, V);
    NMul(sV, 2);
  // compute series length
    D := NLn(R);
    C := NLn(U) - NLn(V);
    L := 1;
    while D > 0 do
    begin
      D := D + C - Ln(L * (L + L +1));
      L := L + 1;
    end;
    NBinarySplitting(P, Q, Round(L), @DoCOSH);  // P / Q = Sin(U / V)
    NMul(R, P);
    NDiv(R, Q);
  end;
end;

procedure NTan(var R: IInteger; const U,V: IInteger);
// R = R * tan(U / V)
// todo: implement a direct binary splitting
// computation of NTan(R, 1, 1) can be checked with NCFE(.., .., CFETan1, ...);
var
  A: IInteger;