// Copyright (c) 1995-1999 Matra Datavision
// Copyright (c) 1999-2012 OPEN CASCADE SAS
//
// The content of this file is subject to the Open CASCADE Technology Public
// License Version 6.5 (the "License"). You may not use the content of this file
// except in compliance with the License. Please obtain a copy of the License
// at http://www.opencascade.org and read it completely before using this file.
//
// The Initial Developer of the Original Code is Open CASCADE S.A.S., having its
// main offices at: 1, place des Freres Montgolfier, 78280 Guyancourt, France.
//
// The Original Code and all software distributed under the License is
// distributed on an "AS IS" basis, without warranty of any kind, and the
// Initial Developer hereby disclaims all such warranties, including without
// limitation, any warranties of merchantability, fitness for a particular
// purpose or non-infringement. Please see the License for the specific terms
// and conditions governing the rights and limitations under the License.

//Modified by MPS : 21-05-97   PRO 7598 
//                  Le nombre de solutions rendu est mySqDist.Length() et non
//                  mySqDist.Length()/2
//                  ajout des valeurs absolues dans le test d'orthogonalite de
//                  GetStateNumber()

/*-----------------------------------------------------------------------------
 Fonctions permettant de rechercher une distance extremale entre 2 courbes
C1 et C2 (en partant de points approches C1(u0) et C2(v0)).
 Cette classe herite de math_FunctionSetWithDerivatives et est utilisee par
l'algorithme math_FunctionSetRoot.
 Si on note Du et Dv, les derivees en u et v, les 2 fonctions a annuler sont:
  { F1(u,v) = (C2(v)-C1(u)).Du(u)/||Du|| }
  { F2(u,v) = (C2(v)-C1(u)).Dv(v)||Dv|| }
 Si on note Duu et Dvv, les derivees secondes de C1 et C2, les derivees de F1
et F2 sont egales a:
  { Duf1(u,v) = -||Du|| + C1C2.Duu/||Du||- F1(u,v)*Duu*Du/||Du||**2
  { Dvf1(u,v) = Dv.Du/||Du|| 
  { Duf2(u,v) = -Du.Dv/||Dv||
  { Dvf2(u,v) = ||Dv|| + C2C1.Dvv/||Dv||- F2(u,v)*Dv*Dvv/||Dv||**2

----------------------------------------------------------------------------*/

#include <Precision.hxx>


static const Standard_Real MinTol    = 1.e-20;
static const Standard_Real TolFactor = 1.e-12;
static const Standard_Real MinStep   = 1e-7;

static const Standard_Integer MaxOrder = 3;



//=============================================================================
Standard_Real Extrema_FuncExtCC::SearchOfTolerance(const Standard_Address C)
  {
  const Standard_Integer NPoint = 10;
  Standard_Real aStartParam, anEndParam;
  
  if(C==myC1)
    {
    aStartParam = myUinfium;
    anEndParam = myUsupremum;
    }
  else if(C==myC2)
    {
    aStartParam = myVinfium;
    anEndParam = myVsupremum;
    }
  else
    {
#ifdef DEB
    cout << "+++ Function Extrema_FuncExtCC::SearchOfTolerance(...)" << endl;
    cout << "Warning: No curve for tolerance computing!---"<<endl;
#endif
    return MinTol;
    }
    
  const Standard_Real aStep = (anEndParam - aStartParam)/(Standard_Real)NPoint;
    
  Standard_Integer aNum = 0;
  Standard_Real aMax = -Precision::Infinite();  //Maximum value of 1st derivative
                                                //(it is computed with using NPoint point)
    
  do
    {
    Standard_Real u = aStartParam + aNum*aStep;    //parameter for every point
    if(u > anEndParam)
      u = anEndParam;
    
    Pnt Ptemp;  //empty point (is not used below)
    Vec VDer;   // 1st derivative vector
    Tool1::D1(*((Curve1*)C), u, Ptemp, VDer);
    Standard_Real vm = VDer.Magnitude();
    if(vm > aMax)
      aMax = vm;
    }
  while(++aNum < NPoint+1);
  
  return Max(aMax*TolFactor,MinTol);
  }

//=============================================================================
Extrema_FuncExtCC::Extrema_FuncExtCC(const Standard_Real thetol) : myC1 (0), myC2 (0), myTol (thetol)
  {
  math_Vector V1(1,2), V2(1,2);
  V1(1) = 0.0;
  V2(1) = 0.0;
  V1(2) = 0.0;
  V2(2) = 0.0;
  SubIntervalInitialize(V1, V2);
  myMaxDerivOrderC1 = 0;
  myTolC1=MinTol;
  myMaxDerivOrderC2 = 0;
  myTolC2=MinTol;
  }

//=============================================================================
Extrema_FuncExtCC::Extrema_FuncExtCC (const Curve1& C1,
				      const Curve2& C2,
				      const Standard_Real thetol) :
				        myC1 ((Standard_Address)&C1), myC2 ((Standard_Address)&C2),
				        myTol (thetol)
  {
  math_Vector V1(1,2), V2(1,2);
  V1(1) = Tool1::FirstParameter(*((Curve1*)myC1));
  V2(1) = Tool1::LastParameter(*((Curve1*)myC1));
  V1(2) = Tool2::FirstParameter(*((Curve2*)myC2));
  V2(2) = Tool2::LastParameter(*((Curve2*)myC2));
  SubIntervalInitialize(V1, V2);
                        
  switch(Tool1::GetType(*((Curve1*)myC1)))
    {
    case GeomAbs_BezierCurve:
    case GeomAbs_BSplineCurve:
    case GeomAbs_OtherCurve:
      myMaxDerivOrderC1 = MaxOrder;
      myTolC1 = SearchOfTolerance((Standard_Address)&C1);
      break;
    default:
      myMaxDerivOrderC1 = 0;
      myTolC1=MinTol;
      break;
    }
    
  switch(Tool2::GetType(*((Curve2*)myC2)))
    {
    case GeomAbs_BezierCurve:
    case GeomAbs_BSplineCurve:
    case GeomAbs_OtherCurve:
      myMaxDerivOrderC2 = MaxOrder;
      myTolC2 = SearchOfTolerance((Standard_Address)&C2);
      break;
    default:
      myMaxDerivOrderC2 = 0;
      myTolC2=MinTol;
      break;
    }
  }

//=============================================================================
void Extrema_FuncExtCC::SetCurve (const Standard_Integer theRank, const Curve1& C)
  {
  Standard_OutOfRange_Raise_if (theRank < 1 || theRank > 2, "Extrema_FuncExtCC::SetCurve()")
  
  if (theRank == 1)
    {
    myC1 = (Standard_Address)&C;
    switch(/*Tool1::GetType(*((Curve1*)myC1))*/ C.GetType())
      {
      case GeomAbs_BezierCurve:
      case GeomAbs_BSplineCurve:
      case GeomAbs_OtherCurve:
        myMaxDerivOrderC1 = MaxOrder;
        myTolC1 = SearchOfTolerance((Standard_Address)&C);
        break;
      default:
        myMaxDerivOrderC1 = 0;
        myTolC1=MinTol;
        break;
      }
    }
  else if (theRank == 2)
    {
    myC2 = (Standard_Address)&C;
    switch(/*Tool2::GetType(*((Curve2*)myC2))*/C.GetType())
      {
      case GeomAbs_BezierCurve:
      case GeomAbs_BSplineCurve:
      case GeomAbs_OtherCurve:
        myMaxDerivOrderC2 = MaxOrder;
        myTolC2 = SearchOfTolerance((Standard_Address)&C);
        break;
      default:
        myMaxDerivOrderC2 = 0;
        myTolC2=MinTol;
        break;
      }
    }
  }

//=============================================================================
Standard_Boolean Extrema_FuncExtCC::Value (const math_Vector& UV, math_Vector& F)
  {
  myU = UV(1);
  myV = UV(2);
  Tool1::D1(*((Curve1*)myC1), myU,myP1,myDu);
  Tool2::D1(*((Curve2*)myC2), myV,myP2,myDv);

  Vec P1P2 (myP1,myP2);

  Standard_Real Ndu = myDu.Magnitude();


  if(myMaxDerivOrderC1 != 0)
    {
    if (Ndu <= myTolC1)
      {
    //Derivative is approximated by Taylor-series
      const Standard_Real DivisionFactor = 1.e-3;
      Standard_Real du;
      if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
        du = 0.0;
      else
        du = myUsupremum-myUinfium;
        
      const Standard_Real aDelta = Max(du*DivisionFactor,MinStep);
      
      Standard_Integer n = 1; //Derivative order
      Vec V;
      Standard_Boolean IsDeriveFound;
      
      do
        {
        V = Tool1::DN(*((Curve1*)myC1),myU,++n);
        Ndu = V.Magnitude();
        IsDeriveFound = (Ndu > myTolC1);
        }
      while(!IsDeriveFound && n < myMaxDerivOrderC1);
      
      if(IsDeriveFound)
        {
        Standard_Real u;
        
        if(myU-myUinfium < aDelta)
          u = myU+aDelta;
        else
          u = myU-aDelta;
          
        Pnt P1, P2;
        Tool1::D0(*((Curve1*)myC1),Min(myU, u),P1);
        Tool1::D0(*((Curve1*)myC1),Max(myU, u),P2);
        
        Vec V1(P1,P2);
        Standard_Real aDirFactor = V.Dot(V1);
        
        if(aDirFactor < 0.0)
          myDu = -V;
        else
          myDu = V;
        }//if(IsDeriveFound)
      else
        {
  //Derivative is approximated by three points

        Pnt Ptemp; //(0,0,0)-coordinate
        Pnt P1, P2, P3;
        Standard_Boolean IsParameterGrown;
                  
        if(myU-myUinfium < 2*aDelta)
          {
          Tool1::D0(*((Curve1*)myC1),myU,P1);
          Tool1::D0(*((Curve1*)myC1),myU+aDelta,P2);
          Tool1::D0(*((Curve1*)myC1),myU+2*aDelta,P3);
          IsParameterGrown = Standard_True;
          }
        else
          {
          Tool1::D0(*((Curve1*)myC1),myU-2*aDelta,P1);
          Tool1::D0(*((Curve1*)myC1),myU-aDelta,P2);
          Tool1::D0(*((Curve1*)myC1),myU,P3);
          IsParameterGrown = Standard_False;
          }
          
        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
        
        if(IsParameterGrown)
          myDu=-3*V1+4*V2-V3;
        else
          myDu=V1-4*V2+3*V3;
        }//else of if(IsDeriveFound)
      Ndu = myDu.Magnitude();      
      }//if (Ndu <= myTolC1) condition
    }//if(myMaxDerivOrder != 0)

  if (Ndu <= MinTol)
    {
#ifdef DEB
    cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
    cout << "Warning: 1st derivative of C1 is equal to zero!---"<<endl;
#endif
    return Standard_False;
    }


  Standard_Real Ndv = myDv.Magnitude();

  if(myMaxDerivOrderC2 != 0)
    {
    if (Ndv <= myTolC2)
      { 
      const Standard_Real DivisionFactor = 1.e-3;
      Standard_Real dv;
      if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
        dv = 0.0;
      else
        dv = myVsupremum-myVinfium;
        
      const Standard_Real aDelta = Max(dv*DivisionFactor,MinStep);

//Derivative is approximated by Taylor-series
      
      Standard_Integer n = 1; //Derivative order
      Vec V;
      Standard_Boolean IsDeriveFound;
      
      do
        {
        V = Tool2::DN(*((Curve2*)myC2),myV,++n);
        Ndv = V.Magnitude();
        IsDeriveFound = (Ndv > myTolC2);
        }
      while(!IsDeriveFound && n < myMaxDerivOrderC2);
      
      if(IsDeriveFound)
        {
        Standard_Real v;
        
        if(myV-myVinfium < aDelta)
          v = myV+aDelta;
        else
          v = myV-aDelta;
          
        Pnt P1, P2;
        Tool2::D0(*((Curve2*)myC2),Min(myV, v),P1);
        Tool2::D0(*((Curve2*)myC2),Max(myV, v),P2);
        
        Vec V1(P1,P2);
        Standard_Real aDirFactor = V.Dot(V1);
        
        if(aDirFactor < 0.0)
          myDv = -V;
        else
          myDv = V;
        }//if(IsDeriveFound)
      else
        {
  //Derivative is approximated by three points

        Pnt Ptemp; //(0,0,0)-coordinate
        Pnt P1, P2, P3;
        Standard_Boolean IsParameterGrown;
                  
        if(myV-myVinfium < 2*aDelta)
          {
          Tool2::D0(*((Curve2*)myC2),myV,P1);
          Tool2::D0(*((Curve2*)myC2),myV+aDelta,P2);
          Tool2::D0(*((Curve2*)myC2),myV+2*aDelta,P3);
          IsParameterGrown = Standard_True;
          }
        else
          {
          Tool2::D0(*((Curve2*)myC2),myV-2*aDelta,P1);
          Tool2::D0(*((Curve2*)myC2),myV-aDelta,P2);
          Tool2::D0(*((Curve2*)myC2),myV,P3);
          IsParameterGrown = Standard_False;
          }
          
        Vec V1(Ptemp,P1), V2(Ptemp,P2), V3(Ptemp,P3);
        
        if(IsParameterGrown)
          myDv=-3*V1+4*V2-V3;
        else
          myDv=V1-4*V2+3*V3;
        }//else of if(IsDeriveFound)
        
        Ndv = myDv.Magnitude();
      }//if (Ndv <= myTolC2)
    }//if(myMaxDerivOrder != 0)

  if (Ndv <= MinTol)
    {
#ifdef DEB
    cout << "+++Function Extrema_FuncExtCC::Value(...)." << endl;
    cout << "1st derivative of C2 is equal to zero!---"<<endl;
#endif
    return Standard_False;
    }

  F(1) = P1P2.Dot(myDu)/Ndu;
  F(2) = P1P2.Dot(myDv)/Ndv;
  return Standard_True;
  }
//=============================================================================

Standard_Boolean Extrema_FuncExtCC::Derivatives (const math_Vector& UV, 
						 math_Matrix& Df)
  {
  math_Vector F(1,2);
  return Values(UV,F,Df);
  }
//=============================================================================

Standard_Boolean Extrema_FuncExtCC::Values (const math_Vector& UV, 
                                            math_Vector& F,
                                            math_Matrix& Df)
  {
  myU = UV(1);
  myV = UV(2);

  if(Value(UV, F) == Standard_False) //Computes F, myDu, myDv
    {
#ifdef DEB
    cout << "+++Standard_Boolean Extrema_FuncExtCC::Values(...)." << endl;
    cout << "Warning: No function value found!---"<<endl;
#endif
    return Standard_False;
    }//if(Value(UV, F) == Standard_False)

  Vec Du, Dv, Duu, Dvv;
  Tool1::D2(*((Curve1*)myC1), myU,myP1,Du,Duu);
  Tool2::D2(*((Curve2*)myC2), myV,myP2,Dv,Dvv);
  
//Calling of "Value(...)" function change class member values. 
//After running it is necessary to return to previous values.  
  const Standard_Real myU_old  = myU,    myV_old = myV;
  const Pnt           myP1_old = myP1,  myP2_old = myP2;
  const Vec           myDu_old = myDu,  myDv_old = myDv;
  

//Attention:  aDelta value must be greater than same value for "Value(...)"
//            function to avoid of points' collisions.

  const Standard_Real DivisionFactor = 0.01;
  
  Standard_Real du;
  if((myUsupremum >= RealLast()) || (myUinfium <= RealFirst()))
    du = 0.0;
  else
    du = myUsupremum-myUinfium;
    
  const Standard_Real aDeltaU = Max(du*DivisionFactor,MinStep);
  
  Standard_Real dv;
  if((myVsupremum >= RealLast()) || (myVinfium <= RealFirst()))
    dv = 0.0;
  else
    dv = myVsupremum-myVinfium;
    
  const Standard_Real aDeltaV = Max(dv*DivisionFactor,MinStep); 

  Vec P1P2 (myP1,myP2);
  
  if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
    {
//Derivative is approximated by three points

    math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
    Standard_Real F1, F2, F3;

/////////////////////////// Search of DF1_u derivative (begin) ///////////////////
    if(myU-myUinfium < 2*aDeltaU)
      {
      F1=F(1);
      math_Vector UV2(1,2), UV3(1,2);
      UV2(1)=myU+aDeltaU;
      UV2(2)=myV;
      UV3(1)=myU+2*aDeltaU;
      UV3(2)=myV;
      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F2 = FF2(1);
      F3 = FF3(1);
      
      Df(1,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
      }//if(myU-myUinfium < 2*aDeltaU)
    else
      {
      F3 = F(1);
      math_Vector UV2(1,2), UV1(1,2);
      UV2(1)=myU-aDeltaU;
      UV2(2)=myV;
      UV1(1)=myU-2*aDeltaU;
      UV1(2)=myV;

      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F1 = FF1(1);
      F2 = FF2(1);

      Df(1,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
      }//else of if(myU-myUinfium < 2*aDeltaU) condition
/////////////////////////// Search of DF1_u derivative (end) ///////////////////
      
      //Return to previous values
      myU = myU_old;
      myV = myV_old;

/////////////////////////// Search of DF1_v derivative (begin) ///////////////////
    if(myV-myVinfium < 2*aDeltaV)
      {
      F1=F(1);
      math_Vector UV2(1,2), UV3(1,2);
      UV2(1)=myU;
      UV2(2)=myV+aDeltaV;
      UV3(1)=myU;
      UV3(2)=myV+2*aDeltaV;

      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
      F2 = FF2(1);
      F3 = FF3(1);
      
      Df(1,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
      }//if(myV-myVinfium < 2*aDeltaV)
    else
      {
      F3 = F(1);
      math_Vector UV2(1,2), UV1(1,2);
      UV2(1)=myU;
      UV2(2)=myV-aDeltaV;
      UV1(1)=myU;
      UV1(2)=myV-2*aDeltaV;
      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F1 = FF1(1);
      F2 = FF2(1);

      Df(1,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
      }//else of if(myV-myVinfium < 2*aDeltaV)
/////////////////////////// Search of DF1_v derivative (end) ///////////////////
    
    //Return to previous values
    myU = myU_old;
    myV = myV_old;
    myP1 = myP1_old,  myP2 = myP2_old;
    myDu = myDu_old,  myDv = myDv_old;
    }//if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
  else
    {
    const Standard_Real Ndu = myDu.Magnitude();
    Df(1,1) = - Ndu + (P1P2.Dot(Duu)/Ndu) - F(1)*(myDu.Dot(Duu)/(Ndu*Ndu));
    Df(1,2) = myDv.Dot(myDu)/Ndu;
    }//else of if((myMaxDerivOrderC1 != 0) && (Du.Magnitude() <= myTolC1))
  
  if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
    {
//Derivative is approximated by three points

    math_Vector FF1(1,2), FF2(1,2), FF3(1,2);
    Standard_Real F1, F2, F3;

/////////////////////////// Search of DF2_v derivative (begin) ///////////////////
    if(myV-myVinfium < 2*aDeltaV)
      {
      F1=F(2);
      math_Vector UV2(1,2), UV3(1,2);
      UV2(1)=myU;
      UV2(2)=myV+aDeltaV;
      UV3(1)=myU;
      UV3(2)=myV+2*aDeltaV;

      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F2 = FF2(2);
      F3 = FF3(2);
      
      Df(2,2) = (-3*F1+4*F2-F3)/(2.0*aDeltaV);
      
      }//if(myV-myVinfium < 2*aDeltaV)
    else
      {
      F3 = F(2);
      math_Vector UV2(1,2), UV1(1,2);
      UV2(1)=myU;
      UV2(2)=myV-aDeltaV;
      UV1(1)=myU;
      UV1(2)=myV-2*aDeltaV;

      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F1 = FF1(2);
      F2 = FF2(2);

      Df(2,2) = (F1-4*F2+3*F3)/(2.0*aDeltaV);
      }//else of if(myV-myVinfium < 2*aDeltaV)
/////////////////////////// Search of DF2_v derivative (end) ///////////////////

      //Return to previous values
      myU = myU_old;
      myV = myV_old;

/////////////////////////// Search of DF2_u derivative (begin) ///////////////////
    if(myU-myUinfium < 2*aDeltaU)
      {
      F1=F(2);
      math_Vector UV2(1,2), UV3(1,2);
      UV2(1)=myU+aDeltaU;
      UV2(2)=myV;
      UV3(1)=myU+2*aDeltaU;
      UV3(2)=myV;
      if(!((Value(UV2,FF2)) && (Value(UV3,FF3))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F2 = FF2(2);
      F3 = FF3(2);
      
      Df(2,1) = (-3*F1+4*F2-F3)/(2.0*aDeltaU);
      }//if(myU-myUinfium < 2*aDelta)
    else
      {
      F3 = F(2);
      math_Vector UV2(1,2), UV1(1,2);
      UV2(1)=myU-aDeltaU;
      UV2(2)=myV;
      UV1(1)=myU-2*aDeltaU;
      UV1(2)=myV;
            
      if(!((Value(UV2,FF2)) && (Value(UV1,FF1))))
        {
#ifdef DEB
        cout << "+++ Function Extrema_FuncExtCC::Values(...)" << endl;
        cout << "There are many points close to singularity points "
            "and which have zero-derivative." << endl;
        cout << "Try to decrease aDelta variable's value. ---" << endl;
#endif
        return Standard_False;
        }
        
      F1 = FF1(2);
      F2 = FF2(2);

      Df(2,1) = (F1-4*F2+3*F3)/(2.0*aDeltaU);
      }//else of if(myU-myUinfium < 2*aDeltaU)
/////////////////////////// Search of DF2_u derivative (end) ///////////////////
    
    //Return to previous values
    myU = myU_old;
    myV = myV_old;
    myP1 = myP1_old;
    myP2 = myP2_old;
    myDu = myDu_old;
    myDv = myDv_old;
    }//if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
  else
    {
    Standard_Real Ndv = myDv.Magnitude();
    Df(2,2) = Ndv + (P1P2.Dot(Dvv)/Ndv) - F(2)*(myDv.Dot(Dvv)/(Ndv*Ndv));
    Df(2,1) = -myDu.Dot(myDv)/Ndv;    
    }//else of if((myMaxDerivOrderC2 != 0) && (Dv.Magnitude() <= myTolC2))
  
  return Standard_True;

  }//end of function
//=============================================================================

Standard_Integer Extrema_FuncExtCC::GetStateNumber ()
{
  Vec Du (myDu), Dv (myDv);
  Vec P1P2 (myP1, myP2);

  Standard_Real mod = Du.Magnitude();
  if(mod > myTolC1) {
    Du /= mod;
  }
  mod = Dv.Magnitude();
  if(mod > myTolC2) {
    Dv /= mod;
  }

  if (Abs(P1P2.Dot(Du)) <= myTol && Abs(P1P2.Dot(Dv)) <= myTol) {
    mySqDist.Append(myP1.SquareDistance(myP2));
    myPoints.Append(POnC(myU, myP1));
    myPoints.Append(POnC(myV, myP2));
  }
  return 0;
}
//=============================================================================

void Extrema_FuncExtCC::Points (const Standard_Integer N,
				POnC& P1,
				POnC& P2) const
{
  P1 = myPoints.Value(2*N-1);
  P2 = myPoints.Value(2*N);
}
//=============================================================================

void Extrema_FuncExtCC::SubIntervalInitialize(const math_Vector& theInfBound,
                                              const math_Vector& theSupBound)
  {
  myUinfium = theInfBound(1);
  myUsupremum = theSupBound(1);
  myVinfium = theInfBound(2);
  myVsupremum = theSupBound(2);
  }
//=============================================================================