There are mistakes in the formula from Novembre and Slatkin 2008. I changed the modification of S and the weighted average. Now it gives perfectly similar values than with the multinormal density

Working now when assuming that the internal nodes are weighted averages from the two descendants. Crazy speed-ups vs fitBM...

adding some function to compute the likelihood of a continuous character under BM using the pruning algorithm