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

Merge branch 'day24-analytical'

This commit is contained in:
Stefan Müller 2024-05-27 02:53:03 +02:00
parent eca6b8f3f9 b27b14a153
commit 859a5db921
12 changed files with 3313 additions and 35 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
AdventOfCode.lpi
View File

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

763
UBigInt.pas Normal file
View File

@ -0,0 +1,763 @@
{
Solutions to the Advent Of Code.
Copyright (C) 2022-2024 Stefan Müller
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation, either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
}
unit UBigInt;
{$mode ObjFPC}{$H+}
interface
uses
Classes, SysUtils, Math;
type
TDigits = array of Cardinal;
{ TBigInt }
// This is an abbreviated reimplementation in Freepascal of a C# class created in 2022.
TBigInt = object
private
FDigits: TDigits;
FIsNegative: Boolean;
function GetSign: Integer;
// Copies consecutive digits from this BigInt to create a new one. The result will be positive. Leading zeros are
// removed from the result, but AIndex + ACount must not exceed the number of digits of this BigInt.
// AIndex is the first (least significant) digit to be taken. The digit with this index will become the 0th digit of
// the new BigInt.
// ACount is the number of consecutive digits to be taken, and the number of digits of the result.
function GetSegment(const AIndex, ACount: Integer): TBigInt;
// Compares the absolute value of this TBigInt object to the absolute value of another one. Returns -1 if this
// object is less than AOther, 1 if this object is greater than AOther, and 0 if they are equal.
function CompareToAbsoluteValues(constref AOther: TBigInt): Integer;
class function GetZero: TBigInt; static;
class function GetOne: TBigInt; static;
// Adds A and B, ignoring their signs and using ReturnNegative instead. The result is
// Sign * (Abs(A) + Abs(B)),
// where Sign is 1 for ReturnNegative = False and -1 otherwise.
class function AddAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt; static;
// Subtracts B from A, ignoring their signs. However, the result might be negative, and the sign can be reversed by
// setting ReturnNegative to True. The result is
// Sign * (Abs(A) - Abs(B)),
// where Sign is 1 for ReturnNegative = False and -1 otherwise.
class function SubtractAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt; static;
// Multiplies A and B, ignoring their signs and using ReturnNegative instead. This multiplication uses a recursive
// implementation of the Karatsuba algorithm. See
// https://www.geeksforgeeks.org/karatsuba-algorithm-for-fast-multiplication-using-divide-and-conquer-algorithm/
// The result is
// Sign * (Abs(a) * Abs(b))
// where Sign is 1 for ReturnNegative = False and -1 otherwise.
class function MultiplyAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt; static;
class function FromHexOrBinString(const AValue: string; const AFromBase: Integer): TBigInt; static;
class function ConvertDigitBlock(const AValue: string; var AStartIndex: Integer; const ACharBlockSize, AFromBase:
Integer): Cardinal;
public
property IsNegative: Boolean read FIsNegative;
property Sign: Integer read GetSign;
class property Zero: TBigInt read GetZero;
class property One: TBigInt read GetOne;
// Returns the index of the most significant bit, i.e. returns integer k, where 2^k is the largest power of 2 that
// is less than or equal to the absolute value of the number itself. Returns -1 if the given number is 0.
function GetMostSignificantBitIndex: Int64;
function CompareTo(constref AOther: TBigInt): Integer;
function TryToInt64(out AOutput: Int64): Boolean;
// TODO: ToString is currently for debug output only.
function ToString: string;
class function FromInt64(const AValue: Int64): TBigInt; static;
class function FromHexadecimalString(const AValue: string): TBigInt; static;
class function FromBinaryString(const AValue: string): TBigInt; static;
end;
TBigIntArray = array of TBigInt;
{ Operators }
operator := (const A: Int64): TBigInt;
operator - (const A: TBigInt): TBigInt;
operator + (const A, B: TBigInt): TBigInt;
operator - (const A, B: TBigInt): TBigInt;
operator * (const A, B: TBigInt): TBigInt;
operator shl (const A: TBigInt; const B: Integer): TBigInt;
operator shr (const A: TBigInt; const B: Integer): TBigInt;
operator = (const A, B: TBigInt): Boolean;
operator <> (const A, B: TBigInt): Boolean;
operator < (const A, B: TBigInt): Boolean;
operator <= (const A, B: TBigInt): Boolean;
operator > (const A, B: TBigInt): Boolean;
operator >= (const A, B: TBigInt): Boolean;
implementation
const
CBase = Cardinal.MaxValue + 1;
CMaxDigit = Cardinal.MaxValue;
CDigitSize = SizeOf(Cardinal);
CBitsPerDigit = CDigitSize * 8;
CHalfBits = CBitsPerDigit >> 1;
CHalfDigitMax = (1 << CHalfBits) - 1;
CZero: TBigInt = (FDigits: (0); FIsNegative: False);
COne: TBigInt = (FDigits: (1); FIsNegative: False);
{ TBigInt }
function TBigInt.GetSign: Integer;
begin
if FIsNegative then
Result := -1
else if (Length(FDigits) > 1) or (FDigits[0] <> 0) then
Result := 1
else
Result := 0;
end;
function TBigInt.GetSegment(const AIndex, ACount: Integer): TBigInt;
var
trimmedCount: Integer;
begin
trimmedCount := ACount;
while (trimmedCount > 1) and (FDigits[AIndex + trimmedCount - 1] = 0) do
Dec(trimmedCount);
Result.FDigits := Copy(FDigits, AIndex, trimmedCount);
Result.FIsNegative := False;
end;
function TBigInt.CompareToAbsoluteValues(constref AOther: TBigInt): Integer;
var
i: Integer;
begin
Result := Length(FDigits) - Length(AOther.FDigits);
if Result = 0 then
begin
for i := High(FDigits) downto 0 do
if FDigits[i] < AOther.FDigits[i] then
begin
Result := -1;
Break;
end
else if FDigits[i] > AOther.FDigits[i] then
begin
Result := 1;
Break;
end;
end;
end;
class function TBigInt.GetZero: TBigInt;
begin
Result := CZero;
end;
class function TBigInt.GetOne: TBigInt;
begin
Result := COne;
end;
class function TBigInt.AddAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
var
i, j, lenA, lenB, len, shorter: Integer;
carry: Cardinal;
begin
lenA := Length(AA.FDigits);
lenB := Length(AB.FDigits);
// Initializes the digits array of the result, with a simple test to try to predict a carry-over into a new digit. The
// result could still carry into new digit depending on lower digits (carry over multiple digits), which would be
// handled at the end.
if lenA = lenB then
if CMaxDigit - AA.FDigits[lenA - 1] < AB.FDigits[lenB - 1] then
// Result will carry into new digit.
SetLength(Result.FDigits, lenA + 1)
else
SetLength(Result.FDigits, lenA)
else
SetLength(Result.FDigits, Max(lenA, lenB));
len := Length(Result.FDigits);
// Calculates the new digits from less to more significant until the end of the shorter operand is reached.
shorter := Min(Length(AA.FDigits), Length(AB.FDigits));
i := 0;
carry := 0;
while i < shorter do
begin
if (AB.FDigits[i] = CMaxDigit) and (carry > 0) then
begin
Result.FDigits[i] := AA.FDigits[i];
carry := 1;
end
else
if CMaxDigit - AA.FDigits[i] < AB.FDigits[i] + carry then
begin
Result.FDigits[i] := AB.FDigits[i] + carry - 1 - (CMaxDigit - AA.FDigits[i]);
carry := 1;
end
else begin
Result.FDigits[i] := AA.FDigits[i] + AB.FDigits[i] + carry;
carry := 0;
end;
Inc(i);
end;
// Copies the missing unchanged digits from the longer operand to the result, if any, before applying remaining
// carry-overs. This avoids additional tests for finding the shorter digit array.
if (i < lenA) or (i < lenB) then
if lenA >= lenB then
for j := i to len - 1 do
Result.FDigits[j] := AA.FDigits[j]
else
for j := i to len - 1 do
Result.FDigits[j] := AB.FDigits[j];
// Applies the remaining carry-overs until the end of the prepared result digit array.
while (carry > 0) and (i < len) do
begin
if Result.FDigits[i] = CMaxDigit then
Result.FDigits[i] := 0
else begin
Inc(Result.FDigits[i]);
carry := 0;
end;
Inc(i);
end;
// Applies the carry-over into a new digit that was not anticipated in the initialization at the top (carry over
// multiple digits).
if carry > 0 then
Insert(1, Result.FDigits, len);
Result.FIsNegative := AReturnNegative;
end;
class function TBigInt.SubtractAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
var
a, b: TBigInt;
carry: Cardinal;
compare, i, j, lastNonZeroDigitIndex, len: Integer;
begin
// Establishes the operand order, such that Abs(a) is not less than Abs(b).
compare := AA.CompareToAbsoluteValues(AB);
if compare = 0 then
begin
Result := Zero;
Exit;
end;
if compare > 0 then
begin
a := AA;
b := AB;
Result.FIsNegative := AReturnNegative;
end
else begin
a := AB;
b := AA;
Result.FIsNegative := not AReturnNegative;
end;
// Calculates the new digits from less to more significant until the end of the shorter operand is reached and all
// carry-overs have been applied.
len := Length(a.FDigits);
SetLength(Result.FDigits, len);
carry := 0;
// Tracks leading zeros for the trim below.
lastNonZeroDigitIndex := 0;
i := 0;
while i < Length(b.FDigits) do
begin
if (a.FDigits[i] = b.FDigits[i]) and (carry > 0) then
begin
Result.FDigits[i] := CMaxDigit;
carry := 1;
lastNonZeroDigitIndex := i;
end
else begin
if a.FDigits[i] < b.FDigits[i] then
begin
Result.FDigits[i] := CMaxDigit - (b.FDigits[i] - a.FDigits[i]) + 1 - carry;
carry := 1;
end
else begin
Result.FDigits[i] := a.FDigits[i] - b.FDigits[i] - carry;
carry := 0;
end;
if Result.FDigits[i] > 0 then
lastNonZeroDigitIndex := i;
end;
Inc(i);
end;
while carry > 0 do
begin
if a.FDigits[i] = 0 then
begin
Result.FDigits[i] := CMaxDigit;
lastNonZeroDigitIndex := i;
end
else begin
Result.FDigits[i] := a.FDigits[i] - carry;
carry := 0;
if Result.FDigits[i] > 0 then
lastNonZeroDigitIndex := i;
end;
Inc(i);
end;
// Copies the missing unchanged digits from the longer operand to the result, if any. If there are none, then no trim
// needs to occur because the most significant digit is not zero.
if i < len then
for j := i to len - 1 do
Result.FDigits[j] := a.FDigits[j]
else if (lastNonZeroDigitIndex + 1 < len) then
// Trims leading zeros from the digits array.
Delete(Result.FDigits, lastNonZeroDigitIndex + 1, len - lastNonZeroDigitIndex - 1);
end;
class function TBigInt.MultiplyAbsoluteValues(constref AA, AB: TBigInt; const AReturnNegative: Boolean): TBigInt;
var
lenA, lenB, lenMax, floorHalfLength, ceilHalfLength: Integer;
a1, a0, b1, b0, a1b1, a0b0: Cardinal;
am, bm, middle, biga1, biga0, bigb1, bigb0, biga1b1, biga0b0: TBigInt;
begin
lenA := Length(AA.FDigits);
lenB := Length(AB.FDigits);
if (lenA <= 1) and (lenB <= 1) then
begin
if (AA.FDigits[0] <= CHalfDigitMax) and (AB.FDigits[0] <= CHalfDigitMax) then
if (AA.FDigits[0] = 0) or (AB.FDigits[0] = 0) then
Result := Zero
else begin
Result.FDigits := TDigits.Create(AA.FDigits[0] * AB.FDigits[0]);
Result.FIsNegative := AReturnNegative;
end
else begin
// a1, a0, b1, b0 use only the lower (less significant) half of the bits of a digit, so the product of any two of
// these fits in one digit.
a1 := AA.FDigits[0] >> CHalfBits;
a0 := AA.FDigits[0] and CHalfDigitMax;
b1 := AB.FDigits[0] >> CHalfBits;
b0 := AB.FDigits[0] and CHalfDigitMax;
a1b1 := a1 * b1;
a0b0 := a0 * b0;
if a1b1 > 0 then
Result.FDigits := TDigits.Create(a0b0, a1b1)
else
Result.FDigits := TDigits.Create(a0b0);
Result.FIsNegative := AReturnNegative;
// The result of (a1 + a0) * (b1 + b0) might not fit in one digit, so one last recursion step is necessary.
am := FromInt64(a1 + a0);
bm := FromInt64(b1 + b0);
middle := (MultiplyAbsoluteValues(am, bm, False) - a1b1 - a0b0) << CHalfBits;
if AReturnNegative then
Result := Result - middle
else
Result := Result + middle;
end;
end
else begin
// Calculates where to split the two numbers.
lenMax := Max(lenA, lenB);
floorHalfLength := lenMax >> 1;
ceilHalfLength := lenMax - floorHalfLength;
// Performs one recursion step.
if ceilHalfLength < lenA then
begin
biga1 := AA.GetSegment(ceilHalfLength, lenA - ceilHalfLength);
biga0 := AA.GetSegment(0, ceilHalfLength);
if ceilHalfLength < lenB then
begin
bigb1 := AB.GetSegment(ceilHalfLength, lenB - ceilHalfLength);
bigb0 := AB.GetSegment(0, ceilHalfLength);
biga1b1 := MultiplyAbsoluteValues(biga1, bigb1, AReturnNegative);
biga0b0 := MultiplyAbsoluteValues(biga0, bigb0, AReturnNegative);
Result := (biga1b1 << (2 * ceilHalfLength * CBitsPerDigit))
+ ((MultiplyAbsoluteValues(biga1 + biga0, bigb1 + bigb0, AReturnNegative) - biga1b1 - biga0b0) << (ceilHalfLength * CBitsPerDigit))
+ biga0b0;
end
else begin
biga0b0 := MultiplyAbsoluteValues(biga0, AB, AReturnNegative);
Result := ((MultiplyAbsoluteValues(biga1 + biga0, AB, AReturnNegative) - biga0b0) << (ceilHalfLength * CBitsPerDigit)) + biga0b0;
end;
end
else begin
bigb1 := AB.GetSegment(ceilHalfLength, lenB - ceilHalfLength);
bigb0 := AB.GetSegment(0, ceilHalfLength);
biga0b0 := MultiplyAbsoluteValues(AA, bigb0, AReturnNegative);
Result := ((MultiplyAbsoluteValues(AA, bigb1 + bigb0, AReturnNegative) - biga0b0) << (ceilHalfLength * CBitsPerDigit)) + biga0b0;
end;
end;
end;
class function TBigInt.FromHexOrBinString(const AValue: string; const AFromBase: Integer): TBigInt;
var
charBlockSize, offset, i, j, k, remainder: Integer;
d: Cardinal;
begin
// 2 ^ (32 / charBlockSize) = AFromBase
case AFromBase of
2: charBlockSize := 32;
16: charBlockSize := 8;
end;
if AValue[1] = '-' then
begin
offset := 2;
Result.FIsNegative := True;
end
else begin
offset := 1;
Result.FIsNegative := False;
end;
// Calculates the first (most significant) digit d of the result.
DivMod(AValue.Length - offset + 1, charBlockSize, i, remainder);
k := offset;
d := 0;
// Checks the first block of chars that is not a full block.
if remainder > 0 then
d := ConvertDigitBlock(AValue, k, remainder, AFromBase);
// Checks full blocks of chars for first digit.
while (d = 0) and (i > 0) do
begin
Dec(i);
d := ConvertDigitBlock(AValue, k, charBlockSize, AFromBase);
end;
// Checks for zero.
if (d = 0) and (i = 0) then
Result := Zero
else begin
// Initializes the array of digits.
SetLength(Result.FDigits, i + 1);
Result.FDigits[i] := d;
// Calculates the other digits.
for j := i - 1 downto 0 do
Result.FDigits[j] := ConvertDigitBlock(AValue, k, charBlockSize, AFromBase);
end;
end;
class function TBigInt.ConvertDigitBlock(const AValue: string; var AStartIndex: Integer; const ACharBlockSize,
AFromBase: Integer): Cardinal;
var
part: string;
begin
part := Copy(AValue, AStartIndex, ACharBlockSize);
Inc(AStartIndex, ACharBlockSize);
case AFromBase of
2: part := '%' + part;
8: part := '&' + part;
16: part := '$' + part;
end;
Result := StrToDWord(part);
end;
function TBigInt.GetMostSignificantBitIndex: Int64;
var
high, i: Integer;
begin
high := Length(FDigits) - 1;
if (high = 0) and (FDigits[0] = 0) then
Result := -1
else begin
i := CBitsPerDigit - 1;
while ((1 << i) and FDigits[high]) = 0 do
Dec(i);
Result := high * CBitsPerDigit + i;
end;
end;
function TBigInt.CompareTo(constref AOther: TBigInt): Integer;
begin
if FIsNegative = AOther.FIsNegative then
Result := CompareToAbsoluteValues(AOther)
else
Result := 1;
if FIsNegative then
Result := -Result;
end;
function TBigInt.TryToInt64(out AOutput: Int64): Boolean;
begin
AOutput := 0;
Result := False;
case Length(FDigits) of
0: Result := True;
1: begin
AOutput := FDigits[0];
if FIsNegative then
AOutput := -AOutput;
Result := True;
end;
2: begin
if FDigits[1] <= Integer.MaxValue then
begin
AOutput := FDigits[1] * CBase + FDigits[0];
if FIsNegative then
AOutput := -AOutput;
Result := True;
end
else if (FDigits[1] = Integer.MaxValue + 1) and (FDigits[0] = 0) and FIsNegative then
begin
AOutput := Int64.MinValue;
Result := True;
end;
end;
end;
end;
function TBigInt.ToString: string;
var
i: Integer;
begin
if FIsNegative then
Result := '(-'
else
Result := '';
for i := 0 to Length(FDigits) - 2 do
Result := Result + '(' + IntToStr(FDigits[i]) + ' + 2^32 * ';
Result := Result + IntToStr(FDigits[Length(FDigits) - 1]);
for i := 0 to Length(FDigits) - 2 do
Result := Result + ')';
if FIsNegative then
Result := Result + ')'
end;
class function TBigInt.FromInt64(const AValue: Int64): TBigInt;
var
absVal: Int64;
begin
if AValue <> Int64.MinValue then
begin
absVal := Abs(AValue);
if absVal >= CBase then
Result.FDigits := TDigits.Create(absVal mod CBase, absVal div CBase)
else
Result.FDigits := TDigits.Create(absVal);
Result.FIsNegative := AValue < 0;
end
else begin
Result.FDigits := TDigits.Create(0, 1 << 31);
Result.FIsNegative := True;
end;
end;
class function TBigInt.FromHexadecimalString(const AValue: string): TBigInt;
begin
Result := FromHexOrBinString(AValue, 16);
end;
class function TBigInt.FromBinaryString(const AValue: string): TBigInt;
begin
Result := FromHexOrBinString(AValue, 2);
end;
{ Operators }
operator := (const A: Int64): TBigInt;
begin
Result := TBigInt.FromInt64(A);
end;
operator - (const A: TBigInt): TBigInt;
var
len: Integer;
begin
len := Length(A.FDigits);
if (len > 1) or (A.FDigits[0] > 0) then
begin
Result.FDigits := Copy(A.FDigits, 0, len);
Result.FIsNegative := not A.FIsNegative;
end
else
Result := TBigInt.Zero;
end;
operator + (const A, B: TBigInt): TBigInt;
begin
if A.IsNegative = B.IsNegative then
Result := TBigInt.AddAbsoluteValues(A, B, A.IsNegative)
else
Result := TBigInt.SubtractAbsoluteValues(A, B, A.IsNegative);
end;
operator - (const A, B: TBigInt): TBigInt;
begin
if A.IsNegative = B.IsNegative then
Result := TBigInt.SubtractAbsoluteValues(A, B, A.IsNegative)
else
Result := TBigInt.AddAbsoluteValues(A, B, A.IsNegative);
end;
operator * (const A, B: TBigInt): TBigInt;
begin
if (A = TBigInt.Zero) or (B = TBigInt.Zero) then
Result := TBigInt.Zero
else
Result := TBigInt.MultiplyAbsoluteValues(A, B, A.IsNegative <> B.IsNegative);
end;
operator shl(const A: TBigInt; const B: Integer): TBigInt;
var
i, j, digitShifts, bitShifts, reverseShift, len, newLength: Integer;
lastDigit: Cardinal;
begin
// Handles shift of zero.
if A = 0 then
begin
Result := TBigInt.Zero;
Exit;
end;
// Determines full digit shifts and bit shifts.
DivMod(B, CBitsPerDigit, digitShifts, bitShifts);
if bitShifts > 0 then
begin
reverseShift := CBitsPerDigit - bitShifts;
len := Length(A.FDigits);
lastDigit := A.FDigits[len - 1] >> reverseShift;
newLength := len + digitShifts;
if lastDigit = 0 then
SetLength(Result.FDigits, newLength)
else
SetLength(Result.FDigits, newLength + 1);
// Performs full digit shifts by shifting the access index j for A.FDigits.
Result.FDigits[digitShifts] := A.FDigits[0] << bitShifts;
j := 0;
for i := digitShifts + 1 to newLength - 1 do
begin
// Performs bit shifts.
Result.FDigits[i] := A.FDigits[j] >> reverseShift;
Inc(j);
Result.FDigits[i] := Result.FDigits[i] or (A.FDigits[j] << bitShifts);
end;
if lastDigit > 0 then
Result.FDigits[newLength] := lastDigit;
end
else begin
// Performs full digit shifts by copy if there are no bit shifts.
len := Length(A.FDigits);
SetLength(Result.FDigits, len + digitShifts);
for i := 0 to digitShifts - 1 do
Result.FDigits[i] := 0;
for i := 0 to len - 1 do
Result.FDigits[i + digitShifts] := A.FDigits[i];
end;
Result.FIsNegative := A.IsNegative;
end;
operator shr(const A: TBigInt; const B: Integer): TBigInt;
var
i, j, digitShifts, bitShifts, reverseShift, len, newLength: Integer;
lastDigit: Cardinal;
begin
// Handles shift of zero.
if A = 0 then
begin
Result := TBigInt.Zero;
Exit;
end;
// Determines full digit shifts and bit shifts.
DivMod(B, CBitsPerDigit, digitShifts, bitShifts);
// Handles shift to zero.
if digitShifts >= Length(A.FDigits) then
begin
Result := TBigInt.Zero;
Exit;
end;
if bitShifts > 0 then
begin
reverseShift := CBitsPerDigit - bitShifts;
len := Length(A.FDigits);
lastDigit := A.FDigits[len - 1] >> bitShifts;
newLength := len - digitShifts;
if lastDigit = 0 then
SetLength(Result.FDigits, newLength - 1)
else
SetLength(Result.FDigits, newLength);
// Performs full digit shifts by shifting the access index j for A.FDigits.
j := digitShifts;
for i := 0 to newLength - 2 do
begin
// Performs bit shifts.
Result.FDigits[i] := A.FDigits[j] >> bitShifts;
Inc(j);
Result.FDigits[i] := Result.FDigits[i] or (A.FDigits[j] << reverseShift);
end;
if lastDigit > 0 then
Result.FDigits[newLength - 1] := lastDigit;
end
else
// Performs full digit shifts by copy if there are no bit shifts.
Result.FDigits := Copy(A.FDigits, digitShifts, Length(A.FDigits) - digitShifts);
Result.FIsNegative := A.IsNegative;
end;
operator = (const A, B: TBigInt): Boolean;
begin
Result := A.CompareTo(B) = 0;
end;
operator <> (const A, B: TBigInt): Boolean;
begin
Result := A.CompareTo(B) <> 0;
end;
operator < (const A, B: TBigInt): Boolean;
begin
Result := A.CompareTo(B) < 0;
end;
operator <= (const A, B: TBigInt): Boolean;
begin
Result := A.CompareTo(B) <= 0;
end;
operator > (const A, B: TBigInt): Boolean;
begin
Result := A.CompareTo(B) > 0;
end;
operator >= (const A, B: TBigInt): Boolean;
begin
Result := A.CompareTo(B) >= 0;
end;
end.

220
UNumberTheory.pas
View File

@ -22,7 +22,7 @@ unit UNumberTheory;
interface interface
uses uses
Classes, SysUtils; Classes, SysUtils, Generics.Collections, Math;
type type
@ -34,6 +34,52 @@ type
class function LeastCommonMultiple(AValue1, AValue2: Int64): Int64; class function LeastCommonMultiple(AValue1, AValue2: Int64): Int64;
end; end;
TInt64Array = array of Int64;
{ TIntegerFactor }
TIntegerFactor = record
Factor: Int64;
Exponent: Byte;
end;
TIntegerFactors = specialize TList<TIntegerFactor>;
{ TIntegerFactorization }
TIntegerFactorization = class
public
class function PollardsRhoAlgorithm(const AValue: Int64): TInt64Array;
class function GetNormalized(constref AIntegerFactorArray: TInt64Array): TIntegerFactors;
end;
{ TDividersEnumerator }
TDividersEnumerator = class
private
FFactors: TIntegerFactors;
FCurrentExponents: array of Byte;
function GetCount: Integer;
public
constructor Create(constref AIntegerFactorArray: TInt64Array);
destructor Destroy; override;
function GetCurrent: Int64;
function MoveNext: Boolean;
procedure Reset;
property Current: Int64 read GetCurrent;
property Count: Integer read GetCount;
end;
{ TDividers }
TDividers = class
private
FFactorArray: TInt64Array;
public
constructor Create(constref AIntegerFactorArray: TInt64Array);
function GetEnumerator: TDividersEnumerator;
end;
implementation implementation
{ TNumberTheory } { TNumberTheory }
@ -58,5 +104,177 @@ begin
Result := (Abs(AValue1) div GreatestCommonDivisor(AValue1, AValue2)) * Abs(AValue2); Result := (Abs(AValue1) div GreatestCommonDivisor(AValue1, AValue2)) * Abs(AValue2);
end; end;
{ TIntegerFactorization }
// https://en.wikipedia.org/wiki/Pollard%27s_rho_algorithm
class function TIntegerFactorization.PollardsRhoAlgorithm(const AValue: Int64): TInt64Array;
var
primes: specialize TList<Int64>;
composites: specialize TStack<Int64>;
factor, n: Int64;
i: Integer;
function G(const AX, AC: Int64): Int64;
begin
Result := (AX * AX + AC) mod n;
end;
function FindFactor(const AStartValue, AC: Int64): Int64;
var
x, y, d: Int64;
begin
x := AStartValue;
y := x;
d := 1;
while d = 1 do
begin
x := G(x, AC);
y := G(G(y, AC), AC);
d := TNumberTheory.GreatestCommonDivisor(Abs(x - y), n);
end;
Result := d;
end;
begin
primes := specialize TList<Int64>.Create;
composites := specialize TStack<Int64>.Create;
n := Abs(AValue);
while (n and 1) = 0 do
begin
primes.Add(2);
n := n shr 1;
end;
composites.Push(n);
while composites.Count > 0 do
begin
n := composites.Pop;
i := 0;
repeat
factor := FindFactor(2 + (i + 1) div 2, 1 - i div 2);
if factor < n then
begin
composites.Push(factor);
composites.Push(n div factor);
end;
Inc(i);
until (factor < n) or (i > 3);
if factor = n then
primes.Add(factor);
end;
Result := primes.ToArray;
primes.Free;
composites.Free;
end;
class function TIntegerFactorization.GetNormalized(constref AIntegerFactorArray: TInt64Array): TIntegerFactors;
var
i: Integer;
factor: Int64;
normal: TIntegerFactor;
found: Boolean;
begin
Result := TIntegerFactors.Create;
for factor in AIntegerFactorArray do
begin
found := False;
for i := 0 to Result.Count - 1 do
if Result[i].Factor = factor then
begin
found := True;
normal := Result[i];
Inc(normal.Exponent);
Result[i] := normal;
Break;
end;
if not found then
begin
normal.Factor := factor;
normal.Exponent := 1;
Result.Add(normal);
end;
end;
end;
{ TDividersEnumerator }
function TDividersEnumerator.GetCount: Integer;
var
factor: TIntegerFactor;
begin
if FFactors.Count > 0 then
begin
Result := 1;
for factor in FFactors do
Result := Result * factor.Exponent;
Dec(Result);
end
else
Result := 0;
end;
constructor TDividersEnumerator.Create(constref AIntegerFactorArray: TInt64Array);
begin
FFactors := TIntegerFactorization.GetNormalized(AIntegerFactorArray);
SetLength(FCurrentExponents, FFactors.Count);
end;
destructor TDividersEnumerator.Destroy;
begin
FFactors.Free;
end;
function TDividersEnumerator.GetCurrent: Int64;
var
i: Integer;
begin
Result := 1;
for i := Low(FCurrentExponents) to High(FCurrentExponents) do
if FCurrentExponents[i] > 0 then
Result := Result * Round(Power(FFactors[i].Factor, FCurrentExponents[i]));
end;
function TDividersEnumerator.MoveNext: Boolean;
var
i: Integer;
begin
Result := False;
i := 0;
while (i <= High(FCurrentExponents)) and (FCurrentExponents[i] >= FFactors[i].Exponent) do
begin
FCurrentExponents[i] := 0;
Inc(i);
end;
if i <= High(FCurrentExponents) then
begin
Inc(FCurrentExponents[i]);
Result := True;
end;
end;
procedure TDividersEnumerator.Reset;
var
i: Integer;
begin
for i := Low(FCurrentExponents) to High(FCurrentExponents) do
FCurrentExponents[i] := 0;
end;
{ TDividers }
constructor TDividers.Create(constref AIntegerFactorArray: TInt64Array);
begin
FFactorArray := AIntegerFactorArray;
end;
function TDividers.GetEnumerator: TDividersEnumerator;
begin
Result := TDividersEnumerator.Create(FFactorArray);
end;
end. end.

297
UPolynomial.pas Normal file
View File

@ -0,0 +1,297 @@
{
Solutions to the Advent Of Code.
Copyright (C) 2024 Stefan Müller
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation, either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
}
unit UPolynomial;
{$mode ObjFPC}{$H+}
interface
uses
Classes, SysUtils, UBigInt;
type
TInt64Array = array of Int64;
{ TBigIntPolynomial }
TBigIntPolynomial = object
private
FCoefficients: array of TBigInt;
function GetDegree: Integer;
function GetCoefficient(const AIndex: Integer): TBigInt;
public
property Degree: Integer read GetDegree;
property Coefficient[const AIndex: Integer]: TBigInt read GetCoefficient;
function CalcValueAt(const AX: Int64): TBigInt;
function CalcSignVariations: Integer;
// Returns 2^n * f(x), given a polynomial f(x) and exponent n.
function ScaleByPowerOfTwo(const AExponent: Cardinal): TBigIntPolynomial;
// Returns f(s * x), given a polynomial f(x) and scale factor s.
function ScaleVariable(const AScaleFactor: TBigInt): TBigIntPolynomial;
// Returns f(2^n * x), given a polynomial f(x) and an exponent n.
function ScaleVariableByPowerOfTwo(const AExponent: Cardinal): TBigIntPolynomial;
// Returns f(x / 2), given a polynomial f(x).
function ScaleVariableByHalf: TBigIntPolynomial;
// Returns f(x + 1), given a polynomial f(x).
function TranslateVariableByOne: TBigIntPolynomial;
// Returns a polynomial with the reverse order of coefficients, i.e. the polynomial
// a_0 * x^n + a_1 * x^(n - 1) + ... + a_(n - 1) * x + a_n,
// given a polynomial
// a_n * x^n + a_(n - 1) * x^(n - 1) + ... + a_1 * x + a_0.
function RevertOrderOfCoefficients: TBigIntPolynomial;
// Returns a polynomial with all coefficents shifted down one position, and the constant term removed. This should
// only be used when the constant term is zero and is then equivalent to a division of polynomial f(x) by x.
function DivideByVariable: TBigIntPolynomial;
function IsEqualTo(const AOther: TBigIntPolynomial): Boolean;
function ToString: string;
class function Create(const ACoefficients: array of TBigInt): TBigIntPolynomial; static;
end;
{ Operators }
operator = (const A, B: TBigIntPolynomial): Boolean;
operator <> (const A, B: TBigIntPolynomial): Boolean;
implementation
{ TBigIntPolynomial }
function TBigIntPolynomial.GetDegree: Integer;
begin
Result := Length(FCoefficients) - 1;
end;
function TBigIntPolynomial.GetCoefficient(const AIndex: Integer): TBigInt;
begin
Result := FCoefficients[AIndex];
end;
function TBigIntPolynomial.CalcValueAt(const AX: Int64): TBigInt;
var
i: Integer;
begin
Result := TBigInt.Zero;
for i := High(FCoefficients) downto 0 do
Result := Result * AX + FCoefficients[i];
end;
function TBigIntPolynomial.CalcSignVariations: Integer;
var
current, last, i: Integer;
begin
Result := 0;
last := 0;
for i := 0 to Length(FCoefficients) - 1 do
begin
current := FCoefficients[i].Sign;
if (current <> 0) and (last <> current) then
begin
if last <> 0 then
Inc(Result);
last := current
end;
end;
end;
function TBigIntPolynomial.ScaleByPowerOfTwo(const AExponent: Cardinal): TBigIntPolynomial;
var
len, i: Integer;
begin
len := Length(FCoefficients);
SetLength(Result.FCoefficients, len);
for i := 0 to len - 1 do
Result.FCoefficients[i] := FCoefficients[i] << AExponent;
end;
function TBigIntPolynomial.ScaleVariable(const AScaleFactor: TBigInt): TBigIntPolynomial;
var
len, i: Integer;
factor: TBigInt;
begin
if AScaleFactor <> TBigInt.Zero then
begin
len := Length(FCoefficients);
SetLength(Result.FCoefficients, len);
Result.FCoefficients[0] := FCoefficients[0];
factor := AScaleFactor;
for i := 1 to len - 1 do begin
Result.FCoefficients[i] := FCoefficients[i] * factor;
factor := factor * AScaleFactor;
end;
end
else begin
SetLength(Result.FCoefficients, 1);
Result.FCoefficients[0] := TBigInt.Zero;
end;
end;
function TBigIntPolynomial.ScaleVariableByPowerOfTwo(const AExponent: Cardinal): TBigIntPolynomial;
var
len, i: Integer;
shift: Cardinal;
begin
len := Length(FCoefficients);
SetLength(Result.FCoefficients, len);
Result.FCoefficients[0] := FCoefficients[0];
shift := AExponent;
for i := 1 to len - 1 do begin
Result.FCoefficients[i] := FCoefficients[i] << shift;
Inc(shift, AExponent);
end;
end;
function TBigIntPolynomial.ScaleVariableByHalf: TBigIntPolynomial;
var
len, i: Integer;
begin
len := Length(FCoefficients);
SetLength(Result.FCoefficients, len);
Result.FCoefficients[0] := FCoefficients[0];
for i := 1 to len - 1 do
Result.FCoefficients[i] := FCoefficients[i] >> i;
end;
function TBigIntPolynomial.TranslateVariableByOne: TBigIntPolynomial;
var
len, i, j: Integer;
factors: array of Cardinal;
begin
len := Length(FCoefficients);
SetLength(Result.FCoefficients, len);
SetLength(factors, len);
for i := 0 to len - 1 do
begin
Result.FCoefficients[i] := TBigInt.Zero;
factors[i] := 1;
end;
// Calculates new coefficients.
for i := 0 to len - 1 do
begin
for j := 0 to len - i - 1 do
begin
if (i <> 0) and (j <> 0) then
factors[j] := factors[j] + factors[j - 1];
Result.FCoefficients[i] := Result.FCoefficients[i] + factors[j] * FCoefficients[j + i];
end;
end;
end;
function TBigIntPolynomial.RevertOrderOfCoefficients: TBigIntPolynomial;
var
len, skip, i: Integer;
begin
// Counts the trailing zeros to skip.
len := Length(FCoefficients);
skip := 0;
while (skip < len) and (FCoefficients[skip] = 0) do
Inc(skip);
// Copies the other coefficients in reverse order.
SetLength(Result.FCoefficients, len - skip);
for i := skip to len - 1 do
Result.FCoefficients[len - i - 1] := FCoefficients[i];
end;
function TBigIntPolynomial.DivideByVariable: TBigIntPolynomial;
var
len: Integer;
begin
len := Length(FCoefficients);
if len > 1 then
Result.FCoefficients := Copy(FCoefficients, 1, len - 1)
else begin
SetLength(Result.FCoefficients, 1);
Result.FCoefficients[0] := TBigInt.Zero;
end;
end;
function TBigIntPolynomial.IsEqualTo(const AOther: TBigIntPolynomial): Boolean;
var
i: Integer;
begin
if Length(FCoefficients) = Length(AOther.FCoefficients) then
begin
Result := True;
for i := 0 to Length(FCoefficients) - 1 do
if FCoefficients[i] <> AOther.FCoefficients[i] then
begin
Result := False;
Break;
end;
end
else
Result := False;
end;
function TBigIntPolynomial.ToString: string;
var
i: Integer;
begin
Result := FCoefficients[0].ToString;
for i := 1 to Length(FCoefficients) - 1 do
if i > 1 then
Result := Result + ' + ' + FCoefficients[i].ToString + ' * x^' + IntToStr(i)
else
Result := Result + ' + ' + FCoefficients[i].ToString + ' * x';
end;
class function TBigIntPolynomial.Create(const ACoefficients: array of TBigInt): TBigIntPolynomial;
var
high, i: integer;
begin
high := -1;
for i := Length(ACoefficients) - 1 downto 0 do
if ACoefficients[i] <> 0 then
begin
high := i;
Break;
end;
if high >= 0 then
begin
SetLength(Result.FCoefficients, high + 1);
for i := 0 to high do
Result.FCoefficients[i] := ACoefficients[i];
end
else begin
SetLength(Result.FCoefficients, 1);
Result.FCoefficients[0] := TBigInt.Zero;
end;
end;
{ Operators }
operator = (const A, B: TBigIntPolynomial): Boolean;
begin
Result := A.IsEqualTo(B);
end;
operator <> (const A, B: TBigIntPolynomial): Boolean;
begin
Result := not A.IsEqualTo(B);
end;
end.

205
UPolynomialRoots.pas Normal file
View File

@ -0,0 +1,205 @@
{
Solutions to the Advent Of Code.
Copyright (C) 2024 Stefan Müller
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation, either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
}
unit UPolynomialRoots;
{$mode ObjFPC}{$H+}
interface
uses
Classes, SysUtils, Generics.Collections, UPolynomial, UBigInt;
type
{ TIsolatingInterval }
// Represents an isolating interval of the form [C / 2^K, (C + H) / 2^K] in respect to [0, 1] or [A, B] in respect to
// [0, 2^boundexp], with A = C * 2^boundexp / 2^K and B = (C + H) * 2^boundexp / 2^K.
TIsolatingInterval = record
C: TBigInt;
K, H, BoundExp: Cardinal;
A, B: TBigInt;
end;
TIsolatingIntervals = specialize TList<TIsolatingInterval>;
TIsolatingIntervalArray = array of TIsolatingInterval;
{ TPolynomialRoots }
TPolynomialRoots = class
private
// Returns the exponent (base two) of an upper bound for the roots of the given polynomial, i.e. all real roots of
// the given polynomial are less or equal than 2^b, where b is the returned positive integer.
class function CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): Cardinal;
class function CreateIsolatingInterval(constref AC: TBigInt; const AK, AH: Cardinal; constref ABoundExp: Cardinal):
TIsolatingInterval;
public
// Returns root-isolating intervals for non-negative, non-multiple roots.
class function BisectIsolation(constref APolynomial: TBigIntPolynomial): TIsolatingIntervalArray;
// Returns root-isolating intervals for non-multiple roots in the interval [0, 2^boundexp].
class function BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal;
const AFindIntegers: Boolean = False): TIsolatingIntervalArray;
// Returns non-negative, non-multiple, integer roots in the interval [0, 2^boundexp].
class function BisectInteger(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal):
TBigIntArray;
end;
implementation
{ TPolynomialRoots }
class function TPolynomialRoots.CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): Cardinal;
var
i, sign: Integer;
an, ai, max: TBigInt;
numeratorBit, denominatorBit: Int64;
begin
// We need a_n > 0 here, so we use -sign(a_n) instead of actually flipping the polynomial.
// Sign is not 0 because a_n is not 0.
an := APolynomial.Coefficient[APolynomial.Degree];
sign := -an.Sign;
// This is a simplification of Cauchy's bound to avoid division and make it a power of two.
// https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Bounds_of_positive_real_roots
max := TBigInt.Zero;
for i := 0 to APolynomial.Degree - 1 do begin
ai := sign * APolynomial.Coefficient[i];
if max < ai then
max := ai;
end;
numeratorBit := max.GetMostSignificantBitIndex + 1;
denominatorBit := an.GetMostSignificantBitIndex;
Result := numeratorBit - denominatorBit;
end;
class function TPolynomialRoots.CreateIsolatingInterval(constref AC: TBigInt; const AK, AH: Cardinal;
constref ABoundExp: Cardinal): TIsolatingInterval;
begin
Result.C := AC;
Result.K := AK;
Result.H := AH;
Result.BoundExp := ABoundExp;
if ABoundExp >= AK then
begin
Result.A := AC << (ABoundExp - AK);
Result.B := (AC + AH) << (ABoundExp - AK);
end
else begin
Result.A := AC << (ABoundExp - AK);
Result.B := (AC + AH) << (ABoundExp - AK);
end;
end;
class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial): TIsolatingIntervalArray;
var
boundExp: Cardinal;
begin
boundExp := CalcUpperRootBound(APolynomial);
Result := BisectIsolation(APolynomial, boundExp);
end;
// This is adapted from https://en.wikipedia.org/wiki/Real-root_isolation#Bisection_method
class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal;
const AFindIntegers: Boolean): TIsolatingIntervalArray;
type
TWorkItem = record
C: TBigInt;
K: Cardinal;
P: TBigIntPolynomial;
end;
TWorkStack = specialize TStack<TWorkItem>;
var
item: TWorkItem;
stack: TWorkStack;
n, v: Integer;
varq: TBigIntPolynomial;
iso: TIsolatingIntervals;
begin
iso := TIsolatingIntervals.Create;
stack := TWorkStack.Create;
item.C := 0;
item.K := 0;
item.P := APolynomial.ScaleVariableByPowerOfTwo(ABoundExp);
stack.Push(item);
n := item.P.Degree;
while stack.Count > 0 do
begin
item := stack.Pop;
if item.P.Coefficient[0] = TBigInt.Zero then
begin
// Found an integer root at 0.
item.P := item.P.DivideByVariable;
Dec(n);
iso.Add(CreateIsolatingInterval(item.C, item.K, 0, ABoundExp));
end;
varq := item.P.RevertOrderOfCoefficients.TranslateVariableByOne;
v := varq.CalcSignVariations;
if (v > 1)
or ((v = 1) and AFindIntegers and (item.K < ABoundExp)) then
begin
// Bisects, first new work item is (2c, k + 1, 2^n * q(x/2)).
item.C := item.C << 1;
Inc(item.K);
item.P := item.P.ScaleVariableByHalf.ScaleByPowerOfTwo(n);
stack.Push(item);
// ... second new work item is (2c + 1, k + 1, 2^n * q((x+1)/2)).
item.C := item.C + 1;
item.P := item.P.TranslateVariableByOne;
stack.Push(item);
end
else if v = 1 then
begin
// Found isolating interval.
iso.Add(CreateIsolatingInterval(item.C, item.K, 1, ABoundExp));
end;
end;
Result := iso.ToArray;
iso.Free;
stack.Free;
end;
class function TPolynomialRoots.BisectInteger(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal):
TBigIntArray;
var
intervals: TIsolatingIntervalArray;
i: TIsolatingInterval;
r: specialize TList<TBigInt>;
value: Int64;
begin
// Calculates isolating intervals.
intervals := BisectIsolation(APolynomial, ABoundExp, True);
r := specialize TList<TBigInt>.Create;
for i in intervals do
if i.H = 0 then
r.Add(i.A)
else if i.A.TryToInt64(value) and (APolynomial.CalcValueAt(value) = 0) then
r.Add(value)
else if i.B.TryToInt64(value) and (APolynomial.CalcValueAt(value) = 0) then
r.Add(value);
Result := r.ToArray;
r.Free;
end;
end.

463
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,37 @@ unit UNeverTellMeTheOdds;
interface interface
uses uses
Classes, SysUtils, Generics.Collections, Math, USolver; Classes, SysUtils, Generics.Collections, Math, USolver, UNumberTheory, UBigInt, UPolynomial, UPolynomialRoots;
type type
{ THailstone } { THailstone }
THailstone = record THailstone = class
X, Y, Z: Int64; public
VX, VY, VZ: Integer; P0, P1, P2: Int64;
V0, V1, V2: Integer;
constructor Create(const ALine: string);
constructor Create;
end; end;
THailstones = specialize TList<THailstone>; THailstones = specialize TObjectList<THailstone>;
TInt64Array = array of Int64;
{ TNeverTellMeTheOdds } { TNeverTellMeTheOdds }
TNeverTellMeTheOdds = class(TSolver) TNeverTellMeTheOdds = class(TSolver)
private private
FMin, FMax: Int64; FMin, FMax: Int64;
FHailStones: THailstones; FHailstones: THailstones;
function AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean; 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 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;
@ -53,6 +64,26 @@ type
implementation 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 } { TNeverTellMeTheOdds }
function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean; function TNeverTellMeTheOdds.AreIntersecting(constref AHailstone1, AHailstone2: THailstone): Boolean;
@ -60,58 +91,432 @@ var
m1, m2, x, y: Double; m1, m2, x, y: Double;
begin begin
Result := False; Result := False;
m1 := AHailstone1.VY / AHailstone1.VX; m1 := AHailstone1.V1 / AHailstone1.V0;
m2 := AHailstone2.VY / AHailstone2.VX; m2 := AHailstone2.V1 / AHailstone2.V0;
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.P1 - m2 * AHailstone2.P0
- AHailstone1.P1 + m1 * AHailstone1.P0)
/ (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.V0) >= AHailstone1.P0 * Sign(AHailstone1.V0))
and (x * Sign(AHailstone2.VX) >= AHailstone2.X * Sign(AHailstone2.VX)) then and (x * Sign(AHailstone2.V0) >= AHailstone2.P0 * Sign(AHailstone2.V0))
then
begin begin
y := m1 * (x - AHailstone1.X) + AHailstone1.Y; y := m1 * (x - AHailstone1.P0) + AHailstone1.P1;
if (FMin <= y) and (y <= FMax) then if (FMin <= y) and (y <= FMax) then
Result := True Result := True
end; end;
end; 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..74] 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
// = k_14 * t_0 + k_15
// 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
// = k_16 * t_0^2 + k_17 * t_0 + k_18
// 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
// = k_19 * t_0^2 + k_20 * 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;
k[14] := k[11] * k[9] + k[8] * k[12];
k[15] := k[4] * k[9] + k[8] * k[6];
k[16] := k[11] * k[13] - k[10] * k[12];
k[17] := k[11] * k[7] + k[4] * k[13] + k[8] * k[5] - k[2] * k[9] - k[10] * k[6] - k[3] * k[12];
k[18] := k[4] * k[7] - k[3] * k[6];
k[19] := k[10] * k[5] + k[2] * k[13];
k[20] := 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_21 * t_0 + k_22
// 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_23 * t_0^2 + k_24 * t_0 + k_25
// 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_26 * t_0
// f_1^2
// = (k_16 * t_0^2 + k_17 * t_0 + k_18)^2
// = k_16^2 * t_0^4 + k_17^2 * t_0^2 + k_18^2 + 2 * k_16 * t_0^2 * k_17 * t_0 + 2 * k_16 * t_0^2 * k_18 + 2 * k_17 * t_0 * k_18
// = k_16^2 * t_0^4 + 2 * k_16 * k_17 * t_0^3 + (k_17^2 + 2 * k_16 * k_18) * t_0^2 + 2 * k_17 * k_18 * t_0 + k_18^2
// = k_16^2 * t_0^4 + k_27 * t_0^3 + k_29 * t_0^2 + k_30 * t_0 + k_18^2
// f_2^2
// = (k_14 * t_0 + k_15)^2
// = k_14^2 * t_0^2 + 2 * k_14 * k_15 * t_0 + k_15^2
// = k_14^2 * t_0^2 + k_31 * t_0 + k_15^2
// f_0 * f_2
// = (k_19 * t_0^2 + k_20 * t_0) * (k_14 * t_0 + k_15)
// = k_19 * k_14 * t_0^3 + (k_19 * k_15 + k_20 * k_14) * t_0^2 + k_20 * k_15 * t_0
// = k_33 * t_0^3 + k_34 * t_0^2 + k_35 * t_0
// f_1^2 + 4 * f_0 * f_2
// = k_16^2 * t_0^4 + k_27 * t_0^3 + k_29 * t_0^2 + k_30 * t_0 + k_18^2 + 4 * (k_33 * t_0^3 + k_34 * t_0^2 + k_35 * t_0)
// = k_37 * t_0^4 + k_75 * t_0^3 + k_76 * t_0^2 + k_77 * t_0 + k_59
// f_1^2 + 2 * f_0 * f_2
// = k_16^2 * t_0^4 + k_27 * t_0^3 + k_29 * t_0^2 + k_30 * t_0 + k_18^2 + 2 * (k_33 * t_0^3 + k_34 * t_0^2 + k_35 * t_0)
// = k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59
k[21] := AHailstone0.V0 * k[9] - AHailstone2.V0 * k[12];
k[22] := k[0] * k[9] - AHailstone2.V0 * k[6];
k[23] := AHailstone0.V0 * k[13];
k[24] := AHailstone0.V0 * k[7] + k[0] * k[13] - AHailstone2.V0 * k[5];
k[25] := k[0] * k[7];
k[26] := k[5] - k[7];
k[27] := 2 * k[16] * k[17];
k[28] := k[17] * k[17];
k[29] := k[28] + 2 * k[16] * k[18];
k[30] := 2 * k[17] * k[18];
k[31] := 2 * k[14] * k[15];
k[32] := k[14] * k[14];
k[33] := k[19] * k[14];
k[34] := k[19] * k[15] + k[20] * k[14];
k[35] := k[20] * k[15];
k[36] := k[15] * k[15];
k[37] := k[16] * k[16];
k[38] := k[27] + 2 * k[33];
k[39] := k[29] + 2 * k[34];
k[40] := k[30] + 2 * k[35];
k[41] := k[21] + AHailstone1.V0 * (k[12] - k[9]);
k[42] := k[22] + AHailstone1.V0 * k[6];
k[43] := k[23] - k[21] - AHailstone1.V0 * k[13] - (k[12] - k[9]) * AHailstone0.V0;
k[44] := k[24] - k[22] + AHailstone1.V0 * k[26] - (k[12] - k[9]) * k[1] - k[6] * AHailstone0.V0;
k[45] := k[25] - k[6] * k[1];
k[46] := k[43] * k[14] - k[41] * k[16];
k[47] := k[43] * k[15] + k[44] * k[14] - k[42] * k[16] - k[41] * k[17];
k[48] := k[44] * k[15] + k[45] * k[14] - k[42] * k[17] - k[41] * k[18];
k[49] := k[45] * k[15] - k[42] * k[18];
k[50] := k[43] * k[16];
k[51] := k[41] * k[37] - k[50] * k[14];
k[52] := k[43] * k[17] + k[44] * k[16];
k[53] := k[41] * k[38] + k[42] * k[37] - k[52] * k[14] - k[50] * k[15];
k[54] := k[43] * k[18] + k[44] * k[17] + k[45] * k[16];
k[55] := k[41] * k[39] + k[42] * k[38] - k[54] * k[14] - k[52] * k[15];
k[56] := k[44] * k[18] + k[45] * k[17];
k[57] := k[41] * k[40] + k[42] * k[39] - k[56] * k[14] - k[54] * k[15];
k[58] := k[45] * k[18];
k[59] := k[18] * k[18];
k[60] := k[41] * k[59] + k[42] * k[40] - k[58] * k[14] - k[56] * k[15];
k[61] := k[42] * k[59] - k[58] * k[15];
k[62] := k[13] * AHailstone0.V0 - k[23];
k[63] := 2 * k[32] * k[62];
k[64] := k[13] * k[1] - k[26] * AHailstone0.V0 - k[24];
k[65] := 2 * (k[31] * k[62] + k[32] * k[64]);
k[66] := - k[26] * k[1] - k[25];
k[67] := 2 * (k[36] * k[62] + k[31] * k[64] + k[32] * k[66]);
k[68] := 2 * (k[36] * k[64] + k[31] * k[66]);
k[69] := 2 * k[36] * k[66];
k[70] := k[51] + k[63];
k[71] := k[53] + k[65];
k[72] := k[55] + k[67];
k[73] := k[57] + k[68];
k[74] := k[60] + k[69];
// Unused, they are part of the polynomial inside the square root.
//k[75] := k[27] + 4 * k[33];
//k[76] := k[29] + 4 * k[34];
//k[77] := k[30] + 4 * k[35];
// Continuing calculations for equation 5.
// 0 = (k_21 * t_0 + k_22 + V_10 * ((k_12 - k_9) * t_0 + k_6)) * (k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59 -+ f_1 * sqrt(f_1^2 + 4 * f_0 * f_2))
// + (k_23 * t_0^2 + k_24 * t_0 + k_25 - t_0 * (k_21 * t_0 + k_22) - ((k_12 - k_9) * t_0 + k_6) * a_2 - V_10 * (k_13 * t_0^2 - k_26 * t_0)) * (- f_1 * f_2 +- f_2 * sqrt(f_1^2 + 4 * f_0 * f_2))
// - 2 * t_0 * (k_23 * t_0^2 + k_24 * t_0 + k_25) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2) + 2 * (k_13 * t_0^2 - k_26 * t_0) * a_2 * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2)
// 0 = (k_41 * t_0 + k_42) * (k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59 -+ f_1 * sqrt(f_1^2 + 4 * f_0 * f_2))
// + ((k_23 - k_21 - V_10 * k_13 - (k_12 - k_9) * V_00) * t_0^2 + (k_24 - k_22 + V_10 * k_26 - (k_12 - k_9) * k_1 - k_6 * V_00) * t_0 + k_25 - k_6 * k_1) * (- f_1 * f_2 +- f_2 * sqrt(f_1^2 + 4 * f_0 * f_2))
// - 2 * t_0 * (k_23 * t_0^2 + k_24 * t_0 + k_25) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2) + 2 * (k_13 * t_0^2 - k_26 * t_0) * a_2 * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2)
// 0 = (k_41 * t_0 + k_42) * (k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59)
// -+ (k_41 * t_0 + k_42) * f_1 * sqrt(f_1^2 + 4 * f_0 * f_2)
// + (k_43 * t_0^2 + k_44 * t_0 + k_45) * (- f_1 * f_2 +- f_2 * sqrt(f_1^2 + 4 * f_0 * f_2))
// - 2 * t_0 * (k_23 * t_0^2 + k_24 * t_0 + k_25) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2) + 2 * (k_13 * t_0^2 - k_26 * t_0) * a_2 * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2)
// 0 = (k_41 * t_0 + k_42) * (k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59)
// -+ (k_41 * t_0 + k_42) * f_1 * sqrt(f_1^2 + 4 * f_0 * f_2)
// - (k_43 * t_0^2 + k_44 * t_0 + k_45) * f_1 * f_2
// +- (k_43 * t_0^2 + k_44 * t_0 + k_45) * f_2 * sqrt(f_1^2 + 4 * f_0 * f_2)
// - 2 * t_0 * (k_23 * t_0^2 + k_24 * t_0 + k_25) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2) + 2 * (k_13 * t_0^2 - k_26 * t_0) * a_2 * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2)
// 0 = +- ((k_43 * t_0^2 + k_44 * t_0 + k_45) * f_2 - (k_41 * t_0 + k_42) * f_1) * sqrt(f_1^2 + 4 * f_0 * f_2)
// + (k_41 * t_0 + k_42) * (k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59)
// - (k_43 * t_0^2 + k_44 * t_0 + k_45) * f_1 * f_2
// - 2 * t_0 * (k_23 * t_0^2 + k_24 * t_0 + k_25) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2) + 2 * (k_13 * t_0^2 - k_26 * t_0) * a_2 * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2)
// 0 = +- ((k_43 * t_0^2 + k_44 * t_0 + k_45) * (k_14 * t_0 + k_15) - (k_41 * t_0 + k_42) * (k_16 * t_0^2 + k_17 * t_0 + k_18)) * sqrt(f_1^2 + 4 * f_0 * f_2)
// + (k_41 * t_0 + k_42) * (k_37 * t_0^4 + k_38 * t_0^3 + k_39 * t_0^2 + k_40 * t_0 + k_59)
// - (k_43 * t_0^2 + k_44 * t_0 + k_45) * (k_16 * t_0^2 + k_17 * t_0 + k_18) * (k_14 * t_0 + k_15)
// - 2 * t_0 * (k_23 * t_0^2 + k_24 * t_0 + k_25) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2) + 2 * (k_13 * t_0^2 - k_26 * t_0) * (V_00 * t_0 + k_1) * (k_14^2 * t_0^2 + k_31 * t_0 + k_15^2)
// 0 = +- (
// (k_43 * k_14 - k_41 * k_16) * t_0^3
// + (k_43 * k_15 + k_44 * k_14 - k_42 * k_16 - k_41 * k_17) * t_0^2
// + (k_44 * k_15 + k_45 * k_14 - k_42 * k_17 - k_41 * k_18) * t_0
// + k_45 * k_15 - k_42 * k_18
// ) * sqrt(f_1^2 + 4 * f_0 * f_2)
// + (k_41 * k_37 - k_43 * k_16 * k_14) * t_0^5
// + (k_41 * k_38 + k_42 * k_37 - k_43 * k_17 * k_14 - k_44 * k_16 * k_14 - k_43 * k_16 * k_15) * t_0^4
// + (k_41 * k_39 + k_42 * k_38 - k_43 * k_18 * k_14 - k_44 * k_17 * k_14 - k_45 * k_16 * k_14 - k_43 * k_17 * k_15 - k_44 * k_16 * k_15) * t_0^3
// + (k_41 * k_40 + k_42 * k_39 - k_44 * k_18 * k_14 - k_45 * k_17 * k_14 - k_43 * k_18 * k_15 - k_44 * k_17 * k_15 - k_45 * k_16 * k_15) * t_0^2
// + (k_41 * k_59 + k_42 * k_40 - k_45 * k_18 * k_14 - k_44 * k_18 * k_15 - k_45 * k_17 * k_15) * t_0
// + k_42 * k_59 - k_45 * k_18 * k_15
// + 2 * (k_13 * V_00 * k_14^2 - k_23 * k_14^2) * t_0^5
// + 2 * (k_13 * V_00 * k_31 + k_13 * k_1 * k_14^2 - k_26 * V_00 * k_14^2 - k_23 * k_31 - k_24 * k_14^2) * t_0^4
// + 2 * (k_13 * V_00 * k_15^2 + k_13 * k_1 * k_31 - k_26 * V_00 * k_31 - k_26 * k_1 * k_14^2 - k_23 * k_15^2 - k_24 * k_31 - k_25 * k_14^2) * t_0^3
// + 2 * (k_13 * k_1 * k_15^2 - k_26 * V_00 * k_15^2 - k_26 * k_1 * k_31 - k_24 * k_15^2 - k_25 * k_31) * t_0^2
// + 2 * (- k_26 * k_1 * k_15^2 - k_25 * k_15^2) * t_0
// 0 = +- (k_46 * t_0^3 + k_47 * t_0^2 + k_48 * t_0 + k_49) * sqrt(f_1^2 + 4 * f_0 * f_2)
// + (k_51 + k_63) * t_0^5 + (k_53 + k_65) * t_0^4 + (k_55 + k_67) * t_0^3 + (k_57 + k_68) * t_0^2 + (k_60 + k_69) * t_0 + k_61
// 0 = +- (k_46 * t_0^3 + k_47 * t_0^2 + k_48 * t_0 + k_49) * sqrt(k_37 * t_0^4 + k_75 * t_0^3 + k_76 * t_0^2 + k_77 * t_0 + k_59)
// + k_70 * t_0^5 + k_71 * t_0^4 + k_72 * t_0^3 + k_73 * t_0^2 + k_74 * t_0 + k_61
OPolynomial0 := TBigIntPolynomial.Create([k[61], k[74], k[73], k[72], k[71], k[70]]);
OPolynomial1 := TBigIntPolynomial.Create([k[49], k[48], k[47], k[46]]);
// 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_46 * t_0^3 + k_47 * t_0^2 + k_48 * t_0 + k_49) * sqrt(k_37 * t_0^4 + k_75 * t_0^3 + k_76 * t_0^2 + k_77 * t_0 + k_59) =
// k_70 * t_0^5 + k_71 * t_0^4 + k_72 * t_0^3 + k_73 * t_0^2 + k_74 * t_0 + k_61
// (k_46 * t_0^3 + k_47 * t_0^2 + k_48 * t_0 + k_49)^2 * (k_37 * t_0^4 + k_75 * t_0^3 + k_76 * t_0^2 + k_77 * t_0 + k_59) =
// (k_70 * t_0^5 + k_71 * t_0^4 + k_72 * t_0^3 + k_73 * t_0^2 + k_74 * t_0 + k_61)^2
// 0 =
// (k_46^2 * t_0^6
// + 2 * k_46 * k_47 * t_0^5
// + k_47^2 * t_0^4 + 2 * k_46 * k_48 * t_0^4
// + 2 * k_46 * k_49 * t_0^3 + 2 * k_47 * k_48 * t_0^3
// + k_48^2 * t_0^2 + 2 * k_47 * k_49 * t_0^2
// + 2 * k_48 * k_49 * t_0
// + k_49^2
// ) * (k_37 * t_0^4 + k_75 * t_0^3 + k_76 * t_0^2 + k_77 * t_0 + k_59)
// - k_70^2 * t_0^10
// - 2 * k_70 * k_71 * t_0^9
// - (k_71^2 + 2 * k_70 * k_72) * t_0^8
// - 2 * (k_70 * k_73 + k_71 * k_72) * t_0^7
// - (k_72^2 + 2 * k_70 * k_74 + 2 * k_71 * k_73) * t_0^6
// - 2 * (k_70 * k_61 + k_71 * k_74 + k_72 * k_73) * t_0^5
// - (k_73^2 + 2 * k_71 * k_61 + 2 * k_72 * k_74) * t_0^4
// - 2 * (k_72 * k_61 + k_73 * k_74) * t_0^3
// - (k_74^2 + 2 * k_73 * k_61) * t_0^2
// - 2 * k_74 * k_61 * t_0
// - k_61^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[78] := k[46] * k[46];
//k[79] := 2 * k[46] * k[47];
//k[80] := k[47] * k[47] + 2 * k[46] * k[48];
//k[81] := 2 * (k[46] * k[49] + k[47] * k[48]);
//k[82] := k[48] * k[48] + 2 * k[47] * k[49];
//k[83] := 2 * k[48] * k[49];
//k[84] := k[49] * k[49];
//ak[0] := k[59] * k[84] - k[61] * k[61];
//ak[1] := k[77] * k[84] + k[59] * k[83] - 2 * k[74] * k[61];
//ak[2] := k[76] * k[84] + k[77] * k[83] + k[59] * k[82] - k[74] * k[74] - 2 * k[73] * k[61];
//ak[3] := k[76] * k[83] + k[77] * k[82] + k[59] * k[81] + k[75] * k[84]
// - 2 * (k[72] * k[61] + k[73] * k[74]);
//ak[4] := k[37] * k[84] + k[76] * k[82] + k[77] * k[81] + k[59] * k[80] + k[75] * k[83] - k[73] * k[73]
// - 2 * (k[71] * k[61] + k[72] * k[74]);
//ak[5] := k[37] * k[83] + k[76] * k[81] + k[77] * k[80] + k[59] * k[79] + k[75] * k[82]
// - 2 * (k[70] * k[61] + k[71] * k[74] + k[72] * k[73]);
//ak[6] := k[37] * k[82] + k[76] * k[80] + k[77] * k[79] + k[59] * k[78] + k[75] * k[81] - k[72] * k[72]
// - 2 * (k[70] * k[74] + k[71] * k[73]);
//ak[7] := k[37] * k[81] + k[76] * k[79] + k[77] * k[78] + k[75] * k[80] - 2 * (k[70] * k[73] + k[71] * k[72]);
//ak[8] := k[37] * k[80] + k[75] * k[79] + k[76] * k[78] - k[71] * k[71] - 2 * k[70] * k[72];
//ak[9] := k[37] * k[79] + k[75] * k[78] - 2 * k[70] * k[71];
//ak[10] := k[37] * k[78] - k[70] * k[70];
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); 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: 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);
if FHailstones.Count >= 3 then
FPart2 := FindRockThrow(0, 1, 2);
end; end;
function TNeverTellMeTheOdds.GetDataFileName: string; function TNeverTellMeTheOdds.GetDataFileName: string;

16
tests/AdventOfCodeFPCUnit.lpi
View File

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

2
tests/AdventOfCodeFPCUnit.lpr
View File

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

1033
tests/UBigIntTestCases.pas Normal file
View File

File diff suppressed because it is too large Load Diff

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(856642398547748, 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;

138
tests/UPolynomialRootsTestCases.pas Normal file
View File

@ -0,0 +1,138 @@
{
Solutions to the Advent Of Code.
Copyright (C) 2024 Stefan Müller
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation, either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
}
unit UPolynomialRootsTestCases;
{$mode ObjFPC}{$H+}
interface
uses
Classes, SysUtils, fpcunit, testregistry, UPolynomial, UPolynomialRoots, UBigInt;
type
{ TPolynomialRootsTestCase }
TPolynomialRootsTestCase = class(TTestCase)
private
procedure AssertBisectIntervals(AIsolatingIntervals: TIsolatingIntervalArray; constref AExpectedRoots:
array of Cardinal);
procedure AssertBisectIntegers(ARoots: TBigIntArray; constref AExpectedRoots: array of Cardinal);
published
procedure TestBisectNoBound;
procedure TestBisectWithBound;
procedure TestBisectInteger;
end;
implementation
{ TPolynomialRootsTestCase }
procedure TPolynomialRootsTestCase.AssertBisectIntervals(AIsolatingIntervals: TIsolatingIntervalArray;
constref AExpectedRoots: array of Cardinal);
var
exp: Cardinal;
found: Boolean;
i, foundIndex: Integer;
begin
AssertEquals('Unexpected number of isolating intervals.', Length(AExpectedRoots), Length(AIsolatingIntervals));
for exp in AExpectedRoots do
begin
found := False;
for i := 0 to Length(AIsolatingIntervals) - 1 do
if (AIsolatingIntervals[i].A <= exp) and (exp <= AIsolatingIntervals[i].B) then
begin
found := True;
foundIndex := i;
Break;
end;
AssertTrue('No isolating interval for expected root ' + IntToStr(exp) + ' found.', found);
Delete(AIsolatingIntervals, foundIndex, 1);
end;
end;
procedure TPolynomialRootsTestCase.AssertBisectIntegers(ARoots: TBigIntArray; constref AExpectedRoots:
array of Cardinal);
var
exp: Cardinal;
found: Boolean;
i, foundIndex: Integer;
begin
AssertEquals('Unexpected number of integer roots.', Length(AExpectedRoots), Length(ARoots));
for exp in AExpectedRoots do
begin
found := False;
for i := 0 to Length(ARoots) - 1 do
if ARoots[i] = exp then
begin
found := True;
foundIndex := i;
Break;
end;
AssertTrue('Expected root ' + IntToStr(exp) + ' not found.', found);
Delete(ARoots, foundIndex, 1);
end;
end;
procedure TPolynomialRootsTestCase.TestBisectNoBound;
const
expRoots: array of Cardinal = (34000, 23017, 5);
var
a: TBigIntPolynomial;
r: TIsolatingIntervalArray;
begin
// y = 3 * (x - 34000) * (x - 23017) * (x - 5) * (x^2 - 19) * (x + 112)
// = 3 * x^6 - 170730 * x^5 + 2329429920 * x^4 + 251300082690 * x^3 - 1270471872603 * x^2 + 4774763204640 * x - 24979889760000
a := TBigIntPolynomial.Create([-24979889760000, 4774763204640, -1270471872603, 251300082690, 2329429920, -170730, 3]);
r := TPolynomialRoots.BisectIsolation(a);
AssertBisectIntervals(r, expRoots);
end;
procedure TPolynomialRootsTestCase.TestBisectWithBound;
const
expRoots: array of Cardinal = (23017, 5);
var
a: TBigIntPolynomial;
r: TIsolatingIntervalArray;
begin
// y = 3 * (x - 34000) * (x - 23017) * (x - 5) * (x^2 - 19) * (x + 112)
// = 3 * x^6 - 170730 * x^5 + 2329429920 * x^4 + 251300082690 * x^3 - 1270471872603 * x^2 + 4774763204640 * x - 24979889760000
a := TBigIntPolynomial.Create([-24979889760000, 4774763204640, -1270471872603, 251300082690, 2329429920, -170730, 3]);
r := TPolynomialRoots.BisectIsolation(a, 15);
AssertBisectIntervals(r, expRoots);
end;
procedure TPolynomialRootsTestCase.TestBisectInteger;
const
expRoots: array of Cardinal = (23017, 5);
var
a: TBigIntPolynomial;
r: TBigIntArray;
begin
// y = 3 * (x - 34000) * (x - 23017) * (x - 5) * (x^2 - 19) * (x + 112)
// = 3 * x^6 - 170730 * x^5 + 2329429920 * x^4 + 251300082690 * x^3 - 1270471872603 * x^2 + 4774763204640 * x - 24979889760000
a := TBigIntPolynomial.Create([-24979889760000, 4774763204640, -1270471872603, 251300082690, 2329429920, -170730, 3]);
r := TPolynomialRoots.BisectInteger(a, 15);
AssertBisectIntegers(r, expRoots);
end;
initialization
RegisterTest(TPolynomialRootsTestCase);
end.

187
tests/UPolynomialTestCases.pas Normal file
View File

@ -0,0 +1,187 @@
{
Solutions to the Advent Of Code.
Copyright (C) 2024 Stefan Müller
This program is free software: you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation, either version 3 of the License, or (at your option) any later
version.
This program is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with
this program. If not, see <http://www.gnu.org/licenses/>.
}
unit UPolynomialTestCases;
{$mode ObjFPC}{$H+}
interface
uses
Classes, SysUtils, fpcunit, testregistry, UPolynomial, UBigInt;
type
{ TBigIntPolynomialTestCase }
TBigIntPolynomialTestCase = class(TTestCase)
private
procedure TestCreateWithDegree(const ACoefficients: array of TBigInt; const ADegree: Integer);
published
procedure TestCreate;
procedure TestCreateDegreeZero;
procedure TestEqual;
procedure TestUnequalSameLength;
procedure TestUnequalDifferentLength;
procedure TestTrimLeadingZeros;
procedure TestCalcValueAt;
procedure TestSignVariations;
procedure TestScaleByPowerOfTwo;
procedure TestScaleVariable;
procedure TestScaleVariableByPowerOfTwo;
procedure TestTranslateVariableByOne;
procedure TestRevertOrderOfCoefficients;
procedure TestDivideByVariable;
end;
implementation
{ TBigIntPolynomialTestCase }
procedure TBigIntPolynomialTestCase.TestCreateWithDegree(const ACoefficients: array of TBigInt; const ADegree: Integer);
var
a: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create(ACoefficients);
AssertEquals('Degree of created polynomial incorrect.', ADegree, a.Degree);
end;
procedure TBigIntPolynomialTestCase.TestCreate;
begin
TestCreateWithDegree([992123, 7, 20, 4550022], 3);
end;
procedure TBigIntPolynomialTestCase.TestCreateDegreeZero;
begin
TestCreateWithDegree([4007], 0);
TestCreateWithDegree([], 0);
TestCreateWithDegree([0], 0);
end;
procedure TBigIntPolynomialTestCase.TestEqual;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([10, 7, 5, 1034]);
b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestUnequalSameLength;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([103, 7, 5, 10]);
b := TBigIntPolynomial.Create([1034, 7, 5, 10]);
AssertTrue('Polynomials are equal.', a <> b);
end;
procedure TBigIntPolynomialTestCase.TestUnequalDifferentLength;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([40000, 10, 7, 5, 1034]);
b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
AssertTrue('Polynomials are equal.', a <> b);
end;
procedure TBigIntPolynomialTestCase.TestTrimLeadingZeros;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([10, 7, 5, 1034, 0, 0]);
b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestCalcValueAt;
var
a: TBigIntPolynomial;
exp: TBigInt;
begin
a := TBigIntPolynomial.Create([80, 892477222, 0, 921556, 7303]);
exp:= TBigInt.FromInt64(16867124285);
AssertTrue('Polynomial evaluation unexpected.', a.CalcValueAt(15) = exp);
end;
procedure TBigIntPolynomialTestCase.TestSignVariations;
var
a: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([-10, 15, 0, 10, -20, -15, 0, 0, 5, 10, -10]);
AssertEquals(4, a.CalcSignVariations);
end;
procedure TBigIntPolynomialTestCase.TestScaleByPowerOfTwo;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([10, 7, 5, 1034]).ScaleByPowerOfTwo(7);
b := TBigIntPolynomial.Create([128 * 10, 128 * 7, 128 * 5, 128 * 1034]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestScaleVariable;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([10, 7, 5, 1034]).ScaleVariable(TBigInt.FromInt64(10));
b := TBigIntPolynomial.Create([10, 70, 500, 1034000]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestScaleVariableByPowerOfTwo;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([10, 7, 5, 1034]).ScaleVariableByPowerOfTwo(5);
b := TBigIntPolynomial.Create([10, 7 * 32, 5 * 32 * 32, 1034 * 32 * 32 * 32]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestTranslateVariableByOne;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([10, 7, 5, 1034]).TranslateVariableByOne;
b := TBigIntPolynomial.Create([1056, 3119, 3107, 1034]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestRevertOrderOfCoefficients;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([0, 10, 7, 5, 1034]).RevertOrderOfCoefficients;
b := TBigIntPolynomial.Create([1034, 5, 7, 10]);
AssertTrue('Polynomials are not equal.', a = b);
end;
procedure TBigIntPolynomialTestCase.TestDivideByVariable;
var
a, b: TBigIntPolynomial;
begin
a := TBigIntPolynomial.Create([0, 10, 7, 5, 1034]).DivideByVariable;
b := TBigIntPolynomial.Create([10, 7, 5, 1034]);
AssertTrue('Polynomials are not equal.', a = b);
end;
initialization
RegisterTest(TBigIntPolynomialTestCase);
end.