AdventOfCode2023/UNumberTheory.pas

{
  Solutions to the Advent Of Code.
  Copyright (C) 2023-2024  Stefan Müller

  This program is free software: you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation, either version 3 of the License, or (at your option) any later
  version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.

  You should have received a copy of the GNU General Public License along with
  this program.  If not, see <http://www.gnu.org/licenses/>.
}

unit UNumberTheory;

{$mode ObjFPC}{$H+}

interface

uses
  Classes, SysUtils, Generics.Collections, Math;

type

  { TNumberTheory }

  TNumberTheory = class
  public
    class function GreatestCommonDivisor(AValue1, AValue2: Int64): Int64;
    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 }

class function TNumberTheory.GreatestCommonDivisor(AValue1, AValue2: Int64): Int64;
var
  t: Int64;
begin
  AValue1 := Abs(AValue1);
  AValue2 := Abs(AValue2);
  while AValue2 <> 0 do
  begin
    t := AValue2;
    AValue2 := AValue1 mod AValue2;
    AValue1 := t;
  end;
  Result := AValue1;
end;

class function TNumberTheory.LeastCommonMultiple(AValue1, AValue2: Int64): Int64;
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.