I took a look at the R source code for the tukey test (HSD test) from the agricolae library http://cran.r-project.org/web/packages/agricolae/index.html.  The interesting thing is that they purposefully call round:

round(1-ptukey(abs(dif[k])*sqrt(2)/sdtdif,ntr,DFerror),6)

Also, they assume lm or aov.  So, we know that they don’t like big numbers. However, the nice thing is that the R file shows how to use the ptukey function!  Which is where I was a bit iffy.  So, I think I could get away with writing it in R and just instead of requiring lm or aov I just want the data then I can perform the requisite calculations.

I think I’ve convinced myself that I’ll trust mathematica and I’ll deal with how long it takes to get the data.  I don’t like how ineffiecient it is, but I can sacrifice time if it means that I get accurate results.  I don’t think that I’m qualified or have the time to create such a library.  It is interesting though that we don’t have many papers on knowing how many significicant digits are necessary for accurate tests.  Or that the R libraries don’t report the errors in their calculations.  This must be why people use Mathematica or Matlab because we trust them to do it right to enough sig figs that it doesn’t matter the error because it is insignificant.

So, I just need a nice desktop.