unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    procedure Button1Click(Sender: TObject);
    procedure FormCreate(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation
uses Math, VarCmplx;

{$R *.dfm}

var
  _zero, _i: Variant;

const
  kEpsilon = 1E-6;

function laguerre(degree: integer;
         const a: array of Variant; var zp: Variant): boolean;
const
  f: array[0 .. 7] of double =
    (1.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88);
var
  z: Variant;
  sn, snm: double;
  j, k: longint;

  ca, caca, cb, c1, dz, cc, cp, cm, zz: Variant;
  alpha, beta, gamma: Variant;

  kk: integer;
begin
  z := zp;

  sn := degree;
  snm := degree - 1;

  for k := 0 to pred(10000) do begin // Not more than ten thousand iterations.
    alpha := a[degree];
    beta := VarComplexCreate(0.0, 0.0);
    gamma := VarComplexCreate(0.0, 0.0);

    for j := degree - 1 downto 0 do begin // Horner's scheme.
      gamma := (z * gamma) + beta;
      beta  := (z * beta)  + alpha;
      alpha := (z * alpha) + a[j];
    end;

    if alpha = _zero then begin
        zp := z;
        result := false; exit;
    end;

    (* done
    if (alpha == 0.0)
        {
        }
    *)
    ca   := (-1.0 * beta) / alpha;
    caca := ca * ca;
    cb   := caca - (2.0 * gamma / alpha);
    c1   := VarComplexSqrt(snm * ((sn * cb) - caca));
    cp   := ca + c1;
    cm   := ca - c1;
    if (VarComplexAbsSqr(cm) > VarComplexAbsSqr(cp)) then      // norm() is faster than abs().
        cp := cm;

    if (cp = _zero) then                // Re-init trick from 'Numerical Recipes'.
        dz := (1.0 + VarComplexAbs(z)) * VarComplexFromPolar(1.0, k)
    else
        dz := sn / cp;

    // short kk = k % 9;               // Once every 9 iterations.
    kk := k mod 9;
    if (kk = 0) then
        dz := dz * f[(k div 9) and 7];       // Cycle breaking trick (but f[0]=1.0).

    z := z + dz;
    if ((kk <> 0) and (VarComplexAbsSqr(dz) < 1e-31)) then begin
      // OK, converged: almost no change.
      // THIS 1e-31 IS RELATED TO kEPSILON !
      zp := z;
      result := false; exit;
    end;

    if (k = 2000) then begin
      z := z * -1.0;           // Kick z around when it
    end
    else
      if (k = 4000) then begin
        z := VarComplexConjugate(z);     // is taking too long.
      end
      else
        if (k = 6000) then begin
          z := -1 * VarComplexConjugate(z);
        end
        else
          if (k = 8000) then begin
            z := z * _i;
          end;
  end;

  laguerre := true;
end;

type
  parrtype = ^arrtype;
  arrtype = array[0 .. maxint div (2 * sizeof(variant))] of variant;


function factor(coeffs_size: integer;
         const coeffs: array of variant;
         var roots: array of variant;
         var num_roots: integer): integer;
const
  n_attempts = 4;
var
  guess: array[0 .. pred(n_attempts)] of variant;
  p, z, a2, d: Variant;
  attempt: integer;
  magnitude, angle, angle_inc: double;
  degree, degreeH, e, n: integer;
  coeff: integer;

  cc: parrtype;
  c: integer;
begin

  guess[0] := VarComplexCreate(8.0, 0.0);
  guess[1] := VarComplexCreate(-6.0, 0.0);
  guess[2] := VarComplexCreate(0.0, 6.0);
  guess[3] := VarComplexCreate(0.0, -8.0);

  attempt := 0;
    coeff := 0;
    degree := coeffs_size - 1;
    num_roots := 0;
    e := 0;

    // Find the ACTUAL degree of the polynomial. This guarentees last coeff is non-zero.

    while((degree > 0) and (coeffs[degree] = _zero)) do // Use kEPSILON instead of 0.0?
        dec(degree);

    if (degree < 1) then begin                // There must at least be 2 coefficients, of which
        factor := 1; exit;                    // only the first may be zero.
    end;

    // Factor-out roots at the origin.
    while ((degree > 0) and (coeffs[coeff] = _zero)) do begin // Use kEPSILON instead of 0.0?
        roots[num_roots] := _zero;       // Insert root at 0.
        inc(num_roots);

        inc(coeff);                    // Never reuse this coeff (INC ITERATOR).
        dec(degree);
    end;

    degreeH := degree; // Remember for retries.
    // From here on, first and last coeff are non-zero.
    if (degree > 0) then begin        // Copy input array so we can write in it.
                                      // (Only necessary for Laguerre, actually).

        // cc = c = new complex<double>[degree + 1]; // MAY THROW AN EXCEPTION.
        getmem(cc, sizeof(variant) * (degree + 1));
        c := 0;
        for n := 0 to degree do begin
          cc^[n] := coeffs[coeff + n];
        end;
        (*
        for (int n = 0; n <= degree; n++)
            c[n] = coeff[n];
        *)

    end; // if

    while (degree > 2) do begin
        // It is perhaps a bit stupid to investigate this EACH time: the chance that
        // such a simple form appears in between Laguerre iterations is almost nil(?).
        // See whether all coeffs in between are zero.
        n := 1;
        while (n < degree) and (cc^[c + n] = _zero) do begin // Use kEPSILON instead of 0.0?
          inc(n);
        end;

                                        // Then we have a complex    88
        if (n = degree) then begin      // binomial in this form:   x  = -a0 / a88
                                        // Now we have 88 solutions all at once:
            p := -1 * cc^[c] / cc^[degree];
            magnitude := VarComplexPower(VarComplexAbsSqr(p), 0.5 / degree);
            angle     := VarComplexAngle(p) / degree;
            angle_inc := 8.0 * arctan(1.0) / degree;
            repeat
              roots[num_roots] := VarComplexFromPolar(magnitude, angle);
              inc(num_roots);
              angle := angle + angle_inc;
              dec(degree);
            until degree <> 0; // degree = 0 returns immediately.
        end // if

        else begin

                                      // Solve numerically.
            z := guess[attempt];       // Select a starting point.
            if (laguerre(degree,      // Number of coefficients - 1.
                         cc^[c],         // Ptr to array of complex coefficients.
                         z)) then     // Starting point and found root.
                e := 2                // 2 means Laguerre did not converge.
            else begin
              // Synthetic division (and check if z is really a root).
	            for n := degree - 1 downto 0 do // c[degree] stays the same.
                      cc^[c + n] := cc^[c + n] + cc^[c + n + 1] * z;

                if (VarComplexAbs(cc^[c]) < kEPSILON) then begin   // Must be (almost) zero.
                    dec(degree);
                    inc(c);                              // Remove first coefficient.
                    roots[num_roots] := z;               // Store the divisor.
                    inc(num_roots);
                    e := 0;
                end
                else
                    e := 3; // 3 means synthetic division left non-zero remainder.
            end; // else

            if (e > 0) then begin
                inc(attempt);
                if (attempt < n_attempts) then begin // Try all over again.
                    num_roots := num_roots - (degreeH - degree);    // Remove only the non-
                    degree := degreeH;                  // zero roots found.
                    c := 0;                             // Restore coeffs.
                    for n := 0 to degree do
                        cc^[n] := coeffs[coeff + n];
                end // if
                else
                  degree := 0; // Exit the loop, e=2 or e=3.
            end; // if
        end; // else
    end; // while

    if (degree = 2) then begin                // Solve quadratic equation analytically:
        a2 := cc^[c + 2] * 2.0;               //        --------
        D  := VarComplexSqrt(cc^[c + 1]*cc^[c+1]          //       / 2
                           -4.0 * cc^[c+2]*cc^[c+0]);           // D = \/ b - 4ac
        roots[num_roots] := ((-1 * cc^[c+1]) + D) / a2;          // x1 = (-b+D)/2a
        inc(num_roots);
        roots[num_roots] := ((-1 * cc^[c+1]) - D) / a2;          // x2 = (-b-D)/2a
        inc(num_roots);
    end // if

    else
      if (degree = 1) then begin            // Solve linear equation analytically:
                                            // a0 + a1 x = 0  <=>  x = -a0 / a1.
        roots[num_roots] := (-1 * cc^[c+0]) / cc^[c+1];    // Division by zero can never occur.
        inc(num_roots);

      end; // if

    freemem(cc, sizeof(variant) * (degree + 1));
    factor := e;                               // 0 = Ok.
end;

function test(degree: integer): integer;
var
  coeffs, roots: parrtype;
  i, e, n_roots: integer;
begin
  getmem(coeffs, sizeof(variant) * (degree + 1));

  coeffs^[2] := VarComplexCreate(2.0, 1.0);
  coeffs^[1] := -1 * VarComplexCreate(7.0, 1.0);
  coeffs^[0] := VarComplexCreate(15.0, 10.0);

  getmem(roots, sizeof(variant) * degree);
  if factor(degree + 1, coeffs^, roots^, n_roots) = 0 then begin

    for i := 0 to n_roots - 1 do begin
      Form1.Memo1.Lines.Add(format('root #%2d = (%s)', [i + 1, string(roots^[i])]));
    end;

  end;
  test := 0;
end;


procedure TForm1.Button1Click(Sender: TObject);
begin
  test(2);
end;

procedure TForm1.FormCreate(Sender: TObject);
begin
  _zero := VarComplexCreate(0.0, 0.0);
  _i := VarComplexCreate(0.0, 1.0);
end;

end.