Lemma of Neyman Pearson (how to get critical region) We have the Random Variable X, which is $\Gamma(p, \lambda)$ distributed with the density: 
$$f_{p,\lambda}(x)=\frac{\lambda^p}{\Gamma(p)}\cdot x^{p-1} \cdot e^{- \lambda x}$$
with $p = 10$ and $H_0: \lambda=2$ or $H_1: \lambda=4$ and $\alpha = 0.001$.
I want to apply the Lemma of Neyman Pearson which states:
Be $c > 0$ fixed and chosen in the way that $A(c) = \{x \in B: \frac{f_0(x)}{f_1(x)} \geq c \}$, such that $\mathbb{P}_{H_0}(X \in A(c)) = \alpha$.
Then the test with the region $A(c)$ among all tests with significance level $\alpha$ is the most powerful.
I am now trying to calculate $A(c)$, but got stuck. I have:
$$\int_{A(c)} f_0(x)dx = \int_{A(c)} \frac{\lambda^p}{\Gamma(p)}\cdot x^{p-1} \cdot e^{- \lambda x}dx = \alpha.$$ 
But I don't know how to get $A(c)$ from this integral...
Also $\frac{f_0(x)}{f_1(x)} = \frac{1}{1024}\cdot e^{2x}$, if that helps...
 A: For the hypotheses $H_0: \lambda = \lambda_0$ and $H_1: \lambda = \lambda_1$ you have:
$$\ln f_0(x) - \ln f_1(x) = p (\ln \lambda_0 - \ln \lambda_1) - x (\lambda_0 - \lambda_1).$$
Hence, we have:
$$\frac{d}{dx} (\ln f_0(x) - \ln f_1(x)) = \lambda_1 - \lambda_0.$$
If $\lambda_1 > \lambda_0$ then the likelihood ratio is an increasing function of $x$ so the region $A(c)$ is of the form $[x_c, \infty)$.  (Alternatively, if $\lambda_1 < \lambda_0$ then the log-likelihood ratio is a decreasing function of $x$ so the region $A(c)$ is of the form $[0, x_c]$.)  In the former case we have:
$$\begin{equation} \begin{aligned}
\alpha = \mathbb{P}(X \in A(c)|H_0) 
&= \frac{\lambda_0^p}{\Gamma(p)} \int \limits_{x_c}^\infty x^{p-1} e^{-\lambda_0 x} \ dx  \\[6pt]
&= \frac{\lambda_0}{\Gamma(p)} \int \limits_{\lambda_0 x_c}^\infty t^{p-1} e^{-t} \ dt  \\[6pt]
&= \lambda_0 \cdot \frac{\Gamma(p, \lambda_0 x_c)}{\Gamma(p)}, \\[6pt]
\end{aligned} \end{equation}$$
where the last line in this equation uses the incomplete gamma function.  Hence, the boundary value $x_c$ is the value that solves:
$$\Gamma(p, \lambda_0 x_c) = \frac{\alpha}{\lambda_0} \cdot \Gamma(p).$$
This value can be obtained numerically from the incomplete gamma function.

Application to your example: In your question you have $p=10$, $\lambda_0=2$ and $\lambda_1=4$.  Since $\lambda_1 > \lambda_0$ the likelihood ratio is an increasing function of $x$ and so you have $A(c) = [x_c, \infty)$ with the boundary value $x_c$ solving:
$$\Gamma(10, 2 x_c) = \frac{\alpha}{2} \cdot \Gamma(10).$$
