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;