"Used to" is sort of an overstatement, IMO. The use of Taylor series here is remarkable because they need to generate coefficients on the fly. Chebyshev and Remez have been around for quite a long time, since the early 1930's, and people doing numerics in a serious fashion have generally used them for polynomial approximation since the birth of computers (except in unique circumstances like the paper you cited).
The new engineering that has come recently with Sollya and similar solvers is explicit design around floating point operations. Chebyshev and Remez use the real numbers, but Sollya uses floats in its minimax function, and the coefficients you get are actually somewhat different due to rounding. Fast LP solvers have also enabled this approach more generally.
It should be said that Sollya doesn't _really_ use floats. It restricts its coefficients to rationals that are (mostly) representable by floats, but the minimax is still run in full precision. Which means you can often beat it by brute force or similar.
Yes, that is a good clarification - it uses wider precision internally, but it takes floating point restrictions and operations into account. If you otherwise used Remez, you would have to just quantize the coefficients blindly and then tweak (probably manually) if something was off.
Shamelessly plugging, you can sort of see the old-school process here (using integers and mixed-precision fixed point is harder with Sollya): https://specbranch.com/posts/faster-div8/
The coefficients and constants in both cases are substantially higher than what you get from Remez, to allow for truncation. When you have quantized operands, they are quite a bit higher. The same goes for approximations generated by Sollya - the optimal coefficients are relatively far from what Remez would tell you to do because of the error you get from rounding.
The new engineering that has come recently with Sollya and similar solvers is explicit design around floating point operations. Chebyshev and Remez use the real numbers, but Sollya uses floats in its minimax function, and the coefficients you get are actually somewhat different due to rounding. Fast LP solvers have also enabled this approach more generally.