fun pwr_of_two 0 = 1 | pwr_of_two n = 2 * pwr_of_two (n-1);

fun sum_over(from, to, f) = if from > to then 0.0 else f(from) + sum_over(from+1, to, f);

fun trapzStep(f, a:real, b:real, n:int, prevS:real) = if n = 1 then 0.5*(b-a)*(f a + f b) else let
         val n_extra:int = pwr_of_two (n-1);
         val step:real = (b-a) / real(n_extra);
         val sum:real = sum_over(1, n_extra, fn i => f(a + (real i + 0.5)*step));
    in 0.5 * (prevS + (b-a) * sum / real n_extra) end;
         
fun integrate(f, a, b, eps) = 
  let fun inner(prev,n) = 
     let val next = trapzStep(f,a,b,n+1,prev) in
       if abs(next - prev) < eps then next else inner(next,n+1) 
     end
  in inner(0.0,0) end;

fun t_pdf(x,a) = gamma((a+1)/2)/gamma(a/2)/sqrt(a*PI) * Math.pow((1+x*x/a),(-(a+1)/2))