AdventOfCode2023/UBigInt.pas

{
  Solutions to the Advent Of Code.
  Copyright (C) 2022-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 UBigInt;

{$mode ObjFPC}{$H+}

interface

uses
  Classes, SysUtils, Math;

type
  TDigits = array of Cardinal;

  { TBigInt }

  // This is an abbreviated reimplementation of a C# class created in 2022.
  TBigInt = object
  private
    FDigits: TDigits;
    FIsNegative: Boolean;
    // Adds A and B, ignoring their signs and using ReturnNegative instead. The result is
    //   Sign * (Abs(A) + Abs(B)),
    // where Sign is 1 for ReturnNegative = False and -1 otherwise.
    class function AddAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt; static;
    // Subtracts B from A, ignoring their signs. However, the result might be negative, and the sign can be reversed by
    // setting ReturnNegative to True. The result is
    //   Sign * (Abs(A) - Abs(B)),
    // where Sign is 1 for ReturnNegative = False and -1 otherwise.
    class function SubtractAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt; static;
    // Multiplies A and B, ignoring their signs and using ReturnNegative instead. This multiplication uses a recursive
    // implementation of the Karatsuba algorithm. See
    // https://www.geeksforgeeks.org/karatsuba-algorithm-for-fast-multiplication-using-divide-and-conquer-algorithm/
    // The result is
    //   Sign * (Abs(a) * Abs(b))
    // where Sign is 1 for ReturnNegative = False and -1 otherwise.
    class function MultiplyAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt; static;
    // Copies consecutive digits from this BigInt to create a new one. The result will be positive. Leading zeros are
    // removed from the result, but AIndex + ACount must not exceed the number of digits of this BigInt.
    // AIndex is the first (least significant) digit to be taken. The digit with this index will become the 0th digit of
    // the new BigInt.
    // ACount is the number of consecutive digits to be taken, and the number of digits of the result.
    function GetSegment(const AIndex, ACount: Integer): TBigInt;
    // Compares the absolute value of this TBigInt object to the absolute value of another one. Returns -1 if this
    // object is less than AOther, 1 if this object is greater than AOther, and 0 if they are equal.
    function CompareToAbsoluteValues(constref AOther: TBigInt): Integer;
  public
    property IsNegative: Boolean read FIsNegative;
    constructor InitZero;
    constructor Init(const AValue: Int64);
    destructor Done;
    function CompareTo(constref AOther: TBigInt): Integer;
  end;

  operator := (const A: Int64): TBigInt;
  operator + (const A, B: TBigInt): TBigInt;
  operator - (const A, B: TBigInt): TBigInt;
  operator * (const A: TBigInt; const B: Int64): TBigInt;
  operator shl (const A: TBigInt; const B: Integer): TBigInt;
  operator = (const A: TBigInt; const B: Int64): Boolean;

implementation

const
  CBase = Cardinal.MaxValue + 1;
  CMaxDigit = Cardinal.MaxValue;
  CDigitSize = SizeOf(Cardinal);
  CBitsPerDigit = CDigitSize * 8;
  CHalfBits = CBitsPerDigit >> 1;
  CHalfDigitMax = (1 << CHalfBits) - 1;

{ TBigInt }

class function TBigInt.AddAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
var
  i, lenA, lenB, len, shorter: Integer;
  carry: Cardinal;
begin
  lenA := Length(AA.FDigits);
  lenB := Length(AB.FDigits);

  // Initializes the digits array of the result, with a simple test to try to predict a carry-over into a new digit. The
  // result could still carry into new digit depending on lower digits (carry over multiple digits), which would be
  // handled at the end.
  if lenA = lenB then
    if CMaxDigit - AA.FDigits[lenA - 1] < AB.FDigits[lenB - 1] then
      // Result will carry into new digit.
      SetLength(Result.FDigits, lenA + 1)
    else
      SetLength(Result.FDigits, lenA)
  else
    SetLength(Result.FDigits, Max(lenA, lenB));
  len := Length(Result.FDigits);

  // Calculates the new digits from less to more significant until the end of the shorter operand is reached.
  shorter := Min(Length(AA.FDigits), Length(AB.FDigits));
  i := 0;
  carry := 0;
  while i < shorter do
  begin
    if (AB.FDigits[i] = CMaxDigit) and (carry > 0) then
    begin
      Result.FDigits[i] := AA.FDigits[i];
      carry := 1;
    end
    else
      if CMaxDigit - AA.FDigits[i] < AB.FDigits[i] + carry then
      begin
        Result.FDigits[i] := AB.FDigits[i] + carry - 1 - (CMaxDigit - AA.FDigits[i]);
        carry := 1;
      end
      else begin
        Result.FDigits[i] := AA.FDigits[i] + AB.FDigits[i] + carry;
        carry := 0;
      end;
    Inc(i);
  end;

  // Copies the missing unchanged digits from the longer operand to the result, if any, before applying remaining
  // carry-overs. This avoids additional tests for finding the shorter digit array.
  if (i < lenA) or (i < lenB) then
    if lenA >= lenB then
      Move(AA.FDigits[i], Result.FDigits[i], CDigitSize * (len - i))
    else
      Move(AB.FDigits[i], Result.FDigits[i], CDigitSize * (len - i));

  // Applies the remaining carry-overs until the end of the prepared result digit array.
  while (carry > 0) and (i < len) do
  begin
    if Result.FDigits[i] = CMaxDigit then
      Result.FDigits[i] := 0
    else begin
      Inc(Result.FDigits[i]);
      carry := 0;
    end;
    Inc(i);
  end;

  // Applies the carry-over into a new digit that was not anticipated in the initialization at the top (carry over
  // multiple digits).
  if carry > 0 then
    Insert(1, Result.FDigits, 0);
  Result.FIsNegative := AReturnNegative;
end;

class function TBigInt.SubtractAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
var
  a, b: TBigInt;
  carry: Cardinal;
  i, lastNonZeroDigitIndex, len: Integer;
begin
  // Establishes the operand order, such that Abs(a) is not less than Abs(b).
  if (AA.CompareToAbsoluteValues(AB) >= 0) then
  begin
    a := AA;
    b := AB;
    Result.FIsNegative := AReturnNegative;
  end
  else begin
    a := AB;
    b := AA;
    Result.FIsNegative := not AReturnNegative;
  end;

  // Calculates the new digits from less to more significant until the end of the shorter operand is reached and all
  // carry-overs have been applied.
  len := Length(a.FDigits);
  SetLength(Result.FDigits, len);
  carry := 0;
  // Tracks leading zeros for the trim below.
  lastNonZeroDigitIndex := 0;
  i := 0;
  while i < Length(b.FDigits) do
  begin
    if (a.FDigits[i] = b.FDigits[i]) and (carry > 0) then
    begin
      Result.FDigits[i] := CMaxDigit;
      carry := 1;
      lastNonZeroDigitIndex := i;
    end
    else begin
      if a.FDigits[i] < b.FDigits[i] then
      begin
        Result.FDigits[i] := CMaxDigit - (b.FDigits[i] - a.FDigits[i]) + 1 - carry;
        carry := 1;
      end
      else begin
        Result.FDigits[i] := a.FDigits[i] - b.FDigits[i] - carry;
        carry := 0;
      end;
      if Result.FDigits[i] > 0 then
        lastNonZeroDigitIndex := i;
    end;
    Inc(i);
  end;
  while carry > 0 do
  begin
    if a.FDigits[i] = 0 then
    begin
      Result.FDigits[i] := CMaxDigit;
      lastNonZeroDigitIndex := i;
    end
    else begin
      Result.FDigits[i] := a.FDigits[i] - carry;
      carry := 0;
      if Result.FDigits[i] > 0 then
        lastNonZeroDigitIndex := i;
    end;
    Inc(i);
  end;

  // Copies the missing unchanged digits from the longer operand to the result, if any. If there are none, then no trim
  // needs to occur because the most significant digit is not zero.
  if i < len then
    Move(a.FDigits[i], Result.FDigits[i], CDigitSize * (len - i))
  else if (lastNonZeroDigitIndex + 1 < len) then
    // Trims leading zeros from the digits array.
    Delete(Result.FDigits, lastNonZeroDigitIndex + 1, len - lastNonZeroDigitIndex - 1);
end;

class function TBigInt.MultiplyAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
var
  lenA, lenB, lenMax, floorHalfLength, ceilHalfLength: Integer;
  a1, a0, b1, b0, a1b1, a0b0: Cardinal;
  am, bm, middle, biga1, biga0, bigb1, bigb0, biga1b1, biga0b0: TBigInt;
begin
  lenA := Length(AA.FDigits);
  lenB := Length(AB.FDigits);
  if (lenA <= 1) and (lenB <= 1) then
  begin
    if (AA.FDigits[0] <= CHalfDigitMax) and (AB.FDigits[0] <= CHalfDigitMax) then
      if (AA.FDigits[0] = 0) or (AB.FDigits[0] = 0) then
        Result.InitZero
      else begin
        Result.FDigits := TDigits.Create(AA.FDigits[0] * AB.FDigits[0]);
        Result.FIsNegative := AReturnNegative;
      end
    else begin
      // a1, a0, b1, b0 use only the lower (less significant) half of the bits of a digit, so the product of any two of
      // these fits in one digit.
      a1 := AA.FDigits[0] >> CHalfBits;
      a0 := AA.FDigits[0] and CHalfDigitMax;
      b1 := AB.FDigits[0] >> CHalfBits;
      b0 := AB.FDigits[0] and CHalfDigitMax;
      a1b1 := a1 * b1;
      a0b0 := a0 * b0;

      if a1b1 > 0 then
        Result.FDigits := TDigits.Create(a0b0, a1b1)
      else
        Result.FDigits := TDigits.Create(a0b0);
      Result.FIsNegative := AReturnNegative;

      // The result of (a1 + a0) * (b1 + b0) might not fit in one digit, so one last recursion step is necessary.
      am.Init(a1 + a0);
      bm.Init(b1 + b0);
      middle := (MultiplyAbsoluteValues(am, bm, False) - a1b1 - a0b0) << CHalfBits;
      if AReturnNegative then
        Result := Result - middle
      else
        Result := Result + middle;
    end;
  end
  else begin
    // Calculates where to split the two numbers.
    lenMax := Max(lenA, lenB);
    floorHalfLength := lenMax >> 1;
    ceilHalfLength := lenMax - floorHalfLength;

    // Performs one recursion step.
    if ceilHalfLength < lenA then
    begin
      biga1 := AA.GetSegment(ceilHalfLength, lenA - ceilHalfLength);
      biga0 := AA.GetSegment(0, ceilHalfLength);

      if ceilHalfLength < lenB then
      begin
        bigb1 := AB.GetSegment(ceilHalfLength, lenB - ceilHalfLength);
        bigb0 := AB.GetSegment(0, ceilHalfLength);
        biga1b1 := MultiplyAbsoluteValues(biga1, bigb1, AReturnNegative);
        biga0b0 := MultiplyAbsoluteValues(biga0, bigb0, AReturnNegative);
        Result := (biga1b1 << (2 * ceilHalfLength * CBitsPerDigit))
          + ((MultiplyAbsoluteValues(biga1 + biga0, bigb1 + bigb0, AReturnNegative) - biga1b1 - biga0b0) << (ceilHalfLength * CBitsPerDigit))
          + biga0b0;
      end
      else begin
        biga0b0 := MultiplyAbsoluteValues(biga0, AB, AReturnNegative);
        Result := ((MultiplyAbsoluteValues(biga1 + biga0, AB, AReturnNegative) - biga0b0) << (ceilHalfLength * CBitsPerDigit)) + biga0b0;
      end;
    end
    else begin
      bigb1 := AB.GetSegment(ceilHalfLength, lenB - ceilHalfLength);
      bigb0 := AB.GetSegment(0, ceilHalfLength);
      biga0b0 := MultiplyAbsoluteValues(AA, bigb0, AReturnNegative);
      Result := ((MultiplyAbsoluteValues(AA, bigb1 + bigb0, AReturnNegative) - biga0b0) << (ceilHalfLength * CBitsPerDigit)) + biga0b0;
    end;
  end;
end;

function TBigInt.GetSegment(const AIndex, ACount: Integer): TBigInt;
var
  trimmedCount: Integer;
begin
  trimmedCount := ACount;
  while (trimmedCount > 1) and (FDigits[AIndex + trimmedCount - 1] = 0) do
    Dec(trimmedCount);
  SetLength(Result.FDigits, trimmedCount);
  Move(FDigits[AIndex], Result.FDigits[0], CDigitSize * trimmedCount);
  Result.FIsNegative := False;
end;

function TBigInt.CompareToAbsoluteValues(constref AOther: TBigInt): Integer;
var
  i: Integer;
begin
  if Length(FDigits) < Length(AOther.FDigits) then
    Result := -1
  else if Length(FDigits) > Length(AOther.FDigits) then
    Result := 1
  else begin
    Result := 0;
    for i := High(FDigits) downto 0 do
      if FDigits[i] < AOther.FDigits[i] then
      begin
       Result := -1;
       Break;
      end
      else if FDigits[i] > AOther.FDigits[i] then
      begin
       Result := 1;
       Break;
      end;
  end;
end;

constructor TBigInt.InitZero;
begin
  FIsNegative := False;
  FDigits := TDigits.Create(0);
end;

constructor TBigInt.Init(const AValue: Int64);
var
  absVal: Int64;
begin
  FIsNegative := AValue < 0;
  if AValue <> Int64.MinValue then
  begin
    absVal := Abs(AValue);
    if absVal >= CBase then
      FDigits := TDigits.Create(absVal mod CBase, absVal div CBase)
    else
      FDigits := TDigits.Create(absVal);
  end
  else begin
    FIsNegative := True;
    FDigits := TDigits.Create(0, 1 << 31);
  end;
end;

destructor TBigInt.Done;
begin
  SetLength(FDigits, 0);
end;

function TBigInt.CompareTo(constref AOther: TBigInt): Integer;
begin
  if IsNegative = AOther.IsNegative then
    Result := CompareToAbsoluteValues(AOther)
  else
    Result := 1;
  if IsNegative then
    Result := -Result;
end;

operator := (const A: Int64): TBigInt;
begin
  Result.Done;
  Result.Init(A);
end;

operator + (const A, B: TBigInt): TBigInt;
begin
  if A.IsNegative = B.IsNegative then
    Result := TBigInt.AddAbsoluteValues(A, B, A.IsNegative)
  else
    Result := TBigInt.SubtractAbsoluteValues(A, B, A.IsNegative);
end;

operator - (const A, B: TBigInt): TBigInt;
begin
  if A.IsNegative = B.IsNegative then
    Result := TBigInt.SubtractAbsoluteValues(A, B, A.IsNegative)
  else
    Result := TBigInt.AddAbsoluteValues(A, B, A.IsNegative);
end;

operator * (const A: TBigInt; const B: Int64): TBigInt;
begin
  if  (a = 0) or (b = 0) then
    Result.InitZero
  else
    Result := TBigInt.MultiplyAbsoluteValues(A, B, A.IsNegative = (B > 0));
end;

operator shl(const A: TBigInt; const B: Integer): TBigInt;
var
  i, j, digitShifts, bitShifts, reverseShift, len, newLength: Integer;
  lastDigit: Cardinal;
begin
  // Handles shift of zero.
  if A = 0 then
    Result.InitZero
  else begin
    // Determines full digit shifts and bit shifts.
    DivMod(B, CBitsPerDigit, digitShifts, bitShifts);

    if bitShifts > 0 then
    begin
      reverseShift := CBitsPerDigit - bitShifts;
      len := Length(A.FDigits);
      lastDigit := A.FDigits[len - 1] >> reverseShift;
      newLength := len + digitShifts;

      if lastDigit = 0 then
        SetLength(Result.FDigits, newLength)
      else
        SetLength(Result.FDigits, newLength + 1);

      // Performs full digit shifts by shifting the access index j for A.FDigits.
      Result.FDigits[digitShifts] := A.FDigits[0] << bitShifts;
      j := 0;
      for i := digitShifts + 1 to newLength - 1 do
      begin
        // Performs bit shifts.
        Result.FDigits[i] := A.FDigits[j] >> reverseShift;
        Inc(j);
        Result.FDigits[i] := Result.FDigits[i] or (A.FDigits[j] << bitShifts);
      end;

      if Length(Result.FDigits) > newLength then
        Result.FDigits[newLength] := lastDigit;
    end
    else begin
      // Performs full digit shifts by copy if there are no bit shifts.
      len := Length(A.FDigits);
      SetLength(Result.FDigits, len + digitShifts);
      Move(A.FDigits[0], Result.FDigits[digitShifts], Length(A.FDigits));
    end;

    Result.FIsNegative := A.IsNegative;
  end;
end;

operator = (const A: TBigInt; const B: Int64): Boolean;
begin
  Result := A.CompareTo(B) = 0;
end;

end.