Re: phypot - Pygmy Hippotause ?

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(at)postgresql(dot)org)
> To make changes to your subscription:
> http://www.postgresql.org/mailpref/pgsql-hackers

--
Bruce Momjian <bruce(at)momjian(dot)us> http://momjian.us
EnterpriseDB http://enterprisedb.com
PG East: http://www.enterprisedb.com/community/nav-pg-east-2010.do
+ If your life is a hard drive, Christ can be your backup. +