diff --git a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
index 8fe2d216a3..5f1977d3aa 100644
--- a/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
+++ b/src/IntPatch/IntPatch_ImpImpIntersection_4.gxx
@@ -2691,7 +2691,62 @@ static IntPatch_ImpImpIntersection::IntStatus
     if ((aRangeS1.Delta() > aMaxV1Range) || (aRangeS2.Delta() > aMaxV2Range))
       return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve;
   }
+  //
+  Standard_Boolean isGoodIntersection = Standard_False;
+  Standard_Real anOptdu = 0.;
+  do
+  {
+    //Checking parameters of cylinders in order to define "good intersection"
+    //"Good intersection" means that axes of cylinders are almost perpendicular and
+    // one radius is much smaller than the other and small cylinder is "inside" big one.
+    const Standard_Real aToMuchCoeff = 3.;
+    const Standard_Real aCritAngle = M_PI / 18.; // 10 degree
+    Standard_Real anR1 = theCyl1.Radius();
+    Standard_Real anR2 = theCyl2.Radius();
+    Standard_Real anRmin = 0., anRmax = 0.;
+    //Radius criterion
+    if (anR1 > aToMuchCoeff * anR2)
+    {
+      anRmax = anR1; anRmin = anR2;
+    }
+    else if (anR2 > aToMuchCoeff * anR1)
+    {
+      anRmax = anR2; anRmin = anR1;
+    }
+    else
+    {
+      break;
+    }
+    //Angle criterion
+    const gp_Ax1& anAx1 = theCyl1.Axis();
+    const gp_Ax1& anAx2 = theCyl2.Axis();
+    if (!anAx1.IsNormal(anAx2, aCritAngle))
+    {
+      break;
+    }
+    //Placement criterion
+    gp_Lin anL1(anAx1), anL2(anAx2);
+    Standard_Real aDist = anL1.Distance(anL2);
+    if (aDist > anRmax / 2.)
+    {
+      break;
+    }
 
+    isGoodIntersection = Standard_True;
+    //Estimation of "optimal" du
+    //Relative deflection, absolut deflection is Rmin*aDeflection
+    Standard_Real aDeflection = 0.001;
+    Standard_Integer aNbP = 3;
+    if (anRmin * aDeflection > 1.e-3)
+    {
+      Standard_Real anAngle = 1.0e0 - aDeflection;
+      anAngle = 2.0e0 * ACos(anAngle);
+      aNbP = (Standard_Integer)(2. * M_PI / anAngle) + 1;
+    }
+    anOptdu = 2. * M_PI_2 / (Standard_Real)(aNbP - 1);
+
+  } while (0);
+//
   const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
   const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
   const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
@@ -2699,15 +2754,27 @@ static IntPatch_ImpImpIntersection::IntStatus
   const Standard_Real aTol2D = theBW.Get2dTolerance();
   const Standard_Real aTol3D = theBW.Get3dTolerance();
   const Standard_Real aPeriod = 2.0*M_PI;
-  const Standard_Integer aNbMaxPoints = 2000;
-  const Standard_Integer aNbMinPoints = 200;
-  const Standard_Integer aNbPoints = Min(Max(aNbMinPoints,
-                                  RealToInt(20.0*theCyl1.Radius())), aNbMaxPoints);
-  const Standard_Real aStepMin = aTol2D,
-                                  aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
-                                  (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) :
-                                  aUSurf1l - aUSurf1f;
+  Standard_Integer aNbMaxPoints = 1000;
+  Standard_Integer aNbMinPoints = 200;
+  Standard_Real du;
+  if (isGoodIntersection)
+  {
+    du = anOptdu;
+    aNbMaxPoints = 200;
+    aNbMinPoints = 50;
+  }
+  else
+  {
+    du = 2. * M_PI / aNbMaxPoints;
+  }
+  Standard_Integer aNbPts = Min(RealToInt((aUSurf1l - aUSurf1f) / du) + 1, 
+                                RealToInt(20.0*theCyl1.Radius()));
+  const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, aNbPts), aNbMaxPoints);
+  const Standard_Real aStepMin = Max(aTol2D, Precision::PConfusion()), 
+    aStepMax = (aUSurf1l - aUSurf1f > M_PI / 100.0) ?
+    (aUSurf1l - aUSurf1f) / IntToReal(aNbPoints) : aUSurf1l - aUSurf1f;
 
+ 
   //The main idea of the algorithm is to change U1-parameter
   //(U-parameter of theCyl1) from aU1f to aU1l with some step
   //(step is adaptive) and to obtain set of intersection points.
@@ -3279,6 +3346,16 @@ static IntPatch_ImpImpIntersection::IntStatus
               aMinUexp = Min(aMinUexp, anUexpect[i]);
               continue;
             }
+            //
+            if (isGoodIntersection)
+            {
+              //Use constant step
+              anUexpect[i] += aStepMax;
+              aMinUexp = Min(aMinUexp, anUexpect[i]);
+
+              continue;
+            }
+            //
 
             Standard_Real aStepTmp = aStepMax;
 
diff --git a/tests/lowalgos/intss/bug31552 b/tests/lowalgos/intss/bug31552
new file mode 100644
index 0000000000..f5f9a4b076
--- /dev/null
+++ b/tests/lowalgos/intss/bug31552
@@ -0,0 +1,48 @@
+puts "========="
+puts "0031552: Modeling Algorithms - Bad performance of intersection of cylindrical surfaces"
+puts "========="
+puts ""
+
+cpulimit 200
+
+cylinder s1 0 557.800010681152 0 0 0 1 0 1 0  347.850006103516
+trim s1 s1 3.2385 6.28319 -39.957 39.957
+
+cylinder s2 5.45639675193397 189.209310035068 0 -0.254420004160252 0.96709382248213 0 0.96709382248213 0.254420004160252 0 88.5
+trim s2 s2 6.19592 6.37045 0 10.0002
+
+dchrono z reset
+dchrono z start
+
+set ic 1
+while { $ic <= 100 } {
+  intersect res s1 s2 
+  incr ic
+}
+
+dchrono z stop counter Bug31552
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  puts "Error: only one curve"
+} 
+
+set che [whatis res_1]
+set ind [string first "3d curve" $che]
+if {${ind} < 0} {
+  puts "Error no first curve"
+}
+
+set che [whatis res_2]
+set ind [string first "3d curve" $che]
+if {${ind} < 0} {
+  puts "Error no second curve"
+}
+
+smallview AXON
+don res_1 res_2
+fit
+checkview -screenshot -2d -path ${imagedir}/${test_image}.png
+
+
diff --git a/tests/perf/modalg/bug26310_1 b/tests/perf/modalg/bug26310_1
index 1de235a3b8..06280519de 100644
--- a/tests/perf/modalg/bug26310_1
+++ b/tests/perf/modalg/bug26310_1
@@ -6,7 +6,7 @@ puts ""
 # Very slow boolean cut operations on cylinders
 #################################################
 
-set ExpTol 3.32e-07
+set ExpTol 3.05e-07
 
 set GoodNbCurv 4