diff --git a/src/IntAna/IntAna_QuadQuadGeo.cxx b/src/IntAna/IntAna_QuadQuadGeo.cxx
index af70c434c5..e90ee8896c 100644
--- a/src/IntAna/IntAna_QuadQuadGeo.cxx
+++ b/src/IntAna/IntAna_QuadQuadGeo.cxx
@@ -977,6 +977,8 @@ void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
       gp_Dir DirCyl = Cyl1.Position().Direction();
       Standard_Real ProjP2OnDirCyl1=gp_Vec(DirCyl).Dot(gp_Vec(P1,P2t));
       
+      //P2 is a projection the location of the 2nd cylinder on the base
+      //of the 1st cylinder
       P2.SetCoord(P2t.X() - ProjP2OnDirCyl1*DirCyl.X(),
                   P2t.Y() - ProjP2OnDirCyl1*DirCyl.Y(),
                   P2t.Z() - ProjP2OnDirCyl1*DirCyl.Z());
@@ -987,7 +989,7 @@ void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
         typeres=IntAna_Empty;
         nbint=0;
       }
-      else if(DistA1A2>(R1pR2))
+      else if((R1pR2 - DistA1A2) <= RealSmall())
       {
         //-- 1 Tangent line -------------------------------------OK
         typeres=IntAna_Line;
@@ -1005,32 +1007,85 @@ void IntAna_QuadQuadGeo::Perform(const gp_Cylinder& Cyl1,
         typeres=IntAna_Line;
         nbint=2;
         dir1=DirCyl;
-        gp_Vec P1P2(P1,P2);
-        gp_Dir DirA1A2=gp_Dir(P1P2);
-        gp_Dir Ortho_dir1_P1P2 = dir1.Crossed(DirA1A2);
         dir2=dir1;
-        Standard_Real Alpha=0.5*(R1*R1-R2*R2+DistA1A2*DistA1A2)/(DistA1A2);       
 
-        //Standard_Real Beta = Sqrt(R1*R1-Alpha*Alpha);
-        Standard_Real anSqrtArg = R1*R1-Alpha*Alpha;
-        Standard_Real Beta = (anSqrtArg > 0.) ? Sqrt(anSqrtArg) : 0.;
+        const Standard_Real aR1R1 = R1*R1;
+
+        /*
+                      P1
+                      o
+                    * | *
+                  * O1|   *
+              A o-----o-----o B
+                  *   |   *
+                    * | *
+                      o
+                      P2
+
+          Two cylinders have axes collinear. Therefore, problem can be reformulated as
+          to find intersection point of two circles (the bases of the cylinders) on
+          the plane: 1st circle has center P1 and radius R1 (the radius of the
+          1st cylinder) and 2nd circle has center P2 and radius R2 (the radius of the
+          2nd cylinder). The plane is the base of the 1st cylinder. Points A and B
+          are intersection point of these circles. Distance P1P2 is equal to DistA1A2.
+          O1 is the intersection point of P1P2 and AB segments.
+
+          At that, if distance AB < Tol we consider that the circles are tangent and
+          has only one intersection point.
+
+            AB = 2*R1*sin(angle AP1P2).
+          Accordingly, 
+            AB^2 < Tol^2 => 4*R1*R1*sin(angle AP1P2)^2 < Tol*Tol.
+        */
+
         
-        if((Beta+Beta)<Tol)
+        //Cosine and Square of Sine of the A-P1-P2 angle
+        const Standard_Real aCos = 0.5*(aR1R1-R2*R2+DistA1A2*DistA1A2)/(R1*DistA1A2);
+        const Standard_Real aSin2 = 1-aCos*aCos;
+
+        const Standard_Boolean isTangent =((4.0*aR1R1*aSin2) < Tol*Tol);
+
+        //Normalized vector P1P2
+        const gp_Vec DirA1A2((P2.XYZ() - P1.XYZ())/DistA1A2); 
+
+        if(isTangent)
         {
+          //Intercept the segment from P1 point along P1P2 direction
+          //and having |P1O1| length 
           nbint=1;
-          pt1.SetCoord( P1.X() + Alpha*DirA1A2.X()
-                       ,P1.Y() + Alpha*DirA1A2.Y()
-                       ,P1.Z() + Alpha*DirA1A2.Z());
+          pt1.SetXYZ(P1.XYZ() + DirA1A2.XYZ()*R1*aCos);
         }
         else
-        { 
-          pt1.SetCoord( P1.X() + Alpha*DirA1A2.X() + Beta*Ortho_dir1_P1P2.X(),
-                        P1.Y() + Alpha*DirA1A2.Y() + Beta*Ortho_dir1_P1P2.Y(),
-                        P1.Z() + Alpha*DirA1A2.Z() + Beta*Ortho_dir1_P1P2.Z());
-          
-          pt2.SetCoord( P1.X() + Alpha*DirA1A2.X() - Beta*Ortho_dir1_P1P2.X(),
-                        P1.Y() + Alpha*DirA1A2.Y() - Beta*Ortho_dir1_P1P2.Y(),
-                        P1.Z() + Alpha*DirA1A2.Z() - Beta*Ortho_dir1_P1P2.Z());
+        {
+          //Sine of the A-P1-P2 angle (if aSin2 < 0 then isTangent == TRUE =>
+          //go to another branch)
+          const Standard_Real aSin = sqrt(aSin2);
+
+          //1. Rotate P1P2 to the angle A-P1-P2 relative to P1
+          //(clockwise and anticlockwise for getting
+          //two intersection points).
+          //2. Intercept the segment from P1 along direction,
+          //determined in the preview paragraph and having R1 length
+          const gp_Dir  &aXDir = Cyl1.Position().XDirection(),
+                        &aYDir = Cyl1.Position().YDirection();
+          const gp_Vec  aR1Xdir = R1*aXDir.XYZ(),
+                        aR1Ydir = R1*aYDir.XYZ();
+
+          //Source 2D-coordinates of the P1P2 vector normalized
+          //in coordinate system, based on the X- and Y-directions
+          //of the 1st cylinder in the plane of the 1st cylinder base
+          //(P1 is the origin of the coordinate system).
+          const Standard_Real aDx = DirA1A2.Dot(aXDir),
+                              aDy = DirA1A2.Dot(aYDir);
+
+          //New coordinate (after rotation) of the P1P2 vector normalized.
+          Standard_Real aNewDx = aDx*aCos - aDy*aSin,
+                        aNewDy = aDy*aCos + aDx*aSin;
+          pt1.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
+
+          aNewDx = aDx*aCos + aDy*aSin;
+          aNewDy = aDy*aCos - aDx*aSin;
+          pt2.SetXYZ(P1.XYZ() + aNewDx*aR1Xdir.XYZ() + aNewDy*aR1Ydir.XYZ());
         }
       }
       else if(DistA1A2>(RmR-Tol))
diff --git a/tests/bugs/moddata_3/bug25782_1 b/tests/bugs/moddata_3/bug25782_1
new file mode 100755
index 0000000000..af4822b215
--- /dev/null
+++ b/tests/bugs/moddata_3/bug25782_1
@@ -0,0 +1,104 @@
+puts "========"
+puts "OCC25782"
+puts "========"
+puts ""
+######################################################
+# The result of intersection between two cylinders is incorrect
+# Algorithm must find one curves only
+######################################################
+
+set GoodNbCurv 1
+
+restore [locate_data_file bug25782_fz19.brep] b1
+restore [locate_data_file bug25782_fz53.brep] b2
+
+mksurface s1 b1
+mksurface s2 b2
+
+intersect res s1 s2
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  
+  copy res res_1
+}
+
+set ic 1
+set AllowRepeate 1
+while { $AllowRepeate != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeate 0
+  } else {
+    dlog reset
+    dlog on
+    xdistcs res_$ic s1 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    dlog reset
+    dlog on
+    xdistcs res_$ic s2 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    incr ic
+  }
+}
+  
+if {[expr {$ic - 1}] == $GoodNbCurv} {
+  puts "OK: Curve Number is good!"
+} else {
+  puts "Error: Curves Number is bad!"
+}
+
+set log [bopcurves b1 b2]
+
+regexp {Tolerance Reached=+([-0-9.+eE]+)\n+([-0-9.+eE]+)} ${log} full Toler NbCurv
+set MaxTol 1.e-7
+if {${Toler} > ${MaxTol}} {
+   puts "Error: Tolerance is too big!"
+}
+
+if {$NbCurv != $GoodNbCurv} {
+  puts "Error: Curves Number is bad!"
+}
+
+for {set i 1} {$i <= ${NbCurv}} {incr i} {
+  bounds c_$i U1 U2
+  
+  if {[dval U2-U1] < 1.0e-9} {
+    puts "Error: Wrong curve's range!"
+  }
+
+  dlog reset
+  dlog on
+  xdistcs c_$i s1 U1 U2 10
+  set Log2 [dlog get]
+  set List2 [split ${Log2} {TD= \t\n}]
+  set Tolerance 1.0e-7
+  set Limit_Tol 1.0e-7
+  set D_good 0.
+  checkList ${List2} ${Tolerance} ${D_good} ${Limit_Tol}
+
+  dlog reset
+  dlog on
+  xdistcs c_$i s2 U1 U2 10
+  set Log2 [dlog get]
+  set List2 [split ${Log2} {TD= \t\n}]
+  set Tolerance 1.0e-7
+  set Limit_Tol 1.0e-7
+  set D_good 0.
+  checkList ${List2} ${Tolerance} ${D_good} ${Limit_Tol}
+}
diff --git a/tests/bugs/moddata_3/bug25782_2 b/tests/bugs/moddata_3/bug25782_2
new file mode 100755
index 0000000000..e889d35d24
--- /dev/null
+++ b/tests/bugs/moddata_3/bug25782_2
@@ -0,0 +1,68 @@
+puts "========"
+puts "OCC25782"
+puts "========"
+puts ""
+######################################################
+# The result of intersection between two cylinders is incorrect
+######################################################
+
+set GoodNbCurv 2
+
+cylinder s1 0 0 0 12 35 47 5
+cylinder s2 3 2 8 12 35 47 4
+
+set bug_info [intersect res s1 s2]
+
+set che [whatis res]
+set ind [string first "3d curve" $che]
+if {${ind} >= 0} {
+  #Only variable "res" exists
+  
+  copy res res_1
+}
+
+if {[llength ${bug_info}] != $GoodNbCurv} {
+  puts "Error: The result of intersection between two cylinders is incorrect"
+}
+
+set Tolerance 1.e-7
+set D_good 0.
+set Limit_Tol 1.0e-7
+
+set ic 1
+set AllowRepeate 1
+while { $AllowRepeate != 0 } {
+  set che [whatis res_$ic]
+  set ind [string first "3d curve" $che]
+  if {${ind} < 0} {
+    set AllowRepeate 0
+  } else {
+    if { [regexp {\*\*\nLine} [dump res_$ic]] } {
+      #puts "OK : Correct intersection"
+    } else {
+      puts "Error : Bad intersection"
+    }
+  
+    dlog reset
+    dlog on
+    xdistcs res_$ic s1 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    dlog reset
+    dlog on
+    xdistcs res_$ic s2 0 100 10
+    set Log1 [dlog get]
+    set List1 [split ${Log1} {TD= \t\n}]
+    set Tolerance 1.0e-7
+    set Limit_Tol 1.0e-7
+    set D_good 0.
+    checkList ${List1} ${Tolerance} ${D_good} ${Limit_Tol}
+    
+    incr ic
+  }
+}