Interesting q , Please see below:

1 The sin(nx )  and cos(nx)  terms are orthogonal so cross terms will be zero.

2 Consider a simple case f(x)=cos(x) = ( ejx +e-jx )/2

3  cos2(x)  =  (1+cos(2x) )/2 = ( ejx + e-jx)^2 /4 

4  So in general each cos sin term to power of n is binomial expansion of euler and cross product terms are zero because of 1 .it should be possible to obtain an algorithimic solution

Cheers

Ramesh gives the correct hint.

More generally if f is a trigonometric polynomial then it is enough to make a convolution m times with itsself: this can be done in a fast way using Fast Fourier Transforms (if the degree is not trivial).

For a more general (L^1) function, you just use a truncation argument, by reducing to a trigonomatric polynomial and then you use the FFT for convolution as before. One has be careful on the truncation since it depends on how fast the k-th Fourier coefficient of f goes to zero and the latter depends on the regularity of f. For instance if f is periodic and q times differetiable then a_k will go to zero faster than 1/k^q

Best regards 

It the fourier coefficients of f(x) are also unknowns (assumed variables), then is it possible to express the fourier coefficients of (f(x))^m in terms of these assumed variables?

We can try to use convolution. Let F denote the Fourier transform operator.Then

F { f * g } = F {f} · F {g},where · denotes point-wise multiplication and * denotes convolution. It also works the other way around:

F { f · g } = F { f } * F { g }. In particular F(f^2)= F(f)*F(f).

Let me explain the words about convolution explicitly. Assume f and g are square integrable functions on [0, 2\pi] (i.e. elements of L2 on corresponding interval),    \sum_n aneint   and   \sum_n bneint  are their Fourier series (in complex form). Then the Fourier coefficient cn of function fg  can be expressed as follows:

cn = \sum_{k=-\infty}^\infty akbn-k

I agree with Vladimir and a fast computation can be done by truncating and using the Fast Fourier Transform

Best regards Stefano

Why don't use you the convolution formula for Fourier series it is similar as for Fourier transform?