But using trigonometry, square roots and division is suboptimal. *Very* suboptimal for several reasons.
Hint: you can determine the solution of every geometric predicate with a determinant. Often this solution is not optimal either, but at least it is "only" a polynomial.
In this case, substract one point from the others and write the result in a 2x2 matrix. Then the sign of the determinant is your answer. For problems regarding rounding errors see here: http://www.cplusplus.com/forum/articles/3827/
There is an exact, lazy adaptive solution for this problem by Shewchuk: http://www.cs.berkeley.edu/~jrs/papers/robustr.ps
It describes floating point expansions and then uses them for the orientation tests in 2D and 3D.
If you use the formula I described above, you get (for the points p,q,r to avoid confusing them with the coordinates):
...which has an error bound of 22b-48
, where b is the number of bits you require to store ceil(max(px,py,qx,qy,rx,ry)).
So, basically, if you don't want Shewchuk's algorithm (it's a bit complicated), you can do:
int bitsize = calc_bitsize(ceil(max(px,py,qx,qy,rx,ry)));
double errorbound = exp(2, bitsize-48);
double result = (qx-px)*(ry-py)-(rx-px)*(qy-py);
if(result < errorbound)
if(result > errorbound)
return <calculate result using exact arithmetic, or, if you want a potentially wrong answer, return "p, q and r are on a line">
This, of course, implies that neither overflow nor underflow occurs. And if you know the largest number at compile time, you should perhaps precompute it (it is usually better to use, say, 10 for b instead of computing it every time, getting, say, 4-8 depending on the case (if you don't make a case study about degenerate cases))