Representing a continuous function by a polynomial is called a polynomial fit, and the usual method is to choose coefficients that minimize an error function according to a norm, most popular squared error (gauss norm) or absolute error (chebyshev norm). Even spreadsheet calculators like MS Excel or LibreOffice Calc can calculate polynomial fits.
The substitution r=(x-y)/(x+y) has the property of self-normalization, that's why it apparently has been chosen.
Similar "tricks" to approximate transcendent functions have been developped since the beginning of numerical mathematics and can be found in many "practical" text books. Their "why" not always obvious.