diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c index d6ce5ccf6cc..c3383bd6da6 100644 --- a/src/backend/access/gist/gistproc.c +++ b/src/backend/access/gist/gistproc.c @@ -55,7 +55,7 @@ rt_box_union(BOX *n, const BOX *a, const BOX *b) * Size of a BOX for penalty-calculation purposes. * The result can be +Infinity, but not NaN. */ -static double +static float8 size_box(const BOX *box) { /* @@ -76,20 +76,21 @@ size_box(const BOX *box) */ if (isnan(box->high.x) || isnan(box->high.y)) return get_float8_infinity(); - return (box->high.x - box->low.x) * (box->high.y - box->low.y); + return float8_mul(float8_mi(box->high.x, box->low.x), + float8_mi(box->high.y, box->low.y)); } /* * Return amount by which the union of the two boxes is larger than * the original BOX's area. The result can be +Infinity, but not NaN. */ -static double +static float8 box_penalty(const BOX *original, const BOX *new) { BOX unionbox; rt_box_union(&unionbox, original, new); - return size_box(&unionbox) - size_box(original); + return float8_mi(size_box(&unionbox), size_box(original)); } /* @@ -263,7 +264,7 @@ typedef struct /* Index of entry in the initial array */ int index; /* Delta between penalties of entry insertion into different groups */ - double delta; + float8 delta; } CommonEntry; /* @@ -279,13 +280,13 @@ typedef struct bool first; /* true if no split was selected yet */ - double leftUpper; /* upper bound of left interval */ - double rightLower; /* lower bound of right interval */ + float8 leftUpper; /* upper bound of left interval */ + float8 rightLower; /* lower bound of right interval */ float4 ratio; float4 overlap; int dim; /* axis of this split */ - double range; /* width of general MBR projection to the + float8 range; /* width of general MBR projection to the * selected axis */ } ConsiderSplitContext; @@ -294,7 +295,7 @@ typedef struct */ typedef struct { - double lower, + float8 lower, upper; } SplitInterval; @@ -304,7 +305,7 @@ typedef struct static int interval_cmp_lower(const void *i1, const void *i2) { - double lower1 = ((const SplitInterval *) i1)->lower, + float8 lower1 = ((const SplitInterval *) i1)->lower, lower2 = ((const SplitInterval *) i2)->lower; return float8_cmp_internal(lower1, lower2); @@ -316,7 +317,7 @@ interval_cmp_lower(const void *i1, const void *i2) static int interval_cmp_upper(const void *i1, const void *i2) { - double upper1 = ((const SplitInterval *) i1)->upper, + float8 upper1 = ((const SplitInterval *) i1)->upper, upper2 = ((const SplitInterval *) i2)->upper; return float8_cmp_internal(upper1, upper2); @@ -339,14 +340,14 @@ non_negative(float val) */ static inline void g_box_consider_split(ConsiderSplitContext *context, int dimNum, - double rightLower, int minLeftCount, - double leftUpper, int maxLeftCount) + float8 rightLower, int minLeftCount, + float8 leftUpper, int maxLeftCount) { int leftCount, rightCount; float4 ratio, overlap; - double range; + float8 range; /* * Calculate entries distribution ratio assuming most uniform distribution @@ -369,8 +370,7 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum, * Ratio of split - quotient between size of lesser group and total * entries count. */ - ratio = ((float4) Min(leftCount, rightCount)) / - ((float4) context->entriesCount); + ratio = float4_div(Min(leftCount, rightCount), context->entriesCount); if (ratio > LIMIT_RATIO) { @@ -384,11 +384,13 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum, * or less range with same overlap. */ if (dimNum == 0) - range = context->boundingBox.high.x - context->boundingBox.low.x; + range = float8_mi(context->boundingBox.high.x, + context->boundingBox.low.x); else - range = context->boundingBox.high.y - context->boundingBox.low.y; + range = float8_mi(context->boundingBox.high.y, + context->boundingBox.low.y); - overlap = (leftUpper - rightLower) / range; + overlap = float8_div(float8_mi(leftUpper, rightLower), range); /* If there is no previous selection, select this */ if (context->first) @@ -444,20 +446,14 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum, /* * Compare common entries by their deltas. - * (We assume the deltas can't be NaN.) */ static int common_entry_cmp(const void *i1, const void *i2) { - double delta1 = ((const CommonEntry *) i1)->delta, + float8 delta1 = ((const CommonEntry *) i1)->delta, delta2 = ((const CommonEntry *) i2)->delta; - if (delta1 < delta2) - return -1; - else if (delta1 > delta2) - return 1; - else - return 0; + return float8_cmp_internal(delta1, delta2); } /* @@ -531,7 +527,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS) context.first = true; /* nothing selected yet */ for (dim = 0; dim < 2; dim++) { - double leftUpper, + float8 leftUpper, rightLower; int i1, i2; @@ -728,7 +724,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS) */ for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i)) { - double lower, + float8 lower, upper; /* @@ -783,7 +779,7 @@ gist_box_picksplit(PG_FUNCTION_ARGS) * Calculate minimum number of entries that must be placed in both * groups, to reach LIMIT_RATIO. */ - int m = ceil(LIMIT_RATIO * (double) nentries); + int m = ceil(LIMIT_RATIO * nentries); /* * Calculate delta between penalties of join "common entries" to @@ -792,8 +788,8 @@ gist_box_picksplit(PG_FUNCTION_ARGS) for (i = 0; i < commonEntriesCount; i++) { box = DatumGetBoxP(entryvec->vector[commonEntries[i].index].key); - commonEntries[i].delta = Abs(box_penalty(leftBox, box) - - box_penalty(rightBox, box)); + commonEntries[i].delta = Abs(float8_mi(box_penalty(leftBox, box), + box_penalty(rightBox, box))); } /* @@ -1107,10 +1103,10 @@ gist_circle_compress(PG_FUNCTION_ARGS) BOX *r; r = (BOX *) palloc(sizeof(BOX)); - r->high.x = in->center.x + in->radius; - r->low.x = in->center.x - in->radius; - r->high.y = in->center.y + in->radius; - r->low.y = in->center.y - in->radius; + r->high.x = float8_pl(in->center.x, in->radius); + r->low.x = float8_mi(in->center.x, in->radius); + r->high.y = float8_pl(in->center.y, in->radius); + r->low.y = float8_mi(in->center.y, in->radius); retval = (GISTENTRY *) palloc(sizeof(GISTENTRY)); gistentryinit(*retval, PointerGetDatum(r), @@ -1148,10 +1144,10 @@ gist_circle_consistent(PG_FUNCTION_ARGS) * rtree_internal_consistent even at leaf nodes. (This works in part * because the index entries are bounding boxes not circles.) */ - bbox.high.x = query->center.x + query->radius; - bbox.low.x = query->center.x - query->radius; - bbox.high.y = query->center.y + query->radius; - bbox.low.y = query->center.y - query->radius; + bbox.high.x = float8_pl(query->center.x, query->radius); + bbox.low.x = float8_mi(query->center.x, query->radius); + bbox.high.y = float8_pl(query->center.y, query->radius); + bbox.low.y = float8_mi(query->center.y, query->radius); result = rtree_internal_consistent(DatumGetBoxP(entry->key), &bbox, strategy); @@ -1216,10 +1212,10 @@ gist_point_fetch(PG_FUNCTION_ARGS) DatumGetFloat8(DirectFunctionCall2(point_distance, \ PointPGetDatum(p1), PointPGetDatum(p2))) -static double +static float8 computeDistance(bool isLeaf, BOX *box, Point *point) { - double result = 0.0; + float8 result = 0.0; if (isLeaf) { @@ -1237,9 +1233,9 @@ computeDistance(bool isLeaf, BOX *box, Point *point) /* point is over or below box */ Assert(box->low.y <= box->high.y); if (point->y > box->high.y) - result = point->y - box->high.y; + result = float8_mi(point->y, box->high.y); else if (point->y < box->low.y) - result = box->low.y - point->y; + result = float8_mi(box->low.y, point->y); else elog(ERROR, "inconsistent point values"); } @@ -1248,9 +1244,9 @@ computeDistance(bool isLeaf, BOX *box, Point *point) /* point is to left or right of box */ Assert(box->low.x <= box->high.x); if (point->x > box->high.x) - result = point->x - box->high.x; + result = float8_mi(point->x, box->high.x); else if (point->x < box->low.x) - result = box->low.x - point->x; + result = float8_mi(box->low.x, point->x); else elog(ERROR, "inconsistent point values"); } @@ -1258,7 +1254,7 @@ computeDistance(bool isLeaf, BOX *box, Point *point) { /* closest point will be a vertex */ Point p; - double subresult; + float8 subresult; result = point_point_distance(point, &box->low); @@ -1449,7 +1445,7 @@ gist_point_distance(PG_FUNCTION_ARGS) { GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2); - double distance; + float8 distance; StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset; switch (strategyGroup) @@ -1478,11 +1474,11 @@ gist_point_distance(PG_FUNCTION_ARGS) * This is a lower bound estimate of distance from point to indexed geometric * type. */ -static double +static float8 gist_bbox_distance(GISTENTRY *entry, Datum query, StrategyNumber strategy, bool *recheck) { - double distance; + float8 distance; StrategyNumber strategyGroup = strategy / GeoStrategyNumberOffset; /* Bounding box distance is always inexact. */ @@ -1512,7 +1508,7 @@ gist_circle_distance(PG_FUNCTION_ARGS) /* Oid subtype = PG_GETARG_OID(3); */ bool *recheck = (bool *) PG_GETARG_POINTER(4); - double distance; + float8 distance; distance = gist_bbox_distance(entry, query, strategy, recheck); @@ -1528,7 +1524,7 @@ gist_poly_distance(PG_FUNCTION_ARGS) /* Oid subtype = PG_GETARG_OID(3); */ bool *recheck = (bool *) PG_GETARG_POINTER(4); - double distance; + float8 distance; distance = gist_bbox_distance(entry, query, strategy, recheck); diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index bd8e21775db..13dc1ab6e60 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -59,7 +59,7 @@ static inline float8 lseg_sl(LSEG *lseg); static inline float8 lseg_invsl(LSEG *lseg); static bool lseg_interpt_line(Point *result, LSEG *lseg, LINE *line); static bool lseg_interpt_lseg(Point *result, LSEG *l1, LSEG *l2); -static int lseg_crossing(double x, double y, double px, double py); +static int lseg_crossing(float8 x, float8 y, float8 px, float8 py); static bool lseg_contain_point(LSEG *lseg, Point *point); static float8 lseg_closept_point(Point *result, LSEG *lseg, Point *pt); static float8 lseg_closept_line(Point *result, LSEG *lseg, LINE *line); @@ -69,9 +69,9 @@ static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2); static inline void box_construct(BOX *result, Point *pt1, Point *pt2); static void box_cn(Point *center, BOX *box); static bool box_ov(BOX *box1, BOX *box2); -static double box_ar(BOX *box); -static double box_ht(BOX *box); -static double box_wd(BOX *box); +static float8 box_ar(BOX *box); +static float8 box_ht(BOX *box); +static float8 box_wd(BOX *box); static bool box_contain_point(BOX *box, Point *point); static bool box_contain_box(BOX *box1, BOX *box2); static bool box_contain_lseg(BOX *box, LSEG *lseg); @@ -80,7 +80,7 @@ static float8 box_closept_point(Point *result, BOX *box, Point *point); static float8 box_closept_lseg(Point *result, BOX *box, LSEG *lseg); /* Routines for circles */ -static double circle_ar(CIRCLE *circle); +static float8 circle_ar(CIRCLE *circle); /* Routines for polygons */ static void make_bound_box(POLYGON *poly); @@ -91,10 +91,10 @@ static bool plist_same(int npts, Point *p1, Point *p2); static float8 dist_ppoly_internal(Point *pt, POLYGON *poly); /* Routines for encoding and decoding */ -static double single_decode(char *num, char **endptr_p, +static float8 single_decode(char *num, char **endptr_p, const char *type_name, const char *orig_string); static void single_encode(float8 x, StringInfo str); -static void pair_decode(char *str, double *x, double *y, char **endptr_p, +static void pair_decode(char *str, float8 *x, float8 *y, char **endptr_p, const char *type_name, const char *orig_string); static void pair_encode(float8 x, float8 y, StringInfo str); static int pair_count(char *s, char delim); @@ -146,7 +146,7 @@ static char *path_encode(enum path_delim path_delim, int npts, Point *pt); * and restore that order for text output - tgl 97/01/16 */ -static double +static float8 single_decode(char *num, char **endptr_p, const char *type_name, const char *orig_string) { @@ -163,7 +163,7 @@ single_encode(float8 x, StringInfo str) } /* single_encode() */ static void -pair_decode(char *str, double *x, double *y, char **endptr_p, +pair_decode(char *str, float8 *x, float8 *y, char **endptr_p, const char *type_name, const char *orig_string) { bool has_delim; @@ -376,19 +376,19 @@ box_in(PG_FUNCTION_ARGS) char *str = PG_GETARG_CSTRING(0); BOX *box = (BOX *) palloc(sizeof(BOX)); bool isopen; - double x, + float8 x, y; path_decode(str, false, 2, &(box->high), &isopen, NULL, "box", str); /* reorder corners if necessary... */ - if (box->high.x < box->low.x) + if (float8_lt(box->high.x, box->low.x)) { x = box->high.x; box->high.x = box->low.x; box->low.x = x; } - if (box->high.y < box->low.y) + if (float8_lt(box->high.y, box->low.y)) { y = box->high.y; box->high.y = box->low.y; @@ -416,7 +416,7 @@ box_recv(PG_FUNCTION_ARGS) { StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); BOX *box; - double x, + float8 x, y; box = (BOX *) palloc(sizeof(BOX)); @@ -427,13 +427,13 @@ box_recv(PG_FUNCTION_ARGS) box->low.y = pq_getmsgfloat8(buf); /* reorder corners if necessary... */ - if (box->high.x < box->low.x) + if (float8_lt(box->high.x, box->low.x)) { x = box->high.x; box->high.x = box->low.x; box->low.x = x; } - if (box->high.y < box->low.y) + if (float8_lt(box->high.y, box->low.y)) { y = box->high.y; box->high.y = box->low.y; @@ -466,7 +466,7 @@ box_send(PG_FUNCTION_ARGS) static inline void box_construct(BOX *result, Point *pt1, Point *pt2) { - if (pt1->x > pt2->x) + if (float8_gt(pt1->x, pt2->x)) { result->high.x = pt1->x; result->low.x = pt2->x; @@ -476,7 +476,7 @@ box_construct(BOX *result, Point *pt1, Point *pt2) result->high.x = pt2->x; result->low.x = pt1->x; } - if (pt1->y > pt2->y) + if (float8_gt(pt1->y, pt2->y)) { result->high.y = pt1->y; result->low.y = pt2->y; @@ -808,10 +808,10 @@ box_center(PG_FUNCTION_ARGS) /* box_ar - returns the area of the box. */ -static double +static float8 box_ar(BOX *box) { - return box_wd(box) * box_ht(box); + return float8_mul(box_wd(box), box_ht(box)); } @@ -820,28 +820,28 @@ box_ar(BOX *box) static void box_cn(Point *center, BOX *box) { - center->x = (box->high.x + box->low.x) / 2.0; - center->y = (box->high.y + box->low.y) / 2.0; + center->x = float8_div(float8_pl(box->high.x, box->low.x), 2.0); + center->y = float8_div(float8_pl(box->high.y, box->low.y), 2.0); } /* box_wd - returns the width (length) of the box * (horizontal magnitude). */ -static double +static float8 box_wd(BOX *box) { - return box->high.x - box->low.x; + return float8_mi(box->high.x, box->low.x); } /* box_ht - returns the height of the box * (vertical magnitude). */ -static double +static float8 box_ht(BOX *box) { - return box->high.y - box->low.y; + return float8_mi(box->high.y, box->low.y); } @@ -865,10 +865,10 @@ box_intersect(PG_FUNCTION_ARGS) result = (BOX *) palloc(sizeof(BOX)); - result->high.x = Min(box1->high.x, box2->high.x); - result->low.x = Max(box1->low.x, box2->low.x); - result->high.y = Min(box1->high.y, box2->high.y); - result->low.y = Max(box1->low.y, box2->low.y); + result->high.x = float8_min(box1->high.x, box2->high.x); + result->low.x = float8_max(box1->low.x, box2->low.x); + result->high.y = float8_min(box1->high.y, box2->high.y); + result->low.y = float8_max(box1->low.y, box2->low.y); PG_RETURN_BOX_P(result); } @@ -1014,8 +1014,8 @@ line_construct(LINE *result, Point *pt, float8 m) if (m == DBL_MAX) { /* vertical - use "x = C" */ - result->A = -1; - result->B = 0; + result->A = -1.0; + result->B = 0.0; result->C = pt->x; } else @@ -1023,7 +1023,7 @@ line_construct(LINE *result, Point *pt, float8 m) /* use "mx - y + yinter = 0" */ result->A = m; result->B = -1.0; - result->C = pt->y - m * pt->x; + result->C = float8_mi(pt->y, float8_mul(m, pt->x)); /* on some platforms, the preceding expression tends to produce -0 */ if (result->C == 0.0) result->C = 0.0; @@ -1079,7 +1079,8 @@ line_perp(PG_FUNCTION_ARGS) else if (FPzero(l1->B)) PG_RETURN_BOOL(FPzero(l2->A)); - PG_RETURN_BOOL(FPeq(((l1->A * l2->B) / (l1->B * l2->A)), -1.0)); + PG_RETURN_BOOL(FPeq(float8_div(float8_mul(l1->A, l2->B), + float8_mul(l1->B, l2->A)), -1.0)); } Datum @@ -1098,25 +1099,35 @@ line_horizontal(PG_FUNCTION_ARGS) PG_RETURN_BOOL(FPzero(line->A)); } + +/* + * Check whether the two lines are the same + * + * We consider NaNs values to be equal to each other to let those lines + * to be found. + */ Datum line_eq(PG_FUNCTION_ARGS) { LINE *l1 = PG_GETARG_LINE_P(0); LINE *l2 = PG_GETARG_LINE_P(1); - double k; + float8 ratio; - if (!FPzero(l2->A)) - k = l1->A / l2->A; - else if (!FPzero(l2->B)) - k = l1->B / l2->B; - else if (!FPzero(l2->C)) - k = l1->C / l2->C; + if (!FPzero(l2->A) && !isnan(l2->A)) + ratio = float8_div(l1->A, l2->A); + else if (!FPzero(l2->B) && !isnan(l2->B)) + ratio = float8_div(l1->B, l2->B); + else if (!FPzero(l2->C) && !isnan(l2->C)) + ratio = float8_div(l1->C, l2->C); else - k = 1.0; + ratio = 1.0; - PG_RETURN_BOOL(FPeq(l1->A, k * l2->A) && - FPeq(l1->B, k * l2->B) && - FPeq(l1->C, k * l2->C)); + PG_RETURN_BOOL((FPeq(l1->A, float8_mul(ratio, l2->A)) && + FPeq(l1->B, float8_mul(ratio, l2->B)) && + FPeq(l1->C, float8_mul(ratio, l2->C))) || + (float8_eq(l1->A, l2->A) && + float8_eq(l1->B, l2->B) && + float8_eq(l1->C, l2->C))); } @@ -1134,7 +1145,7 @@ line_invsl(LINE *line) return DBL_MAX; if (FPzero(line->B)) return 0.0; - return line->B / line->A; + return float8_div(line->B, line->A); } @@ -1152,7 +1163,7 @@ line_distance(PG_FUNCTION_ARGS) if (line_interpt_line(NULL, l1, l2)) /* intersecting? */ PG_RETURN_FLOAT8(0.0); if (FPzero(l1->B)) /* vertical? */ - PG_RETURN_FLOAT8(fabs(l1->C - l2->C)); + PG_RETURN_FLOAT8(fabs(float8_mi(l1->C, l2->C))); point_construct(&tmp, 0.0, l1->C); result = line_closept_point(NULL, l2, &tmp); PG_RETURN_FLOAT8(result); @@ -1185,11 +1196,15 @@ 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. + * + * If the lines have NaN constants, we will return true, and the intersection + * point would have NaN coordinates. We shouldn't return false in this case + * because that would mean the lines are parallel. */ static bool line_interpt_line(Point *result, LINE *l1, LINE *l2) { - double x, + float8 x, y; if (FPzero(l1->B)) /* l1 vertical? */ @@ -1198,20 +1213,20 @@ line_interpt_line(Point *result, LINE *l1, LINE *l2) return false; x = l1->C; - y = (l2->A * x + l2->C); + y = float8_pl(float8_mul(l2->A, x), l2->C); } else if (FPzero(l2->B)) /* l2 vertical? */ { x = l2->C; - y = (l1->A * x + l1->C); + y = float8_pl(float8_mul(l1->A, x), l1->C); } else { - if (FPeq(l2->A, l1->A * (l2->B / l1->B))) + if (FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B)))) return false; - x = (l1->C - l2->C) / (l2->A - l1->A); - y = (l1->A * x + l1->C); + 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); } if (result != NULL) @@ -1247,7 +1262,7 @@ Datum path_area(PG_FUNCTION_ARGS) { PATH *path = PG_GETARG_PATH_P(0); - double area = 0.0; + float8 area = 0.0; int i, j; @@ -1257,12 +1272,11 @@ path_area(PG_FUNCTION_ARGS) for (i = 0; i < path->npts; i++) { j = (i + 1) % path->npts; - area += path->p[i].x * path->p[j].y; - area -= path->p[i].y * path->p[j].x; + area = float8_pl(area, float8_mul(path->p[i].x, path->p[j].y)); + area = float8_mi(area, float8_mul(path->p[i].y, path->p[j].x)); } - area *= 0.5; - PG_RETURN_FLOAT8(area < 0.0 ? -area : area); + PG_RETURN_FLOAT8(float8_div(fabs(area), 2.0)); } @@ -1532,19 +1546,19 @@ path_inter(PG_FUNCTION_ARGS) b1.high.y = b1.low.y = p1->p[0].y; for (i = 1; i < p1->npts; i++) { - b1.high.x = Max(p1->p[i].x, b1.high.x); - b1.high.y = Max(p1->p[i].y, b1.high.y); - b1.low.x = Min(p1->p[i].x, b1.low.x); - b1.low.y = Min(p1->p[i].y, b1.low.y); + b1.high.x = float8_max(p1->p[i].x, b1.high.x); + b1.high.y = float8_max(p1->p[i].y, b1.high.y); + b1.low.x = float8_min(p1->p[i].x, b1.low.x); + b1.low.y = float8_min(p1->p[i].y, b1.low.y); } b2.high.x = b2.low.x = p2->p[0].x; b2.high.y = b2.low.y = p2->p[0].y; for (i = 1; i < p2->npts; i++) { - b2.high.x = Max(p2->p[i].x, b2.high.x); - b2.high.y = Max(p2->p[i].y, b2.high.y); - b2.low.x = Min(p2->p[i].x, b2.low.x); - b2.low.y = Min(p2->p[i].y, b2.low.y); + b2.high.x = float8_max(p2->p[i].x, b2.high.x); + b2.high.y = float8_max(p2->p[i].y, b2.high.y); + b2.low.x = float8_min(p2->p[i].x, b2.low.x); + b2.low.y = float8_min(p2->p[i].y, b2.low.y); } if (!box_ov(&b1, &b2)) PG_RETURN_BOOL(false); @@ -1634,7 +1648,7 @@ path_distance(PG_FUNCTION_ARGS) statlseg_construct(&seg2, &p2->p[jprev], &p2->p[j]); tmp = lseg_closept_lseg(NULL, &seg1, &seg2); - if (!have_min || tmp < min) + if (!have_min || float8_lt(tmp, min)) { min = tmp; have_min = true; @@ -1673,7 +1687,7 @@ path_length(PG_FUNCTION_ARGS) iprev = path->npts - 1; /* include the closure segment */ } - result += point_dt(&path->p[iprev], &path->p[i]); + result = float8_pl(result, point_dt(&path->p[iprev], &path->p[i])); } PG_RETURN_FLOAT8(result); @@ -1833,10 +1847,18 @@ point_ne(PG_FUNCTION_ARGS) PG_RETURN_BOOL(!point_eq_point(pt1, pt2)); } + +/* + * Check whether the two points are the same + * + * We consider NaNs coordinates to be equal to each other to let those points + * to be found. + */ static inline bool point_eq_point(Point *pt1, Point *pt2) { - return FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y); + return ((FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y)) || + (float8_eq(pt1->x, pt2->x) && float8_eq(pt1->y, pt2->y))); } @@ -1856,7 +1878,7 @@ point_distance(PG_FUNCTION_ARGS) static inline float8 point_dt(Point *pt1, Point *pt2) { - return HYPOT(pt1->x - pt2->x, pt1->y - pt2->y); + return HYPOT(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y)); } Datum @@ -1881,7 +1903,7 @@ point_sl(Point *pt1, Point *pt2) return DBL_MAX; if (FPeq(pt1->y, pt2->y)) return 0.0; - return (pt1->y - pt2->y) / (pt1->x - pt2->x); + return float8_div(float8_mi(pt1->y, pt2->y), float8_mi(pt1->x, pt2->x)); } @@ -1897,7 +1919,7 @@ point_invsl(Point *pt1, Point *pt2) return 0.0; if (FPeq(pt1->y, pt2->y)) return DBL_MAX; - return (pt1->x - pt2->x) / (pt2->y - pt1->y); + return float8_div(float8_mi(pt1->x, pt2->x), float8_mi(pt2->y, pt1->y)); } @@ -2171,8 +2193,8 @@ lseg_center(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); - result->x = (lseg->p[0].x + lseg->p[1].x) / 2.0; - result->y = (lseg->p[0].y + lseg->p[1].y) / 2.0; + result->x = float8_div(float8_pl(lseg->p[0].x, lseg->p[1].x), 2.0); + result->y = float8_div(float8_pl(lseg->p[0].y, lseg->p[1].y), 2.0); PG_RETURN_POINT_P(result); } @@ -2295,7 +2317,7 @@ dist_ppath(PG_FUNCTION_ARGS) statlseg_construct(&lseg, &path->p[iprev], &path->p[i]); tmp = lseg_closept_point(NULL, &lseg, pt); - if (!have_min || tmp < result) + if (!have_min || float8_lt(tmp, result)) { result = tmp; have_min = true; @@ -2325,19 +2347,14 @@ dist_sl(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); LINE *line = PG_GETARG_LINE_P(1); - float8 result, - d2; + float8 result; if (lseg_interpt_line(NULL, lseg, line)) result = 0.0; else - { - result = line_closept_point(NULL, line, &lseg->p[0]); - d2 = line_closept_point(NULL, line, &lseg->p[1]); /* XXX shouldn't we take the min not max? */ - if (d2 > result) - result = d2; - } + result = float8_max(line_closept_point(NULL, line, &lseg->p[0]), + line_closept_point(NULL, line, &lseg->p[1])); PG_RETURN_FLOAT8(result); } @@ -2384,11 +2401,10 @@ dist_cpoly(PG_FUNCTION_ARGS) float8 result; /* calculate distance to center, and subtract radius */ - result = dist_ppoly_internal(&circle->center, poly); - - result -= circle->radius; - if (result < 0) - result = 0; + result = float8_mi(dist_ppoly_internal(&circle->center, poly), + circle->radius); + if (result < 0.0) + result = 0.0; PG_RETURN_FLOAT8(result); } @@ -2414,7 +2430,7 @@ dist_polyp(PG_FUNCTION_ARGS) PG_RETURN_FLOAT8(dist_ppoly_internal(point, poly)); } -static double +static float8 dist_ppoly_internal(Point *pt, POLYGON *poly) { float8 result; @@ -2440,7 +2456,7 @@ dist_ppoly_internal(Point *pt, POLYGON *poly) seg.p[1].x = poly->p[i + 1].x; seg.p[1].y = poly->p[i + 1].y; d = lseg_closept_point(NULL, &seg, pt); - if (d < result) + if (float8_lt(d, result)) result = d; } @@ -2543,7 +2559,8 @@ close_pl(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); - line_closept_point(result, line, pt); + if (isnan(line_closept_point(result, line, pt))) + PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } @@ -2600,8 +2617,8 @@ static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2) { Point point; - double dist; - double d; + float8 dist, + d; d = lseg_closept_point(NULL, l1, &l2->p[0]); dist = d; @@ -2609,20 +2626,20 @@ lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2) *result = l2->p[0]; d = lseg_closept_point(NULL, l1, &l2->p[1]); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) *result = l2->p[1]; } - if (lseg_closept_point(&point, l2, &l1->p[0]) < dist) + if (float8_lt(lseg_closept_point(&point, l2, &l1->p[0]), dist)) d = lseg_closept_point(result, l1, &point); - if (lseg_closept_point(&point, l2, &l1->p[1]) < dist) + if (float8_lt(lseg_closept_point(&point, l2, &l1->p[1]), dist)) d = lseg_closept_point(result, l1, &point); - if (d < dist) + if (float8_lt(d, dist)) dist = d; return dist; @@ -2637,7 +2654,8 @@ close_lseg(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); - lseg_closept_lseg(result, l2, l1); + if (isnan(lseg_closept_lseg(result, l2, l1))) + PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } @@ -2652,11 +2670,11 @@ close_lseg(PG_FUNCTION_ARGS) static float8 box_closept_point(Point *result, BOX *box, Point *pt) { - LSEG lseg; + float8 dist, + d; Point point, closept; - double dist, - d; + LSEG lseg; if (box_contain_point(box, pt)) { @@ -2674,7 +2692,7 @@ box_closept_point(Point *result, BOX *box, Point *pt) statlseg_construct(&lseg, &box->high, &point); d = lseg_closept_point(&closept, &lseg, pt); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) @@ -2685,7 +2703,7 @@ box_closept_point(Point *result, BOX *box, Point *pt) point.y = box->low.y; statlseg_construct(&lseg, &box->low, &point); d = lseg_closept_point(&closept, &lseg, pt); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) @@ -2694,7 +2712,7 @@ box_closept_point(Point *result, BOX *box, Point *pt) statlseg_construct(&lseg, &box->high, &point); d = lseg_closept_point(&closept, &lseg, pt); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) @@ -2713,7 +2731,8 @@ close_pb(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); - box_closept_point(result, box, pt); + if (isnan(box_closept_point(result, box, pt))) + PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } @@ -2745,7 +2764,7 @@ close_sl(PG_FUNCTION_ARGS) d1 = line_closept_point(NULL, line, &lseg->p[0]); d2 = line_closept_point(NULL, line, &lseg->p[1]); - if (d1 < d2) + if (float8_lt(d1, d2)) *result = lseg->p[0]; else *result = lseg->p[1]; @@ -2809,7 +2828,8 @@ close_ls(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); - lseg_closept_line(result, lseg, line); + if (isnan(lseg_closept_line(result, lseg, line))) + PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } @@ -2824,11 +2844,11 @@ close_ls(PG_FUNCTION_ARGS) static float8 box_closept_lseg(Point *result, BOX *box, LSEG *lseg) { + float8 dist, + d; Point point, closept; LSEG bseg; - double dist, - d; if (box_interpt_lseg(result, box, lseg)) return 0.0; @@ -2841,7 +2861,7 @@ box_closept_lseg(Point *result, BOX *box, LSEG *lseg) statlseg_construct(&bseg, &box->high, &point); d = lseg_closept_lseg(&closept, &bseg, lseg); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) @@ -2852,7 +2872,7 @@ box_closept_lseg(Point *result, BOX *box, LSEG *lseg) point.y = box->low.y; statlseg_construct(&bseg, &box->low, &point); d = lseg_closept_lseg(&closept, &bseg, lseg); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) @@ -2861,7 +2881,7 @@ box_closept_lseg(Point *result, BOX *box, LSEG *lseg) statlseg_construct(&bseg, &box->high, &point); d = lseg_closept_lseg(&closept, &bseg, lseg); - if (d < dist) + if (float8_lt(d, dist)) { dist = d; if (result != NULL) @@ -2880,7 +2900,8 @@ close_sb(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); - box_closept_lseg(result, box, lseg); + if (isnan(box_closept_lseg(result, box, lseg))) + PG_RETURN_NULL(); PG_RETURN_POINT_P(result); } @@ -2913,7 +2934,9 @@ close_lb(PG_FUNCTION_ARGS) static bool line_contain_point(LINE *line, Point *point) { - return FPzero(line->A * point->x + line->B * point->y + line->C); + return FPzero(float8_pl(float8_pl(float8_mul(line->A, point->x), + float8_mul(line->B, point->y)), + line->C)); } Datum @@ -2994,7 +3017,7 @@ on_ppath(PG_FUNCTION_ARGS) PATH *path = PG_GETARG_PATH_P(1); int i, n; - double a, + float8 a, b; /*-- OPEN --*/ @@ -3005,8 +3028,7 @@ on_ppath(PG_FUNCTION_ARGS) for (i = 0; i < n; i++) { b = point_dt(pt, &path->p[i + 1]); - if (FPeq(a + b, - point_dt(&path->p[i], &path->p[i + 1]))) + if (FPeq(float8_pl(a, b), point_dt(&path->p[i], &path->p[i + 1]))) PG_RETURN_BOOL(true); a = b; } @@ -3092,10 +3114,10 @@ box_interpt_lseg(Point *result, BOX *box, LSEG *lseg) LSEG bseg; Point point; - lbox.low.x = Min(lseg->p[0].x, lseg->p[1].x); - lbox.low.y = Min(lseg->p[0].y, lseg->p[1].y); - lbox.high.x = Max(lseg->p[0].x, lseg->p[1].x); - lbox.high.y = Max(lseg->p[0].y, lseg->p[1].y); + lbox.low.x = float8_min(lseg->p[0].x, lseg->p[1].x); + lbox.low.y = float8_min(lseg->p[0].y, lseg->p[1].y); + lbox.high.x = float8_max(lseg->p[0].x, lseg->p[1].x); + lbox.high.y = float8_max(lseg->p[0].y, lseg->p[1].y); /* nothing close to overlap? then not going to intersect */ if (!box_ov(&lbox, box)) @@ -3202,7 +3224,7 @@ static void make_bound_box(POLYGON *poly) { int i; - double x1, + float8 x1, y1, x2, y2; @@ -3213,13 +3235,13 @@ make_bound_box(POLYGON *poly) y2 = y1 = poly->p[0].y; for (i = 1; i < poly->npts; i++) { - if (poly->p[i].x < x1) + if (float8_lt(poly->p[i].x, x1)) x1 = poly->p[i].x; - if (poly->p[i].x > x2) + if (float8_gt(poly->p[i].x, x2)) x2 = poly->p[i].x; - if (poly->p[i].y < y1) + if (float8_lt(poly->p[i].y, y1)) y1 = poly->p[i].y; - if (poly->p[i].y > y2) + if (float8_gt(poly->p[i].y, y2)) y2 = poly->p[i].y; } @@ -3732,8 +3754,8 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) * if X-intersection wasn't found then check central point of tested * segment. In opposite case we already check all subsegments */ - p.x = (t.p[0].x + t.p[1].x) / 2.0; - p.y = (t.p[0].y + t.p[1].y) / 2.0; + p.x = float8_div(float8_pl(t.p[0].x, t.p[1].x), 2.0); + p.y = float8_div(float8_pl(t.p[0].y, t.p[1].y), 2.0); res = point_inside(&p, poly->npts, poly->p); } @@ -3873,8 +3895,8 @@ static inline void point_add_point(Point *result, Point *pt1, Point *pt2) { point_construct(result, - pt1->x + pt2->x, - pt1->y + pt2->y); + float8_pl(pt1->x, pt2->x), + float8_pl(pt1->y, pt2->y)); } Datum @@ -3896,8 +3918,8 @@ static inline void point_sub_point(Point *result, Point *pt1, Point *pt2) { point_construct(result, - pt1->x - pt2->x, - pt1->y - pt2->y); + float8_mi(pt1->x, pt2->x), + float8_mi(pt1->y, pt2->y)); } Datum @@ -3919,8 +3941,10 @@ static inline void point_mul_point(Point *result, Point *pt1, Point *pt2) { point_construct(result, - (pt1->x * pt2->x) - (pt1->y * pt2->y), - (pt1->x * pt2->y) + (pt1->y * pt2->x)); + float8_mi(float8_mul(pt1->x, pt2->x), + float8_mul(pt1->y, pt2->y)), + float8_pl(float8_mul(pt1->x, pt2->y), + float8_mul(pt1->y, pt2->x))); } Datum @@ -3941,18 +3965,15 @@ point_mul(PG_FUNCTION_ARGS) static inline void point_div_point(Point *result, Point *pt1, Point *pt2) { - double div; + float8 div; - div = (pt2->x * pt2->x) + (pt2->y * pt2->y); - - if (div == 0.0) - ereport(ERROR, - (errcode(ERRCODE_DIVISION_BY_ZERO), - errmsg("division by zero"))); + div = float8_pl(float8_mul(pt2->x, pt2->x), float8_mul(pt2->y, pt2->y)); point_construct(result, - ((pt1->x * pt2->x) + (pt1->y * pt2->y)) / div, - ((pt2->x * pt1->y) - (pt2->y * pt1->x)) / div); + float8_div(float8_pl(float8_mul(pt1->x, pt2->x), + float8_mul(pt1->y, pt2->y)), div), + float8_div(float8_mi(float8_mul(pt1->y, pt2->x), + float8_mul(pt1->x, pt2->y)), div)); } Datum @@ -4089,10 +4110,10 @@ boxes_bound_box(PG_FUNCTION_ARGS) container = (BOX *) palloc(sizeof(BOX)); - container->high.x = Max(box1->high.x, box2->high.x); - container->low.x = Min(box1->low.x, box2->low.x); - container->high.y = Max(box1->high.y, box2->high.y); - container->low.y = Min(box1->low.y, box2->low.y); + container->high.x = float8_max(box1->high.x, box2->high.x); + container->low.x = float8_min(box1->low.x, box2->low.x); + container->high.y = float8_max(box1->high.y, box2->high.y); + container->low.y = float8_min(box1->low.y, box2->low.y); PG_RETURN_BOX_P(container); } @@ -4412,7 +4433,8 @@ circle_in(PG_FUNCTION_ARGS) s++; circle->radius = single_decode(s, &s, "circle", str); - if (circle->radius < 0) + /* We have to accept NaN. */ + if (circle->radius < 0.0) ereport(ERROR, (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), errmsg("invalid input syntax for type %s: \"%s\"", @@ -4479,7 +4501,8 @@ circle_recv(PG_FUNCTION_ARGS) circle->center.y = pq_getmsgfloat8(buf); circle->radius = pq_getmsgfloat8(buf); - if (circle->radius < 0) + /* We have to accept NaN. */ + if (circle->radius < 0.0) ereport(ERROR, (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION), errmsg("invalid radius in external \"circle\" value"))); @@ -4510,6 +4533,9 @@ circle_send(PG_FUNCTION_ARGS) *---------------------------------------------------------*/ /* circles identical? + * + * We consider NaNs values to be equal to each other to let those circles + * to be found. */ Datum circle_same(PG_FUNCTION_ARGS) @@ -4517,7 +4543,8 @@ circle_same(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPeq(circle1->radius, circle2->radius) && + PG_RETURN_BOOL(((isnan(circle1->radius) && isnan(circle1->radius)) || + FPeq(circle1->radius, circle2->radius)) && point_eq_point(&circle1->center, &circle2->center)); } @@ -4530,7 +4557,7 @@ circle_overlap(PG_FUNCTION_ARGS) CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), - circle1->radius + circle2->radius)); + float8_pl(circle1->radius, circle2->radius))); } /* circle_overleft - is the right edge of circle1 at or left of @@ -4542,8 +4569,8 @@ circle_overleft(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPle((circle1->center.x + circle1->radius), - (circle2->center.x + circle2->radius))); + PG_RETURN_BOOL(FPle(float8_pl(circle1->center.x, circle1->radius), + float8_pl(circle2->center.x, circle2->radius))); } /* circle_left - is circle1 strictly left of circle2? @@ -4554,8 +4581,8 @@ circle_left(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPlt((circle1->center.x + circle1->radius), - (circle2->center.x - circle2->radius))); + PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.x, circle1->radius), + float8_mi(circle2->center.x, circle2->radius))); } /* circle_right - is circle1 strictly right of circle2? @@ -4566,8 +4593,8 @@ circle_right(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPgt((circle1->center.x - circle1->radius), - (circle2->center.x + circle2->radius))); + PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.x, circle1->radius), + float8_pl(circle2->center.x, circle2->radius))); } /* circle_overright - is the left edge of circle1 at or right of @@ -4579,8 +4606,8 @@ circle_overright(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPge((circle1->center.x - circle1->radius), - (circle2->center.x - circle2->radius))); + PG_RETURN_BOOL(FPge(float8_mi(circle1->center.x, circle1->radius), + float8_mi(circle2->center.x, circle2->radius))); } /* circle_contained - is circle1 contained by circle2? @@ -4592,7 +4619,7 @@ circle_contained(PG_FUNCTION_ARGS) CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), - circle2->radius - circle1->radius)); + float8_mi(circle2->radius, circle1->radius))); } /* circle_contain - does circle1 contain circle2? @@ -4604,7 +4631,7 @@ circle_contain(PG_FUNCTION_ARGS) CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), - circle1->radius - circle2->radius)); + float8_mi(circle1->radius, circle2->radius))); } @@ -4616,8 +4643,8 @@ circle_below(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPlt((circle1->center.y + circle1->radius), - (circle2->center.y - circle2->radius))); + PG_RETURN_BOOL(FPlt(float8_pl(circle1->center.y, circle1->radius), + float8_mi(circle2->center.y, circle2->radius))); } /* circle_above - is circle1 strictly above circle2? @@ -4628,8 +4655,8 @@ circle_above(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPgt((circle1->center.y - circle1->radius), - (circle2->center.y + circle2->radius))); + PG_RETURN_BOOL(FPgt(float8_mi(circle1->center.y, circle1->radius), + float8_pl(circle2->center.y, circle2->radius))); } /* circle_overbelow - is the upper edge of circle1 at or below @@ -4641,8 +4668,8 @@ circle_overbelow(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPle((circle1->center.y + circle1->radius), - (circle2->center.y + circle2->radius))); + PG_RETURN_BOOL(FPle(float8_pl(circle1->center.y, circle1->radius), + float8_pl(circle2->center.y, circle2->radius))); } /* circle_overabove - is the lower edge of circle1 at or above @@ -4654,8 +4681,8 @@ circle_overabove(PG_FUNCTION_ARGS) CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); - PG_RETURN_BOOL(FPge((circle1->center.y - circle1->radius), - (circle2->center.y - circle2->radius))); + PG_RETURN_BOOL(FPge(float8_mi(circle1->center.y, circle1->radius), + float8_mi(circle2->center.y, circle2->radius))); } @@ -4768,7 +4795,7 @@ circle_mul_pt(PG_FUNCTION_ARGS) result = (CIRCLE *) palloc(sizeof(CIRCLE)); point_mul_point(&result->center, &circle->center, point); - result->radius = circle->radius * HYPOT(point->x, point->y); + result->radius = float8_mul(circle->radius, HYPOT(point->x, point->y)); PG_RETURN_CIRCLE_P(result); } @@ -4783,7 +4810,7 @@ circle_div_pt(PG_FUNCTION_ARGS) result = (CIRCLE *) palloc(sizeof(CIRCLE)); point_div_point(&result->center, &circle->center, point); - result->radius = circle->radius / HYPOT(point->x, point->y); + result->radius = float8_div(circle->radius, HYPOT(point->x, point->y)); PG_RETURN_CIRCLE_P(result); } @@ -4807,7 +4834,7 @@ circle_diameter(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); - PG_RETURN_FLOAT8(2 * circle->radius); + PG_RETURN_FLOAT8(float8_mul(circle->radius, 2.0)); } @@ -4832,10 +4859,11 @@ circle_distance(PG_FUNCTION_ARGS) CIRCLE *circle2 = PG_GETARG_CIRCLE_P(1); float8 result; - result = point_dt(&circle1->center, &circle2->center) - - (circle1->radius + circle2->radius); - if (result < 0) - result = 0; + result = float8_mi(point_dt(&circle1->center, &circle2->center), + float8_pl(circle1->radius, circle2->radius)); + if (result < 0.0) + result = 0.0; + PG_RETURN_FLOAT8(result); } @@ -4845,7 +4873,7 @@ circle_contain_pt(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); Point *point = PG_GETARG_POINT_P(1); - double d; + float8 d; d = point_dt(&circle->center, point); PG_RETURN_BOOL(d <= circle->radius); @@ -4857,7 +4885,7 @@ pt_contained_circle(PG_FUNCTION_ARGS) { Point *point = PG_GETARG_POINT_P(0); CIRCLE *circle = PG_GETARG_CIRCLE_P(1); - double d; + float8 d; d = point_dt(&circle->center, point); PG_RETURN_BOOL(d <= circle->radius); @@ -4874,9 +4902,11 @@ dist_pc(PG_FUNCTION_ARGS) CIRCLE *circle = PG_GETARG_CIRCLE_P(1); float8 result; - result = point_dt(point, &circle->center) - circle->radius; - if (result < 0) - result = 0; + result = float8_mi(point_dt(point, &circle->center), + circle->radius); + if (result < 0.0) + result = 0.0; + PG_RETURN_FLOAT8(result); } @@ -4890,9 +4920,10 @@ dist_cpoint(PG_FUNCTION_ARGS) Point *point = PG_GETARG_POINT_P(1); float8 result; - result = point_dt(point, &circle->center) - circle->radius; - if (result < 0) - result = 0; + result = float8_mi(point_dt(point, &circle->center), circle->radius); + if (result < 0.0) + result = 0.0; + PG_RETURN_FLOAT8(result); } @@ -4914,10 +4945,10 @@ circle_center(PG_FUNCTION_ARGS) /* circle_ar - returns the area of the circle. */ -static double +static float8 circle_ar(CIRCLE *circle) { - return M_PI * (circle->radius * circle->radius); + return float8_mul(float8_mul(circle->radius, circle->radius), M_PI); } @@ -4946,16 +4977,16 @@ circle_box(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); BOX *box; - double delta; + float8 delta; box = (BOX *) palloc(sizeof(BOX)); - delta = circle->radius / sqrt(2.0); + delta = float8_div(circle->radius, sqrt(2.0)); - box->high.x = circle->center.x + delta; - box->low.x = circle->center.x - delta; - box->high.y = circle->center.y + delta; - box->low.y = circle->center.y - delta; + box->high.x = float8_pl(circle->center.x, delta); + box->low.x = float8_mi(circle->center.x, delta); + box->high.y = float8_pl(circle->center.y, delta); + box->low.y = float8_mi(circle->center.y, delta); PG_RETURN_BOX_P(box); } @@ -4971,8 +5002,8 @@ box_circle(PG_FUNCTION_ARGS) circle = (CIRCLE *) palloc(sizeof(CIRCLE)); - circle->center.x = (box->high.x + box->low.x) / 2; - circle->center.y = (box->high.y + box->low.y) / 2; + circle->center.x = float8_div(float8_pl(box->high.x, box->low.x), 2.0); + circle->center.y = float8_div(float8_pl(box->high.y, box->low.y), 2.0); circle->radius = point_dt(&circle->center, &box->high); @@ -4989,8 +5020,8 @@ circle_poly(PG_FUNCTION_ARGS) int base_size, size; int i; - double angle; - double anglestep; + float8 angle; + float8 anglestep; if (FPzero(circle->radius)) ereport(ERROR, @@ -5015,13 +5046,16 @@ circle_poly(PG_FUNCTION_ARGS) SET_VARSIZE(poly, size); poly->npts = npts; - anglestep = (2.0 * M_PI) / npts; + anglestep = float8_div(2.0 * M_PI, npts); for (i = 0; i < npts; i++) { - angle = i * anglestep; - poly->p[i].x = circle->center.x - (circle->radius * cos(angle)); - poly->p[i].y = circle->center.y + (circle->radius * sin(angle)); + angle = float8_mul(anglestep, i); + + poly->p[i].x = float8_mi(circle->center.x, + float8_mul(circle->radius, cos(angle))); + poly->p[i].y = float8_pl(circle->center.y, + float8_mul(circle->radius, sin(angle))); } make_bound_box(poly); @@ -5050,12 +5084,13 @@ poly_to_circle(CIRCLE *result, POLYGON *poly) for (i = 0; i < poly->npts; i++) point_add_point(&result->center, &result->center, &poly->p[i]); - result->center.x /= poly->npts; - result->center.y /= poly->npts; + result->center.x = float8_div(result->center.x, poly->npts); + result->center.y = float8_div(result->center.y, poly->npts); for (i = 0; i < poly->npts; i++) - result->radius += point_dt(&poly->p[i], &result->center); - result->radius /= poly->npts; + result->radius = float8_pl(result->radius, + point_dt(&poly->p[i], &result->center)); + result->radius = float8_div(result->radius, poly->npts); } Datum @@ -5094,12 +5129,12 @@ poly_circle(PG_FUNCTION_ARGS) static int point_inside(Point *p, int npts, Point *plist) { - double x0, + float8 x0, y0; - double prev_x, + float8 prev_x, prev_y; int i = 0; - double x, + float8 x, y; int cross, total_cross = 0; @@ -5107,8 +5142,8 @@ point_inside(Point *p, int npts, Point *plist) Assert(npts > 0); /* compute first polygon point relative to single point */ - x0 = plist[0].x - p->x; - y0 = plist[0].y - p->y; + x0 = float8_mi(plist[0].x, p->x); + y0 = float8_mi(plist[0].y, p->y); prev_x = x0; prev_y = y0; @@ -5116,8 +5151,8 @@ point_inside(Point *p, int npts, Point *plist) for (i = 1; i < npts; i++) { /* compute next polygon point relative to single point */ - x = plist[i].x - p->x; - y = plist[i].y - p->y; + x = float8_mi(plist[i].x, p->x); + y = float8_mi(plist[i].y, p->y); /* compute previous to current point crossing */ if ((cross = lseg_crossing(x, y, prev_x, prev_y)) == POINT_ON_POLYGON) @@ -5149,9 +5184,9 @@ point_inside(Point *p, int npts, Point *plist) */ static int -lseg_crossing(double x, double y, double prev_x, double prev_y) +lseg_crossing(float8 x, float8 y, float8 prev_x, float8 prev_y) { - double z; + float8 z; int y_sign; if (FPzero(y)) @@ -5162,42 +5197,47 @@ lseg_crossing(double x, double y, double prev_x, double prev_y) { /* x > 0 */ if (FPzero(prev_y)) /* y and prev_y are zero */ /* prev_x > 0? */ - return FPgt(prev_x, 0) ? 0 : POINT_ON_POLYGON; - return FPlt(prev_y, 0) ? 1 : -1; + return FPgt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON; + return FPlt(prev_y, 0.0) ? 1 : -1; } else { /* x < 0, x not on positive X axis */ if (FPzero(prev_y)) /* prev_x < 0? */ - return FPlt(prev_x, 0) ? 0 : POINT_ON_POLYGON; + return FPlt(prev_x, 0.0) ? 0 : POINT_ON_POLYGON; return 0; } } else { /* y != 0 */ /* compute y crossing direction from previous point */ - y_sign = FPgt(y, 0) ? 1 : -1; + y_sign = FPgt(y, 0.0) ? 1 : -1; if (FPzero(prev_y)) /* previous point was on X axis, so new point is either off or on */ - return FPlt(prev_x, 0) ? 0 : y_sign; - else if (FPgt(y_sign * prev_y, 0)) + return FPlt(prev_x, 0.0) ? 0 : y_sign; + else if ((y_sign < 0 && FPlt(prev_y, 0.0)) || + (y_sign > 0 && FPgt(prev_y, 0.0))) /* both above or below X axis */ return 0; /* same sign */ else { /* y and prev_y cross X-axis */ - if (FPge(x, 0) && FPgt(prev_x, 0)) + if (FPge(x, 0.0) && FPgt(prev_x, 0.0)) /* both non-negative so cross positive X-axis */ return 2 * y_sign; - if (FPlt(x, 0) && FPle(prev_x, 0)) + if (FPlt(x, 0.0) && FPle(prev_x, 0.0)) /* both non-positive so do not cross positive X-axis */ return 0; /* x and y cross axises, see URL above point_inside() */ - z = (x - prev_x) * y - (y - prev_y) * x; + z = float8_mi(float8_mul(float8_mi(x, prev_x), y), + float8_mul(float8_mi(y, prev_y), x)); if (FPzero(z)) return POINT_ON_POLYGON; - return FPgt((y_sign * z), 0) ? 0 : 2 * y_sign; + if ((y_sign < 0 && FPlt(z, 0.0)) || + (y_sign > 0 && FPgt(z, 0.0))) + return 0; + return 2 * y_sign; } } } @@ -5265,10 +5305,11 @@ plist_same(int npts, Point *p1, Point *p2) * case of hypot(inf,nan) results in INF, and not NAN. *----------------------------------------------------------------------- */ -double -pg_hypot(double x, double y) +float8 +pg_hypot(float8 x, float8 y) { - double yx; + float8 yx, + result; /* Handle INF and NaN properly */ if (isinf(x) || isinf(y)) @@ -5284,7 +5325,7 @@ pg_hypot(double x, double y) /* Swap x and y if needed to make x the larger one */ if (x < y) { - double temp = x; + float8 temp = x; x = y; y = temp; @@ -5300,5 +5341,9 @@ pg_hypot(double x, double y) /* Determine the hypotenuse */ yx = y / x; - return x * sqrt(1.0 + (yx * yx)); + result = x * sqrt(1.0 + (yx * yx)); + + check_float8_val(result, false, false); + + return result; } diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c index fea36f361ae..4aff973ef37 100644 --- a/src/backend/utils/adt/geo_spgist.c +++ b/src/backend/utils/adt/geo_spgist.c @@ -84,14 +84,14 @@ * Comparator for qsort * * We don't need to use the floating point macros in here, because this - * is going only going to be used in a place to effect the performance + * is only going to be used in a place to effect the performance * of the index, not the correctness. */ static int compareDoubles(const void *a, const void *b) { - double x = *(double *) a; - double y = *(double *) b; + float8 x = *(float8 *) a; + float8 y = *(float8 *) b; if (x == y) return 0; @@ -100,8 +100,8 @@ compareDoubles(const void *a, const void *b) typedef struct { - double low; - double high; + float8 low; + float8 high; } Range; typedef struct @@ -175,7 +175,7 @@ static RectBox * initRectBox(void) { RectBox *rect_box = (RectBox *) palloc(sizeof(RectBox)); - double infinity = get_float8_infinity(); + float8 infinity = get_float8_infinity(); rect_box->range_box_x.left.low = -infinity; rect_box->range_box_x.left.high = infinity; @@ -418,10 +418,10 @@ spg_box_quad_picksplit(PG_FUNCTION_ARGS) BOX *centroid; int median, i; - double *lowXs = palloc(sizeof(double) * in->nTuples); - double *highXs = palloc(sizeof(double) * in->nTuples); - double *lowYs = palloc(sizeof(double) * in->nTuples); - double *highYs = palloc(sizeof(double) * in->nTuples); + float8 *lowXs = palloc(sizeof(float8) * in->nTuples); + float8 *highXs = palloc(sizeof(float8) * in->nTuples); + float8 *lowYs = palloc(sizeof(float8) * in->nTuples); + float8 *highYs = palloc(sizeof(float8) * in->nTuples); /* Calculate median of all 4D coordinates */ for (i = 0; i < in->nTuples; i++) @@ -434,10 +434,10 @@ spg_box_quad_picksplit(PG_FUNCTION_ARGS) highYs[i] = box->high.y; } - qsort(lowXs, in->nTuples, sizeof(double), compareDoubles); - qsort(highXs, in->nTuples, sizeof(double), compareDoubles); - qsort(lowYs, in->nTuples, sizeof(double), compareDoubles); - qsort(highYs, in->nTuples, sizeof(double), compareDoubles); + qsort(lowXs, in->nTuples, sizeof(float8), compareDoubles); + qsort(highXs, in->nTuples, sizeof(float8), compareDoubles); + qsort(lowYs, in->nTuples, sizeof(float8), compareDoubles); + qsort(highYs, in->nTuples, sizeof(float8), compareDoubles); median = in->nTuples / 2; diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h index 0e066894cd1..9f8505804e8 100644 --- a/src/include/utils/geo_decls.h +++ b/src/include/utils/geo_decls.h @@ -8,9 +8,6 @@ * * src/include/utils/geo_decls.h * - * NOTE - * These routines do *not* use the float types from adt/. - * * XXX These routines were not written by a numerical analyst. * * XXX I have made some attempt to flesh out the operators @@ -21,14 +18,14 @@ #ifndef GEO_DECLS_H #define GEO_DECLS_H -#include - #include "fmgr.h" /*-------------------------------------------------------------------- * Useful floating point utilities and constants. - *-------------------------------------------------------------------*/ - + *------------------------------------------------------------------- + * + * XXX: They are not NaN-aware. + */ #define EPSILON 1.0E-06 @@ -57,7 +54,7 @@ *-------------------------------------------------------------------*/ typedef struct { - double x, + float8 x, y; } Point; @@ -89,7 +86,7 @@ typedef struct *-------------------------------------------------------------------*/ typedef struct { - double A, + float8 A, B, C; } LINE; @@ -124,7 +121,7 @@ typedef struct typedef struct { Point center; - double radius; + float8 radius; } CIRCLE; /* @@ -178,6 +175,6 @@ typedef struct * in geo_ops.c */ -extern double pg_hypot(double x, double y); +extern float8 pg_hypot(float8 x, float8 y); #endif /* GEO_DECLS_H */