Updated bisection root finding algorithm and test case

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<TIsolatingInterval>;
+
   { 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<TWorkItem>;
 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.

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