warrence
/
AdventOfCode2023
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Added WIP analytical solution attempt

This commit is contained in:
Stefan Müller 2024-02-14 12:07:12 +01:00
parent 356cc2ad5e
commit 44c2c845e0
2 changed files with 560 additions and 30 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

578
solvers/UNeverTellMeTheOdds.pas
View File

@ -1,6 +1,6 @@
{ {
Solutions to the Advent Of Code. 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 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 the terms of the GNU General Public License as published by the Free Software
@ -22,26 +22,45 @@ unit UNeverTellMeTheOdds;
interface interface
uses uses
Classes, SysUtils, Generics.Collections, Math, USolver; Classes, SysUtils, Generics.Collections, Math, matrix, USolver, UNumberTheory, UBigInt;
type type
{ THailstone } { THailstone }
THailstone = record THailstone = class
X, Y, Z: Int64; public
VX, VY, VZ: Integer; Position, Velocity: Tvector3_extended;
constructor Create(const ALine: string);
constructor Create(const APosition, AVelocity: Tvector3_extended);
end; 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 }
TNeverTellMeTheOdds = class(TSolver) TNeverTellMeTheOdds = class(TSolver)
private private
FMin, FMax: Int64; 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; function AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
procedure FindRockThrow(const AIndex1, AIndex2, AIndex3: Integer);
public public
constructor Create(const AMin: Int64 = 200000000000000; const AMax: Int64 = 400000000000000); constructor Create(const AMin: Int64 = 200000000000000; const AMax: Int64 = 400000000000000);
destructor Destroy; override; destructor Destroy; override;
@ -51,8 +70,452 @@ type
function GetPuzzleName: string; override; function GetPuzzleName: string; override;
end; end;
const
CIterationThreshold = 0.00001;
CEpsilon = 0.0000000001;
implementation 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 } { TNeverTellMeTheOdds }
function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean; function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
@ -60,58 +523,113 @@ var
m1, m2, x, y: Double; m1, m2, x, y: Double;
begin begin
Result := False; Result := False;
m1 := AHailstone1.VY / AHailstone1.VX; m1 := AHailstone1.Velocity.data[1] / AHailstone1.Velocity.data[0];
m2 := AHailstone2.VY / AHailstone2.VX; m2 := AHailstone2.Velocity.data[1] / AHailstone2.Velocity.data[0];
if m1 <> m2 then if m1 <> m2 then
begin 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) if (FMin <= x) and (x <= FMax)
and (x * Sign(AHailstone1.VX) >= AHailstone1.X * Sign(AHailstone1.VX)) and (x * Sign(AHailstone1.Velocity.data[0]) >= AHailstone1.Position.data[0] * Sign(AHailstone1.Velocity.data[0]))
and (x * Sign(AHailstone2.VX) >= AHailstone2.X * Sign(AHailstone2.VX)) then and (x * Sign(AHailstone2.Velocity.data[0]) >= AHailstone2.Position.data[0] * Sign(AHailstone2.Velocity.data[0]))
then
begin 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 if (FMin <= y) and (y <= FMax) then
Result := True Result := True
end; end;
end; 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); constructor TNeverTellMeTheOdds.Create(const AMin: Int64; const AMax: Int64);
begin begin
FMin := AMin; FMin := AMin;
FMax := AMax; FMax := AMax;
FHailStones := THailstones.Create; FHailstones := THailstones.Create;
end; end;
destructor TNeverTellMeTheOdds.Destroy; destructor TNeverTellMeTheOdds.Destroy;
begin begin
FHailStones.Free; FHailstones.Free;
inherited Destroy; inherited Destroy;
end; end;
procedure TNeverTellMeTheOdds.ProcessDataLine(const ALine: string); procedure TNeverTellMeTheOdds.ProcessDataLine(const ALine: string);
var
split: TStringArray;
hailstone: THailstone;
begin begin
split := ALine.Split([',', '@']); FHailstones.Add(THailstone.Create(ALine));
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);
end; end;
procedure TNeverTellMeTheOdds.Finish; procedure TNeverTellMeTheOdds.Finish;
var var
i, j: Integer; i, j, k: Integer;
begin begin
for i := 0 to FHailStones.Count - 2 do for i := 0 to FHailstones.Count - 2 do
for j := i + 1 to FHailStones.Count - 1 do for j := i + 1 to FHailstones.Count - 1 do
if AreIntersecting(FHailStones[i], FHailStones[j]) then if AreIntersecting(FHailstones[i], FHailstones[j]) then
Inc(FPart1); 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; end;
function TNeverTellMeTheOdds.GetDataFileName: string; function TNeverTellMeTheOdds.GetDataFileName: string;

12
tests/UNeverTellMeTheOddsTestCases.pas
View File

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