From 53e3922654b2267a20377353630d56f572ba7058 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20M=C3=BCller?= Date: Fri, 24 May 2024 20:20:28 +0200 Subject: [PATCH] Updated bisection root finding algorithm and test case --- UPolynomialRoots.pas | 104 ++++++++++++++++++++++++---- tests/UPolynomialRootsTestCases.pas | 20 +++++- 2 files changed, 109 insertions(+), 15 deletions(-) diff --git a/UPolynomialRoots.pas b/UPolynomialRoots.pas index 71ee652..acf15f1 100644 --- a/UPolynomialRoots.pas +++ b/UPolynomialRoots.pas @@ -22,17 +22,29 @@ unit UPolynomialRoots; interface uses - Classes, SysUtils, UPolynomial, UBigInt; + 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 / 2^K, B / 2^K] in + // respect to [0, bound], with A = C * bound and B = (C + H) * bound. + TIsolatingInterval = record + C, K, H: Cardinal; + Bound, A, B: TBigInt; + end; + + TIsolatingIntervals = specialize TList; + { TRootIsolation } TRootIsolation = class private function CalcSimpleRootBound(constref APolynomial: TBigIntPolynomial): TBigInt; + function GetIsolatingInterval(const AC, AK, AH: Cardinal; constref ABound: TBigInt): TIsolatingInterval; public - function Bisect(constref APolynomial: TBigIntPolynomial): Int64; + function Bisect(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals; end; implementation @@ -42,30 +54,96 @@ implementation function TRootIsolation.CalcSimpleRootBound(constref APolynomial: TBigIntPolynomial): TBigInt; var i, sign: Integer; - a: TBigInt; + 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. - sign := -APolynomial.Coefficient[APolynomial.Degree].Sign; + an := APolynomial.Coefficient[APolynomial.Degree]; + sign := -an.Sign; - // This is a simplification of Cauchy's bound to avoid division. + // 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 - Result := TBigInt.Zero; + max := TBigInt.Zero; for i := 0 to APolynomial.Degree - 1 do begin - a := sign * APolynomial.Coefficient[i]; - if Result < a then - Result := a; + ai := sign * APolynomial.Coefficient[i]; + if max < ai then + max := ai; end; - Result := Result + 1; + numeratorBit := max.GetMostSignificantBitIndex + 1; + denominatorBit := an.GetMostSignificantBitIndex; + Result := TBigInt.One << (numeratorBit - denominatorBit); end; -function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial): Int64; +function TRootIsolation.GetIsolatingInterval(const AC, AK, AH: Cardinal; constref ABound: TBigInt): TIsolatingInterval; +begin + Result.C := AC; + Result.K := AK; + Result.H := AH; + Result.Bound := ABound; + Result.A := AC * ABound; + Result.B := (AC + AH) * ABound; +end; + +// This is adapted from +// https://en.wikipedia.org/wiki/Real-root_isolation#Bisection_method +function TRootIsolation.Bisect(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals; +type + TWorkItem = record + C, K: Cardinal; + P: TBigIntPolynomial; + end; + TWorkStack = specialize TStack; var bound: TBigInt; - p: TBigIntPolynomial; + item: TWorkItem; + stack: TWorkStack; + n, v: Integer; + varq: TBigIntPolynomial; begin + Result := TIsolatingIntervals.Create; + stack := TWorkStack.Create; + bound := CalcSimpleRootBound(APolynomial); - p := APolynomial.ScaleVariable(bound); + n := item.P.Degree; + + item.C := 0; + item.K := 0; + item.P := APolynomial.ScaleVariable(bound); + stack.Push(item); + + 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); + Result.Add(GetIsolatingInterval(item.C, item.K, 0, bound)); + end; + + varq := item.P.RevertOrderOfCoefficients.TranslateVariableByOne; + v := varq.CalcSignVariations; + if v = 1 then + begin + // Found isolating interval. + Result.Add(GetIsolatingInterval(item.C, item.K, 1, bound)); + end + else if v > 1 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; + end; + stack.Free; end; end. diff --git a/tests/UPolynomialRootsTestCases.pas b/tests/UPolynomialRootsTestCases.pas index c7ae57c..81cb8fb 100644 --- a/tests/UPolynomialRootsTestCases.pas +++ b/tests/UPolynomialRootsTestCases.pas @@ -54,15 +54,31 @@ begin end; procedure TPolynomialRootsTestCase.TestBisectionRootIsolation; +const + expRoots: array of Cardinal = (34000, 23017, 5); var + exp: Cardinal; a: TBigIntPolynomial; - r: Int64; + 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(0, r); + 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; end; initialization