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}

const kEPSILON = 1e-6;

type
  TComplex = record
    re, im: Double;
  end;

function ComplexValue(ARe, AIm: Double): TComplex;
begin
  result.re := ARe; result.im := AIm;
end;

function ComplexMult(const a, b: TComplex): tcomplex; overload;
begin
  result.re := a.re * b.re - a.im * b.im;
  result.im := a.re * b.im + a.im * b.re;
end;
function ComplexMult(const a: TComplex; f: double): tcomplex; overload;
begin
  result.re := a.re * f;
  result.im := a.im * f;
end;

function ComplexDiv(const a, b: tcomplex): tcomplex; overload;
begin
  if abs(b.re) >= abs(b.im) then begin
    result.re := (a.re + a.im * b.im / b.re) / (b.re + sqr(b.im) / b.re);
    result.im := (a.im - a.re * b.im / b.re) / (b.re + sqr(b.im) / b.re);
  end
  else begin
    result.re := (a.im + a.re * b.re / b.im) / (b.im + sqr(b.re) / b.im);
    result.im := (a.im * b.re / b.im - a.re) / (b.im + sqr(b.re) / b.im);
  end;
end;
function ComplexDiv(r: double; zden: TComplex): TComplex; overload;
{ division : z := r / zden }
var denom : real;
begin
  denom := (zden.re * zden.re) + (zden.im * zden.im);
  { generates a fpu exception if denom=0 as for reals }
  result.re := (r * zden.re) / denom;
  result.im := - (r * zden.im) / denom;
end;


function ComplexAdd(const a, b: TComplex): tcomplex;
begin
  result.re := a.re + b.re;
  result.im := a.im + b.im;
end;
function ComplexSub(const a, b: tcomplex): tcomplex;
begin
  result.re := a.re - b.re;
  result.im := a.im - b.im;
end;

function ComplexMinus(const a: Tcomplex): tcomplex;
begin
  result := ComplexMult(a, -1);
end;

function ComplexIsZero(const a: TComplex): Boolean;
begin
  result := (abs(a.re) < kEpsilon) and (abs(a.im) < kEpsilon);
end;

function ComplexModule(const z: TComplex): Double;
{ module : r = |z| }
begin
  result := sqrt((z.re * z.re) + (z.im * z.im));
end;

function ComplexArg(const z: tcomplex): double;
begin
  result := arctan2(z.im, z.re);
end;


function ComplexSqrt(const z: TComplex): tcomplex;
{ square root : r := sqrt(z) }
var
  root, q: double;
begin
  if (z.re <> 0.0) or (z.im <> 0.0) then begin
    root := sqrt(0.5 * (abs(z.re) + ComplexModule(z)));
    q := z.im / (2.0 * root);
    if z.re >= 0.0 then begin
      result.re := root;
      result.im := q;
    end
    else
      if z.im < 0.0 then begin
        result.re := - q;
        result.im := - root
      end
      else begin
        result.re :=  q;
        result.im :=  root
      end
  end
  else result := z;
end;

function ComplexExp(const z: tcomplex): tcomplex;
var expz: double;
begin
  expz := exp(z.re);
  result.re := expz * cos(z.im);
  result.im := expz * sin(z.im);
end;
function ComplexLn(const z: tcomplex) : tcomplex;
begin
  result.re := ln(ComplexModule(z));
  result.im := arctan2(z.im, z.re);
end;


function ComplexPow(z1: tcomplex; r: double): tcomplex;
begin
  result := ComplexExp(ComplexMult(ComplexLn(z1), r));
end;



function ComplexPolar(magnitude, angle: double): TComplex;
begin
  result.re := magnitude * Cos(angle);
  result.im := magnitude * Sin(angle);
end;
function ComplexConj(z : tcomplex): tcomplex;
begin
  result.re := z.re;
  result.im := - z.im;
end;




function laguerre(degree: integer;             // Number of coefficients-1.
                  start: integer;
                  const a: array of TComplex;  // Array of coefficients.
                  var z: TComplex): boolean;   // Starting and end point.

const
  f: array[0 .. 7] of double = (
    1.0, 0.5, 0.25, 0.75, 0.13, 0.38, 0.62, 0.88
  );

var
  sn, snm: double;
  k: longint;
  j, kk: integer;
  ca, caca, cb, c1, dz, cc, cp, cm, zz: TComplex;
  alpha, beta, gamma: TComplex;

begin
  result := false;

  sn := degree; // degree = sn = number of coefficients - 1.
  snm := degree - 1;
  for k := 0 to pred(10000) do begin  // Not more than ten thousand iterations.
    alpha := a[start + degree];
    beta := ComplexValue(0.0, 0.0);
    gamma := ComplexValue(0.0, 0.0);

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

    if ComplexIsZero(alpha) then exit;

    ca   := ComplexDiv(ComplexMinus(beta), alpha);
    caca := ComplexMult(ca, ca);
    cb   := ComplexSub(caca, ComplexDiv(ComplexMult(gamma, 2.0), alpha));
    c1   := ComplexSqrt(ComplexMult(ComplexSub(ComplexMult(cb, sn), caca), snm));
    cp   := ComplexAdd(ca, c1);
    cm   := ComplexSub(ca, c1);

    if (ComplexModule(cm) > ComplexModule(cp)) then      // norm() is faster than abs().
        cp := cm;

    if ComplexIsZero(cp) then // Re-init trick from 'Numerical Recipes'.
        dz := ComplexMult(
        ComplexPolar(1.0, k), (1.0 + ComplexModule(z))
        )
    else
        dz := ComplexDiv(sn, cp);

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

    z := ComplexAdd(z, dz);
    if (kk <> 0) and (ComplexModule(dz) < 1e-31) then exit;
      // OK, converged: almost no change.
      // THIS 1e-31 IS RELATED TO kEPSILON !

    if k = 2000 then z := ComplexMinus(z)           // Kick z around when it
    else
      if k = 4000 then z := ComplexConj(z)          // is taking too long.
      else
        if k = 6000 then z := ComplexMinus(ComplexConj(z))
        else
          if k = 8000 then z := ComplexMult(z, ComplexValue(0.0, 1.0));
  end;
  result := true;
end;



type
  parrtype = ^arrtype;
  arrtype = array[0 .. pred(4096)] of tcomplex;

(*
    In most cases, the number of roots will equal the number of input coeffi-
    cients minus one. Only if the last (or latter) coefficient(s) is (are) zero,
    the number of resulting roots may be less than expected (coeffs.size() - 1).
    Coefficients expected in ascending order of degree.
*)
function factor(coeffs_size: integer;
         const coeffs: array of tcomplex;
         var roots: array of tcomplex;
         var num_roots: integer): integer;
const
  n_attempts = 4;
var
  p, z, a2, D: tcomplex;
  attempt, degree, degreeH, e, coeff, n: integer;
  guess: array[0 .. pred(n_attempts)] of tcomplex;
  magnitude, angle, angle_inc: double;

  cc: parrtype;
  c: integer;
begin
  attempt := 0;
  guess[0] := ComplexValue(8.0, 0.0);
  guess[1] := ComplexValue(-6.0, 0.0);
  guess[2] := ComplexValue(0.0, 6.0);
  guess[3] := ComplexValue(0.0, -8.0);

  (*
    complex<double> *cc = NULL, // Tmp copy of coefficients.
                    *c  = NULL; // Itarator within it.
    vector<complex<double> >::const_iterator  coeff  = coeffs.begin();
  *)
  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 ComplexIsZero(coeffs[degree]) do dec(degree);

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

  (*
    if (roots.size() != static_cast<unsigned int>(degree)) // Prepare output vector.
        roots.resize(degree);
  *)

  // Factor-out roots at the origin.
  while (degree > 0) and ComplexIsZero(coeffs[coeff + 0]) do begin
    roots[num_roots] := ComplexValue(0.0, 0.0); // Insert root at 0.
    inc(num_roots);
    inc(coeff);
    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).
    getmem(cc, sizeof(tcomplex) * (degree + 1));
    c := 0;
    // cc = c = new complex<double>[degree + 1]; // MAY THROW AN EXCEPTION.
    for n := 0 to degree do cc^[c + n] := coeffs[coeff + n];

  end;

  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 ComplexIsZero(cc^[c + n]) do inc(n);
                                        // 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 := ComplexDiv(ComplexMinus(cc^[c + 0]), cc^[c + degree]);
            magnitude := power(ComplexModule(p), 0.5 / degree);
            angle     := ComplexArg(p) / degree;
            angle_inc := 8.0 * arctan(1.0) / degree;
            repeat
                roots[num_roots] := ComplexPolar(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.
                        c,      // Offset
                        cc^,    // Ptr to array of complex coefficients.
                        z       // Starting point and found root.
               )
            then 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] := ComplexAdd(cc^[c + n], ComplexMult(cc^[c + n + 1], z));
                if (ComplexModule(cc^[c + 0]) < 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^[c + n] := coeffs[coeff + n];
                end
                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 := ComplexMult(cc^[c + 2], 2.0);   //        --------
        D  := ComplexSqrt(
                ComplexSub(
                  ComplexMult(cc^[c + 1], cc^[c + 1]),
                  ComplexMult(
                    ComplexMult(cc^[c + 2], cc^[c + 0]),
                    4.0
                  )
                )
              );
        roots[num_roots] := ComplexDiv(ComplexAdd(ComplexMinus(cc^[c + 1]), D), a2);          // x1 = (-b+D)/2a
        inc(num_roots);
        roots[num_roots] := ComplexDiv(ComplexSub(ComplexMinus(cc^[c + 1]), D), a2);          // x2 = (-b-D)/2a
        inc(num_roots);
    end
    else
      if (degree = 1) then begin            // Solve linear equation analytically:
                                            // a0 + a1 x = 0  <=>  x = -a0 / a1.
        roots[num_roots] := ComplexDiv(ComplexMinus(cc^[c + 0]), cc^[c + 1]);    // Division by zero can never occur.
        inc(num_roots);

      end;
      // else (if (degree == 0)) we're done.

      freemem(cc, sizeof(tcomplex) * (degree + 1)); // Release temporary workspace.
    result := e;                               // 0 = Ok.
end;

procedure test(degree: integer);
const
  max_degree = 64;
var
  coeffs: array[0 .. max_degree] of tcomplex;
  roots: array[0 .. pred(max_degree)] of tcomplex;

  i, e, n_roots: integer;
begin
    //vector<complex<double> >    coeffs(degree + 1);         // INPUT.

    coeffs[3] := ComplexValue(1.0, 0.0);
    coeffs[2] := ComplexMinus(ComplexValue(6.0, 0.0));
    coeffs[1] := ComplexValue(21.0, 0.0);
    coeffs[0] := ComplexMinus(ComplexValue(52.0, 0.0));
    (*
    coeffs[2] := ComplexValue(2.0, 1.0);
    coeffs[1] := ComplexMinus(ComplexValue(7.0, 1.0));
    coeffs[0] := ComplexValue(15.0, 10.0);
    *)

    // vector<complex<double> >    roots(degree);              // OUT/INPUT.
    // vector<complex<double> >    coeffs_back(degree + 1);    // OUTPUT.

    e := factor(degree + 1,
                coeffs,           // Input.
                roots, n_roots);  // Output.

    for i := 0 to pred(n_roots) do begin
      Form1.Memo1.Lines.Add(format('root #%2d = (%6.2f, %6.2f)', [i + 1, roots[i].re, roots[i].im]));
      // writeln('root #', i + 1, ' = (', roots[i].re:6:2, ' , ', roots[i].im: 6:2, ')')
    end;
end;


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

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

end.