Is there a way to find the intersection of boundaries of two cylinders without the user having to struggle with equations and, so to speak, reinvent the wheel? (I am particularly interested in intersection of infinite cylinders so that there is no problem with dealing with cylinder bases)
reg does not seem to be of any use - it can't be plotted, RandomPoint does not work on it, it produces wrong results like RegionMember giving False for point {1, 0, 1/905 (590 - Sqrt[469551])} which is in fact member of the region.
c1 = Cylinder[2 {{0, 0, -2}, {0, 0, 2}}];
c2 = Cylinder[{-2 {1, 1, 1}, 2 {1, 1 - 1/5, 1}}, 4/5];
Graphics3D[{c1, c2}]
reg=RegionIntersection @@ RegionBoundary /@ {c1, c2}
