## Rational Computational Geometry in Python

Article 2 of 2 in the series Robust Geometric Computation in Python

In the previous article, we looked at how a standard technique for determining the collinearity of points, based on computing the sign of the area of the triangle formed by two points on the line and a third query point. We discovered, that when used with Python's float type [1] the routine was unreliable in a region close to the line. This shortcoming has nothing to do with Python specifically and everything to do with the finite precision of the float number type. This time, we'll examine the behaviour of the algorithm more systematically using the following program:

def sign(x):
"""Determine the sign of x.

Returns:
-1 if x is negative, +1 if x is positive or 0 if x is zero.
"""
return (x > 0) - (x < 0)

def orientation(p, q, r):
"""Determine the orientation of three points in the plane.

Args:
p, q, r: Two-tuples representing coordinate pairs of three points.

Returns:
-1 if p, q, r is a turn to the right, +1 if p, q, r is a turn to the
left, otherwise 0 if p, q, and r are collinear.
"""
d = (q&#91;0&#93; - p&#91;0&#93;) * (r&#91;1&#93; - p&#91;1&#93;) - (q&#91;1&#93; - p&#91;1&#93;) * (r&#91;0&#93; - p&#91;0&#93;)
return sign(d)

def main():
"""
Test whether points immediately above and below the point (0.5, 0.5)
lie above, on, or below the line through (12.0, 12.0) and (24.0, 24.0).
"""
px = 0.5

pys = 0.49999999999999,
0.49999999999999006,
0.4999999999999901,
0.4999999999999902,
0.49999999999999023,
0.4999999999999903,
0.49999999999999034,
0.4999999999999904,
0.49999999999999045,
0.4999999999999905,
0.49999999999999056,
0.4999999999999906,
0.4999999999999907,
0.49999999999999073,
0.4999999999999908,
0.49999999999999084,
0.4999999999999909,
0.49999999999999095,
0.499999999999991,
0.49999999999999106,
0.4999999999999911,
0.4999999999999912,
0.49999999999999123,
0.4999999999999913,
0.49999999999999134,
0.4999999999999914,
0.49999999999999145,
0.4999999999999915,
0.49999999999999156,
0.4999999999999916,
0.4999999999999917,
0.49999999999999173,
0.4999999999999918,
0.49999999999999184,
0.4999999999999919,
0.49999999999999195,
0.499999999999992,
0.49999999999999206,
0.4999999999999921,
0.4999999999999922,
0.49999999999999223,
0.4999999999999923,
0.49999999999999234,
0.4999999999999924,
0.49999999999999245,
0.4999999999999925,
0.49999999999999256,
0.4999999999999926,
0.4999999999999927,
0.49999999999999273,
0.4999999999999928,
0.49999999999999284,
0.4999999999999929,
0.49999999999999295,
0.499999999999993,
0.49999999999999306,
0.4999999999999931,
0.49999999999999317,
0.4999999999999932,
0.4999999999999933,
0.49999999999999334,
0.4999999999999934,
0.49999999999999345,
0.4999999999999935,
0.49999999999999356,
0.4999999999999936,
0.49999999999999367,
0.4999999999999937,
0.4999999999999938,
0.49999999999999384,
0.4999999999999939,
0.49999999999999395,
0.499999999999994,
0.49999999999999406,
0.4999999999999941,
0.49999999999999417,
0.4999999999999942,
0.4999999999999943,
0.49999999999999434,
0.4999999999999944,
0.49999999999999445,
0.4999999999999945,
0.49999999999999456,
0.4999999999999946,
0.49999999999999467,
0.4999999999999947,
0.4999999999999948,
0.49999999999999484,
0.4999999999999949,
0.49999999999999495,
0.499999999999995,
0.49999999999999506,
0.4999999999999951,
0.49999999999999517,
0.4999999999999952,
0.4999999999999953,
0.49999999999999534,
0.4999999999999954,
0.49999999999999545,
0.4999999999999955,
0.49999999999999556,
0.4999999999999956,
0.49999999999999567,
0.4999999999999957,
0.4999999999999958,
0.49999999999999584,
0.4999999999999959,
0.49999999999999595,
0.499999999999996,
0.49999999999999606,
0.4999999999999961,
0.49999999999999617,
0.4999999999999962,
0.4999999999999963,
0.49999999999999634,
0.4999999999999964,
0.49999999999999645,
0.4999999999999965,
0.49999999999999656,
0.4999999999999966,
0.49999999999999667,
0.4999999999999967,
0.4999999999999968,
0.49999999999999684,
0.4999999999999969,
0.49999999999999695,
0.499999999999997,
0.49999999999999706,
0.4999999999999971,
0.49999999999999717,
0.4999999999999972,
0.4999999999999973,
0.49999999999999734,
0.4999999999999974,
0.49999999999999745,
0.4999999999999975,
0.49999999999999756,
0.4999999999999976,
0.49999999999999767,
0.4999999999999977,
0.4999999999999978,
0.49999999999999784,
0.4999999999999979,
0.49999999999999795,
0.499999999999998,
0.49999999999999806,
0.4999999999999981,
0.49999999999999817,
0.4999999999999982,
0.4999999999999983,
0.49999999999999833,
0.4999999999999984,
0.49999999999999845,
0.4999999999999985,
0.49999999999999856,
0.4999999999999986,
0.49999999999999867,
0.4999999999999987,
0.4999999999999988,
0.49999999999999883,
0.4999999999999989,
0.49999999999999895,
0.499999999999999,
0.49999999999999906,
0.4999999999999991,
0.49999999999999917,
0.4999999999999992,
0.4999999999999993,
0.49999999999999933,
0.4999999999999994,
0.49999999999999944,
0.4999999999999995,
0.49999999999999956,
0.4999999999999996,
0.49999999999999967,
0.4999999999999997,
0.4999999999999998,
0.49999999999999983,
0.4999999999999999,
0.49999999999999994,  # The previous representable float < 0.5
0.5,
0.5000000000000001,   # The next representable float > 0.5
0.5000000000000002,
0.5000000000000003,
0.5000000000000004,
0.5000000000000006,
0.5000000000000007,
0.5000000000000008,
0.5000000000000009,
0.500000000000001,
0.5000000000000011,
0.5000000000000012,
0.5000000000000013,
0.5000000000000014,
0.5000000000000016,
0.5000000000000017,
0.5000000000000018,
0.5000000000000019,
0.500000000000002,
0.5000000000000021,
0.5000000000000022,
0.5000000000000023,
0.5000000000000024,
0.5000000000000026,
0.5000000000000027,
0.5000000000000028,
0.5000000000000029,
0.500000000000003,
0.5000000000000031,
0.5000000000000032,
0.5000000000000033,
0.5000000000000034,
0.5000000000000036,
0.5000000000000037,
0.5000000000000038,
0.5000000000000039,
0.500000000000004,
0.5000000000000041,
0.5000000000000042,
0.5000000000000043,
0.5000000000000044,
0.5000000000000046,
0.5000000000000047,
0.5000000000000048,
0.5000000000000049,
0.500000000000005,
0.5000000000000051,
0.5000000000000052,
0.5000000000000053,
0.5000000000000054,
0.5000000000000056,
0.5000000000000057,
0.5000000000000058,
0.5000000000000059,
0.500000000000006,
0.5000000000000061,
0.5000000000000062,
0.5000000000000063,
0.5000000000000064,
0.5000000000000066,
0.5000000000000067,
0.5000000000000068,
0.5000000000000069,
0.500000000000007,
0.5000000000000071,
0.5000000000000072,
0.5000000000000073,
0.5000000000000074,
0.5000000000000075,
0.5000000000000077,
0.5000000000000078,
0.5000000000000079,
0.500000000000008,
0.5000000000000081,
0.5000000000000082,
0.5000000000000083,
0.5000000000000084,
0.5000000000000085,
0.5000000000000087,
0.5000000000000088,
0.5000000000000089,
0.500000000000009,
0.5000000000000091,
0.5000000000000092,
0.5000000000000093,
0.5000000000000094,
0.5000000000000095,
0.5000000000000097,
0.5000000000000098,
0.5000000000000099,
0.50000000000001]

q = (12.0, 12.0)
r = (24.0, 24.0)

for py in pys:
p = (px, py)
o = orientation(p, q, r)
print("orientation(({p[0]:>3}, {p[1]:<19}) q, r) -> {o:>2}".format(
p=p, o=o))

if __name__  == '__main__':
main()


The program includes definitions of our sign() and orientation() functions, together with a main() function which runs the test. The main function includes a list of the 271 nearest representable $$y$$-coordinate values to 0.5. We haven't included the code to generate these values successive float values because it's somewhat besides the point; we're referenced the necessary technique in the previous article.

The program iterates over these py values and performs the orientation test each time, printing the result. The complex format string is used to get readable output which lines up in columns. When we look at that output we see an intricate pattern of results emerge, which isn't even symmetrical around the central 0.5 value:

orientation((0.5, 0.50000000000001   ) q, r) ->  1
orientation((0.5, 0.5000000000000099 ) q, r) ->  1
orientation((0.5, 0.5000000000000098 ) q, r) ->  1
orientation((0.5, 0.5000000000000097 ) q, r) ->  1
orientation((0.5, 0.5000000000000095 ) q, r) ->  1
orientation((0.5, 0.5000000000000094 ) q, r) ->  1
orientation((0.5, 0.5000000000000093 ) q, r) ->  1
orientation((0.5, 0.5000000000000092 ) q, r) ->  1
orientation((0.5, 0.5000000000000091 ) q, r) ->  1
orientation((0.5, 0.500000000000009  ) q, r) ->  1
orientation((0.5, 0.5000000000000089 ) q, r) ->  1
orientation((0.5, 0.5000000000000088 ) q, r) ->  1
orientation((0.5, 0.5000000000000087 ) q, r) ->  1
orientation((0.5, 0.5000000000000085 ) q, r) ->  1
orientation((0.5, 0.5000000000000084 ) q, r) ->  1
orientation((0.5, 0.5000000000000083 ) q, r) ->  1
orientation((0.5, 0.5000000000000082 ) q, r) ->  1
orientation((0.5, 0.5000000000000081 ) q, r) ->  1
orientation((0.5, 0.500000000000008  ) q, r) ->  1
orientation((0.5, 0.5000000000000079 ) q, r) ->  1
orientation((0.5, 0.5000000000000078 ) q, r) ->  1
orientation((0.5, 0.5000000000000077 ) q, r) ->  1
orientation((0.5, 0.5000000000000075 ) q, r) ->  1
orientation((0.5, 0.5000000000000074 ) q, r) ->  1
orientation((0.5, 0.5000000000000073 ) q, r) ->  1
orientation((0.5, 0.5000000000000072 ) q, r) ->  1
orientation((0.5, 0.5000000000000071 ) q, r) ->  1
orientation((0.5, 0.500000000000007  ) q, r) ->  1
orientation((0.5, 0.5000000000000069 ) q, r) ->  1
orientation((0.5, 0.5000000000000068 ) q, r) ->  1
orientation((0.5, 0.5000000000000067 ) q, r) ->  1
orientation((0.5, 0.5000000000000066 ) q, r) ->  1
orientation((0.5, 0.5000000000000064 ) q, r) ->  1
orientation((0.5, 0.5000000000000063 ) q, r) ->  1
orientation((0.5, 0.5000000000000062 ) q, r) ->  1
orientation((0.5, 0.5000000000000061 ) q, r) ->  1
orientation((0.5, 0.500000000000006  ) q, r) ->  1
orientation((0.5, 0.5000000000000059 ) q, r) ->  1
orientation((0.5, 0.5000000000000058 ) q, r) ->  1
orientation((0.5, 0.5000000000000057 ) q, r) ->  1
orientation((0.5, 0.5000000000000056 ) q, r) ->  1
orientation((0.5, 0.5000000000000054 ) q, r) ->  1
orientation((0.5, 0.5000000000000053 ) q, r) ->  1
orientation((0.5, 0.5000000000000052 ) q, r) ->  1
orientation((0.5, 0.5000000000000051 ) q, r) ->  1
orientation((0.5, 0.500000000000005  ) q, r) ->  1
orientation((0.5, 0.5000000000000049 ) q, r) ->  1
orientation((0.5, 0.5000000000000048 ) q, r) ->  1
orientation((0.5, 0.5000000000000047 ) q, r) ->  1
orientation((0.5, 0.5000000000000046 ) q, r) ->  1
orientation((0.5, 0.5000000000000044 ) q, r) ->  0
orientation((0.5, 0.5000000000000043 ) q, r) ->  0
orientation((0.5, 0.5000000000000042 ) q, r) ->  0
orientation((0.5, 0.5000000000000041 ) q, r) ->  0
orientation((0.5, 0.500000000000004  ) q, r) ->  0
orientation((0.5, 0.5000000000000039 ) q, r) ->  0
orientation((0.5, 0.5000000000000038 ) q, r) ->  0
orientation((0.5, 0.5000000000000037 ) q, r) ->  0
orientation((0.5, 0.5000000000000036 ) q, r) ->  0
orientation((0.5, 0.5000000000000034 ) q, r) ->  0
orientation((0.5, 0.5000000000000033 ) q, r) ->  0
orientation((0.5, 0.5000000000000032 ) q, r) ->  0
orientation((0.5, 0.5000000000000031 ) q, r) ->  0
orientation((0.5, 0.500000000000003  ) q, r) ->  0
orientation((0.5, 0.5000000000000029 ) q, r) ->  0
orientation((0.5, 0.5000000000000028 ) q, r) ->  0
orientation((0.5, 0.5000000000000027 ) q, r) ->  0
orientation((0.5, 0.5000000000000026 ) q, r) ->  0
orientation((0.5, 0.5000000000000024 ) q, r) ->  0
orientation((0.5, 0.5000000000000023 ) q, r) ->  0
orientation((0.5, 0.5000000000000022 ) q, r) ->  0
orientation((0.5, 0.5000000000000021 ) q, r) ->  0
orientation((0.5, 0.500000000000002  ) q, r) ->  0
orientation((0.5, 0.5000000000000019 ) q, r) ->  0
orientation((0.5, 0.5000000000000018 ) q, r) ->  1
orientation((0.5, 0.5000000000000017 ) q, r) ->  1
orientation((0.5, 0.5000000000000016 ) q, r) ->  1
orientation((0.5, 0.5000000000000014 ) q, r) ->  1
orientation((0.5, 0.5000000000000013 ) q, r) ->  1
orientation((0.5, 0.5000000000000012 ) q, r) ->  1
orientation((0.5, 0.5000000000000011 ) q, r) ->  1
orientation((0.5, 0.500000000000001  ) q, r) ->  1
orientation((0.5, 0.5000000000000009 ) q, r) ->  0
orientation((0.5, 0.5000000000000008 ) q, r) ->  0
orientation((0.5, 0.5000000000000007 ) q, r) ->  0
orientation((0.5, 0.5000000000000006 ) q, r) ->  0
orientation((0.5, 0.5000000000000004 ) q, r) ->  0
orientation((0.5, 0.5000000000000003 ) q, r) ->  0
orientation((0.5, 0.5000000000000002 ) q, r) ->  0
orientation((0.5, 0.5000000000000001 ) q, r) ->  0
orientation((0.5, 0.5                ) q, r) ->  0
orientation((0.5, 0.49999999999999994) q, r) ->  0
orientation((0.5, 0.4999999999999999 ) q, r) ->  0
orientation((0.5, 0.49999999999999983) q, r) ->  0
orientation((0.5, 0.4999999999999998 ) q, r) ->  0
orientation((0.5, 0.4999999999999997 ) q, r) ->  0
orientation((0.5, 0.49999999999999967) q, r) ->  0
orientation((0.5, 0.4999999999999996 ) q, r) ->  0
orientation((0.5, 0.49999999999999956) q, r) ->  0
orientation((0.5, 0.4999999999999995 ) q, r) ->  0
orientation((0.5, 0.49999999999999944) q, r) ->  0
orientation((0.5, 0.4999999999999994 ) q, r) ->  0
orientation((0.5, 0.49999999999999933) q, r) ->  0
orientation((0.5, 0.4999999999999993 ) q, r) ->  0
orientation((0.5, 0.4999999999999992 ) q, r) ->  0
orientation((0.5, 0.49999999999999917) q, r) ->  0
orientation((0.5, 0.4999999999999991 ) q, r) ->  0
orientation((0.5, 0.49999999999999906) q, r) -> -1
orientation((0.5, 0.499999999999999  ) q, r) -> -1
orientation((0.5, 0.49999999999999895) q, r) -> -1
orientation((0.5, 0.4999999999999989 ) q, r) -> -1
orientation((0.5, 0.49999999999999883) q, r) -> -1
orientation((0.5, 0.4999999999999988 ) q, r) -> -1
orientation((0.5, 0.4999999999999987 ) q, r) -> -1
orientation((0.5, 0.49999999999999867) q, r) -> -1
orientation((0.5, 0.4999999999999986 ) q, r) -> -1
orientation((0.5, 0.49999999999999856) q, r) -> -1
orientation((0.5, 0.4999999999999985 ) q, r) -> -1
orientation((0.5, 0.49999999999999845) q, r) -> -1
orientation((0.5, 0.4999999999999984 ) q, r) -> -1
orientation((0.5, 0.49999999999999833) q, r) -> -1
orientation((0.5, 0.4999999999999983 ) q, r) -> -1
orientation((0.5, 0.4999999999999982 ) q, r) -> -1
orientation((0.5, 0.49999999999999817) q, r) ->  0
orientation((0.5, 0.4999999999999981 ) q, r) ->  0
orientation((0.5, 0.49999999999999806) q, r) ->  0
orientation((0.5, 0.499999999999998  ) q, r) ->  0
orientation((0.5, 0.49999999999999795) q, r) ->  0
orientation((0.5, 0.4999999999999979 ) q, r) ->  0
orientation((0.5, 0.49999999999999784) q, r) ->  0
orientation((0.5, 0.4999999999999978 ) q, r) ->  0
orientation((0.5, 0.4999999999999977 ) q, r) ->  0
orientation((0.5, 0.49999999999999767) q, r) ->  0
orientation((0.5, 0.4999999999999976 ) q, r) ->  0
orientation((0.5, 0.49999999999999756) q, r) ->  0
orientation((0.5, 0.4999999999999975 ) q, r) ->  0
orientation((0.5, 0.49999999999999745) q, r) ->  0
orientation((0.5, 0.4999999999999974 ) q, r) ->  0
orientation((0.5, 0.49999999999999734) q, r) ->  0
orientation((0.5, 0.4999999999999973 ) q, r) ->  0
orientation((0.5, 0.4999999999999972 ) q, r) ->  0
orientation((0.5, 0.49999999999999717) q, r) ->  0
orientation((0.5, 0.4999999999999971 ) q, r) ->  0
orientation((0.5, 0.49999999999999706) q, r) ->  0
orientation((0.5, 0.499999999999997  ) q, r) ->  0
orientation((0.5, 0.49999999999999695) q, r) ->  0
orientation((0.5, 0.4999999999999969 ) q, r) ->  0
orientation((0.5, 0.49999999999999684) q, r) ->  0
orientation((0.5, 0.4999999999999968 ) q, r) ->  0
orientation((0.5, 0.4999999999999967 ) q, r) ->  0
orientation((0.5, 0.49999999999999667) q, r) ->  0
orientation((0.5, 0.4999999999999966 ) q, r) ->  0
orientation((0.5, 0.49999999999999656) q, r) ->  0
orientation((0.5, 0.4999999999999965 ) q, r) ->  0
orientation((0.5, 0.49999999999999645) q, r) ->  0
orientation((0.5, 0.4999999999999964 ) q, r) ->  0
orientation((0.5, 0.49999999999999634) q, r) ->  0
orientation((0.5, 0.4999999999999963 ) q, r) ->  0
orientation((0.5, 0.4999999999999962 ) q, r) ->  0
orientation((0.5, 0.49999999999999617) q, r) ->  0
orientation((0.5, 0.4999999999999961 ) q, r) ->  0
orientation((0.5, 0.49999999999999606) q, r) ->  0
orientation((0.5, 0.499999999999996  ) q, r) ->  0
orientation((0.5, 0.49999999999999595) q, r) ->  0
orientation((0.5, 0.4999999999999959 ) q, r) ->  0
orientation((0.5, 0.49999999999999584) q, r) ->  0
orientation((0.5, 0.4999999999999958 ) q, r) ->  0
orientation((0.5, 0.4999999999999957 ) q, r) ->  0
orientation((0.5, 0.49999999999999567) q, r) ->  0
orientation((0.5, 0.4999999999999956 ) q, r) ->  0
orientation((0.5, 0.49999999999999556) q, r) ->  0
orientation((0.5, 0.4999999999999955 ) q, r) -> -1
orientation((0.5, 0.49999999999999545) q, r) -> -1
orientation((0.5, 0.4999999999999954 ) q, r) -> -1
orientation((0.5, 0.49999999999999534) q, r) -> -1
orientation((0.5, 0.4999999999999953 ) q, r) -> -1
orientation((0.5, 0.4999999999999952 ) q, r) -> -1
orientation((0.5, 0.49999999999999517) q, r) -> -1
orientation((0.5, 0.4999999999999951 ) q, r) -> -1
orientation((0.5, 0.49999999999999506) q, r) -> -1
orientation((0.5, 0.499999999999995  ) q, r) -> -1
orientation((0.5, 0.49999999999999495) q, r) -> -1
orientation((0.5, 0.4999999999999949 ) q, r) -> -1
orientation((0.5, 0.49999999999999484) q, r) -> -1
orientation((0.5, 0.4999999999999948 ) q, r) -> -1
orientation((0.5, 0.4999999999999947 ) q, r) -> -1
orientation((0.5, 0.49999999999999467) q, r) -> -1
orientation((0.5, 0.4999999999999946 ) q, r) -> -1
orientation((0.5, 0.49999999999999456) q, r) -> -1
orientation((0.5, 0.4999999999999945 ) q, r) -> -1
orientation((0.5, 0.49999999999999445) q, r) -> -1
orientation((0.5, 0.4999999999999944 ) q, r) -> -1
orientation((0.5, 0.49999999999999434) q, r) -> -1
orientation((0.5, 0.4999999999999943 ) q, r) -> -1
orientation((0.5, 0.4999999999999942 ) q, r) -> -1
orientation((0.5, 0.49999999999999417) q, r) -> -1
orientation((0.5, 0.4999999999999941 ) q, r) -> -1
orientation((0.5, 0.49999999999999406) q, r) -> -1
orientation((0.5, 0.499999999999994  ) q, r) -> -1
orientation((0.5, 0.49999999999999395) q, r) -> -1
orientation((0.5, 0.4999999999999939 ) q, r) -> -1
orientation((0.5, 0.49999999999999384) q, r) -> -1
orientation((0.5, 0.4999999999999938 ) q, r) -> -1
orientation((0.5, 0.4999999999999937 ) q, r) -> -1
orientation((0.5, 0.49999999999999367) q, r) -> -1
orientation((0.5, 0.4999999999999936 ) q, r) -> -1
orientation((0.5, 0.49999999999999356) q, r) -> -1
orientation((0.5, 0.4999999999999935 ) q, r) -> -1
orientation((0.5, 0.49999999999999345) q, r) -> -1
orientation((0.5, 0.4999999999999934 ) q, r) -> -1
orientation((0.5, 0.49999999999999334) q, r) -> -1
orientation((0.5, 0.4999999999999933 ) q, r) -> -1
orientation((0.5, 0.4999999999999932 ) q, r) -> -1
orientation((0.5, 0.49999999999999317) q, r) -> -1
orientation((0.5, 0.4999999999999931 ) q, r) -> -1
orientation((0.5, 0.49999999999999306) q, r) -> -1
orientation((0.5, 0.499999999999993  ) q, r) -> -1
orientation((0.5, 0.49999999999999295) q, r) -> -1
orientation((0.5, 0.4999999999999929 ) q, r) -> -1
orientation((0.5, 0.49999999999999284) q, r) -> -1
orientation((0.5, 0.4999999999999928 ) q, r) -> -1
orientation((0.5, 0.49999999999999273) q, r) -> -1
orientation((0.5, 0.4999999999999927 ) q, r) -> -1
orientation((0.5, 0.4999999999999926 ) q, r) -> -1
orientation((0.5, 0.49999999999999256) q, r) -> -1
orientation((0.5, 0.4999999999999925 ) q, r) -> -1
orientation((0.5, 0.49999999999999245) q, r) -> -1
orientation((0.5, 0.4999999999999924 ) q, r) -> -1
orientation((0.5, 0.49999999999999234) q, r) -> -1
orientation((0.5, 0.4999999999999923 ) q, r) -> -1
orientation((0.5, 0.49999999999999223) q, r) -> -1
orientation((0.5, 0.4999999999999922 ) q, r) -> -1
orientation((0.5, 0.4999999999999921 ) q, r) -> -1
orientation((0.5, 0.49999999999999206) q, r) -> -1
orientation((0.5, 0.499999999999992  ) q, r) -> -1
orientation((0.5, 0.49999999999999195) q, r) -> -1
orientation((0.5, 0.4999999999999919 ) q, r) -> -1
orientation((0.5, 0.49999999999999184) q, r) -> -1
orientation((0.5, 0.4999999999999918 ) q, r) -> -1
orientation((0.5, 0.49999999999999173) q, r) -> -1
orientation((0.5, 0.4999999999999917 ) q, r) -> -1
orientation((0.5, 0.4999999999999916 ) q, r) -> -1
orientation((0.5, 0.49999999999999156) q, r) -> -1
orientation((0.5, 0.4999999999999915 ) q, r) -> -1
orientation((0.5, 0.49999999999999145) q, r) -> -1
orientation((0.5, 0.4999999999999914 ) q, r) -> -1
orientation((0.5, 0.49999999999999134) q, r) -> -1
orientation((0.5, 0.4999999999999913 ) q, r) -> -1
orientation((0.5, 0.49999999999999123) q, r) -> -1
orientation((0.5, 0.4999999999999912 ) q, r) -> -1
orientation((0.5, 0.4999999999999911 ) q, r) -> -1
orientation((0.5, 0.49999999999999106) q, r) -> -1
orientation((0.5, 0.499999999999991  ) q, r) -> -1
orientation((0.5, 0.49999999999999095) q, r) -> -1
orientation((0.5, 0.4999999999999909 ) q, r) -> -1
orientation((0.5, 0.49999999999999084) q, r) -> -1
orientation((0.5, 0.4999999999999908 ) q, r) -> -1
orientation((0.5, 0.49999999999999073) q, r) -> -1
orientation((0.5, 0.4999999999999907 ) q, r) -> -1
orientation((0.5, 0.4999999999999906 ) q, r) -> -1
orientation((0.5, 0.49999999999999056) q, r) -> -1
orientation((0.5, 0.4999999999999905 ) q, r) -> -1
orientation((0.5, 0.49999999999999045) q, r) -> -1
orientation((0.5, 0.4999999999999904 ) q, r) -> -1
orientation((0.5, 0.49999999999999034) q, r) -> -1
orientation((0.5, 0.4999999999999903 ) q, r) -> -1
orientation((0.5, 0.49999999999999023) q, r) -> -1
orientation((0.5, 0.4999999999999902 ) q, r) -> -1
orientation((0.5, 0.4999999999999901 ) q, r) -> -1
orientation((0.5, 0.49999999999999006) q, r) -> -1
orientation((0.5, 0.49999999999999   ) q, r) -> -1



The colour coding (added later) represents whether the algorithm reckons the points are above the line (in blue), on the line (in yellow) or below the line (in red). The only point which is actually on the line is in green.

By this point you should at least be wary of using floating point arithmetic for geometric computation. Lest you think this can easily be solved by introducing a tolerance value, or some other clunky solution, we'll save you the bother by pointing out that doing do merely moves these fringing effects to the edge of the tolerance zone.

What to do? Fortunately, as we alluded to at the beginning of this tale, Python gives us a solution into the form of the rational numbers, implemented as the Fraction type.

Let's make a small change to our program, converting all numbers to Fractions before proceeding with the computation. We'll do this by modifying the orientation() to convert each of its three arguments from a tuple containing a pair of numeric objects into a pair of Fractions. The Fraction constructor accepts a selection of numeric types, including float:

def orientation(p, q, r):
"""Determine the orientation of three points in the plane.

Args:
p, q, r: Two-tuples representing coordinate pairs of three points.

Returns:
-1 if p, q, r is a turn to the right, +1 if p, q, r is a turn to the
left, otherwise 0 if p, q, and r are collinear.
"""
p = (Fraction(p[0]), Fraction(p[1]))
q = (Fraction(q[0]), Fraction(q[1]))
r = (Fraction(r[0]), Fraction(r[1]))

d = (q[0] - p[0]) * (r[1] - p[1]) - (q[1] - p[1]) * (r[0] - p[0])
return sign(d)


The variable d will now also be a Fraction and the sign() function will work as expected with this type since it only uses comparison to zero.

Let's run our modified example:

orientation((0.5, 0.49999999999999   ) q, r) -> -1
orientation((0.5, 0.49999999999999006) q, r) -> -1
orientation((0.5, 0.4999999999999901 ) q, r) -> -1
orientation((0.5, 0.4999999999999902 ) q, r) -> -1
orientation((0.5, 0.49999999999999023) q, r) -> -1
orientation((0.5, 0.4999999999999903 ) q, r) -> -1
orientation((0.5, 0.49999999999999034) q, r) -> -1
orientation((0.5, 0.4999999999999904 ) q, r) -> -1
orientation((0.5, 0.49999999999999045) q, r) -> -1
orientation((0.5, 0.4999999999999905 ) q, r) -> -1
orientation((0.5, 0.49999999999999056) q, r) -> -1
orientation((0.5, 0.4999999999999906 ) q, r) -> -1
orientation((0.5, 0.4999999999999907 ) q, r) -> -1
orientation((0.5, 0.49999999999999073) q, r) -> -1
orientation((0.5, 0.4999999999999908 ) q, r) -> -1
orientation((0.5, 0.49999999999999084) q, r) -> -1
orientation((0.5, 0.4999999999999909 ) q, r) -> -1
orientation((0.5, 0.49999999999999095) q, r) -> -1
orientation((0.5, 0.499999999999991  ) q, r) -> -1
orientation((0.5, 0.49999999999999106) q, r) -> -1
orientation((0.5, 0.4999999999999911 ) q, r) -> -1
orientation((0.5, 0.4999999999999912 ) q, r) -> -1
orientation((0.5, 0.49999999999999123) q, r) -> -1
orientation((0.5, 0.4999999999999913 ) q, r) -> -1
orientation((0.5, 0.49999999999999134) q, r) -> -1
orientation((0.5, 0.4999999999999914 ) q, r) -> -1
orientation((0.5, 0.49999999999999145) q, r) -> -1
orientation((0.5, 0.4999999999999915 ) q, r) -> -1
orientation((0.5, 0.49999999999999156) q, r) -> -1
orientation((0.5, 0.4999999999999916 ) q, r) -> -1
orientation((0.5, 0.4999999999999917 ) q, r) -> -1
orientation((0.5, 0.49999999999999173) q, r) -> -1
orientation((0.5, 0.4999999999999918 ) q, r) -> -1
orientation((0.5, 0.49999999999999184) q, r) -> -1
orientation((0.5, 0.4999999999999919 ) q, r) -> -1
orientation((0.5, 0.49999999999999195) q, r) -> -1
orientation((0.5, 0.499999999999992  ) q, r) -> -1
orientation((0.5, 0.49999999999999206) q, r) -> -1
orientation((0.5, 0.4999999999999921 ) q, r) -> -1
orientation((0.5, 0.4999999999999922 ) q, r) -> -1
orientation((0.5, 0.49999999999999223) q, r) -> -1
orientation((0.5, 0.4999999999999923 ) q, r) -> -1
orientation((0.5, 0.49999999999999234) q, r) -> -1
orientation((0.5, 0.4999999999999924 ) q, r) -> -1
orientation((0.5, 0.49999999999999245) q, r) -> -1
orientation((0.5, 0.4999999999999925 ) q, r) -> -1
orientation((0.5, 0.49999999999999256) q, r) -> -1
orientation((0.5, 0.4999999999999926 ) q, r) -> -1
orientation((0.5, 0.4999999999999927 ) q, r) -> -1
orientation((0.5, 0.49999999999999273) q, r) -> -1
orientation((0.5, 0.4999999999999928 ) q, r) -> -1
orientation((0.5, 0.49999999999999284) q, r) -> -1
orientation((0.5, 0.4999999999999929 ) q, r) -> -1
orientation((0.5, 0.49999999999999295) q, r) -> -1
orientation((0.5, 0.499999999999993  ) q, r) -> -1
orientation((0.5, 0.49999999999999306) q, r) -> -1
orientation((0.5, 0.4999999999999931 ) q, r) -> -1
orientation((0.5, 0.49999999999999317) q, r) -> -1
orientation((0.5, 0.4999999999999932 ) q, r) -> -1
orientation((0.5, 0.4999999999999933 ) q, r) -> -1
orientation((0.5, 0.49999999999999334) q, r) -> -1
orientation((0.5, 0.4999999999999934 ) q, r) -> -1
orientation((0.5, 0.49999999999999345) q, r) -> -1
orientation((0.5, 0.4999999999999935 ) q, r) -> -1
orientation((0.5, 0.49999999999999356) q, r) -> -1
orientation((0.5, 0.4999999999999936 ) q, r) -> -1
orientation((0.5, 0.49999999999999367) q, r) -> -1
orientation((0.5, 0.4999999999999937 ) q, r) -> -1
orientation((0.5, 0.4999999999999938 ) q, r) -> -1
orientation((0.5, 0.49999999999999384) q, r) -> -1
orientation((0.5, 0.4999999999999939 ) q, r) -> -1
orientation((0.5, 0.49999999999999395) q, r) -> -1
orientation((0.5, 0.499999999999994  ) q, r) -> -1
orientation((0.5, 0.49999999999999406) q, r) -> -1
orientation((0.5, 0.4999999999999941 ) q, r) -> -1
orientation((0.5, 0.49999999999999417) q, r) -> -1
orientation((0.5, 0.4999999999999942 ) q, r) -> -1
orientation((0.5, 0.4999999999999943 ) q, r) -> -1
orientation((0.5, 0.49999999999999434) q, r) -> -1
orientation((0.5, 0.4999999999999944 ) q, r) -> -1
orientation((0.5, 0.49999999999999445) q, r) -> -1
orientation((0.5, 0.4999999999999945 ) q, r) -> -1
orientation((0.5, 0.49999999999999456) q, r) -> -1
orientation((0.5, 0.4999999999999946 ) q, r) -> -1
orientation((0.5, 0.49999999999999467) q, r) -> -1
orientation((0.5, 0.4999999999999947 ) q, r) -> -1
orientation((0.5, 0.4999999999999948 ) q, r) -> -1
orientation((0.5, 0.49999999999999484) q, r) -> -1
orientation((0.5, 0.4999999999999949 ) q, r) -> -1
orientation((0.5, 0.49999999999999495) q, r) -> -1
orientation((0.5, 0.499999999999995  ) q, r) -> -1
orientation((0.5, 0.49999999999999506) q, r) -> -1
orientation((0.5, 0.4999999999999951 ) q, r) -> -1
orientation((0.5, 0.49999999999999517) q, r) -> -1
orientation((0.5, 0.4999999999999952 ) q, r) -> -1
orientation((0.5, 0.4999999999999953 ) q, r) -> -1
orientation((0.5, 0.49999999999999534) q, r) -> -1
orientation((0.5, 0.4999999999999954 ) q, r) -> -1
orientation((0.5, 0.49999999999999545) q, r) -> -1
orientation((0.5, 0.4999999999999955 ) q, r) -> -1
orientation((0.5, 0.49999999999999556) q, r) -> -1
orientation((0.5, 0.4999999999999956 ) q, r) -> -1
orientation((0.5, 0.49999999999999567) q, r) -> -1
orientation((0.5, 0.4999999999999957 ) q, r) -> -1
orientation((0.5, 0.4999999999999958 ) q, r) -> -1
orientation((0.5, 0.49999999999999584) q, r) -> -1
orientation((0.5, 0.4999999999999959 ) q, r) -> -1
orientation((0.5, 0.49999999999999595) q, r) -> -1
orientation((0.5, 0.499999999999996  ) q, r) -> -1
orientation((0.5, 0.49999999999999606) q, r) -> -1
orientation((0.5, 0.4999999999999961 ) q, r) -> -1
orientation((0.5, 0.49999999999999617) q, r) -> -1
orientation((0.5, 0.4999999999999962 ) q, r) -> -1
orientation((0.5, 0.4999999999999963 ) q, r) -> -1
orientation((0.5, 0.49999999999999634) q, r) -> -1
orientation((0.5, 0.4999999999999964 ) q, r) -> -1
orientation((0.5, 0.49999999999999645) q, r) -> -1
orientation((0.5, 0.4999999999999965 ) q, r) -> -1
orientation((0.5, 0.49999999999999656) q, r) -> -1
orientation((0.5, 0.4999999999999966 ) q, r) -> -1
orientation((0.5, 0.49999999999999667) q, r) -> -1
orientation((0.5, 0.4999999999999967 ) q, r) -> -1
orientation((0.5, 0.4999999999999968 ) q, r) -> -1
orientation((0.5, 0.49999999999999684) q, r) -> -1
orientation((0.5, 0.4999999999999969 ) q, r) -> -1
orientation((0.5, 0.49999999999999695) q, r) -> -1
orientation((0.5, 0.499999999999997  ) q, r) -> -1
orientation((0.5, 0.49999999999999706) q, r) -> -1
orientation((0.5, 0.4999999999999971 ) q, r) -> -1
orientation((0.5, 0.49999999999999717) q, r) -> -1
orientation((0.5, 0.4999999999999972 ) q, r) -> -1
orientation((0.5, 0.4999999999999973 ) q, r) -> -1
orientation((0.5, 0.49999999999999734) q, r) -> -1
orientation((0.5, 0.4999999999999974 ) q, r) -> -1
orientation((0.5, 0.49999999999999745) q, r) -> -1
orientation((0.5, 0.4999999999999975 ) q, r) -> -1
orientation((0.5, 0.49999999999999756) q, r) -> -1
orientation((0.5, 0.4999999999999976 ) q, r) -> -1
orientation((0.5, 0.49999999999999767) q, r) -> -1
orientation((0.5, 0.4999999999999977 ) q, r) -> -1
orientation((0.5, 0.4999999999999978 ) q, r) -> -1
orientation((0.5, 0.49999999999999784) q, r) -> -1
orientation((0.5, 0.4999999999999979 ) q, r) -> -1
orientation((0.5, 0.49999999999999795) q, r) -> -1
orientation((0.5, 0.499999999999998  ) q, r) -> -1
orientation((0.5, 0.49999999999999806) q, r) -> -1
orientation((0.5, 0.4999999999999981 ) q, r) -> -1
orientation((0.5, 0.49999999999999817) q, r) -> -1
orientation((0.5, 0.4999999999999982 ) q, r) -> -1
orientation((0.5, 0.4999999999999983 ) q, r) -> -1
orientation((0.5, 0.49999999999999833) q, r) -> -1
orientation((0.5, 0.4999999999999984 ) q, r) -> -1
orientation((0.5, 0.49999999999999845) q, r) -> -1
orientation((0.5, 0.4999999999999985 ) q, r) -> -1
orientation((0.5, 0.49999999999999856) q, r) -> -1
orientation((0.5, 0.4999999999999986 ) q, r) -> -1
orientation((0.5, 0.49999999999999867) q, r) -> -1
orientation((0.5, 0.4999999999999987 ) q, r) -> -1
orientation((0.5, 0.4999999999999988 ) q, r) -> -1
orientation((0.5, 0.49999999999999883) q, r) -> -1
orientation((0.5, 0.4999999999999989 ) q, r) -> -1
orientation((0.5, 0.49999999999999895) q, r) -> -1
orientation((0.5, 0.499999999999999  ) q, r) -> -1
orientation((0.5, 0.49999999999999906) q, r) -> -1
orientation((0.5, 0.4999999999999991 ) q, r) -> -1
orientation((0.5, 0.49999999999999917) q, r) -> -1
orientation((0.5, 0.4999999999999992 ) q, r) -> -1
orientation((0.5, 0.4999999999999993 ) q, r) -> -1
orientation((0.5, 0.49999999999999933) q, r) -> -1
orientation((0.5, 0.4999999999999994 ) q, r) -> -1
orientation((0.5, 0.49999999999999944) q, r) -> -1
orientation((0.5, 0.4999999999999995 ) q, r) -> -1
orientation((0.5, 0.49999999999999956) q, r) -> -1
orientation((0.5, 0.4999999999999996 ) q, r) -> -1
orientation((0.5, 0.49999999999999967) q, r) -> -1
orientation((0.5, 0.4999999999999997 ) q, r) -> -1
orientation((0.5, 0.4999999999999998 ) q, r) -> -1
orientation((0.5, 0.49999999999999983) q, r) -> -1
orientation((0.5, 0.4999999999999999 ) q, r) -> -1
orientation((0.5, 0.49999999999999994) q, r) -> -1
orientation((0.5, 0.5                ) q, r) ->  0
orientation((0.5, 0.5000000000000001 ) q, r) ->  1
orientation((0.5, 0.5000000000000002 ) q, r) ->  1
orientation((0.5, 0.5000000000000003 ) q, r) ->  1
orientation((0.5, 0.5000000000000004 ) q, r) ->  1
orientation((0.5, 0.5000000000000006 ) q, r) ->  1
orientation((0.5, 0.5000000000000007 ) q, r) ->  1
orientation((0.5, 0.5000000000000008 ) q, r) ->  1
orientation((0.5, 0.5000000000000009 ) q, r) ->  1
orientation((0.5, 0.500000000000001  ) q, r) ->  1
orientation((0.5, 0.5000000000000011 ) q, r) ->  1
orientation((0.5, 0.5000000000000012 ) q, r) ->  1
orientation((0.5, 0.5000000000000013 ) q, r) ->  1
orientation((0.5, 0.5000000000000014 ) q, r) ->  1
orientation((0.5, 0.5000000000000016 ) q, r) ->  1
orientation((0.5, 0.5000000000000017 ) q, r) ->  1
orientation((0.5, 0.5000000000000018 ) q, r) ->  1
orientation((0.5, 0.5000000000000019 ) q, r) ->  1
orientation((0.5, 0.500000000000002  ) q, r) ->  1
orientation((0.5, 0.5000000000000021 ) q, r) ->  1
orientation((0.5, 0.5000000000000022 ) q, r) ->  1
orientation((0.5, 0.5000000000000023 ) q, r) ->  1
orientation((0.5, 0.5000000000000024 ) q, r) ->  1
orientation((0.5, 0.5000000000000026 ) q, r) ->  1
orientation((0.5, 0.5000000000000027 ) q, r) ->  1
orientation((0.5, 0.5000000000000028 ) q, r) ->  1
orientation((0.5, 0.5000000000000029 ) q, r) ->  1
orientation((0.5, 0.500000000000003  ) q, r) ->  1
orientation((0.5, 0.5000000000000031 ) q, r) ->  1
orientation((0.5, 0.5000000000000032 ) q, r) ->  1
orientation((0.5, 0.5000000000000033 ) q, r) ->  1
orientation((0.5, 0.5000000000000034 ) q, r) ->  1
orientation((0.5, 0.5000000000000036 ) q, r) ->  1
orientation((0.5, 0.5000000000000037 ) q, r) ->  1
orientation((0.5, 0.5000000000000038 ) q, r) ->  1
orientation((0.5, 0.5000000000000039 ) q, r) ->  1
orientation((0.5, 0.500000000000004  ) q, r) ->  1
orientation((0.5, 0.5000000000000041 ) q, r) ->  1
orientation((0.5, 0.5000000000000042 ) q, r) ->  1
orientation((0.5, 0.5000000000000043 ) q, r) ->  1
orientation((0.5, 0.5000000000000044 ) q, r) ->  1
orientation((0.5, 0.5000000000000046 ) q, r) ->  1
orientation((0.5, 0.5000000000000047 ) q, r) ->  1
orientation((0.5, 0.5000000000000048 ) q, r) ->  1
orientation((0.5, 0.5000000000000049 ) q, r) ->  1
orientation((0.5, 0.500000000000005  ) q, r) ->  1
orientation((0.5, 0.5000000000000051 ) q, r) ->  1
orientation((0.5, 0.5000000000000052 ) q, r) ->  1
orientation((0.5, 0.5000000000000053 ) q, r) ->  1
orientation((0.5, 0.5000000000000054 ) q, r) ->  1
orientation((0.5, 0.5000000000000056 ) q, r) ->  1
orientation((0.5, 0.5000000000000057 ) q, r) ->  1
orientation((0.5, 0.5000000000000058 ) q, r) ->  1
orientation((0.5, 0.5000000000000059 ) q, r) ->  1
orientation((0.5, 0.500000000000006  ) q, r) ->  1
orientation((0.5, 0.5000000000000061 ) q, r) ->  1
orientation((0.5, 0.5000000000000062 ) q, r) ->  1
orientation((0.5, 0.5000000000000063 ) q, r) ->  1
orientation((0.5, 0.5000000000000064 ) q, r) ->  1
orientation((0.5, 0.5000000000000066 ) q, r) ->  1
orientation((0.5, 0.5000000000000067 ) q, r) ->  1
orientation((0.5, 0.5000000000000068 ) q, r) ->  1
orientation((0.5, 0.5000000000000069 ) q, r) ->  1
orientation((0.5, 0.500000000000007  ) q, r) ->  1
orientation((0.5, 0.5000000000000071 ) q, r) ->  1
orientation((0.5, 0.5000000000000072 ) q, r) ->  1
orientation((0.5, 0.5000000000000073 ) q, r) ->  1
orientation((0.5, 0.5000000000000074 ) q, r) ->  1
orientation((0.5, 0.5000000000000075 ) q, r) ->  1
orientation((0.5, 0.5000000000000077 ) q, r) ->  1
orientation((0.5, 0.5000000000000078 ) q, r) ->  1
orientation((0.5, 0.5000000000000079 ) q, r) ->  1
orientation((0.5, 0.500000000000008  ) q, r) ->  1
orientation((0.5, 0.5000000000000081 ) q, r) ->  1
orientation((0.5, 0.5000000000000082 ) q, r) ->  1
orientation((0.5, 0.5000000000000083 ) q, r) ->  1
orientation((0.5, 0.5000000000000084 ) q, r) ->  1
orientation((0.5, 0.5000000000000085 ) q, r) ->  1
orientation((0.5, 0.5000000000000087 ) q, r) ->  1
orientation((0.5, 0.5000000000000088 ) q, r) ->  1
orientation((0.5, 0.5000000000000089 ) q, r) ->  1
orientation((0.5, 0.500000000000009  ) q, r) ->  1
orientation((0.5, 0.5000000000000091 ) q, r) ->  1
orientation((0.5, 0.5000000000000092 ) q, r) ->  1
orientation((0.5, 0.5000000000000093 ) q, r) ->  1
orientation((0.5, 0.5000000000000094 ) q, r) ->  1
orientation((0.5, 0.5000000000000095 ) q, r) ->  1
orientation((0.5, 0.5000000000000097 ) q, r) ->  1
orientation((0.5, 0.5000000000000098 ) q, r) ->  1
orientation((0.5, 0.5000000000000099 ) q, r) ->  1
orientation((0.5, 0.50000000000001   ) q, r) ->  1



Using Fractions internally, our orientation() function gets the full benefit of exact arithmetic with effectively infinite precision and consequently produces an exact result with only one position of p being reported as collinear with q and r.

In the next article, we'll more fully explore the behaviour of the non-robust float-based version of this function based graphically, to get an impression of how lines are 'seen' by floating-point geometric functions.

 [1] Python's float is an IEEE-754 double precision 64-bit float.
Category: misc
Tags: Numerics, Python