ncombi.pas

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

{Copyright:      Hagen Reddmann  HaReddmann at T-Online dot de
 Author:         Hagen Reddmann
 Remarks:        public domain
 known Problems: none
 Version:        5.1, Delphi 5, 9.3.2002
 Description:    different Algorithms to compute the Factorial N! to arbitary value
                 and other usefull combinatoric algorithm

 Comments:       Thanks to Peter Luschny
                 Some factorial functions are based on the JAVA Sources at
                 http://www.luschny.de/math/factorial/FastFactorialFunctions.htm

 * 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 NCombi;

interface

uses NInts;

procedure NFactorial_Moessner(var A: IInteger; N: Cardinal);
procedure NFactorial_Naive(var A: IInteger; N: Cardinal);
procedure NFactorial_Recursive(var A: IInteger; N: Cardinal);
procedure NFactorial_DivideAndConquer(var A: IInteger; N: Cardinal);
procedure NFactorial_Binomial(var A: IInteger; N: Cardinal);
procedure NFactorial_Jason_GMP(var A: IInteger; N: Cardinal);

procedure NFactorial(var A: IInteger; N: Cardinal; const M: IInteger = nil); overload;
procedure NComporial(var A: IInteger; N: Cardinal; const M: IInteger = nil); overload;
procedure NBinomial(var A: IInteger; N,K: Cardinal; const M: IInteger = nil); overload;
procedure NProduct(var A: IInteger; N,K: Cardinal; const M: IInteger = nil); overload;
procedure NPermutation(var A: IInteger; N,K: Cardinal; const M: IInteger = nil); overload;
function  NOddFactorial(var A: IInteger; N: Cardinal; const M: IInteger = nil): Cardinal; overload;
procedure NPrimorial(var A: IInteger; K,N: Cardinal; const M: IInteger = nil); overload;

procedure NHalfFactorial(var A: IInteger; N: Cardinal; const M: IInteger = nil); overload;
procedure NHalfComporial(var A: IInteger; N: Cardinal; const M: IInteger = nil); overload;
procedure NHalfPrimorial(var A: IInteger; N: Cardinal; const M: IInteger = nil); overload;
procedure NHalfBinomial(var A: IInteger; N: Cardinal; const M: IInteger = nil); overload;

function  NFactorialTrailingZeros(const N: Int64; Base: Cardinal): Int64;

procedure NTestCombi(N: Cardinal = 200); overload;

// internal most used stuff
type
  TPowerTable = array of record
    B: Cardinal; // base
    E: Cardinal; // exponent
  end;

function NPowerTable(var T: TPowerTable; N: Cardinal; L: Cardinal = 0; K: Cardinal = 0): Cardinal; overload;
function NPowerTable(var A: IInteger; N: Cardinal; L: Cardinal = 0; K: Cardinal = 0;
                       Shift: Boolean = False; const M: IInteger = nil): Cardinal; overload;
function NPrd(var A: IInteger; const T: TPowerTable; E: Cardinal; const M: IInteger = nil): Boolean; overload;
function NPrd(var A: IInteger; const T: TPowerTable; const M: IInteger = nil): Boolean; overload;

implementation

uses SysUtils, Prime, NMath, Math;

// follow four different Factorial procedures, based on Peter Luschny JAVA source
procedure NFactorial_Moessner(var A: IInteger; N: Cardinal);
var
  S: array of IInteger;
  M,K,I: Cardinal;
begin
  SetLength(S, N +1);
  NSet(S[0], 1);
  for M := 1 to N do
    for K := M downto 1 do
      for I := 1 to K do
         NAdd(S[I], S[I -1]);
  NSwp(A, S[N]);
end;

procedure NFactorial_Naive(var A: IInteger; N: Cardinal);
var
  I: Cardinal;
begin
  NSet(A, 1);
  for I := 2 to N do
    NMul(A, I);
end;

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

  procedure NSplit(var P: IInteger; var K: Cardinal; N: Cardinal);
  var
    M: Cardinal;
    Q: IInteger;
  begin
    M := N shr 1;
    if M = 0 then
    begin
      NSet(P, K);
      Inc(K);
    end else
      if N = 2 then
      begin
        NSet(P, K);
        NMul(P, K +1);
        Inc(K, 2);
      end else
      begin
        NSplit(P, K, M);
        NSplit(Q, K, N - M);
        NMul(P, Q);
      end;
  end;

var
  K: Cardinal;
begin
  K := 1;
  if N < 2 then NSet(A, 1)
    else NSplit(A, K, N);
end;

procedure NFactorial_DivideAndConquer(var A: IInteger; N: Cardinal);
var
  BreakEven: Cardinal;
  
  procedure NSplit(var P: IInteger; var K: Cardinal; N: Cardinal);
  var
    M: Cardinal;
    Q: IInteger;
  begin
    if N < BreakEven then
    begin
      NSet(P, K);
      Inc(K, 2);
      while N > 1 do
      begin
        NMul(P, K);
        Inc(K, 2);
        Dec(N);
      end;
    end else
    begin
      M := N shr 1;
      NSplit(P, K, M);
      NSplit(Q, K, N - M);
      NMul(P, Q);
    end;
  end;

var
  I,J,K: Cardinal;
  L: Integer;
  Len: array of Cardinal;
  P,Q: IInteger;
begin
  if N >= 2 then
  begin
    BreakEven := NBreakEven(1);   // Karatsuba Multiplication BreakEven
    J := NLog2(N);
    SetLength(Len, J);
    Dec(J);
    K := N;
    for I := J downto 0 do
    begin
      L := K + K and 1 -1;
      K := K shr 1;
      Dec(L, K + K and 1 +1);
      if L >= 0 then
        Len[I] := 1 + L shr 1;
    end;
    K := 3;
    NSet(A, 1);
    NSet(Q, 1);
    for I := 0 to J do
      if Len[I] > 0 then
      begin
        NSplit(P, K, Len[I]);
        NMul(Q, P);
        NMul(A, Q);
      end;
    NShl(A, N - Cardinal(NBitWeight(N)));
  end else NSet(A, 1);
end;

procedure NFactorial_Binomial(var A: IInteger; N: Cardinal);
// my Method, this is equal fast as plain Sch鰊hage and are the same group of algos.
var
  I: Integer;
  K: Cardinal;
  F,B: IInteger;
begin
  NSet(F, 1);
  I := NLog2(N);
  K := 1;
  while I > 0 do
  begin
    Dec(I);
    Inc(K, K + N shr I and 1);
    NPowerTable(B, K, K div 2, K div 2);  // OddBinomial(K, [K/2])
    NSqr(F);
    NMul(F, B);
  end;
  NShl(A, F, N - Cardinal(NBitWeight(N)));
end;

procedure NFactorial_Jason_GMP(var A: IInteger; N: Cardinal);
// Factorial contributed to GMP MP Project by Jason at http://217.35.81.229/

  procedure NOddProduct(Low, High: Cardinal; var S: IIntegerArray);

    procedure NSmallOddProduct(var A: IInteger; Start, Step, Count: Cardinal);
    // this procedure is very simplified compared to Jason's one,
    // because they runningtime is never important for big factorials
    // of course with my underlaying mul.
    begin
      NSet(A, Start);
      while Count > 1 do
      begin
        Dec(Count);
        Inc(Start, Step);
        NMul(A, Start);
      end;
    end;

  var
    A,Z,Y,N,Mask,Stc,Stn: Cardinal;
  begin
    Inc(Low);
    if not Odd(Low) then Inc(Low);
    if High = 0 then Inc(High) else
      if not Odd(High) then Dec(High);
    if High < Low then
    begin
      NSet(S[0], 1);
      Exit;
    end;
    if High = Low then
    begin
      NSet(S[0], Low);
      Exit;
    end;
    High := (High - Low) div 2 +1;
    if High <= 1 shl 5 then
    begin
      NSmallOddProduct(S[0], Low, 2, High);
      Exit;
    end;
    N := NLog2(High) - 5;
    Mask := 1 shl N;
    A := Mask + Mask;
    Dec(Mask);
    stn := 0;
    stc := 1;
    for Z := Mask downto 0 do
    begin
      Y := NBitSwap(Z) shr (32 - N);
      NSmallOddProduct(S[stn], Low + 2 * ((not Y) and Mask), A, (High + Y) shr N);
      Inc(Stn);
      Y := Stc;
      Inc(Stc);
      while not Odd(Y) do
      begin
        NMul(S[stn -2], S[stn -1]);
        Dec(stn);
        Y := Y shr 1;
      end;
    end;
  end;

var
  B: IInteger;
  S: IIntegerArray;
  I,J,Z: Integer;
begin
  SetLength(S, 33);
  Z := NLog2(N div 3) +1;
  NSet(A, 1);

{$IFDEF Old_Jason}
// in fact, because of the clever used asymmetric splitting methods in all
// my multiplication stuff like COMBA, Karatsuba, TOOM-3 and Sch鰊hage FFT
// this old used method is equal fast as below code
  NSet(B, 1);
  for I := Z downto 1 do
  begin
    NOddProduct(N shr I, N shr (I -1), S);
    NMul(B, S[0]);
    NMul(A, B);
  end;
{$ELSE}

// Jason's proposed way with squaring by binary powering
  J := 16;
  while J > 0 do
  begin
    NSet(B, 1);
    I := 32 - J;
    while I >= J do
    begin
      if Z >= I then
      begin
        NOddProduct(N shr I, N shr (I -1), S);
        if I <> J then
          NPow(S[0], S[0], I div J);
        NMul(B, S[0]);
      end;
      Dec(I, 2 * J);
    end;
    if (Z >= J) and (J <> 1) then
    begin
      NMul(A, B);
      NSqr(A);
    end;
    J := J shr 1;
  end;
  NMul(A, B);
{$ENDIF}

  NShl(A, N - Cardinal(NBitWeight(N)));
end;

// follow four procedure are the essential of Sch鰊hages Factorial trick

function NPowerTable(var T: TPowerTable; N: Cardinal; L: Cardinal; K: Cardinal): Cardinal;
// compute primepowers in T[] where ((N! / L! / K!) / 2^Result)

  function FindEP(K,L,N,P: Cardinal): Cardinal;
  // extract exponents to base P, where P is prime
  var
    E: Cardinal;
  begin
    E := 0;
    while K >= P do
    begin
      N := N div P;
      L := L div P;
      K := K div P;
      Inc(E, N - L - K);
    end;
    while L >= P do
    begin
      N := N div P;
      L := L div P;
      Inc(E, N - L);
    end;
    while N >= P do
    begin
      N := N div P;
      Inc(E, N);
    end;
    Result := E;
  end;

resourcestring
  sNPowerTable = 'NPowerTable(), requiere K <= N and L <= N and L + K <= N';
var
  P: Cardinal;
  E,I,C: Cardinal;
  S: TSmallPrimeSieve;
begin
  if K > L then
  begin
    I := K;
    K := L;
    L := I;
  end;
  if (K > N) or (L > N) or (L + K > N) then NRaise(@sNPowerTable);
// extract power to base 2
  Result := (N - L - K) - Cardinal(NBitWeight(N) - NBitWeight(L) - NBitWeight(K));
// extract powers to each prime <= N
  S := Primes;
  NAutoRelease(S);
  C := 0;
  I := 1;  // start with prime 3
  repeat
    Inc(I);
    P := S[I];
    if P > N then Break;
    E := FindEP(K, L, N, P);
    if E <> 0 then
    begin
      if C mod 2048 = 0 then SetLength(T, C + 2048);
      T[C].E := E;
      T[C].B := P;
      Inc(C);
    end;
  until False;
  SetLength(T, C);
end;

function NPowerTable(var A: IInteger; N: Cardinal; L: Cardinal; K: Cardinal; Shift: Boolean; const M: IInteger): Cardinal;
// A = ((N! / L! / K!) / 2^Result) mod M if M <> nil
// "divison" by 2^result only if shift = false
var
  T: TPowerTable;
  R: IInteger;
begin
  Result := NPowerTable(T, N, L, K);