Thread: phypot - Pygmy Hippotause ?
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 ----
Paul Matthews <plm@netspace.net.au> wrote: > Feedback appreciated. + /* 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; The first test seems unnecessary if you have the second. x >= 0, so x can't be zero unless y is, too. Returning x on y == 0.0 will return 0.0 whenever x == 0.0. -Kevin
Kevin Grittner wrote:
>
> The first test seems unnecessary if you have the second.
> x >= 0, so x can't be zero unless y is, too.
> Returning x on y == 0.0 will return 0.0 whenever x == 0.0.
>
> -Kevin
>
Wish granted. :-)
--
--
Fools ignore complexity. Pragmatists suffer it.
Some can avoid it. Geniuses remove it.
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 29 Aug 2009 02:47:14 -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,5465 ----
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;
+ }
+
+ /* If x is also 0.0, then 0.0 is returned. This is the correct result
+ * and avoids division by zero. When x > 0.0 we can trivially avoid
+ * the calculation of the hypotenuse. */
+ 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 29 Aug 2009 02:47:19 -0000
***************
*** 50,56 ****
#define FPge(A,B) ((A) >= (B))
#endif
- #define HYPOT(A, B) sqrt((A) * (A) + (B) * (B))
/*---------------------------------------------------------------------
* Point - (x,y)
--- 50,55 ----
I assume this is not something we are supposed to apply.
---------------------------------------------------------------------------
Paul Matthews wrote:
> Kevin Grittner wrote:
> >
> > The first test seems unnecessary if you have the second.
> > x >= 0, so x can't be zero unless y is, too.
> > Returning x on y == 0.0 will return 0.0 whenever x == 0.0.
> >
> > -Kevin
> >
> Wish granted. :-)
>
> --
> --
> Fools ignore complexity. Pragmatists suffer it.
> Some can avoid it. Geniuses remove it.
>
> 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 29 Aug 2009 02:47:14 -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,5465 ----
>
> 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;
> + }
> +
> + /* If x is also 0.0, then 0.0 is returned. This is the correct result
> + * and avoids division by zero. When x > 0.0 we can trivially avoid
> + * the calculation of the hypotenuse. */
> + 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 29 Aug 2009 02:47:19 -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 (pgsql-hackers@postgresql.org)
> To make changes to your subscription:
> http://www.postgresql.org/mailpref/pgsql-hackers
-- Bruce Momjian <bruce@momjian.us> http://momjian.us EnterpriseDB
http://enterprisedb.comPG East: http://www.enterprisedb.com/community/nav-pg-east-2010.do + If your life is a hard
drive,Christ can be your backup. +
Bruce Momjian <bruce@momjian.us> wrote: > I assume this is not something we are supposed to apply. While it appears to improve conformance with the IEEE Std 1003.1 and expand the range of numbers which are correctly handled, it does more calculations. I wouldn't want to see it get in without performance testing and confirmation that existing apps don't rely on the non-standard behavior. I don't remember either of those happening. -Kevin
On Tue, Feb 23, 2010 at 4:03 PM, Kevin Grittner <Kevin.Grittner@wicourts.gov> wrote: > Bruce Momjian <bruce@momjian.us> wrote: >> I assume this is not something we are supposed to apply. > > While it appears to improve conformance with the IEEE Std 1003.1 and > expand the range of numbers which are correctly handled, it does > more calculations. I wouldn't want to see it get in without > performance testing and confirmation that existing apps don't rely > on the non-standard behavior. I don't remember either of those > happening. Perhaps we should add it to the next CommitFest? ...Robert
Robert Haas <robertmhaas@gmail.com> wrote: > Perhaps we should add it to the next CommitFest? Sounds like the right course of action to me. If nobody objects or beats me to it, I'll do that. -Kevin
"Kevin Grittner" wrote: > Robert Haas wrote:> >> Perhaps we should add it to the next CommitFest?> > Sounds like the right course of action to me. If nobody objects > or beats me to it, I'll do that. Done. -Kevin