The symmetric tensor rank approximation problem (STA) consists in computing the best low rank approximation of a symmetric tensor. We describe a Riemannian Newton iteration with trust region scheme for the STA problem. We formulate this problem as a Riemannian optimization problem by parameterizing the constraint set as the Cartesian product of Veronese manifolds. We present an explicit and exact formula for the gradient vector and the Hessian matrix of the method, in terms of the weights and points of the low rank approximation and the symmetric tensor to approximate, by exploiting the properties of the apolar product. We introduce a retraction operator on the Veronese manifold. The Newton Riemannian iterations are performed for best low rank approximation over the real or complex numbers. Numerical experiments are implemented to show the numerical behavior of the new method first against perturbation, to compute the best rank-1 approximation and the spectral norm of a symmetric tensor, and to compare with some existing state-of-the-art methods.