I’ve stumbled upon a little trick to compute jenks breaks faster than with the classInt package, just be sure to use n+1 instead of n as the breaks are computed a little bit differently. That is to say, if you want 5 breaks, n=6, no biggie there.
For more on the Bayesian Analysis of Macroevolutionary Mixtures see BAMMtools library
> user system elapsed
> 0.970 0.001 0.971
On the other hand this takes way more time with large datasets
system.time(classIntervals(mydata$myvar, n=5, style="jenks"))
> Timing stopped at: 1081.894 1.345 1083.511