In the particular case you are having problems with is special: polygon = [(0,0), (0,10), (10,0), (10,10)]
Changing the order of points in a polygon can have significant impact on algorithms.
If you draw your polygon on a graph you'll see you have a horizontal hourglass shape. The polygon border overlaps itself. In geospatial analysis this overlap is not allowed because visually and logically you now have two closed polygons with a common intersection point. By the way most geospatial software doesn't deal well with triangles either.
In this case the point at 9,9 will trick the ray casting algorithm used in your method above because it can easily cross the doubled-over polygon boundary twice.
Please run the following code to see what is going on. (9,9) is on the line and this algorithm doesn't account for it. (5,8) is way outside:
import turtle as t
polygon = [(0,0), (0,100), (100,0), (100,100)]
t.goto(0,0)
fp = None
for p in polygon:
t.goto(p)
if not fp: fp=p
t.goto(fp)
t.up()
t.goto(90,90)
t.write("90,90")
t.dot(10)
t.goto(50,80)
t.write("50,80")
t.dot(10)
t.done()
This code handles the (9,9) edge case: http://geospatialpython.com/2011/08/point-in-polygon-2-on-line.html