The Cephes implementation of log1p is described in the documentation as
Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
where P(x) and Q(x) are two polynomials (with Q(x) a monic polynomial.
This is partially the Maclaurin series for log(1 + x), but it's not immediately clear how they get P(x) and Q(x).
Is there some name for this form of approximation or resources on how to replicate the procedure to get the coefficients?
This is a particular form of a rational approximation, sometimes referred to as the Cody & Waite form. The Cody & Waite form is a form of constrained approximation that can be used with both polynomial and rational approximations and is advantageous for numerical reasons.
Approximation accuracy is improved when the leading coefficient is exactly representable in a given floating-point format. So instead of using an approximation such as sinh(x) ≈ a0x + a1x3 + ... + anx2n+1, one forces the leading coefficient to 1: sinh(x) ≈ x + a1x3 + ... + anx2n+1, since the function argument x is considered exact for the purposes of constructing math library implementations.
This arrangement is called the Cody & Waite form because it was used quite extensively in the following book which was a standard reference for authors of math libraries throughout much of the 1980s and 1990s:
William J. Cody and William Waite, Software Manual for the Elementary Functions, Englewood Cliffs, NJ: Prentice Hall 1980.
However, initial uses of this arrangement can be traced back to Hirondo Kuki, who used it for implementations of asin() and sinh(), for example. See:
Hirondo Kuki, Mathematical Functions. A Descriptions of the Center's 7094 Fortran II Mathematical Function Library, University of Chicago Computation Center Report, February 1966.
The author of the Cephes library, Stephen L. Moshier, used minimax approximations pretty much exclusively. These are approximations that minimize the maximum error over the approximation interval: approximations optimized according to the L∞ norm. For a polynomial minimax approximation Pn(x) of degree n the graph of the error function d(x) = P(x) - f(x) has characteristic properties: It wiggles back and forth across the x-axis, and on the interval of approximation reaches maxima at exactly n+2 points, where all the maxima have exactly the same magnitude but adjacent maxima alternate in their sign (equioscillation property).
This leads to a set of n+2 equations in n+2 unknowns P(x) - f(x) ± d = 0. In 1934 the Soviet mathematician Evgeny Yakovlevich Remez devised an iterative numerical algorithm for minimizing d commonly called the Remez exchange algorithm, where each iteration requires solving the system of linear equations:
E. Remes, "Sur un procédé convergent d'approximations successives pour déterminer les polynomes d'approximation, Comptes rendus hebdomadaires des séances de l'Académie des sciences, 198 (1934), pp. 2063-2065 (scan online).
This algorithm is still in use today and has been extended for use with minimax rational approximations, where the error function d(x) = Pm(x)/Qn(x) - f(x) has m+n+2 maxima, again of like magnitude but alternating signs. Moshier wrote a book describing the design underlying the Cephes library, which includes an implementation of the Remez algorithm:
Stephen L. B. Moshier, Methods and Programs for Mathematical Functions, Chichester: Ellis Horwood 1989.
As the book points out, a rational approximation in Cody & Waite form, such as f(x) ≈ x + x3 P(x2) / Q(x2), leads to a system of equations of the form x3 P(x2) - (f(x) - x ± d) Q(x2) = 0. The specific form of the log1p(x) core approximation in the question extends this by constraining the leading two coefficients.
While it is certainly possible to construct one's own implementation of the Remez algorithm, creating a robust implementation is not trivial in my experience. If one does not want to start from scratch, one might consider using Moshier's Remez implementation in the Netlib repository as a potential starting point for ones efforts. From the file remesf.c one can see that the specific form of the core approximation for log1p() is handled as a special case.
The free open-source Sollya tool has excellent facilities for generating minimax polynomial approximations, but last I checked it could not generate minimax rational approximations. One may want to look to commercial products like Mathematica and Maple for that (I have used neither).