Add benchmark for incomplete gamma
Move coefficient for Gauss-Legendre integration to the top level
Use integration by quadrature for large p in incomplete gamma

Normal approximation have poor precision for moderately large p
(around 2e-6 for p=1000). Approximation which uses Gauss-Legendre
integration gives better precision. It's however deteriorate
as `p' grows so for p>2e5 we switch to the normal approximation