From 89b5669a201de0cb597e3facda51eccbc3884e8e Mon Sep 17 00:00:00 2001 From: Emre Hasegeli Date: Sat, 27 May 2017 16:17:42 -0400 Subject: [PATCH 3/6] geo-float-v08 Use the built-in float datatype to implement geometric types The changes can be summarised as: * Check for underflow and overflow * Check for division by zero * Let the objects with NaN values to be the same * Return NULL when the distance is NaN for all closest point operators * Don't round the slope to DBL_MAX * Use NaN as the slope of two equal points * Favour not-NaN over NaN where it makes sense The patch also replaces all occurrences of "double" as "float8". They are the same, but were randomly spread on the same file. --- src/backend/access/gist/gistproc.c | 91 +++---- src/backend/utils/adt/geo_ops.c | 494 ++++++++++++++++++++----------------- src/backend/utils/adt/geo_spgist.c | 28 +-- src/include/utils/geo_decls.h | 17 +- 4 files changed, 331 insertions(+), 299 deletions(-) diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c index 4bbfdadf9f..4b6233f212 100644 --- a/src/backend/access/gist/gistproc.c +++ b/src/backend/access/gist/gistproc.c @@ -16,20 +16,21 @@ *------------------------------------------------------------------------- */ #include "postgres.h" #include #include #include "access/gist.h" #include "access/stratnum.h" #include "utils/builtins.h" +#include "utils/float.h" #include "utils/geo_decls.h" static bool gist_box_leaf_consistent(BOX *key, BOX *query, StrategyNumber strategy); static bool rtree_internal_consistent(BOX *key, BOX *query, StrategyNumber strategy); /* Minimum accepted ratio of split */ #define LIMIT_RATIO 0.3 @@ -48,55 +49,56 @@ rt_box_union(BOX *n, const BOX *a, const BOX *b) n->high.x = float8_max(a->high.x, b->high.x); n->high.y = float8_max(a->high.y, b->high.y); n->low.x = float8_min(a->low.x, b->low.x); n->low.y = float8_min(a->low.y, b->low.y); } /* * 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) { /* * Check for zero-width cases. Note that we define the size of a zero- * by-infinity box as zero. It's important to special-case this somehow, * as naively multiplying infinity by zero will produce NaN. * * The less-than cases should not happen, but if they do, say "zero". */ if (float8_le(box->high.x, box->low.x) || float8_le(box->high.y, box->low.y)) return 0.0; /* * We treat NaN as larger than +Infinity, so any distance involving a NaN * and a non-NaN is infinite. Note the previous check eliminated the * possibility that the low fields are NaNs. */ 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)); } /* * The GiST Consistent method for boxes * * Should return false if for all data items x below entry, * the predicate x op query must be false, where op is the oper * corresponding to strategy in the pg_amop table. */ Datum @@ -256,74 +258,74 @@ fallbackSplit(GistEntryVector *entryvec, GIST_SPLITVEC *v) /* * Represents information about an entry that can be placed to either group * without affecting overlap over selected axis ("common entry"). */ 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; /* * Context for g_box_consider_split. Contains information about currently * selected split and some general information. */ typedef struct { int entriesCount; /* total number of entries being split */ BOX boundingBox; /* minimum bounding box across all entries */ /* Information about currently selected split follows */ 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; /* * Interval represents projection of box to axis. */ typedef struct { - double lower, + float8 lower, upper; } SplitInterval; /* * Interval comparison function by lower bound of the interval; */ 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); } /* * Interval comparison function by upper bound of the interval; */ 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); } /* * Replace negative (or NaN) value with zero. */ static inline float non_negative(float val) @@ -332,28 +334,28 @@ non_negative(float val) return val; else return 0.0f; } /* * Consider replacement of currently selected split with the better one. */ 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 * of common entries. */ if (minLeftCount >= (context->entriesCount + 1) / 2) { leftCount = minLeftCount; } else @@ -362,40 +364,41 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum, leftCount = maxLeftCount; else leftCount = context->entriesCount / 2; } rightCount = context->entriesCount - leftCount; /* * 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) { bool selectthis = false; /* * The ratio is acceptable, so compare current split with previously * selected one. Between splits of one dimension we search for minimal * overlap (allowing negative values) and minimal ration (between same * overlaps. We switch dimension if find less overlap (non-negative) * 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) selectthis = true; else if (context->dim == dimNum) { /* * Within the same dimension, choose the new split if it has a * smaller overlap, or same overlap but better ratio. */ @@ -524,21 +527,21 @@ gist_box_picksplit(PG_FUNCTION_ARGS) else adjustBox(&context.boundingBox, box); } /* * Iterate over axes for optimal split searching. */ context.first = true; /* nothing selected yet */ for (dim = 0; dim < 2; dim++) { - double leftUpper, + float8 leftUpper, rightLower; int i1, i2; /* Project each entry as an interval on the selected axis. */ for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i)) { box = DatumGetBoxP(entryvec->vector[i].key); if (dim == 0) { @@ -721,21 +724,21 @@ gist_box_picksplit(PG_FUNCTION_ARGS) *rightBox = *(box); \ v->spl_right[v->spl_nright++] = off; \ } while(0) /* * Distribute entries which can be distributed unambiguously, and collect * common entries. */ for (i = FirstOffsetNumber; i <= maxoff; i = OffsetNumberNext(i)) { - double lower, + float8 lower, upper; /* * Get upper and lower bounds along selected axis. */ box = DatumGetBoxP(entryvec->vector[i].key); if (context.dim == 0) { lower = box->low.x; upper = box->high.x; @@ -776,31 +779,31 @@ gist_box_picksplit(PG_FUNCTION_ARGS) /* * Distribute "common entries", if any. */ if (commonEntriesCount > 0) { /* * 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 * (float8) nentries); /* * Calculate delta between penalties of join "common entries" to * different groups. */ 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))); } /* * Sort "common entries" by calculated deltas in order to distribute * the most ambiguous entries first. */ qsort(commonEntries, commonEntriesCount, sizeof(CommonEntry), common_entry_cmp); /* * Distribute "common entries" between groups. @@ -1100,24 +1103,24 @@ gist_circle_compress(PG_FUNCTION_ARGS) { GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); GISTENTRY *retval; if (entry->leafkey) { CIRCLE *in = DatumGetCircleP(entry->key); 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), entry->rel, entry->page, entry->offset, false); } else retval = entry; PG_RETURN_POINTER(retval); } @@ -1141,24 +1144,24 @@ gist_circle_consistent(PG_FUNCTION_ARGS) *recheck = true; if (DatumGetBoxP(entry->key) == NULL || query == NULL) PG_RETURN_BOOL(false); /* * Since the operators require recheck anyway, we can just use * 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); PG_RETURN_BOOL(result); } /************************************************** * Point ops **************************************************/ @@ -1209,63 +1212,63 @@ gist_point_fetch(PG_FUNCTION_ARGS) entry->offset, false); PG_RETURN_POINTER(retval); } #define point_point_distance(p1,p2) \ 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) { /* simple point to point distance */ result = point_point_distance(point, &box->low); } else if (point->x <= box->high.x && point->x >= box->low.x && point->y <= box->high.y && point->y >= box->low.y) { /* point inside the box */ result = 0.0; } else if (point->x <= box->high.x && point->x >= box->low.x) { /* 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"); } else if (point->y <= box->high.y && point->y >= box->low.y) { /* 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"); } else { /* closest point will be a vertex */ Point p; - double subresult; + float8 subresult; result = point_point_distance(point, &box->low); subresult = point_point_distance(point, &box->high); if (result > subresult) result = subresult; p.x = box->low.x; p.y = box->high.y; subresult = point_point_distance(point, &p); @@ -1442,21 +1445,21 @@ gist_point_consistent(PG_FUNCTION_ARGS) } PG_RETURN_BOOL(result); } Datum 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) { case PointStrategyNumberGroup: distance = computeDistance(GIST_LEAF(entry), DatumGetBoxP(entry->key), PG_GETARG_POINT_P(1)); break; default: @@ -1471,25 +1474,25 @@ gist_point_distance(PG_FUNCTION_ARGS) /* * The inexact GiST distance method for geometric types that store bounding * boxes. * * Compute lossy distance from point to index entries. The result is inexact * because index entries are bounding boxes, not the exact shapes of the * indexed geometric types. We use distance from point to MBR of index entry. * 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. */ *recheck = true; switch (strategyGroup) { case PointStrategyNumberGroup: distance = computeDistance(false, DatumGetBoxP(entry->key), @@ -1505,32 +1508,32 @@ gist_bbox_distance(GISTENTRY *entry, Datum query, Datum gist_circle_distance(PG_FUNCTION_ARGS) { GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); Datum query = PG_GETARG_DATUM(1); StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2); /* 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); PG_RETURN_FLOAT8(distance); } Datum gist_poly_distance(PG_FUNCTION_ARGS) { GISTENTRY *entry = (GISTENTRY *) PG_GETARG_POINTER(0); Datum query = PG_GETARG_DATUM(1); StrategyNumber strategy = (StrategyNumber) PG_GETARG_UINT16(2); /* 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); PG_RETURN_FLOAT8(distance); } diff --git a/src/backend/utils/adt/geo_ops.c b/src/backend/utils/adt/geo_ops.c index 16bda3087c..2334da10fe 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -49,56 +49,56 @@ static int point_inside(Point *p, int npts, Point *plist); static inline void line_construct_pm(LINE *result, Point *pt, float8 m); static inline void line_construct_pts(LINE *result, Point *pt1, Point *pt2); static bool line_interpt_line(Point *result, LINE *l1, LINE *l2); static bool line_contain_point(LINE *line, Point *point); static float8 line_closept_point(Point *result, LINE *line, Point *pt); /* Routines for two-dimensional line segments */ static inline void statlseg_construct(LSEG *lseg, Point *pt1, Point *pt2); 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); static float8 lseg_closept_lseg(Point *result, LSEG *l1, LSEG *l2); /* Routines for two-dimensional boxes */ 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); static bool box_interpt_lseg(Point *result, BOX *box, LSEG *lseg); static float8 box_closept_point(Point *result, BOX *box, Point *point); static float8 box_closept_lseg(Point *result, BOX *box, LSEG *lseg); /* Routines for two-dimensional circles */ -static double circle_ar(CIRCLE *circle); +static float8 circle_ar(CIRCLE *circle); /* Routines for two-dimensional polygons */ static void make_bound_box(POLYGON *poly); static void poly_to_circle(CIRCLE *result, POLYGON *poly); static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start); static bool poly_contain_poly(POLYGON *polya, POLYGON *polyb); 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); static void path_decode(char *str, bool opentype, int npts, Point *p, bool *isopen, char **endptr_p, const char *type_name, const char *orig_string); static char *path_encode(enum path_delim path_delim, int npts, Point *pt); /* @@ -136,38 +136,38 @@ static char *path_encode(enum path_delim path_delim, int npts, Point *pt); * * For boxes, the points are opposite corners with the first point at the top right. * For closed paths and polygons, the points should be reordered to allow * fast and correct equality comparisons. * * XXX perhaps points in complex shapes should be reordered internally * to allow faster internal operations, but should keep track of input order * 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) { return float8in_internal(num, endptr_p, type_name, orig_string); } /* single_decode() */ static void single_encode(float8 x, StringInfo str) { char *xstr = float8out_internal(x); appendStringInfoString(str, xstr); pfree(xstr); } /* 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; while (isspace((unsigned char) *str)) str++; if ((has_delim = (*str == LDELIM))) str++; *x = float8in_internal(str, &str, type_name, orig_string); @@ -366,33 +366,33 @@ pair_count(char *s, char delim) * External format: (two corners of box) * "(f8, f8), (f8, f8)" * also supports the older style "(f8, f8, f8, f8)" */ Datum 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; box->low.y = y; } PG_RETURN_BOX_P(box); } /* box_out - convert a box to external form. @@ -406,38 +406,38 @@ box_out(PG_FUNCTION_ARGS) } /* * box_recv - converts external binary format to box */ Datum box_recv(PG_FUNCTION_ARGS) { StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); BOX *box; - double x, + float8 x, y; box = (BOX *) palloc(sizeof(BOX)); box->high.x = pq_getmsgfloat8(buf); box->high.y = pq_getmsgfloat8(buf); box->low.x = pq_getmsgfloat8(buf); 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; box->low.y = y; } PG_RETURN_BOX_P(box); } /* @@ -456,31 +456,31 @@ box_send(PG_FUNCTION_ARGS) pq_sendfloat8(&buf, box->low.y); PG_RETURN_BYTEA_P(pq_endtypsend(&buf)); } /* box_construct - fill in a new box. */ 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; } else { 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; } else { result->high.y = pt2->y; result->low.y = pt1->y; } } @@ -798,54 +798,54 @@ box_center(PG_FUNCTION_ARGS) Point *result = (Point *) palloc(sizeof(Point)); box_cn(result, box); PG_RETURN_POINT_P(result); } /* 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)); } /* box_cn - stores the centerpoint of the box into *center. */ 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); } /*---------------------------------------------------------- * Funky operations. *---------------------------------------------------------*/ /* box_intersect - * returns the overlapping portion of two boxes, * or NULL if they do not intersect. @@ -855,24 +855,24 @@ box_intersect(PG_FUNCTION_ARGS) { BOX *box1 = PG_GETARG_BOX_P(0); BOX *box2 = PG_GETARG_BOX_P(1); BOX *result; if (!box_ov(box1, box2)) PG_RETURN_NULL(); 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); } /* box_diagonal - * returns a line segment which happens to be the * positive-slope diagonal of "box". */ Datum @@ -1004,37 +1004,37 @@ line_send(PG_FUNCTION_ARGS) /* * Fill already-allocated LINE struct from the point and the slope */ static inline void line_construct_pm(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 if (m == 0.0) { /* horizontal - use "y = C" */ result->A = 0.0; result->B = -1.0; result->C = pt->y; } else { /* 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)); } } /* * Fill already-allocated LINE struct from two points on the line */ static inline void line_construct_pts(LINE *result, Point *pt1, Point *pt2) { float8 m; @@ -1073,35 +1073,36 @@ line_intersect(PG_FUNCTION_ARGS) Datum line_parallel(PG_FUNCTION_ARGS) { LINE *l1 = PG_GETARG_LINE_P(0); LINE *l2 = PG_GETARG_LINE_P(1); if (FPzero(l1->B)) PG_RETURN_BOOL(FPzero(l2->B)); - PG_RETURN_BOOL(FPeq(l2->A, l1->A * (l2->B / l1->B))); + PG_RETURN_BOOL(FPeq(l2->A, float8_mul(l1->A, float8_div(l2->B, l1->B)))); } 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)) 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 line_vertical(PG_FUNCTION_ARGS) { LINE *line = PG_GETARG_LINE_P(0); PG_RETURN_BOOL(FPzero(line->B)); } @@ -1111,34 +1112,37 @@ line_horizontal(PG_FUNCTION_ARGS) LINE *line = PG_GETARG_LINE_P(0); 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); - 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(((isnan(l1->A) && isnan(l2->A)) || + FPeq(l1->A, float8_mul(ratio, l2->A))) && + ((isnan(l1->B) && isnan(l2->B)) || + FPeq(l1->B, float8_mul(ratio, l2->B))) && + ((isnan(l1->C) && isnan(l2->C)) || + FPeq(l1->C, float8_mul(ratio, l2->C)))); } /*---------------------------------------------------------- * Line arithmetic routines. *---------------------------------------------------------*/ /* line_distance() * Distance between two lines. */ @@ -1146,21 +1150,21 @@ Datum line_distance(PG_FUNCTION_ARGS) { LINE *l1 = PG_GETARG_LINE_P(0); LINE *l2 = PG_GETARG_LINE_P(1); float8 result; Point tmp; if (line_interpt_line(NULL, l1, l2)) 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); } /* line_interpt() * Point where two lines l1, l2 intersect (if any) */ Datum line_interpt(PG_FUNCTION_ARGS) @@ -1179,47 +1183,51 @@ line_interpt(PG_FUNCTION_ARGS) /* * Internal version of line_interpt * * This returns true if two lines intersect (they do, if they are not * parallel), false if they do not. This also sets the intersection point * to *result, if it is not NULL. * * 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? */ { if (FPzero(l2->B)) /* l2 vertical? */ 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) point_construct(result, x, y); return true; } /*********************************************************************** ** @@ -1240,36 +1248,35 @@ line_interpt_line(Point *result, LINE *l1, LINE *l2) * "(xcoord, ycoord),... " * "[xcoord, ycoord,... ]" * Also support older format: * "(closed, npts, xcoord, ycoord,... )" *---------------------------------------------------------*/ Datum path_area(PG_FUNCTION_ARGS) { PATH *path = PG_GETARG_PATH_P(0); - double area = 0.0; + float8 area = 0.0; int i, j; if (!path->closed) PG_RETURN_NULL(); 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)); } Datum path_in(PG_FUNCTION_ARGS) { char *str = PG_GETARG_CSTRING(0); PATH *path; bool isopen; char *s; @@ -1525,33 +1532,33 @@ path_inter(PG_FUNCTION_ARGS) j; LSEG seg1, seg2; Assert(p1->npts > 0 && p2->npts > 0); b1.high.x = b1.low.x = p1->p[0].x; 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); /* pairwise check lseg intersections */ for (i = 0; i < p1->npts; i++) { int iprev; if (i > 0) @@ -1627,21 +1634,21 @@ path_distance(PG_FUNCTION_ARGS) { if (!p2->closed) continue; jprev = p2->npts - 1; /* include the closure segment */ } statlseg_construct(&seg1, &p1->p[iprev], &p1->p[i]); 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; } } } if (!have_min) PG_RETURN_NULL(); @@ -1666,21 +1673,21 @@ path_length(PG_FUNCTION_ARGS) if (i > 0) iprev = i - 1; else { if (!path->closed) continue; 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); } /*********************************************************************** ** ** Routines for 2D points. ** ***********************************************************************/ @@ -1832,41 +1839,42 @@ point_ne(PG_FUNCTION_ARGS) { Point *pt1 = PG_GETARG_POINT_P(0); Point *pt2 = PG_GETARG_POINT_P(1); PG_RETURN_BOOL(!point_eq_point(pt1, pt2)); } static inline bool point_eq_point(Point *pt1, Point *pt2) { - return FPeq(pt1->x, pt2->x) && FPeq(pt1->y, pt2->y); + return ((isnan(pt1->x) && isnan(pt2->x)) || FPeq(pt1->x, pt2->x)) && + ((isnan(pt1->y) && isnan(pt2->y)) || FPeq(pt1->y, pt2->y)); } /*---------------------------------------------------------- * "Arithmetic" operators on points. *---------------------------------------------------------*/ Datum point_distance(PG_FUNCTION_ARGS) { Point *pt1 = PG_GETARG_POINT_P(0); Point *pt2 = PG_GETARG_POINT_P(1); PG_RETURN_FLOAT8(point_dt(pt1, pt2)); } static 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 point_slope(PG_FUNCTION_ARGS) { Point *pt1 = PG_GETARG_POINT_P(0); Point *pt2 = PG_GETARG_POINT_P(1); PG_RETURN_FLOAT8(slope(pt1->x, pt2->x, pt1->y, pt2->y)); } @@ -1875,27 +1883,29 @@ point_slope(PG_FUNCTION_ARGS) /* * Return slope of two points * * This function accepts x and y coordinates separately to let it be used * to calculate inverse slope. To achieve that, pass the values in * (y1, y2, x2, x1) order. */ static inline float8 slope(float8 x1, float8 x2, float8 y1, float8 y2) { - /* - * XXX This should be exact checking. The callers are using the macros - * when necessary. - */ - if (FPeq(x1, x2)) + if (x1 == x2) + { + if (y1 == y2) + return NAN; + return DBL_MAX; - return (y1 - y2) / (x1 - x2); + } + + return float8_div(float8_mi(y1, y2), float8_mi(x1, x2)); } /*********************************************************************** ** ** Routines for 2D line segments. ** ***********************************************************************/ /*---------------------------------------------------------- @@ -2035,21 +2045,21 @@ lseg_parallel(PG_FUNCTION_ARGS) * '((0,0),(0,1))'::lseg ?-| '((0,0),(1,0))'::lseg * So, modified it to check explicitly for slope of vertical line * returned by slope() and the results seem better. * - thomas 1998-01-31 */ Datum lseg_perp(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); - double m1, + float8 m1, m2; m1 = slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y); m2 = slope(l2->p[0].x, l2->p[1].x, l2->p[0].y, l2->p[1].y); #ifdef GEODEBUG printf("lseg_perp- slopes are %g and %g\n", m1, m2); #endif if (FPzero(m1)) PG_RETURN_BOOL(FPeq(m2, DBL_MAX)); @@ -2157,22 +2167,22 @@ lseg_distance(PG_FUNCTION_ARGS) Datum lseg_center(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); Point *result; 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); } /* * Find the intersection point of two segments (if any). * * This returns true if two line segments intersect, false if they do not. * This also sets the intersection point to *result, if it is not NULL. @@ -2281,21 +2291,21 @@ dist_ppath(PG_FUNCTION_ARGS) iprev = i - 1; else { if (!path->closed) continue; iprev = path->npts - 1; /* Include the closure segment */ } 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; } } PG_RETURN_FLOAT8(result); } /* @@ -2311,33 +2321,28 @@ dist_pb(PG_FUNCTION_ARGS) } /* * Distance from a lseg to a line */ Datum 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); } /* * Distance from a lseg to a box */ Datum dist_sb(PG_FUNCTION_ARGS) { @@ -2370,25 +2375,24 @@ dist_lb(PG_FUNCTION_ARGS) * Distance from a circle to a polygon */ Datum dist_cpoly(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); POLYGON *poly = PG_GETARG_POLYGON_P(1); 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); } /* * Distance from a point to a polygon */ Datum dist_ppoly(PG_FUNCTION_ARGS) { @@ -2437,21 +2441,21 @@ dist_ppoly_internal(Point *pt, POLYGON *poly) for (i = 0; i < poly->npts - 1; i++) { seg.p[0].x = poly->p[i].x; seg.p[0].y = poly->p[i].y; 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); #ifdef GEODEBUG printf("dist_ppoly_internal- segment %d distance is %f\n", (i + 1), d); #endif - if (d < result) + if (float8_lt(d, result)) result = d; } return result; } /*--------------------------------------------------------------------- * interpt_ * Intersection point of objects. @@ -2532,36 +2536,37 @@ line_closept_point(Point *result, LINE *line, Point *point) Datum close_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); Point *result; 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); } /* * Closest point on line segment to specified point. * * This sets the closest point to the *result if it is not NULL and returns * the distance to the closest point. */ static float8 lseg_closept_point(Point *result, LSEG *lseg, Point *pt) { - double invm; + float8 invm; Point closept; LINE tmp; /* * To find the closest point, we draw a perpendicular line from the point * to the line segment. */ invm = slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x); line_construct_pm(&tmp, pt, invm); lseg_closept_line(&closept, lseg, &tmp); @@ -2594,134 +2599,136 @@ close_ps(PG_FUNCTION_ARGS) * This sets the closest point to the *result if it is not NULL and returns * the distance to the closest point. * * XXX This function is wrong. If must never set the *result to a point on * the second segment. */ 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; if (result != NULL) *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; } Datum close_lseg(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); Point *result; 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); } /* * Closest point on or in box to specified point. * * This sets the closest point to the *result if it is not NULL and returns * the distance to the closest point. */ 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)) { if (result != NULL) *result = *pt; return 0.0; } /* pairwise check lseg distances */ point.x = box->low.x; point.y = box->high.y; statlseg_construct(&lseg, &box->low, &point); dist = lseg_closept_point(result, &lseg, 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) *result = closept; } point.x = box->high.x; 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) *result = closept; } 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) *result = closept; } return dist; } Datum close_pb(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); BOX *box = PG_GETARG_BOX_P(1); Point *result; 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); } /* close_sl() * Closest point on line to line segment. * * XXX THIS CODE IS WRONG * The code is actually calculating the point on the line segment @@ -2739,21 +2746,21 @@ close_sl(PG_FUNCTION_ARGS) float8 d1, d2; result = (Point *) palloc(sizeof(Point)); if (lseg_interpt_line(result, lseg, line)) PG_RETURN_POINT_P(result); 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]; PG_RETURN_POINT_P(result); #endif ereport(ERROR, (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), errmsg("function \"close_sl\" not implemented"))); @@ -2803,92 +2810,94 @@ lseg_closept_line(Point *result, LSEG *lseg, LINE *line) Datum close_ls(PG_FUNCTION_ARGS) { LINE *line = PG_GETARG_LINE_P(0); LSEG *lseg = PG_GETARG_LSEG_P(1); Point *result; 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); } /* * Closest point on or in box to line segment. * * This sets the closest point to the *result if it is not NULL and returns * the distance to the closest point. */ 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; /* pairwise check lseg distances */ point.x = box->low.x; point.y = box->high.y; statlseg_construct(&bseg, &box->low, &point); dist = lseg_closept_lseg(result, &bseg, 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) *result = closept; } point.x = box->high.x; 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) *result = closept; } 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) *result = closept; } return dist; } Datum close_sb(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); BOX *box = PG_GETARG_BOX_P(1); Point *result; 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); } Datum close_lb(PG_FUNCTION_ARGS) { #ifdef NOT_USED LINE *line = PG_GETARG_LINE_P(0); @@ -2907,21 +2916,23 @@ close_lb(PG_FUNCTION_ARGS) * on_ * Whether one object lies completely within another. *-------------------------------------------------------------------*/ /* * Does the point satisfy the equation? */ 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 on_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); PG_RETURN_BOOL(line_contain_point(line, pt)); } @@ -2989,33 +3000,32 @@ box_contain_pt(PG_FUNCTION_ARGS) * but not cross. * (we can do p-in-p in lg(n), but it takes preprocessing) */ Datum on_ppath(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); PATH *path = PG_GETARG_PATH_P(1); int i, n; - double a, + float8 a, b; /*-- OPEN --*/ if (!path->closed) { n = path->npts - 1; a = point_dt(pt, &path->p[0]); 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; } PG_RETURN_BOOL(false); } /*-- CLOSED --*/ PG_RETURN_BOOL(point_inside(pt, path->npts, path->p) != 0); } @@ -3087,24 +3097,24 @@ inter_sl(PG_FUNCTION_ARGS) * Optimize for non-intersection by checking for box intersection first. * - thomas 1998-01-30 */ static bool box_interpt_lseg(Point *result, BOX *box, LSEG *lseg) { BOX lbox; 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)) return false; if (result != NULL) { box_cn(&point, box); lseg_closept_point(result, lseg, &point); } @@ -3197,38 +3207,38 @@ inter_lb(PG_FUNCTION_ARGS) * make_bound_box - create the bounding box for the input polygon *------------------------------------------------------------------*/ /*--------------------------------------------------------------------- * Make the smallest bounding box for the given polygon. *---------------------------------------------------------------------*/ static void make_bound_box(POLYGON *poly) { int i; - double x1, + float8 x1, y1, x2, y2; Assert(poly->npts > 0); x1 = x2 = poly->p[0].x; 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; } poly->boundbox.low.x = x1; poly->boundbox.high.x = x2; poly->boundbox.low.y = y1; poly->boundbox.high.y = y2; } /*------------------------------------------------------------------ @@ -3726,22 +3736,22 @@ lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start) } if (res && !intersection) { Point p; /* * 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); } return res; } /*----------------------------------------------------------------- * Determine if polygon A contains polygon B. *-----------------------------------------------------------------*/ @@ -3867,22 +3877,22 @@ construct_point(PG_FUNCTION_ARGS) point_construct(result, x, y); PG_RETURN_POINT_P(result); } 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 point_add(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; result = (Point *) palloc(sizeof(Point)); @@ -3890,22 +3900,22 @@ point_add(PG_FUNCTION_ARGS) point_add_point(result, p1, p2); PG_RETURN_POINT_P(result); } 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 point_sub(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; result = (Point *) palloc(sizeof(Point)); @@ -3913,54 +3923,53 @@ point_sub(PG_FUNCTION_ARGS) point_sub_point(result, p1, p2); PG_RETURN_POINT_P(result); } 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 point_mul(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; result = (Point *) palloc(sizeof(Point)); point_mul_point(result, p1, p2); PG_RETURN_POINT_P(result); } 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 point_div(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; result = (Point *) palloc(sizeof(Point)); @@ -4083,24 +4092,24 @@ point_box(PG_FUNCTION_ARGS) */ Datum boxes_bound_box(PG_FUNCTION_ARGS) { BOX *box1 = PG_GETARG_BOX_P(0), *box2 = PG_GETARG_BOX_P(1), *container; 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); } /*********************************************************************** ** ** Routines for 2D paths. ** ***********************************************************************/ @@ -4406,21 +4415,22 @@ circle_in(PG_FUNCTION_ARGS) if (*cp == LDELIM) s = cp; } pair_decode(s, &circle->center.x, &circle->center.y, &s, "circle", str); if (*s == DELIM) s++; circle->radius = single_decode(s, &s, "circle", str); - if (circle->radius < 0) + /* We have to accept NaN as well. */ + if (circle->radius < 0.0) ereport(ERROR, (errcode(ERRCODE_INVALID_TEXT_REPRESENTATION), errmsg("invalid input syntax for type %s: \"%s\"", "circle", str))); while (depth > 0) { if ((*s == RDELIM) || ((*s == RDELIM_C) && (depth == 1))) { depth--; @@ -4473,21 +4483,21 @@ circle_recv(PG_FUNCTION_ARGS) { StringInfo buf = (StringInfo) PG_GETARG_POINTER(0); CIRCLE *circle; circle = (CIRCLE *) palloc(sizeof(CIRCLE)); circle->center.x = pq_getmsgfloat8(buf); circle->center.y = pq_getmsgfloat8(buf); circle->radius = pq_getmsgfloat8(buf); - if (circle->radius < 0) + if (circle->radius <= 0.0) ereport(ERROR, (errcode(ERRCODE_INVALID_BINARY_REPRESENTATION), errmsg("invalid radius in external \"circle\" value"))); PG_RETURN_CIRCLE_P(circle); } /* * circle_send - converts circle to binary format */ @@ -4511,159 +4521,160 @@ circle_send(PG_FUNCTION_ARGS) *---------------------------------------------------------*/ /* circles identical? */ Datum 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)); } /* circle_overlap - does circle1 overlap circle2? */ Datum circle_overlap(PG_FUNCTION_ARGS) { CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); 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 * the right edge of circle2? */ Datum 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? */ Datum 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? */ Datum 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 * the left edge of circle2? */ Datum 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? */ Datum circle_contained(PG_FUNCTION_ARGS) { CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); 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? */ Datum circle_contain(PG_FUNCTION_ARGS) { CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); 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))); } /* circle_below - is circle1 strictly below circle2? */ Datum 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? */ Datum 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 * the upper edge of circle2? */ Datum 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 * the lower edge of circle2? */ Datum 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))); } /* circle_relop - is area(circle1) relop area(circle2), within * our accuracy constraint? */ Datum circle_eq(PG_FUNCTION_ARGS) { CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); @@ -4762,36 +4773,36 @@ circle_sub_pt(PG_FUNCTION_ARGS) Datum circle_mul_pt(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); Point *point = PG_GETARG_POINT_P(1); CIRCLE *result; 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); } Datum circle_div_pt(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); Point *point = PG_GETARG_POINT_P(1); CIRCLE *result; 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); } /* circle_area - returns the area of the circle. */ Datum circle_area(PG_FUNCTION_ARGS) { @@ -4801,21 +4812,21 @@ circle_area(PG_FUNCTION_ARGS) } /* circle_diameter - returns the diameter of the circle. */ Datum 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)); } /* circle_radius - returns the radius of the circle. */ Datum circle_radius(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); @@ -4826,81 +4837,85 @@ circle_radius(PG_FUNCTION_ARGS) /* circle_distance - returns the distance between * two circles. */ Datum circle_distance(PG_FUNCTION_ARGS) { CIRCLE *circle1 = PG_GETARG_CIRCLE_P(0); 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); } Datum 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); } Datum 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); } /* dist_pc - returns the distance between * a point and a circle. */ Datum dist_pc(PG_FUNCTION_ARGS) { Point *point = PG_GETARG_POINT_P(0); 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); } /* * Distance from a circle to a point */ Datum dist_cpoint(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); 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); } /* circle_center - returns the center point of the circle. */ Datum circle_center(PG_FUNCTION_ARGS) { CIRCLE *circle = PG_GETARG_CIRCLE_P(0); Point *result; @@ -4908,24 +4923,24 @@ circle_center(PG_FUNCTION_ARGS) result = (Point *) palloc(sizeof(Point)); result->x = circle->center.x; result->y = circle->center.y; PG_RETURN_POINT_P(result); } /* 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); } /*---------------------------------------------------------- * Conversion operators. *---------------------------------------------------------*/ Datum cr_circle(PG_FUNCTION_ARGS) { @@ -4940,65 +4955,65 @@ cr_circle(PG_FUNCTION_ARGS) result->radius = radius; PG_RETURN_CIRCLE_P(result); } Datum 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); } /* box_circle() * Convert a box to a circle. */ Datum box_circle(PG_FUNCTION_ARGS) { BOX *box = PG_GETARG_BOX_P(0); CIRCLE *circle; 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); PG_RETURN_CIRCLE_P(circle); } Datum circle_poly(PG_FUNCTION_ARGS) { int32 npts = PG_GETARG_INT32(0); CIRCLE *circle = PG_GETARG_CIRCLE_P(1); POLYGON *poly; int base_size, size; int i; - double angle; - double anglestep; + float8 angle; + float8 anglestep; if (FPzero(circle->radius)) ereport(ERROR, (errcode(ERRCODE_FEATURE_NOT_SUPPORTED), errmsg("cannot convert circle with radius zero to polygon"))); if (npts < 2) ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("must request at least 2 points"))); @@ -5009,27 +5024,30 @@ circle_poly(PG_FUNCTION_ARGS) /* Check for integer overflow */ if (base_size / npts != sizeof(poly->p[0]) || size <= base_size) ereport(ERROR, (errcode(ERRCODE_PROGRAM_LIMIT_EXCEEDED), errmsg("too many points requested"))); poly = (POLYGON *) palloc0(size); /* zero any holes */ 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); PG_RETURN_POLYGON_P(poly); } /* * Convert polygon to circle * @@ -5044,26 +5062,27 @@ poly_to_circle(CIRCLE *result, POLYGON *poly) int i; Assert(poly->npts > 0); result->center.x = 0; result->center.y = 0; result->radius = 0; 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 poly_circle(PG_FUNCTION_ARGS) { POLYGON *poly = PG_GETARG_POLYGON_P(0); CIRCLE *result; result = (CIRCLE *) palloc(sizeof(CIRCLE)); @@ -5088,44 +5107,44 @@ poly_circle(PG_FUNCTION_ARGS) * http://hopf.math.northwestern.edu/index.html * Description of algorithm: http://www.linuxjournal.com/article/2197 * http://www.linuxjournal.com/article/2029 */ #define POINT_ON_POLYGON INT_MAX 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; 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; /* loop over polygon points and aggregate total_cross */ 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) return 2; total_cross += cross; prev_x = x; prev_y = y; } @@ -5143,69 +5162,74 @@ point_inside(Point *p, int npts, Point *plist) /* lseg_crossing() * Returns +/-2 if line segment crosses the positive X-axis in a +/- direction. * Returns +/-1 if one point is on the positive X-axis. * Returns 0 if both points are on the positive X-axis, or there is no crossing. * Returns POINT_ON_POLYGON if the segment contains (0,0). * Wow, that is one confusing API, but it is used above, and when summed, * can tell is if a point is in a polygon. */ 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)) { /* y == 0, on X axis */ if (FPzero(x)) /* (x,y) is (0,0)? */ return POINT_ON_POLYGON; else if (FPgt(x, 0)) { /* 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; } } } static bool plist_same(int npts, Point *p1, Point *p2) { int i, ii, @@ -5275,47 +5299,53 @@ plist_same(int npts, Point *p1, Point *p2) * = 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. *----------------------------------------------------------------------- */ -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)) return get_float8_infinity(); if (isnan(x) || isnan(y)) return get_float8_nan(); /* Else, drop any minus signs */ x = fabs(x); y = fabs(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; } /* * If y is zero, the hypotenuse is x. This test saves a few cycles in * such cases, but more importantly it also protects against * divide-by-zero errors, since now x >= y. */ if (y == 0.0) return x; /* 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); + Assert(result >= 0.0); + + return result; } diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c index 096e59f817..b4008d7423 100644 --- a/src/backend/utils/adt/geo_spgist.c +++ b/src/backend/utils/adt/geo_spgist.c @@ -76,38 +76,38 @@ #include "access/spgist.h" #include "access/stratnum.h" #include "catalog/pg_type.h" #include "utils/fmgrprotos.h" #include "utils/geo_decls.h" /* * 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; return (x > y) ? 1 : -1; } typedef struct { - double low; - double high; + float8 low; + float8 high; } Range; typedef struct { Range left; Range right; } RangeBox; typedef struct { @@ -167,21 +167,21 @@ getRangeBox(BOX *box) /* * Initialize the traversal value * * In the beginning, we don't have any restrictions. We have to * initialize the struct to cover the whole 4D space. */ 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; rect_box->range_box_x.right.low = -infinity; rect_box->range_box_x.right.high = infinity; rect_box->range_box_y.left.low = -infinity; rect_box->range_box_y.left.high = infinity; @@ -410,40 +410,40 @@ spg_box_quad_choose(PG_FUNCTION_ARGS) * point as the median of the coordinates of the boxes. */ Datum spg_box_quad_picksplit(PG_FUNCTION_ARGS) { spgPickSplitIn *in = (spgPickSplitIn *) PG_GETARG_POINTER(0); spgPickSplitOut *out = (spgPickSplitOut *) PG_GETARG_POINTER(1); 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++) { BOX *box = DatumGetBoxP(in->datums[i]); lowXs[i] = box->low.x; highXs[i] = box->high.x; lowYs[i] = box->low.y; 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; centroid = palloc(sizeof(BOX)); centroid->low.x = lowXs[median]; centroid->high.x = highXs[median]; centroid->low.y = lowYs[median]; centroid->high.y = highYs[median]; diff --git a/src/include/utils/geo_decls.h b/src/include/utils/geo_decls.h index aeb31515e4..d14e8c8f72 100644 --- a/src/include/utils/geo_decls.h +++ b/src/include/utils/geo_decls.h @@ -1,42 +1,41 @@ /*------------------------------------------------------------------------- * * geo_decls.h - Declarations for various 2D constructs. * * * Portions Copyright (c) 1996-2018, PostgreSQL Global Development Group * Portions Copyright (c) 1994, Regents of the University of California * * 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 * and data types. There are still some more to do. - tgl 97/04/19 * *------------------------------------------------------------------------- */ #ifndef GEO_DECLS_H #define GEO_DECLS_H #include #include "fmgr.h" #include "utils/float.h" /*-------------------------------------------------------------------- * Useful floating point utilities and constants. - *-------------------------------------------------------------------*/ - + *------------------------------------------------------------------- + * + * XXX They are not NaN-aware. + */ #define EPSILON 1.0E-06 #ifdef EPSILON #define FPzero(A) (fabs(A) <= EPSILON) #define FPeq(A,B) (fabs((A) - (B)) <= EPSILON) #define FPne(A,B) (fabs((A) - (B)) > EPSILON) #define FPlt(A,B) ((B) - (A) > EPSILON) #define FPle(A,B) ((A) - (B) <= EPSILON) #define FPgt(A,B) ((A) - (B) > EPSILON) @@ -51,21 +50,21 @@ #define FPge(A,B) ((A) >= (B)) #endif #define HYPOT(A, B) pg_hypot(A, B) /*--------------------------------------------------------------------- * Point - (x,y) *-------------------------------------------------------------------*/ typedef struct { - double x, + float8 x, y; } Point; /*--------------------------------------------------------------------- * LSEG - A straight line, specified by endpoints. *-------------------------------------------------------------------*/ typedef struct { Point p[2]; @@ -83,21 +82,21 @@ typedef struct int32 dummy; /* padding to make it double align */ Point p[FLEXIBLE_ARRAY_MEMBER]; } PATH; /*--------------------------------------------------------------------- * LINE - Specified by its general equation (Ax+By+C=0). *-------------------------------------------------------------------*/ typedef struct { - double A, + float8 A, B, C; } LINE; /*--------------------------------------------------------------------- * BOX - Specified by two corner points, which are * sorted to save calculation time later. *-------------------------------------------------------------------*/ typedef struct @@ -118,21 +117,21 @@ typedef struct BOX boundbox; Point p[FLEXIBLE_ARRAY_MEMBER]; } POLYGON; /*--------------------------------------------------------------------- * CIRCLE - Specified by a center point and radius. *-------------------------------------------------------------------*/ typedef struct { Point center; - double radius; + float8 radius; } CIRCLE; /* * fmgr interface macros * * Path and Polygon are toastable varlena types, the others are just * fixed-size pass-by-reference types. */ #define DatumGetPointP(X) ((Point *) DatumGetPointer(X)) @@ -172,13 +171,13 @@ typedef struct #define DatumGetCircleP(X) ((CIRCLE *) DatumGetPointer(X)) #define CirclePGetDatum(X) PointerGetDatum(X) #define PG_GETARG_CIRCLE_P(n) DatumGetCircleP(PG_GETARG_DATUM(n)) #define PG_RETURN_CIRCLE_P(x) return CirclePGetDatum(x) /* * in geo_ops.c */ -extern double pg_hypot(double x, double y); +extern float8 pg_hypot(float8 x, float8 y); #endif /* GEO_DECLS_H */ -- 2.14.3 (Apple Git-98)