Let $f:\mathcal{R}^+\to \mathcal{R}^+$ a continuous distribution of the type $$ f(x|a_0,a_1,\cdots,a_n)=\frac{1}{k}(a_0+\sum_{i=1} a_i x^i), $$ where $a_0,a_1,\cdots,a_n$ are given parameters and $k$ the normalizing constant, that is $$ k=\int_0^{\infty} f(x|a_0,a_1,\cdots,a_n)dx $$
My questions is, what is the most efficient way to simulate from this distribution?
My first thought: since this is a polynomial, the c.d.f. is given by the integral $F(x|a_0,a_1,\cdots,a_n)=\int_{0}^x f(t|a_0,a_1,\cdots,a_n)dt$, which might be computationally fast (?), so we could use the inverse transform sampling (https://en.wikipedia.org/wiki/Inverse_transform_sampling). However, thinking of thousands of values, this might become (time) inefficient.
are there other alternatives? thanks.