{
  Solutions to the Advent Of Code.
  Copyright (C) 2024  Stefan Müller

  This program is free software: you can redistribute it and/or modify it under
  the terms of the GNU General Public License as published by the Free Software
  Foundation, either version 3 of the License, or (at your option) any later
  version.

  This program is distributed in the hope that it will be useful, but WITHOUT
  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.

  You should have received a copy of the GNU General Public License along with
  this program.  If not, see <http://www.gnu.org/licenses/>.
}

unit UPolynomialRoots;

{$mode ObjFPC}{$H+}

interface

uses
  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, B] in respect to
  // [0, 2^boundexp], with A = C * 2^boundexp / 2^K and B = (C + H) * 2^boundexp / 2^K.
  TIsolatingInterval = record
    C, K, H, BoundExp: Cardinal;
    A, B: TBigInt;
  end;

  TIsolatingIntervals = specialize TList<TIsolatingInterval>;

  TIsolatingIntervalArray = array of TIsolatingInterval;

  { TPolynomialRoots }

  TPolynomialRoots = class
  private
    // 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): TIsolatingIntervalArray;
    // Returns root-isolating intervals for non-multiple roots in the interval [0, 2^boundexp].
    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

{ TPolynomialRoots }

class function TPolynomialRoots.CalcUpperRootBound(constref APolynomial: TBigIntPolynomial): Cardinal;
var
  i, sign: Integer;
  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.
  an := APolynomial.Coefficient[APolynomial.Degree];
  sign := -an.Sign;

  // 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
  max := TBigInt.Zero;
  for i := 0 to APolynomial.Degree - 1 do begin
    ai := sign * APolynomial.Coefficient[i];
    if max < ai then
      max := ai;
  end;
  numeratorBit := max.GetMostSignificantBitIndex + 1;
  denominatorBit := an.GetMostSignificantBitIndex;
  Result := numeratorBit - denominatorBit;
end;

class function TPolynomialRoots.CreateIsolatingInterval(const AC, AK, AH: Cardinal; constref ABoundExp: Cardinal):
  TIsolatingInterval;
begin
  Result.C := AC;
  Result.K := AK;
  Result.H := AH;
  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): TIsolatingIntervalArray;
var
  boundExp: Cardinal;
begin
  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 ABoundExp: Cardinal;
  const AFindIntegers: Boolean): TIsolatingIntervalArray;
type
  TWorkItem = record
    C, K: Cardinal;
    P: TBigIntPolynomial;
  end;
  TWorkStack = specialize TStack<TWorkItem>;
var
  item: TWorkItem;
  stack: TWorkStack;
  n, v: Integer;
  varq: TBigIntPolynomial;
  iso: TIsolatingIntervals;
begin
  iso := TIsolatingIntervals.Create;
  stack := TWorkStack.Create;

  item.C := 0;
  item.K := 0;
  item.P := APolynomial.ScaleVariableByPowerOfTwo(ABoundExp);
  stack.Push(item);
  n := item.P.Degree;

  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);
      iso.Add(CreateIsolatingInterval(item.C, item.K, 0, ABoundExp));
    end;

    varq := item.P.RevertOrderOfCoefficients.TranslateVariableByOne;
    v := varq.CalcSignVariations;
    //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;
      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
    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;
  iso.Free;
  stack.Free;
end;

class function TPolynomialRoots.BisectInteger(constref APolynomial: TBigIntPolynomial; constref ABoundExp: Cardinal):
  TBigIntArray;
var
  intervals: TIsolatingIntervalArray;
  i: TIsolatingInterval;
  r: specialize TList<TBigInt>;
  value: Int64;
begin
  // Calculates isolating intervals.
  intervals := BisectIsolation(APolynomial, ABoundExp, True);
  r := specialize TList<TBigInt>.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.