BUG: special: Fix asymptotic behavior of rgamma and reimplement poch

The rgamma function in Cephes underflowed too early to 0, 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