  It is now possible to compute linear in mass-ratio terms in the
post-Newtonian (PN) expansion for compact binaries to very high orders using
black hole perturbation theory applied to various invariants. For instance, a
computation of the redshift invariant of a point particle in a circular orbit
about a black hole in linear perturbation theory gives the linear-in-mass-ratio
portion of the binding energy of a circular binary with arbitrary mass ratio.
This binding energy, in turn, encodes the system's conservative dynamics. We
give a method for extracting the analytic forms of these PN coefficients from
high-accuracy numerical data using experimental mathematics techniques, notably
an integer relation algorithm. Such methods should be particularly important
when the calculations progress to the considerably more difficult case of
perturbations of the Kerr metric. As an example, we apply this method to the
redshift invariant in Schwarzschild. Here we obtain analytic coefficients to
12.5PN, and higher-order terms in mixed analytic-numerical form to 21.5PN,
including analytic forms for the complete 13.5PN coefficient, and all the
logarithmic terms at 13PN. At these high orders, an individual coefficient can
have over 30 terms, including a wide variety of transcendental numbers, when
written out in full. We are still able to obtain analytic forms for such
coefficients from the numerical data through a careful study of the structure
of the expansion. The structure we find also allows us to predict certain
"leading logarithm"-type contributions to all orders. The additional terms in
the expansion we obtain improve the accuracy of the PN series for the redshift
observable, even in the very strong-field regime inside the innermost stable
circular orbit, particularly when combined with exponential resummation.
