From 53827acf9b86ed04155d95f30fececb7a2829e2a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20M=C3=BCller?= Date: Thu, 4 Apr 2024 20:34:57 +0200 Subject: [PATCH] Added new unit for polynomial root finding algorithm --- AdventOfCode.lpi | 4 ++ UPolynomialRoots.pas | 72 +++++++++++++++++++++++++++++ tests/AdventOfCodeFPCUnit.lpi | 4 ++ tests/AdventOfCodeFPCUnit.lpr | 2 +- tests/UPolynomialRootsTestCases.pas | 72 +++++++++++++++++++++++++++++ tests/UPolynomialTestCases.pas | 11 ----- 6 files changed, 153 insertions(+), 12 deletions(-) create mode 100644 UPolynomialRoots.pas create mode 100644 tests/UPolynomialRootsTestCases.pas diff --git a/AdventOfCode.lpi b/AdventOfCode.lpi index 5174e60..3061d84 100644 --- a/AdventOfCode.lpi +++ b/AdventOfCode.lpi @@ -145,6 +145,10 @@ + + + + diff --git a/UPolynomialRoots.pas b/UPolynomialRoots.pas new file mode 100644 index 0000000..71ee652 --- /dev/null +++ b/UPolynomialRoots.pas @@ -0,0 +1,72 @@ +{ + Solutions to the Advent Of Code. + Copyright (C) 2024 Stefan Müller + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 3 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +} + +unit UPolynomialRoots; + +{$mode ObjFPC}{$H+} + +interface + +uses + Classes, SysUtils, UPolynomial, UBigInt; + +type + + { TRootIsolation } + + TRootIsolation = class + private + function CalcSimpleRootBound(constref APolynomial: TBigIntPolynomial): TBigInt; + public + function Bisect(constref APolynomial: TBigIntPolynomial): Int64; + end; + +implementation + +{ TRootIsolation } + +function TRootIsolation.CalcSimpleRootBound(constref APolynomial: TBigIntPolynomial): TBigInt; +var + i, sign: Integer; + a: TBigInt; +begin + // We need a_n > 0 here, so we use -sign(a_n) instead of actually flipping the polynomial. + // Sign is not 0 because a_n is not 0. + sign := -APolynomial.Coefficient[APolynomial.Degree].Sign; + + // This is a simplification of Cauchy's bound to avoid division. + // https://en.wikipedia.org/wiki/Geometrical_properties_of_polynomial_roots#Bounds_of_positive_real_roots + Result := TBigInt.Zero; + for i := 0 to APolynomial.Degree - 1 do begin + a := sign * APolynomial.Coefficient[i]; + if Result < a then + Result := a; + end; + Result := Result + 1; +end; + +function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial): Int64; +var + bound: TBigInt; + p: TBigIntPolynomial; +begin + bound := CalcSimpleRootBound(APolynomial); + p := APolynomial.ScaleVariable(bound); +end; + +end. + diff --git a/tests/AdventOfCodeFPCUnit.lpi b/tests/AdventOfCodeFPCUnit.lpi index dcef8d9..6d98efa 100644 --- a/tests/AdventOfCodeFPCUnit.lpi +++ b/tests/AdventOfCodeFPCUnit.lpi @@ -144,6 +144,10 @@ + + + + diff --git a/tests/AdventOfCodeFPCUnit.lpr b/tests/AdventOfCodeFPCUnit.lpr index eba70d8..eddfbea 100644 --- a/tests/AdventOfCodeFPCUnit.lpr +++ b/tests/AdventOfCodeFPCUnit.lpr @@ -9,7 +9,7 @@ uses UHotSpringsTestCases, UPointOfIncidenceTestCases, UParabolicReflectorDishTestCases, ULensLibraryTestCases, UFloorWillBeLavaTestCases, UClumsyCrucibleTestCases, ULavaductLagoonTestCases, UAplentyTestCases, UPulsePropagationTestCases, UStepCounterTestCases, USandSlabsTestCases, ULongWalkTestCases, - UNeverTellMeTheOddsTestCases, UPolynomialTestCases; + UNeverTellMeTheOddsTestCases, UPolynomialTestCases, UPolynomialRootsTestCases; {$R *.res} diff --git a/tests/UPolynomialRootsTestCases.pas b/tests/UPolynomialRootsTestCases.pas new file mode 100644 index 0000000..c7ae57c --- /dev/null +++ b/tests/UPolynomialRootsTestCases.pas @@ -0,0 +1,72 @@ +{ + Solutions to the Advent Of Code. + Copyright (C) 2024 Stefan Müller + + This program is free software: you can redistribute it and/or modify it under + the terms of the GNU General Public License as published by the Free Software + Foundation, either version 3 of the License, or (at your option) any later + version. + + This program is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS + FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along with + this program. If not, see . +} + +unit UPolynomialRootsTestCases; + +{$mode ObjFPC}{$H+} + +interface + +uses + Classes, SysUtils, fpcunit, testregistry, UPolynomial, UPolynomialRoots, UBigInt; + +type + + { TPolynomialRootsTestCase } + + TPolynomialRootsTestCase = class(TTestCase) + protected + FRootIsolation: TRootIsolation; + procedure SetUp; override; + procedure TearDown; override; + published + procedure TestBisectionRootIsolation; + end; + +implementation + +{ TPolynomialRootsTestCase } + +procedure TPolynomialRootsTestCase.SetUp; +begin + inherited SetUp; + FRootIsolation := TRootIsolation.Create; +end; + +procedure TPolynomialRootsTestCase.TearDown; +begin + FRootIsolation.Free; + inherited TearDown; +end; + +procedure TPolynomialRootsTestCase.TestBisectionRootIsolation; +var + a: TBigIntPolynomial; + r: Int64; +begin + // y = 3 * (x - 34000) * (x - 23017) * (x - 5) * (x^2 - 19) * (x + 112) + // = 3 * x^6 - 170730 * x^5 + 2329429920 * x^4 + 251300082690 * x^3 - 1270471872603 * x^2 + 4774763204640 * x - 24979889760000 + a := TBigIntPolynomial.Create([-24979889760000, 4774763204640, -1270471872603, 251300082690, 2329429920, -170730, 3]); + r := FRootIsolation.Bisect(a); + AssertEquals(0, r); +end; + +initialization + + RegisterTest(TPolynomialRootsTestCase); +end. + diff --git a/tests/UPolynomialTestCases.pas b/tests/UPolynomialTestCases.pas index 848caf3..b4d75a0 100644 --- a/tests/UPolynomialTestCases.pas +++ b/tests/UPolynomialTestCases.pas @@ -40,7 +40,6 @@ type procedure TestUnequalSameLength; procedure TestUnequalDifferentLength; procedure TestTrimLeadingZeros; - procedure TestBisectionRootIsolation; end; implementation @@ -111,16 +110,6 @@ begin AssertTrue('Polynomials are not equal.', a = b); end; -procedure TBigIntPolynomialTestCase.TestBisectionRootIsolation; -var - a: TBigIntPolynomial; -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]); - Fail('Not implemented'); -end; - initialization RegisterTest(TBigIntPolynomialTestCase);