// Copyright (c) 1997-1999 Matra Datavision
// Copyright (c) 1999-2014 OPEN CASCADE SAS
//
// This file is part of Open CASCADE Technology software library.
//
// This library is free software; you can redistribute it and/or modify it under
// the terms of the GNU Lesser General Public License version 2.1 as published
// by the Free Software Foundation, with special exception defined in the file
// OCCT_LGPL_EXCEPTION.txt. Consult the file LICENSE_LGPL_21.txt included in OCCT
// distribution for complete text of the license and disclaimer of any warranty.
//
// Alternatively, this file may be used under the terms of Open CASCADE
// commercial license or contractual agreement.

// #ifndef OCCT_DEBUG
#define No_Standard_RangeError
#define No_Standard_OutOfRange
#define No_Standard_DimensionError

// #endif

#include <math_DirectPolynomialRoots.hxx>
#include <math_FunctionRoots.hxx>
#include <math_FunctionWithDerivative.hxx>
#include <math_BracketedRoot.hxx>
#include <TColStd_Array1OfReal.hxx>

#define ITMAX 100
#define EPS 1e-14
#define EPSEPS 2e-14
#define MAXBIS 100

#ifdef OCCT_DEBUG
static Standard_Boolean myDebug = 0;
static Standard_Integer nbsolve = 0;
#endif

class DerivFunction : public math_Function
{
  math_FunctionWithDerivative* myF;

public:
  DerivFunction(math_FunctionWithDerivative& theF)
      : myF(&theF)
  {
  }

  virtual Standard_Boolean Value(const Standard_Real theX, Standard_Real& theFval)
  {
    return myF->Derivative(theX, theFval);
  }
};

static void AppendRoot(TColStd_SequenceOfReal&      Sol,
                       TColStd_SequenceOfInteger&   NbStateSol,
                       const Standard_Real          X,
                       math_FunctionWithDerivative& F,
                       //			const Standard_Real K,
                       const Standard_Real,
                       const Standard_Real dX)
{

  Standard_Integer n = Sol.Length();
  Standard_Real    t;
#ifdef OCCT_DEBUG
  if (myDebug)
  {
    std::cout << "   Ajout de la solution numero : " << n + 1 << std::endl;
    std::cout << "   Valeur de la racine :" << X << std::endl;
  }
#endif

  if (n == 0)
  {
    Sol.Append(X);
    F.Value(X, t);
    NbStateSol.Append(F.GetStateNumber());
  }
  else
  {
    Standard_Integer i  = 1;
    Standard_Integer pl = n + 1;
    while (i <= n)
    {
      t = Sol.Value(i);
      if (t >= X)
      {
        pl = i;
        i  = n;
      }
      if (Abs(X - t) <= dX)
      {
        pl = 0;
        i  = n;
      }
      i++;
    } //-- while
    if (pl > n)
    {
      Sol.Append(X);
      F.Value(X, t);
      NbStateSol.Append(F.GetStateNumber());
    }
    else if (pl > 0)
    {
      Sol.InsertBefore(pl, X);
      F.Value(X, t);
      NbStateSol.InsertBefore(pl, F.GetStateNumber());
    }
  }
}

static void Solve(math_FunctionWithDerivative& F,
                  const Standard_Real          K,
                  const Standard_Real          x1,
                  const Standard_Real          y1,
                  const Standard_Real          x2,
                  const Standard_Real          y2,
                  const Standard_Real          tol,
                  const Standard_Real          dX,
                  TColStd_SequenceOfReal&      Sol,
                  TColStd_SequenceOfInteger&   NbStateSol)
{
#ifdef OCCT_DEBUG
  if (myDebug)
  {
    std::cout << "--> Resolution :" << ++nbsolve << std::endl;
    std::cout << "   x1 =" << x1 << " y1 =" << y1 << std::endl;
    std::cout << "   x2 =" << x2 << " y2 =" << y2 << std::endl;
  }
#endif

  Standard_Integer iter  = 0;
  Standard_Real    tols2 = 0.5 * tol;
  Standard_Real    a, b, c, d = 0, e = 0, fa, fb, fc, p, q, r, s, tol1, xm, min1, min2;
  a = x1;
  b = c = x2;
  fa    = y1;
  fb = fc = y2;
  for (iter = 1; iter <= ITMAX; iter++)
  {
    if ((fb > 0.0 && fc > 0.0) || (fb < 0.0 && fc < 0.0))
    {
      c  = a;
      fc = fa;
      e = d = b - a;
    }
    if (Abs(fc) < Abs(fb))
    {
      a  = b;
      b  = c;
      c  = a;
      fa = fb;
      fb = fc;
      fc = fa;
    }
    tol1 = EPSEPS * Abs(b) + tols2;
    xm   = 0.5 * (c - b);
    if (Abs(xm) < tol1 || fb == 0)
    {
      //-- On tente une iteration de newton
      Standard_Real    Xp, Yp, Dp;
      Standard_Integer itern = 5;
      Standard_Boolean Ok;
      Xp = b;
      do
      {
        Ok = F.Values(Xp, Yp, Dp);
        if (Ok)
        {
          Ok = Standard_False;
          if (Dp > 1e-10 || Dp < -1e-10)
          {
            Xp = Xp - (Yp - K) / Dp;
          }
          if (Xp <= x2 && Xp >= x1)
          {
            F.Value(Xp, Yp);
            Yp -= K;
            if (Abs(Yp) < Abs(fb))
            {
              b  = Xp;
              fb = Yp;
              Ok = Standard_True;
            }
          }
        }
      } while (Ok && --itern >= 0);

      AppendRoot(Sol, NbStateSol, b, F, K, dX);
      return;
    }
    if (Abs(e) >= tol1 && Abs(fa) > Abs(fb))
    {
      s = fb / fa;
      if (a == c)
      {
        p = xm * s;
        p += p;
        q = 1.0 - s;
      }
      else
      {
        q = fa / fc;
        r = fb / fc;
        p = s * ((xm + xm) * q * (q - r) - (b - a) * (r - 1.0));
        q = (q - 1.0) * (r - 1.0) * (s - 1.0);
      }
      if (p > 0.0)
      {
        q = -q;
      }
      p    = Abs(p);
      min1 = 3.0 * xm * q - Abs(tol1 * q);
      min2 = Abs(e * q);
      if ((p + p) < ((min1 < min2) ? min1 : min2))
      {
        e = d;
        d = p / q;
      }
      else
      {
        d = xm;
        e = d;
      }
    }
    else
    {
      d = xm;
      e = d;
    }
    a  = b;
    fa = fb;
    if (Abs(d) > tol1)
    {
      b += d;
    }
    else
    {
      if (xm >= 0)
        b += Abs(tol1);
      else
        b += -Abs(tol1);
    }
    F.Value(b, fb);
    fb -= K;
  }
#ifdef OCCT_DEBUG
  std::cout << " Non Convergence dans math_FunctionRoots.cxx " << std::endl;
#endif
}

#define NEWSEQ 1

#define MATH_FUNCTIONROOTS_NEWCODE // Nv Traitement

// #define MATH_FUNCTIONROOTS_OLDCODE // Ancien
// #define MATH_FUNCTIONROOTS_CHECK // Check

math_FunctionRoots::math_FunctionRoots(math_FunctionWithDerivative& F,
                                       const Standard_Real          A,
                                       const Standard_Real          B,
                                       const Standard_Integer       NbSample,
                                       const Standard_Real          _EpsX,
                                       const Standard_Real          EpsF,
                                       const Standard_Real          EpsNull,
                                       const Standard_Real          K)
{
#ifdef OCCT_DEBUG
  if (myDebug)
  {
    std::cout << "---- Debut de math_FunctionRoots ----" << std::endl;
    nbsolve = 0;
  }
#endif

#if NEWSEQ
  TColStd_SequenceOfReal StaticSol;
#endif
  Sol.Clear();
  NbStateSol.Clear();
#ifdef MATH_FUNCTIONROOTS_NEWCODE
  {
    Done                = Standard_True;
    Standard_Real    X0 = A;
    Standard_Real    XN = B;
    Standard_Integer N  = NbSample;
    //-- ------------------------------------------------------------
    //-- Verifications de bas niveau
    if (B < A)
    {
      X0 = B;
      XN = A;
    }
    N *= 2;
    if (N < 20)
    {
      N = 20;
    }
    //--  On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
    Standard_Real EpsX   = _EpsX;
    Standard_Real DeltaU = Abs(X0) + Abs(XN);
    Standard_Real NEpsX  = 0.0000000001 * DeltaU;
    if (EpsX < NEpsX)
    {
      EpsX = NEpsX;
    }

    //-- recherche d un intervalle ou F(xi) et F(xj) sont de signes differents
    //-- A .............................................................. B
    //-- X0   X1   X2 ........................................  Xn-1      Xn
    Standard_Integer     i;
    Standard_Real        X = X0;
    Standard_Boolean     Ok;
    double               dx = (XN - X0) / N;
    TColStd_Array1OfReal ptrval(0, N);
    Standard_Integer     Nvalid = -1;
    Standard_Real        aux    = 0;
    for (i = 0; i <= N; i++, X += dx)
    {
      if (X > XN)
        X = XN;
      Ok = F.Value(X, aux);
      if (Ok)
        ptrval(++Nvalid) = aux - K;
      //      ptrval(i)-=K;
    }
    //-- Toute la fonction est nulle ?

    if (Nvalid < N)
    {
      Done = Standard_False;
      return;
    }

    AllNull = Standard_True;
    //    for(i=0;AllNull && i<=N;i++) {
    for (i = 0; AllNull && i <= N; i++)
    {
      if (ptrval(i) > EpsNull || ptrval(i) < -EpsNull)
      {
        AllNull = Standard_False;
      }
    }
    if (AllNull)
    {
      //-- tous les points echantillons sont dans la tolerance
    }
    else
    {
      //-- Il y a des points hors tolerance
      //-- on detecte les changements de signes STRICTS
      Standard_Integer ip1;
      //      Standard_Boolean chgtsign=Standard_False;
      Standard_Real tol = EpsX;
      Standard_Real X2;
      for (i = 0, ip1 = 1, X = X0; i < N; i++, ip1++, X += dx)
      {
        X2 = X + dx;
        if (X2 > XN)
          X2 = XN;
        if (ptrval(i) < 0.0)
        {
          if (ptrval(ip1) > 0.0)
          {
            //-- --------------------------------------------------
            //-- changement de signe dans Xi Xi+1
            Solve(F, K, X, ptrval(i), X2, ptrval(ip1), tol, NEpsX, Sol, NbStateSol);
          }
        }
        else
        {
          if (ptrval(ip1) < 0.0)
          {
            //-- --------------------------------------------------
            //-- changement de signe dans Xi Xi+1
            Solve(F, K, X, ptrval(i), X2, ptrval(ip1), tol, NEpsX, Sol, NbStateSol);
          }
        }
      }
      //-- On detecte les cas ou la fct s annule sur des Xi et est
      //-- non nulle au voisinage de Xi
      //--
      //-- On prend 2 points u0,u1 au voisinage de Xi
      //-- Si (F(u0)-K)*(F(u1)-K) <0   on lance une recherche
      //-- Sinon si (F(u0)-K)*(F(u1)-K) !=0 on insere le point X
      for (i = 0; i <= N; i++)
      {
        if (ptrval(i) == 0)
        {
          //	  Standard_Real Val,Deriv;
          X = X0 + i * dx;
          if (X > XN)
            X = XN;
          Standard_Real u0, u1;
          u0 = dx * 0.5;
          u1 = X + u0;
          u0 += X;
          if (u0 < X0)
            u0 = X0;
          if (u0 > XN)
            u0 = XN;
          if (u1 < X0)
            u1 = X0;
          if (u1 > XN)
            u1 = XN;

          Standard_Real y0, y1;
          F.Value(u0, y0);
          y0 -= K;
          F.Value(u1, y1);
          y1 -= K;
          if (y0 * y1 < 0.0)
          {
            Solve(F, K, u0, y0, u1, y1, tol, NEpsX, Sol, NbStateSol);
          }
          else
          {
            if (y0 != 0.0 || y1 != 0.0)
            {
              AppendRoot(Sol, NbStateSol, X, F, K, NEpsX);
            }
          }
        }
      }
      //-- --------------------------------------------------------------------------------
      //-- Il faut traiter differement le cas des points en bout :
      if (ptrval(0) <= EpsF && ptrval(0) >= -EpsF)
      {
        AppendRoot(Sol, NbStateSol, X0, F, K, NEpsX);
      }
      if (ptrval(N) <= EpsF && ptrval(N) >= -EpsF)
      {
        AppendRoot(Sol, NbStateSol, XN, F, K, NEpsX);
      }

      //-- --------------------------------------------------------------------------------
      //-- --------------------------------------------------------------------------------
      //-- On detecte les zones ou on a sur les points echantillons un minimum avec f(x)>0
      //--                                                          un maximum avec f(x)<0
      //-- On reprend une discretisation plus fine au voisinage de ces extremums
      //--
      //-- Recherche d un minima positif
      Standard_Real    xm, ym, dym, xm1, xp1;
      Standard_Real    majdx = 5.0 * dx;
      Standard_Boolean Rediscr;
      //      Standard_Real ptrvalbis[MAXBIS];
      Standard_Integer im1 = 0;
      ip1                  = 2;
      for (i = 1, xm = X0 + dx; i < N; xm += dx, i++, im1++, ip1++)
      {
        Rediscr = Standard_False;
        if (xm > XN)
          xm = XN;
        if (ptrval(i) > 0.0)
        {
          if ((ptrval(im1) > ptrval(i)) && (ptrval(ip1) > ptrval(i)))
          {
            //-- Peut on traverser l axe Ox
            //-- -------------- Estimation a partir de Xim1
            xm1 = xm - dx;
            if (xm1 < X0)
              xm1 = X0;
            F.Values(xm1, ym, dym);
            ym -= K;
            if (dym < -1e-10 || dym > 1e-10)
            {                             // normalement dym < 0
              Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
              if (t < majdx && t > -majdx)
              {
                Rediscr = Standard_True;
              }
            }
            //-- -------------- Estimation a partir de Xip1
            if (Rediscr == Standard_False)
            {
              xp1 = xm + dx;
              if (xp1 > XN)
                xp1 = XN;
              F.Values(xp1, ym, dym);
              ym -= K;
              if (dym < -1e-10 || dym > 1e-10)
              {                             // normalement dym > 0
                Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
                if (t < majdx && t > -majdx)
                {
                  Rediscr = Standard_True;
                }
              }
            }
          }
        }
        else if (ptrval(i) < 0.0)
        {
          if ((ptrval(im1) < ptrval(i)) && (ptrval(ip1) < ptrval(i)))
          {
            //-- Peut on traverser l axe Ox
            //-- -------------- Estimation a partir de Xim1
            xm1 = xm - dx;
            if (xm1 < X0)
              xm1 = X0;
            F.Values(xm1, ym, dym);
            ym -= K;
            if (dym > 1e-10 || dym < -1e-10)
            {                             // normalement dym > 0
              Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
              if (t < majdx && t > -majdx)
              {
                Rediscr = Standard_True;
              }
            }
            //-- -------------- Estimation a partir de Xim1
            if (Rediscr == Standard_False)
            {
              xm1 = xm - dx;
              if (xm1 < X0)
                xm1 = X0;
              F.Values(xm1, ym, dym);
              ym -= K;
              if (dym > 1e-10 || dym < -1e-10)
              {                             // normalement dym < 0
                Standard_Real t = ym / dym; //-- t=xm-x* = (ym-0)/dym
                if (t < majdx && t > -majdx)
                {
                  Rediscr = Standard_True;
                }
              }
            }
          }
        }
        if (Rediscr)
        {
          Standard_Real x0 = xm - dx;
          Standard_Real x3 = xm + dx;
          if (x0 < X0)
            x0 = X0;
          if (x3 > XN)
            x3 = XN;
          Standard_Real    aSolX1 = 0., aSolX2 = 0.;
          Standard_Real    aVal1 = 0., aVal2 = 0.;
          Standard_Real    aDer1 = 0., aDer2 = 0.;
          Standard_Boolean isSol1 = Standard_False;
          Standard_Boolean isSol2 = Standard_False;
          //-- ----------------------------------------------------
          //-- Find minimum of the function |F| between x0 and x3
          //-- by searching for the zero of the function derivative
          DerivFunction      aDerF(F);
          math_BracketedRoot aBR(aDerF, x0, x3, _EpsX);
          if (aBR.IsDone())
          {
            aSolX1 = aBR.Root();
            F.Value(aSolX1, aVal1);
            aVal1 = Abs(aVal1);
            if (aVal1 < EpsF)
            {
              isSol1 = Standard_True;
              aDer1  = aBR.Value();
            }
          }

          //-- --------------------------------------------------
          //-- On recherche un extrema entre x0 et x3
          //-- x1 et x2 sont tels que x0<x1<x2<x3
          //-- et |f(x0)| > |f(x1)|   et |f(x3)| > |f(x2)|
          //--
          //-- En entree : a=xm-dx  b=xm c=xm+dx
          Standard_Real x1, x2, f0, f3;
          Standard_Real R                    = 0.61803399;
          Standard_Real C                    = 1.0 - R;
          Standard_Real tolCR                = NEpsX * 10.0;
          f0                                 = ptrval(im1);
          f3                                 = ptrval(ip1);
          Standard_Boolean recherche_minimum = (f0 > 0.0);

          if (Abs(x3 - xm) > Abs(x0 - xm))
          {
            x1 = xm;
            x2 = xm + C * (x3 - xm);
          }
          else
          {
            x2 = xm;
            x1 = xm - C * (xm - x0);
          }
          Standard_Real f1, f2;
          F.Value(x1, f1);
          f1 -= K;
          F.Value(x2, f2);
          f2 -= K;
          //-- printf("\n *************** RECHERCHE MINIMUM **********\n");
          Standard_Real tolX = 0.001 * NEpsX;
          while (Abs(x3 - x0) > tolCR * (Abs(x1) + Abs(x2)) && (Abs(x1 - x2) > tolX))
          {
            //-- printf("\n (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) (%10.5g,%10.5g) ",
            //--    x0,f0,x1,f1,x2,f2,x3,f3);
            if (recherche_minimum)
            {
              if (f2 < f1)
              {
                x0 = x1;
                x1 = x2;
                x2 = R * x1 + C * x3;
                f0 = f1;
                f1 = f2;
                F.Value(x2, f2);
                f2 -= K;
              }
              else
              {
                x3 = x2;
                x2 = x1;
                x1 = R * x2 + C * x0;
                f3 = f2;
                f2 = f1;
                F.Value(x1, f1);
                f1 -= K;
              }
            }
            else
            {
              if (f2 > f1)
              {
                x0 = x1;
                x1 = x2;
                x2 = R * x1 + C * x3;
                f0 = f1;
                f1 = f2;
                F.Value(x2, f2);
                f2 -= K;
              }
              else
              {
                x3 = x2;
                x2 = x1;
                x1 = R * x2 + C * x0;
                f3 = f2;
                f2 = f1;
                F.Value(x1, f1);
                f1 -= K;
              }
            }
            //-- On ne fait pas que chercher des extremas. Il faut verifier
            //-- si on ne tombe pas sur une racine
            if (f1 * f0 < 0.0)
            {
              //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x0,f0,x1,f1);
              Solve(F, K, x0, f0, x1, f1, tol, NEpsX, Sol, NbStateSol);
            }
            if (f2 * f3 < 0.0)
            {
              //-- printf("\n Recherche entre  (%10.5g,%10.5g) (%10.5g,%10.5g) ",x2,f2,x3,f3);
              Solve(F, K, x2, f2, x3, f3, tol, NEpsX, Sol, NbStateSol);
            }
          }
          if ((recherche_minimum && f1 < f2) || (!recherche_minimum && f1 > f2))
          {
            //-- x1,f(x1) minimum
            if (Abs(f1) < EpsF)
            {
              isSol2 = Standard_True;
              aSolX2 = x1;
              aVal2  = Abs(f1);
            }
          }
          else
          {
            //-- x2.f(x2) minimum
            if (Abs(f2) < EpsF)
            {
              isSol2 = Standard_True;
              aSolX2 = x2;
              aVal2  = Abs(f2);
            }
          }
          // Choose the best solution between aSolX1, aSolX2
          if (isSol1 && isSol2)
          {
            if (aVal2 - aVal1 > EpsF)
              AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
            else if (aVal1 - aVal2 > EpsF)
              AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
            else
            {
              aDer1 = Abs(aDer1);
              F.Derivative(aSolX2, aDer2);
              aDer2 = Abs(aDer2);
              if (aDer1 < aDer2)
                AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
              else
                AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
            }
          }
          else if (isSol1)
            AppendRoot(Sol, NbStateSol, aSolX1, F, K, NEpsX);
          else if (isSol2)
            AppendRoot(Sol, NbStateSol, aSolX2, F, K, NEpsX);
        } //-- Recherche d un extrema
      } //-- for
    }

  #if NEWSEQ
    #ifdef MATH_FUNCTIONROOTS_CHECK
    {
      StaticSol.Clear();
      Standard_Integer n = Sol.Length();
      for (Standard_Integer ii = 1; ii <= n; ii++)
      {
        StaticSol.Append(Sol.Value(ii));
      }
      Sol.Clear();
      NbStateSol.Clear();
    }
    #endif
  #endif
#endif
  }
#ifdef MATH_FUNCTIONROOTS_OLDCODE
  {
    //-- ********************************************************************************
    //--                              ANCIEN TRAITEMENT
    //-- ********************************************************************************

    // calculate all the real roots of a function within the range
    // A..B. without condition on A and B
    // a solution X is found when
    //   abs(Xi - Xi-1) <= EpsX and abs(F(Xi)-K) <= Epsf.
    // The function is considered as null between A and B if
    // abs(F-K) <= EpsNull within this range.
    Standard_Real EpsX = _EpsX; //-- Cas ou le parametre va de 100000000 a 1000000001
                                //-- Il ne faut pas EpsX = 0.000...001  car dans ce cas
                                //-- U + Nn*EpsX     ==     U
    Standard_Real    Lowr, Upp;
    Standard_Real    Increment;
    Standard_Real    Null2;
    Standard_Real    FLowr, FUpp, DFLowr, DFUpp;
    Standard_Real    U, Xu;
    Standard_Real    Fxu, DFxu, FFxu, DFFxu;
    Standard_Real    Fyu, DFyu, FFyu, DFFyu;
    Standard_Boolean Finish;
    Standard_Real    FFi;
    Standard_Integer Nbiter = 30;
    Standard_Integer Iter;
    Standard_Real    Ambda, T;
    Standard_Real    AA, BB, CC;
    Standard_Integer Nn;
    Standard_Real    Alfa1              = 0, Alfa2;
    Standard_Real    OldDF              = RealLast();
    Standard_Real    Standard_Underflow = 1e-32; //-- RealSmall();
    Standard_Boolean Ok;

    Done = Standard_False;

    StdFail_NotDone_Raise_if(NbSample <= 0, " ");

    // initialisation

    if (A > B)
    {
      Lowr = B;
      Upp  = A;
    }
    else
    {
      Lowr = A;
      Upp  = B;
    }

    Increment = (Upp - Lowr) / NbSample;
    StdFail_NotDone_Raise_if(Increment < EpsX, " ");
    Done = Standard_True;
    //--  On teste si EpsX est trop petit (ie : U+Nn*EpsX == U )
    Standard_Real DeltaU = Abs(Upp) + Abs(Lowr);
    Standard_Real NEpsX  = 0.0000000001 * DeltaU;
    if (EpsX < NEpsX)
    {
      EpsX = NEpsX;
      //-- std::cout<<" \n EpsX Init = "<<_EpsX<<" devient : (deltaU : "<<DeltaU<<" )   EpsX =
      //"<<EpsX<<std::endl;
    }
    //--
    Null2 = EpsNull * EpsNull;

    Ok = F.Values(Lowr, FLowr, DFLowr);

    if (!Ok)
    {
      Done = Standard_False;
      return;
    }

    FLowr = FLowr - K;

    Ok = F.Values(Upp, FUpp, DFUpp);

    if (!Ok)
    {
      Done = Standard_False;
      return;
    }

    FUpp = FUpp - K;

    // Calcul sur U

    U     = Lowr - EpsX;
    Fyu   = FLowr - EpsX * DFLowr; // extrapolation lineaire
    DFyu  = DFLowr;
    FFyu  = Fyu * Fyu;
    DFFyu = Fyu * DFyu;
    DFFyu += DFFyu;
    AllNull = (FFyu <= Null2);

    while (U < Upp)
    {

      Xu    = U;
      Fxu   = Fyu;
      DFxu  = DFyu;
      FFxu  = FFyu;
      DFFxu = DFFyu;

      U = Xu + Increment;
      if (U <= Lowr)
      {
        Fyu  = FLowr + (U - Lowr) * DFLowr;
        DFyu = DFLowr;
      }
      else if (U >= Upp)
      {
        Fyu  = FUpp + (U - Upp) * DFUpp;
        DFyu = DFUpp;
      }
      else
      {
        Ok = F.Values(U, Fyu, DFyu);

        if (!Ok)
        {
          Done = Standard_False;
          return;
        }

        Fyu = Fyu - K;
      }
      FFyu  = Fyu * Fyu;
      DFFyu = Fyu * DFyu;
      DFFyu += DFFyu; //-- DFFyu = 2.*Fyu*DFyu;

      if (!AllNull || (FFyu > Null2 && U <= Upp))
      {

        if (AllNull)
        { // rechercher les vraix zeros depuis le debut

          AllNull = Standard_False;
          Xu      = Lowr - EpsX;
          Fxu     = FLowr - EpsX * DFLowr;
          DFxu    = DFLowr;
          FFxu    = Fxu * Fxu;
          DFFxu   = Fxu * DFxu;
          DFFxu += DFFxu; //-- DFFxu = 2.*Fxu*DFxu;
          U  = Xu + Increment;
          Ok = F.Values(U, Fyu, DFyu);

          if (!Ok)
          {
            Done = Standard_False;
            return;
          }

          Fyu   = Fyu - K;
          FFyu  = Fyu * Fyu;
          DFFyu = Fyu * DFyu;
          DFFyu += DFFyu; //-- DFFyu = 2.*Fyu*DFyu;
        }
        Standard_Real FxuFyu = Fxu * Fyu;

        if ((DFFyu > 0. && DFFxu <= 0.) || (DFFyu < 0. && FFyu >= FFxu && DFFxu <= 0.)
            || (DFFyu > 0. && FFyu <= FFxu && DFFxu >= 0.) || (FxuFyu <= 0.))
        {
          // recherche d 1 minimun possible
          Finish = Standard_False;
          Ambda  = Increment;
          T      = 0.;
          Iter   = 0;
          FFi    = 0.;

          if (FxuFyu > 0.)
          {
            // chercher si f peut s annuler pour eviter
            //  des iterations inutiles
            if (Fxu * (Fxu + 2. * DFxu * Increment) > 0.
                && Fyu * (Fyu - 2. * DFyu * Increment) > 0.)
            {

              Finish = Standard_True;
              FFi    = Min(FFxu, FFyu); // pour ne pas recalculer yu
            }
            else if ((DFFxu <= Standard_Underflow && -DFFxu <= Standard_Underflow)
                     || (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow))
            {

              Finish = Standard_True;
              FFxu   = 0.0;
              FFi    = FFyu; // pour recalculer yu
            }
            else if ((DFFyu <= Standard_Underflow && -DFFyu <= Standard_Underflow)
                     || (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow))
            {

              Finish = Standard_True;
              FFyu   = 0.0;
              FFi    = FFxu; // pour recalculer U
            }
          }
          else if (FFxu <= Standard_Underflow && -FFxu <= Standard_Underflow)
          {

            Finish = Standard_True;
            FFxu   = 0.0;
            FFi    = FFyu;
          }
          else if (FFyu <= Standard_Underflow && -FFyu <= Standard_Underflow)
          {

            Finish = Standard_True;
            FFyu   = 0.0;
            FFi    = FFxu;
          }
          while (!Finish)
          {

            // calcul des 2 solutions annulant la derivee de l interpolation cubique
            //    Ambda*t=(U-Xu)  F(t)=aa*t*t*t/3+bb*t*t+cc*t+d
            //    df=aa*t*t+2*bb*t+cc

            AA = 3. * (Ambda * (DFFxu + DFFyu) + 2. * (FFxu - FFyu));
            BB = -2 * (Ambda * (DFFyu + 2. * DFFxu) + 3. * (FFxu - FFyu));
            CC = Ambda * DFFxu;

            if (Abs(AA) < 1e-14 && Abs(BB) < 1e-14 && Abs(CC) < 1e-14)
            {
              AA = BB = CC = 0;
            }

            math_DirectPolynomialRoots Solv(AA, BB, CC);
            if (!Solv.InfiniteRoots())
            {
              Nn = Solv.NbSolutions();
              if (Nn <= 0)
              {
                if (Fxu * Fyu < 0.)
                {
                  Alfa1 = 0.5;
                }
                else
                  Finish = Standard_True;
              }
              else
              {
                Alfa1 = Solv.Value(1);
                if (Nn == 2)
                {
                  Alfa2 = Solv.Value(2);
                  if (Alfa1 > Alfa2)
                  {
                    AA    = Alfa1;
                    Alfa1 = Alfa2;
                    Alfa2 = AA;
                  }
                  if (Alfa1 > 1. || Alfa2 < 0.)
                  {
                    // resolution par dichotomie
                    if (Fxu * Fyu < 0.)
                      Alfa1 = 0.5;
                    else
                      Finish = Standard_True;
                  }
                  else if (Alfa1 < 0. || (DFFxu > 0. && DFFyu >= 0.))
                  {
                    // si 2 derivees >0
                    //(cas changement de signe de la distance signee sans
                    // changement de signe de la derivee:
                    // cas de 'presque'tangence avec 2
                    // solutions proches) ,on prend la plus grane racine
                    if (Alfa2 > 1.)
                    {
                      if (Fxu * Fyu < 0.)
                        Alfa1 = 0.5;
                      else
                        Finish = Standard_True;
                    }
                    else
                      Alfa1 = Alfa2;
                  }
                }
              }
            }
            else if (Fxu * Fyu < -1e-14)
              Alfa1 = 0.5;
            //-- else if (Fxu*Fyu < 0.) Alfa1 = 0.5;
            else
              Finish = Standard_True;

            if (!Finish)
            {
              // petits tests pour diminuer le nombre d iterations
              if (Alfa1 <= EpsX)
              {
                Alfa1 += Alfa1;
              }
              else if (Alfa1 >= (1. - EpsX))
              {
                Alfa1 = Alfa1 + Alfa1 - 1.;
              }
              Alfa1 = Ambda * (1. - Alfa1);
              U     = U + T - Alfa1;
              if (U <= Lowr)
              {
                AA = FLowr + (U - Lowr) * DFLowr;
                BB = DFLowr;
              }
              else if (U >= Upp)
              {
                AA = FUpp + (U - Upp) * DFUpp;
                BB = DFUpp;
              }
              else
              {
                Ok = F.Values(U, AA, BB);

                if (!Ok)
                {
                  Done = Standard_False;
                  return;
                }

                AA = AA - K;
              }
              FFi = AA * AA;
              CC  = AA * BB;
              CC += CC;
              if (((CC < 0. && FFi < FFxu) || DFFxu > 0.) && AA * Fxu > 0.)
              {
                FFxu  = FFi;
                DFFxu = CC;
                Fxu   = AA;
                DFxu  = BB;
                T     = Alfa1;
                if (Alfa1 > Ambda * 0.5)
                {
                  // remarque (1)
                  // determination d 1 autre borne pour diviser
                  // le nouvel intervalle par 2 au -
                  Xu = U + Alfa1 * 0.5;
                  if (Xu <= Lowr)
                  {
                    AA = FLowr + (Xu - Lowr) * DFLowr;
                    BB = DFLowr;
                  }
                  else if (Xu >= Upp)
                  {
                    AA = FUpp + (Xu - Upp) * DFUpp;
                    BB = DFUpp;
                  }
                  else
                  {
                    Ok = F.Values(Xu, AA, BB);

                    if (!Ok)
                    {
                      Done = Standard_False;
                      return;
                    }

                    AA = AA - K;
                  }
                  FFi = AA * AA;
                  CC  = AA * BB;
                  CC += CC;
                  if (((CC >= 0. || FFi >= FFxu) && DFFxu < 0.) || Fxu * AA < 0.)
                  {
                    Fyu   = AA;
                    DFyu  = BB;
                    FFyu  = FFi;
                    DFFyu = CC;
                    T     = Alfa1 * 0.5;
                    Ambda = Alfa1 * 0.5;
                  }
                  else if (AA * Fyu < 0. && AA * Fxu > 0.)
                  {
                    // changement de signe sur l intervalle u,U
                    Fxu   = AA;
                    DFxu  = BB;
                    FFxu  = FFi;
                    DFFxu = CC;
                    FFi   = Min(FFxu, FFyu);
                    T     = Alfa1 * 0.5;
                    Ambda = Alfa1 * 0.5;
                    U     = Xu;
                  }
                  else
                  {
                    Ambda = Alfa1;
                  }
                }
                else
                {
                  Ambda = Alfa1;
                }
              }
              else
              {
                Fyu   = AA;
                DFyu  = BB;
                FFyu  = FFi;
                DFFyu = CC;
                if ((Ambda - Alfa1) > Ambda * 0.5)
                {
                  // meme remarque (1)
                  Xu = U - (Ambda - Alfa1) * 0.5;
                  if (Xu <= Lowr)
                  {
                    AA = FLowr + (Xu - Lowr) * DFLowr;
                    BB = DFLowr;
                  }
                  else if (Xu >= Upp)
                  {
                    AA = FUpp + (Xu - Upp) * DFUpp;
                    BB = DFUpp;
                  }
                  else
                  {
                    Ok = F.Values(Xu, AA, BB);

                    if (!Ok)
                    {
                      Done = Standard_False;
                      return;
                    }

                    AA = AA - K;
                  }
                  FFi = AA * AA;
                  CC  = AA * BB;
                  CC += CC;
                  if (AA * Fyu <= 0. && AA * Fxu > 0.)
                  {
                    FFxu  = FFi;
                    DFFxu = CC;
                    Fxu   = AA;
                    DFxu  = BB;
                    Ambda = (Ambda - Alfa1) * 0.5;
                    T     = 0.;
                  }
                  else if (AA * Fxu < 0. && AA * Fyu > 0.)
                  {
                    FFyu  = FFi;
                    DFFyu = CC;
                    Fyu   = AA;
                    DFyu  = BB;
                    Ambda = (Ambda - Alfa1) * 0.5;
                    T     = 0.;
                    FFi   = Min(FFxu, FFyu);
                    U     = Xu;
                  }
                  else
                  {
                    T     = 0.;
                    Ambda = Ambda - Alfa1;
                  }
                }
                else
                {
                  T     = 0.;
                  Ambda = Ambda - Alfa1;
                }
              }
              // tests d arrets
              if (Abs(FFxu) <= Standard_Underflow
                  || (Abs(DFFxu) <= Standard_Underflow && Fxu * Fyu > 0.))
              {
                Finish = Standard_True;
                if (Abs(FFxu) <= Standard_Underflow)
                {
                  FFxu = 0.0;
                }
                FFi = FFyu;
              }
              else if (Abs(FFyu) <= Standard_Underflow
                       || (Abs(DFFyu) <= Standard_Underflow && Fxu * Fyu > 0.))
              {
                Finish = Standard_True;
                if (Abs(FFyu) <= Standard_Underflow)
                {
                  FFyu = 0.0;
                }
                FFi = FFxu;
              }
              else
              {
                Iter = Iter + 1;
                Finish =
                  Iter >= Nbiter || (Ambda <= EpsX && (Fxu * Fyu >= 0. || FFi <= EpsF * EpsF));
              }
            }
          } // fin interpolation cubique

          // restitution du meilleur resultat

          if (FFxu < FFi)
          {
            U = U + T - Ambda;
          }
          else if (FFyu < FFi)
          {
            U = U + T;
          }

          if (U >= (Lowr - EpsX) && U <= (Upp + EpsX))
          {
            U   = Max(Lowr, Min(U, Upp));
            Ok  = F.Value(U, FFi);
            FFi = FFi - K;
            if (Abs(FFi) < EpsF)
            {
              // coherence
              if (Abs(Fxu) <= Standard_Underflow)
              {
                AA = DFxu;
              }
              else if (Abs(Fyu) <= Standard_Underflow)
              {
                AA = DFyu;
              }
              else if (Fxu * Fyu > 0.)
              {
                AA = 0.;
              }
              else
              {
                AA = Fyu - Fxu;
              }
              if (!Sol.IsEmpty())
              {
                if (Abs(Sol.Last() - U) > 5. * EpsX || (OldDF != RealLast() && OldDF * AA < 0.))
                {
                  Sol.Append(U);
                  NbStateSol.Append(F.GetStateNumber());
                }
              }
              else
              {
                Sol.Append(U);
                NbStateSol.Append(F.GetStateNumber());
              }
              OldDF = AA;
            }
          }
          DFFyu = 0.;
          Fyu   = 0.;
          Nn    = 1;
          while (Nn < 1000000 && DFFyu <= 0.)
          {
            // on repart d 1 pouyem plus loin
            U = U + Nn * EpsX;
            if (U <= Lowr)
            {
              AA = FLowr + (U - Lowr) * DFLowr;
              BB = DFLowr;
            }
            else if (U >= Upp)
            {
              AA = FUpp + (U - Upp) * DFUpp;
              BB = DFUpp;
            }
            else
            {
              Ok = F.Values(U, AA, BB);
              AA = AA - K;
            }
            if (AA * Fyu < 0.)
            {
              U  = U - Nn * EpsX;
              Nn = 1000001;
            }
            else
            {
              Fyu   = AA;
              DFyu  = BB;
              FFyu  = AA * AA;
              DFFyu = 2. * DFyu * Fyu;
              Nn    = 2 * Nn;
            }
          }
        }
      }
    }
  #if NEWSEQ
    #ifdef MATH_FUNCTIONROOTS_CHECK
    {
      Standard_Integer n1 = StaticSol.Length();
      Standard_Integer n2 = Sol.Length();
      if (n1 != n2)
      {
        printf("\n mathFunctionRoots : n1=%d  n2=%d EpsF=%g EpsX=%g\n", n1, n2, EpsF, NEpsX);
        for (Standard_Integer x1 = 1; x1 <= n1; x1++)
        {
          Standard_Real v;
          F.Value(StaticSol(x1), v);
          v -= K;
          printf(" (%+13.8g:%+13.8g) ", StaticSol(x1), v);
        }
        printf("\n");
        for (Standard_Integer x2 = 1; x2 <= n2; x2++)
        {
          Standard_Real v;
          F.Value(Sol(x2), v);
          v -= K;
          printf(" (%+13.8g:%+13.8g) ", Sol(x2), v);
        }
        printf("\n");
      }
      Standard_Integer n = n1;
      if (n1 > n2)
        n = n2;
      for (Standard_Integer i = 1; i <= n; i++)
      {
        Standard_Real t = Sol(i) - StaticSol(i);
        if (Abs(t) > NEpsX)
        {
          printf("\n mathFunctionRoots : i:%d/%d  delta: %g", i, n, t);
        }
      }
    }
    #endif
  #endif
  }
#endif
}

void math_FunctionRoots::Dump(Standard_OStream& o) const
{

  o << "math_FunctionRoots ";
  if (Done)
  {
    o << " Status = Done \n";
    o << " Number of solutions = " << Sol.Length() << std::endl;
    for (Standard_Integer i = 1; i <= Sol.Length(); i++)
    {
      o << " Solution Number " << i << "= " << Sol.Value(i) << std::endl;
    }
  }
  else
  {
    o << " Status = not Done \n";
  }
}