BUG: special: Fix asymptotic expansion in rgamma and reimplement poch

The rgamma function in Cephes overflowed too early to inf, fix this by using lgam to evaluate the tail.
The behavior at large negative values is OK, the function is beyond float range (except at integers).

Reimplement the pochhammer symbol. The naive definition via gamma functions suffers from over/underflows and from hitting poles of the gamma functions which cancel out in several common cases.

Fixes gh-1437