From 3e3e1d45d336a2bcd2312b2595b120bb6bd4a25f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20M=C3=BCller?= Date: Mon, 27 May 2024 02:29:49 +0200 Subject: [PATCH] Added solution "Day 24: Never Tell Me The Odds", part 2 --- solvers/UNeverTellMeTheOdds.pas | 773 +++++++++++-------------- tests/UNeverTellMeTheOddsTestCases.pas | 2 +- 2 files changed, 332 insertions(+), 443 deletions(-) diff --git a/solvers/UNeverTellMeTheOdds.pas b/solvers/UNeverTellMeTheOdds.pas index b4b2e82..793a7aa 100644 --- a/solvers/UNeverTellMeTheOdds.pas +++ b/solvers/UNeverTellMeTheOdds.pas @@ -22,7 +22,7 @@ unit UNeverTellMeTheOdds; interface uses - Classes, SysUtils, Generics.Collections, Math, matrix, USolver, UNumberTheory, UBigInt; + Classes, SysUtils, Generics.Collections, Math, USolver, UNumberTheory, UBigInt, UPolynomial, UPolynomialRoots; type @@ -30,26 +30,15 @@ type THailstone = class public - Position, Velocity: Tvector3_extended; + P0, P1, P2: Int64; + V0, V1, V2: Integer; constructor Create(const ALine: string); - constructor Create(const APosition, AVelocity: Tvector3_extended); + constructor Create; end; THailstones = specialize TObjectList; - { 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; + TInt64Array = array of Int64; { TNeverTellMeTheOdds } @@ -57,10 +46,15 @@ type private FMin, FMax: Int64; FHailstones: THailstones; - FA: array[0..10] of Int64; - FH: array[0..6] of Int64; + FA: array[0..10] of TBigInt; + FH: array[0..6] of TBigInt; function AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean; - procedure FindRockThrow(const AIndex1, AIndex2, AIndex3: Integer); + 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; @@ -70,10 +64,6 @@ type function GetPuzzleName: string; override; end; -const - CIterationThreshold = 0.00001; - CEpsilon = 0.0000000001; - implementation { THailstone } @@ -83,69 +73,71 @@ 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]))); + 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(const APosition, AVelocity: Tvector3_extended); +constructor THailstone.Create; begin - Position := APosition; - Velocity := AVelocity; + end; -{ TFirstCollisionPolynomial } +{ TNeverTellMeTheOdds } -procedure TFirstCollisionPolynomial.NormalizeCoefficients; +function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean; var - shift: Integer; - i: Low(FA)..High(FA); - //gcd: TBigInt; + m1, m2, x, y: Double; begin - // Eliminates zero constant term. - shift := 0; - while (shift <= High(FA)) and (FA[shift] = 0) do - Inc(shift); - - if shift <= High(FA) then + Result := False; + m1 := AHailstone1.V1 / AHailstone1.V0; + m2 := AHailstone2.V1 / AHailstone2.V0; + if m1 <> m2 then begin - if shift > 0 then + 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 - 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; + y := m1 * (x - AHailstone1.P0) + AHailstone1.P1; + if (FMin <= y) and (y <= FMax) then + Result := True 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); +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 - 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 @@ -155,55 +147,88 @@ begin // 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) + // 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) + // 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) + // 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: 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 + // 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. - 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]); + // 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 - 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; + // 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]; @@ -213,10 +238,45 @@ begin 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]; + // 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]; @@ -224,34 +284,19 @@ begin 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[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[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[28] := FH[1] * FH[1]; 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[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]; @@ -269,9 +314,9 @@ begin 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[61] := k[13] * AHailstone0.V0 - k[16]; k[62] := 2 * k[25] * k[61]; - k[63] := k[13] * k[1] - k[19] * V_00 - k[17]; + 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]); @@ -282,319 +327,167 @@ begin 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]; + // 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]; - 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]; + // 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 - // 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], ')'); + 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]]); - 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], - ' '); + // 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 - WriteLn('debug96: ', EvaluateAt(t_0) = 0); - - NormalizeCoefficients; - - WriteLn('debug99: ', EvaluateAt(t_0) = 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 TFirstCollisionPolynomial.EvaluateAt(const AT0: Int64): TBigInt; +function TNeverTellMeTheOdds.CalcRockThrowCollisionOptions(constref AHailstone0, AHailstone1, AHailstone2: THailstone): + TInt64Array; var - i: Low(FA)..High(FA); + a0, a1: TBigIntPolynomial; + a0Roots, a1Roots: TBigIntArray; + options: specialize TList; + i, j: TBigInt; + val: Int64; begin - Result := TBigInt.Zero; - for i := High(FA) downto Low(FA) do - Result := Result * AT0 + FA[i]; + CalcCollisionPolynomials(AHailstone0, AHailstone1, AHailstone2, a0, a1); + a0Roots := TPolynomialRoots.BisectInteger(a0, 64); + a1Roots := TPolynomialRoots.BisectInteger(a1, 64); + + options := specialize TList.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 TFirstCollisionPolynomial.CalcPositiveIntegerRoot: Int64; +function TNeverTellMeTheOdds.ValidateRockThrow(constref AHailstone0, AHailstone1, AHailstone2: THailstone; const AT0, + AT1: Int64): Int64; var - dividers: TDividers; - factors: TInt64Array; - divider: Int64; + divisor, t: Int64; + rock: THailstone; 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; + // 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; -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; + // 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; -{ TNeverTellMeTheOdds } + Result := rock.P0 + rock.P1 + rock.P2; -function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean; -var - m1, m2, x, y: Double; -begin - Result := False; - 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.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.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.Position.data[0]) + AHailstone1.Position.data[1]; - if (FMin <= y) and (y <= FMax) then - Result := True - end; + // 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; -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; + rock.Free; end; constructor TNeverTellMeTheOdds.Create(const AMin: Int64; const AMax: Int64); @@ -617,19 +510,15 @@ end; procedure TNeverTellMeTheOdds.Finish; var - i, j, k: Integer; + 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); - 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); + if FHailstones.Count >= 3 then + FPart2 := FindRockThrow(0, 1, 2); end; function TNeverTellMeTheOdds.GetDataFileName: string; diff --git a/tests/UNeverTellMeTheOddsTestCases.pas b/tests/UNeverTellMeTheOddsTestCases.pas index a2f1235..2f56a79 100644 --- a/tests/UNeverTellMeTheOddsTestCases.pas +++ b/tests/UNeverTellMeTheOddsTestCases.pas @@ -81,7 +81,7 @@ end; procedure TNeverTellMeTheOddsFullDataTestCase.TestPart2; begin - AssertEquals(-1, FSolver.GetResultPart2); + AssertEquals(856642398547748, FSolver.GetResultPart2); end; { TNeverTellMeTheOddsExampleTestCase }