diff --git a/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx b/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx
index cea67e17a5..65a2de26ee 100644
--- a/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx
+++ b/src/ChFi3d/ChFi3d_Builder_CnCrn.cxx
@@ -437,7 +437,7 @@ static void CurveHermite (const TopOpeBRepDS_DataStructure& DStr,
     GeomAdaptor_Curve L (Bezier);
     Extrema_ExtCC ext (C,L);
     if (ext.IsDone()){ 
-      if (ext.NbExt()!=0){ 
+      if (!ext.IsParallel() && ext.NbExt()!=0){ 
         Extrema_POnCurv POnC, POnL;
         ext.Points(1, POnC, POnL);
         if (POnC.Value().Distance(POnL.Value()) < Precision::Confusion())
diff --git a/src/math/math_BFGS.cxx b/src/math/math_BFGS.cxx
index 6a5c450b34..64758dbb86 100644
--- a/src/math/math_BFGS.cxx
+++ b/src/math/math_BFGS.cxx
@@ -243,8 +243,8 @@ static Standard_Boolean MinimizeDirection(math_Vector&   P,
       // Make direction to go along the border
       for (Standard_Integer anIdx = 1; anIdx <= theLeft.Upper(); anIdx++)
       {
-        if (Abs(P(anIdx) - theRight(anIdx)) < Precision::PConfusion() ||
-          Abs(P(anIdx) - theLeft(anIdx)) < Precision::PConfusion())
+        if ((Abs(P(anIdx) - theRight(anIdx)) < Precision::PConfusion() && Dir(anIdx) > 0.0) ||
+            (Abs(P(anIdx) - theLeft(anIdx))  < Precision::PConfusion() && Dir(anIdx) < 0.0))
         {
           Dir(anIdx) = 0.0;
         }
diff --git a/src/math/math_GlobOptMin.cxx b/src/math/math_GlobOptMin.cxx
index 3aae44446a..d908eeff97 100644
--- a/src/math/math_GlobOptMin.cxx
+++ b/src/math/math_GlobOptMin.cxx
@@ -69,7 +69,8 @@ math_GlobOptMin::math_GlobOptMin(math_MultipleVarFunction* theFunc,
   myMaxV(1, myN),
   myCellSize(0, myN - 1),
   myFilter(theFunc->NbVariables()),
-  myCont(2)
+  myCont(2),
+  myF(Precision::Infinite())
 {
   Standard_Integer i;
 
@@ -382,17 +383,6 @@ void math_GlobOptMin::computeInitialValues()
     myC = Max(aLipConst * aMinEps, aMinLC);
   else if (aLipConst > myC * aMaxEps)
     myC = Min(myC * aMaxEps, aMaxLC);
-
-  // Clear all solutions except one.
-  if (myY.Size() != myN)
-  {
-    for(i = 1; i <= myN; i++)
-      aBestPnt(i) = myY(i);
-    myY.Clear();
-    for(i = 1; i <= myN; i++)
-      myY.Append(aBestPnt(i));
-  }
-  mySolCount = 1;
 }
 
 //=======================================================================
@@ -453,37 +443,8 @@ void math_GlobOptMin::computeGlobalExtremum(Standard_Integer j)
       aStepBestValue = (isInside && (val < d))? val : d;
       aStepBestPoint = (isInside && (val < d))? myTmp : myX;
 
-      // Solutions are close to each other 
-      // and it is allowed to have more than one solution.
-      if (Abs(aStepBestValue - myF) < mySameTol * 0.01 &&
-          !myIsFindSingleSolution)
-      {
-        if (!isStored(aStepBestPoint))
-        {
-          if ((aStepBestValue - myF) * myZ > 0.0)
-            myF = aStepBestValue;
-          for(i = 1; i <= myN; i++)
-            myY.Append(aStepBestPoint(i));
-          mySolCount++;
-        }
-      }
-
-      // New best solution:
-      // new point is out of (mySameTol * 0.01) surrounding or
-      // new point is better than old + single point search.
-      Standard_Real aFunctionalDelta = (aStepBestValue - myF) * myZ;
-      if (aFunctionalDelta > mySameTol * 0.01 ||
-         (aFunctionalDelta > 0.0 && myIsFindSingleSolution))
-      {
-        mySolCount = 0;
-        myF = aStepBestValue;
-        myY.Clear();
-        for(i = 1; i <= myN; i++)
-          myY.Append(aStepBestPoint(i));
-        mySolCount++;
-
-        isFirstCellFilterInvoke = Standard_True;
-      }
+      // Check point and value on the current step to be optimal.
+      checkAddCandidate(aStepBestPoint, aStepBestValue);
 
       if (CheckFunctionalStopCriteria())
         return; // Best possible value is obtained.
@@ -639,34 +600,59 @@ Standard_Boolean math_GlobOptMin::CheckFunctionalStopCriteria()
 //=======================================================================
 void math_GlobOptMin::ComputeInitSol()
 {
-  Standard_Real aCurrVal, aBestVal;
-  math_Vector aCurrPnt(1, myN);
-  math_Vector aBestPnt(1, myN);
-  math_Vector aParamStep(1, myN);
-  // Check functional value in midpoint, lower and upper border points and
-  // in each point try to perform local optimization.
-  aBestPnt = (myGlobA + myGlobB) * 0.5;
-  myFunc->Value(aBestPnt, aBestVal);
+  Standard_Real aVal;
+  math_Vector aPnt(1, myN);
 
-  Standard_Integer i;
-  for(i = 1; i <= 3; i++)
+  // Check functional value in midpoint. It is necessary since local optimization
+  // algorithm may fail and return nothing. This is a protection from uninitialized
+  // variables.
+  aPnt = (myGlobA + myGlobB) * 0.5;
+  myFunc->Value(aPnt, aVal);
+  checkAddCandidate(aPnt, aVal);
+
+  // Run local optimization from lower corner, midpoint, and upper corner.
+  for(Standard_Integer i = 1; i <= 3; i++)
   {
-    aCurrPnt = myA + (myB - myA) * (i - 1) / 2.0;
+    aPnt = myA + (myB - myA) * (i - 1) / 2.0;
 
-    if(computeLocalExtremum(aCurrPnt, aCurrVal, aCurrPnt))
+    if(computeLocalExtremum(aPnt, aVal, aPnt))
+      checkAddCandidate(aPnt, aVal);
+  }
+}
+
+//=======================================================================
+//function : checkAddCandidate
+//purpose  :
+//=======================================================================
+void math_GlobOptMin::checkAddCandidate(const math_Vector&  thePnt,
+                                        const Standard_Real theValue)
+{
+  if (Abs(theValue - myF) < mySameTol * 0.01 && // Value in point is close to optimal value.
+      !myIsFindSingleSolution)                  // Several optimal solutions are allowed.
+  {
+    if (!isStored(thePnt))
     {
-      // Local search tries to find better solution than current point.
-      if (aCurrVal < aBestVal)
-      {
-        aBestVal = aCurrVal;
-        aBestPnt = aCurrPnt;
-      }
+      if ((theValue - myF) * myZ > 0.0)
+        myF = theValue;
+      for (Standard_Integer j = 1; j <= myN; j++)
+        myY.Append(thePnt(j));
+      mySolCount++;
     }
   }
 
-  myF = aBestVal;
-  myY.Clear();
-  for(i = 1; i <= myN; i++)
-    myY.Append(aBestPnt(i));
-  mySolCount = 1;
+  // New best solution:
+  // new point is out of (mySameTol * 0.01) surrounding or
+  // new point is better than old and single point search.
+  Standard_Real aDelta = (theValue - myF) * myZ;
+  if (aDelta > mySameTol * 0.01 ||
+     (aDelta > 0.0 && myIsFindSingleSolution))
+  {
+    myF = theValue;
+    myY.Clear();
+    for (Standard_Integer j = 1; j <= myN; j++)
+      myY.Append(thePnt(j));
+    mySolCount = 1;
+
+    isFirstCellFilterInvoke = Standard_True;
+  }
 }
diff --git a/src/math/math_GlobOptMin.hxx b/src/math/math_GlobOptMin.hxx
index e98681181e..930301dc50 100644
--- a/src/math/math_GlobOptMin.hxx
+++ b/src/math/math_GlobOptMin.hxx
@@ -232,6 +232,13 @@ private:
   //! Check presence of thePnt in GlobOpt sequence.
   Standard_Boolean isStored(const math_Vector &thePnt);
 
+  //! Check and add candidate point into solutions.
+  //! @param thePnt   Candidate point.
+  //! @param theValue Function value in the candidate point.
+  void checkAddCandidate(const math_Vector&  thePnt,
+                         const Standard_Real theValue);
+
+
   // Input.
   math_MultipleVarFunction* myFunc;
   Standard_Integer myN;
diff --git a/src/math/math_NewtonMinimum.cxx b/src/math/math_NewtonMinimum.cxx
index 01a6178c46..b407501ab9 100644
--- a/src/math/math_NewtonMinimum.cxx
+++ b/src/math/math_NewtonMinimum.cxx
@@ -195,8 +195,8 @@ void math_NewtonMinimum::Perform(math_MultipleVarFunctionWithHessian& F,
           // Nullify corresponding TheStep indexes.
           for(Standard_Integer anIdx = 1; anIdx <= myLeft.Upper(); anIdx++)
           {
-            if (Abs(precedent->Value(anIdx) - myRight(anIdx)) < Precision::PConfusion() ||
-                Abs(precedent->Value(anIdx) - myLeft(anIdx) ) < Precision::PConfusion())
+            if ((Abs(precedent->Value(anIdx) - myRight(anIdx)) < Precision::PConfusion() && TheStep(anIdx) < 0.0) ||
+                (Abs(precedent->Value(anIdx) - myLeft(anIdx) ) < Precision::PConfusion() && TheStep(anIdx) > 0.0) )
             {
               TheStep(anIdx) = 0.0;
             }
diff --git a/tests/lowalgos/extcc/bug29858_01 b/tests/lowalgos/extcc/bug29858_01
new file mode 100644
index 0000000000..173d4f77ae
--- /dev/null
+++ b/tests/lowalgos/extcc/bug29858_01
@@ -0,0 +1,21 @@
+puts "========"
+puts "OCC29858"
+puts "========"
+puts ""
+#################################################
+# Regression in GeomAPI_ExtremaCurveCurve
+#################################################
+
+# Read input
+restore [locate_data_file bug29858_01_e1.brep] e1
+restore [locate_data_file bug29858_01_e2.brep] e2
+
+# Extract geometry from topology
+mkcurve c1 e1
+mkcurve c2 e2
+
+# Run extrema
+extrema c1 c2
+
+# Check result
+checklength ext_1 -l 9.4224516867962155e-10 -eps 0.01
\ No newline at end of file
diff --git a/tests/lowalgos/extcc/bug29858_02 b/tests/lowalgos/extcc/bug29858_02
new file mode 100644
index 0000000000..c150f970d3
--- /dev/null
+++ b/tests/lowalgos/extcc/bug29858_02
@@ -0,0 +1,28 @@
+puts "========"
+puts "OCC29858"
+puts "========"
+puts ""
+#################################################
+# Regression in GeomAPI_ExtremaCurveCurve
+#################################################
+
+# Read input
+restore [locate_data_file bug29858_02_e1.brep] e1
+restore [locate_data_file bug29858_02_e2.brep] e2
+
+# Extract geometry from topology
+mkcurve c1 e1
+mkcurve c2 e2
+
+# Run extrema
+set info [extrema c1 c2]
+
+# Check result
+regexp {Extrema 1 is point : +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $info full x y z
+# Point check
+set good_x 0.001
+set good_y 0.0073371170841145797
+set good_z -0.083369169789921913
+checkreal "Intersection point x:" ${x} ${good_x} 0.01 0.01
+checkreal "Intersection point y:" ${y} ${good_y} 0.01 0.01
+checkreal "Intersection point z:" ${z} ${good_z} 0.01 0.01
\ No newline at end of file
diff --git a/tests/lowalgos/extcc/bug29858_03 b/tests/lowalgos/extcc/bug29858_03
new file mode 100644
index 0000000000..7058727e53
--- /dev/null
+++ b/tests/lowalgos/extcc/bug29858_03
@@ -0,0 +1,28 @@
+puts "========"
+puts "OCC29858"
+puts "========"
+puts ""
+#################################################
+# Regression in GeomAPI_ExtremaCurveCurve
+#################################################
+
+# Read input
+restore [locate_data_file bug29858_03_e1.brep] e1
+restore [locate_data_file bug29858_03_e2.brep] e2
+
+# Extract geometry from topology
+mkcurve c1 e1
+mkcurve c2 e2
+
+# Run extrema
+set info [extrema c1 c2]
+
+# Check result
+regexp {Extrema 1 is point : +([-0-9.+eE]+) +([-0-9.+eE]+) +([-0-9.+eE]+)} $info full x y z
+# Point check
+set good_x 0.0
+set good_y 0.070710562195021642
+set good_z -0.65305318986891325
+checkreal "Intersection point x:" ${x} ${good_x} 0.01 0.01
+checkreal "Intersection point y:" ${y} ${good_y} 0.01 0.01
+checkreal "Intersection point z:" ${z} ${good_z} 0.01 0.01
\ No newline at end of file
diff --git a/tests/lowalgos/extcc/bug29858_04 b/tests/lowalgos/extcc/bug29858_04
new file mode 100644
index 0000000000..5afc248b21
--- /dev/null
+++ b/tests/lowalgos/extcc/bug29858_04
@@ -0,0 +1,27 @@
+puts "========"
+puts "OCC29858"
+puts "========"
+puts ""
+#################################################
+# Regression in GeomAPI_ExtremaCurveCurve
+#################################################
+
+# Read input
+restore [locate_data_file bug29858_04_e1.brep] e1
+restore [locate_data_file bug29858_04_e2.brep] e2
+
+# Extract geometry from topology
+mkcurve c1 e1
+mkcurve c2 e2
+
+# Run extrema
+set info [extrema c1 c2]
+
+# Amount Check
+if { [llength $info] != 2 } {
+   puts "Error: Invalid extrema number in extrema output"
+}
+
+# Check result
+checklength ext_1 -l 2.2204460492503131e-16 -eps 0.01
+checklength ext_2 -l 9.693032344962995e-16  -eps 0.01
\ No newline at end of file