Merge branch 'bigint' into day24-analytical

- Fixed some uses of Move
- Added Sign, string initializers (hexadecimal and binary), explicit
converter to Int64, comparison operators

AdventOfCode.lpi

@@ -137,6 +137,18 @@
         <Filename Value="solvers\UNeverTellMeTheOdds.pas"/>
         <IsPartOfProject Value="True"/>
       </Unit>
+      <Unit>
+        <Filename Value="UBigInt.pas"/>
+        <IsPartOfProject Value="True"/>
+      </Unit>
+      <Unit>
+        <Filename Value="UPolynomial.pas"/>
+        <IsPartOfProject Value="True"/>
+      </Unit>
+      <Unit>
+        <Filename Value="UPolynomialRoots.pas"/>
+        <IsPartOfProject Value="True"/>
+      </Unit>
     </Units>
   </ProjectOptions>
   <CompilerOptions>

UBigInt.pas

@@ -0,0 +1,708 @@
+{
+  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 in Freepascal of a C# class created in 2022.
+  TBigInt = object
+  private
+    FDigits: TDigits;
+    FIsNegative: Boolean;
+
+    function GetSign: Integer;
+
+    // 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;
+
+    class function GetZero: TBigInt; static;
+
+    // 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;
+
+    class function FromHexOrBinString(const AValue: string; const AFromBase: Integer): TBigInt; static;
+    class function ConvertDigitBlock(const AValue: string; var AStartIndex: Integer; const ACharBlockSize, AFromBase:
+      Integer): Cardinal;
+  public
+    property IsNegative: Boolean read FIsNegative;
+    property Sign: Integer read GetSign;
+    class property Zero: TBigInt read GetZero;
+
+    // Returns the index of the most significant bit, i.e. returns integer k, where 2^k is the largest power of 2 that
+    // is less than or equal to the absolute value of the number itself. Returns -1 if the given number is 0.
+    function GetMostSignificantBitIndex: Int64;
+    function CompareTo(constref AOther: TBigInt): Integer;
+    function TryToInt64(out AOutput: Int64): Boolean;
+    // TODO: ToString is currently for debug output only.
+    function ToString: string;
+    class function FromInt64(const AValue: Int64): TBigInt; static;
+    class function FromHexadecimalString(const AValue: string): TBigInt; static;
+    class function FromBinaryString(const AValue: string): TBigInt; static;
+  end;
+
+  { Operators }
+
+  operator := (const A: Int64): TBigInt;
+  operator - (const A: TBigInt): TBigInt;
+  operator + (const A, B: TBigInt): TBigInt;
+  operator - (const A, B: TBigInt): TBigInt;
+  operator * (const A, B: TBigInt): TBigInt;
+  operator shl (const A: TBigInt; const B: Integer): TBigInt;
+  operator = (const A, B: TBigInt): Boolean;
+  operator <> (const A, B: TBigInt): Boolean;
+  operator < (const A, B: TBigInt): Boolean;
+  operator <= (const A, B: TBigInt): Boolean;
+  operator > (const A, B: TBigInt): Boolean;
+  operator >= (const A, B: TBigInt): Boolean;
+
+implementation
+
+const
+  CBase = Cardinal.MaxValue + 1;
+  CMaxDigit = Cardinal.MaxValue;
+  CDigitSize = SizeOf(Cardinal);
+  CBitsPerDigit = CDigitSize * 8;
+  CHalfBits = CBitsPerDigit >> 1;
+  CHalfDigitMax = (1 << CHalfBits) - 1;
+
+  CZero: TBigInt = (FDigits: (0); FIsNegative: False);
+
+{ TBigInt }
+
+function TBigInt.GetSign: Integer;
+begin
+  if FIsNegative then
+    Result := -1
+  else if (Length(FDigits) > 1) or (FDigits[0] <> 0) then
+    Result := 1
+  else
+    Result := 0;
+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);
+  Result.FDigits := Copy(FDigits, AIndex, trimmedCount);
+  Result.FIsNegative := False;
+end;
+
+function TBigInt.CompareToAbsoluteValues(constref AOther: TBigInt): Integer;
+var
+  i: Integer;
+begin
+  Result := Length(FDigits) - Length(AOther.FDigits);
+  if Result = 0 then
+  begin
+    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;
+
+class function TBigInt.GetZero: TBigInt;
+begin
+  Result := CZero;
+end;
+
+class function TBigInt.AddAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
+var
+  i, j, 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
+      for j := i to len - 1 do
+        Result.FDigits[j] := AA.FDigits[j]
+    else
+      for j := i to len - 1 do
+        Result.FDigits[j] := AB.FDigits[j];
+
+  // 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, len);
+  Result.FIsNegative := AReturnNegative;
+end;
+
+class function TBigInt.SubtractAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
+var
+  a, b: TBigInt;
+  carry: Cardinal;
+  compare, i, j, lastNonZeroDigitIndex, len: Integer;
+begin
+  // Establishes the operand order, such that Abs(a) is not less than Abs(b).
+  compare := AA.CompareToAbsoluteValues(AB);
+  if compare = 0 then
+  begin
+    Result := Zero;
+    Exit;
+  end;
+
+  if compare > 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
+    for j := i to len - 1 do
+      Result.FDigits[j] := a.FDigits[j]
+  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 := Zero
+      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 := FromInt64(a1 + a0);
+      bm := FromInt64(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;
+
+class function TBigInt.FromHexOrBinString(const AValue: string; const AFromBase: Integer): TBigInt;
+var
+  charBlockSize, offset, i, j, k, remainder: Integer;
+  d: Cardinal;
+begin
+  // 2 ^ (32 / charBlockSize) = AFromBase
+  case AFromBase of
+    2: charBlockSize := 32;
+    16: charBlockSize := 8;
+  end;
+
+  if AValue[1] = '-' then
+  begin
+    offset := 2;
+    Result.FIsNegative := True;
+  end
+  else begin
+    offset := 1;
+    Result.FIsNegative := False;
+  end;
+
+  // Calculates the first (most significant) digit d of the result.
+  DivMod(AValue.Length - offset + 1, charBlockSize, i, remainder);
+  k := offset;
+  d := 0;
+  // Checks the first block of chars that is not a full block.
+  if remainder > 0 then
+    d := ConvertDigitBlock(AValue, k, remainder, AFromBase);
+  // Checks full blocks of chars for first digit.
+  while (d = 0) and (i > 0) do
+  begin
+    Dec(i);
+    d := ConvertDigitBlock(AValue, k, charBlockSize, AFromBase);
+  end;
+
+  // Checks for zero.
+  if (d = 0) and (i = 0) then
+    Result := Zero
+  else begin
+    // Initializes the array of digits.
+    SetLength(Result.FDigits, i + 1);
+    Result.FDigits[i] := d;
+
+    // Calculates the other digits.
+    for j := i - 1 downto 0 do
+      Result.FDigits[j] := ConvertDigitBlock(AValue, k, charBlockSize, AFromBase);
+  end;
+end;
+
+class function TBigInt.ConvertDigitBlock(const AValue: string; var AStartIndex: Integer; const ACharBlockSize,
+  AFromBase: Integer): Cardinal;
+var
+  part: string;
+begin
+  part := Copy(AValue, AStartIndex, ACharBlockSize);
+  Inc(AStartIndex, ACharBlockSize);
+  case AFromBase of
+    2: part := '%' + part;
+    8: part := '&' + part;
+    16: part := '$' + part;
+  end;
+  Result := StrToDWord(part);
+end;
+
+function TBigInt.GetMostSignificantBitIndex: Int64;
+var
+  high, i: Integer;
+begin
+  high := Length(FDigits) - 1;
+  if (high = 0) and (FDigits[0] = 0) then
+    Result := -1
+  else begin
+    i := CBitsPerDigit - 1;
+    while ((1 << i) and FDigits[high]) = 0 do
+      Dec(i);
+    Result := high * CBitsPerDigit + i;
+  end;
+end;
+
+function TBigInt.CompareTo(constref AOther: TBigInt): Integer;
+begin
+  if FIsNegative = AOther.FIsNegative then
+    Result := CompareToAbsoluteValues(AOther)
+  else
+    Result := 1;
+  if FIsNegative then
+    Result := -Result;
+end;
+
+function TBigInt.TryToInt64(out AOutput: Int64): Boolean;
+begin
+  AOutput := 0;
+  Result := False;
+  case Length(FDigits) of
+    0: Result := True;
+    1: begin
+      AOutput := FDigits[0];
+      if FIsNegative then
+        AOutput := -AOutput;
+      Result := True;
+    end;
+    2: begin
+      if FDigits[1] <= Integer.MaxValue then
+      begin
+        AOutput := FDigits[1] * CBase + FDigits[0];
+        if FIsNegative then
+          AOutput := -AOutput;
+        Result := True;
+      end
+      else if (FDigits[1] = Integer.MaxValue + 1) and (FDigits[0] = 0) and FIsNegative then
+      begin
+        AOutput := Int64.MinValue;
+        Result := True;
+      end;
+    end;
+  end;
+end;
+
+function TBigInt.ToString: string;
+var
+  i: Integer;
+begin
+  if FIsNegative then
+    Result := '-'
+  else
+    Result := '';
+  for i := 0 to Length(FDigits) - 2 do
+    Result := Result + '(' + IntToStr(FDigits[i]) + ' + 2^32 * ';
+  Result := Result + IntToStr(FDigits[Length(FDigits) - 1]);
+  for i := 0 to Length(FDigits) - 2 do
+    Result := Result + ')';
+end;
+
+class function TBigInt.FromInt64(const AValue: Int64): TBigInt;
+var
+  absVal: Int64;
+begin
+  if AValue <> Int64.MinValue then
+  begin
+    absVal := Abs(AValue);
+    if absVal >= CBase then
+      Result.FDigits := TDigits.Create(absVal mod CBase, absVal div CBase)
+    else
+      Result.FDigits := TDigits.Create(absVal);
+    Result.FIsNegative := AValue < 0;
+  end
+  else begin
+    Result.FDigits := TDigits.Create(0, 1 << 31);
+    Result.FIsNegative := True;
+  end;
+end;
+
+class function TBigInt.FromHexadecimalString(const AValue: string): TBigInt;
+begin
+  Result := FromHexOrBinString(AValue, 16);
+end;
+
+class function TBigInt.FromBinaryString(const AValue: string): TBigInt;
+begin
+  Result := FromHexOrBinString(AValue, 2);
+end;
+
+{ Operators }
+
+operator := (const A: Int64): TBigInt;
+begin
+  Result := TBigInt.FromInt64(A);
+  //WriteLn(':=a op: ', Result.ToString);
+end;
+
+operator - (const A: TBigInt): TBigInt;
+var
+  len: Integer;
+begin
+  len := Length(A.FDigits);
+  if (len > 1) or (A.FDigits[0] > 0) then
+  begin
+    Result.FDigits := Copy(A.FDigits, 0, len);
+    Result.FIsNegative := not A.FIsNegative;
+  end
+  else
+    Result := TBigInt.Zero;
+  //WriteLn('            a: ', A.ToString);
+  //WriteLn('-a  op: ', Result.ToString);
+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);
+  //WriteLn('            a: ', A.ToString);
+  //WriteLn('            b: ', B.ToString);
+  //WriteLn('a+b op: ', Result.ToString);
+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);
+  //WriteLn('            a: ', A.ToString);
+  //WriteLn('            b: ', B.ToString);
+  //WriteLn('a-b op: ', Result.ToString);
+end;
+
+operator * (const A, B: TBigInt): TBigInt;
+begin
+  if  (A = TBigInt.Zero) or (B = TBigInt.Zero) then
+    Result := TBigInt.Zero
+  else
+    Result := TBigInt.MultiplyAbsoluteValues(A, B, A.IsNegative <> B.IsNegative);
+  //WriteLn('            a: ', A.ToString);
+  //WriteLn('            b: ', B.ToString);
+  //WriteLn('a*b op: ', Result.ToString);
+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 := TBigInt.Zero
+  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);
+      for i := 0 to digitShifts - 1 do
+        Result.FDigits[i] := 0;
+      for i := 0 to len - 1 do
+        Result.FDigits[i + digitShifts] := A.FDigits[i];
+    end;
+
+    Result.FIsNegative := A.IsNegative;
+  end;
+  //WriteLn('            a: ', A.ToString);
+  //WriteLn('a<< op: ', Result.ToString);
+end;
+
+operator = (const A, B: TBigInt): Boolean;
+begin
+  Result := A.CompareTo(B) = 0;
+end;
+
+operator <> (const A, B: TBigInt): Boolean;
+begin
+  Result := A.CompareTo(B) <> 0;
+end;
+
+operator < (const A, B: TBigInt): Boolean;
+begin
+  Result := A.CompareTo(B) < 0;
+end;
+
+operator <= (const A, B: TBigInt): Boolean;
+begin
+  Result := A.CompareTo(B) <= 0;
+end;
+
+operator > (const A, B: TBigInt): Boolean;
+begin
+  Result := A.CompareTo(B) > 0;
+end;
+
+operator >= (const A, B: TBigInt): Boolean;
+begin
+  Result := A.CompareTo(B) >= 0;
+end;
+
+end.
+

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.
 

UPolynomial.pas

@@ -0,0 +1,157 @@
+{
+  Solutions to the Advent Of Code.
+  Copyright (C) 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 UPolynomial;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, UBigInt;
+
+type
+  TInt64Array = array of Int64;
+
+  { TBigIntPolynomial }
+
+  TBigIntPolynomial = object
+  private
+    FCoefficients: array of TBigInt;
+    function GetDegree: Integer;
+    function GetCoefficient(const AIndex: Integer): TBigInt;
+  public
+    property Degree: Integer read GetDegree;
+    property Coefficient[const AIndex: Integer]: TBigInt read GetCoefficient;
+    function CalcValueAt(const AX: Int64): TBigInt;
+    function IsEqualTo(const AOther: TBigIntPolynomial): Boolean;
+    function ScaleVariable(const AScaleFactor: TBigInt): TBigIntPolynomial;
+    function ToString: string;
+    class function Create(const ACoefficients: array of TBigInt): TBigIntPolynomial; static;
+  end;
+
+  { Operators }
+
+  operator = (const A, B: TBigIntPolynomial): Boolean;
+  operator <> (const A, B: TBigIntPolynomial): Boolean;
+
+implementation
+
+{ TBigIntPolynomial }
+
+function TBigIntPolynomial.GetDegree: Integer;
+begin
+  Result := Length(FCoefficients) - 1;
+end;
+
+function TBigIntPolynomial.GetCoefficient(const AIndex: Integer): TBigInt;
+begin
+  Result := FCoefficients[AIndex];
+end;
+
+function TBigIntPolynomial.CalcValueAt(const AX: Int64): TBigInt;
+var
+  i: Integer;
+begin
+  Result := TBigInt.Zero;
+  for i := High(FCoefficients) downto 0 do
+    Result := Result * AX + FCoefficients[i];
+end;
+
+function TBigIntPolynomial.IsEqualTo(const AOther: TBigIntPolynomial): Boolean;
+var
+  i: Integer;
+begin
+  if Length(FCoefficients) = Length(AOther.FCoefficients) then
+  begin
+    Result := True;
+    for i := 0 to Length(FCoefficients) - 1 do
+      if FCoefficients[i] <> AOther.FCoefficients[i] then
+      begin
+        Result := False;
+        Break;
+      end;
+  end
+  else
+    Result := False;
+end;
+
+function TBigIntPolynomial.ScaleVariable(const AScaleFactor: TBigInt): TBigIntPolynomial;
+var
+  len, i: Integer;
+  factor: TBigInt;
+begin
+  if AScaleFactor <> TBigInt.Zero then
+  begin
+    len := Length(FCoefficients);
+    SetLength(Result.FCoefficients, len);
+    Result.FCoefficients[0] := FCoefficients[0];
+    factor := AScaleFactor;
+    for i := 1 to len - 1 do begin
+      Result.FCoefficients[i] := FCoefficients[i] * factor;
+      factor := factor * AScaleFactor;
+    end;
+  end
+  else
+    SetLength(Result.FCoefficients, 0);
+end;
+
+function TBigIntPolynomial.ToString: string;
+var
+  i: Integer;
+begin
+  Result := FCoefficients[0].ToString;
+  for i := 1 to Length(FCoefficients) - 1 do
+    if i > 1 then
+      Result := Result + ' + ' + FCoefficients[i].ToString + ' * x^' + IntToStr(i)
+    else
+      Result := Result + ' + ' + FCoefficients[i].ToString + ' * x';
+end;
+
+class function TBigIntPolynomial.Create(const ACoefficients: array of TBigInt): TBigIntPolynomial;
+var
+  high, i: integer;
+begin
+  high := -1;
+  for i := Length(ACoefficients) - 1 downto 0 do
+    if ACoefficients[i] <> 0 then
+    begin
+      high := i;
+      Break;
+    end;
+  if high >= 0 then
+  begin
+    SetLength(Result.FCoefficients, high + 1);
+    for i := 0 to high do
+      Result.FCoefficients[i] := ACoefficients[i];
+  end;
+end;
+
+{ Operators }
+
+operator = (const A, B: TBigIntPolynomial): Boolean;
+begin
+  Result := A.IsEqualTo(B);
+end;
+
+operator <> (const A, B: TBigIntPolynomial): Boolean;
+begin
+  Result := not A.IsEqualTo(B);
+end;
+
+end.
+

UPolynomialRoots.pas

@@ -0,0 +1,72 @@
+{
+  Solutions to the Advent Of Code.
+  Copyright (C) 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 UPolynomialRoots;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, UPolynomial, UBigInt;
+
+type
+
+  { TRootIsolation }
+
+  TRootIsolation = class
+  private
+    function CalcSimpleRootBound(constref APolynomial: TBigIntPolynomial): TBigInt;
+  public
+    function Bisect(constref APolynomial: TBigIntPolynomial): Int64;
+  end;
+
+implementation
+
+{ TRootIsolation }
+
+function TRootIsolation.CalcSimpleRootBound(constref APolynomial: TBigIntPolynomial): TBigInt;
+var
+  i, sign: Integer;
+  a: TBigInt;
+begin
+  // We need a_n > 0 here, so we use -sign(a_n) instead of actually flipping the polynomial.
+  // Sign is not 0 because a_n is not 0.
+  sign := -APolynomial.Coefficient[APolynomial.Degree].Sign;
+
+  // This is a simplification of Cauchy's bound to avoid division.
+  // https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Bounds_of_positive_real_roots
+  Result := TBigInt.Zero;
+  for i := 0 to APolynomial.Degree - 1 do begin
+    a := sign * APolynomial.Coefficient[i];
+    if Result < a then
+      Result := a;
+  end;
+  Result := Result + 1;
+end;
+
+function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial): Int64;
+var
+  bound: TBigInt;
+  p: TBigIntPolynomial;
+begin
+  bound := CalcSimpleRootBound(APolynomial);
+  p := APolynomial.ScaleVariable(bound);
+end;
+
+end.
+

solvers/UNeverTellMeTheOdds.pas

@@ -1,6 +1,6 @@
 {
   Solutions to the Advent Of Code.
-  Copyright (C) 2023  Stefan Müller
+  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
@@ -22,26 +22,45 @@ unit UNeverTellMeTheOdds;
 interface
 
 uses
-  Classes, SysUtils, Generics.Collections, Math, USolver;
+  Classes, SysUtils, Generics.Collections, Math, matrix, USolver, UNumberTheory, UBigInt;
 
 type
 
   { THailstone }
 
-  THailstone = record
-    X, Y, Z: Int64;
-    VX, VY, VZ: Integer;
+  THailstone = class
+  public
+    Position, Velocity: Tvector3_extended;
+    constructor Create(const ALine: string);
+    constructor Create(const APosition, AVelocity: Tvector3_extended);
   end;
 
-  THailstones = specialize TList<THailstone>;
+  THailstones = specialize TObjectList<THailstone>;
+
+  { TFirstCollisionPolynomial }
+
+  TFirstCollisionPolynomial = class
+  private
+    FA: array[0..10] of TBigInt;
+    FH: array[0..6] of TBigInt;
+    procedure NormalizeCoefficients;
+  public
+    procedure Init(constref AHailstone1, AHailstone2, AHailstone3: THailstone; const t_0, t_1, t_2: Int64);
+    function EvaluateAt(const AT0: Int64): TBigInt;
+    function CalcPositiveIntegerRoot: Int64;
+    function CalcT1(const AT0: Int64): Int64;
+  end;
 
   { TNeverTellMeTheOdds }
 
   TNeverTellMeTheOdds = class(TSolver)
   private
     FMin, FMax: Int64;
-    FHailStones: THailstones;
+    FHailstones: THailstones;
+    FA: array[0..10] of Int64;
+    FH: array[0..6] of Int64;
     function AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
+    procedure FindRockThrow(const AIndex1, AIndex2, AIndex3: Integer);
   public
     constructor Create(const AMin: Int64 = 200000000000000; const AMax: Int64 = 400000000000000);
     destructor Destroy; override;
@@ -51,8 +70,452 @@ type
     function GetPuzzleName: string; override;
   end;
 
+const
+  CIterationThreshold = 0.00001;
+  CEpsilon = 0.0000000001;
+
 implementation
 
+{ THailstone }
+
+constructor THailstone.Create(const ALine: string);
+var
+  split: TStringArray;
+begin
+  split := ALine.Split([',', '@']);
+  Position.init(
+    StrToFloat(Trim(split[0])),
+    StrToFloat(Trim(split[1])),
+    StrToFloat(Trim(split[2])));
+  Velocity.init(
+    StrToFloat(Trim(split[3])),
+    StrToFloat(Trim(split[4])),
+    StrToFloat(Trim(split[5])));
+end;
+
+constructor THailstone.Create(const APosition, AVelocity: Tvector3_extended);
+begin
+  Position := APosition;
+  Velocity := AVelocity;
+end;
+
+{ TFirstCollisionPolynomial }
+
+procedure TFirstCollisionPolynomial.NormalizeCoefficients;
+var
+  shift: Integer;
+  i: Low(FA)..High(FA);
+  //gcd: TBigInt;
+begin
+  // Eliminates zero constant term.
+  shift := 0;
+  while (shift <= High(FA)) and (FA[shift] = 0) do
+    Inc(shift);
+
+  if shift <= High(FA) then
+  begin
+    if shift > 0 then
+    begin
+      for i := Low(FA) to High(FA) - shift do
+        FA[i] := FA[i + shift];
+      for i := High(FA) - shift + 1 to High(FA) do
+        FA[i] := 0;
+    end;
+
+    //// Finds GCD of all coefficients.
+    //gcd := FA[Low(FA)];
+    //for i := Low(FA) + 1 to High(FA) do
+    //  if FA[i] <> 0 then
+    //    gcd := TNumberTheory.GreatestCommonDivisor(gcd, FA[i]);
+    //WriteLn('GCD: ', gcd);
+    //
+    //for i := Low(FA) to High(FA) do
+    //  FA[i] := FA[i] div gcd;
+  end;
+
+  //WriteLn('(', FA[10], ') * x^10 + (', FA[9], ') * x^9 + (', FA[8], ') * x^8 + (', FA[7], ') * x^7 + (',
+  //  FA[6], ') * x^6 + (', FA[5], ') * x^5 + (', FA[4], ') * x^4 + (', FA[3], ') * x^3 + (', FA[2], ') * x^2 + (',
+  //  FA[1], ') * x + (', FA[0], ')');
+end;
+
+procedure TFirstCollisionPolynomial.Init(constref AHailstone1, AHailstone2, AHailstone3: THailstone; const t_0, t_1,
+  t_2: Int64);
+var
+  P_00, P_01, P_02, P_10, P_11, P_12, P_20, P_21, P_22,
+  V_00, V_01, V_02, V_10, V_11, V_12, V_20, V_21, V_22: Int64;
+  k: array[0..139] of TBigInt;
+  // For debug calculations
+  act, a_1, a_2, b_0, b_1, c_0, c_1, d_0, d_1, e_0, e_1, f_0, f_1, f_2: Int64;
+begin
+  // Solving this non-linear equation system, with velocities V_i and start positions P_i:
+  //   V_0 * t_0 + P_0 = V_x * t_0 + P_x
+  //   V_1 * t_1 + P_1 = V_x * t_1 + P_x
+  //   V_2 * t_2 + P_2 = V_x * t_2 + P_x
+  // Which gives:
+  //   P_x = (V_0 - V_x) * t_0 + P_0
+  //   V_x = (V_0 * t_0 - V_1 * t_1 + P_0 - P_1) / (t_0 - t_1)
+  // And with vertex components:
+  //   1: 0 = (t_1 - t_0) * (V_00 * t_0 - V_20 * t_2 + P_00 - P_20) - (t_2 - t_0) * (V_00 * t_0 - V_10 * t_1 + P_00 - P_10)
+  //   2: t_1 = (((V_01 - V_21) * t_2 + P_11 - P_21) * t_0 + (P_01 - P_11) * t_2)
+  //            / ((V_01 - V_11) * t_0 + (V_11 - V_21) * t_2 + P_01 - P_21)
+  //   3: t_2 = (((V_02 - V_12) * t_1 + P_22 - P_12) * t_0 + (P_02 - P_22) * t_1)
+  //            / ((V_02 - V_22) * t_0 + (V_22 - V_12) * t_1 + P_02 - P_12)
+  //   for t_0, t_1, t_2 not pairwise equal.
+  // With some substitutions depending only on t_0 this gives
+  //   1: 0 = (t_1 - t_0) * (f_2 - V_20 * t_2) - (t_2 - t_0) * (f_1 - V_10 * t_1)
+  //   2: t_1 = (b_0 + b_1 * t_2) / (c_0 + c_1 * t_2)
+  //   3: t_2 = (d_0 + d_1 * t_1) / (e_0 + e_1 * t_1)
+  // And 3 in 2 gives:
+  //   4: g_2 * t_1^2 - g_1 * t_1 - g_0 = 0
+  // Then, with 4 and 3 in 1 and lengthy calculations with many substitutions (see constants k below, now independent of
+  // t_0), the following polynomial can be constructed, with t_0 being a positive integer root of this polynomial.
+  //   y = a_10 * x^10 + a_9 * x^9 + ... + a_0
+
+  P_00 := Round(AHailstone1.Position.data[0]);
+  P_01 := Round(AHailstone1.Position.data[1]);
+  P_02 := Round(AHailstone1.Position.data[2]);
+  P_10 := Round(AHailstone2.Position.data[0]);
+  P_11 := Round(AHailstone2.Position.data[1]);
+  P_12 := Round(AHailstone2.Position.data[2]);
+  P_20 := Round(AHailstone3.Position.data[0]);
+  P_21 := Round(AHailstone3.Position.data[1]);
+  P_22 := Round(AHailstone3.Position.data[2]);
+  V_00 := Round(AHailstone1.Velocity.data[0]);
+  V_01 := Round(AHailstone1.Velocity.data[1]);
+  V_02 := Round(AHailstone1.Velocity.data[2]);
+  V_10 := Round(AHailstone2.Velocity.data[0]);
+  V_11 := Round(AHailstone2.Velocity.data[1]);
+  V_12 := Round(AHailstone2.Velocity.data[2]);
+  V_20 := Round(AHailstone3.Velocity.data[0]);
+  V_21 := Round(AHailstone3.Velocity.data[1]);
+  V_22 := Round(AHailstone3.Velocity.data[2]);
+
+  k[0] := P_00 - P_20;
+  k[1] := P_00 - P_10;
+  k[2] := P_11 - P_21;
+  k[3] := P_01 - P_11;
+  k[4] := P_01 - P_21;
+  k[5] := P_22 - P_12;
+  k[6] := P_02 - P_22;
+  k[7] := P_02 - P_12;
+  k[8] := V_11 - V_21;
+  k[9] := V_22 - V_12;
+  k[10] := V_01 - V_21;
+  k[11] := V_01 - V_11;
+  k[12] := V_02 - V_12;
+  k[13] := V_02 - V_22;
+
+  FH[0] := k[11] * k[9] + k[8] * k[12];
+  FH[1] := k[4] * k[9] + k[8] * k[6];
+  FH[2] := k[11] * k[13] - k[10] * k[12];
+  FH[3] := k[11] * k[7] + k[4] * k[13] + k[8] * k[5] - k[2] * k[9] - k[10] * k[6] - k[3] * k[12];
+  FH[4] := k[4] * k[7] - k[3] * k[6];
+  FH[5] := k[10] * k[5] + k[2] * k[13];
+  FH[6] := k[3] * k[5] + k[2] * k[7];
+
+  k[14] := V_00 * k[9] - V_20 * k[12];
+  k[15] := k[0] * k[9] - V_20 * k[6];
+  k[16] := V_00 * k[13];
+  k[17] := V_00 * k[7] + k[0] * k[13] - V_20 * k[5];
+  k[18] := k[0] * k[7];
+  k[19] := k[5] - k[7];
+  k[20] := 2 * FH[2] * FH[3];
+  k[21] := FH[3] * FH[3];
+  k[22] := k[21] + 2 * FH[2] * FH[4];
+  k[23] := 2 * FH[3] * FH[4];
+  k[24] := 2 * FH[0] * FH[1];
+  k[25] := FH[0] * FH[0];      // KILL?
+  k[26] := FH[5] * k[25];    // KILL?
+  k[126] := FH[5] * FH[0];
+  k[127] := FH[5] * FH[1] + FH[6] * FH[0];
+  k[128] := FH[6] * FH[1];
+  k[27] := FH[5] * k[24] + FH[6] * k[25];   // KILL?
+  k[28] := FH[1] * FH[1];                // KILL?
+  k[29] := FH[5] * k[28] + FH[6] * k[24];   // KILL?
+  k[30] := FH[6] * k[28];        // KILL?
+  k[31] := FH[2] * FH[2];
+  k[132] := k[20] + 4 * k[126];
+  k[133] := k[22] + 4 * k[127];
+  k[134] := k[23] + 4 * k[128];
+  k[32] := k[31] + 4 * k[26];      // KILL?
+  k[33] := k[20] + 4 * k[27];       // KILL?
+  k[34] := k[22] + 4 * k[29];      // KILL?
+  k[35] := k[23] + 4 * k[30];     // KILL?
+  k[36] := k[31] + 2 * k[26];    // KILL?
+  k[37] := k[20] + 2 * k[27];    // KILL?
+  k[38] := k[22] + 2 * k[29];    // KILL?
+  k[39] := k[23] + 2 * k[30];      // KILL?
+  k[137] := k[20] + 2 * k[126];
+  k[138] := k[22] + 2 * k[127];
+  k[139] := k[23] + 2 * k[128];
+  k[40] := k[14] + V_10 * (k[12] - k[9]);
+  k[41] := k[15] + V_10 * k[6];
+  k[42] := k[16] - k[14] - V_10 * k[13] - (k[12] - k[9]) * V_00;
+  k[43] := k[17] - k[15] + V_10 * k[19] - (k[12] - k[9]) * k[1] - k[6] * V_00;
+  k[44] := k[18] - k[6] * k[1];
+  k[45] := k[42] * FH[0] - k[40] * FH[2];
+  k[46] := k[42] * FH[1] + k[43] * FH[0] - k[41] * FH[2] - k[40] * FH[3];
+  k[47] := k[43] * FH[1] + k[44] * FH[0] - k[41] * FH[3] - k[40] * FH[4];
+  k[48] := k[44] * FH[1] - k[41] * FH[4];
+  k[49] := k[42] * FH[2];
+  k[50] := k[40] * k[31] - k[49] * FH[0];
+  k[51] := k[42] * FH[3] + k[43] * FH[2];
+  k[52] := k[40] * k[137] + k[41] * k[31] - k[51] * FH[0] - k[49] * FH[1];
+  k[53] := k[42] * FH[4] + k[43] * FH[3] + k[44] * FH[2];
+  k[54] := k[40] * k[138] + k[41] * k[137] - k[53] * FH[0] - k[51] * FH[1];
+  k[55] := k[43] * FH[4] + k[44] * FH[3];
+  k[56] := k[40] * k[139] + k[41] * k[138] - k[55] * FH[0] - k[53] * FH[1];
+  k[57] := k[44] * FH[4];
+  k[58] := FH[4] * FH[4];
+  k[59] := k[40] * k[58] + k[41] * k[139] - k[57] * FH[0] - k[55] * FH[1];
+  k[60] := k[41] * k[58] - k[57] * FH[1];
+  k[61] := k[13] * V_00 - k[16];
+  k[62] := 2 * k[25] * k[61];
+  k[63] := k[13] * k[1] - k[19] * V_00 - k[17];
+  k[64] := 2 * (k[24] * k[61] + k[25] * k[63]);
+  k[65] := - k[19] * k[1] - k[18];
+  k[66] := 2 * (k[28] * k[61] + k[24] * k[63] + k[25] * k[65]);
+  k[67] := 2 * (k[28] * k[63] + k[24] * k[65]);
+  k[68] := 2 * k[28] * k[65];
+  k[69] := k[50] + k[62];
+  k[70] := k[52] + k[64];
+  k[71] := k[54] + k[66];
+  k[72] := k[56] + k[67];
+  k[73] := k[59] + k[68];
+  k[74] := k[45] * k[45];
+  k[75] := 2 * k[45] * k[46];
+  k[76] := k[46] * k[46] + 2 * k[45] * k[47];
+  k[77] := 2 * (k[45] * k[48] + k[46] * k[47]);
+  k[78] := k[47] * k[47] + 2 * k[46] * k[48];
+  k[79] := 2 * k[47] * k[48];
+  k[80] := k[48] * k[48];
+
+  FA[0] := k[58] * k[80] - k[60] * k[60];
+  FA[1] := k[134] * k[80] + k[58] * k[79] - 2 * k[73] * k[60];
+  FA[2] := k[133] * k[80] + k[134] * k[79] + k[58] * k[78] - k[73] * k[73] - 2 * k[72] * k[60];
+  FA[3] := k[133] * k[79] + k[134] * k[78] + k[58] * k[77] + k[132] * k[80]
+    - 2 * (k[71] * k[60] + k[72] * k[73]);
+  FA[4] := k[31] * k[80] + k[133] * k[78] + k[134] * k[77] + k[58] * k[76] + k[132] * k[79] - k[72] * k[72]
+    - 2 * (k[70] * k[60] + k[71] * k[73]);
+  FA[5] := k[31] * k[79] + k[133] * k[77] + k[134] * k[76] + k[58] * k[75] + k[132] * k[78]
+    - 2 * (k[69] * k[60] + k[70] * k[73] + k[71] * k[72]);
+  FA[6] := k[31] * k[78] + k[133] * k[76] + k[134] * k[75] + k[58] * k[74] + k[132] * k[77] - k[71] * k[71]
+    - 2 * (k[69] * k[73] + k[70] * k[72]);
+  FA[7] := k[31] * k[77] + k[133] * k[75] + k[134] * k[74] + k[132] * k[76] - 2 * (k[69] * k[72] + k[70] * k[71]);
+  FA[8] := k[31] * k[76] + k[132] * k[75] + k[133] * k[74] - k[70] * k[70] - 2 * k[69] * k[71];
+  FA[9] := k[31] * k[75] + k[132] * k[74] - 2 * k[69] * k[70];
+  FA[10] := k[31] * k[74] - k[69] * k[69];
+
+  // Debug calculations
+  //a_1 := V_00 * t_0 + P_00 - P_20;
+  //a_2 := V_00 * t_0 + P_00 - P_10;
+  //b_0 := (P_11 - P_21) * t_0;
+  //b_1 := (V_01 - V_21) * t_0 + P_01 - P_11;
+  //c_0 := (V_01 - V_11) * t_0 + P_01 - P_21;
+  //c_1 := V_11 - V_21;
+  //d_0 := (P_22 - P_12) * t_0;
+  //d_1 := (V_02 - V_12) * t_0 + P_02 - P_22;
+  //e_0 := (V_02 - V_22) * t_0 + P_02 - P_12;
+  //e_1 := V_22 - V_12;
+  //f_2 := c_0 * e_1 + c_1 * d_1;
+  //f_1 := c_0 * e_0 + c_1 * d_0 - b_0 * e_1 - b_1 * d_1;
+  //f_0 := b_1 * d_0 + b_0 * e_0;
+  //
+  //act := f_2 * t_1 * t_1 + f_1 * t_1 - f_0;
+  //Write('debug10: ', 0 = act, '    ');
+  //
+  //if f_2 <> 0 then
+  //begin
+  //  act := Round(- f_1 / (2 * f_2) + Sqrt((f_1 / (2 * f_2)) * (f_1 / (2 * f_2)) + f_0 / f_2));
+  //  Write('debug15: ', t_1 = act);
+  //  act := Round(- f_1 / (2 * f_2) - Sqrt((f_1 / (2 * f_2)) * (f_1 / (2 * f_2)) + f_0 / f_2));
+  //  Write(' OR ', t_1 = act, '    ');
+  //end;
+  //
+  //act := (e_0 + e_1 * t_1) * t_2 - (d_0 + d_1 * t_1);
+  //Write('debug20: ', 0 = act, '    ');
+  //
+  //act := (a_1 * e_1 - V_20 * d_1 + V_10 * (d_1 - e_1 * t_0)) * t_1 * t_1
+  //  + (a_1 * e_0 - V_20 * d_0 - t_0 * (a_1 * e_1  - V_20 * d_1) - (d_1 - e_1 * t_0) * a_2 + V_10 * (d_0 - e_0 * t_0)) * t_1
+  //  + t_0 * (V_20 * d_0 - a_1 * e_0) + (e_0 * t_0 - d_0) * a_2;
+  //Write('debug30: ', 0 = act, '    ');
+  //
+  //act := Round((a_1 * e_1 - V_20 * d_1 + V_10 * (d_1 - e_1 * t_0)) * (f_1 * f_1 + 2 * f_0 * f_2 - f_1 * Sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //  + (a_1 * e_0 - V_20 * d_0 - t_0 * (a_1 * e_1 - V_20 * d_1) - (d_1 - e_1 * t_0) * a_2 + V_10 * (d_0 - e_0 * t_0)) * (- f_1 * f_2 + f_2 * Sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //  + t_0 * (V_20 * d_0 - a_1 * e_0) * 2 * f_2 * f_2 + (e_0 * t_0 - d_0) * a_2 * 2 * f_2 * f_2);
+  //Write('debug40: ', 0 = act, '    ');
+  //
+  //Write('debug41: ',
+  //  a_1 * k[9] - V_20 * d_1
+  //  = k[14] * t_0 + k[15], '    ');
+  //Write('debug42: ',
+  //  d_1 - k[9] * t_0
+  //  = (k[12] - k[9]) * t_0 + k[6], '    ');
+  //Write('debug43: ',
+  //  a_1 * e_0 - V_20 * d_0
+  //  = k[16] * t_0 * t_0 + k[17] * t_0 + k[18], '    ');
+  //Write('debug44: ',
+  //  d_0 - e_0 * t_0
+  //  = - k[13] * t_0 * t_0 + k[19] * t_0, '    ');
+  //Write('debug45: ',
+  //  f_1 * f_1
+  //  = FH[2] * FH[2] * t_0 * t_0 * t_0 * t_0 + k[20] * t_0 * t_0 * t_0 + k[22] * t_0 * t_0 + k[23] * t_0 + FH[4] * FH[4], '    ');
+  //Write('debug46: ',
+  //  f_2 * f_2
+  //  = FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1], '    ');
+  //Write('debug47: ',
+  //  f_0 * f_2
+  //  = k[126] * t_0 * t_0 * t_0 + k[127] * t_0 * t_0 + k[128] * t_0, '    ');
+  //Write('debug48: ',
+  //  f_1 * f_1 + 4 * f_0 * f_2
+  //  = k[31] * t_0 * t_0 * t_0 * t_0 + k[132] * t_0 * t_0 * t_0 + k[133] * t_0 * t_0 + k[134] * t_0 + k[58], '    ');
+  //Write('debug49: ',
+  //  f_1 * f_1 + 2 * f_0 * f_2
+  //  = k[31] * t_0 * t_0 * t_0 * t_0 + k[137] * t_0 * t_0 * t_0 + k[138] * t_0 * t_0 + k[139] * t_0 + k[58], '    ');
+  //
+  //act := Round((k[14] * t_0 + k[15] + V_10 * ((k[12] - k[9]) * t_0 + k[6])) * (k[31] * t_0 * t_0 * t_0 * t_0 + k[137] * t_0 * t_0 * t_0 + k[138] * t_0 * t_0 + k[139] * t_0 + k[58] - f_1 * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //  + (k[16] * t_0 * t_0 + k[17] * t_0 + k[18] - t_0 * (k[14] * t_0 + k[15]) - ((k[12] - k[9]) * t_0 + k[6]) * a_2 - V_10 * (k[13] * t_0 * t_0 - k[19] * t_0)) * (- f_1 * f_2 + f_2 * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //  - 2 * t_0 * (k[16] * t_0 * t_0 + k[17] * t_0 + k[18]) * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1]) + 2 * (k[13] * t_0 * t_0 - k[19] * t_0) * a_2 * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1]));
+  //Write('debug50: ', 0 = act, '    ');
+  //
+  //Write('debug53: ',
+  //  0 = Round((k[40] * t_0 + k[41]) * (k[31] * t_0 * t_0 * t_0 * t_0 + k[137] * t_0 * t_0 * t_0 + k[138] * t_0 * t_0 + k[139] * t_0 + k[58] - f_1 * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //    + ((k[16] - k[14] - V_10 * k[13] - (k[12] - k[9]) * V_00) * t_0 * t_0 + (k[17] - k[15] + V_10 * k[19] - (k[12] - k[9]) * k[1] - k[6] * V_00) * t_0 + k[18] - k[6] * k[1]) * (- f_1 * f_2 + f_2 * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //    - 2 * t_0 * (k[16] * t_0 * t_0 + k[17] * t_0 + k[18]) * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1]) + 2 * (k[13] * t_0 * t_0 - k[19] * t_0) * a_2 * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1])),
+  //    '    ');
+  //
+  //Write('debug55: ',
+  //  0 = Round((k[40] * t_0 + k[41]) * (k[31] * t_0 * t_0 * t_0 * t_0 + k[137] * t_0 * t_0 * t_0 + k[138] * t_0 * t_0 + k[139] * t_0 + k[58])
+  //    - (k[40] * t_0 + k[41]) * f_1 * sqrt(f_1 * f_1 + 4 * f_0 * f_2)
+  //    + (k[42] * t_0 * t_0 + k[43] * t_0 + k[44]) * (- f_1 * f_2 + f_2 * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //    - 2 * t_0 * (k[16] * t_0 * t_0 + k[17] * t_0 + k[18]) * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1]) + 2 * (k[13] * t_0 * t_0 - k[19] * t_0) * a_2 * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1])),
+  //    '    ');
+  //
+  //Write('debug70: ',
+  //  0 = Round(((k[42] * t_0 * t_0 + k[43] * t_0 + k[44]) * (FH[0] * t_0 + FH[1]) - (k[40] * t_0 + k[41]) * (FH[2] * t_0 * t_0 + FH[3] * t_0 + FH[4])) * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+  //    + (k[40] * t_0 + k[41]) * (k[31] * t_0 * t_0 * t_0 * t_0 + k[137] * t_0 * t_0 * t_0 + k[138] * t_0 * t_0 + k[139] * t_0 + k[58])
+  //    - (k[42] * t_0 * t_0 + k[43] * t_0 + k[44]) * (FH[2] * t_0 * t_0 + FH[3] * t_0 + FH[4]) * (FH[0] * t_0 + FH[1])
+  //    - 2 * t_0 * (k[16] * t_0 * t_0 + k[17] * t_0 + k[18]) * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1]) + 2 * (k[13] * t_0 * t_0 - k[19] * t_0) * (V_00 * t_0 + k[1]) * (FH[0] * FH[0] * t_0 * t_0 + k[24] * t_0 + FH[1] * FH[1]),
+  //    '    ');
+//
+//  Write('debug73: ',
+//    0 = Round((
+//        (k[42] * FH[0] - k[40] * FH[2]) * t_0 * t_0 * t_0
+//        + (k[42] * FH[1] + k[43] * FH[0] - k[41] * FH[2] - k[40] * FH[3]) * t_0 * t_0
+//        + (k[43] * FH[1] + k[44] * FH[0] - k[41] * FH[3] - k[40] * FH[4]) * t_0
+//        + k[44] * FH[1] - k[41] * FH[4]
+//    ) * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+//    + (k[40] * k[31] - k[42] * FH[2] * FH[0]) * t_0 * t_0 * t_0 * t_0 * t_0
+//    + (k[40] * k[137] + k[41] * k[31] - k[42] * FH[3] * FH[0] - k[43] * FH[2] * FH[0] - k[42] * FH[2] * FH[1]) * t_0 * t_0 * t_0 * t_0
+//    + (k[40] * k[138] + k[41] * k[137] - k[42] * FH[4] * FH[0] - k[43] * FH[3] * FH[0] - k[44] * FH[2] * FH[0] - k[42] * FH[3] * FH[1] - k[43] * FH[2] * FH[1]) * t_0 * t_0 * t_0
+//    + (k[40] * k[139] + k[41] * k[138] - k[43] * FH[4] * FH[0] - k[44] * FH[3] * FH[0] - k[42] * FH[4] * FH[1] - k[43] * FH[3] * FH[1] - k[44] * FH[2] * FH[1]) * t_0 * t_0
+//    + (k[40] * k[58] + k[41] * k[139] - k[44] * FH[4] * FH[0] - k[43] * FH[4] * FH[1] - k[44] * FH[3] * FH[1]) * t_0
+//    + k[41] * k[58] - k[44] * FH[4] * FH[1]
+//    + 2 * (k[13] * V_00 * FH[0] * FH[0] - k[16] * FH[0] * FH[0]) * t_0 * t_0 * t_0 * t_0 * t_0
+//    + 2 * (k[13] * V_00 * k[24] + k[13] * k[1] * FH[0] * FH[0] - k[19] * V_00 * FH[0] * FH[0] - k[16] * k[24] - k[17] * FH[0] * FH[0]) * t_0 * t_0 * t_0 * t_0
+//    + 2 * (k[13] * V_00 * FH[1] * FH[1] + k[13] * k[1] * k[24] - k[19] * V_00 * k[24] - k[19] * k[1] * FH[0] * FH[0] - k[16] * FH[1] * FH[1] - k[17] * k[24] - k[18] * FH[0] * FH[0]) * t_0 * t_0 * t_0
+//    + 2 * (k[13] * k[1] * FH[1] * FH[1] - k[19] * V_00 * FH[1] * FH[1] - k[19] * k[1] * k[24] - k[17] * FH[1] * FH[1] - k[18] * k[24]) * t_0 * t_0
+//    + 2 * (- k[19] * k[1] * FH[1] * FH[1] - k[18] * FH[1] * FH[1]) * t_0,
+//    '    ');
+//
+//  Write('debug78: ',
+//    0 = Round((k[45] * t_0 * t_0 * t_0 + k[46] * t_0 * t_0 + k[47] * t_0 + k[48]) * sqrt(f_1 * f_1 + 4 * f_0 * f_2))
+//      + (k[50] + k[62]) * t_0 * t_0 * t_0 * t_0 * t_0 + (k[52] + k[64]) * t_0 * t_0 * t_0 * t_0 + (k[54] + k[66]) * t_0 * t_0 * t_0 + (k[56] + k[67]) * t_0 * t_0 + (k[59] + k[68]) * t_0 + k[60],
+//      '    ');
+//
+//  Write('debug80: ',
+//    0 = Round((k[45] * t_0 * t_0 * t_0 + k[46] * t_0 * t_0 + k[47] * t_0 + k[48]) * sqrt(k[31] * t_0 * t_0 * t_0 * t_0 + k[132] * t_0 * t_0 * t_0 + k[133] * t_0 * t_0 + k[134] * t_0 + k[58])
+//      + k[69] * t_0 * t_0 * t_0 * t_0 * t_0 + k[70] * t_0 * t_0 * t_0 * t_0 + k[71] * t_0 * t_0 * t_0 + k[72] * t_0 * t_0 + k[73] * t_0 + k[60]),
+//      '    ');
+//  WriteLn;
+//  WriteLn('  0 = ((', k[45], ') * x^3 + (', k[46], ') * x^2 + (', k[47], ') * x + (', k[48], ')) * sqrt((', k[31], ') * x^4 + (', k[132], ') * x^3 + (', k[133], ') * x^2 + (', k[134], ') * x + (', k[58], ')) + (',
+//    k[69], ') * x^5 + (', k[70], ') * x^4 + (', k[71], ') * x^3 + (', k[72], ') * x^2 + (', k[73], ') * x + (', k[60], ')');
+
+  Write('debug83: ',
+    (k[45] * t_0 * t_0 * t_0 + k[46] * t_0 * t_0 + k[47] * t_0 + k[48]) * (k[45] * t_0 * t_0 * t_0 + k[46] * t_0 * t_0 + k[47] * t_0 + k[48]) * (k[31] * t_0 * t_0 * t_0 * t_0 + k[132] * t_0 * t_0 * t_0 + k[133] * t_0 * t_0 + k[134] * t_0 + k[58]) =
+      (k[69] * t_0 * t_0 * t_0 * t_0 * t_0 + k[70] * t_0 * t_0 * t_0 * t_0 + k[71] * t_0 * t_0 * t_0 + k[72] * t_0 * t_0 + k[73] * t_0 + k[60]) * (k[69] * t_0 * t_0 * t_0 * t_0 * t_0 + k[70] * t_0 * t_0 * t_0 * t_0 + k[71] * t_0 * t_0 * t_0 + k[72] * t_0 * t_0 + k[73] * t_0 + k[60]),
+      '    ');
+  Write('debug85: ',
+    0 =
+      (
+          k[45] * k[45] * t_0 * t_0 * t_0 * t_0 * t_0 * t_0
+          + 2 * k[45] * k[46] * t_0 * t_0 * t_0 * t_0 * t_0
+          + k[46] * k[46] * t_0 * t_0 * t_0 * t_0
+          + 2 * k[45] * k[47] * t_0 * t_0 * t_0 * t_0
+          + 2 * k[45] * k[48] * t_0 * t_0 * t_0
+          + 2 * k[46] * k[47] * t_0 * t_0 * t_0
+          + k[47] * k[47] * t_0 * t_0
+          + 2 * k[46] * k[48] * t_0 * t_0
+          + 2 * k[47] * k[48] * t_0
+          + k[48] * k[48]
+      ) * (k[31] * t_0 * t_0 * t_0 * t_0 + k[132] * t_0 * t_0 * t_0 + k[133] * t_0 * t_0 + k[134] * t_0 + k[58])
+      - k[69] * k[69] * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0
+      - 2 * k[69] * k[70] * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0
+      - (k[70] * k[70] + 2 * k[69] * k[71]) * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0
+      - 2 * (k[69] * k[72] + k[70] * k[71]) * t_0 * t_0 * t_0 * t_0 * t_0 * t_0 * t_0
+      - (k[71] * k[71] + 2 * k[69] * k[73] + 2 * k[70] * k[72]) * t_0 * t_0 * t_0 * t_0 * t_0 * t_0
+      - 2 * (k[69] * k[60] + k[70] * k[73] + k[71] * k[72]) * t_0 * t_0 * t_0 * t_0 * t_0
+      - (k[72] * k[72] + 2 * k[70] * k[60] + 2 * k[71] * k[73]) * t_0 * t_0 * t_0 * t_0
+      - 2 * (k[71] * k[60] + k[72] * k[73]) * t_0 * t_0 * t_0
+      - (k[73] * k[73] + 2 * k[72] * k[60]) * t_0 * t_0
+      - 2 * k[73] * k[60] * t_0
+      - k[60] * k[60],
+      '    ');
+
+  WriteLn('debug96: ', EvaluateAt(t_0) = 0);
+
+  NormalizeCoefficients;
+
+  WriteLn('debug99: ', EvaluateAt(t_0) = 0, '    ');
+end;
+
+function TFirstCollisionPolynomial.EvaluateAt(const AT0: Int64): TBigInt;
+var
+  i: Low(FA)..High(FA);
+begin
+  Result := TBigInt.Zero;
+  for i := High(FA) downto Low(FA) do
+    Result := Result * AT0 + FA[i];
+end;
+
+function TFirstCollisionPolynomial.CalcPositiveIntegerRoot: Int64;
+var
+  dividers: TDividers;
+  factors: TInt64Array;
+  divider: Int64;
+begin
+  Result := 0;
+  //factors := TIntegerFactorization.PollardsRhoAlgorithm(FA[0]);
+  //dividers := TDividers.Create(factors);
+  //
+  //try
+  //for divider in dividers do
+  //begin
+  //  //WriteLn('Check if ', divider, ' is a root...');
+  //  if EvaluateAt(divider) = 0 then
+  //  begin
+  //    Result := divider;
+  //    Break;
+  //  end;
+  //end;
+  //
+  //finally
+  //  dividers.Free;
+  //end;
+end;
+
+function TFirstCollisionPolynomial.CalcT1(const AT0: Int64): Int64;
+var
+  g_0, g_1, g_2: Int64;
+  g: Extended;
+begin
+  //g_2 := FH[0] * AT0 + FH[1];
+  //g_1 := FH[2] * AT0 * AT0 + FH[3] * AT0 + FH[4];
+  //g_0 := FH[5] * AT0 * AT0 + FH[6] * AT0;
+  //g := - g_1 / (2 * g_2);
+  //Result := Round(g + sqrt(g * g + g_0));
+end;
+
 { TNeverTellMeTheOdds }
 
 function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
@@ -60,58 +523,113 @@ var
   m1, m2, x, y: Double;
 begin
   Result := False;
-  m1 := AHailstone1.VY / AHailstone1.VX;
-  m2 := AHailstone2.VY / AHailstone2.VX;
+  m1 := AHailstone1.Velocity.data[1] / AHailstone1.Velocity.data[0];
+  m2 := AHailstone2.Velocity.data[1] / AHailstone2.Velocity.data[0];
   if m1 <> m2 then
   begin
-    x := (AHailstone2.Y - m2 * AHailstone2.X - AHailstone1.Y + m1 * AHailstone1.X) / (m1 - m2);
+    x := (AHailstone2.Position.data[1] - m2 * AHailstone2.Position.data[0]
+      - AHailstone1.Position.data[1] + m1 * AHailstone1.Position.data[0])
+      / (m1 - m2);
     if (FMin <= x) and (x <= FMax)
-    and (x * Sign(AHailstone1.VX) >= AHailstone1.X * Sign(AHailstone1.VX))
-    and (x * Sign(AHailstone2.VX) >= AHailstone2.X * Sign(AHailstone2.VX)) then
+    and (x * Sign(AHailstone1.Velocity.data[0]) >= AHailstone1.Position.data[0] * Sign(AHailstone1.Velocity.data[0]))
+    and (x * Sign(AHailstone2.Velocity.data[0]) >= AHailstone2.Position.data[0] * Sign(AHailstone2.Velocity.data[0]))
+    then
     begin
-      y := m1 * (x - AHailstone1.X) + AHailstone1.Y;
+      y := m1 * (x - AHailstone1.Position.data[0]) + AHailstone1.Position.data[1];
       if (FMin <= y) and (y <= FMax) then
         Result := True
     end;
   end;
 end;
 
+// For debug calculations:
+Const
+  T : array[0..4] of Byte = (5, 3, 4, 6, 1);
+
+procedure TNeverTellMeTheOdds.FindRockThrow(const AIndex1, AIndex2, AIndex3: Integer);
+var
+  //i, j, k: Integer;
+  //x0, x1, x2: Extended;
+  f: TFirstCollisionPolynomial;
+  t0, t1: Int64;
+  p, v: Tvector3_extended;
+  test: TBigInt;
+begin
+  WriteLn;
+  WriteLn(AIndex1, ' ', AIndex2, ' ', AIndex3);
+  f := TFirstCollisionPolynomial.Create;
+  f.Init(FHailstones[AIndex1], FHailstones[AIndex2], FHailstones[AIndex3], T[AIndex1], T[AIndex2], T[AIndex3]);
+  //t0 := f.CalcPositiveIntegerRoot;
+  //WriteLn('t0: ', t0, ' ', t0 = T[AIndex1]);
+  //t1 := f.CalcT1(t0);
+  //WriteLn(', t1: ', t1);
+  f.Free;
+
+  //// V_x = (V_0 * t_0 - V_1 * t_1 + P_0 - P_1) / (t_0 - t_1)
+  //v := (FHailstones[AIndex1].Velocity * t0 - FHailstones[AIndex2].Velocity * t1
+  //  + FHailstones[AIndex1].Position - FHailstones[AIndex2].Position) / (t0 - t1);
+  //// P_x = (V_0 - V_x) * t_0 + P_0
+  //p := (FHailstones[AIndex1].Velocity - v) * t0 + FHailstones[AIndex1].Position;
+  //FPart2 := Round(p.data[0]) + Round(p.data[1]) + Round(p.data[2]);
+
+  //for i := 0 to FHailstones.Count - 3 do
+  //  for j := i + 1 to FHailstones.Count - 2 do
+  //    for k:= j + 1 to FHailstones.Count - 1 do
+  //    begin
+  //      WriteLn(i, j, k);
+  //      solver := TRockThrowSolver.Create(FHailstones[i], FHailstones[j], FHailstones[k], 0);
+  //      case i of
+  //        0: x0 := 5;
+  //        1: x0 := 3;
+  //        2: x0 := 4;
+  //      end;
+  //      f := solver.CalcValue(x0);
+  //      solver.Free;
+  //    end;
+
+  //for i := 80 to 120 do
+  //begin
+  //  solver := TRockThrowSolver.Create(FHailstones[0], FHailstones[1], FHailstones[2], 0);
+  //  x0 := i / 20;
+  //  f := solver.CalcValue(x0);
+  //  WriteLn(x0, ' ', f.Valid, ' ', f.Value);
+  //  solver.Free;
+  //end;
+end;
+
 constructor TNeverTellMeTheOdds.Create(const AMin: Int64; const AMax: Int64);
 begin
   FMin := AMin;
   FMax := AMax;
-  FHailStones := THailstones.Create;
+  FHailstones := THailstones.Create;
 end;
 
 destructor TNeverTellMeTheOdds.Destroy;
 begin
-  FHailStones.Free;
+  FHailstones.Free;
   inherited Destroy;
 end;
 
 procedure TNeverTellMeTheOdds.ProcessDataLine(const ALine: string);
-var
-  split: TStringArray;
-  hailstone: THailstone;
 begin
-  split := ALine.Split([',', '@']);
-  hailstone.X := StrToInt64(Trim(split[0]));
-  hailstone.Y := StrToInt64(Trim(split[1]));
-  hailstone.Z := StrToInt64(Trim(split[2]));
-  hailstone.VX := StrToInt(Trim(split[3]));
-  hailstone.VY := StrToInt(Trim(split[4]));
-  hailstone.VZ := StrToInt(Trim(split[5]));
-  FHailStones.Add(hailstone);
+  FHailstones.Add(THailstone.Create(ALine));
 end;
 
 procedure TNeverTellMeTheOdds.Finish;
 var
-  i, j: Integer;
+  i, j, k: Integer;
 begin
-  for i := 0 to FHailStones.Count - 2 do
-    for j := i + 1 to FHailStones.Count - 1 do
-      if AreIntersecting(FHailStones[i], FHailStones[j]) then
+  for i := 0 to FHailstones.Count - 2 do
+    for j := i + 1 to FHailstones.Count - 1 do
+      if AreIntersecting(FHailstones[i], FHailstones[j]) then
         Inc(FPart1);
+
+  for i := 0 to FHailstones.Count - 1 do
+    for j := 0 to FHailstones.Count - 1 do
+      for k := 0 to FHailstones.Count - 1 do
+        if (i <> j) and (i <> k) and (j <> k) then
+          FindRockThrow(i, j, k);
+  //FindRockThrow(0, 1, 2);
 end;
 
 function TNeverTellMeTheOdds.GetDataFileName: string;

tests/AdventOfCodeFPCUnit.lpi

@@ -40,10 +40,6 @@
         <Filename Value="UGearRatiosTestCases.pas"/>
         <IsPartOfProject Value="True"/>
       </Unit>
-      <Unit>
-        <Filename Value="..\USolver.pas"/>
-        <IsPartOfProject Value="True"/>
-      </Unit>
       <Unit>
         <Filename Value="UBaseTestCases.pas"/>
         <IsPartOfProject Value="True"/>
@@ -140,6 +136,18 @@
         <Filename Value="UNeverTellMeTheOddsTestCases.pas"/>
         <IsPartOfProject Value="True"/>
       </Unit>
+      <Unit>
+        <Filename Value="UPolynomialTestCases.pas"/>
+        <IsPartOfProject Value="True"/>
+      </Unit>
+      <Unit>
+        <Filename Value="UPolynomialRootsTestCases.pas"/>
+        <IsPartOfProject Value="True"/>
+      </Unit>
+      <Unit>
+        <Filename Value="UBigIntTestCases.pas"/>
+        <IsPartOfProject Value="True"/>
+      </Unit>
     </Units>
   </ProjectOptions>
   <CompilerOptions>

tests/AdventOfCodeFPCUnit.lpr

@@ -9,7 +9,7 @@ uses
   UHotSpringsTestCases, UPointOfIncidenceTestCases, UParabolicReflectorDishTestCases, ULensLibraryTestCases,
   UFloorWillBeLavaTestCases, UClumsyCrucibleTestCases, ULavaductLagoonTestCases, UAplentyTestCases,
   UPulsePropagationTestCases, UStepCounterTestCases, USandSlabsTestCases, ULongWalkTestCases,
-  UNeverTellMeTheOddsTestCases;
+  UNeverTellMeTheOddsTestCases, UBigIntTestCases, UPolynomialTestCases, UPolynomialRootsTestCases;
 
 {$R *.res}
 

tests/UBigIntTestCases.pas

@@ -0,0 +1,960 @@
+{
+  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 UBigIntTestCases;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, UBigInt;
+
+type
+  // TODO: TBigIntAbsTestCase
+  // TODO: TBigIntFromDecTestCase
+  //   TestCase('80091110002223289436355210101231765016693209176113648', 10, '80091110002223289436355210101231765016693209176113648');
+  //   TestCase('-3201281944789146858325672846129872035678639439423746384198412894', 10, '-3201281944789146858325672846129872035678639439423746384198412894');
+  //   TestCase('000000000000000000000000000000000000000080091', 10, '80091');
+  //   TestCase('0000000000000000000000000000000000000000000000', 10, '0');
+  //   TestCase('-000000000000000000000000000000004687616873215464763146843215486643215', 10, '-4687616873215464763146843215486643215');
+  //   TestTryFromDec
+  // TODO: TBigIntFromOctTestCase
+  //   OctValue, 8, DecValue:
+  //   TestCase('3132521', 8, '832849');
+  //   TestCase('772350021110002226514', 8, '9123449598017875276');
+  //   TestCase('-247515726712721', 8, '-11520970560977');
+  //   TestCase('000000000000000000000000000000000000000072', 8, '58');
+  //   TestCase('000000000000000000000000000000000000000000000000000000000', 8, '0');
+  //   OctValue, HexValue:
+  //   TestCase('7440737076307076026772', '3C83BE3E638F82DFA');
+  //   TestCase('1336625076203401614325664', '16F6547C8380E31ABB4');
+  //   TestCase('3254410316166274457257121050201413373225', '6AC84338ECBC97ABCA22840C2DF695');
+  //   TestCase('1441330072562106272767270117307274151276215', '642D81D5C88CBAFBAE09EC75E1A57C8D');
+  //   TestCase('3533665174054677131163436413756431236774257271133', '3ADED4F82CDF964E71E85FBA329EFE2BD725B');
+  //   TestCase('57346415317074055631627141161054073014006076155710653', '5EE686B3C782DCCE5CC271160EC18061F1B791AB');
+  //   TestCase('000000000000000000000000000000245745251354033134', 8, '5838773085550172');
+  //   TestTryFromOct
+  // TODO: TBigIntUnaryPlusTestCase
+  // TODO: TBigIntBitwiseLogicalTestCase
+  // TODO: TBigIntComplementTestCase
+  // TODO: TBigIntConversionTestCase
+  // TODO: TBigIntIncrementDecrementTestCase
+  // TODO: TBigIntQuotientTestCase
+  // TODO: TBigIntShiftRightTestCase
+
+  { TBigIntSignTestCase }
+
+  TBigIntSignTestCase = class(TTestCase)
+  private
+    procedure Test(const AHexValue: string; const AExpectedSign: Integer);
+  published
+    procedure TestZero;
+    procedure TestShortPositive;
+    procedure TestShortNegative;
+    procedure TestLongPositive;
+    procedure TestLongNegative;
+  end;
+
+  { TBigIntMostSignificantBitIndexTestCase }
+
+  TBigIntMostSignificantBitIndexTestCase = class(TTestCase)
+  private
+    procedure Test(const ABinValue: string; const AExpectedIndex: Int64);
+  published
+    procedure TestZero;
+    procedure TestShort;
+    procedure TestLong;
+  end;
+
+  { TBigIntFromInt64TestCase }
+
+  TBigIntFromInt64TestCase = class(TTestCase)
+  private
+    procedure Test(const AValue: Int64);
+  published
+    procedure TestShortPositive;
+    procedure TestShortNegative;
+    procedure TestLongPositive;
+    procedure TestLongNegative;
+    procedure TestZero;
+  end;
+
+  { TBigIntFromHexTestCase }
+
+  TBigIntFromHexTestCase = class(TTestCase)
+  private
+    procedure TestShort(const AHexValue: string; const ADecValue: Int64);
+  published
+    procedure TestPositive;
+    procedure TestNegative;
+    procedure TestZero;
+    procedure TestLeadingZeros;
+    // TODO: TestTryFromHex
+  end;
+
+  { TBigIntFromBinTestCase }
+
+  TBigIntFromBinTestCase = class(TTestCase)
+  private
+    procedure TestShort(const ABinValue: string; const ADecValue: Int64);
+  published
+    procedure TestPositive;
+    procedure TestNegative;
+    procedure TestLeadingZeros;
+    // TODO: TestTryFromBin
+  end;
+
+  { TBigIntUnaryMinusTestCase }
+
+  TBigIntUnaryMinusTestCase = class(TTestCase)
+  private
+    procedure Test(const AValue: Int64);
+    procedure TestTwice(const AValue: Int64);
+  published
+    procedure TestZero;
+    procedure TestPositive;
+    procedure TestNegative;
+    procedure TestPositiveTwice;
+    procedure TestNegativeTwice;
+  end;
+
+  { TBigIntSumTestCase }
+
+  TBigIntSumTestCase = class(TTestCase)
+  private
+    procedure Test(const AHexValueLeft, AHexValueRight, AHexValueSum: string);
+  published
+    procedure TestShort;
+    procedure TestPositivePlusPositive;
+    procedure TestNegativePlusNegative;
+    procedure TestLargePositivePlusSmallNegative;
+    procedure TestSmallPositivePlusLargeNegative;
+    procedure TestLargeNegativePlusSmallPositive;
+    procedure TestSmallNegativePlusLargePositive;
+    procedure TestZeroPlusPositive;
+    procedure TestPositivePlusZero;
+    procedure TestZero;
+    procedure TestSumShorterLeft;
+    procedure TestSumShorterRight;
+    procedure TestSumEndsInCarry;
+    procedure TestSumEndsInCarryNewDigit;
+    procedure TestSumMultiCarry;
+  end;
+
+  { TBigIntDifferenceTestCase }
+
+  TBigIntDifferenceTestCase = class(TTestCase)
+  private
+    procedure Test(const AHexValueMinuend, AHexValueSubtrahend, AHexValueDifference: string);
+  published
+    procedure TestShort;
+    procedure TestLargePositiveMinusSmallPositive;
+    procedure TestSmallPositiveMinusLargePositive;
+    procedure TestLargeNegativeMinusSmallNegative;
+    procedure TestSmallNegativeMinusLargeNegative;
+    procedure TestNegativeMinusPositive;
+    procedure TestPositiveMinusNegative;
+    procedure TestLargeOperands;
+    procedure TestPositiveMinusZero;
+    procedure TestZeroMinusPositive;
+    procedure TestZero;
+    procedure TestDifferenceShorterLeft;
+    procedure TestDifferenceShorterRight;
+    procedure TestDifferenceEndsInCarry;
+    procedure TestDifferenceEndsInCarryLosingDigit;
+    procedure TestDifferenceMultiCarry;
+  end;
+
+  { TBigIntProductTestCase }
+
+  TBigIntProductTestCase = class(TTestCase)
+  private
+    procedure Test(const AHexValueLeft, AHexValueRight, AHexValueProduct: string);
+  published
+    procedure TestShort;
+    procedure TestLongPositiveTimesPositive;
+    procedure TestLongPositiveTimesNegative;
+    procedure TestLongNegativeTimesPositive;
+    procedure TestLongNegativeTimesNegative;
+    procedure TestZeroTimesPositive;
+    procedure TestZeroTimesNegative;
+    procedure TestZero;
+  end;
+
+  { TBigIntShiftLeftTestCase }
+
+  TBigIntShiftLeftTestCase = class(TTestCase)
+  private
+    procedure Test(const AHexValueOperand: string; const AShift: Integer; const AHexValueResult: string);
+  published
+    procedure TestShort;
+    procedure TestShortWithCarry;
+    procedure TestLongWithCarry;
+    procedure TestLongWithMultiDigitCarry;
+    procedure TestWithAlignedDigits;
+    procedure TestZero;
+  end;
+
+  { TBigIntEqualityTestCase }
+
+  TBigIntEqualityTestCase = class(TTestCase)
+  private
+    procedure TestEqual(const AValue: Int64);
+    procedure TestEqualHex(const AHexValue: string);
+    procedure TestNotEqualHex(const AHexValueLeft, AHexValueRight: string);
+  published
+    procedure TestShortEqual;
+    procedure TestLongEqual;
+    procedure TestZeroEqual;
+    procedure TestShortNotEqual;
+    procedure TestLongNotEqualSign;
+    procedure TestZeroNotEqual;
+  end;
+
+  { TBigIntComparisonTestCase }
+
+  TBigIntComparisonTestCase = class(TTestCase)
+  private
+    procedure TestLessThan(const AHexValueLeft, AHexValueRight: string);
+    procedure TestGreaterThan(const AHexValueLeft, AHexValueRight: string);
+    procedure TestEqual(const AHexValue: string);
+  published
+    procedure TestLessSameLength;
+    procedure TestLessShorterLeft;
+    procedure TestLessNegativeShorterRight;
+    procedure TestGreaterSameLength;
+    procedure TestGreaterShorterRight;
+    procedure TestGreaterNegativeShorterLeft;
+    procedure TestEqualPositive;
+    procedure TestEqualNegative;
+    procedure TestEqualZero;
+  end;
+
+implementation
+
+{ TBigIntSignTestCase }
+
+procedure TBigIntSignTestCase.Test(const AHexValue: string; const AExpectedSign: Integer);
+var
+  a: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValue);
+  AssertEquals(AExpectedSign, a.Sign);
+end;
+
+procedure TBigIntSignTestCase.TestZero;
+begin
+  Test('0', 0);
+end;
+
+procedure TBigIntSignTestCase.TestShortPositive;
+begin
+  Test('36A', 1);
+end;
+
+procedure TBigIntSignTestCase.TestShortNegative;
+begin
+  Test('-29B7', -1);
+end;
+
+procedure TBigIntSignTestCase.TestLongPositive;
+begin
+  Test('1240168AB09ABDF8283B77124', 1);
+end;
+
+procedure TBigIntSignTestCase.TestLongNegative;
+begin
+  Test('-192648F1305DD04757A24C873F29F60B6184B5', -1);
+end;
+
+{ TBigIntMostSignificantBitIndexTestCase }
+
+procedure TBigIntMostSignificantBitIndexTestCase.Test(const ABinValue: string; const AExpectedIndex: Int64);
+var
+  a: TBigInt;
+begin
+  a := TBigInt.FromBinaryString(ABinValue);
+  AssertEquals('BigInt from binary string ''' + ABinValue + ''' did not have its most significant bit at index '''
+    + IntToStr(AExpectedIndex) + '''.',
+    AExpectedIndex, a.GetMostSignificantBitIndex);
+end;
+
+procedure TBigIntMostSignificantBitIndexTestCase.TestZero;
+begin
+  Test('0', -1);
+end;
+
+procedure TBigIntMostSignificantBitIndexTestCase.TestShort;
+begin
+  Test('1', 0);
+  Test('11111101111001', 13);
+  Test('10010111101100001110110101111000', 31);
+end;
+
+procedure TBigIntMostSignificantBitIndexTestCase.TestLong;
+begin
+  Test('111100010101010100011101010100011', 32);
+  Test('11101001101010111101000101010001010101010101111111001010101010010100101000101011111001000111001001100011', 103);
+  Test('111101100011110110111011010111100000000001010111101110101101101100101010110111101011010101001100', 95);
+end;
+
+{ TBigIntFromInt64TestCase }
+
+procedure TBigIntFromInt64TestCase.Test(const AValue: Int64);
+var
+  a: TBigInt;
+  act: Int64;
+begin
+  a := TBigInt.FromInt64(AValue);
+  AssertTrue('BigInt from ''' + IntToStr(AValue) + ''' could not be converted back to an Int64 value.',
+    a.TryToInt64(act));
+  AssertEquals('BigInt from ''' + IntToStr(AValue) + ''' converted back to Int64 was not equal to initial value.',
+    AValue, act);
+end;
+
+procedure TBigIntFromInt64TestCase.TestShortPositive;
+begin
+  Test(17);
+  Test(4864321);
+  Test(Cardinal.MaxValue);
+end;
+
+procedure TBigIntFromInt64TestCase.TestShortNegative;
+begin
+  Test(-54876);
+  Test(Integer.MinValue);
+end;
+
+procedure TBigIntFromInt64TestCase.TestLongPositive;
+begin
+  Test(5800754643214654);
+  Test(Int64.MaxValue);
+end;
+
+procedure TBigIntFromInt64TestCase.TestLongNegative;
+begin
+  Test(-94136445555883);
+  Test(Int64.MinValue);
+end;
+
+procedure TBigIntFromInt64TestCase.TestZero;
+begin
+  Test(0);
+end;
+
+{ TBigIntFromHexTestCase }
+
+procedure TBigIntFromHexTestCase.TestShort(const AHexValue: string; const ADecValue: Int64);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValue);
+  b := TBigInt.FromInt64(ADecValue);
+  AssertTrue('BigInt from hexadecimal string ''' + AHexValue + ''' was not equal to ''' + IntToStr(ADecValue) + '''.',
+    a = b);
+end;
+
+procedure TBigIntFromHexTestCase.TestPositive;
+begin
+  TestShort('91', 145);
+  TestShort('800E111000222', 2252766466540066);
+  TestShort('DE0B227802AB233', 999995000000000563);
+  TestShort('19C0048155000', 453000000000000);
+  TestShort('3B9ACA00', 1000000000);
+end;
+
+procedure TBigIntFromHexTestCase.TestNegative;
+begin
+  TestShort('-47', -71);
+  TestShort('-800E111000222', -2252766466540066);
+  TestShort('-4512FF3', -72429555);
+end;
+
+procedure TBigIntFromHexTestCase.TestZero;
+begin
+  TestShort('0', 0);
+  TestShort('-0', 0);
+end;
+
+procedure TBigIntFromHexTestCase.TestLeadingZeros;
+begin
+  TestShort('000000000000000000000000000000004B', 75);
+  TestShort('-00000000000000000000000000000000A2', -162);
+  TestShort('00000000000000000000000000000000FF452849DB01', 280672493755137);
+  TestShort('000000000000000000000000000000000000', 0);
+  TestShort('-000000000000000000000000000000000000', 0);
+end;
+
+{ TBigIntFromBinTestCase }
+
+procedure TBigIntFromBinTestCase.TestShort(const ABinValue: string; const ADecValue: Int64);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromBinaryString(ABinValue);
+  b := TBigInt.FromInt64(ADecValue);
+  AssertTrue('BigInt from binary string ''' + ABinValue + ''' was not equal to ''' + IntToStr(ADecValue) + '''.',
+    a = b);
+end;
+
+procedure TBigIntFromBinTestCase.TestPositive;
+begin
+  TestShort('110101010101101101010010110000101010100101001001010100110', 120109162101379750);
+end;
+
+procedure TBigIntFromBinTestCase.TestNegative;
+begin
+  TestShort('-11100100111010100111000111100011110000100110110111100101010010', -4123780452057839954);
+end;
+
+procedure TBigIntFromBinTestCase.TestLeadingZeros;
+begin
+  TestShort('0000000000000000000000000000000000000000000000000000000000000000000000111', 7);
+  TestShort('0000000000000000000000000000000000000000000000000000000000000000000000000000000', 0);
+end;
+
+{ TBigIntUnaryMinusTestCase }
+
+procedure TBigIntUnaryMinusTestCase.Test(const AValue: Int64);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromInt64(AValue);
+  b := TBigInt.FromInt64(-AValue);
+  AssertTrue('Negative BigInt from ''' + IntToStr(AValue) + ''' was not equal to the BigInt from the negative value.',
+    b = -a);
+end;
+
+procedure TBigIntUnaryMinusTestCase.TestTwice(const AValue: Int64);
+var
+  a: TBigInt;
+begin
+  a := TBigInt.FromInt64(AValue);
+  AssertTrue('BigInt from ''' + IntToStr(AValue) + '''was not equal to the double negative of itself.',
+    a = -(-a));
+end;
+
+procedure TBigIntUnaryMinusTestCase.TestZero;
+begin
+  Test(0);
+end;
+
+procedure TBigIntUnaryMinusTestCase.TestPositive;
+begin
+  Test(234972358233);
+end;
+
+procedure TBigIntUnaryMinusTestCase.TestNegative;
+begin
+  Test(-999214927400);
+end;
+
+procedure TBigIntUnaryMinusTestCase.TestPositiveTwice;
+begin
+  TestTwice(8647613456601);
+end;
+
+procedure TBigIntUnaryMinusTestCase.TestNegativeTwice;
+begin
+  TestTwice(-38600421308534);
+end;
+
+{ TBigIntSumTestCase }
+
+procedure TBigIntSumTestCase.Test(const AHexValueLeft, AHexValueRight, AHexValueSum: string);
+var
+  a, b, s: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueLeft);
+  b := TBigInt.FromHexadecimalString(AHexValueRight);
+  s := TBigInt.FromHexadecimalString(AHexValueSum);
+  AssertTrue('Sum of BigInt from ''' + AHexValueLeft + ''' and from ''' + AHexValueRight
+    + ''' was not equal to the BigInt from ''' + AHexValueSum + '''.',
+    s = a + b);
+end;
+
+procedure TBigIntSumTestCase.TestShort;
+begin
+  Test('5', '7', 'C');
+end;
+
+procedure TBigIntSumTestCase.TestPositivePlusPositive;
+begin
+  Test('7000822000111', '9000911000333', '10001133000444');
+end;
+
+procedure TBigIntSumTestCase.TestNegativePlusNegative;
+begin
+  Test('-129B92A84643D608141', '-4574DBA206951ECFE', '-12E10783E84A6B26E3F');
+end;
+
+procedure TBigIntSumTestCase.TestLargePositivePlusSmallNegative;
+begin
+  Test('87BAD26984', '-8DCB20461', '7EDE206523');
+end;
+
+procedure TBigIntSumTestCase.TestSmallPositivePlusLargeNegative;
+begin
+  Test('A58301E4006', '-9851DA0FD433', '-8DF9A9F1942D');
+end;
+
+procedure TBigIntSumTestCase.TestLargeNegativePlusSmallPositive;
+begin
+  Test('-1FDB60CB5698870', '99CB1E00DE', '-1FDB572EA4B8792');
+end;
+
+procedure TBigIntSumTestCase.TestSmallNegativePlusLargePositive;
+begin
+  Test('-1ED598BBFEC2', '59CD4F02ECB56', '57DFF5772CC94');
+end;
+
+procedure TBigIntSumTestCase.TestZeroPlusPositive;
+begin
+  Test('0', '9BB000911FF5A000333', '9BB000911FF5A000333');
+end;
+
+procedure TBigIntSumTestCase.TestPositivePlusZero;
+begin
+  Test('23009605A16BCBB294A1FD', '0', '23009605A16BCBB294A1FD');
+end;
+
+procedure TBigIntSumTestCase.TestZero;
+begin
+  Test('0', '0', '0');
+end;
+
+procedure TBigIntSumTestCase.TestSumShorterLeft;
+begin
+  Test('3FFFF', '9000911000222', '9000911040221');
+end;
+
+procedure TBigIntSumTestCase.TestSumShorterRight;
+begin
+  Test('9000911000555', '3000EEEE', '900094100F443');
+end;
+
+procedure TBigIntSumTestCase.TestSumEndsInCarry;
+begin
+  Test('4000444000220', 'FFFFFFFF', '400054400021F');
+end;
+
+procedure TBigIntSumTestCase.TestSumEndsInCarryNewDigit;
+begin
+  Test('99990000', '99990055', '133320055');
+end;
+
+procedure TBigIntSumTestCase.TestSumMultiCarry;
+begin
+  Test('FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF', 'FFFFFFFF', '1000000000000000000000000FFFFFFFE');
+end;
+
+{ TBigIntDifferenceTestCase }
+
+procedure TBigIntDifferenceTestCase.Test(const AHexValueMinuend, AHexValueSubtrahend, AHexValueDifference: string);
+var
+  a, b, s: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueMinuend);
+  b := TBigInt.FromHexadecimalString(AHexValueSubtrahend);
+  s := TBigInt.FromHexadecimalString(AHexValueDifference);
+  AssertTrue('Difference of BigInt from ''' + AHexValueMinuend + ''' and from ''' + AHexValueSubtrahend
+    + ''' was not equal to the BigInt from ''' + AHexValueDifference + '''.',
+    s = a - b);
+end;
+
+procedure TBigIntDifferenceTestCase.TestShort;
+begin
+  Test('230', '6D', '1C3');
+end;
+
+procedure TBigIntDifferenceTestCase.TestLargePositiveMinusSmallPositive;
+begin
+  Test('910AAD86E5002455', '72DE020F932A1', '91037FA6C406F1B4');
+end;
+
+procedure TBigIntDifferenceTestCase.TestSmallPositiveMinusLargePositive;
+begin
+  Test('3127541F4C0AA', '3818CD9FBE652B', '-3506585DC9A481');
+end;
+
+procedure TBigIntDifferenceTestCase.TestLargeNegativeMinusSmallNegative;
+begin
+  Test('-B12BE1E20098', '-148F3137CED396', '13DE0555ECD2FE');
+end;
+
+procedure TBigIntDifferenceTestCase.TestSmallNegativeMinusLargeNegative;
+begin
+  Test('-AF3FF1EC626908C', '-18295', '-AF3FF1EC6250DF7');
+end;
+
+procedure TBigIntDifferenceTestCase.TestNegativeMinusPositive;
+begin
+  Test('-E493506B19', '20508ED255', '-104E3DF3D6E');
+end;
+
+procedure TBigIntDifferenceTestCase.TestPositiveMinusNegative;
+begin
+  Test('114EEC66851AFD98', '-100AA4308C5249FBBFADEB89CD6A7D9CC', '100AA4308C5249FBC0C2DA5035BC2D764');
+end;
+
+procedure TBigIntDifferenceTestCase.TestLargeOperands;
+begin
+  Test('1069FC8EA3D99C39E1C07AA078B77B5CC679CB448563345A878C603D32FB0F8FEFE02AD9574A5EB6',
+    '1069FC8EA3D99C39E1C07AA078B77B5CC679CB448563345A878C603D32FB0F8FEFE023C522B87F8C',
+    '7143491DF2A');
+end;
+
+procedure TBigIntDifferenceTestCase.TestPositiveMinusZero;
+begin
+  Test('8ABB000911FF5A000333', '0', '8ABB000911FF5A000333');
+end;
+
+procedure TBigIntDifferenceTestCase.TestZeroMinusPositive;
+begin
+  Test('0', '243961982DDD64F81B2', '-243961982DDD64F81B2');
+end;
+
+procedure TBigIntDifferenceTestCase.TestZero;
+begin
+  Test('0', '0', '0');
+end;
+
+procedure TBigIntDifferenceTestCase.TestDifferenceShorterLeft;
+begin
+  Test('3FFFF', '9000911040221', '-9000911000222');
+end;
+
+procedure TBigIntDifferenceTestCase.TestDifferenceShorterRight;
+begin
+  Test('900094100F443', '3000EEEE', '9000911000555');
+end;
+
+procedure TBigIntDifferenceTestCase.TestDifferenceEndsInCarry;
+begin
+  Test('400054400021F', 'FFFFFFFF', '4000444000220');
+end;
+
+procedure TBigIntDifferenceTestCase.TestDifferenceEndsInCarryLosingDigit;
+begin
+  Test('133320055', '99990000', '99990055');
+end;
+
+procedure TBigIntDifferenceTestCase.TestDifferenceMultiCarry;
+begin
+  Test('1000000000000000000000000FFFFFFFE', 'FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF', 'FFFFFFFF');
+end;
+
+{ TBigIntProductTestCase }
+
+procedure TBigIntProductTestCase.Test(const AHexValueLeft, AHexValueRight, AHexValueProduct: string);
+var
+  a, b, s: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueLeft);
+  b := TBigInt.FromHexadecimalString(AHexValueRight);
+  s := TBigInt.FromHexadecimalString(AHexValueProduct);
+  AssertTrue('Product of BigInt from ''' + AHexValueLeft + ''' and from ''' + AHexValueRight
+    + ''' was not equal to the BigInt from ''' + AHexValueProduct + '''.',
+    s = a * b);
+end;
+
+procedure TBigIntProductTestCase.TestShort;
+begin
+  Test('9', 'B', '63');
+  Test('29E531A', 'DF5F299', '248E3915103E8A');
+end;
+
+procedure TBigIntProductTestCase.TestLongPositiveTimesPositive;
+begin
+  Test('FFFFFFFF', 'FFFFFFFF', 'FFFFFFFE00000001');
+  Test('4873B23699D07741F1544862221B1966AA84512', '4DE0013', '160A322C656DEB1DD6D721D36E35F2E29B4D2A18192056');
+  Test('74FD3E6988116762', '22DB271AFC4941', 'FEDC8CD51DEE46BE83C283B5E31E2');
+end;
+
+procedure TBigIntProductTestCase.TestLongPositiveTimesNegative;
+begin
+  Test('23401834190D12FF3210F0B0129123AA', '-A4C0530234', '-16AF8B019CA1436BBFD1F1FB08494FFC9EF7E09288');
+end;
+
+procedure TBigIntProductTestCase.TestLongNegativeTimesPositive;
+begin
+  Test('-3ACB78882923810', 'F490B8022A4BCBFF34E01', '-382B2B9851BC93CB0308B502C3B036D71810');
+end;
+
+procedure TBigIntProductTestCase.TestLongNegativeTimesNegative;
+begin
+  Test('-554923FB201', '-9834FDC032', '32B514C1BA1E774EE8432');
+end;
+
+procedure TBigIntProductTestCase.TestZeroTimesPositive;
+begin
+  Test('0', '1AF5D0039B888AC00F299', '0');
+end;
+
+procedure TBigIntProductTestCase.TestZeroTimesNegative;
+begin
+  Test('0', '-1AF5D0039B888AC00F299', '0');
+end;
+
+procedure TBigIntProductTestCase.TestZero;
+begin
+  Test('0', '0', '0');
+end;
+
+{ TBigIntShiftLeftTestCase }
+
+procedure TBigIntShiftLeftTestCase.Test(const AHexValueOperand: string; const AShift: Integer; const AHexValueResult: string);
+var
+  a, s: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueOperand);
+  s := TBigInt.FromHexadecimalString(AHexValueResult);
+  AssertTrue('BigInt from hexadecimal string ''' + AHexValueOperand + ''' shifted left by ''' + IntToStr(AShift)
+    + ''' was not equal to BigInt from hexadecimal string ''' + AHexValueResult + '''.',
+    s = (a << AShift));
+end;
+
+procedure TBigIntShiftLeftTestCase.TestShort;
+begin
+  // BIN 101
+  // BIN 101000
+  Test('5', 3, '28');
+end;
+
+procedure TBigIntShiftLeftTestCase.TestShortWithCarry;
+begin
+  // BIN 1110 1101 0111 0001 1010 1010
+  // BIN 1 1101 1010 1110 0011 0101 0100 0000 0000 0000
+  Test('ED71AA', 13, '1DAE354000');
+end;
+
+procedure TBigIntShiftLeftTestCase.TestLongWithCarry;
+begin
+  // BIN 1 0011 0000 1011 0010 0011 1110 0100 1100 1111 1001 0101 0100 0100 0101
+  // BIN 10 0110 0001 0110 0100 0111 1100 1001 1001 1111 0010 1010 1000 1000 1010 0000 0000
+  Test('130B23E4CF95445', 9, '261647C99F2A88A00');
+end;
+
+procedure TBigIntShiftLeftTestCase.TestLongWithMultiDigitCarry;
+begin
+  Test('C99F12A735A3B83901BF92011', 140, 'C99F12A735A3B83901BF9201100000000000000000000000000000000000');
+end;
+
+procedure TBigIntShiftLeftTestCase.TestWithAlignedDigits;
+begin
+  // Shifts the left operand by a multiple of the length of full TBigInt digits, so the digits will be shifted, but not
+  // changed.
+  Test('10F0F39C5E', 32 * 4, '10F0F39C5E00000000000000000000000000000000');
+end;
+
+procedure TBigIntShiftLeftTestCase.TestZero;
+begin
+  Test('0', 119, '0');
+end;
+
+{ TBigIntEqualityTestCase }
+
+procedure TBigIntEqualityTestCase.TestEqual(const AValue: Int64);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromInt64(AValue);
+  b := TBigInt.FromInt64(AValue);
+  AssertTrue('Two BigInt from ''' + IntToStr(AValue) + ''' were not equal.', a = b);
+  AssertFalse('Two BigInt from ''' + IntToStr(AValue) + ''' were un-equal.', a <> b);
+end;
+
+procedure TBigIntEqualityTestCase.TestEqualHex(const AHexValue: string);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValue);
+  b := TBigInt.FromHexadecimalString(AHexValue);
+  AssertTrue('Two BigInt from ''' + AHexValue + ''' were not equal.', a = b);
+  AssertFalse('Two BigInt from ''' + AHexValue + ''' were un-equal.', a <> b);
+end;
+
+procedure TBigIntEqualityTestCase.TestNotEqualHex(const AHexValueLeft, AHexValueRight: string);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueLeft);
+  b := TBigInt.FromHexadecimalString(AHexValueRight);
+  AssertTrue('BigInt from ''' + AHexValueLeft + ''' and from ''' + AHexValueRight + ''' were not un-equal.',
+    a <> b);
+  AssertFalse('BigInt from ''' + AHexValueLeft + ''' and from ''' + AHexValueRight + ''' were equal.',
+    a = b);
+end;
+
+procedure TBigIntEqualityTestCase.TestShortEqual;
+begin
+  TestEqual(14);
+  TestEqualHex('8000111000111');
+end;
+
+procedure TBigIntEqualityTestCase.TestLongEqual;
+begin
+  TestEqualHex('5AB60FF292A014BF1DD0A');
+end;
+
+procedure TBigIntEqualityTestCase.TestZeroEqual;
+begin
+  TestEqual(0);
+end;
+
+procedure TBigIntEqualityTestCase.TestShortNotEqual;
+begin
+  TestNotEqualHex('9FF99920', '100');
+  TestNotEqualHex('20001110002111', '70001110007111');
+end;
+
+procedure TBigIntEqualityTestCase.TestLongNotEqualSign;
+begin
+  TestNotEqualHex('48843AB320FF0041123A', '-48843AB320FF0041123A');
+  TestNotEqualHex('-B13F79D842A30957DD09523667', 'B13F79D842A30957DD09523667');
+end;
+
+procedure TBigIntEqualityTestCase.TestZeroNotEqual;
+begin
+  TestNotEqualHex('0', 'F');
+  TestNotEqualHex('568F7', '0');
+end;
+
+{ TBigIntComparisonTestCase }
+
+procedure TBigIntComparisonTestCase.TestLessThan(const AHexValueLeft, AHexValueRight: string);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueLeft);
+  b := TBigInt.FromHexadecimalString(AHexValueRight);
+  AssertTrue('BigInt from ''' + AHexValueLeft + ''' was greater than BigInt from ''' + AHexValueRight + '''.',
+    a < b);
+  AssertTrue('BigInt from ''' + AHexValueLeft + ''' was greater or equal than BigInt from ''' + AHexValueRight + '''.',
+    a <= b);
+  AssertFalse('BigInt from ''' + AHexValueLeft + ''' was greater than BigInt from ''' + AHexValueRight + '''.',
+    a > b);
+  AssertFalse('BigInt from ''' + AHexValueLeft + ''' was greater or equal than BigInt from ''' + AHexValueRight + '''.',
+    a >= b);
+end;
+
+procedure TBigIntComparisonTestCase.TestGreaterThan(const AHexValueLeft, AHexValueRight: string);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValueLeft);
+  b := TBigInt.FromHexadecimalString(AHexValueRight);
+  AssertFalse('BigInt from ''' + AHexValueLeft + ''' was less than BigInt from ''' + AHexValueRight + '''.',
+    a < b);
+  AssertFalse('BigInt from ''' + AHexValueLeft + ''' was less or equal than BigInt from ''' + AHexValueRight + '''.',
+    a <= b);
+  AssertTrue('BigInt from ''' + AHexValueLeft + ''' was less than BigInt from ''' + AHexValueRight + '''.',
+    a > b);
+  AssertTrue('BigInt from ''' + AHexValueLeft + ''' was less or equal than BigInt from ''' + AHexValueRight + '''.',
+    a >= b);
+end;
+
+procedure TBigIntComparisonTestCase.TestEqual(const AHexValue: string);
+var
+  a, b: TBigInt;
+begin
+  a := TBigInt.FromHexadecimalString(AHexValue);
+  b := TBigInt.FromHexadecimalString(AHexValue);
+  AssertFalse('First BigInt from ''' + AHexValue + ''' was less than the second.',
+    a < b);
+  AssertTrue('First BigInt from ''' + AHexValue + ''' was greater than the second.',
+    a <= b);
+  AssertFalse('First BigInt from ''' + AHexValue + ''' was greater than the second.',
+    a > b);
+  AssertTrue('First BigInt from ''' + AHexValue + ''' was less than the second.',
+    a >= b);
+end;
+
+procedure TBigIntComparisonTestCase.TestLessSameLength;
+begin
+  TestLessThan('104234FF2B29100C012', '234867AB261F1003429103C');
+  TestLessThan('-9812FB2964AC7632865238BBD3', '294625DF51B2A842582244C18163490');
+  TestLessThan('-12834653A2856DF8', '-A946C2BF89401B8');
+  TestLessThan('-2F846', '0');
+  TestLessThan('0', 'FF67B');
+end;
+
+procedure TBigIntComparisonTestCase.TestLessShorterLeft;
+begin
+  TestLessThan('34FF2B29100C012', '234867AB261F1003429103C');
+  TestLessThan('0', 'BFF008112BA00012');
+end;
+
+procedure TBigIntComparisonTestCase.TestLessNegativeShorterRight;
+begin
+  TestLessThan('-9B72844AC', '1F3486B');
+  TestLessThan('-BB884F022111190', '0');
+end;
+
+procedure TBigIntComparisonTestCase.TestGreaterSameLength;
+begin
+  TestGreaterThan('B042104234FF2B29100C012', '23867AB261F1003429103C');
+  TestGreaterThan('1294B77', '-611F3B93');
+  TestGreaterThan('-782163498326593562D548AAF715B4DC9143E42C68F39A29BB2', '-FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF');
+  TestGreaterThan('783BDA0', '0');
+  TestGreaterThan('0', '-99B6D');
+end;
+
+procedure TBigIntComparisonTestCase.TestGreaterShorterRight;
+begin
+  TestGreaterThan('102343234423FF67278B11ADD345F6BB21232C9129100C012', '234867AB261F1003429103C');
+  TestGreaterThan('249D00F63BBAA8668B23', '0');
+end;
+
+procedure TBigIntComparisonTestCase.TestGreaterNegativeShorterLeft;
+begin
+  TestGreaterThan('81F2A7B78B812', '-23490D7F19F247F91A365B1893BB701');
+  TestGreaterThan('0', '-80F88242B34127');
+end;
+
+procedure TBigIntComparisonTestCase.TestEqualPositive;
+begin
+  TestEqual('A44B80191059CA123318921A219BB');
+end;
+
+procedure TBigIntComparisonTestCase.TestEqualNegative;
+begin
+  TestEqual('-28912798DD1246DAC9FB4269908012D496896812FF3A8D071B32');
+end;
+
+procedure TBigIntComparisonTestCase.TestEqualZero;
+begin
+  TestEqual('0');
+end;
+
+initialization
+
+  RegisterTest(TBigIntSignTestCase);
+  RegisterTest(TBigIntMostSignificantBitIndexTestCase);
+  RegisterTest(TBigIntFromInt64TestCase);
+  RegisterTest(TBigIntFromHexTestCase);
+  RegisterTest(TBigIntFromBinTestCase);
+  RegisterTest(TBigIntUnaryMinusTestCase);
+  RegisterTest(TBigIntSumTestCase);
+  RegisterTest(TBigIntDifferenceTestCase);
+  RegisterTest(TBigIntProductTestCase);
+  RegisterTest(TBigIntShiftLeftTestCase);
+  RegisterTest(TBigIntEqualityTestCase);
+  RegisterTest(TBigIntComparisonTestCase);
+end.
+

tests/UNeverTellMeTheOddsTestCases.pas

@@ -33,6 +33,7 @@ type
     function CreateSolver: ISolver; override;
   published
     procedure TestPart1;
+    procedure TestPart2;
   end;
 
   { TNeverTellMeTheOddsExampleTestCase }
@@ -42,6 +43,7 @@ type
     function CreateSolver: ISolver; override;
   published
     procedure TestPart1;
+    procedure TestPart2;
   end;
 
   { TNeverTellMeTheOddsTestCase }
@@ -77,6 +79,11 @@ begin
   AssertEquals(15107, FSolver.GetResultPart1);
 end;
 
+procedure TNeverTellMeTheOddsFullDataTestCase.TestPart2;
+begin
+  AssertEquals(-1, FSolver.GetResultPart2);
+end;
+
 { TNeverTellMeTheOddsExampleTestCase }
 
 function TNeverTellMeTheOddsExampleTestCase.CreateSolver: ISolver;
@@ -89,6 +96,11 @@ begin
   AssertEquals(2, FSolver.GetResultPart1);
 end;
 
+procedure TNeverTellMeTheOddsExampleTestCase.TestPart2;
+begin
+  AssertEquals(47, FSolver.GetResultPart2);
+end;
+
 { TNeverTellMeTheOddsTestCase }
 
 function TNeverTellMeTheOddsTestCase.CreateSolver: ISolver;

tests/UPolynomialRootsTestCases.pas

@@ -0,0 +1,72 @@
+{
+  Solutions to the Advent Of Code.
+  Copyright (C) 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 UPolynomialRootsTestCases;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, UPolynomial, UPolynomialRoots, UBigInt;
+
+type
+
+  { TPolynomialRootsTestCase }
+
+  TPolynomialRootsTestCase = class(TTestCase)
+  protected
+    FRootIsolation: TRootIsolation;
+    procedure SetUp; override;
+    procedure TearDown; override;
+  published
+    procedure TestBisectionRootIsolation;
+  end;
+
+implementation
+
+{ TPolynomialRootsTestCase }
+
+procedure TPolynomialRootsTestCase.SetUp;
+begin
+  inherited SetUp;
+  FRootIsolation := TRootIsolation.Create;
+end;
+
+procedure TPolynomialRootsTestCase.TearDown;
+begin
+  FRootIsolation.Free;
+  inherited TearDown;
+end;
+
+procedure TPolynomialRootsTestCase.TestBisectionRootIsolation;
+var
+  a: TBigIntPolynomial;
+  r: Int64;
+begin
+  // y = 3 * (x - 34000) * (x - 23017) * (x - 5) * (x^2 - 19) * (x + 112)
+  //   = 3 * x^6 - 170730 * x^5 + 2329429920 * x^4 + 251300082690 * x^3 - 1270471872603 * x^2 + 4774763204640 * x - 24979889760000
+  a := TBigIntPolynomial.Create([-24979889760000, 4774763204640, -1270471872603, 251300082690, 2329429920, -170730, 3]);
+  r := FRootIsolation.Bisect(a);
+  AssertEquals(0, r);
+end;
+
+initialization
+
+  RegisterTest(TPolynomialRootsTestCase);
+end.
+

tests/UPolynomialTestCases.pas

@@ -0,0 +1,117 @@
+{
+  Solutions to the Advent Of Code.
+  Copyright (C) 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 UPolynomialTestCases;
+
+{$mode ObjFPC}{$H+}
+
+interface
+
+uses
+  Classes, SysUtils, fpcunit, testregistry, UPolynomial, UBigInt;
+
+type
+
+  { TBigIntPolynomialTestCase }
+
+  TBigIntPolynomialTestCase = class(TTestCase)
+  private
+    procedure TestCreateWithDegree(const ACoefficients: array of TBigInt; const ADegree: Integer);
+  published
+    procedure TestCreate;
+    procedure TestCreateDegreeOne;
+    procedure TestCreateDegreeZero;
+    procedure TestCreateDegreeZeroTrim;
+    procedure TestEqual;
+    procedure TestUnequalSameLength;
+    procedure TestUnequalDifferentLength;
+    procedure TestTrimLeadingZeros;
+  end;
+
+implementation
+
+{ TBigIntPolynomialTestCase }
+
+procedure TBigIntPolynomialTestCase.TestCreateWithDegree(const ACoefficients: array of TBigInt; const ADegree: Integer);
+var
+  a: TBigIntPolynomial;
+begin
+  a := TBigIntPolynomial.Create(ACoefficients);
+  AssertEquals('Degree of created polynomial incorrect.', ADegree, a.Degree);
+end;
+
+procedure TBigIntPolynomialTestCase.TestCreate;
+begin
+  TestCreateWithDegree([992123, 7, 20, 4550022], 3);
+end;
+
+procedure TBigIntPolynomialTestCase.TestCreateDegreeOne;
+begin
+  TestCreateWithDegree([4007], 0);
+end;
+
+procedure TBigIntPolynomialTestCase.TestCreateDegreeZero;
+begin
+  TestCreateWithDegree([], -1);
+end;
+
+procedure TBigIntPolynomialTestCase.TestCreateDegreeZeroTrim;
+begin
+  TestCreateWithDegree([0], -1);
+end;
+
+procedure TBigIntPolynomialTestCase.TestEqual;
+var
+  a, b: TBigIntPolynomial;
+begin
+  a := TBigIntPolynomial.Create([10, 7, 5, 1034]);
+  b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
+  AssertTrue('Polynomials are not equal.', a = b);
+end;
+
+procedure TBigIntPolynomialTestCase.TestUnequalSameLength;
+var
+  a, b: TBigIntPolynomial;
+begin
+  a := TBigIntPolynomial.Create([103, 7, 5, 10]);
+  b := TBigIntPolynomial.Create([1034, 7, 5, 10]);
+  AssertTrue('Polynomials are equal.', a <> b);
+end;
+
+procedure TBigIntPolynomialTestCase.TestUnequalDifferentLength;
+var
+  a, b: TBigIntPolynomial;
+begin
+  a := TBigIntPolynomial.Create([40000, 10, 7, 5, 1034]);
+  b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
+  AssertTrue('Polynomials are equal.', a <> b);
+end;
+
+procedure TBigIntPolynomialTestCase.TestTrimLeadingZeros;
+var
+  a, b: TBigIntPolynomial;
+begin
+  a := TBigIntPolynomial.Create([10, 7, 5, 1034, 0, 0]);
+  b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
+  AssertTrue('Polynomials are not equal.', a = b);
+end;
+
+initialization
+
+  RegisterTest(TBigIntPolynomialTestCase);
+end.
+