see here: https://gist.github.com/skariel/da85943803a6f57a52fd
it implements fast and robust 2D and 3D orientation and in-circle tests according to the algorithms described in this paper: http://arxiv.org/abs/0901.4107 i.e. calculate using regular Floats while constraining the error. If in danger of getting a bad result repeat using exact calculation. this is good to implement fast Voronoy/Delaunay tessellations. See some usage in the recent biggest cosmological simulation Illustris: http://www.illustris-project.org/ These should be faster than say- http://www.cs.cmu.edu/~quake/robust.html which use adaptive precision I intend to implement Voronoy algorithms on top of this, some day :) Feedback welcome!
