From 8d4a5c2ed8b1fed523d5481e6bd7dc5a953cdee4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Stefan=20M=C3=BCller?= Date: Sun, 26 May 2024 17:05:02 +0200 Subject: [PATCH] Changed root finder to use base-2 exponent of bound Since the root isolation algorithm uses a power of two as bound anyway, it makes sense to use it's exponent in method interfaces and throughout the algorithm. This simplifies multiplications to cheap shifts and will make it easier to detect when the isolating interval size is 1. Also added some method documentation. --- UPolynomialRoots.pas | 50 ++++++++++++++++++----------- tests/UPolynomialRootsTestCases.pas | 2 +- 2 files changed, 32 insertions(+), 20 deletions(-) diff --git a/UPolynomialRoots.pas b/UPolynomialRoots.pas index 2a522f0..6853b80 100644 --- a/UPolynomialRoots.pas +++ b/UPolynomialRoots.pas @@ -29,10 +29,10 @@ type { TIsolatingInterval } // Represents an isolating interval of the form [C / 2^K, (C + H) / 2^K] in respect to [0, 1] or [A, B] in respect to - // [0, bound], with A = C * bound / 2^K and B = (C + H) * bound / 2^K. + // [0, 2^boundexp], with A = C * 2^boundexp / 2^K and B = (C + H) * 2^boundexp / 2^K. TIsolatingInterval = record - C, K, H: Cardinal; - Bound, A, B: TBigInt; + C, K, H, BoundExp: Cardinal; + A, B: TBigInt; end; TIsolatingIntervals = specialize TList; @@ -41,11 +41,16 @@ type TPolynomialRoots = class private - class function CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): TBigInt; - class function CreateIsolatingInterval(const AC, AK, AH: Cardinal; constref ABound: TBigInt): TIsolatingInterval; + // Returns the exponent (base two) of an upper bound for the roots of the given polynomial, i.e. all real roots of + // the given polynomial are less or equal than 2^b, where b is the returned positive integer. + class function CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): Cardinal; + class function CreateIsolatingInterval(const AC, AK, AH: Cardinal; constref ABoundExp: Cardinal): + TIsolatingInterval; public + // Returns root-isolating intervals for non-negative, non-multiple roots. class function BisectIsolation(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals; - class function BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABound: TBigInt): + // Returns root-isolating intervals for non-multiple roots in the interval [0, 2^boundexp]. + class function BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal): TIsolatingIntervals; end; @@ -53,7 +58,7 @@ implementation { TPolynomialRoots } -class function TPolynomialRoots.CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): TBigInt; +class function TPolynomialRoots.CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): Cardinal; var i, sign: Integer; an, ai, max: TBigInt; @@ -74,31 +79,38 @@ begin end; numeratorBit := max.GetMostSignificantBitIndex + 1; denominatorBit := an.GetMostSignificantBitIndex; - Result := TBigInt.One << (numeratorBit - denominatorBit); + Result := numeratorBit - denominatorBit; end; -class function TPolynomialRoots.CreateIsolatingInterval(const AC, AK, AH: Cardinal; constref ABound: TBigInt): +class function TPolynomialRoots.CreateIsolatingInterval(const AC, AK, AH: Cardinal; constref ABoundExp: Cardinal): TIsolatingInterval; begin Result.C := AC; Result.K := AK; Result.H := AH; - Result.Bound := ABound; - Result.A := (AC * ABound) >> AK; - Result.B := ((AC + AH) * ABound) >> AK; + Result.BoundExp := ABoundExp; + if ABoundExp >= AK then + begin + Result.A := AC << (ABoundExp - AK); + Result.B := (AC + AH) << (ABoundExp - AK); + end + else begin + Result.A := AC << (ABoundExp - AK); + Result.B := (AC + AH) << (ABoundExp - AK); + end; end; class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial): TIsolatingIntervals; var - bound: TBigInt; + boundExp: Cardinal; begin - bound := CalcUpperRootBound(APolynomial); - Result := BisectIsolation(APolynomial, bound); + boundExp := CalcUpperRootBound(APolynomial); + Result := BisectIsolation(APolynomial, boundExp); end; // This is adapted from // https://en.wikipedia.org/wiki/Real-root_isolation#Bisection_method -class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABound: TBigInt): +class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal): TIsolatingIntervals; type TWorkItem = record @@ -117,7 +129,7 @@ begin item.C := 0; item.K := 0; - item.P := APolynomial.ScaleVariable(ABound); + item.P := APolynomial.ScaleVariableByPowerOfTwo(ABoundExp); stack.Push(item); n := item.P.Degree; @@ -129,7 +141,7 @@ begin // Found an integer root at 0. item.P := item.P.DivideByVariable; Dec(n); - Result.Add(CreateIsolatingInterval(item.C, item.K, 0, ABound)); + Result.Add(CreateIsolatingInterval(item.C, item.K, 0, ABoundExp)); end; varq := item.P.RevertOrderOfCoefficients.TranslateVariableByOne; @@ -137,7 +149,7 @@ begin if v = 1 then begin // Found isolating interval. - Result.Add(CreateIsolatingInterval(item.C, item.K, 1, ABound)); + Result.Add(CreateIsolatingInterval(item.C, item.K, 1, ABoundExp)); end else if v > 1 then begin diff --git a/tests/UPolynomialRootsTestCases.pas b/tests/UPolynomialRootsTestCases.pas index 03faa50..cfd9e47 100644 --- a/tests/UPolynomialRootsTestCases.pas +++ b/tests/UPolynomialRootsTestCases.pas @@ -89,7 +89,7 @@ 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 := TPolynomialRoots.BisectIsolation(a, TBigInt.One << 15); + r := TPolynomialRoots.BisectIsolation(a, 15); AssertBisectResult(r, expRoots); r.Free; end;