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

Added bisection root finding algorithm with custom upper bound

This commit is contained in:
Stefan Müller 2024-05-24 20:47:52 +02:00
parent 53e3922654
commit baa1f8f31f
2 changed files with 57 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
UPolynomialRoots.pas
View File

@ -45,6 +45,7 @@ type
function GetIsolatingInterval(const AC, AK, AH: Cardinal; constref ABound: TBigInt): TIsolatingInterval;
public
function Bisect(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals;
function Bisect(constref APolynomial: TBigIntPolynomial; constref ABound: TBigInt): TIsolatingIntervals;
end;
implementation
@ -85,9 +86,17 @@ begin
Result.B := (AC + AH) * ABound;
end;
function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals;
var
bound: TBigInt;
begin
bound := CalcSimpleRootBound(APolynomial);
Result := Bisect(APolynomial, bound);
end;
// This is adapted from
// https://en.wikipedia.org/wiki/Real-root_isolation#Bisection_method
function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals;
function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial; constref ABound: TBigInt): TIsolatingIntervals;
type
TWorkItem = record
C, K: Cardinal;
@ -95,7 +104,6 @@ type
end;
TWorkStack = specialize TStack<TWorkItem>;
var
bound: TBigInt;
item: TWorkItem;
stack: TWorkStack;
n, v: Integer;
@ -104,12 +112,10 @@ begin
Result := TIsolatingIntervals.Create;
stack := TWorkStack.Create;
bound := CalcSimpleRootBound(APolynomial);
n := item.P.Degree;
item.C := 0;
item.K := 0;
item.P := APolynomial.ScaleVariable(bound);
item.P := APolynomial.ScaleVariable(ABound);
stack.Push(item);
while stack.Count > 0 do
@ -120,7 +126,7 @@ begin
// Found an integer root at 0.
item.P := item.P.DivideByVariable;
Dec(n);
Result.Add(GetIsolatingInterval(item.C, item.K, 0, bound));
Result.Add(GetIsolatingInterval(item.C, item.K, 0, ABound));
end;
varq := item.P.RevertOrderOfCoefficients.TranslateVariableByOne;
@ -128,7 +134,7 @@ begin
if v = 1 then
begin
// Found isolating interval.
Result.Add(GetIsolatingInterval(item.C, item.K, 1, bound));
Result.Add(GetIsolatingInterval(item.C, item.K, 1, ABound));
end
else if v > 1 then
begin

61
tests/UPolynomialRootsTestCases.pas
View File

@ -29,18 +29,43 @@ type
{ TPolynomialRootsTestCase }
TPolynomialRootsTestCase = class(TTestCase)
private
procedure AssertBisectResult(constref AIsolatingIntervals: TIsolatingIntervals; constref AExpectedRoots:
array of Cardinal);
protected
FRootIsolation: TRootIsolation;
procedure SetUp; override;
procedure TearDown; override;
published
procedure TestBisectionRootIsolation;
procedure TestBisectNoBound;
procedure TestBisectWithBound;
end;
implementation
{ TPolynomialRootsTestCase }
procedure TPolynomialRootsTestCase.AssertBisectResult(constref AIsolatingIntervals: TIsolatingIntervals; constref
AExpectedRoots: array of Cardinal);
var
exp: Cardinal;
ri: TIsolatingInterval;
found: Boolean;
begin
AssertEquals('Unexpected number of isolating intervals.', Length(AExpectedRoots), AIsolatingIntervals.Count);
for exp in AExpectedRoots do
begin
found := False;
for ri in AIsolatingIntervals do
if (ri.A <= exp) and (exp <= ri.B) then
begin
found := True;
Break;
end;
AssertTrue('No isolating interval for expected root ' + IntToStr(exp) + ' found.', found);
end;
end;
procedure TPolynomialRootsTestCase.SetUp;
begin
inherited SetUp;
@ -53,32 +78,34 @@ begin
inherited TearDown;
end;
procedure TPolynomialRootsTestCase.TestBisectionRootIsolation;
procedure TPolynomialRootsTestCase.TestBisectNoBound;
const
expRoots: array of Cardinal = (34000, 23017, 5);
var
exp: Cardinal;
a: TBigIntPolynomial;
r: TIsolatingIntervals;
ri: TIsolatingInterval;
found: Boolean;
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(Length(expRoots), r.Count);
for exp in expRoots do
begin
found := False;
for ri in r do
if (ri.A <= exp) and (exp <= ri.B) then
begin
found := True;
Break;
end;
AssertTrue('No isolating interval for expected root ' + IntToStr(exp) + ' found.', found);
end;
AssertBisectResult(r, expRoots);
r.Free;
end;
procedure TPolynomialRootsTestCase.TestBisectWithBound;
const
expRoots: array of Cardinal = (23017, 5);
var
a: TBigIntPolynomial;
r: TIsolatingIntervals;
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, TBigInt.One << 15);
AssertBisectResult(r, expRoots);
r.Free;
end;
initialization