Robust two lines/segments intersection point in 2D
If you choose a formula with parameter $t$ you'll get the $P_x$ value as a result of a number of operations - additions, multiplications and one division:

$$t = \frac{(x_1-x_3)(y_3-y_4)-(y_1-y_3)(x_3-x_4)}{(x_1-x_2)(y_3-y_4)-(y_1-y_2)(x_3-x_4)}$$
$$P_x = x_1 + t(x_2-x_1)$$

The number $x_2=72057594037954921$ contains 17 significant decimal digits - so it can't be represented exactly using double-precision floating numbers. Apparently, a small error is introduced during some operation (most probably - the division), and this error is propagated to the result. This is a form of the Round-off error.

This bunch of calculations doesn't happen when you choose a formula with parameter $u$, because the numerator there is equal to zero.

Conclusion: if you need to work with numbers, containing a big number of significant digits - use floating numbers with bigger precision, or may be - with arbitrary precision. An example: mpmath.



Another option - use rational numbers. For example, kernels, based on rational numbers, can be found in the computational geometry library CGAL. Developers of this library have invested a lot of efforts to solve this robustness problem - please see this page for more information.

If you want to go deeper into this research area, called Exact Geometric Computation (EGC), you can start from this review. You'll need a good knowledge of Linear Algebra and Linear Systems.