0004-line-fixes-v1.patch

From 79e38d38451e9b6c48183145d2396830e6442fda 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-v1

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 are:

* Make operators more symmetric
* 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

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 | 115 +++++++++++++++++++++++++++-------------
 1 file changed, 77 insertions(+), 38 deletions(-)

diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c
index b1a6d55f8d..681970b0a1 100644
--- a/src/backend/utils/adt/geo_ops.c
+++ b/src/backend/utils/adt/geo_ops.c
@@ -948,20 +948,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);
@@ -1095,25 +1100,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));
 }
 
@@ -1125,20 +1133,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)) &&
@@ -1152,30 +1161,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);
@@ -1198,38 +1213,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
 
@@ -2450,22 +2489,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)
@@ -2646,43 +2684,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.11.0 (Apple Git-81)