Read the problem as: integrate cdfN(f(a,b,x))*pdfGamma(c,x) over the positive reals with f(x) = = a/sqrt(x)+b*sqrt(x), cdfN is the cumulative normal distribution, pdfGamma the density of the gamma distribution cdfGamma (called P in Abramowitz & Stegun). The limits of f depend on the signs of a and b, so the first argument will reach 0 or 1. Or 1/2. Now estimate where to replace it by that constant (and justify why ... but i guess for mj's solution one would have to as well). So for being 'close to 0' one switches to a cdfGamma and for the infinite tail use an adaptive integration (instead of cutting of) which stops if the integrand is decreasing and becomes small. For c = time/nu getting large (i.e. beyond 170, where most implementations will quitt - but then impl volatility is almost constant) one can switch to asymptotics (which i did without extreme care). For a=0 or b=0 one can even write down 'easy' explicite solutions in terms of hypergeometric functions, which i ignored as well as i ignored the numerical case sigma --> 0.