RNG for truncated distributions

[bnicenboim]

@nschiett

real normal_lb_ub_rng(real mu, real sigma, real lb, real ub) {
    real p1 = normal_cdf(lb, mu, sigma);  // cdf with lower bound
    real p2 = normal_cdf(ub, mu, sigma);  // cdf with upper bound
    real u = uniform_rng(p1, p2);
    return (sigma * inv_Phi(u)) + mu;  // inverse cdf 
}

I started to work with this, and one problem that I notice is that the normal cdf overflows quite easily. Is there another alternative? Will there be truncated rngs built-in in stan in the near future?