diff --git a/src/QABugs/QABugs_19.cxx b/src/QABugs/QABugs_19.cxx
index 2744491ca8..7d30c41ece 100755
--- a/src/QABugs/QABugs_19.cxx
+++ b/src/QABugs/QABugs_19.cxx
@@ -1234,6 +1234,69 @@ static Standard_Integer OCC24005 (Draw_Interpretor& theDI, Standard_Integer theN
   return 0;
 }
 
+#include <Extrema_FuncExtPS.hxx>
+#include <math_FunctionSetRoot.hxx>
+#include <math_Vector.hxx>
+#include <BRepBuilderAPI_MakeVertex.hxx>
+static Standard_Integer OCC24137 (Draw_Interpretor& theDI, Standard_Integer theNArg, const char** theArgv) 
+{
+  Standard_Integer anArgIter = 1;
+  if (theNArg < 5)
+    {
+      theDI <<"Usage: " << theArgv[0] << " face vertex U V [N]"<<"\n";
+      return 1;
+    }
+
+  // get target shape
+  Standard_CString aFaceName = theArgv[anArgIter++];
+  Standard_CString aVertName = theArgv[anArgIter++];
+  const TopoDS_Shape aShapeF = DBRep::Get (aFaceName);
+  const TopoDS_Shape aShapeV = DBRep::Get (aVertName);
+  const Standard_Real aUFrom = Atof (theArgv[anArgIter++]);
+  const Standard_Real aVFrom = Atof (theArgv[anArgIter++]);
+  const Standard_Integer aNbIts = (anArgIter < theNArg) ? atol (theArgv[anArgIter++]) : 100;
+  if (aShapeF.IsNull() || aShapeF.ShapeType() != TopAbs_FACE)
+    {
+      std::cout << "Error: " << aFaceName << " shape is null / not a face" << std::endl;
+      return 1;
+    }
+  if (aShapeV.IsNull() || aShapeV.ShapeType() != TopAbs_VERTEX)
+    {
+      std::cout << "Error: " << aVertName << " shape is null / not a vertex" << std::endl;
+      return 1;
+    }
+  const TopoDS_Face   aFace = TopoDS::Face   (aShapeF);
+  const TopoDS_Vertex aVert = TopoDS::Vertex (aShapeV);
+  GeomAdaptor_Surface aSurf (BRep_Tool::Surface (aFace));
+
+  gp_Pnt aPnt = BRep_Tool::Pnt (aVert), aRes;
+
+  Extrema_FuncExtPS    anExtFunc;
+  math_FunctionSetRoot aRoot (anExtFunc, aNbIts);
+
+  math_Vector aTolUV (1, 2), aUVinf  (1, 2), aUVsup  (1, 2), aFromUV (1, 2);
+  aTolUV (1) =  Precision::Confusion(); aTolUV (2) =  Precision::Confusion();
+  aUVinf (1) = -Precision::Infinite();  aUVinf (2) = -Precision::Infinite();
+  aUVsup (1) =  Precision::Infinite();  aUVsup (2) =  Precision::Infinite();
+  aFromUV(1) =  aUFrom; aFromUV(2) = aVFrom;
+
+  anExtFunc.Initialize (aSurf);
+  anExtFunc.SetPoint (aPnt);
+  aRoot.SetTolerance (aTolUV);
+  aRoot.Perform (anExtFunc, aFromUV, aUVinf, aUVsup);
+  if (!aRoot.IsDone())
+    {
+      std::cerr << "No results!\n";
+      return 1;
+    }
+
+  theDI << aRoot.Root()(1) << " " << aRoot.Root()(2) << "\n";
+  
+  aSurf.D0 (aRoot.Root()(1), aRoot.Root()(2), aRes);
+  DBRep::Set ("result", BRepBuilderAPI_MakeVertex (aRes));
+  return 0;
+}
+
 void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   const char *group = "QABugs";
 
@@ -1254,5 +1317,6 @@ void QABugs::Commands_19(Draw_Interpretor& theCommands) {
   theCommands.Add ("OCC24019", "OCC24019 aShape", __FILE__, OCC24019, group);
   theCommands.Add ("OCC11758", "OCC11758", __FILE__, OCC11758, group);
   theCommands.Add ("OCC24005", "OCC24005 result", __FILE__, OCC24005, group);
+  theCommands.Add ("OCC24137", "OCC24137 face vertex U V [N]", __FILE__, OCC24137, group);
   return;
 }
diff --git a/src/math/math_FunctionSetRoot.cxx b/src/math/math_FunctionSetRoot.cxx
index 134ac626a4..e06513dbbd 100755
--- a/src/math/math_FunctionSetRoot.cxx
+++ b/src/math/math_FunctionSetRoot.cxx
@@ -65,56 +65,56 @@
 
 #define FSR_DEBUG(arg)
 // Uncomment the following code to have debug output to cout 
-/* 
-static Standard_Boolean mydebug = Standard_False;
-#undef FSR_DEBUG
-#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
-*/
+//========================================================== 
+//static Standard_Boolean mydebug = Standard_True;
+//#undef FSR_DEBUG
+//#define FSR_DEBUG(arg) {if (mydebug) { cout << arg << endl; }}
+//===========================================================
 
 class MyDirFunction : public math_Function
 {
 
-     math_Vector *P0;
-     math_Vector *Dir;
-     math_Vector *P;
-     math_Vector *FV;
-     math_FunctionSetWithDerivatives *F;
+  math_Vector *P0;
+  math_Vector *Dir;
+  math_Vector *P;
+  math_Vector *FV;
+  math_FunctionSetWithDerivatives *F;
 
 public :
 
-     MyDirFunction(math_Vector& V1, 
-		   math_Vector& V2,
-		   math_Vector& V3,
-		   math_Vector& V4,
-		   math_FunctionSetWithDerivatives& f) ;
-     
-     void Initialize(const math_Vector& p0, const math_Vector& dir) const;
-//For hp :
-     Standard_Boolean Value(const math_Vector& Sol, math_Vector& FF,
-			    math_Matrix& DF, math_Vector& GH, 
-			    Standard_Real& F2, Standard_Real& Gnr1);
-//     Standard_Boolean MyDirFunction::Value(const math_Vector& Sol, math_Vector& FF,
-//					   math_Matrix& DF, math_Vector& GH, 
-//					   Standard_Real& F2, Standard_Real& Gnr1);
-     Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ;
-     
+  MyDirFunction(math_Vector& V1, 
+    math_Vector& V2,
+    math_Vector& V3,
+    math_Vector& V4,
+    math_FunctionSetWithDerivatives& f) ;
+
+  void Initialize(const math_Vector& p0, const math_Vector& dir) const;
+  //For hp :
+  Standard_Boolean Value(const math_Vector& Sol, math_Vector& FF,
+    math_Matrix& DF, math_Vector& GH, 
+    Standard_Real& F2, Standard_Real& Gnr1);
+  //     Standard_Boolean MyDirFunction::Value(const math_Vector& Sol, math_Vector& FF,
+  //					   math_Matrix& DF, math_Vector& GH, 
+  //					   Standard_Real& F2, Standard_Real& Gnr1);
+  Standard_Boolean Value(const Standard_Real x, Standard_Real& fval) ;
+
 };
 
 MyDirFunction::MyDirFunction(math_Vector& V1, 
-			     math_Vector& V2,
-			     math_Vector& V3,
-			     math_Vector& V4,
-			     math_FunctionSetWithDerivatives& f) {
-        
-  P0  = &V1;
-  Dir = &V2;
-  P   = &V3;
-  FV =  &V4;
-  F   = &f;
+                             math_Vector& V2,
+                             math_Vector& V3,
+                             math_Vector& V4,
+                             math_FunctionSetWithDerivatives& f) {
+
+                               P0  = &V1;
+                               Dir = &V2;
+                               P   = &V3;
+                               FV =  &V4;
+                               F   = &f;
 }
 
 void MyDirFunction::Initialize(const math_Vector& p0, 
-				    const math_Vector& dir)  const
+                               const math_Vector& dir)  const
 {
   *P0 = p0;
   *Dir = dir;
@@ -122,7 +122,7 @@ void MyDirFunction::Initialize(const math_Vector& p0,
 
 
 Standard_Boolean MyDirFunction::Value(const Standard_Real x, 
-				      Standard_Real& fval) 
+                                      Standard_Real& fval) 
 {
   Standard_Real p;
   for(Standard_Integer i = P->Lower(); i <= P->Upper(); i++) {
@@ -136,13 +136,13 @@ Standard_Boolean MyDirFunction::Value(const Standard_Real x,
     for(Standard_Integer i = FV->Lower(); i <= FV->Upper(); i++) {
       aVal = FV->Value(i);
       if(aVal < 0.) {
-	if(aVal <= -1.e+100) // Precision::HalfInfinite() later
-//       if(Precision::IsInfinite(Abs(FV->Value(i)))) {
-//	fval = Precision::Infinite();
-	return Standard_False;
+        if(aVal <= -1.e+100) // Precision::HalfInfinite() later
+          //       if(Precision::IsInfinite(Abs(FV->Value(i)))) {
+          //	fval = Precision::Infinite();
+          return Standard_False;
       }
       else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
-	return Standard_False;
+        return Standard_False;
     }
 
     fval = 0.5 * (FV->Norm2());
@@ -152,29 +152,29 @@ Standard_Boolean MyDirFunction::Value(const Standard_Real x,
 }
 
 Standard_Boolean  MyDirFunction::Value(const math_Vector& Sol,
-				       math_Vector& FF,
-				       math_Matrix& DF,
-				       math_Vector& GH,
-				       Standard_Real& F2,
-				       Standard_Real& Gnr1)
+                                       math_Vector& FF,
+                                       math_Matrix& DF,
+                                       math_Vector& GH,
+                                       Standard_Real& F2,
+                                       Standard_Real& Gnr1)
 {
   if( F->Values(Sol, FF, DF) ) {
 
     Standard_Real aVal = 0.;
 
     for(Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
-// modified by NIZHNY-MKK  Mon Oct  3 17:56:50 2005.BEGIN
+      // modified by NIZHNY-MKK  Mon Oct  3 17:56:50 2005.BEGIN
       aVal = FF.Value(i);
       if(aVal < 0.) {
-	if(aVal <= -1.e+100) // Precision::HalfInfinite() later
-//       if(Precision::IsInfinite(Abs(FF.Value(i)))) {
-//	F2 = Precision::Infinite();
-//	Gnr1 = Precision::Infinite();
-	  return Standard_False;
+        if(aVal <= -1.e+100) // Precision::HalfInfinite() later
+          //       if(Precision::IsInfinite(Abs(FF.Value(i)))) {
+          //	F2 = Precision::Infinite();
+          //	Gnr1 = Precision::Infinite();
+          return Standard_False;
       }
       else if(aVal >= 1.e+100) // Precision::HalfInfinite() later
-	return Standard_False;
-// modified by NIZHNY-MKK  Mon Oct  3 17:57:05 2005.END
+        return Standard_False;
+      // modified by NIZHNY-MKK  Mon Oct  3 17:57:05 2005.END
     }
 
 
@@ -196,14 +196,14 @@ Standard_Boolean  MyDirFunction::Value(const math_Vector& Sol,
 
 //--------------------------------------------------------------
 static Standard_Boolean MinimizeDirection(const math_Vector&   P0,
-					  const math_Vector&   P1,
-					  const math_Vector&   P2,
-					  const Standard_Real  F1,
-					  math_Vector&         Delta,
-					  const math_Vector&   Tol,
-					  MyDirFunction& F)
-// Purpose : minimisation a partir de 3 points
-//-------------------------------------------------------
+                                          const math_Vector&   P1,
+                                          const math_Vector&   P2,
+                                          const Standard_Real  F1,
+                                          math_Vector&         Delta,
+                                          const math_Vector&   Tol,
+                                          MyDirFunction& F)
+                                          // Purpose : minimisation a partir de 3 points
+                                          //-------------------------------------------------------
 {
   // (1) Evaluation d'un tolerance parametrique 1D
   Standard_Real tol1d = 2.1 , invnorme, tsol;
@@ -217,11 +217,11 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P0,
   if (tol1d > 1.9) return Standard_False; //Pas la peine de se fatiguer
   tol1d /= 3; 
 
-//JR/Hp :
+  //JR/Hp :
   math_Vector PP0 = P0 ;
   math_Vector PP1 = P1 ;
   Delta = PP1 - PP0;
-//  Delta = P1 - P0;
+  //  Delta = P1 - P0;
   invnorme = Delta.Norm();
   if (invnorme <= Eps) return Standard_False;
   invnorme = ((Standard_Real) 1) / invnorme;
@@ -230,7 +230,7 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P0,
 
   // (2) On minimise
   FSR_DEBUG ("      minimisation dans la direction")
-  ax = -1; bx = 0;
+    ax = -1; bx = 0;
   cx = (P2-P1).Norm()*invnorme;
   if (cx < 1.e-2) return Standard_False;
   math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
@@ -246,15 +246,15 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P0,
 
 //----------------------------------------------------------------------
 static Standard_Boolean MinimizeDirection(const math_Vector&   P,
-					  math_Vector&   Dir,
-					  const Standard_Real& PValue,
-					  const Standard_Real& PDirValue,
-					  const math_Vector&   Gradient,
-					  const math_Vector&   DGradient,
-					  const math_Vector&   Tol,
-					  MyDirFunction& F)
-// Purpose: minimisation a partir de 2 points et une derives
-//----------------------------------------------------------------------
+                                          math_Vector&   Dir,
+                                          const Standard_Real& PValue,
+                                          const Standard_Real& PDirValue,
+                                          const math_Vector&   Gradient,
+                                          const math_Vector&   DGradient,
+                                          const math_Vector&   Tol,
+                                          MyDirFunction& F)
+                                          // Purpose: minimisation a partir de 2 points et une derives
+                                          //----------------------------------------------------------------------
 
 {
   // (0) Evaluation d'un tolerance parametrique 1D
@@ -267,10 +267,10 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
     if (absdir >Eps) tol1d = Min(tol1d, Tol(ii) / absdir);
   }
   if (tol1d > 0.9) return Standard_False;
-    
+
   // (1) On realise une premiere interpolation quadratique
   Standard_Real ax, bx, cx, df1, df2, Delta, tsol, fsol, tsolbis;
-  FSR_DEBUG("     essai d interpolation")
+  FSR_DEBUG("     essai d interpolation");
 
   df1 = Gradient*Dir;
   df2 = DGradient*Dir;
@@ -290,16 +290,16 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
     else { // cas quadratique
       Delta = bx*bx - 4*ax*cx;
       if (Delta > 1.e-9) {
-	// il y a des racines, on prend la plus proche de 0
-	Delta = Sqrt(Delta);
-	tsol = -(bx + Delta);
-	tsolbis = (Delta - bx);
-	if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
-	tsol /= 2*ax;
+        // il y a des racines, on prend la plus proche de 0
+        Delta = Sqrt(Delta);
+        tsol = -(bx + Delta);
+        tsolbis = (Delta - bx);
+        if (Abs(tsolbis) < Abs(tsol)) tsol = tsolbis;
+        tsol /= 2*ax;
       }
       else {
-	// pas ou peu de racine : on "extremise"
-	tsol = -(0.5*bx)/ax;
+        // pas ou peu de racine : on "extremise"
+        tsol = -(0.5*bx)/ax;
       }
     }
   }
@@ -312,26 +312,26 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
   if (fsol<PValue) { 
     good = Standard_True;
     Result = fsol;
-    FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue)
+    FSR_DEBUG("t= "<<tsol<<" F = " << fsol << " OldF = "<<PValue);
   }
 
   // (2) Si l'on a pas assez progresser on realise une recherche 
   //     en bonne et due forme, a partir des inits precedents
   if ((fsol > 0.2*PValue) && (tol1d < 0.5)) {
-    
+
     if (tsol <0) {
       ax = tsol; bx = 0.0; cx = 1.0;
     }
     else {
       ax = 0.0; bx = tsol; cx = 1.0;
     }
-    FSR_DEBUG(" minimisation dans la direction")
+    FSR_DEBUG(" minimisation dans la direction");
     math_BrentMinimum Sol(F, ax, bx, cx, tol1d, 100, tol1d);
     if(Sol.IsDone()) {
       if (Sol.Minimum() <= Result) {
         tsol = Sol.Location();
         good = Standard_True;
-	FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
+        FSR_DEBUG("t= "<<tsol<<" F ="<< Sol.Minimum() << " OldF = "<<Result)
       }
     }
   }
@@ -344,12 +344,12 @@ static Standard_Boolean MinimizeDirection(const math_Vector&   P,
 
 //------------------------------------------------------
 static void SearchDirection(const math_Matrix& DF,
-			    const math_Vector& GH,
-			    const math_Vector& FF,
-			    Standard_Boolean ChangeDirection,
-			    const math_Vector& InvLengthMax, 
-			    math_Vector& Direction,
-			    Standard_Real& Dy)
+                            const math_Vector& GH,
+                            const math_Vector& FF,
+                            Standard_Boolean ChangeDirection,
+                            const math_Vector& InvLengthMax, 
+                            math_Vector& Direction,
+                            Standard_Real& Dy)
 
 {
   Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
@@ -357,12 +357,12 @@ static void SearchDirection(const math_Matrix& DF,
   if (!ChangeDirection) {
     if (Ninc == Neq) {
       for (Standard_Integer i = FF.Lower(); i <= FF.Upper(); i++) {
-	Direction(i) = -FF(i);
+        Direction(i) = -FF(i);
       }
       math_Gauss Solut(DF, 1.e-9);
       if (Solut.IsDone()) Solut.Solve(Direction);
       else { // we have to "forget" singular directions.
-	FSR_DEBUG(" Matrice singuliere : On prend SVD")
+        FSR_DEBUG(" Matrice singuliere : On prend SVD");
         math_SVD SolvebySVD(DF);
         if (SolvebySVD.IsDone()) SolvebySVD.Solve(-1*FF, Direction);
         else ChangeDirection = Standard_True;
@@ -383,9 +383,9 @@ static void SearchDirection(const math_Matrix& DF,
   // Afin de blinder les cas trop mal conditionne
   // PMN 12/05/97 Traitement des singularite dans les conges
   // Sur des surfaces periodiques
-  
+
   Standard_Real ratio = Abs( Direction(Direction.Lower())
-			    *InvLengthMax(Direction.Lower()) );
+    *InvLengthMax(Direction.Lower()) );
   Standard_Integer i;
   for (i = Direction.Lower()+1; i<=Direction.Upper(); i++) {
     ratio = Max(ratio,  Abs( Direction(i)*InvLengthMax(i)) );
@@ -393,7 +393,7 @@ static void SearchDirection(const math_Matrix& DF,
   if (ratio > 1) {
     Direction /= ratio;
   }
-  
+
   Dy = Direction*GH;
   if (Dy >= -Eps) { // newton "ne descend pas" on prend le gradient
     ChangeDirection = Standard_True;
@@ -409,18 +409,18 @@ static void SearchDirection(const math_Matrix& DF,
 
 //=====================================================================
 static void SearchDirection(const math_Matrix& DF, 
-			    const math_Vector& GH,
-			    const math_Vector& FF,
-			    const math_IntegerVector& Constraints,
-//			    const math_Vector& X, // Le point d'init
-			    const math_Vector& , // Le point d'init
-			    Standard_Boolean ChangeDirection,
+                            const math_Vector& GH,
+                            const math_Vector& FF,
+                            const math_IntegerVector& Constraints,
+                            //			    const math_Vector& X, // Le point d'init
+                            const math_Vector& , // Le point d'init
+                            Standard_Boolean ChangeDirection,
                             const math_Vector& InvLengthMax,
-			    math_Vector& Direction,
-			    Standard_Real& Dy)
-//Purpose : Recherche une direction (et un pas si Newton Fonctionne) le long
-//          d'une frontiere
-//=====================================================================
+                            math_Vector& Direction,
+                            Standard_Real& Dy)
+                            //Purpose : Recherche une direction (et un pas si Newton Fonctionne) le long
+                            //          d'une frontiere
+                            //=====================================================================
 {
   Standard_Integer Ninc = DF.ColNumber(), Neq = DF.RowNumber();
   Standard_Integer i, j, k, Cons = 0;
@@ -434,12 +434,12 @@ static void SearchDirection(const math_Matrix& DF,
 
   if (Cons == 0) {
     SearchDirection(DF, GH, FF, ChangeDirection, InvLengthMax, 
-		    Direction, Dy);
+      Direction, Dy);
   }
   else if (Cons == Ninc) { // il n'y a plus rien a faire...
     for(Standard_Integer i = Direction.Lower(); i <= Direction.Upper(); i++) {
-	Direction(i) = 0;
-      }
+      Direction(i) = 0;
+    }
     Dy = 0;
   }
   else { //(1) Cas general : On definit un sous probleme
@@ -450,31 +450,31 @@ static void SearchDirection(const math_Matrix& DF,
 
     for (k=1, i = 1; i <= Ninc; i++) {
       if (Constraints(i) == 0) { 
-	  MyGH(k) = GH(i);
-          MyInvLengthMax(k) = InvLengthMax(i);
-	  MyDirection(k) = Direction(i);
-	  for (j = 1; j <= Neq; j++) {
-	   DF2(j, k) = DF(j, i);
-	 }
-	  k++; //on passe a l'inconnue suivante
-	}
+        MyGH(k) = GH(i);
+        MyInvLengthMax(k) = InvLengthMax(i);
+        MyDirection(k) = Direction(i);
+        for (j = 1; j <= Neq; j++) {
+          DF2(j, k) = DF(j, i);
+        }
+        k++; //on passe a l'inconnue suivante
+      }
     }
     //(2) On le resoud
     SearchDirection(DF2, MyGH, FF, ChangeDirection, MyInvLengthMax, 
-		    MyDirection, Dy);
+      MyDirection, Dy);
 
     // (3) On l'interprete...
     // Reconstruction de Direction:
     for (i = 1, k = 1; i <= Ninc; i++) {
       if (Constraints(i) == 0) {
-	if (!ChangeDirection) {
-	  Direction(i) = MyDirection(k);
-	}
-	else Direction(i) = - GH(i);
-	k++;
+        if (!ChangeDirection) {
+          Direction(i) = MyDirection(k);
+        }
+        else Direction(i) = - GH(i);
+        k++;
       }
       else {
-	Direction(i) = 0.0;
+        Direction(i) = 0.0;
       }
     }
   }
@@ -491,18 +491,18 @@ Standard_Boolean Bounds(const math_Vector&  InfBound,
                         math_IntegerVector& Constraints,
                         math_Vector&        Delta,
                         Standard_Boolean&   theIsNewSol)
-//
-// Purpose: Trims an initial solution Sol to be within a domain defined by
-//   InfBound and SupBound. Delta will contain a distance between final Sol and
-//   SolSave.
-//   IsNewSol returns False, if final Sol fully coincides with SolSave, i.e.
-//   if SolSave already lied on a boundary and initial Sol was fully beyond it
-//======================================================
+                        //
+                        // Purpose: Trims an initial solution Sol to be within a domain defined by
+                        //   InfBound and SupBound. Delta will contain a distance between final Sol and
+                        //   SolSave.
+                        //   IsNewSol returns False, if final Sol fully coincides with SolSave, i.e.
+                        //   if SolSave already lied on a boundary and initial Sol was fully beyond it
+                        //======================================================
 {
   Standard_Boolean Out = Standard_False;
   Standard_Integer i, Ninc = Sol.Length();
   Standard_Real    monratio = 1;
-  
+
   theIsNewSol = Standard_True;
 
   // Calcul du ratio de recadrage
@@ -557,29 +557,29 @@ Standard_Boolean Bounds(const math_Vector&  InfBound,
 
 
 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
-					   const math_Vector& Tolerance,
-					   const Standard_Integer NbIterations) :
-                                           Delta(1, F.NbVariables()),
-                                           Sol(1, F.NbVariables()),
-                                           DF(1, F.NbEquations(),
-					      1, F.NbVariables()),
-                                           Tol(1,F.NbVariables()),
+                                           const math_Vector& Tolerance,
+                                           const Standard_Integer NbIterations) :
+Delta(1, F.NbVariables()),
+Sol(1, F.NbVariables()),
+DF(1, F.NbEquations(),
+   1, F.NbVariables()),
+   Tol(1,F.NbVariables()),
 
-                                           InfBound(1, F.NbVariables()),
-                                           SupBound(1, F.NbVariables()),
+   InfBound(1, F.NbVariables()),
+   SupBound(1, F.NbVariables()),
 
-                                           SolSave(1, F.NbVariables()), 
-					   GH(1, F.NbVariables()), 
-					   DH(1, F.NbVariables()), 
-					   DHSave(1, F.NbVariables()),
-					   FF(1, F.NbEquations()), 
-					   PreviousSolution(1, F.NbVariables()), 
-					   Save(0, NbIterations),
-					   Constraints(1, F.NbVariables()),
-					   Temp1(1, F.NbVariables()), 
-					   Temp2(1, F.NbVariables()), 
-					   Temp3(1, F.NbVariables()), 
-					   Temp4(1, F.NbEquations())
+   SolSave(1, F.NbVariables()), 
+   GH(1, F.NbVariables()), 
+   DH(1, F.NbVariables()), 
+   DHSave(1, F.NbVariables()),
+   FF(1, F.NbEquations()), 
+   PreviousSolution(1, F.NbVariables()), 
+   Save(0, NbIterations),
+   Constraints(1, F.NbVariables()),
+   Temp1(1, F.NbVariables()), 
+   Temp2(1, F.NbVariables()), 
+   Temp3(1, F.NbVariables()), 
+   Temp4(1, F.NbEquations())
 
 {
   for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
@@ -589,28 +589,28 @@ math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
 }
 
 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
-					   const Standard_Integer NbIterations) :
-                                           Delta(1, F.NbVariables()),
-                                           Sol(1, F.NbVariables()),
-                                           DF(1, F.NbEquations(),
-					      1, F.NbVariables()),
-                                           Tol(1, F.NbVariables()),
+                                           const Standard_Integer NbIterations) :
+Delta(1, F.NbVariables()),
+Sol(1, F.NbVariables()),
+DF(1, F.NbEquations(),
+   1, F.NbVariables()),
+   Tol(1, F.NbVariables()),
 
-                                           InfBound(1, F.NbVariables()),
-                                           SupBound(1, F.NbVariables()),
+   InfBound(1, F.NbVariables()),
+   SupBound(1, F.NbVariables()),
 
-                                           SolSave(1, F.NbVariables()), 
-					   GH(1, F.NbVariables()), 
-					   DH(1, F.NbVariables()), 
-					   DHSave(1, F.NbVariables()),
-					   FF(1, F.NbEquations()), 
-					   PreviousSolution(1, F.NbVariables()), 
-					   Save(0, NbIterations),
-					   Constraints(1, F.NbVariables()),
-					   Temp1(1, F.NbVariables()), 
-					   Temp2(1, F.NbVariables()), 
-					   Temp3(1, F.NbVariables()), 
-					   Temp4(1, F.NbEquations())
+   SolSave(1, F.NbVariables()), 
+   GH(1, F.NbVariables()), 
+   DH(1, F.NbVariables()), 
+   DHSave(1, F.NbVariables()),
+   FF(1, F.NbEquations()), 
+   PreviousSolution(1, F.NbVariables()), 
+   Save(0, NbIterations),
+   Constraints(1, F.NbVariables()),
+   Temp1(1, F.NbVariables()), 
+   Temp2(1, F.NbVariables()), 
+   Temp3(1, F.NbVariables()), 
+   Temp4(1, F.NbEquations())
 
 {
   Itermax = NbIterations;
@@ -619,37 +619,37 @@ math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
 
 
 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
-					   const math_Vector& StartingPoint,
-					   const math_Vector& Tolerance,
-					   const math_Vector& infBound,
-					   const math_Vector& supBound,
-					   const Standard_Integer NbIterations,
-					   Standard_Boolean theStopOnDivergent) :
-                                           Delta(1, F.NbVariables()),
-                                           Sol(1, F.NbVariables()),
-                                           DF(1, F.NbEquations(),
-					      1, F.NbVariables()),
-                                           Tol(1,F.NbVariables()),
+                                           const math_Vector& StartingPoint,
+                                           const math_Vector& Tolerance,
+                                           const math_Vector& infBound,
+                                           const math_Vector& supBound,
+                                           const Standard_Integer NbIterations,
+                                           Standard_Boolean theStopOnDivergent) :
+Delta(1, F.NbVariables()),
+Sol(1, F.NbVariables()),
+DF(1, F.NbEquations(),
+   1, F.NbVariables()),
+   Tol(1,F.NbVariables()),
 
 
-                                           InfBound(1, F.NbVariables()),
-                                           SupBound(1, F.NbVariables()),
+   InfBound(1, F.NbVariables()),
+   SupBound(1, F.NbVariables()),
 
-                                           SolSave(1, F.NbVariables()), 
-					   GH(1, F.NbVariables()), 
-					   DH(1, F.NbVariables()), 
-					   DHSave(1, F.NbVariables()),
-					   FF(1, F.NbEquations()), 
-					   PreviousSolution(1, F.NbVariables()), 
-					   Save(0, NbIterations),
-					   Constraints(1, F.NbVariables()),
-					   Temp1(1, F.NbVariables()), 
-					   Temp2(1, F.NbVariables()), 
-					   Temp3(1, F.NbVariables()), 
-					   Temp4(1, F.NbEquations())
+   SolSave(1, F.NbVariables()), 
+   GH(1, F.NbVariables()), 
+   DH(1, F.NbVariables()), 
+   DHSave(1, F.NbVariables()),
+   FF(1, F.NbEquations()), 
+   PreviousSolution(1, F.NbVariables()), 
+   Save(0, NbIterations),
+   Constraints(1, F.NbVariables()),
+   Temp1(1, F.NbVariables()), 
+   Temp2(1, F.NbVariables()), 
+   Temp3(1, F.NbVariables()), 
+   Temp4(1, F.NbEquations())
 
 {
-					  
+
   for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
     Tol(i) =Tolerance(i);
   }
@@ -658,32 +658,32 @@ math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
 }
 
 math_FunctionSetRoot::math_FunctionSetRoot(math_FunctionSetWithDerivatives& F,
-					   const math_Vector& StartingPoint,
-					   const math_Vector& Tolerance,
-					   const Standard_Integer NbIterations) :
-                                           Delta(1, F.NbVariables()),
-                                           Sol(1, F.NbVariables()),
-                                           DF(1, F.NbEquations(),
-					      1, StartingPoint.Length()),
-                                           Tol(1,F.NbVariables()),
+                                           const math_Vector& StartingPoint,
+                                           const math_Vector& Tolerance,
+                                           const Standard_Integer NbIterations) :
+Delta(1, F.NbVariables()),
+Sol(1, F.NbVariables()),
+DF(1, F.NbEquations(),
+   1, StartingPoint.Length()),
+   Tol(1,F.NbVariables()),
 
-                                           InfBound(1, F.NbVariables()),
-                                           SupBound(1, F.NbVariables()),
+   InfBound(1, F.NbVariables()),
+   SupBound(1, F.NbVariables()),
 
-                                           SolSave(1, F.NbVariables()), 
-					   GH(1, F.NbVariables()), 
-					   DH(1, F.NbVariables()), 
-					   DHSave(1, F.NbVariables()),
-					   FF(1, F.NbEquations()), 
-					   PreviousSolution(1, F.NbVariables()), 
-					   Save(0, NbIterations),
-					   Constraints(1, F.NbVariables()),
-					   Temp1(1, F.NbVariables()), 
-					   Temp2(1, F.NbVariables()), 
-					   Temp3(1, F.NbVariables()), 
-					   Temp4(1, F.NbEquations())
+   SolSave(1, F.NbVariables()), 
+   GH(1, F.NbVariables()), 
+   DH(1, F.NbVariables()), 
+   DHSave(1, F.NbVariables()),
+   FF(1, F.NbEquations()), 
+   PreviousSolution(1, F.NbVariables()), 
+   Save(0, NbIterations),
+   Constraints(1, F.NbVariables()),
+   Temp1(1, F.NbVariables()), 
+   Temp2(1, F.NbVariables()), 
+   Temp3(1, F.NbVariables()), 
+   Temp4(1, F.NbEquations())
 
- {
+{
   for (Standard_Integer i = 1; i <= Tol.Length(); i++) {
     Tol(i) = Tolerance(i);
   }
@@ -704,17 +704,17 @@ void math_FunctionSetRoot::SetTolerance(const math_Vector& Tolerance)
 }
 
 void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
-				   const math_Vector& StartingPoint,
-	     		           const math_Vector& InfBound,
-			           const math_Vector& SupBound,
-			           Standard_Boolean theStopOnDivergent)
+                                   const math_Vector& StartingPoint,
+                                   const math_Vector& InfBound,
+                                   const math_Vector& SupBound,
+                                   Standard_Boolean theStopOnDivergent)
 {
   Standard_Integer Ninc = F.NbVariables(), Neq = F.NbEquations();
-  
+
   if ((Neq <= 0)                      ||
-      (StartingPoint.Length()!= Ninc) ||
-      (InfBound.Length() != Ninc)     ||
-      (SupBound.Length() != Ninc))  { Standard_DimensionError:: Raise(); }
+    (StartingPoint.Length()!= Ninc) ||
+    (InfBound.Length() != Ninc)     ||
+    (SupBound.Length() != Ninc))  { Standard_DimensionError:: Raise(); }
 
   Standard_Integer i;
   Standard_Boolean ChangeDirection = Standard_False, Sort = Standard_False, isNewSol = Standard_False;
@@ -726,10 +726,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
   math_Vector InvLengthMax(1, Ninc); // Pour bloquer les pas a 1/4 du domaine
   math_IntegerVector Constraints(1, Ninc); // Pour savoir sur quels bord on se trouve
   for (i = 1; i <= Ninc ; i++) {
-// modified by NIZHNY-MKK  Mon Oct  3 18:03:50 2005
-//      InvLengthMax(i) = 1. / Max(Abs(SupBound(i) - InfBound(i))/4, 1.e-9);
+    // modified by NIZHNY-MKK  Mon Oct  3 18:03:50 2005
+    //      InvLengthMax(i) = 1. / Max(Abs(SupBound(i) - InfBound(i))/4, 1.e-9);
     InvLengthMax(i) = 1. / Max((SupBound(i) - InfBound(i))/4, 1.e-9);
-   }
+  }
 
   MyDirFunction F_Dir(Temp1, Temp2, Temp3, Temp4, F);
   Standard_Integer  DescenteIter;
@@ -770,8 +770,8 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
   // de faire une seconde iteration...
   Save(0) = Max (F2, EpsSqrt);
   Standard_Real aTol_Func = Epsilon(F2);
-  FSR_DEBUG("=== Mode Debug de Function Set Root" << endl)
-  FSR_DEBUG("    F2 Initial = " << F2)
+  FSR_DEBUG("=== Mode Debug de Function Set Root" << endl);
+  FSR_DEBUG("    F2 Initial = " << F2);
 
   if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
     Done = Standard_False;
@@ -825,12 +825,12 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
     }
     //
     Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-		  Constraints, Delta, isNewSol);
+      Constraints, Delta, isNewSol);
+
 
-      
     DHSave = GH;
     if (isNewSol) {
-//      F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+      //      F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
       if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
         Done = Standard_False;
         if (!theStopOnDivergent || !myIsDivergent)
@@ -840,10 +840,10 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
         return;
       }
     }
-    
-    FSR_DEBUG("Kount         = " << Kount)
-    FSR_DEBUG("Le premier F2 = " << F2)
-    FSR_DEBUG("Direction     = " << ChangeDirection)
+
+    FSR_DEBUG("Kount         = " << Kount);
+    FSR_DEBUG("Le premier F2 = " << F2);
+    FSR_DEBUG("Direction     = " << ChangeDirection);
 
     if ((F2 <= Eps) || (Gnr1 <= Eps2)) {
       Done = Standard_False;
@@ -870,74 +870,74 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
       // Traitement standard sans traitement des bords
       // -------------------------------------------------
       if (!Sort) { // si l'on est pas sortie on essaye de progresser en avant
-	while((F2/PreviousMinimum > Progres) && !Stop) {
-	  if (F2 < OldF && (Dy < 0.0)) {
-	    // On essaye de progresser dans cette direction.
-	    FSR_DEBUG(" iteration de descente = " << DescenteIter)
-	    DescenteIter++;
-	    SolSave = Sol;
-	    OldF = F2;
-	    for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
-	      Sol(i) = Sol(i) + Ambda * DH(i);
-	    }
-      //
-      for (i = 1; i <= Ninc; i++)
-      {
-        myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
-      }
-      if (theStopOnDivergent & myIsDivergent)
-      {
-        return;
-      }
-      //
-	    Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave, 
-			  Constraints, Delta, isNewSol);
-	    FSR_DEBUG(" Augmentation de lambda")
-	    Ambda *= 1.7;
-	  }
-	  else {
-	    if ((F2 >= OldF)||(F2 >= PreviousMinimum)) {
-	      Good = Standard_False;
-	      if (DescenteIter == 0) { 
-		// C'est le premier pas qui flanche, on fait une interpolation.
-		// et on minimise si necessaire.
-		DescenteIter++;
-		Good = MinimizeDirection(SolSave, Delta, OldF, F2, DHSave, GH,
-					 Tol, F_Dir);
-	      }
-	      else if (ChangeDirection || (DescenteIter>1) 
-		       || (OldF>PreviousMinimum) ) {
-		// La progression a ete utile, on minimise...
-		DescenteIter++;
-		Good = MinimizeDirection(PreviousSolution, SolSave, Sol, 
-					 OldF, Delta, Tol, F_Dir);
-	      }
-	      if (!Good) {
-		Sol = SolSave;
-		F2 = OldF;
-	      }
-	      else {
-	        Sol = SolSave+Delta;
-          //
-          for (i = 1; i <= Ninc; i++)
-          {
-            myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+        while((F2/PreviousMinimum > Progres) && !Stop) {
+          if (F2 < OldF && (Dy < 0.0)) {
+            // On essaye de progresser dans cette direction.
+            FSR_DEBUG(" iteration de descente = " << DescenteIter);
+            DescenteIter++;
+            SolSave = Sol;
+            OldF = F2;
+            for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
+              Sol(i) = Sol(i) + Ambda * DH(i);
+            }
+            //
+            for (i = 1; i <= Ninc; i++)
+            {
+              myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+            }
+            if (theStopOnDivergent & myIsDivergent)
+            {
+              return;
+            }
+            //
+            Stop = Bounds(InfBound, SupBound, Tol, Sol, SolSave, 
+              Constraints, Delta, isNewSol);
+            FSR_DEBUG(" Augmentation de lambda");
+            Ambda *= 1.7;
           }
-          if (theStopOnDivergent & myIsDivergent)
-          {
-            return;
+          else {
+            if ((F2 >= OldF)||(F2 >= PreviousMinimum)) {
+              Good = Standard_False;
+              if (DescenteIter == 0) { 
+                // C'est le premier pas qui flanche, on fait une interpolation.
+                // et on minimise si necessaire.
+                DescenteIter++;
+                Good = MinimizeDirection(SolSave, Delta, OldF, F2, DHSave, GH,
+                  Tol, F_Dir);
+              }
+              else if (ChangeDirection || (DescenteIter>1) 
+                || (OldF>PreviousMinimum) ) {
+                  // La progression a ete utile, on minimise...
+                  DescenteIter++;
+                  Good = MinimizeDirection(PreviousSolution, SolSave, Sol, 
+                    OldF, Delta, Tol, F_Dir);
+              }
+              if (!Good) {
+                Sol = SolSave;
+                F2 = OldF;
+              }
+              else {
+                Sol = SolSave+Delta;
+                //
+                for (i = 1; i <= Ninc; i++)
+                {
+                  myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+                }
+                if (theStopOnDivergent & myIsDivergent)
+                {
+                  return;
+                }
+                //
+                Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+                  Constraints, Delta, isNewSol);
+              }
+              Sort = Standard_False; // On a rejete le point sur la frontiere
+            }
+            Stop = Standard_True; // et on sort dans tous les cas...
           }
-          //
-	        Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-	          Constraints, Delta, isNewSol);
-	      }
-	      Sort = Standard_False; // On a rejete le point sur la frontiere
-	    }
-	    Stop = Standard_True; // et on sort dans tous les cas...
-	  }
-	  DHSave = GH;
+          DHSave = GH;
           if (isNewSol) {
-//            F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+            //            F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
             if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
               Done = Standard_False;
               if (!theStopOnDivergent || !myIsDivergent)
@@ -947,173 +947,173 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
               return;
             }
           }
-	  Dy = GH*DH;
-	  if (Abs(Dy) <= Eps) {
-	    if (F2 > OldF)
+          Dy = GH*DH;
+          if (Abs(Dy) <= Eps) {
+            if (F2 > OldF)
               Sol = SolSave;
-      Done = Standard_False;
-      if (!theStopOnDivergent || !myIsDivergent)
-      {
-        Done = Standard_True;
+            Done = Standard_False;
+            if (!theStopOnDivergent || !myIsDivergent)
+            {
+              Done = Standard_True;
               ////modified by jgv, 31.08.2011////
               F.Value(Sol, FF); //update F before GetStateNumber
               ///////////////////////////////////
-        State = F.GetStateNumber();
-      }
-	    return;
-	  }
-	  if (DescenteIter >= 10) {
-	    Stop = Standard_True;
-	  }
-	}
-	FSR_DEBUG("--- Sortie du Traitement Standard")
-	FSR_DEBUG("    DescenteIter = "<<DescenteIter << " F2 = " << F2)
+              State = F.GetStateNumber();
+            }
+            return;
+          }
+          if (DescenteIter >= 100) {
+            Stop = Standard_True;
+          }
+        }
+        FSR_DEBUG("--- Sortie du Traitement Standard");
+        FSR_DEBUG("    DescenteIter = "<<DescenteIter << " F2 = " << F2);
       }
       // ------------------------------------
       //  on passe au traitement des bords
       // ------------------------------------
       if (Sort) {
-	Stop = (F2 > 1.001*OldF); // Pour ne pas progresser sur le bord
-	Sortbis = Sort;
-	DescenteIter = 0;
-	while (Sortbis && ((F2<OldF)|| (DescenteIter == 0))
-	       && (!Stop)) {
-	  DescenteIter++;
-	  // On essaye de progresser sur le bord
-	  SolSave = Sol;
-	  OldF = F2;
-	  SearchDirection(DF, GH, FF,  Constraints, Sol, 
-			  ChangeDirection, InvLengthMax, DH, Dy);
-	  FSR_DEBUG(" Conditional Direction = " << ChangeDirection)
-	  if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
-	    if (ChangeDirection) {
+        Stop = (F2 > 1.001*OldF); // Pour ne pas progresser sur le bord
+        Sortbis = Sort;
+        DescenteIter = 0;
+        while (Sortbis && ((F2<OldF)|| (DescenteIter == 0))
+          && (!Stop)) {
+            DescenteIter++;
+            // On essaye de progresser sur le bord
+            SolSave = Sol;
+            OldF = F2;
+            SearchDirection(DF, GH, FF,  Constraints, Sol, 
+              ChangeDirection, InvLengthMax, DH, Dy);
+            FSR_DEBUG(" Conditional Direction = " << ChangeDirection);
+            if (Dy<-Eps) { //Pour eviter des calculs inutiles et des /0...
+              if (ChangeDirection) {
 
-// 	      Ambda = Ambda2 / Sqrt(Abs(Dy));
-	      Ambda = Ambda2 / Sqrt(-Dy);
-	      if (Ambda > 1.0) Ambda = 1.0;
-	    }
-	    else {
-	      Ambda = 1.0;
-	      Ambda2 = 0.5*Ambda/DH.Norm();
-	    }
+                // 	      Ambda = Ambda2 / Sqrt(Abs(Dy));
+                Ambda = Ambda2 / Sqrt(-Dy);
+                if (Ambda > 1.0) Ambda = 1.0;
+              }
+              else {
+                Ambda = 1.0;
+                Ambda2 = 0.5*Ambda/DH.Norm();
+              }
 
-	    for( i = Sol.Lower(); i <= Sol.Upper(); i++) { 
-	      Sol(i) = Sol(i) + Ambda * DH(i);
-	    }
-      //
-      for (i = 1; i <= Ninc; i++)
-      {
-        myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
-      }
-      if (theStopOnDivergent & myIsDivergent)
-      {
-        return;
-      }
-      //
-	    Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-			     Constraints, Delta, isNewSol);
-
-	    DHSave = GH;
-            if (isNewSol) {
-//              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-              if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-                Done = Standard_False;
-                if (!theStopOnDivergent || !myIsDivergent)
-                {
-                  State = F.GetStateNumber();
-                }
+              for( i = Sol.Lower(); i <= Sol.Upper(); i++) { 
+                Sol(i) = Sol(i) + Ambda * DH(i);
+              }
+              //
+              for (i = 1; i <= Ninc; i++)
+              {
+                myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+              }
+              if (theStopOnDivergent & myIsDivergent)
+              {
                 return;
               }
+              //
+              Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+                Constraints, Delta, isNewSol);
+
+              DHSave = GH;
+              if (isNewSol) {
+                //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+                if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+                  Done = Standard_False;
+                  if (!theStopOnDivergent || !myIsDivergent)
+                  {
+                    State = F.GetStateNumber();
+                  }
+                  return;
+                }
+              }
+              Ambda2 = Gnr1;
+              FSR_DEBUG("---  Iteration au bords : " << DescenteIter);
+              FSR_DEBUG("---  F2 = " << F2);
+            }
+            else {
+              Stop = Standard_True;
+            }
+
+            while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
+              DescenteIter++;
+              FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter);
+              if (F2 < OldF && Dy < 0.0) {
+                // On essaye de progresser dans cette direction.
+                SolSave = Sol;
+                OldF = F2;
+                for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
+                  Sol(i) = Sol(i) + Ambda * DH(i);
+                }
+                //
+                for (i = 1; i <= Ninc; i++)
+                {
+                  myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+                }
+                if (theStopOnDivergent & myIsDivergent)
+                {
+                  return;
+                }
+                //
+                Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+                  Constraints, Delta, isNewSol);     
+              }
+              DHSave = GH;
+              if (isNewSol) {
+                //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+                if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+                  Done = Standard_False;
+                  if (!theStopOnDivergent || !myIsDivergent)
+                  {
+                    State = F.GetStateNumber();
+                  }
+                  return;
+                }
+              }
+              Ambda2 = Gnr1;
+              Dy = GH*DH;
+              Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
+            }
+            Stop = ((Dy >=0) || (DescenteIter >= 10));
+        }
+        if (((F2/PreviousMinimum > Progres) &&
+          (F2>=OldF))||(F2>=PreviousMinimum)) {
+            // On minimise par Brent
+            DescenteIter++;
+            Good = MinimizeDirection(SolSave, Delta, OldF, F2,  
+              DHSave, GH, Tol, F_Dir);
+            if (!Good) {
+              Sol = SolSave;
+              Sort = Standard_False;
+            }
+            else {
+              Sol = SolSave + Delta;
+              //
+              for (i = 1; i <= Ninc; i++)
+              {
+                myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
+              }
+              if (theStopOnDivergent & myIsDivergent)
+              {
+                return;
+              }
+              //
+              Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
+                Constraints, Delta, isNewSol);
+              if (isNewSol) {
+                //              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+                if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+                  Done = Standard_False;
+                  if (!theStopOnDivergent || !myIsDivergent)
+                  {
+                    State = F.GetStateNumber();
+                  }
+                  return;
+                }
+              }
             }
-	    Ambda2 = Gnr1;
-	    FSR_DEBUG("---  Iteration au bords : " << DescenteIter)
-            FSR_DEBUG("---  F2 = " << F2)
-	  }
-	  else {
-	    Stop = Standard_True;
-	  }
-
-	  while((F2/PreviousMinimum > Progres) && (F2<OldF) && (!Stop) ) {
-	    DescenteIter++;
-	    FSR_DEBUG("--- Iteration de descente conditionnel = " << DescenteIter)
-	    if (F2 < OldF && Dy < 0.0) {
-	      // On essaye de progresser dans cette direction.
-	      SolSave = Sol;
-	      OldF = F2;
-	      for( i = Sol.Lower(); i <= Sol.Upper(); i++) {
-		Sol(i) = Sol(i) + Ambda * DH(i);
-	      }
-        //
-        for (i = 1; i <= Ninc; i++)
-        {
-          myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
-        }
-        if (theStopOnDivergent & myIsDivergent)
-        {
-          return;
-        }
-        //
-	      Sortbis = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-			       Constraints, Delta, isNewSol);     
-	    }
-	    DHSave = GH;
-            if (isNewSol) {
-//              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-              if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-                Done = Standard_False;
-                if (!theStopOnDivergent || !myIsDivergent)
-                {
-                  State = F.GetStateNumber();
-                }
-                return;
-              }
-            }
-	    Ambda2 = Gnr1;
-	    Dy = GH*DH;
-	    Stop = ((Dy >=0) || (DescenteIter >= 10) || Sortbis);
-	  }
-	  Stop = ((Dy >=0) || (DescenteIter >= 10));
-	}
-	if (((F2/PreviousMinimum > Progres) &&
-	     (F2>=OldF))||(F2>=PreviousMinimum)) {
-	  // On minimise par Brent
-	  DescenteIter++;
-	  Good = MinimizeDirection(SolSave, Delta, OldF, F2,  
-				   DHSave, GH, Tol, F_Dir);
-	  if (!Good) {
-	    Sol = SolSave;
-            Sort = Standard_False;
-	  }
-	  else {
-	    Sol = SolSave + Delta;
-      //
-      for (i = 1; i <= Ninc; i++)
-      {
-        myIsDivergent |= (Sol(i) < InfBound(i)) | (SupBound(i) < Sol(i));
-      }
-      if (theStopOnDivergent & myIsDivergent)
-      {
-        return;
-      }
-      //
-	    Sort = Bounds(InfBound, SupBound, Tol, Sol, SolSave,
-	      Constraints, Delta, isNewSol);
-            if (isNewSol) {
-//              F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-              if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-                Done = Standard_False;
-                if (!theStopOnDivergent || !myIsDivergent)
-                {
-                  State = F.GetStateNumber();
-                }
-                return;
-              }
-	    }
-	  }
-	  Dy = GH*DH;
-	}	
-	FSR_DEBUG("--- Sortie du Traitement des Bords")
-	FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2)
+            Dy = GH*DH;
+        }	
+        FSR_DEBUG("--- Sortie du Traitement des Bords");
+        FSR_DEBUG("--- DescenteIter = "<<DescenteIter << " F2 = " << F2);
       }
     }
 
@@ -1123,56 +1123,56 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
     Save(Kount) = F2; 
     // Est ce la solution ?
     if (ChangeDirection) Verif = Standard_True;
-                        // Gradient : Il faut eviter de boucler
+    // Gradient : Il faut eviter de boucler
     else {
       Verif = Standard_False;
       if (Kount > 1) { // Pour accelerer les cas quasi-quadratique
-	if (Save(Kount-1)<1.e-4*Save(Kount-2)) Verif = Standard_True;
+        if (Save(Kount-1)<1.e-4*Save(Kount-2)) Verif = Standard_True;
       }
       else Verif = (F2 < 1.e-6*Save(0)); //Pour les cas dejas solutions
     }
     if (Verif) {
       for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
-	Delta(i) = PreviousSolution(i) - Sol(i);
+        Delta(i) = PreviousSolution(i) - Sol(i);
       }
-      
+
       if (IsSolutionReached(F)) {
-	if (PreviousMinimum < F2) {
-	  Sol = SolSave;
-	}
-  Done = Standard_False;
-  if (!theStopOnDivergent || !myIsDivergent)
-  {
-    Done = Standard_True;
+        if (PreviousMinimum < F2) {
+          Sol = SolSave;
+        }
+        Done = Standard_False;
+        if (!theStopOnDivergent || !myIsDivergent)
+        {
+          Done = Standard_True;
           ////modified by jgv, 31.08.2011////
           F.Value(Sol, FF); //update F before GetStateNumber
           ///////////////////////////////////
-    State = F.GetStateNumber();
-  }
-	return;
+          State = F.GetStateNumber();
+        }
+        return;
       }
     }
     //fin du test solution
-   
+
     // Analyse de la progression...
     //comparison of current minimum and previous minimum
     if ((F2 - PreviousMinimum) <= aTol_Func){ 
       if (Kount > 5) {
-	// L'historique est il bon ?
-	if (F2 >= 0.95*Save(Kount - 5)) {
-	  if (!ChangeDirection) ChangeDirection = Standard_True;
-	  else 
-    {
-      Done = Standard_False;
-      if (!theStopOnDivergent || !myIsDivergent)
-      {
-        Done = Standard_True;
-        State = F.GetStateNumber();
-      }
-      return; //  si un gain inf a 5% on sort
-    }
-	}
-	else ChangeDirection = Standard_False; //Si oui on recommence
+        // L'historique est il bon ?
+        if (F2 >= 0.95*Save(Kount - 5)) {
+          if (!ChangeDirection) ChangeDirection = Standard_True;
+          else 
+          {
+            Done = Standard_False;
+            if (!theStopOnDivergent || !myIsDivergent)
+            {
+              Done = Standard_True;
+              State = F.GetStateNumber();
+            }
+            return; //  si un gain inf a 5% on sort
+          }
+        }
+        else ChangeDirection = Standard_False; //Si oui on recommence
       }
       else  ChangeDirection = Standard_False; //Pas d'historique on continue
       // Si le gradient ne diminue pas suffisemment par newton on essaie
@@ -1180,47 +1180,47 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
       // paraitre avec NEWTON le gradient de f peut augmenter alors que f 
       // diminue: dans ce cas il faut garder NEWTON)
       if ((Gnr1 > 0.9*Oldgr) && 
-	  (F2 > 0.5*PreviousMinimum)) {
-	ChangeDirection = Standard_True;
+        (F2 > 0.5*PreviousMinimum)) {
+          ChangeDirection = Standard_True;
       }
 
       // Si l'on ne decide pas de changer de strategie, on verifie,
       // si ce n'est dejas fait     
       if ((!ChangeDirection) && (!Verif)) {
-	for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
-	  Delta(i) = PreviousSolution(i) - Sol(i);
-	}
-	if (IsSolutionReached(F)) {
-    Done = Standard_False;
-    if (!theStopOnDivergent || !myIsDivergent)
-    {
-      Done = Standard_True;
+        for(i = Delta.Lower(); i <= Delta.Upper(); i++) {
+          Delta(i) = PreviousSolution(i) - Sol(i);
+        }
+        if (IsSolutionReached(F)) {
+          Done = Standard_False;
+          if (!theStopOnDivergent || !myIsDivergent)
+          {
+            Done = Standard_True;
             ////modified by jgv, 31.08.2011////
             F.Value(Sol, FF); //update F before GetStateNumber
             ///////////////////////////////////
-      State = F.GetStateNumber();
-    }
-	  return;
-	}
+            State = F.GetStateNumber();
+          }
+          return;
+        }
       }
     } 
     else { // Cas de regression
       if (!ChangeDirection) { // On passe au gradient
-	ChangeDirection = Standard_True;
-	Sol = PreviousSolution;
-//	F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
-	if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
-	  Done = Standard_False;
-    if (!theStopOnDivergent || !myIsDivergent)
-    {
-      State = F.GetStateNumber();
-    }
-	  return;
-	}
+        ChangeDirection = Standard_True;
+        Sol = PreviousSolution;
+        //	F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1);
+        if(!F_Dir.Value(Sol, FF, DF, GH, F2, Gnr1)) {
+          Done = Standard_False;
+          if (!theStopOnDivergent || !myIsDivergent)
+          {
+            State = F.GetStateNumber();
+          }
+          return;
+        }
       }
       else 
       {
-        
+
         if (!theStopOnDivergent || !myIsDivergent)
         {
           State = F.GetStateNumber();
@@ -1239,23 +1239,23 @@ void math_FunctionSetRoot::Perform(math_FunctionSetWithDerivatives& F,
 
 
 Standard_Boolean math_FunctionSetRoot::IsSolutionReached(math_FunctionSetWithDerivatives& ) {
-   for(Standard_Integer i = 1; i<= Sol.Length(); i++) {
-     if(Abs(Delta(i)) > Tol(i)) {return Standard_False;}
-   }
-   return Standard_True;
+  for(Standard_Integer i = 1; i<= Sol.Length(); i++) {
+    if(Abs(Delta(i)) > Tol(i)) {return Standard_False;}
+  }
+  return Standard_True;
 }
 
 
 void math_FunctionSetRoot::Dump(Standard_OStream& o) const {
-   o<<" math_FunctionSetRoot";
-   if (Done) {
-     o << " Status = Done\n";
-     o << " Location value = " << Sol << "\n";
-     o << " Number of iterations = " << Kount << "\n";
-   }
-   else {
-     o<<"Status = Not Done\n";
-   }
+  o<<" math_FunctionSetRoot";
+  if (Done) {
+    o << " Status = Done\n";
+    o << " Location value = " << Sol << "\n";
+    o << " Number of iterations = " << Kount << "\n";
+  }
+  else {
+    o<<"Status = Not Done\n";
+  }
 }
 
 
diff --git a/tests/bugs/fclasses/bug24137 b/tests/bugs/fclasses/bug24137
new file mode 100755
index 0000000000..21e7f7256c
--- /dev/null
+++ b/tests/bugs/fclasses/bug24137
@@ -0,0 +1,22 @@
+puts "================"
+puts "OCC24137"
+puts "================"
+puts ""
+#######################################################################
+# math_FunctionSetRoot returns too rough solution
+#######################################################################
+
+pload QAcommands
+
+vertex v 6.65634 -0.201746 2.51477
+restore [locate_data_file bug24137_face.brep] f 
+
+OCC24137 f v 508.326 77.6999
+distmini d v result
+
+regexp {([-0-9.+eE]+)$} [dump d_val] full dist
+set good_dist 0
+set toler 2.0e-06
+if { [expr abs( ${dist} - ${good_dist} )] > ${toler} } {
+    puts "Faulty : the distanse is ${dist}. It is bad value"
+}
diff --git a/tests/bugs/modalg_1/bug10160_4 b/tests/bugs/modalg_1/bug10160_4
index c8283b9efe..c9f61a6334 100755
--- a/tests/bugs/modalg_1/bug10160_4
+++ b/tests/bugs/modalg_1/bug10160_4
@@ -39,7 +39,8 @@ puts "CPU_user_time=${CPU_user_time}"
 set CPU_user_time [expr ${CPU_user_time} / ${NbTests}]
 puts "CPU_user_time=${CPU_user_time}"
 
-set square 3.56087e+07
+#CR24137 set square 3.56087e+07
+set square 3.52471e+07
 
 # Analysis of "nbshapes res"
 set nb_v_good 24
@@ -53,4 +54,3 @@ set nb_compound_good 1
 set nb_shape_good 102
 
 set 2dviewer 0
-
diff --git a/tests/bugs/modalg_1/bug10160_5 b/tests/bugs/modalg_1/bug10160_5
index 87ca7f71e6..b74c12a213 100755
--- a/tests/bugs/modalg_1/bug10160_5
+++ b/tests/bugs/modalg_1/bug10160_5
@@ -39,7 +39,8 @@ puts "CPU_user_time=${CPU_user_time}"
 set CPU_user_time [expr ${CPU_user_time} / ${NbTests}]
 puts "CPU_user_time=${CPU_user_time}"
 
-set square 782201
+#CR24317 set square 782201
+set square 766474
 
 # Analysis of "nbshapes res"
 set nb_v_good 53
diff --git a/tests/de/iges_1/G9 b/tests/de/iges_1/G9
index 0c59d56410..082a390a35 100644
--- a/tests/de/iges_1/G9
+++ b/tests/de/iges_1/G9
@@ -6,11 +6,11 @@ set filename lh93wsddr3370z4.igs
 
 set ref_data {
 DATA        : Faulties = 0  ( 0 )  Warnings = 0  ( 0 )  Summary  = 0  ( 0 )
-TPSTAT      : Faulties = 0  ( 0 )  Warnings = 162  ( 4127 )  Summary  = 162  ( 4127 )
+TPSTAT      : Faulties = 0  ( 0 )  Warnings = 164  ( 4127 )  Summary  = 164  ( 4127 )
 CHECKSHAPE  : Wires    = 0  ( 0 )  Faces    = 0  ( 0 )  Shells   = 0  ( 0 )   Solids   = 0 ( 0 )
-NBSHAPES    : Solid    = 0  ( 0 )  Shell    = 0  ( 0 )  Face     = 448  ( 448 )   Summary  = 7237  ( 7237 )
-STATSHAPE   : Solid    = 0  ( 0 )  Shell    = 0  ( 0 )  Face     = 448  ( 448 )   FreeWire = 0  ( 0 )   FreeEdge  = 0 ( 0 )   SharedEdge = 3192  ( 3192 )
-TOLERANCE   : MaxTol   =  0.08683083502  (  0.04341528762 )  AvgTol   =  0.000870172875  (  0.0007502917279 )
+NBSHAPES    : Solid    = 0  ( 0 )  Shell    = 0  ( 0 )  Face     = 448  ( 448 )   Summary  = 7239  ( 7237 )
+STATSHAPE   : Solid    = 0  ( 0 )  Shell    = 0  ( 0 )  Face     = 448  ( 448 )   FreeWire = 0  ( 0 )   FreeEdge  = 0 ( 0 )   SharedEdge = 3193  ( 3192 )
+TOLERANCE   : MaxTol   =  0.08683083502  (  0.04341528762 )  AvgTol   =  0.0008697846249  (  0.0007558849075 )
 LABELS      : N0Labels = 1  ( 1 )  N1Labels = 0  ( 0 )  N2Labels = 0  ( 0 )   TotalLabels = 1  ( 1 )   NameLabels = 1  ( 1 )   ColorLabels = 0  ( 0 )   LayerLabels = 0  ( 0 )
 PROPS       : Centroid = 0  ( 0 )  Volume   = 0  ( 0 )  Area     = 0  ( 0 )
 NCOLORS     : NColors  = 0  ( 0 )