Another attempt at replacing the current HYPOT macro with a much better
behaved function. I've added comments addressing those sections of code
that tripped people up before. As well as explaining some of the design
decisions. Feedback appreciated.
Index: src/backend/utils/adt/geo_ops.c
===================================================================
RCS file: /projects/cvsroot/pgsql/src/backend/utils/adt/geo_ops.c,v
retrieving revision 1.103
diff -c -r1.103 geo_ops.c
*** src/backend/utils/adt/geo_ops.c 28 Jul 2009 09:47:59 -0000 1.103
--- src/backend/utils/adt/geo_ops.c 28 Aug 2009 08:11:00 -0000
***************
*** 32,37 ****
--- 32,38 ----
* Internal routines
*/
+ static double phypot(double x, double y);
static int point_inside(Point *p, int npts, Point *plist);
static int lseg_crossing(double x, double y, double px, double py);
static BOX *box_construct(double x1, double x2, double y1, double y2);
***************
*** 825,831 ****
box_cn(&a, box1);
box_cn(&b, box2);
! PG_RETURN_FLOAT8(HYPOT(a.x - b.x, a.y - b.y));
}
--- 826,832 ----
box_cn(&a, box1);
box_cn(&b, box2);
! PG_RETURN_FLOAT8(phypot(a.x-b.x, a.y-b.y));
}
***************
*** 1971,1977 ****
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
! PG_RETURN_FLOAT8(HYPOT(pt1->x - pt2->x, pt1->y - pt2->y));
}
double
--- 1972,1978 ----
Point *pt1 = PG_GETARG_POINT_P(0);
Point *pt2 = PG_GETARG_POINT_P(1);
! PG_RETURN_FLOAT8(phypot(pt1->x - pt2->x, pt1->y - pt2->y));
}
double
***************
*** 1979,1987 ****
{
#ifdef GEODEBUG
printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
! pt1->x, pt1->y, pt2->x, pt2->y, HYPOT(pt1->x - pt2->x, pt1->y -
pt2->y));
#endif
! return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y);
}
Datum
--- 1980,1988 ----
{
#ifdef GEODEBUG
printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n",
! pt1->x, pt1->y, pt2->x, pt2->y, phypot(pt1->x - pt2->x, pt1->y -
pt2->y));
#endif
! return phypot(pt1->x - pt2->x, pt1->y - pt2->y);
}
Datum
***************
*** 2444,2450 ****
dist_pl_internal(Point *pt, LINE *line)
{
return (line->A * pt->x + line->B * pt->y + line->C) /
! HYPOT(line->A, line->B);
}
Datum
--- 2445,2451 ----
dist_pl_internal(Point *pt, LINE *line)
{
return (line->A * pt->x + line->B * pt->y + line->C) /
! phypot(line->A, line->B);
}
Datum
***************
*** 4916,4922 ****
PointPGetDatum(point)));
result->center.x = p->x;
result->center.y = p->y;
! result->radius *= HYPOT(point->x, point->y);
PG_RETURN_CIRCLE_P(result);
}
--- 4917,4923 ----
PointPGetDatum(point)));
result->center.x = p->x;
result->center.y = p->y;
! result->radius *= phypot(point->x, point->y);
PG_RETURN_CIRCLE_P(result);
}
***************
*** 4936,4942 ****
PointPGetDatum(point)));
result->center.x = p->x;
result->center.y = p->y;
! result->radius /= HYPOT(point->x, point->y);
PG_RETURN_CIRCLE_P(result);
}
--- 4937,4943 ----
PointPGetDatum(point)));
result->center.x = p->x;
result->center.y = p->y;
! result->radius /= phypot(point->x, point->y);
PG_RETURN_CIRCLE_P(result);
}
***************
*** 5401,5403 ****
--- 5402,5468 ----
return FALSE;
}
+
+
+ /*-------------------------------------------------------------------------
+ * Determine the hypotenuse.
+ *
+ * If required, x and y are swapped to make x the larger number. The
+ * traditional formulae of x^2+y^2 is rearranged to factor x outside the
+ * sqrt. This allows computation of the hypotenuse for significantly
+ * larger values, and with a higher precision than otherwise normally
+ * possible.
+ *
+ * Only arguments of > 1.27e308 are at risk of causing overflow. Whereas
+ * the naive approach limits arguments to < 9.5e153.
+ *
+ * sqrt( x^2 + y^2 ) = sqrt( x^2( 1 + y^2/x^2) )
+ * = x * sqrt( 1 + y^2/x^2 )
+ * = x * sqrt( 1 + y/x * y/x )
+ *
+ * It is expected that this routine will eventually be replaced with the
+ * C99 hypot() function.
+ *
+ * This implementation conforms to IEEE Std 1003.1 and GLIBC, in that the
+ * case of hypot(inf,nan) results in INF, and not NAN. To obtain Single
+ * UNIX Specification behaviour, swap the order of the isnan() and isinf()
+ * test sections.
+ *
+ * The HYPOT macro this function is replacing did not check for, or
+ * signal on overflow. In addition no part of geo_ops checks for overflow,
+ * underflow or NaN's. This function maintains the same behaviour.
+ *-----------------------------------------------------------------------*/
+ double phypot( double x, double y )
+ {
+ double yx;
+
+ /* As per IEEE Std 1003.1 */
+ if( isinf(x) || isinf(y) )
+ return get_float8_infinity();
+
+ /* As per IEEE Std 1003.1 */
+ if( isnan(x) || isnan(y) )
+ return get_float8_nan();
+
+ /* If required, swap x and y */
+ x = fabs(x);
+ y = fabs(y);
+ if (x < y) {
+ double temp = x;
+ x = y;
+ y = temp;
+ }
+
+ /* As x is the larger value, this must be the correct answer. Also
+ * avoids division by zero. */
+ if( x == 0.0 )
+ return 0.0;
+
+ /* Trivial case. */
+ if( y == 0.0 )
+ return x;
+
+ /* Determine the hypotenuse */
+ yx = y/x;
+ return x*sqrt(1.0+yx*yx);
+ }
Index: src/include/utils/geo_decls.h
===================================================================
RCS file: /projects/cvsroot/pgsql/src/include/utils/geo_decls.h,v
retrieving revision 1.55
diff -c -r1.55 geo_decls.h
*** src/include/utils/geo_decls.h 1 Jan 2009 17:24:02 -0000 1.55
--- src/include/utils/geo_decls.h 28 Aug 2009 08:11:05 -0000
***************
*** 50,56 ****
#define FPge(A,B) ((A) >= (B))
#endif
- #define HYPOT(A, B) sqrt((A) * (A) + (B) * (B))
/*---------------------------------------------------------------------
* Point - (x,y)
--- 50,55 ----
--
Sent via pgsql-hackers mailing list ([email protected])
To make changes to your subscription:
http://www.postgresql.org/mailpref/pgsql-hackers
