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

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

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

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

unit UNeverTellMeTheOdds;

{$mode ObjFPC}{$H+}

interface

uses
  Classes, SysUtils, Generics.Collections, Math, USolver, UNumberTheory, UBigInt, UPolynomial, UPolynomialRoots;

type

  { THailstone }

  THailstone = class
  public
    P0, P1, P2: Int64;
    V0, V1, V2: Integer;
    constructor Create(const ALine: string);
    constructor Create;
  end;

  THailstones = specialize TObjectList<THailstone>;

  TInt64Array = array of Int64;

  { TNeverTellMeTheOdds }

  TNeverTellMeTheOdds = class(TSolver)
  private
    FMin, FMax: Int64;
    FHailstones: THailstones;
    FA: array[0..10] of TBigInt;
    FH: array[0..6] of TBigInt;
    function AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
    function FindRockThrow(const AIndex0, AIndex1, AIndex2: Integer): Int64;
    procedure CalcCollisionPolynomials(constref AHailstone0, AHailstone1, AHailstone2: THailstone; out OPolynomial0,
      OPolynomial1: TBigIntPolynomial);
    function CalcRockThrowCollisionOptions(constref AHailstone0, AHailstone1, AHailstone2: THailstone): TInt64Array;
    function ValidateRockThrow(constref AHailstone0, AHailstone1, AHailstone2: THailstone; const AT0, AT1: Int64):
      Int64;
  public
    constructor Create(const AMin: Int64 = 200000000000000; const AMax: Int64 = 400000000000000);
    destructor Destroy; override;
    procedure ProcessDataLine(const ALine: string); override;
    procedure Finish; override;
    function GetDataFileName: string; override;
    function GetPuzzleName: string; override;
  end;

implementation

{ THailstone }

constructor THailstone.Create(const ALine: string);
var
  split: TStringArray;
begin
  split := ALine.Split([',', '@']);
  P0 := StrToInt64(Trim(split[0]));
  P1 := StrToInt64(Trim(split[1]));
  P2 := StrToInt64(Trim(split[2]));
  V0 := StrToInt(Trim(split[3]));
  V1 := StrToInt(Trim(split[4]));
  V2 := StrToInt(Trim(split[5]));
end;

constructor THailstone.Create;
begin

end;

{ TNeverTellMeTheOdds }

function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
var
  m1, m2, x, y: Double;
begin
  Result := False;
  m1 := AHailstone1.V1 / AHailstone1.V0;
  m2 := AHailstone2.V1 / AHailstone2.V0;
  if m1 <> m2 then
  begin
    x := (AHailstone2.P1 - m2 * AHailstone2.P0
      - AHailstone1.P1 + m1 * AHailstone1.P0)
      / (m1 - m2);
    if (FMin <= x) and (x <= FMax)
    and (x * Sign(AHailstone1.V0) >= AHailstone1.P0 * Sign(AHailstone1.V0))
    and (x * Sign(AHailstone2.V0) >= AHailstone2.P0 * Sign(AHailstone2.V0))
    then
    begin
      y := m1 * (x - AHailstone1.P0) + AHailstone1.P1;
      if (FMin <= y) and (y <= FMax) then
        Result := True
    end;
  end;
end;

function TNeverTellMeTheOdds.FindRockThrow(const AIndex0, AIndex1, AIndex2: Integer): Int64;
var
  t0, t1: TInt64Array;
  i, j: Int64;
begin
  t0 := CalcRockThrowCollisionOptions(FHailstones[AIndex0], FHailstones[AIndex1], FHailstones[AIndex2]);
  t1 := CalcRockThrowCollisionOptions(FHailstones[AIndex1], FHailstones[AIndex0], FHailstones[AIndex2]);

  Result := 0;
  for i in t0 do
  begin
    for j in t1 do
    begin
      Result := ValidateRockThrow(FHailstones[AIndex0], FHailstones[AIndex1], FHailstones[AIndex2], i, j);
      if Result > 0 then
        Break;
    end;
    if Result > 0 then
      Break;
  end;
end;

procedure TNeverTellMeTheOdds.CalcCollisionPolynomials(constref AHailstone0, AHailstone1, AHailstone2: THailstone; out
  OPolynomial0, OPolynomial1: TBigIntPolynomial);
var
  k: array[0..139] of TBigInt;
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) * (a_1 - V_20 * t_2) - (t_2 - t_0) * (a_2 - 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:  f_2 * t_1^2 + f_1 * t_1 - f_0 = 0
  // Then, with 4 and 3 in 1 and many substitutions (see constants k below, now independent of t_0), the equation
  //   5:  0 = p_0(t_0) + p_1(t_0) * sqrt(p_2(t_0))
  // can be constructed, where p_0, p_1, and p_2 are polynomials in t_0. Since we are searching for an integer solution,
  // we assume that there is an integer t_0 that is a root of both p_0 and p_1, which would solve the equation.

  // Subsitutions depending on t_0:
  //   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_0 = b_1 * d_0 + b_0 * e_0
  //   f_1 = c_0 * e_0 + c_1 * d_0 - b_0 * e_1 - b_1 * d_1
  //   f_2 = c_0 * e_1 + c_1 * d_1

  // Calculations for equation 5 (4 and 3 in 1).
  // 1:  0 = (t_1 - t_0) * (a_1 - V_20 * t_2) - (t_2 - t_0) * (a_2 - V_10 * t_1)
  // 3:  (e_0 + e_1 * t_1) * t_2 = (d_0 + d_1 * t_1)
  // 0 = (t_1 - t_0) * (a_1 - V_20 * t_2) - (t_2 - t_0) * (a_2 - V_10 * t_1)
  //   = (t_1 - t_0) * (a_1 * (e_0 + e_1 * t_1) - V_20 * (e_0 + e_1 * t_1) * t_2) - ((e_0 + e_1 * t_1) * t_2 - (e_0 + e_1 * t_1) * t_0) * (a_2 - V_10 * t_1)
  //   = (t_1 - t_0) * (a_1 * (e_0 + e_1 * t_1) - V_20 * (d_0 + d_1 * t_1)) - ((d_0 + d_1 * t_1) - (e_0 + e_1 * t_1) * t_0) * (a_2 - V_10 * t_1)
  //   = (t_1 - t_0) * (a_1 * e_0 + a_1 * e_1 * t_1 - V_20 * d_0 - V_20 * d_1 * t_1) - (d_0 + d_1 * t_1 - e_0 * t_0 - e_1 * t_1 * t_0) * (a_2 - V_10 * t_1)
  //   = (a_1 * e_1 - V_20 * d_1) * t_1^2 + (a_1 * e_0 - V_20 * d_0 - t_0 * (a_1 * e_1  - V_20 * d_1)) * t_1 - t_0 * (a_1 * e_0 - V_20 * d_0)
  //     - ( - V_10 * (d_1 - e_1 * t_0) * t_1^2 + ((d_1 - e_1 * t_0) * a_2 - V_10 * (d_0 - e_0 * t_0)) * t_1 + (d_0 - e_0 * t_0) * a_2)
  //   = (a_1 * e_1 - V_20 * d_1 + V_10 * (d_1 - e_1 * t_0)) * t_1^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)) * t_1
  //     + t_0 * (V_20 * d_0 - a_1 * e_0) + (e_0 * t_0 - d_0) * a_2
  // Inserting 4, solved for t_0:  t_1 = - f_1 / (2 * f_2) + sqrt((f_1 / (2 * f_2))^2 + f_0 / f_2)
  //   = (a_1 * e_1 - V_20 * d_1 + V_10 * (d_1 - e_1 * t_0)) * (f_1^2 + 2 * f_0 * f_2 - f_1 * sqrt(f_1^2 + 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^2 + 4 * f_0 * f_2))
  //     + t_0 * (V_20 * d_0 - a_1 * e_0) * 2 * f_2^2 + (e_0 * t_0 - d_0) * a_2 * 2 * f_2^2

  // a_1 = V_00 * t_0 + k_0
  // a_2 = V_00 * t_0 + k_1
  // b_0 = k_2 * t_0
  // b_1 = k_10 * t_0 + k_3
  // c_0 = k_11 * t_0 + k_4
  // d_0 = k_5 * t_0
  // d_1 = k_12 * t_0 + k_6
  // e_0 = k_13 * t_0 + k_7
  // f_2 = (k_11 * t_0 + k_4) * k_9 + k_8 * (k_12 * t_0 + k_6)
  //     = (k_11 * k_9 + k_8 * k_12) * t_0 + k_4 * k_9 + k_8 * k_6
  //     = FH_0 * t_0 + FH_1
  // f_1 = (k_11 * t_0 + k_4) * (k_13 * t_0 + k_7) + k_8 * k_5 * t_0 - k_2 * t_0 * k_9 - (k_10 * t_0 + k_3) * (k_12 * t_0 + k_6)
  //     = (k_11 * k_13 - k_10 * k_12) * t_0^2 + (k_11 * k_7 + k_4 * k_12 + k_8 * k_5 - k_2 * k_9 - k_10 * k_6 - k_3 * k_12) * t_0 + k_4 * k_7 - k_3 * k_6
  //     = FH_2 * t_0^2 + FH_3 * t_0 + FH_4
  // f_0 = (k_10 * t_0 + k_3) * k_5 * t_0 + k_2 * t_0 * (k_13 * t_0 + k_7)
  //     = (k_10 * k_5 + k_2 * k_13) * t_0^2 + (k_3 * k_5 + k_2 * k_7) * t_0
  //     = FH_5 * t_0^2 + FH_6 * t_0

  k[0] := AHailstone0.P0 - AHailstone2.P0;
  k[1] := AHailstone0.P0 - AHailstone1.P0;
  k[2] := AHailstone1.P1 - AHailstone2.P1;
  k[3] := AHailstone0.P1 - AHailstone1.P1;
  k[4] := AHailstone0.P1 - AHailstone2.P1;
  k[5] := AHailstone2.P2 - AHailstone1.P2;
  k[6] := AHailstone0.P2 - AHailstone2.P2;
  k[7] := AHailstone0.P2 - AHailstone1.P2;
  k[8] := AHailstone1.V1 - AHailstone2.V1;
  k[9] := AHailstone2.V2 - AHailstone1.V2;
  k[10] := AHailstone0.V1 - AHailstone2.V1;
  k[11] := AHailstone0.V1 - AHailstone1.V1;
  k[12] := AHailstone0.V2 - AHailstone1.V2;
  k[13] := AHailstone0.V2 - AHailstone2.V2;

  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];

  // Additional substitutions.
  //   a_1 * k_9 - V_20 * d_1
  //       = (V_00 * t_0 + k_0) * k_9 - V_20 * (k_12 * t_0 + k_6)
  //       = (V_00 * k_9 - V_20 * k_12) * t_0 + k_0 * k_9 - V_20 * k_6
  //       = k_14 * t_0 + k_15
  //   d_1 - k_9 * t_0
  //       = k_12 * t_0 + k_6 - k_9 * t_0
  //       = (k_12 - k_9) * t_0 + k_6
  //   a_1 * e_0 - V_20 * d_0
  //       = (V_00 * t_0 + k_0) * (k_13 * t_0 + k_7) - V_20 * k_5 * t_0
  //       = V_00 * k_13 * t_0^2 + (V_00 * k_7 + k_0 * k_13 - V_20 * k_5) * t_0 + k_0 * k_7
  //       = k_16 * t_0^2 + k_17 * t_0 + k_18
  //   d_0 - e_0 * t_0
  //       = k_5 * t_0 - k_13 * t_0^2 - k_7 * t_0
  //       = - k_13 * t_0^2 + k_19 * t_0
  //   f_1^2
  //       = (FH_2 * t_0^2 + FH_3 * t_0 + FH_4)^2
  //       = FH_2^2 * t_0^4 + FH_3^2 * t_0^2 + FH_4^2 + 2 * FH_2 * t_0^2 * FH_3 * t_0 + 2 * FH_2 * t_0^2 * FH_4 + 2 * FH_3 * t_0 * FH_4
  //       = FH_2^2 * t_0^4 + 2 * FH_2 * FH_3 * t_0^3 + (FH_3^2 + 2 * FH_2 * FH_4) * t_0^2 + 2 * FH_3 * FH_4 * t_0 + FH_4^2
  //       = FH_2^2 * t_0^4 + k_20 * t_0^3 + k_22 * t_0^2 + k_23 * t_0 + FH_4^2
  //   f_2^2
  //       = (FH_0 * t_0 + FH_1)^2
  //       = FH_0^2 * t_0^2 + 2 * FH_0 * FH_1 * t_0 + FH_1^2
  //       = FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2
  //   f_0 * f_2
  //       = (FH_5 * t_0^2 + FH_6 * t_0) * (FH_0 * t_0 + FH_1)
  //       = FH_5 * FH_0 * t_0^3 + (FH_5 * FH_1 + FH_6 * FH_0) * t_0^2 + FH_6 * FH_1 * t_0
  //       = k_126 * t_0^3 + k_127 * t_0^2 + k_128 * t_0
  //   f_1^2 + 4 * f_0 * f_2
  //       = FH_2^2 * t_0^4 + k_20 * t_0^3 + k_22 * t_0^2 + k_23 * t_0 + FH_4^2 + 4 * (k_126 * t_0^3 + k_127 * t_0^2 + k_128 * t_0)
  //       = k_31 * t_0^4 + k_132 * t_0^3 + k_133 * t_0^2 + k_134 * t_0 + k_58
  //   f_1^2 + 2 * f_0 * f_2
  //       = FH_2^2 * t_0^4 + k_20 * t_0^3 + k_22 * t_0^2 + k_23 * t_0 + FH_4^2 + 2 * (k_126 * t_0^3 + k_127 * t_0^2 + k_128 * t_0)
  //       = k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58

  k[14] := AHailstone0.V0 * k[9] - AHailstone2.V0 * k[12];
  k[15] := k[0] * k[9] - AHailstone2.V0 * k[6];
  k[16] := AHailstone0.V0 * k[13];
  k[17] := AHailstone0.V0 * k[7] + k[0] * k[13] - AHailstone2.V0 * 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];
  k[126] := FH[5] * FH[0];
  k[127] := FH[5] * FH[1] + FH[6] * FH[0];
  k[128] := FH[6] * FH[1];
  k[28] := FH[1] * FH[1];
  k[31] := FH[2] * FH[2];
  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] + AHailstone1.V0 * (k[12] - k[9]);
  k[41] := k[15] + AHailstone1.V0 * k[6];
  k[42] := k[16] - k[14] - AHailstone1.V0 * k[13] - (k[12] - k[9]) * AHailstone0.V0;
  k[43] := k[17] - k[15] + AHailstone1.V0 * k[19] - (k[12] - k[9]) * k[1] - k[6] * AHailstone0.V0;
  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] * AHailstone0.V0 - k[16];
  k[62] := 2 * k[25] * k[61];
  k[63] := k[13] * k[1] - k[19] * AHailstone0.V0 - 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];
  // Unused, they are part of the polynomial inside the square root.
  //k[132] := k[20] + 4 * k[126];
  //k[133] := k[22] + 4 * k[127];
  //k[134] := k[23] + 4 * k[128];

  // Continuing calculations for equation 5.
  // 0 = (k_14 * t_0 + k_15 + V_10 * ((k_12 - k_9) * t_0 + k_6)) * (k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58 -+ f_1 * sqrt(f_1^2 + 4 * f_0 * f_2))
  //     + (k_16 * t_0^2 + 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^2 - k_19 * t_0)) * (- f_1 * f_2 +- f_2 * sqrt(f_1^2 + 4 * f_0 * f_2))
  //     - 2 * t_0 * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2) + 2 * (k_13 * t_0^2 - k_19 * t_0) * a_2 * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2)
  // 0 = (k_40 * t_0 + k_41) * (k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58 -+ f_1 * sqrt(f_1^2 + 4 * f_0 * f_2))
  //     + ((k_16 - k_14 - V_10 * k_13 - (k_12 - k_9) * V_00) * t_0^2 + (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^2 + 4 * f_0 * f_2))
  //     - 2 * t_0 * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2) + 2 * (k_13 * t_0^2 - k_19 * t_0) * a_2 * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2)
  // 0 = (k_40 * t_0 + k_41) * (k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58)
  //     -+ (k_40 * t_0 + k_41) * f_1 * sqrt(f_1^2 + 4 * f_0 * f_2)
  //     + (k_42 * t_0^2 + k_43 * t_0 + k_44) * (- f_1 * f_2 +- f_2 * sqrt(f_1^2 + 4 * f_0 * f_2))
  //     - 2 * t_0 * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2) + 2 * (k_13 * t_0^2 - k_19 * t_0) * a_2 * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2)
  // 0 = (k_40 * t_0 + k_41) * (k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58)
  //     -+ (k_40 * t_0 + k_41) * f_1 * sqrt(f_1^2 + 4 * f_0 * f_2)
  //     - (k_42 * t_0^2 + k_43 * t_0 + k_44) * f_1 * f_2
  //     +- (k_42 * t_0^2 + k_43 * t_0 + k_44) * f_2 * sqrt(f_1^2 + 4 * f_0 * f_2)
  //     - 2 * t_0 * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2) + 2 * (k_13 * t_0^2 - k_19 * t_0) * a_2 * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2)
  // 0 = +- ((k_42 * t_0^2 + k_43 * t_0 + k_44) * f_2 - (k_40 * t_0 + k_41) * f_1) * sqrt(f_1^2 + 4 * f_0 * f_2)
  //     + (k_40 * t_0 + k_41) * (k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58)
  //     - (k_42 * t_0^2 + k_43 * t_0 + k_44) * f_1 * f_2
  //     - 2 * t_0 * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2) + 2 * (k_13 * t_0^2 - k_19 * t_0) * a_2 * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2)
  // 0 = +- ((k_42 * t_0^2 + k_43 * t_0 + k_44) * (FH_0 * t_0 + FH_1) - (k_40 * t_0 + k_41) * (FH_2 * t_0^2 + FH_3 * t_0 + FH_4)) * sqrt(f_1^2 + 4 * f_0 * f_2)
  //     + (k_40 * t_0 + k_41) * (k_31 * t_0^4 + k_137 * t_0^3 + k_138 * t_0^2 + k_139 * t_0 + k_58)
  //     - (k_42 * t_0^2 + k_43 * t_0 + k_44) * (FH_2 * t_0^2 + FH_3 * t_0 + FH_4) * (FH_0 * t_0 + FH_1)
  //     - 2 * t_0 * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2) + 2 * (k_13 * t_0^2 - k_19 * t_0) * (V_00 * t_0 + k_1) * (FH_0^2 * t_0^2 + k_24 * t_0 + FH_1^2)
  // 0 = +- (
  //         (k_42 * FH_0 - k_40 * FH_2) * t_0^3
  //         + (k_42 * FH_1 + k_43 * FH_0 - k_41 * FH_2 - k_40 * FH_3) * t_0^2
  //         + (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^2 + 4 * f_0 * f_2)
  //     + (k_40 * k_31 - k_42 * FH_2 * FH_0) * t_0^5
  //     + (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^4
  //     + (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^3
  //     + (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^2
  //     + (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^2 - k_16 * FH_0^2) * t_0^5
  //     + 2 * (k_13 * V_00 * k_24 + k_13 * k_1 * FH_0^2 - k_19 * V_00 * FH_0^2 - k_16 * k_24 - k_17 * FH_0^2) * t_0^4
  //     + 2 * (k_13 * V_00 * FH_1^2 + k_13 * k_1 * k_24 - k_19 * V_00 * k_24 - k_19 * k_1 * FH_0^2 - k_16 * FH_1^2 - k_17 * k_24 - k_18 * FH_0^2) * t_0^3
  //     + 2 * (k_13 * k_1 * FH_1^2 - k_19 * V_00 * FH_1^2 - k_19 * k_1 * k_24 - k_17 * FH_1^2 - k_18 * k_24) * t_0^2
  //     + 2 * (- k_19 * k_1 * FH_1^2 - k_18 * FH_1^2) * t_0
  // 0 = +- (k_45 * t_0^3 + k_46 * t_0^2 + k_47 * t_0 + k_48) * sqrt(f_1^2 + 4 * f_0 * f_2)
  //     + (k_50 + k_62) * t_0^5 + (k_52 + k_64) * t_0^4 + (k_54 + k_66) * t_0^3 + (k_56 + k_67) * t_0^2 + (k_59 + k_68) * t_0 + k_60
  // 0 = +- (k_45 * t_0^3 + k_46 * t_0^2 + k_47 * t_0 + k_48) * sqrt(k_31 * t_0^4 + k_132 * t_0^3 + k_133 * t_0^2 + k_134 * t_0 + k_58)
  //     + k_69 * t_0^5 + k_70 * t_0^4 + k_71 * t_0^3 + k_72 * t_0^2 + k_73 * t_0 + k_60

  OPolynomial0 := TBigIntPolynomial.Create([k[60], k[73], k[72], k[71], k[70], k[69]]);
  OPolynomial1 := TBigIntPolynomial.Create([k[48], k[47], k[46], k[45]]);

  // Squaring that formula eliminates the square root, but may lead to a polynomial with all coefficients zero in some
  // cases. Therefore this part is merely included for the interested reader.
  // -+ (k_45 * t_0^3 + k_46 * t_0^2 + k_47 * t_0 + k_48) * sqrt(k_31 * t_0^4 + k_132 * t_0^3 + k_133 * t_0^2 + k_134 * t_0 + k_58) =
  //     k_69 * t_0^5 + k_70 * t_0^4 + k_71 * t_0^3 + k_72 * t_0^2 + k_73 * t_0 + k_60
  // (k_45 * t_0^3 + k_46 * t_0^2 + k_47 * t_0 + k_48)^2 * (k_31 * t_0^4 + k_132 * t_0^3 + k_133 * t_0^2 + k_134 * t_0 + k_58) =
  //     (k_69 * t_0^5 + k_70 * t_0^4 + k_71 * t_0^3 + k_72 * t_0^2 + k_73 * t_0 + k_60)^2
  // 0 =
  //     (k_45^2 * t_0^6
  //         + 2 * k_45 * k_46 * t_0^5
  //         + k_46^2 * t_0^4 + 2 * k_45 * k_47 * t_0^4
  //         + 2 * k_45 * k_48 * t_0^3 + 2 * k_46 * k_47 * t_0^3
  //         + k_47^2 * t_0^2 + 2 * k_46 * k_48 * t_0^2
  //         + 2 * k_47 * k_48 * t_0
  //         + k_48^2
  //     ) * (k_31 * t_0^4 + k_132 * t_0^3 + k_133 * t_0^2 + k_134 * t_0 + k_58)
  //     - k_69^2 * t_0^10
  //     - 2 * k_69 * k_70 * t_0^9
  //     - (k_70^2 + 2 * k_69 * k_71) * t_0^8
  //     - 2 * (k_69 * k_72 + k_70 * k_71) * t_0^7
  //     - (k_71^2 + 2 * k_69 * k_73 + 2 * k_70 * k_72) * t_0^6
  //     - 2 * (k_69 * k_60 + k_70 * k_73 + k_71 * k_72) * t_0^5
  //     - (k_72^2 + 2 * k_70 * k_60 + 2 * k_71 * k_73) * t_0^4
  //     - 2 * (k_71 * k_60 + k_72 * k_73) * t_0^3
  //     - (k_73^2 + 2 * k_72 * k_60) * t_0^2
  //     - 2 * k_73 * k_60 * t_0
  //     - k_60^2
  // 0 = ak_10 * t_0^10 + ak_9 * t_0^9 + ak_8 * t_0^8 + ak_7 * t_0^7 + ak_6 * t_0^6 + ak_5 * t_0^5 + ak_4 * t_0^4 + ak_3 * t_0^3 + ak_2 * t_0^2 + ak_1 * t_0 + ak_0

  //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];
  //ak[0] := k[58] * k[80] - k[60] * k[60];
  //ak[1] := k[134] * k[80] + k[58] * k[79] - 2 * k[73] * k[60];
  //ak[2] := k[133] * k[80] + k[134] * k[79] + k[58] * k[78] - k[73] * k[73] - 2 * k[72] * k[60];
  //ak[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]);
  //ak[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]);
  //ak[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]);
  //ak[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]);
  //ak[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]);
  //ak[8] := k[31] * k[76] + k[132] * k[75] + k[133] * k[74] - k[70] * k[70] - 2 * k[69] * k[71];
  //ak[9] := k[31] * k[75] + k[132] * k[74] - 2 * k[69] * k[70];
  //ak[10] := k[31] * k[74] - k[69] * k[69];
end;

function TNeverTellMeTheOdds.CalcRockThrowCollisionOptions(constref AHailstone0, AHailstone1, AHailstone2: THailstone):
  TInt64Array;
var
  a0, a1: TBigIntPolynomial;
  a0Roots, a1Roots: TBigIntArray;
  options: specialize TList<Int64>;
  i, j: TBigInt;
  val: Int64;
begin
  CalcCollisionPolynomials(AHailstone0, AHailstone1, AHailstone2, a0, a1);
  a0Roots := TPolynomialRoots.BisectInteger(a0, 64);
  a1Roots := TPolynomialRoots.BisectInteger(a1, 64);

  options := specialize TList<Int64>.Create;
  for i in a0Roots do
    for j in a1Roots do
      if (i = j) and i.TryToInt64(val) then
        options.Add(val);
  Result := options.ToArray;
  options.Free;
end;

function TNeverTellMeTheOdds.ValidateRockThrow(constref AHailstone0, AHailstone1, AHailstone2: THailstone; const AT0,
  AT1: Int64): Int64;
var
  divisor, t: Int64;
  rock: THailstone;
begin
  // V_x = (V_0 * t_0 - V_1 * t_1 + P_0 - P_1) / (t_0 - t_1)
  divisor := AT0 - AT1;
  rock := THailstone.Create;
  rock.V0 := (AHailstone0.V0 * AT0 - AHailstone1.V0 * AT1 + AHailstone0.P0 - AHailstone1.P0) div divisor;
  rock.V1 := (AHailstone0.V1 * AT0 - AHailstone1.V1 * AT1 + AHailstone0.P1 - AHailstone1.P1) div divisor;
  rock.V2 := (AHailstone0.V2 * AT0 - AHailstone1.V2 * AT1 + AHailstone0.P2 - AHailstone1.P2) div divisor;

  // P_x = (V_0 - V_x) * t_0 + P_0
  rock.P0 := (AHailstone0.V0 - rock.V0) * AT0 + AHailstone0.P0;
  rock.P1 := (AHailstone0.V1 - rock.V1) * AT0 + AHailstone0.P1;
  rock.P2 := (AHailstone0.V2 - rock.V2) * AT0 + AHailstone0.P2;

  Result := rock.P0 + rock.P1 + rock.P2;

  // Checks collision with the third hailstone.
  if ((AHailstone2.V0 = rock.V0) and (AHailstone2.P0 <> rock.P0))
  or ((AHailstone2.V1 = rock.V1) and (AHailstone2.P1 <> rock.P1))
  or ((AHailstone2.V2 = rock.V2) and (AHailstone2.P2 <> rock.P2)) then
    Result := 0
  else begin
    t := (AHailstone2.P0 - rock.P0) div (rock.V0 - AHailstone2.V0);
    if (t <> (AHailstone2.P1 - rock.P1) div (rock.V1 - AHailstone2.V1))
    or (t <> (AHailstone2.P2 - rock.P2) div (rock.V2 - AHailstone2.V2)) then
      Result := 0;
  end;

  rock.Free;
end;

constructor TNeverTellMeTheOdds.Create(const AMin: Int64; const AMax: Int64);
begin
  FMin := AMin;
  FMax := AMax;
  FHailstones := THailstones.Create;
end;

destructor TNeverTellMeTheOdds.Destroy;
begin
  FHailstones.Free;
  inherited Destroy;
end;

procedure TNeverTellMeTheOdds.ProcessDataLine(const ALine: string);
begin
  FHailstones.Add(THailstone.Create(ALine));
end;

procedure TNeverTellMeTheOdds.Finish;
var
  i, j: 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
        Inc(FPart1);

  if FHailstones.Count >= 3 then
    FPart2 := FindRockThrow(0, 1, 2);
end;

function TNeverTellMeTheOdds.GetDataFileName: string;
begin
  Result := 'never_tell_me_the_odds.txt';
end;

function TNeverTellMeTheOdds.GetPuzzleName: string;
begin
  Result := 'Day 24: Never Tell Me The Odds';
end;

end.