{$n+}

Unit Equation;

Interface

Const
  maxOrder = 12;      { maximum order for a polynomial }


Function solveQuadratic(Var x, y: Array Of Double): Byte;
Function solveCubic(Var x, y: Array Of Double): Byte;

Function solveQuarticAlgebra(Var x, results: Array Of Double): Byte;
Function solveQuarticVieta(Var x, results: Array Of Double): Byte;


Function PolySolve(Const order: Integer;
         Var coeffs, roots: Array Of Double): Integer;


Implementation

Const
  Epsilon     = 1.0E-10;
  maxDistance = 1.0E+07;


{
  Solve the quadratic equation:

  x[0] * x^2 + x[1] * x + x[2] = 0.

  The value returned by this function is the number of real roots.
  The roots themselves are returned in y[0], y[1].
}
Function solveQuadratic(Var x, y: Array Of Double): Byte;
  Var
    a, b, c, Desc: Double;
    T: Double;
  Begin
    a := x[0]; b := -x[1]; c := x[2];
    If a = 0 Then
      Begin
        solveQuadratic := 0;
        If b = 0.0 Then Exit;

        solveQuadratic := 1;
        y[0] := c / b; Exit
      End;

    solveQuadratic := 0;
    Desc := Sqr(b) - 4 * a * c;
    If Desc < 0.0 Then Exit
    Else
      If Abs(Desc) < Epsilon Then
        Begin
          solveQuadratic := 1;
          y[0] := 0.5 * b / a;
          Exit
        End;

    Desc := Sqrt(Desc);
    T := 2 * a;
    y[0] := (b + Desc) / T;
    y[1] := (b - Desc) / T;
    solveQuadratic := 2
  End;


{
  Solve the cubic equation:

  x[0] * x^3 + x[1] * x^2 + x[2] * x + x[3] = 0.

  The result of this function is an integer that tells how many real
  roots exist.  Determination of how many are distinct is up to the
  process that calls this routine.  The roots that exist are stored
  in (y[0], y[1], y[2]).

  Note: this function relies very heavily on trigonometric functions and
  the square root function. If an alternative solution is found that does
  not rely on transcendentals this code should be replaced.
}
Function solveCubic(Var x, y: Array Of Double): Byte;

  Function ArcCos(x: Double): Double;
    Begin
      If x = 0 Then
        ArcCos := Pi / 2
      Else
        ArcCos := ArcTan( Sqrt(1 - Sqr(x)) / x ) + Pi * Byte(x < 0)
    End;

  Function Pow(a, x: Double): Double;
    Begin
      Pow := Exp( x * Ln(a) )
    End;

  Var
    a0, a1, a2, a3, an: Double;
    sq, q, q3, r, r2, sqrA: Double;
    theta, Desc: Double;
  Begin
    a0 := x[0];
    If a0 = 0.0 Then
      Begin
        solveCubic := solveQuadratic(x[1], y);
        Exit
      End
    Else
      Begin
        a1 := x[1] / a0;
        a2 := x[2] / a0;
        a3 := x[3] / a0;
      End;

    sqrA := sqr(a1);
    q := (sqrA - 3.0 * a2) / 9.0;
    r := (2.0 * sqrA * a1 - 9.0 * a1 * a2 + 27.0 * a3) / 54.0;
    q3 := sqr(q) * q;
    r2 := sqr(r);
    Desc := q3 - r2;
    an := a1 / 3.0;
    If Desc >= 0 Then
      Begin
        { Three real roots }
        solveCubic := 3;

        Desc := r / sqrt(q3);
        theta := ArcCos(Desc) / 3;
        sq := -2.0 * Sqrt(q);
        y[0] := sq * Cos(theta) - an;
        y[1] := sq * Cos(theta + (2 * Pi / 3)) - an;
        y[2] := sq * Cos(theta + (4 * Pi / 3)) - an;
        Exit
      End
    Else
      Begin
        solveCubic := 1;

        sq := Pow(Sqrt(r2 - q3) + Abs(r), 1 / 3);
        If r < 0 Then
          y[0] := (sq + q / sq) - an
        Else
          y[0] := -(sq + q / sq) - an;
        Exit
      End
  End;




(*
  Additional functions for solveQuartic.
*)
Const
  fudgeFactorFirst  =  1.0E+12;
  fudgeFactorSecond = -1.0E-05;
  fudgeFactorThird  =  1.0E-07;


Function DifficultCoeffs(n: Integer; Var x: Array Of Double): Boolean;
  Var
    i: Integer;
    biggest: Double;
  Begin
    biggest := Abs(x[0]);
    For i := 1 To n Do
      If Abs(x[i]) > biggest Then
        biggest := Abs(x[i]);

    DifficultCoeffs := False;
    If biggest = 0.0 Then Exit;

    For i := 0 To n Do
      If x[i] <> 0.0 Then
        If Abs(biggest / x[i]) > fudgeFactorFirst Then
          Begin
            DifficultCoeffs := True;
            x[i] := 0.0
          End
  End;


{
  root solver based on the STURM sequences for a polynomial
}
Const
  RelError = 1.0e-14; { smallest relative error we want }
  maxPow   = 32;      { max power of 10 we wish to search to }
  maxIt    = 800;     { max number of iterations }

  {
    a coefficient smaller than SMALL_ENOUGH is
    considered to be zero (0.0).
  }
  SmallEnough = 1.0e-12;

{
  structure type for representing a polynomial
}
Type
  _poly = Record
    ord: Integer;
    coef: Array[0 .. maxOrder] Of Double;
  End;

(*
 * evalpoly
 *
 * evaluate polynomial defined in coef returning its value.
 *)
Function EvalPoly(Const Ord: Integer;
         Var coef: Array Of Double; x: Double): Double;
Var
  f: Double;
  fp: Integer;
Begin
  f := coef[Ord];
  For fp := Pred(Ord) DownTo 0 Do
    f := x * f + coef[fp];
  EvalPoly := f
End;


(*
 * modrf
 *
 * uses the modified regula-falsi method to evaluate the root
 * in interval [a,b] of the polynomial described in coef. The
 * root is returned is returned in *val. The routine returns zero
 * if it can't converge.
 *)
Function modrf(Const Ord: Integer; Var coef: Array Of Double;
         a, b: Double; Var val: Double): Boolean;
Var
  scoef, ecoef, fp, its: Integer;
  fa, fb, x, fx, lfx: Double;
Begin
  scoef := 0;
  ecoef := ord;

  fa := coef[ecoef]; fb := fa;
  For fp := Pred(ecoef) DownTo scoef Do Begin
    fa := a * fa + coef[fp];
    fb := b * fb + coef[fp];
  End;

  (*
   * if there is no sign difference the method won't work
   *)
  If fa * fb > 0.0 Then Begin
    modrf := False; Exit
  End;

  If Abs(fa) < RelError Then Begin
    val := a;
    modrf := True; Exit
  End;

  If Abs(fb) < RelError Then Begin
    val := b;
    modrf := True; Exit
  End;

  lfx := fa;
  For its := 0 To Pred(maxIt) Do Begin
    x := (fb * a - fa * b) / (fb - fa);

    fx := coef[ecoef];
    For fp := Pred(ecoef) DownTo scoef Do
      fx := x * fx + coef[fp];

    If Abs(x) > RelError Then Begin
      If Abs(fx / x) < RelError Then Begin
        val := x;
        modrf := True; Exit
      End
    End
    Else
      If Abs(fx) < RelError Then Begin
        val := x;
        modrf := True; Exit
      End;

    If (fa * fx) < 0 Then Begin
      b := x;
      fb := fx;
      If (lfx * fx) > 0 Then
        fa := fa / 2
    End
    Else Begin
      a := x;
      fa := fx;
      if (lfx * fx) > 0 Then
        fb := fb / 2
    End;

    lfx := fx
  End;

  {$IFDEF DEBUG_VER}
    WriteLn('modrf overflow ', a, ' ', b, ' ', fx);
  {$ENDIF}

  modrf := False
End;


(*
 * modp
 *
 * calculates the modulus of u(x) / v(x) leaving it in r, it
 * returns 0 if r(x) is a constant.
 * note: this function assumes the leading coefficient of v
 * is 1 or -1
 *)
Function modp(Var u, v, r: _poly): Integer;
Var
  k, j: Integer;
  nr, _end, uc: Integer;
Begin

  nr := 0;
  _end := u.ord;

  uc := 0;
  While uc <= _end Do Begin
    r.coef[nr] := u.coef[uc];
    Inc(nr); Inc(uc)
  End;

  If v.coef[v.ord] < 0.0 Then Begin
    k := Pred(u.ord - v.ord);
    While k >= 0 Do Begin
      r.coef[k] := -r.coef[k];
      Dec(k, 2)
    End;

    For k := u.ord - v.ord DownTo 0 Do
      For j := Pred(v.ord + k) DownTo k Do
        r.coef[j] := -r.coef[j] - r.coef[v.ord + k] * v.coef[j - k]
  End
  Else Begin
    For k := u.ord - v.ord DownTo 0 Do
      For j := Pred(v.ord + k) DownTo k Do
        r.coef[j] := r.coef[j] - r.coef[v.ord + k] * v.coef[j - k]
  End;

  k := Pred(v.ord);
  While (k >= 0) and (Abs(r.coef[k]) < SmallEnough) Do Begin
    r.coef[k] := 0.0;
    Dec(k)
  End;

  r.ord := k * Byte(k >= 0);
  modp := r.ord
End;

(*
 * buildsturm
 *
 * build up a sturm sequence for a polynomial in smat, returning
 * the number of polynomials in the sequence
 *)
Function buildsturm(Const Ord: Integer;
         Var sseq: Array Of _poly): Integer;
Var
  f: Double;
  i, fp, fc, sp: Integer;
Begin

  sseq[0].ord := Ord;
  sseq[1].ord := Pred(Ord);

  (*
   * calculate the derivative and normalise the leading
   * coefficient.
   *)
  f := Abs(sseq[0].coef[Ord] * Ord);

  fp := 0;
  fc := 1;
  For i := 1 To Ord Do Begin
    sseq[1].coef[fp] := sseq[0].coef[fc] * i / f;
    Inc(fp); Inc(fc)
  End;

  (*
   * construct the rest of the Sturm sequence
   *)
  sp := 2;
  While modp(sseq[sp-2], sseq[sp-1], sseq[sp]) <> 0 Do Begin
    (*
     * reverse the sign and normalise
     *)
    f := -Abs(sseq[sp].coef[sseq[sp].ord]);
    For fp := sseq[sp].ord DownTo 0 Do
      sseq[sp].coef[fp] := sseq[sp].coef[fp] / f;

    Inc(sp)
  End;

  sseq[sp].coef[0] := -sseq[sp].coef[0]; { reverse the sign }
  buildsturm := sp
End;

(*
 * numroots
 *
 *	return the number of distinct real roots of the polynomial
 * described in sseq.
 *)
Function numroots(np: Integer; Var sseq: Array Of _poly;
         Var atneg, atpos: Integer): Integer;
Var
  s, atposinf, atneginf: Integer;
  f, lf: Double;
Begin

  atposinf := 0; atneginf := 0;
  (*
   * changes at positive infinity
   *)
  lf := sseq[0].coef[sseq[0].ord];

  For s := 1 To np Do Begin

    f := sseq[s].coef[sseq[s].ord];
    If (lf = 0.0) or (lf * f < 0) Then Inc(atposinf);
    lf := f
  End;

  (*
   * changes at negative infinity
   *)
  If (sseq[0].ord and $1) > 0 Then
    lf := -sseq[0].coef[sseq[0].ord]
  Else lf := sseq[0].coef[sseq[0].ord];

  For s := 1 To np Do Begin

    If (sseq[s].ord and $1) > 0 Then
      f := -sseq[s].coef[sseq[s].ord]
    Else f := sseq[s].coef[sseq[s].ord];

    If (lf = 0.0) or (lf * f < 0) Then Inc(atneginf);
    lf := f
  End;

  atneg := atneginf;
  atpos := atposinf;

  numroots := atneginf - atposinf
End;

(*
 * numchanges
 *
 * return the number of sign changes in the Sturm sequence in
 * sseq at the value a.
 *)
Function numchanges(np: Integer; Var sseq: Array Of _poly;
         a: Double): Integer;
Var
  s, changes: Integer;
  f, lf: Double;
Begin

  changes := 0;
  lf := EvalPoly(sseq[0].ord, sseq[0].coef, a);

  For s := 1 To np Do Begin
    f := EvalPoly(sseq[s].ord, sseq[s].coef, a);
    If (lf = 0.0) or (lf * f < 0) Then Inc(changes);
    lf := f
  End;

  numchanges := changes
End;

(*
 * sbisect
 *
 * uses a bisection based on the sturm sequence for the polynomial
 * described in sseq to isolate intervals in which roots occur,
 * the roots are returned in the roots array in order of magnitude.
 *)
Procedure sbisect(np: Integer; Var sseq: Array Of _poly;
          min, max: Double; atmin, atmax: Integer;
          Var roots: Array Of Double);
Var
  n1, n2, nroot, its, atmid: Integer;
  mid: Double;
Begin

  n1 := 0; n2 := 0;
  nroot := atmin - atmax;

  If nroot = 1 Then Begin
    (*
     * first try a less expensive technique.
     *)
    If modrf(sseq[0].ord, sseq[0].coef, min, max, roots[0]) Then Exit;

    (*
     * if we get here we have to evaluate the root the hard
     * way by using the Sturm sequence.
     *)
    For its := 0 To Pred(maxIt) Do Begin

      mid := (min + max) / 2;
      atmid := numchanges(np, sseq, mid);

      If Abs(mid) > RelError Then Begin
	If Abs((max - min) / mid) < RelError Then Begin
	  roots[0] := mid; Exit
        End;
      End
      Else
        If Abs(max - min) < RelError Then Begin
	  roots[0] := mid; Exit
        End;

      If (atmin - atmid) = 0 Then min := mid
      Else max := mid

    End;

    If its = maxIt Then Begin
      {$IFDEF DEBUG_VER}
        WriteLn('sbisect: overflow min ', min, ' max ', max,
                ' diff ', max-min, ' nroot ', nroot, ' n1 ', n1, ' n2 ', n2);
      {$ENDIF}
      roots[0] := mid
    End;
    Exit
  End;

  (*
   * more than one root in the interval, we have to bisect...
   *)
  For its := 0 To Pred(maxIt) Do Begin

    mid := (min + max) / 2;
    atmid := numchanges(np, sseq, mid);

    n1 := atmin - atmid;
    n2 := atmid - atmax;

    If (n1 <> 0) and (n2 <> 0) Then Begin
      sbisect(np, sseq, min, mid, atmin, atmid, roots);
      sbisect(np, sseq, mid, max, atmid, atmax, roots[n1]);
      Break
    End;

    If n1 = 0 Then min := mid
    Else max := mid

  End;

  If its = maxIt Then Begin
    {$IFDEF DEBUG_VER}
      WriteLn('sbisect: roots too close together');
      WriteLn('sbisect: overflow min ', min, ' max ', max,
              ' diff ', max-min, ' nroot ', nroot, ' n1 ', n1, ' n2 ', n2);
    {$ENDIF}
    For n1 := atmax To Pred(atmin) Do
      roots[n1 - atmax] := mid
  End
End;


Var
  sseq: array[0 .. maxOrder] of _poly;

Function PolySolve(Const order: Integer;
         Var coeffs, roots: Array Of Double): Integer;
Var
  i, j, np, nroots, nchanges, atmin, atmax: Integer;
  min, max: double;
Begin
  For i := order DownTo 0 Do Begin
    sseq[0].coef[i] := coeffs[order-i]
  End;

  np := buildsturm(order, sseq);

  {$IFDEF DEBUG_VER}
    WriteLn('Sturm sequence for:');
    For i := order DownTo 0 Do
      Write(sseq[0].coef[i]:7:4, ' ');
    WriteLn;

    For i := 0 To np Do Begin
      For j := sseq[i].ord DownTo 0 Do
        Write(sseq[i].coef[j]:7:4, ' ');
      WriteLn
    End;

    WriteLn;
  {$ENDIF}

  (*
   * get the number of real roots
   *)
  nroots := numroots(np, sseq, atmin, atmax);
  If nroots = 0 Then Begin
    PolySolve := 0; Exit
  End;

  (*
   * calculate the bracket that the roots live in
   *)
  min := -1.0;
  nchanges := numchanges(np, sseq, min);
  i := 0;
  While (nchanges <> atmin) and (i <> maxPow) Do Begin
    min := min * 10.0;
    nchanges := numchanges(np, sseq, min);
    Inc(i)
  End;

  If nchanges <> atmin Then Begin

    {$IFDEF DEBUG_VER}
      WriteLn('solve: unable to bracket all negative roots');
    {$ENDIF}

    atmin := nchanges
  End;

  max := 1.0;
  nchanges := numchanges(np, sseq, max);
  i := 0;
  While (nchanges <> atmax) and (i <> maxPow) Do Begin
    max := max * 10.0;
    nchanges := numchanges(np, sseq, max);
    Inc(i)
  End;

  If nchanges <> atmax Then Begin

    {$IFDEF DEBUG_VER}
      WriteLn('solve: unable to bracket all positive roots');
    {$ENDIF}

    atmax := nchanges
  End;

  nroots := atmin - atmax;

  (*
   * perform the bisection.
   *)
  sbisect(np, sseq, min, max, atmin, atmax, roots);
  PolySolve := nroots
End;



{
  solving quartics algebraically (method of Lodovico Ferrari)
}
Function solveQuarticAlgebra(Var x, results: Array Of Double): Byte;
  Var
    i: Integer;
    a0, a1, y, d1, x1, t1, t2: Double;
    c0, c1, c2, c3, c4, d2, q1, q2: Double;
    cubic: Array[0 .. 3] Of Double;
    roots: Array[0 .. 2] Of Double;
  Begin
    If DifficultCoeffs(4, x) Then
      Begin
        If Abs(x[0]) < Epsilon Then
          If Abs(x[1]) < Epsilon Then
            Begin
              solveQuarticAlgebra := solveQuadratic(x[2], results); Exit
            End
          Else
            Begin
              solveQuarticAlgebra := solveCubic(x[1], results); Exit
            End

        Else
          Begin
            solveQuarticAlgebra := PolySolve(4, x, results); Exit
          End
      End;

    c0 := x[0];
    If Abs(c0) < Epsilon Then
      Begin
        solveQuarticAlgebra := solveCubic(x[1], results); Exit
      End
    Else
      If Abs(x[4]) < Epsilon Then
        Begin
          solveQuarticAlgebra := solveCubic(x, results); Exit
        End
      Else
        Begin
          c1 := x[1] / c0;
          c2 := x[2] / c0;
          c3 := x[3] / c0;
          c4 := x[4] / c0;
        End;

      (*
        The first step is to take the original equation:

	   x^4 + b*x^3 + c*x^2 + d*x + e = 0

        and rewrite it as:

	   x^4 + b*x^3 = -c*x^2 - d*x - e,

        adding (b*x/2)^2 + (x^2 + b*x/2)y + y^2/4 to each side gives a
        perfect square on the lhs:

	   (x^2 + b*x/2 + y/2)^2 = (b^2/4 - c + y)x^2 + (b*y/2 - d)x + y^2/4 - e

        By choosing the appropriate value for y, the rhs can be made a perfect
        square also.  This value is found when the rhs is treated as a quadratic
        in x with the discriminant equal to 0.  This will be true when:

	   (b*y/2 - d)^2 - 4.0 * (b^2/4 - c*y)*(y^2/4 - e) = 0, or

	   y^3 - c*y^2 + (b*d - 4*e)*y - b^2*e + 4*c*e - d^2 = 0.

        This is called the resolvent of the quartic equation
      *)

      a0 := 4.0 * c4;
      cubic[0] := 1.0;
      cubic[1] := -1.0 * c2;
      cubic[2] := c1 * c3 - a0;
      cubic[3] := a0 * c2 - Sqr(c1) * c4 - Sqr(c3);
      i := solveCubic(cubic[0], roots[0]);

      solveQuarticAlgebra := 0;
      If i > 0 Then
        y := roots[0]
      Else
        Exit;

      (*
        What we are left with is a quadratic squared on the lhs and a
        linear term on the right.  The linear term has one of two signs,
        take each and add it to the lhs.  The form of the quartic is now:

	   a' = b^2/4 - c + y,    b' = b*y/2 - d, (from rhs quadritic above)

	   (x^2 + b*x/2 + y/2) = +sqrt(a'*(x + 1/2 * b'/a')^2), and
	   (x^2 + b*x/2 + y/2) = -sqrt(a'*(x + 1/2 * b'/a')^2).

        By taking the linear term from each of the right hand sides and
        adding to the appropriate part of the left hand side, two quadratic
        formulas are created.  By solving each of these the four roots of
        the quartic are determined
      *)

      i := 0;
      a0 := c1 / 2.0;
      a1 := y / 2.0;

      t1 := Sqr(a0) - c2 + y;
      If t1 < 0.0 Then
        Begin
          If t1 > fudgeFactorSecond
            Then t1 := 0.0
          Else
            { First Special case, a' < 0 means all roots are complex }
            Exit
        End;

      If t1 < fudgeFactorThird Then
        Begin
          (*
            Second special case, the "x" term on the right hand side above
	    has vanished.  In this case:
		      (x^2 + b*x/2 + y/2) = +sqrt(y^2/4 - e), and
		      (x^2 + b*x/2 + y/2) = -sqrt(y^2/4 - e).
          *)
          t2 := Sqr(a1) - c4;
          If t2 < 0.0 Then Exit;

          x1 := 0.0;
          d1 := Sqrt(t2)
        End
      Else
        Begin
          x1 := Sqrt(t1);
          d1 := 0.5 * (a0 * y - c3) / x1
        End;

      { Solve the first quadratic }
      q1 := -a0 - x1;
      q2 := a1 + d1;
      d2 := Sqr(q1) - 4.0 * q2;

      If d2 >= 0.0 Then
        Begin
          d2 := Sqrt(d2);
          results[0] := 0.5 * (q1 + d2);
          results[1] := 0.5 * (q1 - d2);
          i := 2
        End;

      { Solve the second quadratic }
      q1 := q1 + x1 + x1;
      q2 := a1 - d1;
      d2 := Sqr(q1) - 4.0 * q2;
      If d2 >= 0.0 Then
        Begin
          d2 := Sqrt(d2);
          results[i] := 0.5 * (q1 + d2); Inc(i);
          results[i] := 0.5 * (q1 - d2); Inc(i)
        End;

    solveQuarticAlgebra := i
  End;


{
  Solve a quartic using the method of Francois Vieta (Circa 1735)
}
Function solveQuarticVieta(Var x, results: Array Of Double): Byte;
  Var
    i: Integer;
    c0, c1, c2, c3, c4: Double;
    c12, z, p, q, q1, q2, r, d1, d2: Double;

    cubic: Array[0 .. 3] Of Double;
    roots: Array[0 .. 2] Of Double;
  Begin
    { Figure out the size difference between coefficients }
    If DifficultCoeffs(4, x) Then

      If Abs(x[0]) < Epsilon Then
        If Abs(x[1]) < Epsilon Then
          Begin
            solveQuarticVieta := solveQuadratic(x[2], results);
            Exit
          End
        Else
          Begin
            solveQuarticVieta := solveCubic(x[1], results);
            Exit;
          End
      Else
        Begin
          solveQuarticVieta := PolySolve(4, x, results);
          Exit
        End;

    { See if the high order term has vanished }
    c0 := x[0];
    If Abs(c0) < Epsilon Then
      Begin
        solveQuarticVieta := solveCubic(x[1], results); Exit
      End;

    { See if the constant term has vanished }
    If Abs(x[4]) < Epsilon Then
      Begin
        solveQuarticVieta := solveCubic(x, results); Exit
      End;

    { Make sure the quartic has a leading coefficient of 1.0 }
    c1 := x[1] / c0;
    c2 := x[2] / c0;
    c3 := x[3] / c0;
    c4 := x[4] / c0;

    { Compute the cubic resolvant }
    c12 := c1 * c1;
    p := -0.375 * c12 + c2;
    q := 0.125 * c12 * c1 - 0.5 * c1 * c2 + c3;
    r := -0.01171875 * c12 * c12 + 0.0625 * c12 * c2 - 0.25 * c1 * c3 + c4;

    cubic[0] := 1.0;
    cubic[1] := -0.5 * p;
    cubic[2] := -r;
    cubic[3] := 0.5 * r * p - 0.125 * Sqr(q);
    i := solveCubic(cubic, roots);

    solveQuarticVieta := 0;
    If i > 0 Then
      z := roots[0]
    Else Exit;

    d1 := 2.0 * z - p;

    If d1 < 0.0 Then
      If d1 > - Epsilon Then
        d1 := 0.0
      Else Exit;

    If d1 < Epsilon Then
      Begin
        d2 := Sqr(z) - r;
        If d2 < 0.0 Then Exit;
        d2 := Sqrt(d2)
      End
    Else
      Begin
        d1 := Sqrt(d1);
        d2 := 0.5 * q / d1
      End;

    { Set up useful values for the quadratic factors }
    q1 := Sqr(d1);
    q2 := -0.25 * c1;
    i := 0;

    { Solve the first quadratic }
    p := q1 - 4.0 * (z - d2);
    If p = 0 Then
      Begin
        results[i] := -0.5 * d1 - q2; Inc(i)
      End
    Else
      If p > 0 Then
        Begin
          p := Sqrt(p);
          results[i] := -0.5 * (d1 + p) + q2; Inc(i);
          results[i] := -0.5 * (d1 - p) + q2; Inc(i);
        End;

    { Solve the second quadratic }
    p := q1 - 4.0 * (z + d2);
    If p = 0 Then
      Begin
        results[i] := 0.5 * d1 - q2; Inc(i)
      End
    Else
      If p > 0 Then
        Begin
          p := Sqrt(p);
          results[i] := 0.5 * (d1 + p) + q2; Inc(i);
          results[i] := 0.5 * (d1 - p) + q2; Inc(i)
        End;

    solveQuarticVieta := i
  End;

End.