0031499: Boolean Operations - Custom fuzzy value corrupts the result of CUT
[occt.git] / src / IntPatch / IntPatch_ImpImpIntersection_4.gxx
index d5cdd42..8fe2d21 100644 (file)
 // Alternatively, this file may be used under the terms of Open CASCADE
 // commercial license or contractual agreement.
 
+#include <algorithm>
+#include <Bnd_Range.hxx>
 #include <IntAna_ListOfCurve.hxx>
-#include <IntAna_ListIteratorOfListOfCurve.hxx>
-#include <IntPatch_WLine.hxx>
+#include <math_Matrix.hxx>
+#include <NCollection_IncAllocator.hxx>
+#include <Standard_DivideByZero.hxx>
 
+//If Abs(a) <= aNulValue then it is considered that a = 0.
+static const Standard_Real aNulValue = 1.0e-11;
+
+static void ShortCosForm( const Standard_Real theCosFactor,
+                          const Standard_Real theSinFactor,
+                          Standard_Real& theCoeff,
+                          Standard_Real& theAngle);
 //
-static 
-  Standard_Boolean ExploreCurve(const gp_Cylinder& aCy,
-                               const gp_Cone& aCo,
-                               IntAna_Curve& aC,
-                               const Standard_Real aTol,
-                               IntAna_ListOfCurve& aLC);
-static
-  Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
-                              const gp_Cylinder& Cy2,
-                              const Standard_Real Tol);
-
-// ------------------------------------------------------------------
-// MinMax : Replaces  theParMIN = MIN(theParMIN, theParMAX),
+static Standard_Boolean ExploreCurve(const gp_Cone& theCo,
+                                     IntAna_Curve& aC,
+                                     const Standard_Real aTol,
+                                     IntAna_ListOfCurve& aLC);
+
+static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
+                                      const Standard_Real theUlTarget,
+                                      Standard_Real& theUGiven,
+                                      const Standard_Real theTol2D,
+                                      const Standard_Real thePeriod,
+                                      const Standard_Boolean theFlForce);
+
+
+class ComputationMethods
+{
+  //Every cylinder can be represented by the following equation in parametric form:
+  //    S(U,V) = L + R*cos(U)*Xd+R*sin(U)*Yd+V*Zd,
+  //where location L, directions Xd, Yd and Zd have type gp_XYZ.
+
+  //Intersection points between two cylinders can be found from the following system:
+  //    S1(U1, V1) = S2(U2, V2)
+  //or
+  //    {X01 + R1*cos(U1)*Xx1 + R1*sin(U1)*Yx1 + V1*Zx1 = X02 + R2*cos(U2)*Xx2 + R2*sin(U2)*Yx2 + V2*Zx2
+  //    {Y01 + R1*cos(U1)*Xy1 + R1*sin(U1)*Yy1 + V1*Zy1 = Y02 + R2*cos(U2)*Xy2 + R2*sin(U2)*Yy2 + V2*Zy2   (1)
+  //    {Z01 + R1*cos(U1)*Xz1 + R1*sin(U1)*Yz1 + V1*Zz1 = Z02 + R2*cos(U2)*Xz2 + R2*sin(U2)*Yz2 + V2*Zz2
+
+  //The formula (1) can be rewritten as follows
+  //    {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1
+  //    {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2                                   (2)
+  //    {C13*V1+C23*V2=A13*cos(U1)+B13*sin(U1)+A23*cos(U2)+B23*sin(U2)+D3
+
+  //Hereafter we consider that in system
+  //    {C11*V1+C21*V2=A11*cos(U1)+B11*sin(U1)+A21*cos(U2)+B21*sin(U2)+D1                                   (3)
+  //    {C12*V1+C22*V2=A12*cos(U1)+B12*sin(U1)+A22*cos(U2)+B22*sin(U2)+D2
+  //variables V1 and V2 can be found unambiguously, i.e. determinant
+  //            |C11 C21|
+  //            |       | != 0
+  //            |C12 C22|
+  //            
+  //In this case, variables V1 and V2 can be found as follows:
+  //    {V1 = K11*sin(U1)+K21*sin(U2)+L11*cos(U1)+L21*cos(U2)+M1 = K1*cos(U1-FIV1)+L1*cos(U2-PSIV1)+M1      (4)
+  //    {V2 = K12*sin(U1)+K22*sin(U2)+L12*cos(U1)+L22*cos(U2)+M2 = K2*cos(U2-FIV2)+L2*cos(U2-PSIV2)+M2
+
+  //Having substituted result of (4) to the 3rd equation of (2), we will obtain equation
+  //    cos(U2-FI2) = B*cos(U1-FI1)+C.                                                                      (5)
+
+  //I.e. when U1 is taken different given values (from domain), correspond U2 value can be computed
+  //from equation (5). After that, V1 and V2 can be computed from the system (4) (see
+  //CylCylComputeParameters(...) methods).
+
+  //It is important to remark that equation (5) (in general) has two solutions: U2=FI2 +/- f(U1).
+  //Therefore, we are getting here two intersection lines.
+
+public:
+  //Stores equations coefficients
+  struct stCoeffsValue
+  {
+    stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
+
+    math_Vector mVecA1;
+    math_Vector mVecA2;
+    math_Vector mVecB1;
+    math_Vector mVecB2;
+    math_Vector mVecC1;
+    math_Vector mVecC2;
+    math_Vector mVecD;
+
+    Standard_Real mK21; //sinU2
+    Standard_Real mK11; //sinU1
+    Standard_Real mL21; //cosU2
+    Standard_Real mL11;  //cosU1
+    Standard_Real mM1;  //Free member
+
+    Standard_Real mK22; //sinU2
+    Standard_Real mK12; //sinU1
+    Standard_Real mL22; //cosU2
+    Standard_Real mL12; //cosU1
+    Standard_Real mM2; //Free member
+
+    Standard_Real mK1;
+    Standard_Real mL1;
+    Standard_Real mK2;
+    Standard_Real mL2;
+
+    Standard_Real mFIV1;
+    Standard_Real mPSIV1;
+    Standard_Real mFIV2;
+    Standard_Real mPSIV2;
+
+    Standard_Real mB;
+    Standard_Real mC;
+    Standard_Real mFI1;
+    Standard_Real mFI2;
+  };
+
+
+  //! Determines, if U2(U1) function is increasing.
+  static Standard_Boolean CylCylMonotonicity(const Standard_Real theU1par,
+                                             const Standard_Integer theWLIndex,
+                                             const stCoeffsValue& theCoeffs,
+                                             const Standard_Real thePeriod,
+                                             Standard_Boolean& theIsIncreasing);
+
+  //! Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
+  //! esimates the tolerance of U2-computing (estimation result is
+  //! assigned to *theDelta value).
+  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+                                                const Standard_Integer theWLIndex,
+                                                const stCoeffsValue& theCoeffs,
+                                                Standard_Real& theU2,
+                                                Standard_Real* const theDelta = 0);
+
+  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1,
+                                                  const Standard_Real theU2,
+                                                  const stCoeffsValue& theCoeffs,
+                                                  Standard_Real& theV1,
+                                                  Standard_Real& theV2);
+
+  static Standard_Boolean CylCylComputeParameters(const Standard_Real theU1par,
+                                                  const Standard_Integer theWLIndex,
+                                                  const stCoeffsValue& theCoeffs,
+                                                  Standard_Real& theU2,
+                                                  Standard_Real& theV1,
+                                                  Standard_Real& theV2);
+
+};
+
+ComputationMethods::stCoeffsValue::stCoeffsValue(const gp_Cylinder& theCyl1,
+                                                 const gp_Cylinder& theCyl2):
+  mVecA1(-theCyl1.Radius()*theCyl1.XAxis().Direction().XYZ()),
+  mVecA2(theCyl2.Radius()*theCyl2.XAxis().Direction().XYZ()),
+  mVecB1(-theCyl1.Radius()*theCyl1.YAxis().Direction().XYZ()),
+  mVecB2(theCyl2.Radius()*theCyl2.YAxis().Direction().XYZ()),
+  mVecC1(theCyl1.Axis().Direction().XYZ()),
+  mVecC2(theCyl2.Axis().Direction().XYZ().Reversed()),
+  mVecD(theCyl2.Location().XYZ() - theCyl1.Location().XYZ())
+{
+  enum CoupleOfEquation
+  {
+    COENONE = 0,
+    COE12 = 1,
+    COE23 = 2,
+    COE13 = 3
+  }aFoundCouple = COENONE;
+
+
+  Standard_Real aDetV1V2 = 0.0;
+
+  const Standard_Real aDelta1 = mVecC1(1)*mVecC2(2)-mVecC1(2)*mVecC2(1); //1-2
+  const Standard_Real aDelta2 = mVecC1(2)*mVecC2(3)-mVecC1(3)*mVecC2(2); //2-3
+  const Standard_Real aDelta3 = mVecC1(1)*mVecC2(3)-mVecC1(3)*mVecC2(1); //1-3
+  const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
+  const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
+  const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
+
+  if(anAbsD1 >= anAbsD2)
+  {
+    if(anAbsD3 > anAbsD1)
+    {
+      aFoundCouple = COE13;
+      aDetV1V2 = aDelta3;
+    }
+    else
+    {
+      aFoundCouple = COE12;
+      aDetV1V2 = aDelta1;
+    }
+  }
+  else
+  {
+    if(anAbsD3 > anAbsD2)
+    {
+      aFoundCouple = COE13;
+      aDetV1V2 = aDelta3;
+    }
+    else
+    {
+      aFoundCouple = COE23;
+      aDetV1V2 = aDelta2;
+    }
+  }
+
+  // In point of fact, every determinant (aDelta1, aDelta2 and aDelta3) is
+  // cross-product between directions (i.e. sine of angle).
+  // If sine is too small then sine is (approx.) equal to angle itself.
+  // Therefore, in this case we should compare sine with angular tolerance. 
+  // This constant is used for check if axes are parallel (see constructor
+  // AxeOperator::AxeOperator(...) in IntAna_QuadQuadGeo.cxx file).
+  if(Abs(aDetV1V2) < Precision::Angular())
+  {
+    throw Standard_Failure("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
+  }
+
+  switch(aFoundCouple)
+  {
+  case COE12:
+    break;
+  case COE23:
+    {
+      math_Vector aVTemp(mVecA1);
+      mVecA1(1) = aVTemp(2);
+      mVecA1(2) = aVTemp(3);
+      mVecA1(3) = aVTemp(1);
+
+      aVTemp = mVecA2;
+      mVecA2(1) = aVTemp(2);
+      mVecA2(2) = aVTemp(3);
+      mVecA2(3) = aVTemp(1);
+
+      aVTemp = mVecB1;
+      mVecB1(1) = aVTemp(2);
+      mVecB1(2) = aVTemp(3);
+      mVecB1(3) = aVTemp(1);
+
+      aVTemp = mVecB2;
+      mVecB2(1) = aVTemp(2);
+      mVecB2(2) = aVTemp(3);
+      mVecB2(3) = aVTemp(1);
+
+      aVTemp = mVecC1;
+      mVecC1(1) = aVTemp(2);
+      mVecC1(2) = aVTemp(3);
+      mVecC1(3) = aVTemp(1);
+
+      aVTemp = mVecC2;
+      mVecC2(1) = aVTemp(2);
+      mVecC2(2) = aVTemp(3);
+      mVecC2(3) = aVTemp(1);
+
+      aVTemp = mVecD;
+      mVecD(1) = aVTemp(2);
+      mVecD(2) = aVTemp(3);
+      mVecD(3) = aVTemp(1);
+
+      break;
+    }
+  case COE13:
+    {
+      math_Vector aVTemp = mVecA1;
+      mVecA1(2) = aVTemp(3);
+      mVecA1(3) = aVTemp(2);
+
+      aVTemp = mVecA2;
+      mVecA2(2) = aVTemp(3);
+      mVecA2(3) = aVTemp(2);
+
+      aVTemp = mVecB1;
+      mVecB1(2) = aVTemp(3);
+      mVecB1(3) = aVTemp(2);
+
+      aVTemp = mVecB2;
+      mVecB2(2) = aVTemp(3);
+      mVecB2(3) = aVTemp(2);
+
+      aVTemp = mVecC1;
+      mVecC1(2) = aVTemp(3);
+      mVecC1(3) = aVTemp(2);
+
+      aVTemp = mVecC2;
+      mVecC2(2) = aVTemp(3);
+      mVecC2(3) = aVTemp(2);
+
+      aVTemp = mVecD;
+      mVecD(2) = aVTemp(3);
+      mVecD(3) = aVTemp(2);
+
+      break;
+    }
+  default:
+    break;
+  }
+
+  //------- For V1 (begin)
+  //sinU2
+  mK21 = (mVecC2(2)*mVecB2(1)-mVecC2(1)*mVecB2(2))/aDetV1V2;
+  //sinU1
+  mK11 = (mVecC2(2)*mVecB1(1)-mVecC2(1)*mVecB1(2))/aDetV1V2;
+  //cosU2
+  mL21 = (mVecC2(2)*mVecA2(1)-mVecC2(1)*mVecA2(2))/aDetV1V2;
+  //cosU1
+  mL11 = (mVecC2(2)*mVecA1(1)-mVecC2(1)*mVecA1(2))/aDetV1V2;
+  //Free member
+  mM1 =  (mVecC2(2)*mVecD(1)-mVecC2(1)*mVecD(2))/aDetV1V2;
+  //------- For V1 (end)
+
+  //------- For V2 (begin)
+  //sinU2
+  mK22 = (mVecC1(1)*mVecB2(2)-mVecC1(2)*mVecB2(1))/aDetV1V2;
+  //sinU1
+  mK12 = (mVecC1(1)*mVecB1(2)-mVecC1(2)*mVecB1(1))/aDetV1V2;
+  //cosU2
+  mL22 = (mVecC1(1)*mVecA2(2)-mVecC1(2)*mVecA2(1))/aDetV1V2;
+  //cosU1
+  mL12 = (mVecC1(1)*mVecA1(2)-mVecC1(2)*mVecA1(1))/aDetV1V2;
+  //Free member
+  mM2 = (mVecC1(1)*mVecD(2)-mVecC1(2)*mVecD(1))/aDetV1V2;
+  //------- For V1 (end)
+
+  ShortCosForm(mL11, mK11, mK1, mFIV1);
+  ShortCosForm(mL21, mK21, mL1, mPSIV1);
+  ShortCosForm(mL12, mK12, mK2, mFIV2);
+  ShortCosForm(mL22, mK22, mL2, mPSIV2);
+
+  const Standard_Real aA1=mVecC1(3)*mK21+mVecC2(3)*mK22-mVecB2(3), //sinU2
+    aA2=mVecC1(3)*mL21+mVecC2(3)*mL22-mVecA2(3), //cosU2
+    aB1=mVecB1(3)-mVecC1(3)*mK11-mVecC2(3)*mK12, //sinU1
+    aB2=mVecA1(3)-mVecC1(3)*mL11-mVecC2(3)*mL12; //cosU1
+
+  mC =mVecD(3) - mVecC1(3)*mM1 -mVecC2(3)*mM2;  //Free
+
+  Standard_Real aA = 0.0;
+
+  ShortCosForm(aB2,aB1,mB,mFI1);
+  ShortCosForm(aA2,aA1,aA,mFI2);
+
+  mB /= aA;
+  mC /= aA;
+}
+
+class WorkWithBoundaries
+{
+public:
+  enum SearchBoundType
+  {
+    SearchNONE = 0,
+    SearchV1 = 1,
+    SearchV2 = 2
+  };
+
+  struct StPInfo
+  {
+    StPInfo()
+    {
+      mySurfID = 0;
+      myU1 = RealLast();
+      myV1 = RealLast();
+      myU2 = RealLast();
+      myV2 = RealLast();
+    }
+
+    //Equal to 0 for 1st surface non-zero for 2nd one.
+    Standard_Integer mySurfID;
+
+    Standard_Real myU1;
+    Standard_Real myV1;
+    Standard_Real myU2;
+    Standard_Real myV2;
+
+    bool operator>(const StPInfo& theOther) const
+    {
+      return myU1 > theOther.myU1;
+    }
+
+    bool operator<(const StPInfo& theOther) const
+    {
+      return myU1 < theOther.myU1;
+    }
+
+    bool operator==(const StPInfo& theOther) const
+    {
+      return myU1 == theOther.myU1;
+    }
+  };
+
+  WorkWithBoundaries(const IntSurf_Quadric& theQuad1,
+                     const IntSurf_Quadric& theQuad2,
+                     const ComputationMethods::stCoeffsValue& theCoeffs,
+                     const Bnd_Box2d& theUVSurf1,
+                     const Bnd_Box2d& theUVSurf2,
+                     const Standard_Integer theNbWLines,
+                     const Standard_Real thePeriod,
+                     const Standard_Real theTol3D,
+                     const Standard_Real theTol2D,
+                     const Standard_Boolean isTheReverse) :
+      myQuad1(theQuad1), myQuad2(theQuad2), myCoeffs(theCoeffs),
+      myUVSurf1(theUVSurf1), myUVSurf2(theUVSurf2), myNbWLines(theNbWLines),
+      myPeriod(thePeriod), myTol3D(theTol3D), myTol2D(theTol2D),
+      myIsReverse(isTheReverse)
+  {
+  };
+
+  // Returns parameters of system solved while finding
+  // intersection line
+  const ComputationMethods::stCoeffsValue &SICoeffs() const
+  {
+    return myCoeffs;
+  }
+
+  // Returns quadric correspond to the index theIdx.
+  const IntSurf_Quadric& GetQSurface(const Standard_Integer theIdx) const
+  {
+    if (theIdx <= 1)
+      return myQuad1;
+
+    return myQuad2;
+  }
+
+  // Returns TRUE in case of reverting surfaces
+  Standard_Boolean IsReversed() const
+  {
+    return myIsReverse;
+  }
+
+  // Returns 2D-tolerance
+  Standard_Real Get2dTolerance() const
+  {
+    return myTol2D;
+  }
+
+  // Returns 3D-tolerance
+  Standard_Real Get3dTolerance() const
+  {
+    return myTol3D;
+  }
+
+  // Returns UV-bounds of 1st surface
+  const Bnd_Box2d& UVS1() const
+  {
+    return myUVSurf1;
+  }
+
+  // Returns UV-bounds of 2nd surface
+  const Bnd_Box2d& UVS2() const
+  {
+    return myUVSurf2;
+  }
+
+  void AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
+                        const Standard_Real theU1,
+                        const Standard_Real theU1Min,
+                        const Standard_Real theU2,
+                        const Standard_Real theV1,
+                        const Standard_Real theV1Prev,
+                        const Standard_Real theV2,
+                        const Standard_Real theV2Prev,
+                        const Standard_Integer theWLIndex,
+                        const Standard_Boolean theFlForce,
+                        Standard_Boolean& isTheFound1,
+                        Standard_Boolean& isTheFound2) const;
+  
+  static Standard_Boolean BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
+                                              const Standard_Real thePeriod,
+                                              Bnd_Range theURange[]);
+  
+  void BoundaryEstimation(const gp_Cylinder& theCy1,
+                          const gp_Cylinder& theCy2,
+                          Bnd_Range& theOutBoxS1,
+                          Bnd_Range& theOutBoxS2) const;
+
+protected:
+
+  //Solves equation (2) (see declaration of ComputationMethods class) in case,
+  //when V1 or V2 (is set by theSBType argument) is known (corresponds to the boundary
+  //and equal to theVzad) but U1 is unknown. Computation is made by numeric methods and
+  //requires initial values (theVInit, theInitU2 and theInitMainVar).
+  Standard_Boolean
+              SearchOnVBounds(const SearchBoundType theSBType,
+                              const Standard_Real theVzad,
+                              const Standard_Real theVInit,
+                              const Standard_Real theInitU2,
+                              const Standard_Real theInitMainVar,
+                              Standard_Real& theMainVariableValue) const;
+
+  const WorkWithBoundaries& operator=(const WorkWithBoundaries&);
+
+private:
+  friend class ComputationMethods;
+
+  const IntSurf_Quadric& myQuad1;
+  const IntSurf_Quadric& myQuad2;
+  const ComputationMethods::stCoeffsValue& myCoeffs;
+  const Bnd_Box2d& myUVSurf1;
+  const Bnd_Box2d& myUVSurf2;
+  const Standard_Integer myNbWLines;
+  const Standard_Real myPeriod;
+  const Standard_Real myTol3D;
+  const Standard_Real myTol2D;
+  const Standard_Boolean myIsReverse;
+};
+
+static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
+                                  const IntSurf_Quadric& theQuad2,
+                                  const Handle(IntSurf_LineOn2S)& theLine,
+                                  const ComputationMethods::stCoeffsValue& theCoeffs,
+                                  const Standard_Integer theWLIndex,
+                                  const Standard_Integer theMinNbPoints,
+                                  const Standard_Integer theStartPointOnLine,
+                                  const Standard_Integer theEndPointOnLine,
+                                  const Standard_Real theTol2D,
+                                  const Standard_Real thePeriodOfSurf2,
+                                  const Standard_Boolean isTheReverse);
+
+//=======================================================================
+//function : MinMax
+//purpose  : Replaces theParMIN = MIN(theParMIN, theParMAX),
 //                    theParMAX = MAX(theParMIN, theParMAX).
-// ------------------------------------------------------------------
+//=======================================================================
 static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
 {
   if(theParMIN > theParMAX)
@@ -44,6 +536,256 @@ static inline void MinMax(Standard_Real& theParMIN, Standard_Real& theParMAX)
   }
 }
 
+//=======================================================================
+//function : ExtremaLineLine
+//purpose  : Computes extrema between the given lines. Returns parameters
+//          on correspond curve (see correspond method for Extrema_ExtElC class). 
+//=======================================================================
+static inline void ExtremaLineLine(const gp_Ax1& theC1,
+                                   const gp_Ax1& theC2,
+                                   const Standard_Real theCosA,
+                                   const Standard_Real theSqSinA,
+                                   Standard_Real& thePar1,
+                                   Standard_Real& thePar2)
+{
+  const gp_Dir &aD1 = theC1.Direction(),
+               &aD2 = theC2.Direction();
+
+  const gp_XYZ aL1L2 = theC2.Location().XYZ() - theC1.Location().XYZ();
+  const Standard_Real aD1L = aD1.XYZ().Dot(aL1L2),
+                      aD2L = aD2.XYZ().Dot(aL1L2);
+
+  thePar1 = (aD1L - theCosA * aD2L) / theSqSinA;
+  thePar2 = (theCosA * aD1L - aD2L) / theSqSinA;
+}
+
+//=======================================================================
+//function : VBoundaryPrecise
+//purpose  : By default, we shall consider, that V1 and V2 will be increased
+//            if U1 is increased. But if it is not, new V1set and/or V2set
+//            must be computed as [V1current - DeltaV1] (analogically
+//            for V2). This function processes this case.
+//=======================================================================
+static void VBoundaryPrecise( const math_Matrix& theMatr,
+                              const Standard_Real theV1AfterDecrByDelta,
+                              const Standard_Real theV2AfterDecrByDelta,
+                              Standard_Real& theV1Set,
+                              Standard_Real& theV2Set)
+{
+  //Now we are going to define if V1 (and V2) increases
+  //(or decreases) when U1 will increase.
+  const Standard_Integer aNbDim = 3;
+  math_Matrix aSyst(1, aNbDim, 1, aNbDim);
+
+  aSyst.SetCol(1, theMatr.Col(1));
+  aSyst.SetCol(2, theMatr.Col(2));
+  aSyst.SetCol(3, theMatr.Col(4));
+
+  //We have the system (see comment to StepComputing(...) function)
+  //    {a11*dV1 + a12*dV2 + a14*dU2 = -a13*dU1
+  //    {a21*dV1 + a22*dV2 + a24*dU2 = -a23*dU1
+  //    {a31*dV1 + a32*dV2 + a34*dU2 = -a33*dU1
+
+  const Standard_Real aDet = aSyst.Determinant();
+
+  aSyst.SetCol(1, theMatr.Col(3));
+  const Standard_Real aDet1 = aSyst.Determinant();
+
+  aSyst.SetCol(1, theMatr.Col(1));
+  aSyst.SetCol(2, theMatr.Col(3));
+
+  const Standard_Real aDet2 = aSyst.Determinant();
+
+  //Now,
+  //    dV1 = -dU1*aDet1/aDet
+  //    dV2 = -dU1*aDet2/aDet
+
+  //If U1 is increased then dU1 > 0.
+  //If (aDet1/aDet > 0) then dV1 < 0 and 
+  //V1 will be decreased after increasing U1.
+
+  //We have analogical situation with V2-parameter. 
+
+  if(aDet*aDet1 > 0.0)
+  {
+    theV1Set = theV1AfterDecrByDelta;
+  }
+
+  if(aDet*aDet2 > 0.0)
+  {
+    theV2Set = theV2AfterDecrByDelta;
+  }
+}
+
+//=======================================================================
+//function : DeltaU1Computing
+//purpose  : Computes new step for U1 parameter.
+//=======================================================================
+static inline 
+        Standard_Boolean DeltaU1Computing(const math_Matrix& theSyst,
+                                          const math_Vector& theFree,
+                                          Standard_Real& theDeltaU1Found)
+{
+  Standard_Real aDet = theSyst.Determinant();
+
+  if(Abs(aDet) > aNulValue)
+  {
+    math_Matrix aSyst1(theSyst);
+    aSyst1.SetCol(2, theFree);
+
+    theDeltaU1Found = Abs(aSyst1.Determinant()/aDet);
+    return Standard_True;
+  }
+
+  return Standard_False;
+}
+
+//=======================================================================
+//function : StepComputing
+//purpose  : 
+//
+//Attention!!!:
+//            theMatr must have 3*5-dimension strictly.
+//            For system
+//                {a11*V1+a12*V2+a13*dU1+a14*dU2=b1; 
+//                {a21*V1+a22*V2+a23*dU1+a24*dU2=b2; 
+//                {a31*V1+a32*V2+a33*dU1+a34*dU2=b3; 
+//            theMatr must be following:
+//                (a11 a12 a13 a14 b1) 
+//                (a21 a22 a23 a24 b2) 
+//                (a31 a32 a33 a34 b3) 
+//=======================================================================
+static Standard_Boolean StepComputing(const math_Matrix& theMatr,
+                                      const Standard_Real theV1Cur,
+                                      const Standard_Real theV2Cur,
+                                      const Standard_Real theDeltaV1,
+                                      const Standard_Real theDeltaV2,
+                                      Standard_Real& theDeltaU1Found/*,
+                                      Standard_Real& theDeltaU2Found,
+                                      Standard_Real& theV1Found,
+                                      Standard_Real& theV2Found*/)
+{
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+  bool flShow = false;
+
+  if(flShow)
+  {
+    printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
+              theMatr(1,1), theMatr(1,2), theMatr(1,3), theMatr(1,4), theMatr(1,5));
+    printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
+              theMatr(2,1), theMatr(2,2), theMatr(2,3), theMatr(2,4), theMatr(2,5));
+    printf("{%+10.20f*V1 + %+10.20f*V2 + %+10.20f*dU1 + %+10.20f*dU2 = %+10.20f\n",
+              theMatr(3,1), theMatr(3,2), theMatr(3,3), theMatr(3,4), theMatr(3,5));
+  }
+#endif
+
+  Standard_Boolean isSuccess = Standard_False;
+  theDeltaU1Found/* = theDeltaU2Found*/ = RealLast();
+  //theV1Found = theV1set;
+  //theV2Found = theV2Set;
+  const Standard_Integer aNbDim = 3;
+
+  math_Matrix aSyst(1, aNbDim, 1, aNbDim);
+  math_Vector aFree(1, aNbDim);
+
+  //By default, increasing V1(U1) and V2(U1) functions is
+  //considered
+  Standard_Real aV1Set = theV1Cur + theDeltaV1,
+                aV2Set = theV2Cur + theDeltaV2;
+
+  //However, what is indeed?
+  VBoundaryPrecise( theMatr, theV1Cur - theDeltaV1,
+                    theV2Cur - theDeltaV2, aV1Set, aV2Set);
+
+  aSyst.SetCol(2, theMatr.Col(3));
+  aSyst.SetCol(3, theMatr.Col(4));
+
+  for(Standard_Integer i = 0; i < 2; i++)
+  {
+    if(i == 0)
+    {//V1 is known
+      aSyst.SetCol(1, theMatr.Col(2));
+      aFree.Set(1, aNbDim, theMatr.Col(5)-aV1Set*theMatr.Col(1));
+    }
+    else
+    {//i==1 => V2 is known
+      aSyst.SetCol(1, theMatr.Col(1));
+      aFree.Set(1, aNbDim, theMatr.Col(5)-aV2Set*theMatr.Col(2));
+    }
+
+    Standard_Real aNewDU = theDeltaU1Found;
+    if(DeltaU1Computing(aSyst, aFree, aNewDU))
+    {
+      isSuccess = Standard_True;
+      if(aNewDU < theDeltaU1Found)
+      {
+        theDeltaU1Found = aNewDU;
+      }
+    }
+  }
+
+  if(!isSuccess)
+  {
+    aFree = theMatr.Col(5) - aV1Set*theMatr.Col(1) - aV2Set*theMatr.Col(2);
+    math_Matrix aSyst1(1, aNbDim, 1, 2);
+    aSyst1.SetCol(1, aSyst.Col(2));
+    aSyst1.SetCol(2, aSyst.Col(3));
+
+    //Now we have overdetermined system.
+
+    const Standard_Real aDet1 = theMatr(1,3)*theMatr(2,4) - theMatr(2,3)*theMatr(1,4);
+    const Standard_Real aDet2 = theMatr(1,3)*theMatr(3,4) - theMatr(3,3)*theMatr(1,4);
+    const Standard_Real aDet3 = theMatr(2,3)*theMatr(3,4) - theMatr(3,3)*theMatr(2,4);
+    const Standard_Real anAbsD1 = Abs(aDet1);
+    const Standard_Real anAbsD2 = Abs(aDet2);
+    const Standard_Real anAbsD3 = Abs(aDet3);
+
+    if(anAbsD1 >= anAbsD2)
+    {
+      if(anAbsD1 >= anAbsD3)
+      {
+        //Det1
+        if(anAbsD1 <= aNulValue)
+          return isSuccess;
+
+        theDeltaU1Found = Abs(aFree(1)*theMatr(2,4) - aFree(2)*theMatr(1,4))/anAbsD1;
+        isSuccess = Standard_True;
+      }
+      else
+      {
+        //Det3
+        if(anAbsD3 <= aNulValue)
+          return isSuccess;
+
+        theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
+        isSuccess = Standard_True;
+      }
+    }
+    else
+    {
+      if(anAbsD2 >= anAbsD3)
+      {
+        //Det2
+        if(anAbsD2 <= aNulValue)
+          return isSuccess;
+
+        theDeltaU1Found = Abs(aFree(1)*theMatr(3,4) - aFree(3)*theMatr(1,4))/anAbsD2;
+        isSuccess = Standard_True;
+      }
+      else
+      {
+        //Det3
+        if(anAbsD3 <= aNulValue)
+          return isSuccess;
+
+        theDeltaU1Found = Abs(aFree(2)*theMatr(3,4) - aFree(3)*theMatr(2,4))/anAbsD3;
+        isSuccess = Standard_True;
+      }
+    }
+  }
+
+  return isSuccess;
+}
 
 //=======================================================================
 //function : ProcessBounds
@@ -89,6 +831,7 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
        if (!procf) {
          d=ptf.Distance(ptsol.Value());
          if (d <= Tol) {
+            ptsol.SetTolerance(Tol);
            if (!ptsol.IsMultiple()) {
              //-- le point ptsol (de aligold) est declare multiple sur aligold
              Multpoint = Standard_True;
@@ -107,6 +850,7 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
        }
        if (!procl) {
          if (ptl.Distance(ptsol.Value()) <= Tol) {
+            ptsol.SetTolerance(Tol);
            if (!ptsol.IsMultiple()) {
              Multpoint = Standard_True;
              ptsol.SetMultiple(Standard_True);
@@ -137,6 +881,8 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
       }
     }
   }
+  
+  ptsol.SetTolerance(Tol);
   if (!procf && !procl) {
     Quad1.Parameters(ptf,U1,V1);
     Quad2.Parameters(ptf,U2,V2);
@@ -184,22 +930,24 @@ void ProcessBounds(const Handle(IntPatch_ALine)& alig,          //-- ligne coura
     alig->SetLastPoint(alig->NbVertex());
   }
 }
+
 //=======================================================================
-//function : IntCyCy
-//purpose  : 
+//function : CyCyAnalyticalIntersect
+//purpose  : Checks if intersection curve is analytical (line, circle, ellipse)
+//            and returns these curves.
 //=======================================================================
-Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
-                        const IntSurf_Quadric& Quad2,
-                        const Standard_Real Tol,
-                        Standard_Boolean& Empty,
-                        Standard_Boolean& Same,
-                        Standard_Boolean& Multpoint,
-                        IntPatch_SequenceOfLine& slin,
-                        IntPatch_SequenceOfPoint& spnt)
-
+Standard_Boolean CyCyAnalyticalIntersect( const IntSurf_Quadric& Quad1,
+                                          const IntSurf_Quadric& Quad2,
+                                          const IntAna_QuadQuadGeo& theInter,
+                                          const Standard_Real Tol,
+                                          Standard_Boolean& Empty,
+                                          Standard_Boolean& Same,
+                                          Standard_Boolean& Multpoint,
+                                          IntPatch_SequenceOfLine& slin,
+                                          IntPatch_SequenceOfPoint& spnt)
 {
   IntPatch_Point ptsol;
-
+  
   Standard_Integer i;
 
   IntSurf_TypeTrans trans1,trans2;
@@ -211,39 +959,34 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
   gp_Cylinder Cy1(Quad1.Cylinder());
   gp_Cylinder Cy2(Quad2.Cylinder());
 
-  IntAna_QuadQuadGeo inter(Cy1,Cy2,Tol);
-
-  if (!inter.IsDone())
-  {
-    return Standard_False;
-  }
-
-  typint = inter.TypeInter();
-  Standard_Integer NbSol = inter.NbSolutions();
+  typint = theInter.TypeInter();
+  Standard_Integer NbSol = theInter.NbSolutions();
   Empty = Standard_False;
   Same  = Standard_False;
 
   switch (typint)
-  {
+      {
   case IntAna_Empty:
     {
       Empty = Standard_True;
-    }
+      }
     break;
 
   case IntAna_Same:
-    {
+      {
       Same  = Standard_True;
-    }
+      }
     break;
-    
+
   case IntAna_Point:
     {
-      gp_Pnt psol(inter.Point(1));
-      Standard_Real U1,V1,U2,V2;
-      Quad1.Parameters(psol,U1,V1);
-      Quad2.Parameters(psol,U2,V2);
+      gp_Pnt psol(theInter.Point(1));
       ptsol.SetValue(psol,Tol,Standard_True);
+
+      Standard_Real U1,V1,U2,V2;
+      Quad1.Parameters(psol, U1, V1);
+      Quad2.Parameters(psol, U2, V2);
+
       ptsol.SetParameters(U1,V1,U2,V2);
       spnt.Append(ptsol);
     }
@@ -254,10 +997,14 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       gp_Pnt ptref;
       if (NbSol == 1)
       { // Cylinders are tangent to each other by line
-        linsol = inter.Line(1);
+        linsol = theInter.Line(1);
         ptref = linsol.Location();
+        
+        //Radius-vectors
         gp_Dir crb1(gp_Vec(ptref,Cy1.Location()));
         gp_Dir crb2(gp_Vec(ptref,Cy2.Location()));
+
+        //outer normal lines
         gp_Vec norm1(Quad1.Normale(ptref));
         gp_Vec norm2(Quad2.Normale(ptref));
         IntSurf_Situation situcyl1;
@@ -265,6 +1012,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
 
         if (crb1.Dot(crb2) < 0.)
         { // centre de courbures "opposes"
+            //ATTENTION!!! 
+            //        Normal and Radius-vector of the 1st(!) cylinder
+            //        is used for judging what the situation of the 2nd(!)
+            //        cylinder is.
+
           if (norm1.Dot(crb1) > 0.)
           {
             situcyl2 = IntSurf_Inside;
@@ -284,7 +1036,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
           }
         }
         else
-        {
+          {
           if (Cy1.Radius() < Cy2.Radius())
           {
             if (norm1.Dot(crb1) > 0.)
@@ -334,9 +1086,11 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       {
         for (i=1; i <= NbSol; i++)
         {
-          linsol = inter.Line(i);
+          linsol = theInter.Line(i);
           ptref = linsol.Location();
           gp_Vec lsd = linsol.Direction();
+
+          //Theoretically, qwe = +/- 1.0.
           Standard_Real qwe = lsd.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
           if (qwe >0.00000001)
           {
@@ -365,31 +1119,34 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       gp_Vec Tgt;
       gp_Pnt ptref;
       IntPatch_Point pmult1, pmult2;
-      
-      elipsol = inter.Ellipse(1);
-      
+
+      elipsol = theInter.Ellipse(1);
+
       gp_Pnt pttang1(ElCLib::Value(0.5*M_PI, elipsol));
       gp_Pnt pttang2(ElCLib::Value(1.5*M_PI, elipsol));
-      
+
       Multpoint = Standard_True;
       pmult1.SetValue(pttang1,Tol,Standard_True);
       pmult2.SetValue(pttang2,Tol,Standard_True);
       pmult1.SetMultiple(Standard_True);
       pmult2.SetMultiple(Standard_True);
-      
+
       Standard_Real oU1,oV1,oU2,oV2;
-      Quad1.Parameters(pttang1,oU1,oV1); 
-      Quad2.Parameters(pttang1,oU2,oV2);
-      pmult1.SetParameters(oU1,oV1,oU2,oV2);
+      Quad1.Parameters(pttang1, oU1, oV1);
+      Quad2.Parameters(pttang1, oU2, oV2);
 
+      pmult1.SetParameters(oU1,oV1,oU2,oV2);
       Quad1.Parameters(pttang2,oU1,oV1); 
       Quad2.Parameters(pttang2,oU2,oV2);
+
       pmult2.SetParameters(oU1,oV1,oU2,oV2);
-      
-      // on traite la premiere ellipse
 
+      // on traite la premiere ellipse
+        
       //-- Calcul de la Transition de la ligne 
       ElCLib::D1(0.,elipsol,ptref,Tgt);
+
+      //Theoretically, qwe = +/- |Tgt|.
       Standard_Real qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
       if (qwe>0.00000001)
       {
@@ -402,7 +1159,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
         trans2 = IntSurf_Out;
       }
       else
-      { 
+      {
         trans1=trans2=IntSurf_Undecided; 
       }
 
@@ -420,6 +1177,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
         //
         Quad1.Parameters(aP, aU1, aV1); 
         Quad2.Parameters(aP, aU2, aV2);
+
         aIP.SetParameters(aU1, aV1, aU2, aV2);
         //
         aIP.SetParameter(0.);
@@ -443,7 +1201,7 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       //-- Transitions calculee au point 0    OK
       //
       // on traite la deuxieme ellipse
-      elipsol = inter.Ellipse(2);
+      elipsol = theInter.Ellipse(2);
 
       Standard_Real param1 = ElCLib::Parameter(elipsol,pttang1);
       Standard_Real param2 = ElCLib::Parameter(elipsol,pttang2);
@@ -462,6 +1220,8 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
       
       //-- Calcul des transitions de ligne pour la premiere ligne 
       ElCLib::D1(parampourtransition,elipsol,ptref,Tgt);      
+
+      //Theoretically, qwe = +/- |Tgt|.
       qwe=Tgt.DotCross(Quad2.Normale(ptref),Quad1.Normale(ptref));
       if(qwe> 0.00000001)
       {
@@ -489,8 +1249,10 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
         aIP.SetValue(aP,Tol,Standard_False);
         aIP.SetMultiple(Standard_False);
         //
-        Quad1.Parameters(aP, aU1, aV1); 
+
+        Quad1.Parameters(aP, aU1, aV1);
         Quad2.Parameters(aP, aU2, aV2);
+
         aIP.SetParameters(aU1, aV1, aU2, aV2);
         //
         aIP.SetParameter(0.);
@@ -509,111 +1271,14 @@ Standard_Boolean IntCyCy(const IntSurf_Quadric& Quad1,
     }
     break;
 
-  case IntAna_NoGeometricSolution:
-    {
-      Standard_Boolean bReverse;
-      Standard_Real U1,V1,U2,V2;
-      gp_Pnt psol;
-      //
-      bReverse=IsToReverse(Cy1, Cy2, Tol);
-      if (bReverse)
-      {
-        Cy2=Quad1.Cylinder();
-        Cy1=Quad2.Cylinder();
-      }
-      //
-      IntAna_IntQuadQuad anaint(Cy1,Cy2,Tol);
-      if (!anaint.IsDone())
-      {
-        return Standard_False;
-      }
-      
-      if (anaint.NbPnt() == 0 && anaint.NbCurve() == 0)
-      {
-        Empty = Standard_True;
-      }
-      else
-      {
-        NbSol = anaint.NbPnt();
-        for (i = 1; i <= NbSol; i++)
-        {
-          psol = anaint.Point(i);
-          Quad1.Parameters(psol,U1,V1);
-          Quad2.Parameters(psol,U2,V2);
-          ptsol.SetValue(psol,Tol,Standard_True);
-          ptsol.SetParameters(U1,V1,U2,V2);
-          spnt.Append(ptsol);
-        }
-
-        gp_Pnt ptvalid, ptf, ptl;
-        gp_Vec tgvalid;
-
-        Standard_Real first,last,para;
-        IntAna_Curve curvsol;
-        Standard_Boolean tgfound;
-        Standard_Integer kount;
-
-        NbSol = anaint.NbCurve();
-        for (i = 1; i <= NbSol; i++)
-        {
-          curvsol = anaint.Curve(i);
-          curvsol.Domain(first,last);
-          ptf = curvsol.Value(first);
-          ptl = curvsol.Value(last);
-
-          para = last;
-          kount = 1;
-          tgfound = Standard_False;
-
-          while (!tgfound)
-          {
-            para = (1.123*first + para)/2.123;
-            tgfound = curvsol.D1u(para,ptvalid,tgvalid);
-            if (!tgfound)
-            {
-              kount ++;
-              tgfound = kount > 5;
-            }
-          }
-
-          Handle(IntPatch_ALine) alig;
-          if (kount <= 5)
-          {
-            Standard_Real qwe = tgvalid.DotCross( Quad2.Normale(ptvalid),
-                                                  Quad1.Normale(ptvalid));
-            if(qwe>0.00000001)
-            { 
-              trans1 = IntSurf_Out;
-              trans2 = IntSurf_In;
-            }
-            else if(qwe<-0.00000001)
-            {
-              trans1 = IntSurf_In;
-              trans2 = IntSurf_Out;
-            }
-            else
-            { 
-              trans1=trans2=IntSurf_Undecided; 
-            }
-            alig = new IntPatch_ALine(curvsol,Standard_False,trans1,trans2);
-          }
-          else
-          {
-            alig = new IntPatch_ALine(curvsol,Standard_False);
-            //-- cout << "Transition indeterminee" << endl;
-          }
-
-          Standard_Boolean TempFalse1 = Standard_False;
-          Standard_Boolean TempFalse2 = Standard_False;
-
-          ProcessBounds(alig,slin,Quad1,Quad2,TempFalse1,ptf,first,
-                        TempFalse2,ptl,last,Multpoint,Tol);
-          slin.Append(alig);
-        }
-      }
-    }
-    break;
+  case IntAna_Parabola:
+  case IntAna_Hyperbola:
+    throw Standard_Failure("IntCyCy(): Wrong intersection type!");
 
+  case IntAna_Circle:
+    // Circle is useful when we will work with trimmed surfaces
+    // (two cylinders can be tangent by their basises, e.g. circle)
+  case IntAna_NoGeometricSolution:
   default:
     return Standard_False;
   }
@@ -637,356 +1302,326 @@ static void ShortCosForm( const Standard_Real theCosFactor,
   if(IsEqual(theCoeff, 0.0))
   {
     theAngle = 0.0;
-    return;
-  }
-
-  theAngle = acos(Abs(theCosFactor/theCoeff));
-
-  if(theSinFactor > 0.0)
-  {
-    if(IsEqual(theCosFactor, 0.0))
-    {
-      theAngle = M_PI/2.0;
-    }
-    else if(theCosFactor < 0.0)
-    {
-      theAngle = M_PI-theAngle;
-    }
-  }
-  else if(IsEqual(theSinFactor, 0.0))
-  {
-    if(theCosFactor < 0.0)
-    {
-      theAngle = M_PI;
-    }
-  }
-  if(theSinFactor < 0.0)
-  {
-    if(theCosFactor > 0.0)
-    {
-      theAngle = 2.0*M_PI-theAngle;
-    }
-    else if(IsEqual(theCosFactor, 0.0))
-    {
-      theAngle = 3.0*M_PI/2.0;
-    }
-    else if(theCosFactor < 0.0)
-    {
-      theAngle = M_PI+theAngle;
-    }
-  }
-}
-
-enum SearchBoundType
-{
-  SearchNONE = 0,
-  SearchV1 = 1,
-  SearchV2 = 2
-};
-
-//Stores equations coefficients
-struct stCoeffsValue
-{
-  stCoeffsValue(const gp_Cylinder&, const gp_Cylinder&);
-
-  gp_Vec mVecA1;
-  gp_Vec mVecA2;
-  gp_Vec mVecB1;
-  gp_Vec mVecB2;
-  gp_Vec mVecC1;
-  gp_Vec mVecC2;
-  gp_Vec mVecD;
-
-  Standard_Real mK21; //sinU2
-  Standard_Real mK11; //sinU1
-  Standard_Real mL21; //cosU2
-  Standard_Real mL11;  //cosU1
-  Standard_Real mM1;  //Free member
-
-  Standard_Real mK22; //sinU2
-  Standard_Real mK12; //sinU1
-  Standard_Real mL22; //cosU2
-  Standard_Real mL12; //cosU1
-  Standard_Real mM2; //Free member
-  
-  Standard_Real mK1;
-  Standard_Real mL1;
-  Standard_Real mK2;
-  Standard_Real mL2;
-
-  Standard_Real mFIV1;
-  Standard_Real mPSIV1;
-  Standard_Real mFIV2;
-  Standard_Real mPSIV2;
-
-  Standard_Real mB;
-  Standard_Real mC;
-  Standard_Real mFI1;
-  Standard_Real mFI2;
-};
-
-stCoeffsValue::stCoeffsValue( const gp_Cylinder& theCyl1,
-                              const gp_Cylinder& theCyl2)
-{
-  const Standard_Real aNulValue = 0.01*Precision::PConfusion();
-
-  mVecA1 = -theCyl1.Radius()*theCyl1.XAxis().Direction();
-  mVecA2 = theCyl2.Radius()*theCyl2.XAxis().Direction();
-
-  mVecB1 = -theCyl1.Radius()*theCyl1.YAxis().Direction();
-  mVecB2 = theCyl2.Radius()*theCyl2.YAxis().Direction();
-
-  mVecC1 = theCyl1.Axis().Direction();
-  mVecC2 = -(theCyl2.Axis().Direction());
-
-  mVecD = theCyl2.Location().XYZ() - theCyl1.Location().XYZ();
-
-  enum CoupleOfEquation
-  {
-    COENONE = 0,
-    COE12 = 1,
-    COE23 = 2,
-    COE13 = 3
-  }aFoundCouple = COENONE;
-
+    return;
+  }
 
-  Standard_Real aDetV1V2 = 0.0;
-  
-  const Standard_Real aDelta1 = mVecC1.X()*mVecC2.Y()-mVecC1.Y()*mVecC2.X(); //1-2
-  const Standard_Real aDelta2 = mVecC1.Y()*mVecC2.Z()-mVecC1.Z()*mVecC2.Y(); //2-3
-  const Standard_Real aDelta3 = mVecC1.X()*mVecC2.Z()-mVecC1.Z()*mVecC2.X(); //1-3
-  const Standard_Real anAbsD1 = Abs(aDelta1); //1-2
-  const Standard_Real anAbsD2 = Abs(aDelta2); //2-3
-  const Standard_Real anAbsD3 = Abs(aDelta3); //1-3
+  theAngle = acos(Abs(theCosFactor/theCoeff));
 
-  if(anAbsD1 >= anAbsD2)
+  if(theSinFactor > 0.0)
   {
-    if(anAbsD3 > anAbsD1)
+    if(IsEqual(theCosFactor, 0.0))
     {
-      aFoundCouple = COE13;
-      aDetV1V2 = aDelta3;
+      theAngle = M_PI/2.0;
     }
-    else
+    else if(theCosFactor < 0.0)
     {
-      aFoundCouple = COE12;
-      aDetV1V2 = aDelta1;
+      theAngle = M_PI-theAngle;
     }
   }
-  else
+  else if(IsEqual(theSinFactor, 0.0))
   {
-    if(anAbsD3 > anAbsD2)
-    {
-      aFoundCouple = COE13;
-      aDetV1V2 = aDelta3;
-    }
-    else
+    if(theCosFactor < 0.0)
     {
-      aFoundCouple = COE23;
-      aDetV1V2 = aDelta2;
+      theAngle = M_PI;
     }
   }
-
-  if(Abs(aDetV1V2) < aNulValue)
+  if(theSinFactor < 0.0)
   {
-    Standard_Failure::Raise("Error. Exception in divide by zerro (IntCyCyTrim)!!!!");
+    if(theCosFactor > 0.0)
+    {
+      theAngle = 2.0*M_PI-theAngle;
+    }
+    else if(IsEqual(theCosFactor, 0.0))
+    {
+      theAngle = 3.0*M_PI/2.0;
+    }
+    else if(theCosFactor < 0.0)
+    {
+      theAngle = M_PI+theAngle;
+    }
   }
+}
 
-  switch(aFoundCouple)
+//=======================================================================
+//function : CylCylMonotonicity
+//purpose  : Determines, if U2(U1) function is increasing.
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylMonotonicity(const Standard_Real theU1par,
+                                                        const Standard_Integer theWLIndex,
+                                                        const stCoeffsValue& theCoeffs,
+                                                        const Standard_Real thePeriod,
+                                                        Standard_Boolean& theIsIncreasing)
+{
+  // U2(U1) = FI2 + (+/-)acos(B*cos(U1 - FI1) + C);
+  //Accordingly,
+  //Func. U2(X1) = FI2 + X1;
+  //Func. X1(X2) = anArccosFactor*X2;
+  //Func. X2(X3) = acos(X3);
+  //Func. X3(X4) = B*X4 + C;
+  //Func. X4(U1) = cos(U1 - FI1).
+  //
+  //Consequently,
+  //U2(X1) is always increasing.
+  //X1(X2) is increasing if anArccosFactor > 0.0 and
+  //is decreasing otherwise.
+  //X2(X3) is always decreasing.
+  //Therefore, U2(X3) is decreasing if anArccosFactor > 0.0 and
+  //is increasing otherwise.
+  //X3(X4) is increasing if B > 0 and is decreasing otherwise.
+  //X4(U1) is increasing if U1 - FI1 in [PI, 2*PI) and
+  //is decreasing U1 - FI1 in [0, PI) or if (U1 - FI1 == 2PI).
+  //After that, we can predict behaviour of U2(U1) function:
+  //if it is increasing or decreasing.
+
+  //For "+/-" sign. If isPlus == TRUE, "+" is chosen, otherwise, we choose "-".
+  Standard_Boolean isPlus = Standard_False;
+
+  switch(theWLIndex)
   {
-  case COE12:
+  case 0: 
+    isPlus = Standard_True;
     break;
-  case COE23:
-    {
-      gp_Vec aVTemp = mVecA1;
-      mVecA1.SetX(aVTemp.Y());
-      mVecA1.SetY(aVTemp.Z());
-      mVecA1.SetZ(aVTemp.X());
+  case 1:
+    isPlus = Standard_False;
+    break;
+  default:
+    //throw Standard_Failure("Error. Range Error!!!!");
+    return Standard_False;
+  }
 
-      aVTemp = mVecA2;
-      mVecA2.SetX(aVTemp.Y());
-      mVecA2.SetY(aVTemp.Z());
-      mVecA2.SetZ(aVTemp.X());
+  Standard_Real aU1Temp = theU1par - theCoeffs.mFI1;
+  InscribePoint(0, thePeriod, aU1Temp, 0.0, thePeriod, Standard_False);
 
-      aVTemp = mVecB1;
-      mVecB1.SetX(aVTemp.Y());
-      mVecB1.SetY(aVTemp.Z());
-      mVecB1.SetZ(aVTemp.X());
-      
-      aVTemp = mVecB2;
-      mVecB2.SetX(aVTemp.Y());
-      mVecB2.SetY(aVTemp.Z());
-      mVecB2.SetZ(aVTemp.X());
+  theIsIncreasing = Standard_True;
 
-      aVTemp = mVecC1;
-      mVecC1.SetX(aVTemp.Y());
-      mVecC1.SetY(aVTemp.Z());
-      mVecC1.SetZ(aVTemp.X());
+  if(((M_PI - aU1Temp) < RealSmall()) && (aU1Temp < thePeriod))
+  {
+    theIsIncreasing = Standard_False;
+  }
 
-      aVTemp = mVecC2;
-      mVecC2.SetX(aVTemp.Y());
-      mVecC2.SetY(aVTemp.Z());
-      mVecC2.SetZ(aVTemp.X());
+  if(theCoeffs.mB < 0.0)
+  {
+    theIsIncreasing = !theIsIncreasing;
+  }
 
-      aVTemp = mVecD;
-      mVecD.SetX(aVTemp.Y());
-      mVecD.SetY(aVTemp.Z());
-      mVecD.SetZ(aVTemp.X());
+  if(!isPlus)
+  {
+    theIsIncreasing = !theIsIncreasing;
+  }  
 
-      break;
-    }
-  case COE13:
-    {
-      gp_Vec aVTemp = mVecA1;
-      mVecA1.SetY(aVTemp.Z());
-      mVecA1.SetZ(aVTemp.Y());
+  return Standard_True;
+}
 
-      aVTemp = mVecA2;
-      mVecA2.SetY(aVTemp.Z());
-      mVecA2.SetZ(aVTemp.Y());
+//=======================================================================
+//function : CylCylComputeParameters
+//purpose  : Computes U2 (U-parameter of the 2nd cylinder) and, if theDelta != 0,
+//            estimates the tolerance of U2-computing (estimation result is
+//            assigned to *theDelta value).
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
+                                                const Standard_Integer theWLIndex,
+                                                const stCoeffsValue& theCoeffs,
+                                                Standard_Real& theU2,
+                                                Standard_Real* const theDelta)
+{
+  //This formula is got from some experience and can be changed.
+  const Standard_Real aTol0 = Min(10.0*Epsilon(1.0)*theCoeffs.mB, aNulValue);
+  const Standard_Real aTol = 1.0 - aTol0;
 
-      aVTemp = mVecB1;
-      mVecB1.SetY(aVTemp.Z());
-      mVecB1.SetZ(aVTemp.Y());
+  if(theWLIndex < 0 || theWLIndex > 1)
+    return Standard_False;
 
-      aVTemp = mVecB2;
-      mVecB2.SetY(aVTemp.Z());
-      mVecB2.SetZ(aVTemp.Y());
+  const Standard_Real aSign = theWLIndex ? -1.0 : 1.0;
 
-      aVTemp = mVecC1;
-      mVecC1.SetY(aVTemp.Z());
-      mVecC1.SetZ(aVTemp.Y());
+  Standard_Real anArg = cos(theU1par - theCoeffs.mFI1);
+  anArg = theCoeffs.mB*anArg + theCoeffs.mC;
 
-      aVTemp = mVecC2;
-      mVecC2.SetY(aVTemp.Z());
-      mVecC2.SetZ(aVTemp.Y());
+  if(anArg >= aTol)
+  {
+    if(theDelta)
+      *theDelta = 0.0;
 
-      aVTemp = mVecD;
-      mVecD.SetY(aVTemp.Z());
-      mVecD.SetZ(aVTemp.Y());
+    anArg = 1.0;
+  }
+  else if(anArg <= -aTol)
+  {
+    if(theDelta)
+      *theDelta = 0.0;
 
-      break;
-    }
-  default:
-    break;
+    anArg = -1.0;
+  }
+  else if(theDelta)
+  {
+    //There is a case, when
+    //  const double aPar = 0.99999999999721167;
+    //  const double aFI2 = 2.3593296083566181e-006;
+
+    //Then
+    //  aPar - cos(aFI2) == -5.10703e-015 ==> cos(aFI2) == aPar. 
+    //Theoreticaly, in this case 
+    //  aFI2 == +/- acos(aPar).
+    //However,
+    //  acos(aPar) - aFI2 == 2.16362e-009.
+    //Error is quite big.
+
+    //This error should be estimated. Let use following way, which is based
+    //on expanding to Taylor series.
+
+    //  acos(p)-acos(p+x) = x/sqrt(1-p*p).
+
+    //If p == (1-d) (when p > 0) or p == (-1+d) (when p < 0) then
+    //  acos(p)-acos(p+x) = x/sqrt(d*(2-d)).
+
+    //Here always aTol0 <= d <= 1. Max(x) is considered (!) to be equal to aTol0.
+    //In this case
+    //  8*aTol0 <= acos(p)-acos(p+x) <= sqrt(2/(2-aTol0)-1),
+    //                                              because 0 < aTol0 < 1.
+    //E.g. when aTol0 = 1.0e-11,
+    //   8.0e-11 <= acos(p)-acos(p+x) < 2.24e-6.
+
+    const Standard_Real aDelta = Min(1.0-anArg, 1.0+anArg);
+    Standard_DivideByZero_Raise_if((aDelta*aDelta < RealSmall()) || (aDelta >= 2.0), 
+          "IntPatch_ImpImpIntersection_4.gxx, CylCylComputeParameters()");
+    *theDelta = aTol0/sqrt(aDelta*(2.0-aDelta));
   }
 
-  //------- For V1 (begin)
-  //sinU2
-  mK21 = (mVecC2.Y()*mVecB2.X()-mVecC2.X()*mVecB2.Y())/aDetV1V2;
-  //sinU1
-  mK11 = (mVecC2.Y()*mVecB1.X()-mVecC2.X()*mVecB1.Y())/aDetV1V2;
-  //cosU2
-  mL21 = (mVecC2.Y()*mVecA2.X()-mVecC2.X()*mVecA2.Y())/aDetV1V2;
-  //cosU1
-  mL11 = (mVecC2.Y()*mVecA1.X()-mVecC2.X()*mVecA1.Y())/aDetV1V2;
-  //Free member
-  mM1 = (mVecC2.Y()*mVecD.X()-mVecC2.X()*mVecD.Y())/aDetV1V2;
-  //------- For V1 (end)
+  theU2 = acos(anArg);
+  theU2 = theCoeffs.mFI2 + aSign*theU2;
 
-  //------- For V2 (begin)
-  //sinU2
-  mK22 = (mVecC1.X()*mVecB2.Y()-mVecC1.Y()*mVecB2.X())/aDetV1V2;
-  //sinU1
-  mK12 = (mVecC1.X()*mVecB1.Y()-mVecC1.Y()*mVecB1.X())/aDetV1V2;
-  //cosU2
-  mL22 = (mVecC1.X()*mVecA2.Y()-mVecC1.Y()*mVecA2.X())/aDetV1V2;
-  //cosU1
-  mL12 = (mVecC1.X()*mVecA1.Y()-mVecC1.Y()*mVecA1.X())/aDetV1V2;
-  //Free member
-  mM2 = (mVecC1.X()*mVecD.Y()-mVecC1.Y()*mVecD.X())/aDetV1V2;
-  //------- For V1 (end)
+  return Standard_True;
+}
 
-  ShortCosForm(mL11, mK11, mK1, mFIV1);
-  ShortCosForm(mL21, mK21, mL1, mPSIV1);
-  ShortCosForm(mL12, mK12, mK2, mFIV2);
-  ShortCosForm(mL22, mK22, mL2, mPSIV2);
+//=======================================================================
+//function : CylCylComputeParameters
+//purpose  : Computes V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1,
+                                                const Standard_Real theU2,
+                                                const stCoeffsValue& theCoeffs,
+                                                Standard_Real& theV1,
+                                                Standard_Real& theV2)
+{
+  theV1 = theCoeffs.mK21 * sin(theU2) + 
+          theCoeffs.mK11 * sin(theU1) +
+          theCoeffs.mL21 * cos(theU2) +
+          theCoeffs.mL11 * cos(theU1) + theCoeffs.mM1;
 
-  const Standard_Real aA1=mVecC1.Z()*mK21+mVecC2.Z()*mK22-mVecB2.Z(), //sinU2
-                      aA2=mVecC1.Z()*mL21+mVecC2.Z()*mL22-mVecA2.Z(), //cosU2
-                      aB1=mVecB1.Z()-mVecC1.Z()*mK11-mVecC2.Z()*mK12, //sinU1
-                      aB2=mVecA1.Z()-mVecC1.Z()*mL11-mVecC2.Z()*mL12; //cosU1
+  theV2 = theCoeffs.mK22 * sin(theU2) +
+          theCoeffs.mK12 * sin(theU1) +
+          theCoeffs.mL22 * cos(theU2) +
+          theCoeffs.mL12 * cos(theU1) + theCoeffs.mM2;
 
-  mC =mVecD.Z() -mVecC1.Z()*mM1 -mVecC2.Z()*mM2;  //Free
+  return Standard_True;
+}
 
-  Standard_Real aA = 0.0;
+//=======================================================================
+//function : CylCylComputeParameters
+//purpose  : Computes U2 (U-parameter of the 2nd cylinder),
+//            V1 and V2 (V-parameters of the 1st and 2nd cylinder respectively).
+//=======================================================================
+Standard_Boolean ComputationMethods::CylCylComputeParameters(const Standard_Real theU1par,
+                                                const Standard_Integer theWLIndex,
+                                                const stCoeffsValue& theCoeffs,
+                                                Standard_Real& theU2,
+                                                Standard_Real& theV1,
+                                                Standard_Real& theV2)
+{
+  if(!CylCylComputeParameters(theU1par, theWLIndex, theCoeffs, theU2))
+    return Standard_False;
 
-  ShortCosForm(aB2,aB1,mB,mFI1);
-  ShortCosForm(aA2,aA1,aA,mFI2);
+  if(!CylCylComputeParameters(theU1par, theU2, theCoeffs, theV1, theV2))
+    return Standard_False;
 
-  mB /= aA;
-  mC /= aA;
+  return Standard_True;
 }
 
 //=======================================================================
 //function : SearchOnVBounds
 //purpose  : 
 //=======================================================================
-static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
-                                        const stCoeffsValue& theCoeffs,
+Standard_Boolean WorkWithBoundaries::
+                        SearchOnVBounds(const SearchBoundType theSBType,
                                         const Standard_Real theVzad,
+                                        const Standard_Real theVInit,
                                         const Standard_Real theInitU2,
                                         const Standard_Real theInitMainVar,
-                                        Standard_Real& theMainVariableValue)
+                                        Standard_Real& theMainVariableValue) const
 {
+  const Standard_Integer aNbDim = 3;
   const Standard_Real aMaxError = 4.0*M_PI; // two periods
   
-  theMainVariableValue = RealLast();
-  const Standard_Real aTol2 = Precision::PConfusion()*Precision::PConfusion();
-  Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVzad;
+  theMainVariableValue = theInitMainVar;
+  const Standard_Real aTol2 = 1.0e-18;
+  Standard_Real aMainVarPrev = theInitMainVar, aU2Prev = theInitU2, anOtherVar = theVInit;
   
+  //Structure of aMatr:
+  //  C_{1}*U_{1} & C_{2}*U_{2} & C_{3}*V_{*},
+  //where C_{1}, C_{2} and C_{3} are math_Vector.
+  math_Matrix aMatr(1, aNbDim, 1, aNbDim);
+
   Standard_Real anError = RealLast();
+  Standard_Real anErrorPrev = anError;
   Standard_Integer aNbIter = 0;
   do
   {
     if(++aNbIter > 1000)
       return Standard_False;
 
-    gp_Vec aVecMainVar = theCoeffs.mVecA1 * sin(aMainVarPrev) - theCoeffs.mVecB1 * cos(aMainVarPrev);
-    gp_Vec aVecFreeMem =  (theCoeffs.mVecA2 * aU2Prev + theCoeffs.mVecB2) * sin(aU2Prev) -
-                          (theCoeffs.mVecB2 * aU2Prev - theCoeffs.mVecA2) * cos(aU2Prev) +
-                          (theCoeffs.mVecA1 * aMainVarPrev + theCoeffs.mVecB1) * sin(aMainVarPrev) -
-                          (theCoeffs.mVecB1 * aMainVarPrev - theCoeffs.mVecA1) * cos(aMainVarPrev) + theCoeffs.mVecD;
+    const Standard_Real aSinU1 = sin(aMainVarPrev),
+                        aCosU1 = cos(aMainVarPrev),
+                        aSinU2 = sin(aU2Prev),
+                        aCosU2 = cos(aU2Prev);
+
+    math_Vector aVecFreeMem = (myCoeffs.mVecA2 * aU2Prev +
+                                              myCoeffs.mVecB2) * aSinU2 -
+                              (myCoeffs.mVecB2 * aU2Prev -
+                                              myCoeffs.mVecA2) * aCosU2 +
+                              (myCoeffs.mVecA1 * aMainVarPrev +
+                                              myCoeffs.mVecB1) * aSinU1 -
+                              (myCoeffs.mVecB1 * aMainVarPrev -
+                                              myCoeffs.mVecA1) * aCosU1 +
+                                                            myCoeffs.mVecD;
 
-    gp_Vec aVecVar1 = theCoeffs.mVecA2 * sin(aU2Prev) - theCoeffs.mVecB2 * cos(aU2Prev);
-    gp_Vec aVecVar2;
+    math_Vector aMSum(1, 3);
 
     switch(theSBType)
     {
     case SearchV1:
-      aVecVar2 = theCoeffs.mVecC2;
-      aVecFreeMem -= theCoeffs.mVecC1 * theVzad;
+      aMatr.SetCol(3, myCoeffs.mVecC2);
+      aMSum = myCoeffs.mVecC1 * theVzad;
+      aVecFreeMem -= aMSum;
+      aMSum += myCoeffs.mVecC2*anOtherVar;
       break;
 
     case SearchV2:
-      aVecVar2 = theCoeffs.mVecC1;
-      aVecFreeMem -= theCoeffs.mVecC2 * theVzad;
+      aMatr.SetCol(3, myCoeffs.mVecC1);
+      aMSum = myCoeffs.mVecC2 * theVzad;
+      aVecFreeMem -= aMSum;
+      aMSum += myCoeffs.mVecC1*anOtherVar;
       break;
 
     default:
       return Standard_False;
     }
 
-    Standard_Real aDetMainSyst =  aVecMainVar.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())-
-                                  aVecMainVar.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+
-                                  aVecMainVar.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X());
+    aMatr.SetCol(1, myCoeffs.mVecA1 * aSinU1 - myCoeffs.mVecB1 * aCosU1);
+    aMatr.SetCol(2, myCoeffs.mVecA2 * aSinU2 - myCoeffs.mVecB2 * aCosU2);
 
-    if(IsEqual(aDetMainSyst, 0.0))
+    Standard_Real aDetMainSyst = aMatr.Determinant();
+
+    if(Abs(aDetMainSyst) < aNulValue)
     {
       return Standard_False;
     }
 
-  
-    Standard_Real aDetMainVar = aVecFreeMem.X()*(aVecVar1.Y()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.Y())-
-                                aVecFreeMem.Y()*(aVecVar1.X()*aVecVar2.Z()-aVecVar1.Z()*aVecVar2.X())+
-                                aVecFreeMem.Z()*(aVecVar1.X()*aVecVar2.Y()-aVecVar1.Y()*aVecVar2.X());
-
-    Standard_Real aDetVar1    = aVecMainVar.X()*(aVecFreeMem.Y()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.Y())-
-                                aVecMainVar.Y()*(aVecFreeMem.X()*aVecVar2.Z()-aVecFreeMem.Z()*aVecVar2.X())+
-                                aVecMainVar.Z()*(aVecFreeMem.X()*aVecVar2.Y()-aVecFreeMem.Y()*aVecVar2.X());
+    math_Matrix aM1(aMatr), aM2(aMatr), aM3(aMatr);
+    aM1.SetCol(1, aVecFreeMem);
+    aM2.SetCol(2, aVecFreeMem);
+    aM3.SetCol(3, aVecFreeMem);
 
-    Standard_Real aDetVar2    = aVecMainVar.X()*(aVecVar1.Y()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.Y())-
-                                aVecMainVar.Y()*(aVecVar1.X()*aVecFreeMem.Z()-aVecVar1.Z()*aVecFreeMem.X())+
-                                aVecMainVar.Z()*(aVecVar1.X()*aVecFreeMem.Y()-aVecVar1.Y()*aVecFreeMem.X());
+    const Standard_Real aDetMainVar = aM1.Determinant();
+    const Standard_Real aDetVar1    = aM2.Determinant();
+    const Standard_Real aDetVar2    = aM3.Determinant();
 
     Standard_Real aDelta = aDetMainVar/aDetMainSyst-aMainVarPrev;
 
@@ -1009,6 +1644,27 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
     aDelta = aDetVar2/aDetMainSyst-anOtherVar;
     anError += aDelta*aDelta;
     anOtherVar += aDelta;
+
+    if(anError > anErrorPrev)
+    {//Method diverges. Keep the best result
+      const Standard_Real aSinU1Last = sin(aMainVarPrev),
+                          aCosU1Last = cos(aMainVarPrev),
+                          aSinU2Last = sin(aU2Prev),
+                          aCosU2Last = cos(aU2Prev);
+      aMSum -= (myCoeffs.mVecA1*aCosU1Last + 
+                myCoeffs.mVecB1*aSinU1Last +
+                myCoeffs.mVecA2*aCosU2Last +
+                myCoeffs.mVecB2*aSinU2Last +
+                myCoeffs.mVecD);
+      const Standard_Real aSQNorm = aMSum.Norm2();
+      return (aSQNorm < aTol2);
+    }
+    else
+    {
+      theMainVariableValue = aMainVarPrev;
+    }
+
+    anErrorPrev = anError;
   }
   while(anError > aTol2);
 
@@ -1019,15 +1675,23 @@ static Standard_Boolean SearchOnVBounds(const SearchBoundType theSBType,
 
 //=======================================================================
 //function : InscribePoint
-//purpose  : If theFlForce==TRUE theUGiven will be adjasted forceful.
+//purpose  : If theFlForce==TRUE theUGiven will be changed forcefully
+//            even if theUGiven is already inscibed in the boundary
+//            (if it is possible; i.e. if new theUGiven is inscribed
+//            in the boundary, too).
 //=======================================================================
-static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
-                                      const Standard_Real theUlTarget,
-                                      Standard_Real& theUGiven,
-                                      const Standard_Real theTol2D,
-                                      const Standard_Real thePeriod,
-                                      const Standard_Boolean theFlForce)
+Standard_Boolean InscribePoint( const Standard_Real theUfTarget,
+                                const Standard_Real theUlTarget,
+                                Standard_Real& theUGiven,
+                                const Standard_Real theTol2D,
+                                const Standard_Real thePeriod,
+                                const Standard_Boolean theFlForce)
 {
+  if(Precision::IsInfinite(theUGiven))
+  {
+    return Standard_False;
+  }
+
   if((theUfTarget - theUGiven <= theTol2D) &&
               (theUGiven - theUlTarget <= theTol2D))
   {//it has already inscribed
@@ -1058,53 +1722,53 @@ static Standard_Boolean InscribePoint(const Standard_Real theUfTarget,
     return Standard_True;
   }
 
-  if(IsEqual(thePeriod, 0.0))
-  {//it cannot be inscribed
-    return Standard_False;
-  }
-
-  Standard_Real aDU = theUGiven - theUfTarget;
-
-  if(aDU > 0.0)
-    aDU = -thePeriod;
-  else
-    aDU = +thePeriod;
+  const Standard_Real aUf = theUfTarget - theTol2D;
+  const Standard_Real aUl = aUf + thePeriod;
 
-  while(((theUGiven - theUfTarget)*aDU < 0.0) && !((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D)))
-  {
-    theUGiven += aDU;
-  }
+  theUGiven = ElCLib::InPeriod(theUGiven, aUf, aUl);
   
-  return ((theUfTarget - theUGiven <= theTol2D) && (theUGiven - theUlTarget <= theTol2D));
+  return ((theUfTarget - theUGiven <= theTol2D) &&
+          (theUGiven - theUlTarget <= theTol2D));
 }
 
 //=======================================================================
 //function : InscribeInterval
-//purpose  : Adjusts theUfGiven and after that fits theUlGiven to result
+//purpose  : Shifts theRange in order to make at least one of its 
+//            boundary in the range [theUfTarget, theUlTarget]
 //=======================================================================
 static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
-                                      const Standard_Real theUlTarget,
-                                      Standard_Real& theUfGiven,
-                                      Standard_Real& theUlGiven,
-                                      const Standard_Real theTol2D,
-                                      const Standard_Real thePeriod)
+                                         const Standard_Real theUlTarget,
+                                         Bnd_Range &theRange,
+                                         const Standard_Real theTol2D,
+                                         const Standard_Real thePeriod)
 {
-  Standard_Real anUpar = theUfGiven;
-  const Standard_Real aDelta = theUlGiven - theUfGiven;
-  if(InscribePoint(theUfTarget, theUlTarget, anUpar,
-                        theTol2D, thePeriod, Standard_False))
+  Standard_Real anUpar = 0.0;
+  if (!theRange.GetMin(anUpar))
+  {
+    return Standard_False;
+  }
+
+  const Standard_Real aDelta = theRange.Delta();
+  if(InscribePoint(theUfTarget, theUlTarget, anUpar, 
+          theTol2D, thePeriod, (Abs(theUlTarget-anUpar) < theTol2D)))
   {
-    theUfGiven = anUpar;
-    theUlGiven = theUfGiven + aDelta;
+    theRange.SetVoid();
+    theRange.Add(anUpar);
+    theRange.Add(anUpar + aDelta);
   }
   else 
   {
-    anUpar = theUlGiven;
+    if (!theRange.GetMax (anUpar))
+    {
+      return Standard_False;
+    }
+
     if(InscribePoint(theUfTarget, theUlTarget, anUpar,
-                        theTol2D, thePeriod, Standard_False))
+          theTol2D, thePeriod, (Abs(theUfTarget-anUpar) < theTol2D)))
     {
-      theUlGiven = anUpar;
-      theUfGiven = theUlGiven - aDelta;
+      theRange.SetVoid();
+      theRange.Add(anUpar);
+      theRange.Add(anUpar - aDelta);
     }
     else
     {
@@ -1116,83 +1780,151 @@ static Standard_Boolean InscribeInterval(const Standard_Real theUfTarget,
 }
 
 //=======================================================================
-//function : InscribeAndSortArray
-//purpose  : Sort from Min to Max value
+//function : ExcludeNearElements
+//purpose  : Checks if theArr contains two almost equal elements.
+//            If it is true then one of equal elements will be excluded
+//            (made infinite).
+//           Returns TRUE if at least one element of theArr has been changed.
+//ATTENTION!!!
+//           1. Every not infinite element of theArr is considered to be 
+//            in [0, T] interval (where T is the period);
+//           2. theArr must be sorted in ascending order.
 //=======================================================================
-static void InscribeAndSortArray( Standard_Real theArr[],
-                                  const Standard_Integer theNOfMembers,
-                                  const Standard_Real theUf,
-                                  const Standard_Real theUl,
-                                  const Standard_Real theTol2D,
-                                  const Standard_Real thePeriod)
+static Standard_Boolean ExcludeNearElements(Standard_Real theArr[],
+                                            const Standard_Integer theNOfMembers,
+                                            const Standard_Real theUSurf1f,
+                                            const Standard_Real theUSurf1l,
+                                            const Standard_Real theTol)
 {
-  for(Standard_Integer i = 0; i < theNOfMembers; i++)
+  Standard_Boolean aRetVal = Standard_False;
+  for(Standard_Integer i = 1; i < theNOfMembers; i++)
   {
-    if(Precision::IsInfinite(theArr[i]))
-    {
-      if(theArr[i] < 0.0)
-        theArr[i] = -theArr[i];
+    Standard_Real &anA = theArr[i],
+                  &anB = theArr[i-1];
 
-      continue;
-    }
+    //Here, anA >= anB
 
-    InscribePoint(theUf, theUl, theArr[i], theTol2D, thePeriod, Standard_False);
+    if(Precision::IsInfinite(anA))
+      break;
 
-    for(Standard_Integer j = i - 1; j >= 0; j--)
+    if((anA-anB) < theTol)
     {
+      if((anB != 0.0) && (anB != theUSurf1f) && (anB != theUSurf1l)) 
+      anA = (anA + anB)/2.0;
+      else
+        anA = anB;
 
-      if(theArr[j+1] < theArr[j])
-      {
-        Standard_Real anUtemp = theArr[j+1];
-        theArr[j+1] = theArr[j];
-        theArr[j] = anUtemp;
-      }
+      //Make this element infinite an forget it
+      //(we will not use it in next iterations).
+      anB = Precision::Infinite();
+      aRetVal = Standard_True;
     }
   }
-}
 
+  return aRetVal;
+}
 
 //=======================================================================
 //function : AddPointIntoWL
 //purpose  : Surf1 is a surface, whose U-par is variable.
+//           If theFlBefore == TRUE then we enable the U1-parameter
+//            of the added point to be less than U1-parameter of
+//           previously added point (in general U1-parameter is
+//           always increased; therefore, this situation is abnormal).
+//           If theOnlyCheck==TRUE then no point will be added to theLine.
 //=======================================================================
 static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
                                         const IntSurf_Quadric& theQuad2,
+                                        const ComputationMethods::stCoeffsValue& theCoeffs,
                                         const Standard_Boolean isTheReverse,
+                                        const Standard_Boolean isThePrecise,
                                         const gp_Pnt2d& thePntOnSurf1,
                                         const gp_Pnt2d& thePntOnSurf2,
                                         const Standard_Real theUfSurf1,
                                         const Standard_Real theUlSurf1,
+                                        const Standard_Real theUfSurf2,
+                                        const Standard_Real theUlSurf2,
+                                        const Standard_Real theVfSurf1,
+                                        const Standard_Real theVlSurf1,
+                                        const Standard_Real theVfSurf2,
+                                        const Standard_Real theVlSurf2,
                                         const Standard_Real thePeriodOfSurf1,
                                         const Handle(IntSurf_LineOn2S)& theLine,
+                                        const Standard_Integer theWLIndex,
                                         const Standard_Real theTol3D,
                                         const Standard_Real theTol2D,
-                                        const Standard_Boolean theFlForce)
+                                        const Standard_Boolean theFlBefore = Standard_False,
+                                        const Standard_Boolean theOnlyCheck = Standard_False)
 {
+  //Check if the point is in the domain or can be inscribed in the domain after adjusting.
+  
   gp_Pnt  aPt1(theQuad1.Value(thePntOnSurf1.X(), thePntOnSurf1.Y())),
           aPt2(theQuad2.Value(thePntOnSurf2.X(), thePntOnSurf2.Y()));
 
-  Standard_Real anUpar = thePntOnSurf1.X();
-  if(!InscribePoint(theUfSurf1, theUlSurf1, anUpar, theTol2D,
-                                    thePeriodOfSurf1, theFlForce))
+  Standard_Real aU1par = thePntOnSurf1.X();
+
+  // aU1par always increases. Therefore, we must reduce its
+  // value in order to continue creation of WLine.
+  if(!InscribePoint(theUfSurf1, theUlSurf1, aU1par, theTol2D,
+                  thePeriodOfSurf1, aU1par > 0.5*(theUfSurf1+theUlSurf1)))
+    return Standard_False;
+
+  if ((theLine->NbPoints() > 0) &&
+      ((theUlSurf1 - theUfSurf1) >= thePeriodOfSurf1) &&      
+      (((aU1par + thePeriodOfSurf1 - theUlSurf1) <= theTol2D) ||
+       ((aU1par - thePeriodOfSurf1 - theUfSurf1) >= theTol2D)))
+  {
+    // aU1par can be adjusted to both theUlSurf1 and theUfSurf1
+    // with equal possibilities. This code fragment allows choosing
+    // correct parameter from these two variants.
+
+    Standard_Real aU1 = 0.0, aV1 = 0.0;
+    if (isTheReverse)
+    {
+      theLine->Value(theLine->NbPoints()).ParametersOnS2(aU1, aV1);
+    }
+    else
+    {
+      theLine->Value(theLine->NbPoints()).ParametersOnS1(aU1, aV1);
+    }
+
+    const Standard_Real aDelta = aU1 - aU1par;
+    if (2.0*Abs(aDelta) > thePeriodOfSurf1)
+    {
+      aU1par += Sign(thePeriodOfSurf1, aDelta);
+    }
+  }
+
+  Standard_Real aU2par = thePntOnSurf2.X();
+  if(!InscribePoint(theUfSurf2, theUlSurf2, aU2par, theTol2D,
+                                    thePeriodOfSurf1, Standard_False))
+    return Standard_False;
+
+  Standard_Real aV1par = thePntOnSurf1.Y();
+  if((aV1par - theVlSurf1 > theTol2D) || (theVfSurf1 - aV1par > theTol2D))
     return Standard_False;
 
+  Standard_Real aV2par = thePntOnSurf2.Y();
+  if((aV2par -  theVlSurf2 > theTol2D) || (theVfSurf2 - aV2par > theTol2D))
+    return Standard_False;
+  
+  //Get intersection point and add it in the WL
   IntSurf_PntOn2S aPnt;
   
   if(isTheReverse)
   {
     aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
-                        thePntOnSurf2.X(), thePntOnSurf2.Y(),
-                        anUpar, thePntOnSurf1.Y());
+                        aU2par, aV2par,
+                        aU1par, aV1par);
   }
   else
   {
     aPnt.SetValue((aPt1.XYZ()+aPt2.XYZ())/2.0,
-                        anUpar, thePntOnSurf1.Y(),
-                        thePntOnSurf2.X(), thePntOnSurf2.Y());
+                        aU1par, aV1par,
+                        aU2par, aV2par);
   }
 
-  const Standard_Integer aNbPnts = theLine->NbPoints();
+  Standard_Integer aNbPnts = theLine->NbPoints();
   if(aNbPnts > 0)
   {
     Standard_Real aUl = 0.0, aVl = 0.0;
@@ -1202,11 +1934,27 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
     else
       aPlast.ParametersOnS1(aUl, aVl);
 
-    if(anUpar <= aUl)
-    {//Parameter value will be always increased.
-      return Standard_False;
+    if(!theFlBefore && (aU1par <= aUl))
+    {
+      //Parameter value must be increased if theFlBefore == FALSE.
+
+      aU1par += thePeriodOfSurf1;
+
+      //The condition is as same as in
+      //InscribePoint(...) function
+      if((theUfSurf1 - aU1par > theTol2D) ||
+         (aU1par - theUlSurf1 > theTol2D))
+      {
+        //New aU1par is out of target interval.
+        //Go back to old value.
+
+        return Standard_False;
+      }
     }
 
+    if (theOnlyCheck)
+      return Standard_True;
+
     //theTol2D is minimal step along parameter changed. 
     //Therefore, if we apply this minimal step two 
     //neighbour points will be always "same". Consequently,
@@ -1218,35 +1966,65 @@ static Standard_Boolean AddPointIntoWL( const IntSurf_Quadric& theQuad1,
     }
   }
 
+  if (theOnlyCheck)
+    return Standard_True;
+
   theLine->Add(aPnt);
+
+  if(!isThePrecise)
+    return Standard_True;
+
+  //Try to precise existing WLine
+  aNbPnts = theLine->NbPoints();
+  if(aNbPnts >= 3)
+  {
+    Standard_Real aU1 = 0.0, aU2 = 0.0, aU3 = 0.0, aV = 0.0;
+    if(isTheReverse)
+    {
+      theLine->Value(aNbPnts).ParametersOnS2(aU3, aV);
+      theLine->Value(aNbPnts-1).ParametersOnS2(aU2, aV);
+      theLine->Value(aNbPnts-2).ParametersOnS2(aU1, aV);
+    }
+    else
+    {
+      theLine->Value(aNbPnts).ParametersOnS1(aU3, aV);
+      theLine->Value(aNbPnts-1).ParametersOnS1(aU2, aV);
+      theLine->Value(aNbPnts-2).ParametersOnS1(aU1, aV);
+    }
+
+    const Standard_Real aStepPrev = aU2-aU1;
+    const Standard_Real aStep = aU3-aU2;
+
+    const Standard_Integer aDeltaStep = RealToInt(aStepPrev/aStep);
+
+    if((1 < aDeltaStep) && (aDeltaStep < 2000))
+    {
+      //Add new points in case of non-uniform distribution of existing points
+      SeekAdditionalPoints( theQuad1, theQuad2, theLine, 
+                            theCoeffs, theWLIndex, aDeltaStep, aNbPnts-2,
+                            aNbPnts-1, theTol2D, thePeriodOfSurf1, isTheReverse);
+    }
+  }
+
   return Standard_True;
 }
 
 //=======================================================================
 //function : AddBoundaryPoint
-//purpose  : 
+//purpose  : Find intersection point on V-boundary.
 //=======================================================================
-static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
-                                          const IntSurf_Quadric& theQuad2,
-                                          const Handle(IntPatch_WLine)& theWL,
-                                          const stCoeffsValue& theCoeffs,
-                                          const Bnd_Box2d& theUVSurf1,
-                                          const Bnd_Box2d& theUVSurf2,
-                                          const Standard_Real theTol3D,
-                                          const Standard_Real theTol2D,
-                                          const Standard_Real thePeriod,
-                                          const Standard_Real theNulValue,
+void WorkWithBoundaries::AddBoundaryPoint(const Handle(IntPatch_WLine)& theWL,
                                           const Standard_Real theU1,
+                                          const Standard_Real theU1Min,
                                           const Standard_Real theU2,
                                           const Standard_Real theV1,
                                           const Standard_Real theV1Prev,
                                           const Standard_Real theV2,
                                           const Standard_Real theV2Prev,
-                                          const Standard_Boolean isTheReverse,
-                                          const Standard_Real theArccosFactor,
+                                          const Standard_Integer theWLIndex,
                                           const Standard_Boolean theFlForce,
                                           Standard_Boolean& isTheFound1,
-                                          Standard_Boolean& isTheFound2)
+                                          Standard_Boolean& isTheFound2) const
 {
   Standard_Real aUSurf1f = 0.0, //const
                 aUSurf1l = 0.0,
@@ -1257,208 +2035,167 @@ static Standard_Boolean AddBoundaryPoint( const IntSurf_Quadric& theQuad1,
                 aVSurf2f = 0.0,
                 aVSurf2l = 0.0;
 
-  theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
-  theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
-
-  SearchBoundType aTS1 = SearchNONE, aTS2 = SearchNONE;
-  Standard_Real aV1zad = aVSurf1f, aV2zad = aVSurf2f;
-
-  Standard_Real anUpar1 = theU1, anUpar2 = theU1;
-  Standard_Real aVf = theV1, aVl = theV1Prev;
-  MinMax(aVf, aVl);
-  if((aVf <= aVSurf1f) && (aVSurf1f <= aVl))
-  {
-    aTS1 = SearchV1;
-    aV1zad = aVSurf1f;
-    isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1f, theU2, theU1, anUpar1);
-  }
-  else if((aVf <= aVSurf1l) && (aVSurf1l <= aVl))
-  {
-    aTS1 = SearchV1;
-    aV1zad = aVSurf1l;
-    isTheFound1 = SearchOnVBounds(aTS1, theCoeffs, aVSurf1l, theU2, theU1, anUpar1);
-  }
-
-  aVf = theV2;
-  aVl = theV2Prev;
-  MinMax(aVf, aVl);
-
-  if((aVf <= aVSurf2f) && (aVSurf2f <= aVl))
-  {
-    aTS2 = SearchV2;
-    aV2zad = aVSurf2f;
-    isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2f, theU2, theU1, anUpar2);
-  }
-  else if((aVf <= aVSurf2l) && (aVSurf2l <= aVl))
-  {
-    aTS2 = SearchV2;
-    aV2zad = aVSurf2l;
-    isTheFound2 = SearchOnVBounds(aTS2, theCoeffs, aVSurf2l, theU2, theU1, anUpar2);
-  }
-
-  if(anUpar1 < anUpar2)
-  {
-    if(isTheFound1)
-    {
-      Standard_Real anArg = theCoeffs.mB * cos(anUpar1 - theCoeffs.mFI1) + theCoeffs.mC;
-
-      if(theNulValue > 1.0 - anArg)
-        anArg = 1.0;
-      if(anArg + 1.0 < theNulValue)
-        anArg = -1.0;
-
-      Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
-      
-      if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
-      {
-        const Standard_Real aV1 = 
-              (aTS1 == SearchV1) ? aV1zad : 
-                                  theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) +
-                                  theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1;
-        const Standard_Real aV2 = 
-              (aTS1 == SearchV2) ? aV2zad : 
-                                  theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) +
-                                  theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2;
-
-        AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
-                          gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
-                          aUSurf1f, aUSurf1l, thePeriod,
-                          theWL->Curve(), theTol3D, theTol2D, theFlForce);
-      }
-      else
-      {
-        isTheFound1 = Standard_False;
-      }
-    }
-
-    if(isTheFound2)
-    {
-      Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC;
+  myUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
+  myUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
+  
+  const Standard_Integer aSize = 4;
+  const Standard_Real anArrVzad[aSize] = {aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l};
 
-      if(theNulValue > 1.0 - anArg)
-        anArg = 1.0;
-      if(anArg + 1.0 < theNulValue)
-        anArg = -1.0;
+  StPInfo aUVPoint[aSize];
 
-      Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
-      if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
-      {
-        const Standard_Real aV1 = 
-              (aTS2 == SearchV1) ? aV1zad : 
-                                  theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) +
-                                  theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1;
-        const Standard_Real aV2 = 
-              (aTS2 == SearchV2) ? aV2zad : 
-                                  theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) +
-                                  theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2;
-
-        AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
-                        gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
-                        aUSurf1f, aUSurf1l, thePeriod,
-                        theWL->Curve(),theTol3D, theTol2D, theFlForce);
-      }
-      else
-      {
-        isTheFound2 = Standard_False;
-      }
-    }
-  }
-  else
+  for(Standard_Integer anIDSurf = 0; anIDSurf < 4; anIDSurf+=2)
   {
-    if(isTheFound2)
+    const Standard_Real aVf = (anIDSurf == 0) ? theV1 : theV2,
+                        aVl = (anIDSurf == 0) ? theV1Prev : theV2Prev;
+
+    const SearchBoundType aTS = (anIDSurf == 0) ? SearchV1 : SearchV2;
+
+    for(Standard_Integer anIDBound = 0; anIDBound < 2; anIDBound++)
     {
-      Standard_Real anArg = theCoeffs.mB * cos(anUpar2 - theCoeffs.mFI1) + theCoeffs.mC;
+      const Standard_Integer anIndex = anIDSurf+anIDBound;
 
-      if(theNulValue > 1.0 - anArg)
-        anArg = 1.0;
-      if(anArg + 1.0 < theNulValue)
-        anArg = -1.0;
+      aUVPoint[anIndex].mySurfID = anIDSurf;
 
-      Standard_Real aU2 = theCoeffs.mFI2 + theArccosFactor * acos(anArg);
-      
-      if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
-      {
-        const Standard_Real aV1 = (aTS2 == SearchV1) ? aV1zad : 
-                                  theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar2) +
-                                  theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar2) + theCoeffs.mM1;
-        const Standard_Real aV2 = (aTS2 == SearchV2) ? aV2zad : 
-                                  theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar2) +
-                                  theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar2) + theCoeffs.mM2;
-
-        AddPointIntoWL(theQuad1, theQuad2, isTheReverse, 
-                        gp_Pnt2d(anUpar2, aV1), gp_Pnt2d(aU2, aV2),
-                        aUSurf1f, aUSurf1l, thePeriod,
-                        theWL->Curve(), theTol3D, theTol2D, theFlForce);
-      }
-      else
+      if((Abs(aVf-anArrVzad[anIndex]) > myTol2D) &&
+          ((aVf-anArrVzad[anIndex])*(aVl-anArrVzad[anIndex]) > 0.0))
       {
-        isTheFound2 = Standard_False;
+        continue;
       }
-    }
 
-    if(isTheFound1)
-    {
-      Standard_Real anArg = theCoeffs.mB*cos(anUpar1-theCoeffs.mFI1)+theCoeffs.mC;
+      //Segment [aVf, aVl] intersects at least one V-boundary (first or last)
+      // (in general, case is possible, when aVf > aVl).
 
-      if(theNulValue > 1.0 - anArg)
-        anArg = 1.0;
-      if(anArg + 1.0 < theNulValue)
-        anArg = -1.0;
+      // Precise intersection point
+      const Standard_Boolean aRes = SearchOnVBounds(aTS, anArrVzad[anIndex],
+                                                    (anIDSurf == 0) ? theV2 : theV1,
+                                                    theU2, theU1,
+                                                    aUVPoint[anIndex].myU1);
 
-      Standard_Real aU2 = theCoeffs.mFI2+theArccosFactor*acos(anArg);
-      if(InscribePoint(aUSurf2f, aUSurf2l, aU2, theTol2D, thePeriod, Standard_False))
+      // aUVPoint[anIndex].myU1 is considered to be nearer to theU1 than
+      // to theU1+/-Period
+      if (!aRes || (aUVPoint[anIndex].myU1 >= theU1) ||
+                              (aUVPoint[anIndex].myU1 < theU1Min))
       {
-        const Standard_Real aV1 = (aTS1 == SearchV1) ? aV1zad :
-                                  theCoeffs.mK21 * sin(aU2) + theCoeffs.mK11 * sin(anUpar1) +
-                                  theCoeffs.mL21 * cos(aU2) + theCoeffs.mL11 * cos(anUpar1) + theCoeffs.mM1;
-        const Standard_Real aV2 = (aTS1 == SearchV2) ? aV2zad : 
-                                  theCoeffs.mK22 * sin(aU2) + theCoeffs.mK12 * sin(anUpar1) +
-                                  theCoeffs.mL22 * cos(aU2) + theCoeffs.mL12 * cos(anUpar1) + theCoeffs.mM2;
-
-        AddPointIntoWL(theQuad1, theQuad2, isTheReverse,
-                        gp_Pnt2d(anUpar1, aV1), gp_Pnt2d(aU2, aV2),
-                        aUSurf1f, aUSurf1l, thePeriod,
-                        theWL->Curve(), theTol3D, theTol2D, theFlForce);
+        //Intersection point is not found or out of the domain
+        aUVPoint[anIndex].myU1 = RealLast();
+        continue;
       }
       else
       {
-        isTheFound1 = Standard_False;
+        //intersection point is found
+
+        Standard_Real &aU1 = aUVPoint[anIndex].myU1,
+                      &aU2 = aUVPoint[anIndex].myU2,
+                      &aV1 = aUVPoint[anIndex].myV1,
+                      &aV2 = aUVPoint[anIndex].myV2;
+        aU2 = theU2;
+        aV1 = theV1;
+        aV2 = theV2;
+
+        if(!ComputationMethods::
+                  CylCylComputeParameters(aU1, theWLIndex, myCoeffs, aU2, aV1, aV2))
+        {
+          // Found point is wrong
+          aU1 = RealLast();
+          continue;
+        }
+
+        //Point on true V-boundary.
+        if(aTS == SearchV1)
+          aV1 = anArrVzad[anIndex];
+        else //if(aTS[anIndex] == SearchV2)
+          aV2 = anArrVzad[anIndex];
       }
     }
   }
 
-  return Standard_True;
+  // Sort with acceding U1-parameter.
+  std::sort(aUVPoint, aUVPoint+aSize);
+    
+  isTheFound1 = isTheFound2 = Standard_False;
+
+  //Adding found points on boundary in the WLine.
+  for(Standard_Integer i = 0; i < aSize; i++)
+  {
+    if(aUVPoint[i].myU1 == RealLast())
+      break;
+
+    if(!AddPointIntoWL(myQuad1, myQuad2, myCoeffs, myIsReverse, Standard_False,
+                        gp_Pnt2d(aUVPoint[i].myU1, aUVPoint[i].myV1),
+                        gp_Pnt2d(aUVPoint[i].myU2, aUVPoint[i].myV2),
+                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, myPeriod,
+                        theWL->Curve(), theWLIndex, myTol3D, myTol2D, theFlForce))
+    {
+      continue;
+    }
+
+    if(aUVPoint[i].mySurfID == 0)
+    {
+      isTheFound1 = Standard_True;
+    }
+    else
+    {
+      isTheFound2 = Standard_True;
+    }
+  }
 }
 
 //=======================================================================
 //function : SeekAdditionalPoints
-//purpose  : 
+//purpose  : Inserts additional intersection points between neighbor points.
+//            This action is bone with several iterations. During every iteration,
+//          new point is inserted in middle of every interval.
+//            The process will be finished if NbPoints >= theMinNbPoints.
 //=======================================================================
 static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
                                   const IntSurf_Quadric& theQuad2,
-                                  const Handle(IntSurf_LineOn2S)& theLile,
-                                  const stCoeffsValue& theCoeffs,
+                                  const Handle(IntSurf_LineOn2S)& theLine,
+                                  const ComputationMethods::stCoeffsValue& theCoeffs,
+                                  const Standard_Integer theWLIndex,
                                   const Standard_Integer theMinNbPoints,
-                                  const Standard_Real theU2f,
-                                  const Standard_Real theU2l,
+                                  const Standard_Integer theStartPointOnLine,
+                                  const Standard_Integer theEndPointOnLine,
                                   const Standard_Real theTol2D,
                                   const Standard_Real thePeriodOfSurf2,
-                                  const Standard_Real theArccosFactor,
                                   const Standard_Boolean isTheReverse)
 {
-  Standard_Integer aNbPoints = theLile->NbPoints();
+  if(theLine.IsNull())
+    return;
+  
+  Standard_Integer aNbPoints = theEndPointOnLine - theStartPointOnLine + 1;
   if(aNbPoints >= theMinNbPoints)
   {
     return;
   }
 
+  Standard_Real aMinDeltaParam = theTol2D;
+
+  {
+    Standard_Real u1 = 0.0, v1 = 0.0, u2 = 0.0, v2 = 0.0;
+
+    if(isTheReverse)
+    {
+      theLine->Value(theStartPointOnLine).ParametersOnS2(u1, v1);
+      theLine->Value(theEndPointOnLine).ParametersOnS2(u2, v2);
+    }
+    else
+    {
+      theLine->Value(theStartPointOnLine).ParametersOnS1(u1, v1);
+      theLine->Value(theEndPointOnLine).ParametersOnS1(u2, v2);
+    }
+    
+    aMinDeltaParam = Max(Abs(u2 - u1)/IntToReal(theMinNbPoints), aMinDeltaParam);
+  }
+
+  Standard_Integer aLastPointIndex = theEndPointOnLine;
   Standard_Real U1prec = 0.0, V1prec = 0.0, U2prec = 0.0, V2prec = 0.0;
 
   Standard_Integer aNbPointsPrev = 0;
   while(aNbPoints < theMinNbPoints && (aNbPoints != aNbPointsPrev))
   {
     aNbPointsPrev = aNbPoints;
-    for(Standard_Integer fp = 1, lp = 2; fp < aNbPoints; fp = lp + 1)
+    for(Standard_Integer fp = theStartPointOnLine, lp = 0; fp < aLastPointIndex; fp = lp + 1)
     {
       Standard_Real U1f = 0.0, V1f = 0.0; //first point in 1st suraface
       Standard_Real U1l = 0.0, V1l = 0.0; //last  point in 1st suraface
@@ -1470,22 +2207,22 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
       
       if(isTheReverse)
       {
-        theLile->Value(fp).ParametersOnS2(U1f, V1f);
-        theLile->Value(lp).ParametersOnS2(U1l, V1l);
+        theLine->Value(fp).ParametersOnS2(U1f, V1f);
+        theLine->Value(lp).ParametersOnS2(U1l, V1l);
 
-        theLile->Value(fp).ParametersOnS1(U2f, V2f);
-        theLile->Value(lp).ParametersOnS1(U2l, V2l);
+        theLine->Value(fp).ParametersOnS1(U2f, V2f);
+        theLine->Value(lp).ParametersOnS1(U2l, V2l);
       }
       else
       {
-        theLile->Value(fp).ParametersOnS1(U1f, V1f);
-        theLile->Value(lp).ParametersOnS1(U1l, V1l);
+        theLine->Value(fp).ParametersOnS1(U1f, V1f);
+        theLine->Value(lp).ParametersOnS1(U1l, V1l);
 
-        theLile->Value(fp).ParametersOnS2(U2f, V2f);
-        theLile->Value(lp).ParametersOnS2(U2l, V2l);
+        theLine->Value(fp).ParametersOnS2(U2f, V2f);
+        theLine->Value(lp).ParametersOnS2(U2l, V2l);
       }
 
-      if(Abs(U1l - U1f) <= theTol2D)
+      if(Abs(U1l - U1f) <= aMinDeltaParam)
       {
         //Step is minimal. It is not necessary to divide it.
         continue;
@@ -1493,33 +2230,23 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
 
       U1prec = 0.5*(U1f+U1l);
       
-      Standard_Real anArg = theCoeffs.mB * cos(U1prec - theCoeffs.mFI1) + theCoeffs.mC;
-      if(anArg > 1.0)
-        anArg = 1.0;
-      if(anArg < -1.0)
-        anArg = -1.0;
-
-      U2prec = theCoeffs.mFI2 + theArccosFactor*acos(anArg);
-      InscribePoint(theU2f, theU2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False);
-
-      gp_Pnt aP1, aP2;
-
-      V1prec = theCoeffs.mK21 * sin(U2prec) + 
-                theCoeffs.mK11 * sin(U1prec) + 
-                theCoeffs.mL21 * cos(U2prec) + 
-                theCoeffs.mL11 * cos(U1prec) + theCoeffs.mM1;
-      V2prec = theCoeffs.mK22 * sin(U2prec) + 
-                theCoeffs.mK12 * sin(U1prec) + 
-                theCoeffs.mL22 * cos(U2prec) + 
-                theCoeffs.mL12 * cos(U1prec) + theCoeffs.mM2;
+      if(!ComputationMethods::
+            CylCylComputeParameters(U1prec, theWLIndex, theCoeffs, U2prec, V1prec, V2prec))
+      {
+        continue;
+      }
 
-      aP1 = theQuad1.Value(U1prec, V1prec);
-      aP2 = theQuad2.Value(U2prec, V2prec);
+      MinMax(U2f, U2l);
+      if(!InscribePoint(U2f, U2l, U2prec, theTol2D, thePeriodOfSurf2, Standard_False))
+      {
+        continue;
+      }
 
-      gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
+      const gp_Pnt aP1(theQuad1.Value(U1prec, V1prec)), aP2(theQuad2.Value(U2prec, V2prec));
+      const gp_Pnt aPInt(0.5*(aP1.XYZ() + aP2.XYZ()));
 
-#ifdef OCCT_DEBUG
-      //cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << endl;
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+      std::cout << "|P1Pi| = " << aP1.SquareDistance(aPInt) << "; |P2Pi| = " << aP2.SquareDistance(aPInt) << std::endl;
 #endif
 
       IntSurf_PntOn2S anIP;
@@ -1532,186 +2259,98 @@ static void SeekAdditionalPoints( const IntSurf_Quadric& theQuad1,
         anIP.SetValue(aPInt, U1prec, V1prec, U2prec, V2prec);
       }
 
-      theLile->InsertBefore(lp, anIP);
+      theLine->InsertBefore(lp, anIP);
 
-      aNbPoints = theLile->NbPoints();
-      if(aNbPoints >= theMinNbPoints)
-      {
-        return;
-      }
+      aNbPoints++;
+      aLastPointIndex++;
     }
-  }
-}
-
-//=======================================================================
-//function : CriticalPointsComputing
-//purpose  : 
-//=======================================================================
-static void CriticalPointsComputing(const stCoeffsValue& theCoeffs,
-                                    const Standard_Real theUSurf1f,
-                                    const Standard_Real theUSurf1l,
-                                    const Standard_Real theUSurf2f,
-                                    const Standard_Real theUSurf2l,
-                                    const Standard_Real thePeriod,
-                                    const Standard_Real theTol2D,
-                                    const Standard_Integer theNbCritPointsMax,
-                                    Standard_Real theU1crit[])
-{
-  theU1crit[0] = 0.0;
-  theU1crit[1] = thePeriod;
-  theU1crit[2] = theUSurf1f;
-  theU1crit[3] = theUSurf1l;
-
-  const Standard_Real aCOS = cos(theCoeffs.mFI2);
-  const Standard_Real aBSB = Abs(theCoeffs.mB);
-  if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
-  {
-    Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
-    if(anArg > 1.0)
-      anArg = 1.0;
-    if(anArg < -1.0)
-      anArg = -1.0;
-
-    theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
-    theU1crit[5] = acos(anArg) + theCoeffs.mFI1;
-  }
-
-  Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
-  Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
-  MinMax(aSf, aSl);
-
-  theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
-  theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ? -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
-  theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?  acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : -Precision::Infinite();
-  theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?  acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 : Precision::Infinite();
 
-  //preparative treatment of array
-  InscribeAndSortArray(theU1crit, theNbCritPointsMax, 0.0, thePeriod, theTol2D, thePeriod);
-  for(Standard_Integer i = 1; i < theNbCritPointsMax; i++)
-  {
-    Standard_Real &a = theU1crit[i],
-                  &b = theU1crit[i-1];
-    const Standard_Real aRemain = fmod(Abs(a - b), thePeriod); // >= 0, because Abs(a - b) >= 0
-    if((Abs(a - b) < theTol2D) || (aRemain < theTol2D) || (Abs(aRemain - thePeriod) < theTol2D))
+    if(aNbPoints >= theMinNbPoints)
     {
-      a = (a + b)/2.0;
-      b = Precision::Infinite();
+      return;
     }
   }
 }
 
 //=======================================================================
-//function : IntCyCyTrim
-//purpose  : 
+//function : BoundariesComputing
+//purpose  : Computes true domain of future intersection curve (allows
+//            avoiding excess iterations)
 //=======================================================================
-Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
-                              const IntSurf_Quadric& theQuad2,
-                              const Standard_Real theTol3D,
-                              const Standard_Real theTol2D,
-                              const Bnd_Box2d& theUVSurf1,
-                              const Bnd_Box2d& theUVSurf2,
-                              const Standard_Boolean isTheReverse,
-                              Standard_Boolean& isTheEmpty,
-                              IntPatch_SequenceOfLine& theSlin,
-                              IntPatch_SequenceOfPoint& theSPnt)
+Standard_Boolean WorkWithBoundaries::
+            BoundariesComputing(const ComputationMethods::stCoeffsValue &theCoeffs,
+                                const Standard_Real thePeriod,
+                                Bnd_Range theURange[])
 {
-  Standard_Real aUSurf1f = 0.0, //const
-                aUSurf1l = 0.0,
-                aVSurf1f = 0.0,
-                aVSurf1l = 0.0;
-  Standard_Real aUSurf2f = 0.0, //const
-                aUSurf2l = 0.0,
-                aVSurf2f = 0.0,
-                aVSurf2l = 0.0;
-
-  theUVSurf1.Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
-  theUVSurf2.Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
-
-  const Standard_Real aNulValue = 0.01*Precision::PConfusion();
-
-  const gp_Cylinder&  aCyl1 = theQuad1.Cylinder(),
-                      aCyl2 = theQuad2.Cylinder();
-
-  IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
+  //All comments to this method is related to the comment
+  //to ComputationMethods class
+  
+  //So, we have the equation
+  //    cos(U2-FI2)=B*cos(U1-FI1)+C
+  //Evidently,
+  //    -1 <= B*cos(U1-FI1)+C <= 1
 
-  if (!anInter.IsDone())
+  if (theCoeffs.mB > 0.0)
   {
-    return Standard_False;
-  }
-
-  IntAna_ResultType aTypInt = anInter.TypeInter();
-
-  if(aTypInt != IntAna_NoGeometricSolution)
-  { //It is not necessary (because result is an analytic curve) or
-    //it is impossible to make Walking-line.
-
-    return Standard_False;
-  }
-    
-  theSlin.Clear();
-  theSPnt.Clear();
-  const Standard_Integer aNbMaxPoints = 2000;
-  const Standard_Integer aNbMinPoints = 200;
-  const Standard_Integer aNbPoints = Min(Max(aNbMinPoints, 
-                      RealToInt(20.0*aCyl1.Radius())), aNbMaxPoints);
-  const Standard_Real aPeriod = 2.0*M_PI;
-  const Standard_Real aStepMin = theTol2D, 
-                      aStepMax = (aUSurf1l-aUSurf1f)/IntToReal(aNbPoints);
-  const Standard_Integer aNbWLines = 2;
-
-  const stCoeffsValue anEquationCoeffs(aCyl1, aCyl2);
+    // -(1+C)/B <= cos(U1-FI1) <= (1-C)/B
 
-  //Boundaries
-  const Standard_Integer aNbOfBoundaries = 2;
-  Standard_Real aU1f[aNbOfBoundaries] = {-Precision::Infinite(), -Precision::Infinite()};
-  Standard_Real aU1l[aNbOfBoundaries] = {Precision::Infinite(), Precision::Infinite()};
-
-  if(anEquationCoeffs.mB > 0.0)
-  {
-    if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) < -1.0)
-    {//There is NOT intersection
-      return Standard_True;
+    if (theCoeffs.mB + Abs(theCoeffs.mC) < -1.0)
+    {
+      //(1-C)/B < -1 or -(1+C)/B > 1  ==> No solution
+      
+      return Standard_False;
     }
-    else if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
-    {//U=[0;2*PI]+aFI1
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
+    else if (theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+    {
+      //(1-C)/B >= 1 and -(1+C)/B <= -1 ==> U=[0;2*PI]+aFI1
+      theURange[0].Add(theCoeffs.mFI1);
+      theURange[0].Add(thePeriod + theCoeffs.mFI1);
     }
-    else if((1 + anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
-            (anEquationCoeffs.mB <= 1 - anEquationCoeffs.mC))
+    else if ((1 + theCoeffs.mC <= theCoeffs.mB) &&
+             (theCoeffs.mB <= 1 - theCoeffs.mC))
     {
-      Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
+      //(1-C)/B >= 1 and -(1+C)/B >= -1 ==> 
+      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
+      //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB)
+
+      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
         anArg = -1.0;
 
       const Standard_Real aDAngle = acos(anArg);
-      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
+      theURange[0].Add(theCoeffs.mFI1);
+      theURange[0].Add(aDAngle + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod + theCoeffs.mFI1);
     }
-    else if((1 - anEquationCoeffs.mC <= anEquationCoeffs.mB) &&
-            (anEquationCoeffs.mB <= 1 + anEquationCoeffs.mC))
+    else if ((1 - theCoeffs.mC <= theCoeffs.mB) &&
+             (theCoeffs.mB <= 1 + theCoeffs.mC))
     {
-      Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
+      //(1-C)/B <= 1 and -(1+C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1
+      //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
+
+      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
         anArg = -1.0;
 
       const Standard_Real aDAngle = acos(anArg);
-      //U=[aDAngle;2*PI-aDAngle]+aFI1
-
-      aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
+      theURange[0].Add(aDAngle + theCoeffs.mFI1);
+      theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
     }
-    else if(anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
+    else if (theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
     {
-      Standard_Real anArg1 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB,
-                    anArg2 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
+      //(1-C)/B <= 1 and -(1+C)/B >= -1 ==>
+      //(U=[aDAngle1;aDAngle2]+aFI1) ||
+      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
+      //where aDAngle1 = acos((1 - myCoeffs.mC) / myCoeffs.mB),
+      //      aDAngle2 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
+
+      Standard_Real anArg1 = (1 - theCoeffs.mC) / theCoeffs.mB,
+                    anArg2 = -(theCoeffs.mC + 1) / theCoeffs.mB;
       if(anArg1 > 1.0)
         anArg1 = 1.0;
       if(anArg1 < -1.0)
@@ -1725,286 +2364,606 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
       const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
       //(U=[aDAngle1;aDAngle2]+aFI1) ||
       //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
-
-      aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
+      theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
+      theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
     }
     else
     {
-      Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
+      return Standard_False;
     }
   }
-  else if(anEquationCoeffs.mB < 0.0)
+  else if (theCoeffs.mB < 0.0)
   {
-    if(anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) > 1.0)
-    {//There is NOT intersection
-      return Standard_True;
+    // (1-C)/B <= cos(U1-FI1) <= -(1+C)/B
+
+    if (theCoeffs.mB + Abs(theCoeffs.mC) > 1.0)
+    {
+      // -(1+C)/B < -1 or (1-C)/B > 1 ==> No solutions
+      return Standard_False;
     }
-    else if(-anEquationCoeffs.mB + Abs(anEquationCoeffs.mC) <= 1.0)
-    {//U=[0;2*PI]+aFI1
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod + anEquationCoeffs.mFI1;
+    else if (-theCoeffs.mB + Abs(theCoeffs.mC) <= 1.0)
+    {
+      //  -(1+C)/B >= 1 and (1-C)/B <= -1 ==> U=[0;2*PI]+aFI1
+      theURange[0].Add(theCoeffs.mFI1);
+      theURange[0].Add(thePeriod + theCoeffs.mFI1);
+    }
+    else if ((-theCoeffs.mC - 1 <= theCoeffs.mB) &&
+             (theCoeffs.mB <= theCoeffs.mC - 1))
+    {
+      //  -(1+C)/B >= 1 and (1-C)/B >= -1 ==> 
+      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1),
+      //where aDAngle = acos((1 - myCoeffs.mC) / myCoeffs.mB)
+
+      Standard_Real anArg = (1 - theCoeffs.mC) / theCoeffs.mB;
+      if(anArg > 1.0)
+        anArg = 1.0;
+      if(anArg < -1.0)
+        anArg = -1.0;
+
+      const Standard_Real aDAngle = acos(anArg);
+      theURange[0].Add(theCoeffs.mFI1);
+      theURange[0].Add(aDAngle + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod - aDAngle + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod + theCoeffs.mFI1);
     }
-    else if((-anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
-            ( anEquationCoeffs.mB <= anEquationCoeffs.mC - 1))
+    else if ((theCoeffs.mC - 1 <= theCoeffs.mB) &&
+             (theCoeffs.mB <= -theCoeffs.mB - 1))
     {
-      Standard_Real anArg = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
+      //  -(1+C)/B <= 1 and (1-C)/B <= -1 ==> U=[aDAngle;2*PI-aDAngle]+aFI1,
+      //where aDAngle = acos(-(myCoeffs.mC + 1) / myCoeffs.mB).
+
+      Standard_Real anArg = -(theCoeffs.mC + 1) / theCoeffs.mB;
       if(anArg > 1.0)
         anArg = 1.0;
       if(anArg < -1.0)
         anArg = -1.0;
 
-      const Standard_Real aDAngle = acos(anArg);
-      //(U=[0;aDAngle]+aFI1) || (U=[2*PI-aDAngle;2*PI]+aFI1)
+      const Standard_Real aDAngle = acos(anArg);
+      theURange[0].Add(aDAngle + theCoeffs.mFI1);
+      theURange[0].Add(thePeriod - aDAngle + theCoeffs.mFI1);
+    }
+    else if (-theCoeffs.mB - Abs(theCoeffs.mC) >= 1.0)
+    {
+      //  -(1+C)/B <= 1 and (1-C)/B >= -1 ==>
+      //(U=[aDAngle1;aDAngle2]+aFI1) || (U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1),
+      //where aDAngle1 = acos(-(myCoeffs.mC + 1) / myCoeffs.mB),
+      //      aDAngle2 = acos((1 - myCoeffs.mC) / myCoeffs.mB)
+
+      Standard_Real anArg1 = -(theCoeffs.mC + 1) / theCoeffs.mB,
+                    anArg2 = (1 - theCoeffs.mC) / theCoeffs.mB;
+      if(anArg1 > 1.0)
+        anArg1 = 1.0;
+      if(anArg1 < -1.0)
+        anArg1 = -1.0;
+
+      if(anArg2 > 1.0)
+        anArg2 = 1.0;
+      if(anArg2 < -1.0)
+        anArg2 = -1.0;
+
+      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
+      theURange[0].Add(aDAngle1 + theCoeffs.mFI1);
+      theURange[0].Add(aDAngle2 + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod - aDAngle2 + theCoeffs.mFI1);
+      theURange[1].Add(thePeriod - aDAngle1 + theCoeffs.mFI1);
+    }
+    else
+    {
+      return Standard_False;
+    }
+  }
+  else
+  {
+    return Standard_False;
+  }
+
+  return Standard_True;
+}
+
+//=======================================================================
+//function : CriticalPointsComputing
+//purpose  : theNbCritPointsMax contains true number of critical points.
+//            It must be initialized correctly before function calling
+//=======================================================================
+static void CriticalPointsComputing(const ComputationMethods::stCoeffsValue& theCoeffs,
+                                    const Standard_Real theUSurf1f,
+                                    const Standard_Real theUSurf1l,
+                                    const Standard_Real theUSurf2f,
+                                    const Standard_Real theUSurf2l,
+                                    const Standard_Real thePeriod,
+                                    const Standard_Real theTol2D,
+                                    Standard_Integer& theNbCritPointsMax,
+                                    Standard_Real theU1crit[])
+{
+  //[0...1] - in these points parameter U1 goes through
+  //          the seam-edge of the first cylinder.
+  //[2...3] - First and last U1 parameter.
+  //[4...5] - in these points parameter U2 goes through
+  //          the seam-edge of the second cylinder.
+  //[6...9] - in these points an intersection line goes through
+  //          U-boundaries of the second surface.
+  //[10...11] - Boundary of monotonicity interval of U2(U1) function
+  //            (see CylCylMonotonicity() function)
+
+  theU1crit[0] = 0.0;
+  theU1crit[1] = thePeriod;
+  theU1crit[2] = theUSurf1f;
+  theU1crit[3] = theUSurf1l;
+
+  const Standard_Real aCOS = cos(theCoeffs.mFI2);
+  const Standard_Real aBSB = Abs(theCoeffs.mB);
+  if((theCoeffs.mC - aBSB <= aCOS) && (aCOS <= theCoeffs.mC + aBSB))
+  {
+    Standard_Real anArg = (aCOS - theCoeffs.mC) / theCoeffs.mB;
+    if(anArg > 1.0)
+      anArg = 1.0;
+    if(anArg < -1.0)
+      anArg = -1.0;
+
+    theU1crit[4] = -acos(anArg) + theCoeffs.mFI1;
+    theU1crit[5] =  acos(anArg) + theCoeffs.mFI1;
+  }
+
+  Standard_Real aSf = cos(theUSurf2f - theCoeffs.mFI2);
+  Standard_Real aSl = cos(theUSurf2l - theCoeffs.mFI2);
+  MinMax(aSf, aSl);
+
+  //In accorance with pure mathematic, theU1crit[6] and [8]
+  //must be -Precision::Infinite() instead of used +Precision::Infinite()
+  theU1crit[6] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+                  -acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+                  Precision::Infinite();
+  theU1crit[7] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+                  -acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+                  Precision::Infinite();
+  theU1crit[8] = Abs((aSf - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+                  acos((aSf - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+                  Precision::Infinite();
+  theU1crit[9] = Abs((aSl - theCoeffs.mC) / theCoeffs.mB) < 1.0 ?
+                  acos((aSl - theCoeffs.mC) / theCoeffs.mB) + theCoeffs.mFI1 :
+                  Precision::Infinite();
+
+  theU1crit[10] = theCoeffs.mFI1;
+  theU1crit[11] = M_PI+theCoeffs.mFI1;
+
+  //preparative treatment of array. This array must have faled to contain negative
+  //infinity number
+
+  for(Standard_Integer i = 0; i < theNbCritPointsMax; i++)
+  {
+    if(Precision::IsInfinite(theU1crit[i]))
+    {
+      continue;
+    }
+
+    theU1crit[i] = fmod(theU1crit[i], thePeriod);
+    if(theU1crit[i] < 0.0)
+      theU1crit[i] += thePeriod;
+  }
+
+  //Here all not infinite elements of theU1crit are in [0, thePeriod) range
+
+  do
+  {
+    std::sort(theU1crit, theU1crit + theNbCritPointsMax);
+  }
+  while(ExcludeNearElements(theU1crit, theNbCritPointsMax,
+                            theUSurf1f, theUSurf1l, theTol2D));
+
+  //Here all not infinite elements in theU1crit are different and sorted.
 
-      aU1f[0] = anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod + anEquationCoeffs.mFI1;
-    }
-    else if((anEquationCoeffs.mC - 1 <= anEquationCoeffs.mB) &&
-            (anEquationCoeffs.mB <= -anEquationCoeffs.mB - 1))
+  while(theNbCritPointsMax > 0)
+  {
+    Standard_Real &anB = theU1crit[theNbCritPointsMax-1];
+    if(Precision::IsInfinite(anB))
     {
-      Standard_Real anArg = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB;
-      if(anArg > 1.0)
-        anArg = 1.0;
-      if(anArg < -1.0)
-        anArg = -1.0;
+      theNbCritPointsMax--;
+      continue;
+    }
 
-      const Standard_Real aDAngle = acos(anArg);
-      //U=[aDAngle;2*PI-aDAngle]+aFI1
+    //1st not infinte element is found
 
-      aU1f[0] = aDAngle + anEquationCoeffs.mFI1;
-      aU1l[0] = aPeriod - aDAngle + anEquationCoeffs.mFI1;
-    }
-    else if(-anEquationCoeffs.mB - Abs(anEquationCoeffs.mC) >= 1.0)
-    {
-      Standard_Real anArg1 = -(anEquationCoeffs.mC + 1) / anEquationCoeffs.mB,
-                    anArg2 = (1 - anEquationCoeffs.mC) / anEquationCoeffs.mB;
-      if(anArg1 > 1.0)
-        anArg1 = 1.0;
-      if(anArg1 < -1.0)
-        anArg1 = -1.0;
+    if(theNbCritPointsMax == 1)
+      break;
 
-      if(anArg2 > 1.0)
-        anArg2 = 1.0;
-      if(anArg2 < -1.0)
-        anArg2 = -1.0;
+    //Here theNbCritPointsMax > 1
 
-      const Standard_Real aDAngle1 = acos(anArg1), aDAngle2 = acos(anArg2);
-      //(U=[aDAngle1;aDAngle2]+aFI1) ||
-      //(U=[2*PI-aDAngle2;2*PI-aDAngle1]+aFI1)
+    Standard_Real &anA = theU1crit[0];
 
-      aU1f[0] = aDAngle1 + anEquationCoeffs.mFI1;
-      aU1l[0] = aDAngle2 + anEquationCoeffs.mFI1;
-      aU1f[1] = aPeriod - aDAngle2 + anEquationCoeffs.mFI1;
-      aU1l[1] = aPeriod - aDAngle1 + anEquationCoeffs.mFI1;
-    }
-    else
-    {
-      Standard_Failure::Raise("Error. Exception. Unhandled case (Range computation)!");
+    //Compare 1st and last significant elements of theU1crit
+    //They may still differs by period.
+
+    if (Abs(anB - anA - thePeriod) < theTol2D)
+    {//E.g. anA == 2.0e-17, anB == (thePeriod-1.0e-18)
+      anA = (anA + anB - thePeriod)/2.0;
+      anB = Precision::Infinite();
+      theNbCritPointsMax--;
     }
+
+    //Out of "while(theNbCritPointsMax > 0)" cycle.
+    break;
   }
-  else
+
+  //Attention! Here theU1crit may be unsorted.
+}
+
+//=======================================================================
+//function : BoundaryEstimation
+//purpose  : Rough estimation of the parameter range.
+//=======================================================================
+void WorkWithBoundaries::BoundaryEstimation(const gp_Cylinder& theCy1,
+                                            const gp_Cylinder& theCy2,
+                                            Bnd_Range& theOutBoxS1,
+                                            Bnd_Range& theOutBoxS2) const
+{
+  const gp_Dir &aD1 = theCy1.Axis().Direction(),
+               &aD2 = theCy2.Axis().Direction();
+  const Standard_Real aR1 = theCy1.Radius(),
+                      aR2 = theCy2.Radius();
+
+  //Let consider a parallelogram. Its edges are parallel to aD1 and aD2.
+  //Its altitudes are equal to 2*aR1 and 2*aR2 (diameters of the cylinders).
+  //In fact, this parallelogram is a projection of the cylinders to the plane
+  //created by the intersected axes aD1 and aD2 (if the axes are skewed then
+  //one axis can be translated by parallel shifting till intersection).
+
+  const Standard_Real aCosA = aD1.Dot(aD2);
+  const Standard_Real aSqSinA = aD1.XYZ().CrossSquareMagnitude(aD2.XYZ());
+
+  //If sine is small then it can be compared with angle.
+  if (aSqSinA < Precision::Angular()*Precision::Angular())
+    return;
+
+  //Half of delta V. Delta V is a distance between 
+  //projections of two opposite parallelogram vertices
+  //(joined by the maximal diagonal) to the cylinder axis.
+  const Standard_Real aSinA = sqrt(aSqSinA);
+  const Standard_Real anAbsCosA = Abs(aCosA);
+  const Standard_Real aHDV1 = (aR1 * anAbsCosA + aR2) / aSinA,
+                      aHDV2 = (aR2 * anAbsCosA + aR1) / aSinA;
+
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+  //The code in this block is created for test only.It is stupidly to create
+  //OCCT-test for the method, which will be changed possibly never.
+  std::cout << "Reference values: aHDV1 = " << aHDV1 << "; aHDV2 = " << aHDV2 << std::endl;
+#endif
+
+  //V-parameters of intersection point of the axes (in case of skewed axes,
+  //see comment above).
+  Standard_Real aV01 = 0.0, aV02 = 0.0;
+  ExtremaLineLine(theCy1.Axis(), theCy2.Axis(), aCosA, aSqSinA, aV01, aV02);
+
+  theOutBoxS1.Add(aV01 - aHDV1);
+  theOutBoxS1.Add(aV01 + aHDV1);
+
+  theOutBoxS2.Add(aV02 - aHDV2);
+  theOutBoxS2.Add(aV02 + aHDV2);
+
+  theOutBoxS1.Enlarge(Precision::Confusion());
+  theOutBoxS2.Enlarge(Precision::Confusion());
+
+  Standard_Real aU1 = 0.0, aV1 = 0.0, aU2 = 0.0, aV2 = 0.0;
+  myUVSurf1.Get(aU1, aV1, aU2, aV2);
+  theOutBoxS1.Common(Bnd_Range(aV1, aV2));
+
+  myUVSurf2.Get(aU1, aV1, aU2, aV2);
+  theOutBoxS2.Common(Bnd_Range(aV1, aV2));
+}
+
+//=======================================================================
+//function : CyCyNoGeometric
+//purpose  : 
+//=======================================================================
+static IntPatch_ImpImpIntersection::IntStatus
+                    CyCyNoGeometric(const gp_Cylinder &theCyl1,
+                                    const gp_Cylinder &theCyl2,
+                                    const WorkWithBoundaries &theBW,
+                                    Bnd_Range theRange[],
+                                    const Standard_Integer theNbOfRanges /*=2*/,
+                                    Standard_Boolean& isTheEmpty,
+                                    IntPatch_SequenceOfLine& theSlin,
+                                    IntPatch_SequenceOfPoint& theSPnt)
+{
+  Standard_Real aUSurf1f = 0.0, aUSurf1l = 0.0,
+                aUSurf2f = 0.0, aUSurf2l = 0.0,
+                aVSurf1f = 0.0, aVSurf1l = 0.0,
+                aVSurf2f = 0.0, aVSurf2l = 0.0;
+
+  theBW.UVS1().Get(aUSurf1f, aVSurf1f, aUSurf1l, aVSurf1l);
+  theBW.UVS2().Get(aUSurf2f, aVSurf2f, aUSurf2l, aVSurf2l);
+
+  Bnd_Range aRangeS1, aRangeS2;
+  theBW.BoundaryEstimation(theCyl1, theCyl2, aRangeS1, aRangeS2);
+  if (aRangeS1.IsVoid() || aRangeS2.IsVoid())
+    return IntPatch_ImpImpIntersection::IntStatus_OK;
+
   {
-    Standard_Failure::Raise("Error. Exception. Unhandled case (B-parameter computation)!");
+    //Quotation of the message from issue #26894 (author MSV):
+    //"We should return fail status from intersector if the result should be an
+    //infinite curve of non-analytical type... I propose to define the limit for the
+    //extent as the radius divided by 1e+2 and multiplied by 1e+7.
+    //Thus, taking into account the number of valuable digits (15), we provide reliable
+    //computations with an error not exceeding R/100."
+    const Standard_Real aF = 1.0e+5;
+    const Standard_Real aMaxV1Range = aF*theCyl1.Radius(), aMaxV2Range = aF*theCyl2.Radius();
+    if ((aRangeS1.Delta() > aMaxV1Range) || (aRangeS2.Delta() > aMaxV2Range))
+      return IntPatch_ImpImpIntersection::IntStatus_InfiniteSectionCurve;
   }
 
-  for(Standard_Integer i = 0; i < aNbOfBoundaries; i++)
+  const ComputationMethods::stCoeffsValue &anEquationCoeffs = theBW.SICoeffs();
+  const IntSurf_Quadric& aQuad1 = theBW.GetQSurface(1);
+  const IntSurf_Quadric& aQuad2 = theBW.GetQSurface(2);
+  const Standard_Boolean isReversed = theBW.IsReversed();
+  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;
+
+  //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.
+
+  for (Standard_Integer i = 0; i < theNbOfRanges; i++)
   {
-    if(Precision::IsInfinite(aU1f[i]) && Precision::IsInfinite(aU1l[i]))
+    if (theRange[i].IsVoid())
       continue;
 
-    InscribeInterval(aUSurf1f, aUSurf1l, aU1f[i], aU1l[i], theTol2D, aPeriod);
+    InscribeInterval(aUSurf1f, aUSurf1l, theRange[i], aTol2D, aPeriod);
   }
 
-  if( !Precision::IsInfinite(aU1f[0]) && !Precision::IsInfinite(aU1f[1]) &&
-      !Precision::IsInfinite(aU1l[0]) && !Precision::IsInfinite(aU1l[1]))
+  if (theRange[0].Union(theRange[1]))
   {
-    if( ((aU1f[1] <= aU1l[0]) || (aU1l[1] <= aU1l[0])) && 
-        ((aU1f[0] <= aU1l[1]) || (aU1l[0] <= aU1l[1])))
-    {//Join all intervals to one
-      aU1f[0] = Min(aU1f[0], aU1f[1]);
-      aU1l[0] = Max(aU1l[0], aU1l[1]);
-
-      aU1f[1] = -Precision::Infinite();
-      aU1l[1] = Precision::Infinite();
-    }
+    // Works only if (theNbOfRanges == 2).
+    theRange[1].SetVoid();
   }
 
-  //Critical points
-  //[0...1] - in these points parameter U1 goes through
-  //          the seam-edge of the first cylinder.
-  //[2...3] - First and last U1 parameter.
-  //[4...5] - in these points parameter U2 goes through
-  //          the seam-edge of the second cylinder.
-  //[6...9] - in these points an intersection line goes through
-  //          U-boundaries of the second surface.
-  const Standard_Integer aNbCritPointsMax = 10;
-  Standard_Real anU1crit[aNbCritPointsMax] = {Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite(),
-                                              Precision::Infinite()};
-
+  //Critical points are the value of U1-parameter in the points
+  //where WL must be decomposed.
+
+  //When U1 goes through critical points its value is set up to this
+  //parameter forcefully and the intersection point is added in the line.
+  //After that, the WL is broken (next U1 value will be correspond to the new WL).
+
+  //See CriticalPointsComputing(...) function to get detail information about this array.
+  const Standard_Integer aNbCritPointsMax = 12;
+  Standard_Real anU1crit[aNbCritPointsMax] = { Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite(),
+                                               Precision::Infinite() };
+
+  //This list of critical points is not full because it does not contain any points
+  //which intersection line goes through V-bounds of cylinders in.
+  //They are computed by numerical methods on - line (during algorithm working).
+  //The moment is caught, when intersection line goes through V-bounds of any cylinder.
+
+  Standard_Integer aNbCritPoints = aNbCritPointsMax;
   CriticalPointsComputing(anEquationCoeffs, aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
-                                      aPeriod, theTol2D, aNbCritPointsMax, anU1crit);
-
+                          aPeriod, aTol2D, aNbCritPoints, anU1crit);
 
   //Getting Walking-line
 
   enum WLFStatus
   {
+    // No points have been added in WL
     WLFStatus_Absent = 0,
-    WLFStatus_Exist  = 1,
+    // WL contains at least one point
+    WLFStatus_Exist = 1,
+    // WL has been finished in some critical point
+    // We should start new line
     WLFStatus_Broken = 2
   };
 
-  for(Standard_Integer aCurInterval = 0; aCurInterval < aNbOfBoundaries; aCurInterval++)
+  const Standard_Integer aNbWLines = 2;
+  for (Standard_Integer aCurInterval = 0; aCurInterval < theNbOfRanges; aCurInterval++)
   {
-    if(Precision::IsInfinite(aU1f[aCurInterval]) && Precision::IsInfinite(aU1l[aCurInterval]))
-      continue;
-
+    //Process every continuous region
     Standard_Boolean isAddedIntoWL[aNbWLines];
-    for(Standard_Integer i = 0; i < aNbWLines; i++)
+    for (Standard_Integer i = 0; i < aNbWLines; i++)
       isAddedIntoWL[i] = Standard_False;
 
-    Standard_Real anUf = aU1f[aCurInterval], anUl = aU1l[aCurInterval];
-    const Standard_Boolean isDeltaPeriod  = IsEqual(anUl-anUf, aPeriod);
+    Standard_Real anUf = 1.0, anUl = 0.0;
+    if (!theRange[aCurInterval].GetBounds(anUf, anUl))
+      continue;
+
+    const Standard_Boolean isDeltaPeriod = IsEqual(anUl - anUf, aPeriod);
 
     //Inscribe and sort critical points
-    InscribeAndSortArray(anU1crit, aNbCritPointsMax, anUf, anUl, theTol2D, aPeriod);
+    for (Standard_Integer i = 0; i < aNbCritPoints; i++)
+    {
+      InscribePoint(anUf, anUl, anU1crit[i], 0.0, aPeriod, Standard_False);
+    }
+
+    std::sort(anU1crit, anU1crit + aNbCritPoints);
 
-    while(anUf < anUl)
+    while (anUf < anUl)
     {
+      //Change value of U-parameter on the 1st surface from anUf to anUl
+      //(anUf will be modified in the cycle body).
+      //Step is computed adaptively (see comments below).
+
       Standard_Real aU2[aNbWLines], aV1[aNbWLines], aV2[aNbWLines];
       WLFStatus aWLFindStatus[aNbWLines];
       Standard_Real aV1Prev[aNbWLines], aV2Prev[aNbWLines];
-      Standard_Real anArccosFactor[aNbWLines] = {1.0, -1.0};
+      Standard_Real anUexpect[aNbWLines];
       Standard_Boolean isAddingWLEnabled[aNbWLines];
 
       Handle(IntSurf_LineOn2S) aL2S[aNbWLines];
       Handle(IntPatch_WLine) aWLine[aNbWLines];
-      for(Standard_Integer i = 0; i < aNbWLines; i++)
+      for (Standard_Integer i = 0; i < aNbWLines; i++)
       {
         aL2S[i] = new IntSurf_LineOn2S();
         aWLine[i] = new IntPatch_WLine(aL2S[i], Standard_False);
+        aWLine[i]->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
         aWLFindStatus[i] = WLFStatus_Absent;
         isAddingWLEnabled[i] = Standard_True;
         aU2[i] = aV1[i] = aV2[i] = 0.0;
         aV1Prev[i] = aV2Prev[i] = 0.0;
+        anUexpect[i] = anUf;
       }
-      
-      Standard_Real anU1 = anUf;
 
-      Standard_Real aCriticalDelta[aNbCritPointsMax];
-      for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
-        aCriticalDelta[i] = anU1 - anU1crit[i];
+      Standard_Real aCriticalDelta[aNbCritPointsMax] = { 0 };
+      for (Standard_Integer aCritPID = 0; aCritPID < aNbCritPoints; aCritPID++)
+      { //We are not interested in elements of aCriticalDelta array
+        //if their index is greater than or equal to aNbCritPoints
 
+        aCriticalDelta[aCritPID] = anUf - anU1crit[aCritPID];
+      }
+
+      Standard_Real anU1 = anUf, aMinCriticalParam = anUf;
       Standard_Boolean isFirst = Standard_True;
 
-      while(anU1 <= anUl)
+      while (anU1 <= anUl)
       {
-        if(isDeltaPeriod)
-        {
-          if(IsEqual(anU1, anUl))
-          {
-            //if isAddedIntoWL* == TRUE WLine contains only one point
-            //(which was end point of previous WLine). If we will
-            //add point found on the current step WLine will contain only
-            //two points. At that both these points will be equal to the
-            //points found earlier. Therefore, new WLine will repeat 
-            //already existing WLine. Consequently, it is necessary 
-            //to forbid building new line in this case.
-
-            for(Standard_Integer i = 0; i < aNbWLines; i++)
-              isAddingWLEnabled[i] = !isAddedIntoWL[i];
-          }
-        }
+        //Change value of U-parameter on the 1st surface from anUf to anUl
+        //(anUf will be modified in the cycle body). However, this cycle
+        //can be broken if WL goes though some critical point.
+        //Step is computed adaptively (see comments below).
 
-        for(Standard_Integer i = 0; i < aNbCritPointsMax; i++)
+        for (Standard_Integer i = 0; i < aNbCritPoints; i++)
         {
-          if((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
+          if ((anU1 - anU1crit[i])*aCriticalDelta[i] < 0.0)
           {
+            //WL has gone through i-th critical point
             anU1 = anU1crit[i];
 
-            for(Standard_Integer i = 0; i < aNbWLines; i++)
-              aWLFindStatus[i] = WLFStatus_Broken;
+            for (Standard_Integer j = 0; j < aNbWLines; j++)
+            {
+              aWLFindStatus[j] = WLFStatus_Broken;
+              anUexpect[j] = anU1;
+            }
 
             break;
           }
         }
 
-        Standard_Real anArg = anEquationCoeffs.mB * 
-              cos(anU1 - anEquationCoeffs.mFI1) + anEquationCoeffs.mC;
-        
-        if(aNulValue > 1.0 - anArg)
-          anArg = 1.0;
-        if(anArg + 1.0 < aNulValue)
-          anArg = -1.0;
-
-        for(Standard_Integer i = 0; i < aNbWLines; i++)
+        if (IsEqual(anU1, anUl))
         {
-          const Standard_Integer aNbPntsWL =  aWLine[i].IsNull() ? 0 :
-                                              aWLine[i]->Curve()->NbPoints();
-          aU2[i] = anEquationCoeffs.mFI2 + anArccosFactor[i]*acos(anArg);
-
-          InscribePoint(aUSurf2f, aUSurf2l, aU2[i], theTol2D, aPeriod, Standard_False);
+          for (Standard_Integer i = 0; i < aNbWLines; i++)
+          {
+            aWLFindStatus[i] = WLFStatus_Broken;
+            anUexpect[i] = anU1;
 
-          if(aNbPntsWL == 0)
-          {//the line has not contained any points yet
-            if(((aUSurf2l - aUSurf2f) >= aPeriod) && 
-                ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
-                  (Abs(aU2[i]-aUSurf2l) < theTol2D)))
+            if (isDeltaPeriod)
+            {
+              //if isAddedIntoWL[i] == TRUE WLine contains only one point
+              //(which was end point of previous WLine). If we will
+              //add point found on the current step WLine will contain only
+              //two points. At that both these points will be equal to the
+              //points found earlier. Therefore, new WLine will repeat 
+              //already existing WLine. Consequently, it is necessary 
+              //to forbid building new line in this case.
+
+              isAddingWLEnabled[i] = (!isAddedIntoWL[i]);
+            }
+            else
             {
-              const Standard_Real anU1Temp = anU1 + aStepMin;
-              Standard_Real anArgTemp = anEquationCoeffs.mB * 
-                                        cos(anU1Temp - anEquationCoeffs.mFI1) +
-                                        anEquationCoeffs.mC;
+              isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
+                                      (aWLFindStatus[i] == WLFStatus_Absent));
+            }
+          }//for(Standard_Integer i = 0; i < aNbWLines; i++)
+        }
+        else
+        {
+          for (Standard_Integer i = 0; i < aNbWLines; i++)
+          {
+            isAddingWLEnabled[i] = ((aTol2D >= (anUexpect[i] - anU1)) ||
+                                    (aWLFindStatus[i] == WLFStatus_Absent));
+          }//for(Standard_Integer i = 0; i < aNbWLines; i++)
+        }
+
+        for (Standard_Integer i = 0; i < aNbWLines; i++)
+        {
+          const Standard_Integer aNbPntsWL = aWLine[i].IsNull() ? 0 :
+            aWLine[i]->Curve()->NbPoints();
 
-              if(aNulValue > 1.0 - anArgTemp)
-                anArgTemp = 1.0;
+          if ((aWLFindStatus[i] == WLFStatus_Broken) ||
+            (aWLFindStatus[i] == WLFStatus_Absent))
+          {//Begin and end of WLine must be on boundary point
+           //or on seam-edge strictly (if it is possible).
 
-              if(anArgTemp + 1.0 < aNulValue)
-                anArgTemp = -1.0;
+            Standard_Real aTol = aTol2D;
+            ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs,
+                                                        aU2[i], &aTol);
+            InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
 
-              Standard_Real aU2Temp = anEquationCoeffs.mFI2 +
-                                      anArccosFactor[i]*acos(anArgTemp);
+            aTol = Max(aTol, aTol2D);
 
-              InscribePoint(aUSurf2f, aUSurf2l, aU2Temp, theTol2D, aPeriod, Standard_False);
+            if (Abs(aU2[i]) <= aTol)
+              aU2[i] = 0.0;
+            else if (Abs(aU2[i] - aPeriod) <= aTol)
+              aU2[i] = aPeriod;
+            else if (Abs(aU2[i] - aUSurf2f) <= aTol)
+              aU2[i] = aUSurf2f;
+            else if (Abs(aU2[i] - aUSurf2l) <= aTol)
+              aU2[i] = aUSurf2l;
+          }
+          else
+          {
+            ComputationMethods::CylCylComputeParameters(anU1, i, anEquationCoeffs, aU2[i]);
+            InscribePoint(aUSurf2f, aUSurf2l, aU2[i], aTol2D, aPeriod, Standard_False);
+          }
 
-              if(2.0*Abs(aU2Temp - aU2[i]) > aPeriod)
+          if (aNbPntsWL == 0)
+          {//the line has not contained any points yet
+            if (((aUSurf2f + aPeriod - aUSurf2l) <= 2.0*aTol2D) &&
+                ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
+                  (Abs(aU2[i] - aUSurf2l) < aTol2D)))
+            {
+              //In this case aU2[i] can have two values: current aU2[i] or
+              //aU2[i]+aPeriod (aU2[i]-aPeriod). It is necessary to choose
+              //correct value.
+
+              Standard_Boolean isIncreasing = Standard_True;
+              ComputationMethods::CylCylMonotonicity(anU1+aStepMin, i, anEquationCoeffs,
+                                                      aPeriod, isIncreasing);
+
+              //If U2(U1) is increasing and U2 is considered to be equal aUSurf2l
+              //then after the next step (when U1 will be increased) U2 will be
+              //increased too. And we will go out of surface boundary.
+              //Therefore, If U2(U1) is increasing then U2 must be equal aUSurf2f.
+              //Analogically, if U2(U1) is decreasing.
+              if (isIncreasing)
               {
-                if(aU2Temp > aU2[i])
-                  aU2[i] += aPeriod;
-                else
-                  aU2[i] -= aPeriod;
+                aU2[i] = aUSurf2f;
+              }
+              else
+              {
+                aU2[i] = aUSurf2l;
               }
             }
           }
-
-          if(aNbPntsWL > 0)
-          {
-            if(((aUSurf2l - aUSurf2f) >= aPeriod) && 
-                ((Abs(aU2[i] - aUSurf2f) < theTol2D) ||
-                  (Abs(aU2[i]-aUSurf2l) < theTol2D)))
+          else
+          {//aNbPntsWL > 0
+            if (((aUSurf2l - aUSurf2f) >= aPeriod) &&
+                ((Abs(aU2[i] - aUSurf2f) < aTol2D) ||
+                  (Abs(aU2[i] - aUSurf2l) < aTol2D)))
             {//end of the line
               Standard_Real aU2prev = 0.0, aV2prev = 0.0;
-              if(isTheReverse)
+              if (isReversed)
                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS1(aU2prev, aV2prev);
               else
                 aWLine[i]->Curve()->Value(aNbPntsWL).ParametersOnS2(aU2prev, aV2prev);
 
-              if(2.0*Abs(aU2prev - aU2[i]) > aPeriod)
+              if (2.0*Abs(aU2prev - aU2[i]) > aPeriod)
               {
-                if(aU2prev > aU2[i])
+                if (aU2prev > aU2[i])
                   aU2[i] += aPeriod;
                 else
                   aU2[i] -= aPeriod;
@@ -2012,29 +2971,10 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
             }
           }
 
-          if(aWLFindStatus[i] == WLFStatus_Broken)
-          {
-            if(Abs(aU2[i]) <= theTol2D)
-              aU2[i] = 0.0;
-            else if(Abs(aU2[i] - aPeriod) <= theTol2D)
-              aU2[i] = aPeriod;
-            else if(Abs(aU2[i] - aUSurf2f) <= theTol2D)
-              aU2[i] = aUSurf2f;
-            else if(Abs(aU2[i] - aUSurf2l) <= theTol2D)
-              aU2[i] = aUSurf2l;
-          }
-
-          aV1[i] =  anEquationCoeffs.mK21 * sin(aU2[i]) + 
-                    anEquationCoeffs.mK11 * sin(anU1) +
-                    anEquationCoeffs.mL21 * cos(aU2[i]) +
-                    anEquationCoeffs.mL11 * cos(anU1) + anEquationCoeffs.mM1;
-
-          aV2[i] =  anEquationCoeffs.mK22 * sin(aU2[i]) +
-                    anEquationCoeffs.mK12 * sin(anU1) +
-                    anEquationCoeffs.mL22 * cos(aU2[i]) +
-                    anEquationCoeffs.mL12 * cos(anU1) + anEquationCoeffs.mM2;
+          ComputationMethods::CylCylComputeParameters(anU1, aU2[i], anEquationCoeffs,
+                                                              aV1[i], aV2[i]);
 
-          if(isFirst)
+          if (isFirst)
           {
             aV1Prev[i] = aV1[i];
             aV2Prev[i] = aV2[i];
@@ -2045,163 +2985,374 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
 
         //Looking for points into WLine
         Standard_Boolean isBroken = Standard_False;
-        for(Standard_Integer i = 0; i < aNbWLines; i++)
+        for (Standard_Integer i = 0; i < aNbWLines; i++)
         {
-          if(!isAddingWLEnabled[i])
+          if (!isAddingWLEnabled[i])
           {
-            aV1Prev[i] = aV1[i];
-            aV2Prev[i] = aV2[i];
+            Standard_Boolean isBoundIntersect = Standard_False;
+            if ((Abs(aV1[i] - aVSurf1f) <= aTol2D) ||
+                ((aV1[i] - aVSurf1f)*(aV1Prev[i] - aVSurf1f) < 0.0))
+            {
+              isBoundIntersect = Standard_True;
+            }
+            else if ((Abs(aV1[i] - aVSurf1l) <= aTol2D) ||
+                    ((aV1[i] - aVSurf1l)*(aV1Prev[i] - aVSurf1l) < 0.0))
+            {
+              isBoundIntersect = Standard_True;
+            }
+            else if ((Abs(aV2[i] - aVSurf2f) <= aTol2D) ||
+                    ((aV2[i] - aVSurf2f)*(aV2Prev[i] - aVSurf2f) < 0.0))
+            {
+              isBoundIntersect = Standard_True;
+            }
+            else if ((Abs(aV2[i] - aVSurf2l) <= aTol2D) ||
+                    ((aV2[i] - aVSurf2l)*(aV2Prev[i] - aVSurf2l) < 0.0))
+            {
+              isBoundIntersect = Standard_True;
+            }
+
+            if (aWLFindStatus[i] == WLFStatus_Broken)
+              isBroken = Standard_True;
+
+            if (!isBoundIntersect)
+            {
+              continue;
+            }
+            else
+            {
+              anUexpect[i] = anU1;
+            }
+          }
+
+          // True if the current point already in the domain
+          const Standard_Boolean isInscribe =
+              ((aUSurf2f - aU2[i]) <= aTol2D) && ((aU2[i] - aUSurf2l) <= aTol2D) &&
+              ((aVSurf1f - aV1[i]) <= aTol2D) && ((aV1[i] - aVSurf1l) <= aTol2D) &&
+              ((aVSurf2f - aV2[i]) <= aTol2D) && ((aV2[i] - aVSurf2l) <= aTol2D);
+
+          //isVIntersect == TRUE if intersection line intersects two (!)
+          //V-bounds of cylinder (1st or 2nd - no matter)
+          const Standard_Boolean isVIntersect =
+              (((aVSurf1f - aV1[i])*(aVSurf1f - aV1Prev[i]) < RealSmall()) &&
+                ((aVSurf1l - aV1[i])*(aVSurf1l - aV1Prev[i]) < RealSmall())) ||
+              (((aVSurf2f - aV2[i])*(aVSurf2f - aV2Prev[i]) < RealSmall()) &&
+                ((aVSurf2l - aV2[i])*(aVSurf2l - aV2Prev[i]) < RealSmall()));
+
+          //isFound1 == TRUE if intersection line intersects V-bounds
+          //  (First or Last - no matter) of the 1st cylynder
+          //isFound2 == TRUE if intersection line intersects V-bounds
+          //  (First or Last - no matter) of the 2nd cylynder
+          Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
+          Standard_Boolean isForce = Standard_False;
+
+          if (aWLFindStatus[i] == WLFStatus_Absent)
+          {
+            if (((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1 - aUSurf1l) < aTol2D))
+            {
+              isForce = Standard_True;
+            }
+          }
+
+          theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
+                                 aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i, isForce,
+                                 isFound1, isFound2);
+
+          const Standard_Boolean isPrevVBound = !isVIntersect &&
+                                                ((Abs(aV1Prev[i] - aVSurf1f) <= aTol2D) ||
+                                                 (Abs(aV1Prev[i] - aVSurf1l) <= aTol2D) ||
+                                                 (Abs(aV2Prev[i] - aVSurf2f) <= aTol2D) ||
+                                                 (Abs(aV2Prev[i] - aVSurf2l) <= aTol2D));
+
+          aV1Prev[i] = aV1[i];
+          aV2Prev[i] = aV2[i];
+
+          if ((aWLFindStatus[i] == WLFStatus_Exist) && (isFound1 || isFound2) && !isPrevVBound)
+          {
+            aWLFindStatus[i] = WLFStatus_Broken; //start a new line
+          }
+          else if (isInscribe)
+          {
+            if ((aWLFindStatus[i] == WLFStatus_Absent) && (isFound1 || isFound2))
+            {
+              aWLFindStatus[i] = WLFStatus_Exist;
+            }
+
+            if ((aWLFindStatus[i] != WLFStatus_Broken) ||
+              (aWLine[i]->NbPnts() >= 1) || IsEqual(anU1, anUl))
+            {
+              if (aWLine[i]->NbPnts() > 0)
+              {
+                Standard_Real aU2p = 0.0, aV2p = 0.0;
+                if (isReversed)
+                  aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
+                else
+                  aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
+
+                const Standard_Real aDelta = aU2[i] - aU2p;
+
+                if (2.0 * Abs(aDelta) > aPeriod)
+                {
+                  if (aDelta > 0.0)
+                  {
+                    aU2[i] -= aPeriod;
+                  }
+                  else
+                  {
+                    aU2[i] += aPeriod;
+                  }
+                }
+              }
+
+              if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
+                                gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
+                                aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+                                aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod,
+                                aWLine[i]->Curve(), i, aTol3D, aTol2D, isForce))
+              {
+                if (aWLFindStatus[i] == WLFStatus_Absent)
+                {
+                  aWLFindStatus[i] = WLFStatus_Exist;
+                }
+              }
+              else if (!isFound1 && !isFound2)
+              {//We do not add any point while doing this iteration
+                if (aWLFindStatus[i] == WLFStatus_Exist)
+                {
+                  aWLFindStatus[i] = WLFStatus_Broken;
+                }
+              }
+            }
+          }
+          else
+          {//We do not add any point while doing this iteration
+            if (aWLFindStatus[i] == WLFStatus_Exist)
+            {
+              aWLFindStatus[i] = WLFStatus_Broken;
+            }
+          }
+
+          if (aWLFindStatus[i] == WLFStatus_Broken)
+            isBroken = Standard_True;
+        }//for(Standard_Integer i = 0; i < aNbWLines; i++)
+
+        if (isBroken)
+        {//current lines are filled. Go to the next lines
+          anUf = anU1;
+
+          Standard_Boolean isAdded = Standard_True;
+
+          for (Standard_Integer i = 0; i < aNbWLines; i++)
+          {
+            if (isAddingWLEnabled[i])
+            {
+              continue;
+            }
+
+            isAdded = Standard_False;
+
+            Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
+
+            theBW.AddBoundaryPoint(aWLine[i], anU1, aMinCriticalParam, aU2[i],
+                                   aV1[i], aV1Prev[i], aV2[i], aV2Prev[i], i,
+                                   Standard_False, isFound1, isFound2);
+
+            if (isFound1 || isFound2)
+            {
+              isAdded = Standard_True;
+            }
+
+            if (aWLine[i]->NbPnts() > 0)
+            {
+              Standard_Real aU2p = 0.0, aV2p = 0.0;
+              if (isReversed)
+                aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS1(aU2p, aV2p);
+              else
+                aWLine[i]->Point(aWLine[i]->NbPnts()).ParametersOnS2(aU2p, aV2p);
+
+              const Standard_Real aDelta = aU2[i] - aU2p;
+
+              if (2 * Abs(aDelta) > aPeriod)
+              {
+                if (aDelta > 0.0)
+                {
+                  aU2[i] -= aPeriod;
+                }
+                else
+                {
+                  aU2[i] += aPeriod;
+                }
+              }
+            }
+
+            if(AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
+                              Standard_True, gp_Pnt2d(anU1, aV1[i]),
+                              gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
+                              aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
+                              aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
+                              i, aTol3D, aTol2D, Standard_False))
+            {
+              isAdded = Standard_True;
+            }
+          }
+
+          if (!isAdded)
+          {
+            //Before breaking WL, we must complete it correctly
+            //(e.g. to prolong to the surface boundary).
+            //Therefore, we take the point last added in some WL
+            //(have maximal U1-parameter) and try to add it in
+            //the current WL.
+            Standard_Real anUmaxAdded = RealFirst();
+
+            {
+              Standard_Boolean isChanged = Standard_False;
+              for (Standard_Integer i = 0; i < aNbWLines; i++)
+              {
+                if ((aWLFindStatus[i] == WLFStatus_Absent) || (aWLine[i]->NbPnts() == 0))
+                  continue;
+
+                Standard_Real aU1c = 0.0, aV1c = 0.0;
+                if (isReversed)
+                  aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS2(aU1c, aV1c);
+                else
+                  aWLine[i]->Curve()->Value(aWLine[i]->NbPnts()).ParametersOnS1(aU1c, aV1c);
+
+                anUmaxAdded = Max(anUmaxAdded, aU1c);
+                isChanged = Standard_True;
+              }
+
+              if (!isChanged)
+              { //If anUmaxAdded were not changed in previous cycle then
+                //we would break existing WLines.
+                break;
+              }
+            }
+
+            for (Standard_Integer i = 0; i < aNbWLines; i++)
+            {
+              if (isAddingWLEnabled[i])
+              {
+                continue;
+              }
+
+              ComputationMethods::CylCylComputeParameters(anUmaxAdded, i, anEquationCoeffs,
+                aU2[i], aV1[i], aV2[i]);
+
+              AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
+                             Standard_True, gp_Pnt2d(anUmaxAdded, aV1[i]),
+                             gp_Pnt2d(aU2[i], aV2[i]), aUSurf1f, aUSurf1l,
+                             aUSurf2f, aUSurf2l, aVSurf1f, aVSurf1l,
+                             aVSurf2f, aVSurf2l, aPeriod, aWLine[i]->Curve(),
+                             i, aTol3D, aTol2D, Standard_False);
+            }
+          }
+
+          break;
+        }
+
+        //Step computing
+
+        {
+          //Step of aU1-parameter is computed adaptively. The algorithm 
+          //aims to provide given aDeltaV1 and aDeltaV2 values (if it is 
+          //possible because the intersection line can go along V-isoline)
+          //in every iteration. It allows avoiding "flying" intersection
+          //points too far each from other (see issue #24915).
 
-            if(aWLFindStatus[i] == WLFStatus_Broken)
-              isBroken = Standard_True;
+          const Standard_Real aDeltaV1 = aRangeS1.Delta() / IntToReal(aNbPoints),
+                              aDeltaV2 = aRangeS2.Delta() / IntToReal(aNbPoints);
 
-            continue;
-          }
+          math_Matrix aMatr(1, 3, 1, 5);
 
-          if( ((aUSurf2f-aU2[i]) <= theTol2D) && ((aU2[i]-aUSurf2l) <= theTol2D) &&
-              ((aVSurf1f - aV1[i]) <= theTol2D) && ((aV1[i] - aVSurf1l) <= theTol2D) &&
-              ((aVSurf2f - aV2[i]) <= theTol2D) && ((aV2[i] - aVSurf2l) <= theTol2D))
+          Standard_Real aMinUexp = RealLast();
+
+          for (Standard_Integer i = 0; i < aNbWLines; i++)
           {
-            Standard_Boolean isForce = Standard_False;
-            if(aWLFindStatus[i] == WLFStatus_Absent)
+            if (aTol2D < (anUexpect[i] - anU1))
             {
-              Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
-
-              if(((aUSurf2l - aUSurf2f) >= aPeriod) && (Abs(anU1-aUSurf1l) < theTol2D))
-              {
-                isForce = Standard_True;
-              }
-
-              AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
-                                theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
-                                aNulValue, anU1, aU2[i], aV1[i], aV1Prev[i],
-                                aV2[i], aV2Prev[i], isTheReverse,
-                                anArccosFactor[i], isForce, isFound1, isFound2);
-              
-              if(isFound1 || isFound2)
-              {
-                aWLFindStatus[i] = WLFStatus_Exist;
-              }
+              continue;
             }
 
-            if(( aWLFindStatus[i] != WLFStatus_Broken) || (aWLine[i]->NbPnts() >= 1))
+            if (aWLFindStatus[i] == WLFStatus_Absent)
             {
-              if(AddPointIntoWL(theQuad1, theQuad2, isTheReverse, 
-                                gp_Pnt2d(anU1, aV1[i]), gp_Pnt2d(aU2[i], aV2[i]),
-                                aUSurf1f, aUSurf1l, aPeriod,
-                                aWLine[i]->Curve(), theTol3D, theTol2D, isForce))
-              {
-                if(aWLFindStatus[i] == WLFStatus_Absent)
-                {
-                  aWLFindStatus[i] = WLFStatus_Exist;
-                }
-              }
+              anUexpect[i] += aStepMax;
+              aMinUexp = Min(aMinUexp, anUexpect[i]);
+              continue;
             }
-          }
-          else
-          {
-            if(aWLFindStatus[i] == WLFStatus_Exist)
-            {
-              Standard_Boolean isFound1 = Standard_False, isFound2 = Standard_False;
-
-              AddBoundaryPoint( theQuad1, theQuad2, aWLine[i], anEquationCoeffs,
-                                theUVSurf1, theUVSurf2, theTol3D, theTol2D, aPeriod,
-                                aNulValue, anU1, aU2[i], aV1[i], aV1Prev[i],
-                                aV2[i], aV2Prev[i], isTheReverse,
-                                anArccosFactor[i], Standard_False, isFound1, isFound2);
 
-              if(isFound1 || isFound2)
-                aWLFindStatus[i] = WLFStatus_Broken; //start a new line
-            }
-          }
+            Standard_Real aStepTmp = aStepMax;
 
-          aV1Prev[i] = aV1[i];
-          aV2Prev[i] = aV2[i];
+            const Standard_Real aSinU1 = sin(anU1),
+                                aCosU1 = cos(anU1),
+                                aSinU2 = sin(aU2[i]),
+                                aCosU2 = cos(aU2[i]);
 
-          if(aWLFindStatus[i] == WLFStatus_Broken)
-            isBroken = Standard_True;
-        }//for(Standard_Integer i = 0; i < aNbWLines; i++)
+            aMatr.SetCol(1, anEquationCoeffs.mVecC1);
+            aMatr.SetCol(2, anEquationCoeffs.mVecC2);
+            aMatr.SetCol(3, anEquationCoeffs.mVecA1*aSinU1 - anEquationCoeffs.mVecB1*aCosU1);
+            aMatr.SetCol(4, anEquationCoeffs.mVecA2*aSinU2 - anEquationCoeffs.mVecB2*aCosU2);
+            aMatr.SetCol(5, anEquationCoeffs.mVecA1*aCosU1 + anEquationCoeffs.mVecB1*aSinU1 +
+                            anEquationCoeffs.mVecA2*aCosU2 + anEquationCoeffs.mVecB2*aSinU2 +
+                            anEquationCoeffs.mVecD);
 
-        if(isBroken)
-        {//current lines are filled. Go to the next lines
-          anUf = anU1;
-          break;
-        }
+            //The main idea is in solving of linearized system (2)
+            //(see description to ComputationMethods class) in order to find new U1-value
+            //to provide new value V1 or V2, which differs from current one by aDeltaV1 or
+            //aDeltaV2 respectively. 
 
-        //Step computing
+            //While linearizing, following Taylor formulas are used:
+            //    cos(x0+dx) = cos(x0) - sin(x0)*dx
+            //    sin(x0+dx) = sin(x0) + cos(x0)*dx
 
-        Standard_Real aStepU1 = aStepMax;
+            //Consequently cos(U1), cos(U2), sin(U1) and sin(U2) in the system (2)
+            //must be substituted by corresponding values.
 
-        for(Standard_Integer i = 0; i < aNbWLines; i++)
-        {
-          Standard_Real aDeltaU1V1 = aStepU1, aDeltaU1V2 = aStepU1;
-
-          Standard_Real aFact1 = !IsEqual(sin(aU2[i] - anEquationCoeffs.mFI2), 0.0) ? 
-                                  anEquationCoeffs.mK1 * sin(anU1 - anEquationCoeffs.mFIV1) +
-                                  anEquationCoeffs.mL1 * anEquationCoeffs.mB * sin(aU2[i] - anEquationCoeffs.mPSIV1) *
-                                  sin(anU1 - anEquationCoeffs.mFI1)/sin(aU2[i]-anEquationCoeffs.mFI2) : 0.0,
-                        aFact2 = !IsEqual(sin(aU2[i]-anEquationCoeffs.mFI2), 0.0) ? 
-                                  anEquationCoeffs.mK2 * sin(anU1 - anEquationCoeffs.mFIV2) + 
-                                  anEquationCoeffs.mL2 * anEquationCoeffs.mB * sin(aU2[i] - anEquationCoeffs.mPSIV2) *
-                                  sin(anU1 - anEquationCoeffs.mFI1)/sin(aU2[i] - anEquationCoeffs.mFI2) : 0.0;
-
-          Standard_Real aDeltaV1 = (aVSurf1l - aVSurf1f)/IntToReal(aNbPoints),
-                        aDeltaV2 = (aVSurf2l - aVSurf2f)/IntToReal(aNbPoints);
-
-          if((aV1[i] < aVSurf1f) && (aFact1 < 0.0))
-          {//Make close to aVSurf1f by increasing anU1
-            aDeltaV1 = Min(aDeltaV1, Abs(aV1[i]-aVSurf1f));
-          }
+            //ATTENTION!!!
+            //The solution is approximate. More over, all requirements to the
+            //linearization must be satisfied in order to obtain quality result.
 
-          if((aV1[i] > aVSurf1l) && (aFact1 > 0.0))
-          {//Make close to aVSurf1l by increasing anU1
-            aDeltaV1 = Min(aDeltaV1, Abs(aV1[i]-aVSurf1l));
-          }
+            if (!StepComputing(aMatr, aV1[i], aV2[i], aDeltaV1, aDeltaV2, aStepTmp))
+            {
+              //To avoid cycling-up
+              anUexpect[i] += aStepMax;
+              aMinUexp = Min(aMinUexp, anUexpect[i]);
 
-          if((aV2[i] < aVSurf2f) && (aFact2 < 0.0))
-          {//Make close to aVSurf2f by increasing anU1
-            aDeltaV2 = Min(aDeltaV2, Abs(aV2[i]-aVSurf2f));
-          }
+              continue;
+            }
 
-          if((aV2[i] > aVSurf2l) && (aFact2 > 0.0))
-          {//Make close to aVSurf2l by increasing anU1
-            aDeltaV2 = Min(aDeltaV2, Abs(aV2[i]-aVSurf1l));
-          }
+            if (aStepTmp < aStepMin)
+              aStepTmp = aStepMin;
 
-          aDeltaU1V1 = !IsEqual(aFact1,0.0)? Abs(aDeltaV1/aFact1) : aStepMax;
-          aDeltaU1V2 = !IsEqual(aFact2,0.0)? Abs(aDeltaV2/aFact2) : aStepMax;
+            if (aStepTmp > aStepMax)
+              aStepTmp = aStepMax;
 
-          if(aDeltaU1V1 < aStepU1)
-            aStepU1 = aDeltaU1V1;
+            anUexpect[i] = anU1 + aStepTmp;
+            aMinUexp = Min(aMinUexp, anUexpect[i]);
+          }
 
-          if(aDeltaU1V2 < aStepU1)
-            aStepU1 = aDeltaU1V2;
+          anU1 = aMinUexp;
         }
 
-        if(aStepU1 < aStepMin)
-          aStepU1 = aStepMin;
-      
-        if(aStepU1 > aStepMax)
-          aStepU1 = aStepMax;
-
-        anU1 += aStepU1;
-
-        const Standard_Real aDiff = anU1 - anUl;
-        if((0.0 < aDiff) && (aDiff < aStepU1-Precision::PConfusion()))
+        if (Precision::PConfusion() >= (anUl - anU1))
           anU1 = anUl;
 
         anUf = anU1;
 
-        for(Standard_Integer i = 0; i < aNbWLines; i++)
+        for (Standard_Integer i = 0; i < aNbWLines; i++)
         {
-          if(aWLine[i]->NbPnts() != 1)
+          if (aWLine[i]->NbPnts() != 1)
             isAddedIntoWL[i] = Standard_False;
+
+          if (anU1 == anUl)
+          {//strictly equal. Tolerance is considered above.
+            anUexpect[i] = anUl;
+          }
         }
       }
 
-      for(Standard_Integer i = 0; i < aNbWLines; i++)
+      for (Standard_Integer i = 0; i < aNbWLines; i++)
       {
-        if((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
+        if ((aWLine[i]->NbPnts() == 1) && (!isAddedIntoWL[i]))
         {
           isTheEmpty = Standard_False;
           Standard_Real u1, v1, u2, v2;
@@ -2209,33 +3360,124 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
           IntPatch_Point aP;
           aP.SetParameter(u1);
           aP.SetParameters(u1, v1, u2, v2);
-          aP.SetTolerance(theTol3D);
+          aP.SetTolerance(aTol3D);
           aP.SetValue(aWLine[i]->Point(1).Value());
 
-          theSPnt.Append(aP);
+          //Check whether the added point exists.
+          //It is enough to check the last point.
+          if (theSPnt.IsEmpty() ||
+              !theSPnt.Last().PntOn2S().IsSame(aP.PntOn2S(), Precision::Confusion()))
+          {
+            theSPnt.Append(aP);
+          }
         }
-        else if(aWLine[i]->NbPnts() > 1)
+        else if (aWLine[i]->NbPnts() > 1)
         {
           Standard_Boolean isGood = Standard_True;
 
-          if(aWLine[i]->NbPnts() == 2)
+          if (aWLine[i]->NbPnts() == 2)
           {
             const IntSurf_PntOn2S& aPf = aWLine[i]->Point(1);
             const IntSurf_PntOn2S& aPl = aWLine[i]->Point(2);
-            
-            if(aPf.IsSame(aPl, Precision::Confusion()))
+
+            if (aPf.IsSame(aPl, Precision::Confusion()))
               isGood = Standard_False;
           }
+          else if (aWLine[i]->NbPnts() > 2)
+          {
+            // Sometimes points of the WLine are distributed
+            // linearly and uniformly. However, such position
+            // of the points does not always describe the real intersection
+            // curve. I.e. real tangents at the ends of the intersection
+            // curve can significantly deviate from this "line" direction.
+            // Here we are processing this case by inserting additional points
+            // to the beginning/end of the WLine to make it more precise.
+            // See description to the issue #30082.
+
+            const Standard_Real aSqTol3D = aTol3D*aTol3D;
+            for (Standard_Integer j = 0; j < 2; j++)
+            {
+              // If j == 0 ==> add point at begin of WLine.
+              // If j == 1 ==> add point at end of WLine.
+
+              for (;;)
+              {
+                if (aWLine[i]->NbPnts() >= aNbMaxPoints)
+                {
+                  break;
+                }
+
+                // Take 1st and 2nd point to compute the "line" direction.
+                // For our convenience, we make 2nd point be the ends of the WLine
+                // because it will be used for computation of the normals 
+                // to the surfaces.
+                const Standard_Integer anIdx1 = j ? aWLine[i]->NbPnts() - 1 : 2;
+                const Standard_Integer anIdx2 = j ? aWLine[i]->NbPnts() : 1;
+
+                const gp_Pnt &aP1 = aWLine[i]->Point(anIdx1).Value();
+                const gp_Pnt &aP2 = aWLine[i]->Point(anIdx2).Value();
+
+                const gp_Vec aDir(aP1, aP2);
+
+                if (aDir.SquareMagnitude() < aSqTol3D)
+                {
+                  break;
+                }
+
+                // Compute tangent in first/last point of the WLine.
+                // We do not take into account the flag "isReversed"
+                // because strict direction of the tangent is not
+                // important here (we are interested in the tangent
+                // line itself and nothing to fear if its direction
+                // is reversed).
+                const gp_Vec aN1 = aQuad1.Normale(aP2);
+                const gp_Vec aN2 = aQuad2.Normale(aP2);
+                const gp_Vec aTg(aN1.Crossed(aN2));
+
+                if (aTg.SquareMagnitude() < Precision::SquareConfusion())
+                {
+                  // Tangent zone
+                  break;
+                }
+
+                // Check of the bending
+                Standard_Real anAngle = aDir.Angle(aTg);
+
+                if (anAngle > M_PI_2)
+                  anAngle -= M_PI;
+
+                if (Abs(anAngle) > 0.25) // ~ 14deg.
+                {
+                  const Standard_Integer aNbPntsPrev = aWLine[i]->NbPnts();
+                  SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
+                                       anEquationCoeffs, i, 3, anIdx1, anIdx2,
+                                       aTol2D, aPeriod, isReversed);
+
+                  if (aWLine[i]->NbPnts() == aNbPntsPrev)
+                  {
+                    // No points have been added. ==> Exit from a loop.
+                    break;
+                  }
+                }
+                else
+                {
+                  // Good result has been achieved. ==> Exit from a loop.
+                  break;
+                }
+              } // for (;;)
+            }
+          }
 
-          if(isGood)
+          if (isGood)
           {
             isTheEmpty = Standard_False;
             isAddedIntoWL[i] = Standard_True;
-            SeekAdditionalPoints( theQuad1, theQuad2, aWLine[i]->Curve(), 
-                                  anEquationCoeffs, aNbPoints, aUSurf2f, aUSurf2l,
-                                  theTol2D, aPeriod, anArccosFactor[i], isTheReverse);
+            SeekAdditionalPoints(aQuad1, aQuad2, aWLine[i]->Curve(),
+                                 anEquationCoeffs, i, aNbPoints, 1,
+                                 aWLine[i]->NbPnts(), aTol2D, aPeriod,
+                                 isReversed);
 
-            aWLine[i]->ComputeVertexParameters(theTol3D);
+            aWLine[i]->ComputeVertexParameters(aTol3D);
             theSlin.Append(aWLine[i]);
           }
         }
@@ -2243,11 +3485,383 @@ Standard_Boolean IntCyCyTrim( const IntSurf_Quadric& theQuad1,
         {
           isAddedIntoWL[i] = Standard_False;
         }
+
+#ifdef INTPATCH_IMPIMPINTERSECTION_DEBUG
+        //aWLine[i]->Dump();
+#endif
       }
     }
   }
 
-  return Standard_True;
+
+  //Delete the points in theSPnt, which
+  //lie at least in one of the line in theSlin.
+  for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
+  {
+    for (Standard_Integer aNbLin = 1; aNbLin <= theSlin.Length(); aNbLin++)
+    {
+      Handle(IntPatch_WLine) aWLine1(Handle(IntPatch_WLine)::
+        DownCast(theSlin.Value(aNbLin)));
+
+      const IntSurf_PntOn2S& aPntFWL1 = aWLine1->Point(1);
+      const IntSurf_PntOn2S& aPntLWL1 = aWLine1->Point(aWLine1->NbPnts());
+
+      const IntSurf_PntOn2S aPntCur = theSPnt.Value(aNbPnt).PntOn2S();
+      if (aPntCur.IsSame(aPntFWL1, aTol3D) ||
+        aPntCur.IsSame(aPntLWL1, aTol3D))
+      {
+        theSPnt.Remove(aNbPnt);
+        aNbPnt--;
+        break;
+      }
+    }
+  }
+
+  //Try to add new points in the neighborhood of existing point
+  for (Standard_Integer aNbPnt = 1; aNbPnt <= theSPnt.Length(); aNbPnt++)
+  {
+    // Standard algorithm (implemented above) could not find any
+    // continuous curve in neighborhood of aPnt2S (e.g. because
+    // this curve is too small; see tests\bugs\modalg_5\bug25292_35 and _36).
+    // Here, we will try to find several new points nearer to aPnt2S.
+
+    // The algorithm below tries to find two points in every
+    // intervals [u1 - aStepMax, u1] and [u1, u1 + aStepMax]
+    // and every new point will be in maximal distance from
+    // u1. If these two points exist they will be joined
+    // by the intersection curve.
+
+    const IntPatch_Point& aPnt2S = theSPnt.Value(aNbPnt);
+
+    Standard_Real u1, v1, u2, v2;
+    aPnt2S.Parameters(u1, v1, u2, v2);
+
+    Handle(IntSurf_LineOn2S) aL2S = new IntSurf_LineOn2S();
+    Handle(IntPatch_WLine) aWLine = new IntPatch_WLine(aL2S, Standard_False);
+    aWLine->SetCreatingWayInfo(IntPatch_WLine::IntPatch_WLImpImp);
+
+    //Define the index of WLine, which lies the point aPnt2S in.
+    Standard_Integer anIndex = 0;
+
+    Standard_Real anUf = 0.0, anUl = 0.0, aCurU2 = 0.0;
+    if (isReversed)
+    {
+      anUf = Max(u2 - aStepMax, aUSurf1f);
+      anUl = Min(u2 + aStepMax, aUSurf1l);
+      aCurU2 = u1;
+    }
+    else
+    {
+      anUf = Max(u1 - aStepMax, aUSurf1f);
+      anUl = Min(u1 + aStepMax, aUSurf1l);
+      aCurU2 = u2;
+    }
+
+    const Standard_Real anUinf = anUf, anUsup = anUl, anUmid = 0.5*(anUf + anUl);
+
+    {
+      //Find the value of anIndex variable.
+      Standard_Real aDelta = RealLast();
+      for (Standard_Integer i = 0; i < aNbWLines; i++)
+      {
+        Standard_Real anU2t = 0.0;
+        if (!ComputationMethods::CylCylComputeParameters(anUmid, i, anEquationCoeffs, anU2t))
+          continue;
+
+        Standard_Real aDU2 = fmod(Abs(anU2t - aCurU2), aPeriod);
+        aDU2 = Min(aDU2, Abs(aDU2 - aPeriod));
+        if (aDU2 < aDelta)
+        {
+          aDelta = aDU2;
+          anIndex = i;
+        }
+      }
+    }
+
+    // Bisection method is used in order to find every new point.
+    // I.e. if we need to find intersection point in the interval [anUinf, anUmid]
+    // we check the point anUC = 0.5*(anUinf+anUmid). If it is an intersection point
+    // we try to find another point in the interval [anUinf, anUC] (because we find the point in
+    // maximal distance from anUmid). If it is not then we try to find another point in the
+    // interval [anUC, anUmid]. Next iterations will be made analogically.
+    // When we find intersection point in the interval [anUmid, anUsup] we try to find
+    // another point in the interval [anUC, anUsup] if anUC is intersection point and
+    // in the interval [anUmid, anUC], otherwise.
+
+    Standard_Real anAddedPar[2] = {isReversed ? u2 : u1, isReversed ? u2 : u1};
+
+    for (Standard_Integer aParID = 0; aParID < 2; aParID++)
+    {
+      if (aParID == 0)
+      {
+        anUf = anUinf;
+        anUl = anUmid;
+      }
+      else // if(aParID == 1)
+      {
+        anUf = anUmid;
+        anUl = anUsup;
+      }
+
+      Standard_Real &aPar1 = (aParID == 0) ? anUf : anUl,
+                    &aPar2 = (aParID == 0) ? anUl : anUf;
+
+      while (Abs(aPar2 - aPar1) > aStepMin)
+      {
+        Standard_Real anUC = 0.5*(anUf + anUl);
+        Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
+        Standard_Boolean isDone = ComputationMethods::
+                CylCylComputeParameters(anUC, anIndex, anEquationCoeffs, aU2, aV1, aV2);
+
+        if (isDone)
+        {
+          if (Abs(aV1 - aVSurf1f) <= aTol2D)
+            aV1 = aVSurf1f;
+
+          if (Abs(aV1 - aVSurf1l) <= aTol2D)
+            aV1 = aVSurf1l;
+
+          if (Abs(aV2 - aVSurf2f) <= aTol2D)
+            aV2 = aVSurf2f;
+
+          if (Abs(aV2 - aVSurf2l) <= aTol2D)
+            aV2 = aVSurf2l;
+
+          isDone = AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed,
+                                  Standard_True, gp_Pnt2d(anUC, aV1), gp_Pnt2d(aU2, aV2),
+                                  aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+                                  aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l,
+                                  aPeriod, aWLine->Curve(), anIndex, aTol3D,
+                                  aTol2D, Standard_False, Standard_True);
+        }
+
+        if (isDone)
+        {
+          anAddedPar[0] = Min(anAddedPar[0], anUC);
+          anAddedPar[1] = Max(anAddedPar[1], anUC);
+          aPar2 = anUC;
+        }
+        else
+        {
+          aPar1 = anUC;
+        }
+      }
+    }
+
+    //Fill aWLine by additional points
+    if (anAddedPar[1] - anAddedPar[0] > aStepMin)
+    {
+      for (Standard_Integer aParID = 0; aParID < 2; aParID++)
+      {
+        Standard_Real aU2 = 0.0, aV1 = 0.0, aV2 = 0.0;
+        ComputationMethods::CylCylComputeParameters(anAddedPar[aParID], anIndex,
+                                                  anEquationCoeffs, aU2, aV1, aV2);
+
+        AddPointIntoWL(aQuad1, aQuad2, anEquationCoeffs, isReversed, Standard_True,
+                        gp_Pnt2d(anAddedPar[aParID], aV1), gp_Pnt2d(aU2, aV2),
+                        aUSurf1f, aUSurf1l, aUSurf2f, aUSurf2l,
+                        aVSurf1f, aVSurf1l, aVSurf2f, aVSurf2l, aPeriod, aWLine->Curve(),
+                        anIndex, aTol3D, aTol2D, Standard_False, Standard_False);
+      }
+
+      SeekAdditionalPoints(aQuad1, aQuad2, aWLine->Curve(),
+                            anEquationCoeffs, anIndex, aNbMinPoints,
+                            1, aWLine->NbPnts(), aTol2D, aPeriod,
+                            isReversed);
+
+      aWLine->ComputeVertexParameters(aTol3D);
+      theSlin.Append(aWLine);
+
+      theSPnt.Remove(aNbPnt);
+      aNbPnt--;
+    }
+  }
+
+  return IntPatch_ImpImpIntersection::IntStatus_OK;
+}
+
+//=======================================================================
+//function : IntCyCy
+//purpose  : 
+//=======================================================================
+IntPatch_ImpImpIntersection::IntStatus IntCyCy(const IntSurf_Quadric& theQuad1,
+                                               const IntSurf_Quadric& theQuad2,
+                                               const Standard_Real theTol3D,
+                                               const Standard_Real theTol2D,
+                                               const Bnd_Box2d& theUVSurf1,
+                                               const Bnd_Box2d& theUVSurf2,
+                                               Standard_Boolean& isTheEmpty,
+                                               Standard_Boolean& isTheSameSurface,
+                                               Standard_Boolean& isTheMultiplePoint,
+                                               IntPatch_SequenceOfLine& theSlin,
+                                               IntPatch_SequenceOfPoint& theSPnt)
+{
+  isTheEmpty = Standard_True;
+  isTheSameSurface = Standard_False;
+  isTheMultiplePoint = Standard_False;
+  theSlin.Clear();
+  theSPnt.Clear();
+
+  const gp_Cylinder aCyl1 = theQuad1.Cylinder(),
+                    aCyl2 = theQuad2.Cylinder();
+
+  IntAna_QuadQuadGeo anInter(aCyl1,aCyl2,theTol3D);
+
+  if (!anInter.IsDone())
+  {
+    return IntPatch_ImpImpIntersection::IntStatus_Fail;
+  }
+
+  if(anInter.TypeInter() != IntAna_NoGeometricSolution)
+  {
+    if (CyCyAnalyticalIntersect(theQuad1, theQuad2, anInter,
+                                theTol3D, isTheEmpty,
+                                isTheSameSurface, isTheMultiplePoint,
+                                theSlin, theSPnt))
+    {
+      return IntPatch_ImpImpIntersection::IntStatus_OK;
+    }
+  }
+  
+  //Here, intersection line is not an analytical curve(line, circle, ellipsis etc.)
+
+  Standard_Real aUSBou[2][2], aVSBou[2][2]; //const
+
+  theUVSurf1.Get(aUSBou[0][0], aVSBou[0][0], aUSBou[0][1], aVSBou[0][1]);
+  theUVSurf2.Get(aUSBou[1][0], aVSBou[1][0], aUSBou[1][1], aVSBou[1][1]);
+
+  const Standard_Real aPeriod = 2.0*M_PI;
+  const Standard_Integer aNbWLines = 2;
+
+  const ComputationMethods::stCoeffsValue anEquationCoeffs1(aCyl1, aCyl2);
+  const ComputationMethods::stCoeffsValue anEquationCoeffs2(aCyl2, aCyl1);
+
+  //Boundaries. 
+  //Intersection result can include two non-connected regions
+  //(see WorkWithBoundaries::BoundariesComputing(...) method).
+  const Standard_Integer aNbOfBoundaries = 2;
+  Bnd_Range anURange[2][aNbOfBoundaries];   //const
+
+  if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs1, aPeriod, anURange[0]))
+    return IntPatch_ImpImpIntersection::IntStatus_OK;
+
+  if (!WorkWithBoundaries::BoundariesComputing(anEquationCoeffs2, aPeriod, anURange[1]))
+    return IntPatch_ImpImpIntersection::IntStatus_OK;
+
+  //anURange[*] can be in different periodic regions in
+  //compare with First-Last surface. E.g. the surface
+  //is full cylinder [0, 2*PI] but anURange is [5, 7].
+  //Trivial common range computation returns [5, 2*PI] and
+  //its summary length is 2*PI-5 == 1.28... only. That is wrong.
+  //This problem can be solved by the following
+  //algorithm:
+  // 1. split anURange[*] by the surface boundary;
+  // 2. shift every new range in order to inscribe it
+  //      in [Ufirst, Ulast] of cylinder;
+  // 3. consider only common ranges between [Ufirst, Ulast]
+  //    and new ranges.
+  //
+  // In above example, we obtain following:
+  // 1. two ranges: [5, 2*PI] and [2*PI, 7];
+  // 2. after shifting: [5, 2*PI] and [0, 7-2*PI];
+  // 3. Common ranges: ([5, 2*PI] and [0, 2*PI]) == [5, 2*PI],
+  //                   ([0, 7-2*PI] and [0, 2*PI]) == [0, 7-2*PI];
+  // 4. Their summary length is (2*PI-5)+(7-2*PI-0)==7-5==2 ==> GOOD.
+
+  Standard_Real aSumRange[2] = { 0.0, 0.0 };
+  Handle(NCollection_IncAllocator) anAlloc = new NCollection_IncAllocator;
+  for (Standard_Integer aCID = 0; aCID < 2; aCID++)
+  {
+    anAlloc->Reset();
+    NCollection_List<Bnd_Range> aListOfRng(anAlloc);
+    
+    aListOfRng.Append(anURange[aCID][0]);
+    aListOfRng.Append(anURange[aCID][1]);
+
+    const Standard_Real aSplitArr[3] = {aUSBou[aCID][0], aUSBou[aCID][1], 0.0};
+
+    NCollection_List<Bnd_Range>::Iterator anITrRng;
+    for (Standard_Integer aSInd = 0; aSInd < 3; aSInd++)
+    {
+      NCollection_List<Bnd_Range> aLstTemp(aListOfRng);
+      aListOfRng.Clear();
+      for (anITrRng.Init(aLstTemp); anITrRng.More(); anITrRng.Next())
+      {
+        Bnd_Range& aRng = anITrRng.Value();
+        aRng.Split(aSplitArr[aSInd], aListOfRng, aPeriod);
+      }
+    }
+
+    anITrRng.Init(aListOfRng);
+    for (; anITrRng.More(); anITrRng.Next())
+    {
+      Bnd_Range& aCurrRange = anITrRng.Value();
+
+      Bnd_Range aBoundR;
+      aBoundR.Add(aUSBou[aCID][0]);
+      aBoundR.Add(aUSBou[aCID][1]);
+
+      if (!InscribeInterval(aUSBou[aCID][0], aUSBou[aCID][1],
+                                          aCurrRange, theTol2D, aPeriod))
+      {
+        //If aCurrRange does not have common block with
+        //[Ufirst, Ulast] of cylinder then we will try
+        //to inscribe [Ufirst, Ulast] in the boundaries of aCurrRange.
+        Standard_Real aF = 1.0, aL = 0.0;
+        if (!aCurrRange.GetBounds(aF, aL))
+          continue;
+
+        if ((aL < aUSBou[aCID][0]))
+        {
+          aCurrRange.Shift(aPeriod);
+        }
+        else if (aF > aUSBou[aCID][1])
+        {
+          aCurrRange.Shift(-aPeriod);
+        }
+      }
+
+      aBoundR.Common(aCurrRange);
+
+      const Standard_Real aDelta = aBoundR.Delta();
+
+      if (aDelta > 0.0)
+      {
+        aSumRange[aCID] += aDelta;
+      }
+    }
+  }
+
+  //The bigger range the bigger number of points in Walking-line (WLine)
+  //we will be able to add and consequently we will obtain more
+  //precise intersection line.
+  //Every point of WLine is determined as function from U1-parameter,
+  //where U1 is U-parameter on 1st quadric.
+  //Therefore, we should use quadric with bigger range as 1st parameter
+  //in IntCyCy() function.
+  //On the other hand, there is no point in reversing in case of
+  //analytical intersection (when result is line, ellipse, point...).
+  //This result is independent of the arguments order.
+  const Standard_Boolean isToReverse = (aSumRange[1] > aSumRange[0]);
+
+  if (isToReverse)
+  {
+    const WorkWithBoundaries aBoundWork(theQuad2, theQuad1, anEquationCoeffs2,
+                                        theUVSurf2, theUVSurf1, aNbWLines,
+                                        aPeriod, theTol3D, theTol2D, Standard_True);
+
+    return CyCyNoGeometric(aCyl2, aCyl1, aBoundWork, anURange[1], aNbOfBoundaries,
+                              isTheEmpty, theSlin, theSPnt);
+  }
+  else
+  {
+    const WorkWithBoundaries aBoundWork(theQuad1, theQuad2, anEquationCoeffs1,
+                                        theUVSurf1, theUVSurf2, aNbWLines,
+                                        aPeriod, theTol3D, theTol2D, Standard_False);
+
+    return CyCyNoGeometric(aCyl1, aCyl2, aBoundWork, anURange[0], aNbOfBoundaries,
+                              isTheEmpty, theSlin, theSPnt);
+  }
 }
 
 //=======================================================================
@@ -2609,7 +4223,7 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
        //curvsol = anaint.Curve(i);
        aC=anaint.Curve(i);
        aLC.Clear();
-       ExploreCurve(Cy, Co, aC, 10.*Tol, aLC);
+       ExploreCurve(Co, aC, 10.*Tol, aLC);
        //
        aIt.Initialize(aLC);
        for (; aIt.More(); aIt.Next()) {
@@ -2658,7 +4272,7 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
            ptvalid = curvsol.Value(para);
            alig = new IntPatch_ALine(curvsol,Standard_False);
            kept = Standard_True;
-           //-- cout << "Transition indeterminee" << endl;
+           //-- std::cout << "Transition indeterminee" << std::endl;
          }
          if (kept) {
            Standard_Boolean Nfirstp = !firstp;
@@ -2682,104 +4296,69 @@ Standard_Boolean IntCyCo(const IntSurf_Quadric& Quad1,
 }
 //=======================================================================
 //function : ExploreCurve
-//purpose  : 
+//purpose  : Splits aC on several curves in the cone apex points.
 //=======================================================================
-Standard_Boolean ExploreCurve(const gp_Cylinder& ,//aCy,
-                             const gp_Cone& aCo,
-                             IntAna_Curve& aC,
-                             const Standard_Real aTol,
-                             IntAna_ListOfCurve& aLC)
-                             
+Standard_Boolean ExploreCurve(const gp_Cone& theCo,
+                              IntAna_Curve& theCrv,
+                              const Standard_Real theTol,
+                              IntAna_ListOfCurve& theLC)
 {
-  Standard_Boolean bFind=Standard_False;
-  Standard_Real aTheta, aT1, aT2, aDst;
-  gp_Pnt aPapx, aPx;
-  //
-  //aC.Dump();
-  //
-  aLC.Clear();
-  aLC.Append(aC);
-  //
-  aPapx=aCo.Apex();
-  //
-  aC.Domain(aT1, aT2);
-  //
-  aPx=aC.Value(aT1);
-  aDst=aPx.Distance(aPapx);
-  if (aDst<aTol) {
-    return bFind;
-  }
-  aPx=aC.Value(aT2);
-  aDst=aPx.Distance(aPapx);
-  if (aDst<aTol) {
-    return bFind;
-  }
+  const Standard_Real aSqTol = theTol*theTol;
+  const gp_Pnt aPapx(theCo.Apex());
+
+  Standard_Real aT1, aT2;
+  theCrv.Domain(aT1, aT2);
+
+  theLC.Clear();
   //
-  bFind=aC.FindParameter(aPapx, aTheta);
-  if (!bFind){
-    return bFind;
+  TColStd_ListOfReal aLParams;
+  theCrv.FindParameter(aPapx, aLParams);
+  if (aLParams.IsEmpty())
+  {
+    theLC.Append(theCrv);
+    return Standard_False;
   }
-  //
-  aPx=aC.Value(aTheta);
-  aDst=aPx.Distance(aPapx);
-  if (aDst>aTol) {
-    return !bFind;
+
+  for (TColStd_ListIteratorOfListOfReal anItr(aLParams); anItr.More(); anItr.Next())
+  {
+    Standard_Real aPrm = anItr.Value();
+
+    if ((aPrm - aT1) < Precision::PConfusion())
+      continue;
+
+    Standard_Boolean isLast = Standard_False;
+    if ((aT2 - aPrm) < Precision::PConfusion())
+    {
+      aPrm = aT2;
+      isLast = Standard_True;
+    }
+
+    const gp_Pnt aP = theCrv.Value(aPrm);
+    const Standard_Real aSqD = aP.SquareDistance(aPapx);
+    if (aSqD < aSqTol)
+    {
+      IntAna_Curve aC1 = theCrv;
+      aC1.SetDomain(aT1, aPrm);
+      aT1 = aPrm;
+      theLC.Append(aC1);
+    }
+
+    if (isLast)
+      break;
   }
-  //
-  // need to be splitted at aTheta
-  IntAna_Curve aC1, aC2;
-  //
-  aC1=aC;
-  aC1.SetDomain(aT1, aTheta);
-  aC2=aC;
-  aC2.SetDomain(aTheta, aT2);
-  //
-  aLC.Clear();
-  aLC.Append(aC1);
-  aLC.Append(aC2);
-  //
-  return bFind;
-}
-//=======================================================================
-//function : IsToReverse
-//purpose  : 
-//=======================================================================
-Standard_Boolean IsToReverse(const gp_Cylinder& Cy1,
-                            const gp_Cylinder& Cy2,
-                            const Standard_Real Tol)
-{
-  Standard_Boolean bRet;
-  Standard_Real aR1,aR2, dR, aSc1, aSc2;
-  //
-  bRet=Standard_False;
-  //
-  aR1=Cy1.Radius();
-  aR2=Cy2.Radius();
-  dR=aR1-aR2;
-  if (dR<0.) {
-    dR=-dR;
+
+  if (theLC.IsEmpty())
+  {
+    theLC.Append(theCrv);
+    return Standard_False;
   }
-  if (dR>Tol) {
-    bRet=aR1>aR2;
+
+  if ((aT2 - aT1) > Precision::PConfusion())
+  {
+    IntAna_Curve aC1 = theCrv;
+    aC1.SetDomain(aT1, aT2);
+    theLC.Append(aC1);
   }
-  else {
-    gp_Dir aDZ(0.,0.,1.);
-    //
-    const gp_Dir& aDir1=Cy1.Axis().Direction();
-    aSc1=aDir1*aDZ;
-    if (aSc1<0.) {
-      aSc1=-aSc1;
-    }
-    //
-    const gp_Dir& aDir2=Cy2.Axis().Direction();
-    aSc2=aDir2*aDZ;
-    if (aSc2<0.) {
-      aSc2=-aSc2;
-    }
-    //
-    if (aSc2<aSc1) {
-      bRet=!bRet;
-    }
-  }//if (dR<Tol) {
-  return bRet;
+
+  return Standard_True;
 }