From a034fbaedcd4941106abe5b44b3e57ba930dea68 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Stefan=20M=C3=BCller?= <stefan.mueller@aksdb.de>
Date: Wed, 31 Jan 2024 18:58:49 +0100
Subject: [PATCH] Added integer factorization and enumeration of dividers

---
 UNumberTheory.pas | 220 +++++++++++++++++++++++++++++++++++++++++++++-
 1 file changed, 219 insertions(+), 1 deletion(-)

diff --git a/UNumberTheory.pas b/UNumberTheory.pas
index 99b7fd5..a247344 100644
--- a/UNumberTheory.pas
+++ b/UNumberTheory.pas
@@ -22,7 +22,7 @@ unit UNumberTheory;
 interface
 
 uses
-  Classes, SysUtils;
+  Classes, SysUtils, Generics.Collections, Math;
 
 type
 
@@ -34,6 +34,52 @@ type
     class function LeastCommonMultiple(AValue1, AValue2: Int64): Int64;
   end;
 
+  TInt64Array = array of Int64;
+
+  { TIntegerFactor }
+
+  TIntegerFactor = record
+    Factor: Int64;
+    Exponent: Byte;
+  end;
+
+  TIntegerFactors = specialize TList<TIntegerFactor>;
+
+  { TIntegerFactorization }
+
+  TIntegerFactorization = class
+  public
+    class function PollardsRhoAlgorithm(const AValue: Int64): TInt64Array;
+    class function GetNormalized(constref AIntegerFactorArray: TInt64Array): TIntegerFactors;
+  end;
+
+  { TDividersEnumerator }
+
+  TDividersEnumerator = class
+  private
+    FFactors: TIntegerFactors;
+    FCurrentExponents: array of Byte;
+    function GetCount: Integer;
+  public
+    constructor Create(constref AIntegerFactorArray: TInt64Array);
+    destructor Destroy; override;
+    function GetCurrent: Int64;
+    function MoveNext: Boolean;
+    procedure Reset;
+    property Current: Int64 read GetCurrent;
+    property Count: Integer read GetCount;
+  end;
+
+  { TDividers }
+
+  TDividers = class
+  private
+    FFactorArray: TInt64Array;
+  public
+    constructor Create(constref AIntegerFactorArray: TInt64Array);
+    function GetEnumerator: TDividersEnumerator;
+  end;
+
 implementation
 
 { TNumberTheory }
@@ -58,5 +104,177 @@ begin
   Result := (Abs(AValue1) div GreatestCommonDivisor(AValue1, AValue2)) * Abs(AValue2);
 end;
 
+{ TIntegerFactorization }
+
+// https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
+class function TIntegerFactorization.PollardsRhoAlgorithm(const AValue: Int64): TInt64Array;
+var
+  primes: specialize TList<Int64>;
+  composites: specialize TStack<Int64>;
+  factor, n: Int64;
+  i: Integer;
+
+  function G(const AX, AC: Int64): Int64;
+  begin
+    Result := (AX * AX + AC) mod n;
+  end;
+
+  function FindFactor(const AStartValue, AC: Int64): Int64;
+  var
+    x, y, d: Int64;
+  begin
+    x := AStartValue;
+    y := x;
+    d := 1;
+    while d = 1 do
+    begin
+      x := G(x, AC);
+      y := G(G(y, AC), AC);
+      d := TNumberTheory.GreatestCommonDivisor(Abs(x - y), n);
+    end;
+    Result := d;
+  end;
+
+begin
+  primes := specialize TList<Int64>.Create;
+  composites := specialize TStack<Int64>.Create;
+
+  n := Abs(AValue);
+  while (n and 1) = 0 do
+  begin
+    primes.Add(2);
+    n := n shr 1;
+  end;
+
+  composites.Push(n);
+  while composites.Count > 0 do
+  begin
+    n := composites.Pop;
+    i := 0;
+    repeat
+      factor := FindFactor(2 + (i + 1) div 2, 1 - i div 2);
+      if factor < n then
+      begin
+        composites.Push(factor);
+        composites.Push(n div factor);
+      end;
+      Inc(i);
+    until (factor < n) or (i > 3);
+    if factor = n then
+      primes.Add(factor);
+  end;
+
+  Result := primes.ToArray;
+
+  primes.Free;
+  composites.Free;
+end;
+
+class function TIntegerFactorization.GetNormalized(constref AIntegerFactorArray: TInt64Array): TIntegerFactors;
+var
+  i: Integer;
+  factor: Int64;
+  normal: TIntegerFactor;
+  found: Boolean;
+begin
+  Result := TIntegerFactors.Create;
+  for factor in AIntegerFactorArray do
+  begin
+    found := False;
+    for i := 0 to Result.Count - 1 do
+      if Result[i].Factor = factor then
+      begin
+        found := True;
+        normal := Result[i];
+        Inc(normal.Exponent);
+        Result[i] := normal;
+        Break;
+      end;
+    if not found then
+    begin
+      normal.Factor := factor;
+      normal.Exponent := 1;
+      Result.Add(normal);
+    end;
+  end;
+end;
+
+{ TDividersEnumerator }
+
+function TDividersEnumerator.GetCount: Integer;
+var
+  factor: TIntegerFactor;
+begin
+  if FFactors.Count > 0 then
+  begin
+    Result := 1;
+    for factor in FFactors do
+      Result := Result * factor.Exponent;
+    Dec(Result);
+  end
+  else
+    Result := 0;
+end;
+
+constructor TDividersEnumerator.Create(constref AIntegerFactorArray: TInt64Array);
+begin
+  FFactors := TIntegerFactorization.GetNormalized(AIntegerFactorArray);
+  SetLength(FCurrentExponents, FFactors.Count);
+end;
+
+destructor TDividersEnumerator.Destroy;
+begin
+  FFactors.Free;
+end;
+
+function TDividersEnumerator.GetCurrent: Int64;
+var
+  i: Integer;
+begin
+  Result := 1;
+  for i := Low(FCurrentExponents) to High(FCurrentExponents) do
+    if FCurrentExponents[i] > 0 then
+      Result := Result * Round(Power(FFactors[i].Factor, FCurrentExponents[i]));
+end;
+
+function TDividersEnumerator.MoveNext: Boolean;
+var
+  i: Integer;
+begin
+  Result := False;
+  i := 0;
+  while (i <= High(FCurrentExponents)) and (FCurrentExponents[i] >= FFactors[i].Exponent) do
+  begin
+    FCurrentExponents[i] := 0;
+    Inc(i);
+  end;
+
+  if i <= High(FCurrentExponents) then
+  begin
+    Inc(FCurrentExponents[i]);
+    Result := True;
+  end;
+end;
+
+procedure TDividersEnumerator.Reset;
+var
+  i: Integer;
+begin
+  for i := Low(FCurrentExponents) to High(FCurrentExponents) do
+    FCurrentExponents[i] := 0;
+end;
+
+{ TDividers }
+
+constructor TDividers.Create(constref AIntegerFactorArray: TInt64Array);
+begin
+  FFactorArray := AIntegerFactorArray;
+end;
+
+function TDividers.GetEnumerator: TDividersEnumerator;
+begin
+  Result := TDividersEnumerator.Create(FFactorArray);
+end;
+
 end.