diff --git a/UBigInt.pas b/UBigInt.pas index 421368d..5d7dd29 100644 --- a/UBigInt.pas +++ b/UBigInt.pas @@ -91,6 +91,8 @@ type class function FromBinaryString(const AValue: string): TBigInt; static; end; + TBigIntArray = array of TBigInt; + { Operators } operator := (const A: Int64): TBigInt; diff --git a/UPolynomialRoots.pas b/UPolynomialRoots.pas index fa1280c..a05c9fa 100644 --- a/UPolynomialRoots.pas +++ b/UPolynomialRoots.pas @@ -52,8 +52,11 @@ type // Returns root-isolating intervals for non-negative, non-multiple roots. class function BisectIsolation(constref APolynomial: TBigIntPolynomial): TIsolatingIntervalArray; // Returns root-isolating intervals for non-multiple roots in the interval [0, 2^boundexp]. - class function BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal): - TIsolatingIntervalArray; + class function BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal; + const AFindIntegers: Boolean = False): TIsolatingIntervalArray; + // Returns non-negative, non-multiple, integer roots in the interval [0, 2^boundexp]. + class function BisectInteger(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal): + TBigIntArray; end; implementation @@ -112,8 +115,8 @@ end; // This is adapted from // https://en.wikipedia.org/wiki/Real-root_isolation#Bisection_method -class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal): - TIsolatingIntervalArray; +class function TPolynomialRoots.BisectIsolation(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal; + const AFindIntegers: Boolean): TIsolatingIntervalArray; type TWorkItem = record C, K: Cardinal; @@ -149,12 +152,11 @@ begin varq := item.P.RevertOrderOfCoefficients.TranslateVariableByOne; v := varq.CalcSignVariations; - if v = 1 then - begin - // Found isolating interval. - iso.Add(CreateIsolatingInterval(item.C, item.K, 1, ABoundExp)); - end - else if v > 1 then + //WriteLn; + //WriteLn('var((x+1)^n*q(1/(x+1))): ', v); + + if (v > 1) + or ((v = 1) and AFindIntegers and (item.K < ABoundExp)) then begin // Bisects, first new work item is (2c, k + 1, 2^n * q(x/2)). item.C := item.C << 1; @@ -165,6 +167,12 @@ begin item.C := item.C + 1; item.P := item.P.TranslateVariableByOne; stack.Push(item); + end + else if v = 1 then + begin + // Found isolating interval. + //WriteLn('Found isolating interval (' + IntToStr(item.C) + ', ' + IntToStr(item.K) + ', 1).'); + iso.Add(CreateIsolatingInterval(item.C, item.K, 1, ABoundExp)); end; end; Result := iso.ToArray; @@ -172,5 +180,29 @@ begin stack.Free; end; +class function TPolynomialRoots.BisectInteger(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal): + TBigIntArray; +var + intervals: TIsolatingIntervalArray; + i: TIsolatingInterval; + r: specialize TList; + value: Int64; +begin + // Calculates isolating intervals. + intervals := BisectIsolation(APolynomial, ABoundExp, True); + r := specialize TList.Create; + + for i in intervals do + if i.H = 0 then + r.Add(i.A) + else if i.A.TryToInt64(value) and (APolynomial.CalcValueAt(value) = 0) then + r.Add(value) + else if i.B.TryToInt64(value) and (APolynomial.CalcValueAt(value) = 0) then + r.Add(value); + + Result := r.ToArray; + r.Free; +end; + end. diff --git a/tests/UPolynomialRootsTestCases.pas b/tests/UPolynomialRootsTestCases.pas index 91c29d3..91ae105 100644 --- a/tests/UPolynomialRootsTestCases.pas +++ b/tests/UPolynomialRootsTestCases.pas @@ -32,9 +32,11 @@ type private procedure AssertBisectIntervals(AIsolatingIntervals: TIsolatingIntervalArray; constref AExpectedRoots: array of Cardinal); + procedure AssertBisectIntegers(ARoots: TBigIntArray; constref AExpectedRoots: array of Cardinal); published procedure TestBisectNoBound; procedure TestBisectWithBound; + procedure TestBisectInteger; end; implementation @@ -64,6 +66,29 @@ begin end; end; +procedure TPolynomialRootsTestCase.AssertBisectIntegers(ARoots: TBigIntArray; constref AExpectedRoots: + array of Cardinal); +var + exp: Cardinal; + found: Boolean; + i, foundIndex: Integer; +begin + AssertEquals('Unexpected number of integer roots.', Length(AExpectedRoots), Length(ARoots)); + for exp in AExpectedRoots do + begin + found := False; + for i := 0 to Length(ARoots) - 1 do + if ARoots[i] = exp then + begin + found := True; + foundIndex := i; + Break; + end; + AssertTrue('Expected root ' + IntToStr(exp) + ' not found.', found); + Delete(ARoots, foundIndex, 1); + end; +end; + procedure TPolynomialRootsTestCase.TestBisectNoBound; const expRoots: array of Cardinal = (34000, 23017, 5); @@ -92,6 +117,20 @@ begin AssertBisectIntervals(r, expRoots); end; +procedure TPolynomialRootsTestCase.TestBisectInteger; +const + expRoots: array of Cardinal = (23017, 5); +var + a: TBigIntPolynomial; + r: TBigIntArray; +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.BisectInteger(a, 15); + AssertBisectIntegers(r, expRoots); +end; + initialization RegisterTest(TPolynomialRootsTestCase);