If $p$ is a polynomial with real coefficients and p(x)>0 on [0,1], then $p(x)=\sum c_{i,j} x^i(1x)^j$ with $c_{i,j}$ positive. I know this is true but but I need a proof/reference. Thanks!

See "Polynomials that are positive on an interval" by myself and Bruce Reznick, Trans. Amer. Math. Soc. 352 (2000), 46774692. We give a brief history of this problem along along with a bound for the minimum $m$ so that $p$ can be written $p = \sum_{i=1}^m c_i x^i (1x)^{mi}$ with $c_i \geq 0$. The bound depends on the minimum of $p$ on $[0,1]$ and the size of the coefficients. The best reference for this (without the bound) is probably P\'olyaSzego. 


The following extension of Keivan Karai's comment proves the result. Consider $q(x):=p(x)c$ where $c=\frac12 \min p([0,1])$. Approximate $q(x)$ as Keivan suggested: let $$ q_n(x) = \sum_{i=0}^n q(i/n) \binom ni x^i(1x)^{ni} . $$ Then $q_n$ is a polynomial of at most the same degree as $q$ (and hence $p$) and $q_n$ converges to $q$ uniformly on $[0,1]$. Since the degree is bounded, it follows that the coefficients of $Q_n$ converge to the coefficients of $q$. In particular, $\deg q_n=\deg q$ for all large enough $n$. Let $a$ be the leading coefficient of $q$ and $a_n$ the leading coefficient of $q_n$. Then $a_n\to a$ and hence the polynomial $$ r_n(x):=\frac a{a_n} q_n(x) $$ converges to $q(x)$ uniformly on $[0,1]$. Choose $n$ so large that $r_nq<c$ on $[0,1]$, then $r_n<p$ on $[0,1]$. Since $r_n$ is already represented in the desired form, it remains to represent $pr_n$. Observe that $\deg(pr_n)<\deg p$ and proceed by induction. Appendix: Why $\deg q_n\le\deg q$. Since $q_n$ is linear in $q$, it suffices to consider the case $q(x)=x^k$. Let $$ F_0(u,v) = (u+v)^n = \sum_{i=0}^n \binom ni u^i v^{ni} $$ and $$ F_m = u\cdot \frac{\partial}{\partial u} F_{m1}, \qquad m=1,2,\dots $$ Then $$ F_k(u,v) = \sum_{i=0}^n i^k \binom in u^iv^{ni}, $$ hence $q_n(x) = n^{n} F_k(x,1x)$. By induction in $k$ one sees that $F_k(u,v)$ is a linear combination of terms of the form $u^j(u+v)^{nj}$ with $j\le k$. Hence $F_k(x,1x)$ has degree at most $k$, and so does $q_n(x)$. Remark. Analyzing the coefficients of the above linear combination, one can see that, as $n\to\infty$ and $k$ fixed, the coefficient at $u^j(u+v)^{nj}$ goes to 1 and all other ones go to 0. This gives an elementary proof of convergence $q_n\to q$ (coefficientwise) avoiding the general theory. 


My first response (which I'm leaving below) was fallacious, but it was based on the idea that this method had to work, given the truth of the assertion. Since it doesn't work ..., I'll show why my interpretation of the question (I thought $n$ was implicitly bounded) is impossible. Note that Sergei Ivanov's answer, based on a comment suggestion of Keivan Karai, solves the problem when $n$ is unbounded.
When the degree $n$ is bounded, the assertion for positive polynomials implies the limiting case, for polynomials that are nonnegative in the interval. But polynomials of the form given in the question do not have zeros on the unit interval (as suggested by commenters). In particular, (x1/2)^2 cannot be represented by in this form, since all coefficients would be forced to be 0. Moreover: there cannot be any positive finite generating set for the convex cone of positive polynomials of degree $> 1$. This cone is not polyhedral, having uncountably many different extremal rays for instance in degree 2, where for each $a$ in the unit interval and $t > 0$, the ray of polynomials $t (xa)^2$ is extremal, not expressible as a convex combination of other positive polynomials. Here is a plot of the convex cone showing coefficients of positive quadratic polynomials. The three triangles bound the convex combinations of basis functions mentioned in the question. The curved side consists of coefficients of nonnegative quadratic functions with a 0 in the interval. These are not positive linear combinations of any finite collection of positive polynomials. The $n$th Bernstein polynomial expressions for quadratic polynomials, as described in Ivanov's answer, map to this space as convex cones. Here is a slice, in the quadratic coefficient space, of the images of cones of positive quadratic polynomials so obtained for $n = 2..6$. The vertices in the pictures are coefficients of approximations of quadratic polynomials that have the form $(xi/n)(x(i+1)/n$. Even though these polynomials are not positive, the coefficients in the Bernstein approximation are postiive, since they only make use of values at the $(i/n)$ points. The cone of positive quadratic $n$th Bernstein polynomials is somewhat bigger than this, but these are enough to illustrate how they converge.
That's a nice, simple representation; I wasn't aware of it, but presumably something this simple must be known. Here's a proof: The case $i+j = n$ is sufficient to get degree $n$ polynomials (perhaps that's what you meant). Here is a plot of this positive basis, for $n = 3$; the picture suggests a working strategy. Since only one of the functions is positive at each endpoint, you know what the coefficients of these functions must be. If you subtract, is it still positive? That follows from
(My proof was fallacious ). From the lemma, it follows that $q(x) = p(x)  p(0) (1x)^n + p(1) x^n$ is still positive. Since $q(x)$ is 0 at the endpoints of the interval, it is divisible by $x(1x)$. Use induction to represent the quotient; this gives the desired representation for $p$. 


Later Still: See Vicki Power's answer for a reference more suitable for your specific question. Alternatively, if you simply want a reference, this is an application of the highly nontrivial Schmüdgen theorem (a.k.a. Schmüdgen's Positivstellensatz) 1. The theorem says that if you have a compact semialgebraic set $K \subset \mathbb{R}^n$ defined by inequalities $\{g_1 \geq 0, \dots, g_s \geq 0\}$, then any positive polynomial $f$ on $K$ belongs to the cone generated by $\{g_1, \dots, g_s\}$, i.e. sums with positive coefficients of products of the $g_i$'s (including powers) together with squares of arbitrary polynomials. In your case, you have $g_1(x)=x$ and $g_2(x)=1x$, and because you're in the one variable case, you should be able to get rid of the arbitrary squares (I might come back to this answer when I have more time). 1 K. Schmüdgen, The Kmoment problem for compact semialgebraic sets. Math. Ann. 289 (1991), no. 2, 203206. (link) Added Later: Note that Schmüdgen's result is only when $fK>0$. If you relax to $fK\geq 0$, it is not true any more. More surprisingly, Stengle points out in the MR that this result does not hold if you replace $\mathbb{R}$ by an arbitrary real closed field. 


As Vicki already pointed out above, PolyaSzegö (Problems and Theorems in Analysis II. Theory of Functions. Zeros. Polynomials. Determinants. Number Theory. Geometry.) is a good reference for this particular result. See §6, Problem 49 in that book: By the Fundamental Theorem of Algebra you can restrict to quadratic polynomials. For lines it is trivial. For parabolas opening down it is easy. The only nontrivial case happens for parabolas opening up. Parabolas opening up whose minimum lies above the interval are also just a little exercise. The only case remaining is parabolas opening up without zero on the real line. Now determine a polynomial $h \in \mathbb Z[K,L,A,B,C]_ 4$ such that for all $a,b,c\in\mathbb R$ and all $\ell\in\mathbb N_{\ge2}$ one has $$aX^2+bX+c=\sum_{k=0}^\ell\frac{(\ell2)!}{k!(\ellk)!}h(k,\ell,a,b,c)X^k(1X)^{\ellk}.$$ Show that each $f\in\mathbb R[X]_ 2$ with $f>0$ on $\mathbb R$ has the desired representation by comparing the discriminants of $f=aX^2+bX+c\in\mathbb R[X]_ 2$ and $h(K,\ell,a,b,c)\in\mathbb Z[K]_ 2$ for $a,b,c\in\mathbb R$ and large $\ell\in\mathbb N_{\ge2}$. This approach unfortunately does not carry over for the numerous generalizations of the theorem. 

