Modified Chebyshev-Picard Iteration methods are presented for solving high precision, long-term orbit propagation problems. Fusing Chebyshev polynomials with the classical Picard iteration method, the proposed methods iteratively refine an orthogonal function approximation of the entire state trajectory, in contrast to traditional, step-wise, forward integration methods. Numerical results demonstrate that for orbit propagation problems, the presented methods are comparable to or superior to a state-of-the-art 12th order Runge-Kutta-Nystrom method in a serial processor as measured by both precision and efficiency. We have found revolutionary long solution arcs with more than eleven digit path approximations over one to three lower-case Earth orbit periods, multiple solution arcs can be patched continuously together to achieve very long-term propagation, leading to more than ten digit accuracy with built-in precise interpolation. Of revolutionary practical promise to much more efficiently solving high precision, long-term orbital trajectory propagation problems is the observation that the presented methods are well suited to massive parallelization because computation of force functions along each path iteration can be rigorously distributed over many parallel cores with negligible cross communication needed.