Convert polynomial to Chebyshev

Let's say you have a polynomial:

f[x_] := 1 - x + 2 x^2 - 5 x^4
n = 5; (* Polynomial order + 1 *)

Now you have to solve a system of linear equations:

a = Transpose @ PadRight @ Table[CoefficientList[ChebyshevT[k - 1, x], x], {k, n}];
c = LinearSolve[a, CoefficientList[f[x], x]]
(* {1/8, -1, -3/2, 0, -5/8} *)

Sum[c[[k]] ChebyshevT[k - 1, x], {k, n}] // Expand
(* 1 - x + 2 x^2 - 5 x^4 *)

If getting the numerical values of coefficients is sufficient, you can use the discrete Fourier transform, which is an extremely fast method.

$\qquad T_n(x) = \cos(n \arccos x)$

c = 2 MapAt[#/2 &, #, 1]/Sqrt[n] & @ FourierDCT @ f @ Cos[π Range[.5, n]/n]
(* {0.125, -1., -1.5, -9.93014*10^-17, -0.625} *)

Sum[c[[k]] ChebyshevT[k - 1, x], {k, n}] // Expand // Chop
(* 1. - 1. x + 2. x^2 - 5. x^4 *)

P.S. The last method can be used for the Chebyshev expansion of an arbitrary function (there can be a problem with Gibbs oscillations, but that is another story).

Another possibility is to use Salzer's algorithm (previously used here) to perform the basis conversion:

poly = 1 - x + 2 x^2 - 5 x^4;
c = CoefficientList[poly, x];
n = Length[c] - 1; Remove[a];
a[0, 2] = c[[n - 1]] + c[[n + 1]]/2;
a[1, 2] = c[[n]]; a[2, 2] = c[[n + 1]]/2;
Do[
   a[0, k + 1] = c[[n - k]] + a[1, k]/2;
   a[1, k + 1] = a[0, k] + a[2, k]/2;
   Do[
      a[m, k + 1] = (a[m + 1, k] + a[m - 1, k])/2,
      {m, 2, k - 1}];
   a[k, k + 1] = a[k - 1, k]/2;
   a[k + 1, k + 1] = a[k, k]/2,
   {k, 2, n - 1}];

ccof = Table[a[m, n], {m, 0, n}]
   {1/8, -1, -3/2, 0, -5/8}

The matrix implicitly used for the basis conversion in ybeltukov's answer actually has a nice closed form. You can use this to perform the conversion directly, without having to use LinearSolve[]:

pcmat[n_Integer?NonNegative] := SparseArray[{{j_, k_} /; j <= k && EvenQ[k - j] :> 
      Binomial[k - 1, Quotient[k - 1, 2] - Quotient[j - 1, 2]]/
      2^(k - Boole[j > 1] - 1)}, {n + 1, n + 1}]

Test:

pcmat[4].CoefficientList[1 - x + 2 x^2 - 5 x^4, x]
   {1/8, -1, -3/2, 0, -5/8}

---

Having the explicit basis change matrix available allows us to study the numerical stability of the conversion. In particular, we can do an analysis similar to the one in this answer, and look at the growth of the condition number of the basis change matrix for various degrees:

Table[Ceiling[Log10[LinearAlgebra`MatrixConditionNumber[pcmat[k]]]], {k, 15}]
   {0, 1, 1, 2, 2, 3, 3, 3, 4, 4, 5, 5, 5, 6, 6}

This says, for instance, that you can lose up to $6$ significant figures when converting a degree-$15$ polynomial in the monomial basis to Chebyshev form.

---

For completeness, here's the routine for the inverse matrix (for converting from the Chebyshev basis to the monomial basis):

cpmat[n_Integer?NonNegative] := 
      SparseArray[{{1, k_} /; OddQ[k] :> (-1)^Quotient[k, 2],
                   {j_, k_} /; 1 < j <= k :> Cos[(k - j) π/2] 2^(j - 2)
                   Binomial[(k + j)/2 - 2, (k - j)/2] (k - 1)/(j - 1)},
                  {n + 1, n + 1}]