diff --git a/UPolynomialRoots.pas b/UPolynomialRoots.pas index acf15f1..a4b029a 100644 --- a/UPolynomialRoots.pas +++ b/UPolynomialRoots.pas @@ -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; 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 diff --git a/tests/UPolynomialRootsTestCases.pas b/tests/UPolynomialRootsTestCases.pas index 81cb8fb..3ae8ffe 100644 --- a/tests/UPolynomialRootsTestCases.pas +++ b/tests/UPolynomialRootsTestCases.pas @@ -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