How is optimal coordinates change chosen for Chebyshev expansion? I'm looking into SLATEC implementation of Bessel function $J_0$ computation (readable in C in GSL), namely at its part for arguments in interval $(0,4)$. There a Chebyshev expansion is used, but the input variable is transformed in non-obvious way. As I've understood, their coefficients are computed using this formula:
$$c_m=\frac2\pi\int_{-1}^1\frac{T_m(x)J_0\left(\sqrt{8(x+1)}\right)}{\sqrt{1-x^2}}dx,$$
and this change of coordinates from $x$ to $\sqrt{8(x+1)}$ is meant to transform from working range of Chebyshev polynomials $(-1,1)$ to the interesting range of Bessel function arguments $(0,4)$. But this change seems somewhat complicated. Naively, I'd rather start from something like $2(x+1)$, which also maps $(-1,1)$ to $(0,4)$. But when I tried it, I found that it gives much worse precision with $13$ terms of the expansion (max absolute error $6\times10^{-11}$ vs $3\times10^{-16}$ in SLATEC).
So I wonder now, how was that transformation, which is actually used in the implementation, found? Can it be proved to be optimal for the precision wanted?
 A: Let $s(x)$ be the substitution, i.e. we are looking for an approximation
$$
J_0(s(x)) \approx P(x) = \frac{c_0}{2}T_0(x) + \sum_{m=1}^{M} c_m T_m(x), 
\quad s(-1) = 0,\, s(1) = 4.
$$
with
$$
c_m = \frac{2}{\pi}\int_{-1}^1 \frac{T_m(x) J_0(s(x))}{\sqrt{1-x^2}} dx.
$$
Let's take such $s(x)$ so $J_0(s(x))$ is a polynomial $Q(x)$ of degree $M$ or less. Now $J_0(s(x))$ is approximated by $P(x)$ exactly. That $s(x)$ would be
$$
s(x) = J_0^{-1}(Q(x)).
$$
Now the problem is not about optimal $s(x)$ (which there are infinitely many) but about simple $s(x)$ of that form. It should not at least involve any special functions.
By taking
$$
J_0(x) = 1 - \frac{x^2}{4} + O(x^4)
$$
and neglecting $O(x^4)$ term we can approximate $J_0^{-1}(v)$ as
$$
J_0^{-1}(v) \approx 2\sqrt{1-v}.
$$
So
$$
s(x) = 2\sqrt{1-Q(x)}\\
s(-1) = 2 \sqrt{1 - Q(-1)} = 0\\
s(1) = 2 \sqrt{1 - Q(1)} = 4\\
Q(-1) = 1, \quad Q(1) = -3.
$$
A linear $Q(x)$ would be
$$
Q(x) = -1 - 2x.
$$
so
$$
s(x) = 2 \sqrt{1 + 1 + 2x} = 2 \sqrt{2 + 2x}.
$$
That's exactly the substitution from GSL.
