0004-line-fixes-v03.patch

From 82fd7f782555f7390210e8c07f811bb1ca419659 Mon Sep 17 00:00:00 2001
From: Emre Hasegeli <emre@hasegeli.com>
Date: Sun, 28 May 2017 11:35:17 +0200
Subject: [PATCH 4/4] line-fixes-v03

Fix obvious problems around the line datatype

I have noticed some line operators retuning wrong results, and Tom Lane
spotted similar problems on more places.  Source history reveals that
during 1990s, the internal format of the line datatype is changed, but
most functions haven't got the hint.  The fixes include:

* Reject invalid specification A=B=0 on receive
* Avoid division by zero on perpendicular operator
* Fix intersection and distance operators when neither A nor B are 1
* Avoid point distance operator crashing on float precision loss
* Fix line segment distance by getting the minimum as the comment suggested

The changes are also aiming make line operators more symmetric.  Changing
results for same lines in different order are surely not expected.

Previous discussion:
https://www.postgresql.org/message-id/flat/CAE2gYzw_-z%3DV2kh8QqFjenu%3D8MJXzOP44wRW%3DAzzeamrmTT1%3DQ%40mail.gmail.com
---
 src/backend/utils/adt/geo_ops.c | 134 ++++++++++++++++++++++++----------------
 1 file changed, 81 insertions(+), 53 deletions(-)

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index dbf1de0ca2..c8714ff78f 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -1003,20 +1003,25 @@ line_recv(PG_FUNCTION_ARGS)
 {
 	StringInfo	buf = (StringInfo) PG_GETARG_POINTER(0);
 	LINE	   *line;
 
 	line = (LINE *) palloc(sizeof(LINE));
 
 	line->A = pq_getmsgfloat8(buf);
 	line->B = pq_getmsgfloat8(buf);
 	line->C = pq_getmsgfloat8(buf);
 
+	if (line->A == 0.0 && line->B == 0.0)
+		ereport(ERROR,
+				(errcode(ERRCODE_INVALID_TEXT_REPRESENTATION),
+				 errmsg("invalid line specification: A and B cannot both be zero")));
+
 	PG_RETURN_LINE_P(line);
 }
 
 /*
  *		line_send			- converts line to binary format
  */
 Datum
 line_send(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
@@ -1123,25 +1128,28 @@ line_parallel(PG_FUNCTION_ARGS)
 }
 
 Datum
 line_perp(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
 
 	if (FPzero(l1->A))
 		PG_RETURN_BOOL(FPzero(l2->B));
-	else if (FPzero(l1->B))
+	if (FPzero(l1->B))
 		PG_RETURN_BOOL(FPzero(l2->A));
+	if (FPzero(l2->A))
+		PG_RETURN_BOOL(FPzero(l1->B));
+	if (FPzero(l2->B))
+		PG_RETURN_BOOL(FPzero(l1->A));
 
-	PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->B),
-								   float8_mul(l1->B, l2->A)), -1.0));
+	PG_RETURN_BOOL(FPeq(float8_mul(l1->A, l2->B), -float8_mul(l1->B, l2->A)));
 }
 
 Datum
 line_vertical(PG_FUNCTION_ARGS)
 {
 	LINE	   *line = PG_GETARG_LINE_P(0);
 
 	PG_RETURN_BOOL(FPzero(line->B));
 }
 
@@ -1153,20 +1161,21 @@ line_horizontal(PG_FUNCTION_ARGS)
 	PG_RETURN_BOOL(FPzero(line->A));
 }
 
 Datum
 line_eq(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
 	float8		ratio;
 
+	/* A and B cannot be both 0, but we never really know with the EPSILON. */
 	if (!FPzero(l2->A))
 		ratio = float8_div(l1->A, l2->A);
 	else if (!FPzero(l2->B))
 		ratio = float8_div(l1->B, l2->B);
 	else if (!FPzero(l2->C))
 		ratio = float8_div(l1->C, l2->C);
 	else
 		ratio = 1.0;
 
 	PG_RETURN_BOOL(FPeq(l1->A, float8_mul(ratio, l2->A)) &&
@@ -1180,30 +1189,36 @@ line_eq(PG_FUNCTION_ARGS)
  *---------------------------------------------------------*/
 
 /* line_distance()
  * Distance between two lines.
  */
 Datum
 line_distance(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
-	float8		result;
-	Point		tmp;
+	float8		ratio;
 
 	if (line_interpt_internal(NULL, l1, l2))
 		PG_RETURN_FLOAT8(0.0);
-	if (FPzero(l1->B))			/* vertical? */
-		PG_RETURN_FLOAT8(fabs(float8_mi(l1->C, l2->C)));
-	point_construct(&tmp, 0.0, l1->C);
-	result = dist_pl_internal(&tmp, l2);
-	PG_RETURN_FLOAT8(result);
+
+	/* A and B cannot be both 0, but we never really know with the EPSILON. */
+	if (!FPzero(l2->A))
+		ratio = float8_div(l1->A, l2->A);
+	else if (!FPzero(l2->B))
+		ratio = float8_div(l1->B, l2->B);
+	else
+		ratio = 1.0;
+
+	PG_RETURN_FLOAT8(float8_div(fabs(float8_mi(l1->C,
+											   float8_mul(ratio, l2->C))),
+								float8_hypot(l1->A, l1->B)));
 }
 
 /* line_interpt()
  * Point where two lines l1, l2 intersect (if any)
  */
 Datum
 line_interpt(PG_FUNCTION_ARGS)
 {
 	LINE	   *l1 = PG_GETARG_LINE_P(0);
 	LINE	   *l2 = PG_GETARG_LINE_P(1);
@@ -1226,38 +1241,62 @@ line_interpt(PG_FUNCTION_ARGS)
  * NOTE: If the lines are identical then we will find they are parallel
  * and report "no intersection".  This is a little weird, but since
  * there's no *unique* intersection, maybe it's appropriate behavior.
  */
 static bool
 line_interpt_internal(Point *result, LINE *l1, LINE *l2)
 {
 	float8		x,
 				y;
 
-	if (FPzero(l1->B))			/* l1 vertical? */
+	if (FPzero(l1->A))			/* l1 horizontal? */
+	{
+		if (FPzero(l2->A))		/* l2 horizontal? */
+			return false;
+
+		y = float8_div(-l1->C, l1->B);
+		x = float8_div(-float8_pl(float8_mul(l2->B, y), l2->C), l2->A);
+	}
+	else if (FPzero(l1->B))		/* l1 vertical? */
 	{
 		if (FPzero(l2->B))		/* l2 vertical? */
 			return false;
 
-		x = l1->C;
-		y = float8_pl(float8_mul(l2->A, x), l2->C);
+		x = float8_div(-l1->C, l1->A);
+		y = float8_div(-float8_pl(float8_mul(l2->A, x), l2->C), l2->B);
+	}
+	else if (FPzero(l2->A))		/* l2 horizontal? */
+	{
+		if (FPzero(l1->A))		/* l1 horizontal? */
+			return false;
+
+		y = float8_div(-l2->C, l2->B);
+		x = float8_div(-float8_pl(float8_mul(l1->B, y), l1->C), l1->A);
+	}
+	else if (FPzero(l2->B))		/* l2 vertical? */
+	{
+		if (FPzero(l1->B))		/* l1 vertical? */
+			return false;
+
+		x = float8_div(-l2->C, l2->A);
+		y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B);
 	}
 	else
 	{
-		if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B))))
+		if (FPeq(float8_mul(l1->A, l2->B), float8_mul(l2->A, l1->B)))
 			return false;
 
-		if (FPzero(l2->B))
-			x = l2->C;
-		else
-			x = float8_div(float8_mi(l1->C, l2->C), float8_mi(l2->A, l1->A));
-		y = float8_pl(float8_mul(l1->A, x), l1->C);
+		x = float8_div(float8_mi(float8_mul(l1->B, l2->C),
+								 float8_mul(l2->B, l1->C)),
+					   float8_mi(float8_mul(l1->A, l2->B),
+								 float8_mul(l2->A, l1->B)));
+		y = float8_div(-float8_pl(float8_mul(l1->A, x), l1->C), l1->B);
 	}
 	if (result != NULL)
 		point_construct(result, x, y);
 
 #ifdef GEODEBUG
 	printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n",
 		   DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C);
 	printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y);
 #endif
 
@@ -2061,46 +2100,35 @@ lseg_parallel(PG_FUNCTION_ARGS)
 
 	m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y);
 	m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y);
 
 	PG_RETURN_BOOL(FPeq(m1, m2));
 }
 
 /* lseg_perp()
  * Determine if two line segments are perpendicular.
  *
- * This code did not get the correct answer for
- *	'((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg
- * So, modified it to check explicitly for slope of vertical line
- *	returned by float8_slope() and the results seem better.
- * - thomas 1998-01-31
+ * We are comparing the slope of the first line segment with the inverse slope
+ * of the second line segment.
  */
 Datum
 lseg_perp(PG_FUNCTION_ARGS)
 {
 	LSEG	   *l1 = PG_GETARG_LSEG_P(0);
 	LSEG	   *l2 = PG_GETARG_LSEG_P(1);
 	float8		m1,
 				m2;
 
 	m1 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y);
-	m2 = float8_slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y);
+	m2 = float8_slope(l2->p[0].y, l2->p[1].y, l2->p[1].x, l2->p[0].x);
 
-#ifdef GEODEBUG
-	printf("lseg_perp- slopes are %g and %g\n", m1, m2);
-#endif
-	if (FPzero(m1))
-		PG_RETURN_BOOL(FPeq(m2, DBL_MAX));
-	else if (FPzero(m2))
-		PG_RETURN_BOOL(FPeq(m1, DBL_MAX));
-
-	PG_RETURN_BOOL(FPeq(float8_div(m1, m2), -1.0));
+	PG_RETURN_BOOL(FPeq(m1, m2));
 }
 
 Datum
 lseg_vertical(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 
 	PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x));
 }
 
@@ -2467,22 +2495,21 @@ dist_pb(PG_FUNCTION_ARGS)
 Datum
 dist_sl(PG_FUNCTION_ARGS)
 {
 	LSEG	   *lseg = PG_GETARG_LSEG_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
 	float8		result;
 
 	if (lseg_interpt_line_internal(NULL, lseg, line))
 		result = 0.0;
 	else
-		/* XXX shouldn't we take the min not max? */
-		result = float8_max(dist_pl_internal(&lseg->p[0], line),
+		result = float8_min(dist_pl_internal(&lseg->p[0], line),
 							dist_pl_internal(&lseg->p[1], line));
 
 	PG_RETURN_FLOAT8(result);
 }
 
 /*
  * Distance from a lseg to a box
  */
 Datum
 dist_sb(PG_FUNCTION_ARGS)
@@ -2663,43 +2690,44 @@ lseg_interpt_line_internal(Point *result, LSEG *lseg, LINE *line)
 /* close_pl -
  *		The intersection point of a perpendicular of the line
  *		through the point.
  */
 Datum
 close_pl(PG_FUNCTION_ARGS)
 {
 	Point	   *pt = PG_GETARG_POINT_P(0);
 	LINE	   *line = PG_GETARG_LINE_P(1);
 	Point	   *result;
-	float8		invm;
-	LINE		tmp;
 
 	result = (Point *) palloc(sizeof(Point));
 
 	if (FPzero(line->B))		/* vertical? */
+		point_construct(result, float8_div(-line->C, line->A), pt->y);
+	else if (FPzero(line->A))	/* horizontal? */
+		point_construct(result, pt->x, float8_div(-line->C, line->B));
+	else
 	{
-		result->x = line->C;
-		result->y = pt->y;
-		PG_RETURN_POINT_P(result);
-	}
-	if (FPzero(line->A))		/* horizontal? */
-	{
-		result->x = pt->x;
-		result->y = line->C;
-		PG_RETURN_POINT_P(result);
-	}
-	/* drop a perpendicular and find the intersection point */
+		LINE		perp;
+
+		/*
+		 * Drop a perpendicular and find the intersection point
+		 *
+		 * We need to invert and flip the sign on the slope to get the
+		 * perpendicular.  We might lose some precision on the division. It
+		 * shouldn't be as much to turn the line.  If it happens anyway, we
+		 * will assume the point it on the line.
+		 */
+		line_construct_pm(&perp, pt, float8_div(line->B, line->A));
+		if (!line_interpt_internal(result, &perp, line))
+			*result = *pt;
+	}
 
-	/* invert and flip the sign on the slope to get a perpendicular */
-	invm = float8_div(line->B, line->A);
-	line_construct_pm(&tmp, pt, invm);
-	line_interpt_internal(result, &tmp, line);
 	PG_RETURN_POINT_P(result);
 }
 
 
 /* close_ps()
  * Closest point on line segment to specified point.
  * Take the closest endpoint if the point is left, right,
  *	above, or below the segment, otherwise find the intersection
  *	point of the segment and its perpendicular through the point.
  *
-- 
2.13.5 (Apple Git-94)