unit dpolynomialf;

{$MODE FPC}
{$MODESWITCH OUT}
{$MODESWITCH RESULT}

interface

function ApplyPolynomialf(P: PSingle; N: PtrInt; X: Single): Single;
//
//  Returns
//
//      P[0] + P[1]*X + P[2]*X^2 + ... + P[N-1]*X^(N-1)
//

function SolveQuadricEquation(A, B, C: Single;
                              X: PSingle): PtrInt;
//
//  Solves the quadric equation for real X:
//
//      A * X^2 + B * X + C = 0
//
//  returns number of roots (0, 1 or 2)
//  and the roots itself in X[0], X[1] if found.
//

function SolveCubicEquation(A, B, C, D: Single;
                            X: PSingle): PtrInt;
//
//  Solves the cubic equation for real X:
//
//      A * X^3 + B * X^2 + C * X + D = 0
//
//  returns number of roots (0, 1, 2 or 3)
//  and the roots itself in X[0], X[1], X[2] if found.
//

function SolveQuarticEquation(A, B, C, D, E: Single;
                              X: PSingle): PtrInt;
//
//  Solves the quartic equation for real X:
//
//      A * X^4 + B * X^3 + C * X^2 + D * X + E = 0
//
//  returns number of roots (0, 1, 2, 3, 4)
//  and the roots itself in X[0], X[1], X[2], X[4] if found.
//

function NewtonPolynomial(P: PSingle;
                          N: PtrInt;
                          Neg, Pos: Single): Single;

procedure DiffPolynomial(P: PSingle;
                         N: PtrInt;
                         Q: PSingle);

function SolvePolynomialf(P: PSingle;
                          N: PtrInt;
                          A, B: Single;
                          X: PSingle): PtrInt;
//
//  Solves for X
//
//      P[0] + P[1]*X + P[2]*X^2 + ... + P[N]*X^N = 0
//
//  with condition A<=X<=B.
//
//  Returns number of roots (0 .. N)
//  and the roots itself in X[0] .. X[Result - 1],
//  so X should point to an array of at least N elements.
//

implementation

const
  EPS = 0.000001;

function ApplyPolynomialf(P: PSingle; N: PtrInt; X: Single): Single;
var
  A: Single;
  E: PSingle;
begin
  Result := 0;
  A := 1;
  E := P + N + 1;
  while P < E do begin
    Result := Result + P^ * A;
    Inc(P);
    A := A * X;
  end;
end;

function SolveQuadricEquation(A, B, C: Single;
                              X: PSingle): PtrInt;
var
  D: Single;
begin
  if Abs(A) < EPS then begin // linear case
    if Abs(B) < EPS then
      Exit(0);
    X[0] := - C / B;
    X[1] := X[0];
    Exit(1);
  end;

  D := B*B - 4*A*C;
  if D < 0 then begin
    Exit(0);
  end else if D < EPS then begin
    X[0] := -B / (2 * A);
    X[1] := X[0];
    Exit(1);
  end else begin
    D := Sqrt(D);
    X[0] := - (B - D) / (2 * A);
    X[1] := - (B + D) / (2 * A);
    Exit(2);
  end;
end;

// TODO i386 optimization
function ArcTan2(Y, X: Single): Single;
begin
    if (x=0) then
      begin
        if y=0 then
          Result:=0.0
        else if y>0 then
          Result:=pi/2
        else {if y<0 then}
          Result:=-pi/2;
      end
    else
      Result:=ArcTan(y/x);
    if x<0.0 then
      Result:=Result+pi;
    if Result>pi then
      Result:=Result-2*pi;
end;

function ArcCos(X: Single): Single;
begin
  if Abs(X) = 1.0 then begin
    if X < 0.0 then begin
      Result := Pi;
	end else
      Result := 0;
  end else
    Result := ArcTan2(Sqrt((1.0 - X) * (1.0 + X)), X);
end;

function SolveCubicEquation(A, B, C, D: Single;
                            X: PSingle): PtrInt;
label
  LAll3RootsAreReal;
var
  F, G, H: Single;
  I, J, K, M, N: Single;
  R, S, T, U: Single;
begin
  if Abs(A) < EPS then
    Exit(SolveQuadricEquation(B, C, D, X));

  //  Trivial normalization to the form X^3 + B * X^2 + C * X + D = 0
  B := B / A;
  C := C / A;
  D := D / A;

  //  Lets solve depressed equation t^3 + F*t + G = 0, where
  //
  //    F = (3c - b^2) / 3
  //    G = (2b^3 - 9bc + 27d) / 27
  //
  //  If t is a root of the depressed equation, then
  //
  //    X = t - b / 3
  //
  //  is a root of the original one.
  //
  F := ((3 * C) - Sqr(B)) / 3;
  G := ((2 * B * Sqr(B)) - (9 * B * C) + (27 * D)) / 27;


  //
  //  H is a discriminant of the equation. It determines number of roots, see
  //  ifs below for details.
  //
  H := (G * G / 4) + (F * F * F / 27);

  if H < 0 then begin // all 3 roots are real
LAll3RootsAreReal:
    //
    //  We gonna use the trigonometric solution.
    //
    I := Sqrt((Sqr(G) / 4) - H);
    J := Exp((1/3) * ln(I)); // = I^(1/3)
    K := ArcCos(- (G / (2 * I)));
    M := Cos(K / 3);
    N := Sqrt(3) * Sin(K / 3);
    X[0] := 2 * J * Cos(K / 3) - (B / 3);
    X[1] := - J * (M + N) - (B / 3);
    X[2] := - J * (M - N) - (B / 3);
    Exit(3);
  end else if H > 0 then begin // only 1 real root
    R := - (G / 2) + Sqrt(H);
    S := Exp((1/3) * ln(R)); // = R^(1/3)
    T := - (G / 2) - Sqrt(H);
    U := Exp((1/3) * ln(T)); // = T^(1/3)
    X[0] := (S + U) - (B / 3);
    X[1] := X[0];
    X[2] := X[0];
    Exit(1);
  end else if (Abs(F) < EPS) and (Abs(G) < EPS) and (Abs(H) < EPS) then begin
    // all 3 roots are real and equal
    X[0] := Exp((1/3) * ln(D)); // = D^(1/3)
    X[1] := X[0];
    X[2] := X[0];
    Exit(1);
  end else begin
    // H=0, but F<>0 or G<>0, the rarest case, but solution is the same as in
    // the H<0 case.
    goto LAll3RootsAreReal;
  end;
end;

function SolveQuarticEquation(A, B, C, D, E: Single;
                              X: PSingle): PtrInt;
var
  P, Q, R: Single;
  S, K: Single;
  D0, D1, D: Single;
begin
  if Abs(A) < EPS then
    Exit(SolveCubicEquation(B, C, D, E, X));
  Exit(0);

  //  Normilize
  B := B / A;
  C := C / A;
  D := D / A;
  E := E / A;

  //
  //  Convert into a depressed quartic form:
  //
  //    Y^4 + P * Y^2 + Q * Y + R = 0
  //
  //  If Y is a root of the depressed quadric equation, then
  //
  //    X = Y - B/4
  //
  //  is a root of the original one.
  //
  P := (8 * C - 3 * B * Sqr(B)) / 8;
  Q := (B * Sqr(B) - 4 * B * C + 8 * D) / 8;
  R := (-3 * Sqr(Sqr(B)) + 256*E + 64*B*D + 16*Sqr(B)*C) / 256;

  //
  //  Compute several helper values.
  //
  D0 := Sqr(C) - 3 * B * D + 12 * E;
  D1 := 2 * C * Sqr(C) - 9 * B * C * D + 27 * Sqr(B) * E + 27 * Sqr(D) - 72 * C * E;
  D := Sqr(D1) - 4 * D0 * Sqr(D0);

  if D < 0 then begin
  end;
  Q := Exp((1/3) * ln((D1 + Sqrt(D)) / 2));
  S := 0.5 * Sqrt((- 2 * P + (K + D0 / K)) / 3);

  //
  //
  //
end;

function NewtonPolynomial(P: PSingle;
                          N: PtrInt;
                          Neg, Pos: Single): Single;
var
  Mid, V: Single;
begin
  while Abs(Pos - Neg) > EPS do begin
    Mid := Neg + (Pos - Neg) / 2;
    V := ApplyPolynomialf(P, N, Mid);
    if V < 0 then begin
      Neg := Mid;
    end else
      Pos := Mid;
  end;
  Exit(Pos);
end;

procedure DiffPolynomial(P: PSingle;
                         N: PtrInt;
                         Q: PSingle);
begin
end;

function SolvePolynomialf(P: PSingle;
                          N: PtrInt;
                          A, B: Single;
                          X: PSingle): PtrInt;
var
  Diffs: PSingle;
  I: LongInt;
begin
//
// First we compute all derivatives of the polynomial down to second degree:
//    P P' P'' P''' ...
//
  Exit(0);
end;

end.
