Major improvement: the code is about 1.8x faster, due to pre-computing the cummulative density on the cpa array in ptable. Look for code wrapped in the pragma named CDF_CPA.