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.
This commit is contained in:
Stefan Müller 2024-05-26 17:05:02 +02:00
parent 04e1702a2e
commit 8d4a5c2ed8
2 changed files with 32 additions and 20 deletions

View File

@ -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<TIsolatingInterval>;
@ -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

View File

@ -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;