Added new unit for polynomial root finding algorithm

This commit is contained in:
Stefan Müller 2024-04-04 20:34:57 +02:00
parent 3f0d13248d
commit eb0eb8b615
6 changed files with 153 additions and 12 deletions

View File

@ -145,6 +145,10 @@
<Filename Value="UPolynomial.pas"/>
<IsPartOfProject Value="True"/>
</Unit>
<Unit>
<Filename Value="UPolynomialRoots.pas"/>
<IsPartOfProject Value="True"/>
</Unit>
</Units>
</ProjectOptions>
<CompilerOptions>

72
UPolynomialRoots.pas Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
}
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.

View File

@ -144,6 +144,10 @@
<Filename Value="UPolynomialTestCases.pas"/>
<IsPartOfProject Value="True"/>
</Unit>
<Unit>
<Filename Value="UPolynomialRootsTestCases.pas"/>
<IsPartOfProject Value="True"/>
</Unit>
</Units>
</ProjectOptions>
<CompilerOptions>

View File

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

View File

@ -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 <http://www.gnu.org/licenses/>.
}
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.

View File

@ -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);