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[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.