From ca7b0db14e8f50bf81c601c66c1b785a977a88c7 Mon Sep 17 00:00:00 2001 From: Emre Hasegeli Date: Sat, 27 May 2017 16:17:42 -0400 Subject: [PATCH 3/4] geo-float-v05 Use the built-in float datatype to implement geometric types This will provide: * Check for underflow and overflow * Check for division by zero * Handle NaNs consistently 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 | 156 +++++------ src/backend/utils/adt/geo_ops.c | 531 ++++++++++++++++++++----------------- src/backend/utils/adt/geo_spgist.c | 36 +-- src/include/utils/geo_decls.h | 41 ++- 4 files changed, 400 insertions(+), 364 deletions(-) diff --git a/src/backend/access/gist/gistproc.c b/src/backend/access/gist/gistproc.c index 13524b4da0..6fa97602c5 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,52 +364,54 @@ 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) + if (float4_gt(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. */ - if (overlap < context->overlap || - (overlap == context->overlap && ratio > context->ratio)) + if (float4_lt(overlap, context->overlap) || + (float4_eq(overlap, context->overlap) && + float4_gt(ratio, context->ratio))) selectthis = true; } else { /* * Across dimensions, choose the new split if it has a smaller * *non-negative* overlap, or same *non-negative* overlap but * bigger range. This condition differs from the one described in * the article. On the datasets where leaf MBRs don't overlap * themselves, non-overlapping splits (i.e. splits which have zero @@ -416,55 +420,49 @@ g_box_consider_split(ConsiderSplitContext *context, int dimNum, * non-overlapping splits (i.e. having lowest negative overlap) * appears to be in the same dimension as in the previous split. * Therefore MBRs appear to be very prolonged along another * dimension, which leads to bad search performance. Using range * as the second split criteria makes MBRs more quadratic. Using * *non-negative* overlap instead of overlap as the first split * criteria gives to range criteria a chance to matter, because * non-overlapping splits are equivalent in this criteria. */ if (non_negative(overlap) < non_negative(context->overlap) || - (range > context->range && + (float4_gt(range, context->range) && non_negative(overlap) <= non_negative(context->overlap))) selectthis = true; } if (selectthis) { /* save information about selected split */ context->first = false; context->ratio = ratio; context->range = range; context->overlap = overlap; context->rightLower = rightLower; context->leftUpper = leftUpper; context->dim = 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); } /* * -------------------------------------------------------------------------- * Double sorting split algorithm. This is used for both boxes and points. * * The algorithm finds split of boxes by considering splits along each axis. * Each entry is first projected as an interval on the X-axis, and different * ways to split the intervals into two groups are considered, trying to * minimize the overlap of the groups. Then the same is repeated for the @@ -524,21 +522,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 +719,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 +774,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 +1098,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 +1139,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,80 +1207,84 @@ 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) + else if (float8_le(point->x, box->high.x) && + float8_ge(point->x, box->low.x) && + float8_le(point->y, box->high.y) && + float8_ge(point->y, box->low.y)) { /* point inside the box */ result = 0.0; } - else if (point->x <= box->high.x && point->x >= box->low.x) + else if (float8_le(point->x, box->high.x) && + float8_ge(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; - else if (point->y < box->low.y) - result = box->low.y - point->y; + Assert(float8_le(box->low.y, box->high.y)); + if (float8_gt(point->y, box->high.y)) + result = float8_mi(point->y, box->high.y); + else if (float8_lt(point->y, box->low.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) + else if (float8_le(point->y, box->high.y) && + float8_ge(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; - else if (point->x < box->low.x) - result = box->low.x - point->x; + Assert(float8_lt(box->low.x, box->high.x)); + if (float8_gt(point->x, box->high.x)) + result = float8_mi(point->x, box->high.x); + else if (float8_lt(point->x, box->low.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) + if (float8_gt(result, subresult)) result = subresult; p.x = box->low.x; p.y = box->high.y; subresult = point_point_distance(point, &p); - if (result > subresult) + if (float8_gt(result, subresult)) result = subresult; p.x = box->high.x; p.y = box->low.y; subresult = point_point_distance(point, &p); - if (result > subresult) + if (float8_gt(result, subresult)) result = subresult; } return result; } static bool gist_point_consistent_internal(StrategyNumber strategy, bool isLeaf, BOX *key, Point *query) { @@ -1363,24 +1365,24 @@ gist_point_consistent(PG_FUNCTION_ARGS) * Instead we write a non-fuzzy overlap test. The same code * will also serve for leaf-page tests, since leaf keys have * high == low. */ BOX *query, *key; query = PG_GETARG_BOX_P(1); key = DatumGetBoxP(entry->key); - result = (key->high.x >= query->low.x && - key->low.x <= query->high.x && - key->high.y >= query->low.y && - key->low.y <= query->high.y); + result = (float8_ge(key->high.x, query->low.x) && + float8_le(key->low.x, query->high.x) && + float8_ge(key->high.y, query->low.y) && + float8_le(key->low.y, query->high.y)); *recheck = false; } break; case PolygonStrategyNumberGroup: { POLYGON *query = PG_GETARG_POLYGON_P(1); result = DatumGetBool(DirectFunctionCall5( gist_poly_consistent, PointerGetDatum(entry), @@ -1389,22 +1391,22 @@ gist_point_consistent(PG_FUNCTION_ARGS) 0, PointerGetDatum(recheck))); if (GIST_LEAF(entry) && result) { /* * We are on leaf page and quick check shows overlapping * of polygon's bounding box and point */ BOX *box = DatumGetBoxP(entry->key); - Assert(box->high.x == box->low.x - && box->high.y == box->low.y); + Assert(float8_eq(box->high.x, box->low.x) && + float8_eq(box->high.y, box->low.y)); result = DatumGetBool(DirectFunctionCall2( poly_contain_pt, PolygonPGetDatum(query), PointPGetDatum(&box->high))); *recheck = false; } } break; case CircleStrategyNumberGroup: { @@ -1418,22 +1420,22 @@ gist_point_consistent(PG_FUNCTION_ARGS) 0, PointerGetDatum(recheck))); if (GIST_LEAF(entry) && result) { /* * We are on leaf page and quick check shows overlapping * of polygon's bounding box and point */ BOX *box = DatumGetBoxP(entry->key); - Assert(box->high.x == box->low.x - && box->high.y == box->low.y); + Assert(float8_eq(box->high.x, box->low.x) && + float8_eq(box->high.y, box->low.y)); result = DatumGetBool(DirectFunctionCall2( circle_contain_pt, CirclePGetDatum(query), PointPGetDatum(&box->high))); *recheck = false; } } break; default: elog(ERROR, "unrecognized strategy number: %d", strategy); @@ -1442,21 +1444,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 +1473,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 +1507,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 2cae70e3df..a4f6c1992e 100644 --- a/src/backend/utils/adt/geo_ops.c +++ b/src/backend/utils/adt/geo_ops.c @@ -29,76 +29,77 @@ * Internal routines */ enum path_delim { PATH_NONE, PATH_OPEN, PATH_CLOSED }; /* Routines for float8 */ static inline float8 float8_slope(float8 x1, float8 x2, float8 y1, float8 y2); +static inline float8 float8_hypot(float8 x, float8 y); /* Routines for two-dimensional points */ static inline void point_construct(Point *result, float8 x, float8 y); static inline void point_add_internal(Point *result, Point *pt1, Point *pt2); static inline void point_sub_internal(Point *result, Point *pt1, Point *pt2); static inline void point_mul_internal(Point *result, Point *pt1, Point *pt2); static inline void point_div_internal(Point *result, Point *pt1, Point *pt2); static inline bool point_eq_internal(Point *pt1, Point *pt2); static float8 point_dt(Point *pt1, Point *pt2); static int point_inside(Point *p, int npts, Point *plist); /* 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); /* Routines for two-dimensional circles */ -static double circle_ar(CIRCLE *circle); +static float8 circle_ar(CIRCLE *circle); /* Routines for two-dimensional lines */ 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_internal(Point *result, LINE *l1, LINE *l2); static inline float8 line_calculate_point(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_internal(Point *result, LSEG *lseg, LINE *line); static bool lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2); static bool lseg_inside_poly(Point *a, Point *b, POLYGON *poly, int start); static float8 lseg_dt(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 on_ps_internal(Point *pt, LSEG *lseg); /* Routines for two-dimensional polygons */ static void make_bound_box(POLYGON *poly); /* Unsorted */ static bool plist_same(int npts, Point *p1, Point *p2); -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); -static double dist_pl_internal(Point *pt, LINE *line); -static double dist_ps_internal(Point *pt, LSEG *lseg); -static double dist_ppoly_internal(Point *pt, POLYGON *poly); +static float8 dist_pl_internal(Point *pt, LINE *line); +static float8 dist_ps_internal(Point *pt, LSEG *lseg); +static float8 dist_ppoly_internal(Point *pt, POLYGON *poly); /* * Delimiters for input and output strings. * LDELIM, RDELIM, and DELIM are left, right, and separator delimiters, respectively. * LDELIM_EP, RDELIM_EP are left and right delimiters for paths with endpoints. */ #define LDELIM '(' #define RDELIM ')' @@ -127,38 +128,38 @@ static double dist_ppoly_internal(Point *pt, POLYGON *poly); * * 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); @@ -358,21 +359,42 @@ pair_count(char *s, char delim) * to calculate inverse slope. To achieve that pass the values in * (y1, y2, x2, x1) order. */ static inline float8 float8_slope(float8 x1, float8 x2, float8 y1, float8 y2) { if (FPeq(x1, x2)) return DBL_MAX; if (FPeq(y1, y2)) return 0.0; - return (y1 - y2) / (x1 - x2); + return float8_div(float8_mi(y1, y2), float8_mi(x1, x2)); +} + + +/* + * Wrapper around libc hypot() + */ +static inline float8 +float8_hypot(float8 x, float8 y) +{ + float8 result; + + errno = 0; + result = hypot(x, y); + if (errno == ERANGE) + ereport(ERROR, + (errcode(ERRCODE_NUMERIC_VALUE_OUT_OF_RANGE), + errmsg("value out of range: overflow"))); + check_float8_val(result, isinf(x) || isinf(y), x == 0.0 && y == 0.0); + Assert(result >= 0.0); + + return result; } /*********************************************************************** ** ** Routines for two-dimensional boxes. ** ***********************************************************************/ /*---------------------------------------------------------- @@ -384,33 +406,33 @@ float8_slope(float8 x1, float8 x2, float8 y1, float8 y2) * 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. @@ -424,38 +446,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); } /* @@ -474,31 +496,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; } } @@ -810,54 +832,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. @@ -867,24 +889,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 @@ -1015,30 +1037,30 @@ 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 { /* 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) { @@ -1066,21 +1088,22 @@ line_construct_pp(PG_FUNCTION_ARGS) /* * Calculate the line equation for a point * * This returns the result of the line equation Ax + By + C. The result * needs to be 0 for the point to be on the line. */ static inline float8 line_calculate_point(LINE *line, Point *pt) { - return line->A * pt->x + line->B * pt->y + line->C; + return float8_pl(float8_pl(float8_mul(line->A, pt->x), + float8_mul(line->B, pt->y)), line->C); } /*---------------------------------------------------------- * Relative position routines. *---------------------------------------------------------*/ Datum line_intersect(PG_FUNCTION_ARGS) { @@ -1103,21 +1126,22 @@ 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)); } @@ -1127,34 +1151,34 @@ 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; + ratio = float8_div(l1->A, l2->A); else if (!FPzero(l2->B)) - k = l1->B / l2->B; + ratio = float8_div(l1->B, l2->B); else if (!FPzero(l2->C)) - k = l1->C / 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))); } /*---------------------------------------------------------- * Line arithmetic routines. *---------------------------------------------------------*/ /* line_distance() * Distance between two lines. */ @@ -1162,21 +1186,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_internal(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 = dist_pl_internal(&tmp, l2); PG_RETURN_FLOAT8(result); } /* line_interpt() * Point where two lines l1, l2 intersect (if any) */ Datum line_interpt(PG_FUNCTION_ARGS) @@ -1199,41 +1223,41 @@ line_interpt(PG_FUNCTION_ARGS) * 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. */ static bool line_interpt_internal(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 (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; if (FPzero(l2->B)) x = l2->C; else - 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); #ifdef GEODEBUG printf("line_interpt- lines are A=%.*g, B=%.*g, C=%.*g, A=%.*g, B=%.*g, C=%.*g\n", DBL_DIG, l1->A, DBL_DIG, l1->B, DBL_DIG, l1->C, DBL_DIG, l2->A, DBL_DIG, l2->B, DBL_DIG, l2->C); printf("line_interpt- lines intersect at (%.*g,%.*g)\n", DBL_DIG, x, DBL_DIG, y); #endif @@ -1260,36 +1284,35 @@ line_interpt_internal(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; @@ -1506,31 +1529,31 @@ path_npoints(PG_FUNCTION_ARGS) PG_RETURN_INT32(path->npts); } Datum path_close(PG_FUNCTION_ARGS) { PATH *path = PG_GETARG_PATH_P_COPY(0); - path->closed = TRUE; + path->closed = true; PG_RETURN_PATH_P(path); } Datum path_open(PG_FUNCTION_ARGS) { PATH *path = PG_GETARG_PATH_P_COPY(0); - path->closed = FALSE; + path->closed = false; PG_RETURN_PATH_P(path); } /* path_inter - * Does p1 intersect p2 at any point? * Use bounding boxes for a quick (O(n)) check, then do a * O(n^2) iterative edge check. */ @@ -1546,33 +1569,33 @@ path_inter(PG_FUNCTION_ARGS) LSEG seg1, seg2; if (p1->npts <= 0 || p2->npts <= 0) PG_RETURN_BOOL(false); 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) @@ -1648,21 +1671,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_dt(&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(); @@ -1687,21 +1710,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. ** ***********************************************************************/ @@ -1875,21 +1898,21 @@ point_distance(PG_FUNCTION_ARGS) Point *pt2 = PG_GETARG_POINT_P(1); PG_RETURN_FLOAT8(point_dt(pt1, pt2)); } static float8 point_dt(Point *pt1, Point *pt2) { float8 result; - result = hypot(pt1->x - pt2->x, pt1->y - pt2->y); + result = float8_hypot(float8_mi(pt1->x, pt2->x), float8_mi(pt1->y, pt2->y)); #ifdef GEODEBUG printf("point_dt- segment (%f,%f),(%f,%f) length is %f\n", pt1->x, pt1->y, pt2->x, pt2->y, result); #endif return result; } Datum @@ -2049,35 +2072,35 @@ 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 float8_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 = float8_slope(l1->p[0].x, l1->p[1].x, l1->p[0].y, l1->p[1].y); m2 = float8_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)); else if (FPzero(m2)) PG_RETURN_BOOL(FPeq(m1, DBL_MAX)); - PG_RETURN_BOOL(FPeq(m1 / m2, -1.0)); + PG_RETURN_BOOL(FPeq(float8_div(m1, m2), -1.0)); } Datum lseg_vertical(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); PG_RETURN_BOOL(FPeq(lseg->p[0].x, lseg->p[1].x)); } @@ -2167,52 +2190,47 @@ lseg_distance(PG_FUNCTION_ARGS) LSEG *l2 = PG_GETARG_LSEG_P(1); PG_RETURN_FLOAT8(lseg_dt(l1, l2)); } /* lseg_dt() * Distance between two line segments. * Must check both sets of endpoints to ensure minimum distance is found. * - thomas 1998-02-01 */ -static double +static float8 lseg_dt(LSEG *l1, LSEG *l2) { - double result, - d; + float8 result; if (lseg_interpt_internal(NULL, l1, l2)) return 0.0; - d = dist_ps_internal(&l1->p[0], l2); - result = d; - d = dist_ps_internal(&l1->p[1], l2); - result = Min(result, d); - d = dist_ps_internal(&l2->p[0], l1); - result = Min(result, d); - d = dist_ps_internal(&l2->p[1], l1); - result = Min(result, d); + result = dist_ps_internal(&l1->p[0], l2); + result = float8_min(result, dist_ps_internal(&l1->p[1], l2)); + result = float8_min(result, dist_ps_internal(&l2->p[0], l1)); + result = float8_min(result, dist_ps_internal(&l2->p[1], l1)); return result; } 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); } static bool lseg_interpt_internal(Point *result, LSEG *l1, LSEG *l2) { Point interpt; LINE tmp1, tmp2; @@ -2296,45 +2314,44 @@ lseg_interpt(PG_FUNCTION_ARGS) */ Datum dist_pl(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LINE *line = PG_GETARG_LINE_P(1); PG_RETURN_FLOAT8(dist_pl_internal(pt, line)); } -static double +static float8 dist_pl_internal(Point *pt, LINE *line) { - return fabs(line_calculate_point(line, pt) / - hypot(line->A, line->B)); + return float8_div(fabs(line_calculate_point(line, pt)), + float8_hypot(line->A, line->B)); } /* * Distance from a point to a lseg */ Datum dist_ps(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LSEG *lseg = PG_GETARG_LSEG_P(1); PG_RETURN_FLOAT8(dist_ps_internal(pt, lseg)); } -static double +static float8 dist_ps_internal(Point *pt, LSEG *lseg) { - double m; /* slope of perp. */ - double result, - tmpdist; + float8 m; /* slope of perp. */ + float8 result; Point interpt; LINE ln; /* * Construct a line perpendicular to the input segment and through the * input point */ m = float8_slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x); line_construct_pm(&ln, pt, m); @@ -2354,24 +2371,22 @@ dist_ps_internal(Point *pt, LSEG *lseg) /* yes, so use distance to the intersection point */ result = point_dt(pt, &interpt); #ifdef GEODEBUG printf("dist_ps- distance is %f to intersection point is (%f,%f)\n", result, tmp.x, tmp.y); #endif } else { /* no, so use distance to the nearer endpoint */ - result = point_dt(pt, &lseg->p[0]); - tmpdist = point_dt(pt, &lseg->p[1]); - if (tmpdist < result) - result = tmpdist; + result = float8_min(point_dt(pt, &lseg->p[0]), + point_dt(pt, &lseg->p[1])); } return result; } /* * Distance from a point to a path */ Datum dist_ppath(PG_FUNCTION_ARGS) @@ -2409,21 +2424,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 = dist_ps_internal(pt, &lseg); - if (!have_min || tmp < result) + if (!have_min || float8_lt(tmp, result)) { result = tmp; have_min = true; } } break; } PG_RETURN_FLOAT8(result); } @@ -2447,33 +2462,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_internal(NULL, lseg, line)) result = 0.0; else - { - result = dist_pl_internal(&lseg->p[0], line); - d2 = dist_pl_internal(&lseg->p[1], line); /* XXX shouldn't we take the min not max? */ - if (d2 > result) - result = d2; - } + result = float8_max(dist_pl_internal(&lseg->p[0], line), + dist_pl_internal(&lseg->p[1], line)); PG_RETURN_FLOAT8(result); } /* * Distance from a lseg to a box */ Datum dist_sb(PG_FUNCTION_ARGS) { @@ -2515,25 +2525,22 @@ 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_min(float8_mi(dist_ppoly_internal(&circle->center, poly), + circle->radius), 0.0); PG_RETURN_FLOAT8(result); } /* * Distance from a point to a polygon */ Datum dist_ppoly(PG_FUNCTION_ARGS) { @@ -2551,21 +2558,21 @@ dist_polyp(PG_FUNCTION_ARGS) { POLYGON *poly = PG_GETARG_POLYGON_P(0); Point *point = PG_GETARG_POINT_P(1); float8 result; result = dist_ppoly_internal(point, poly); PG_RETURN_FLOAT8(result); } -static double +static float8 dist_ppoly_internal(Point *pt, POLYGON *poly) { float8 result; float8 d; int i; LSEG seg; if (point_inside(pt, poly->npts, poly->p) != 0) { #ifdef GEODEBUG @@ -2588,21 +2595,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 = dist_ps_internal(pt, &seg); #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. @@ -2672,21 +2679,21 @@ close_pl(PG_FUNCTION_ARGS) point_construct(result, pt->x, line->C); else { LINE tmp; /* * Drop a perpendicular and find the intersection point * * We need to invert the slope to get the perpendicular. */ - line_construct_pm(&tmp, pt, line->B / line->A); + line_construct_pm(&tmp, pt, float8_div(line->B, line->A)); line_interpt_internal(result, &tmp, line); } PG_RETURN_POINT_P(result); } /* close_ps() * Closest point on line segment to specified point. * Take the closest endpoint if the point is left, right, @@ -2696,21 +2703,21 @@ close_pl(PG_FUNCTION_ARGS) * Some tricky code here, relying on boolean expressions * evaluating to only zero or one to use as an array index. * bug fixes by gthaker@atl.lmco.com; May 1, 1998 */ Datum close_ps(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LSEG *lseg = PG_GETARG_LSEG_P(1); Point *result; - double invm; + float8 invm; int xh, yh; LINE tmp; result = (Point *) palloc(sizeof(Point)); #ifdef GEODEBUG printf("close_sp:pt->x %f pt->y %f\nlseg(0).x %f lseg(0).y %f lseg(1).x %f lseg(1).y %f\n", pt->x, pt->y, lseg->p[0].x, lseg->p[0].y, lseg->p[1].x, lseg->p[1].y); @@ -2755,32 +2762,32 @@ close_ps(PG_FUNCTION_ARGS) } /* * vert. and horiz. cases are down, now check if the closest point is one * of the end points or someplace on the lseg. */ invm = float8_slope(lseg->p[0].y, lseg->p[1].y, lseg->p[1].x, lseg->p[0].x); line_construct_pm(&tmp, &lseg->p[!yh], invm); /* lower edge of the * "band" */ - if (pt->y < (tmp.A * pt->x + tmp.C)) + if (pt->y < float8_pl(float8_mul(tmp.A, pt->x), tmp.C)) { /* we are below the lower edge */ *result = lseg->p[!yh]; /* below the lseg, take lower end pt */ #ifdef GEODEBUG printf("close_ps below: tmp A %f B %f C %f\n", tmp->A, tmp->B, tmp->C); #endif PG_RETURN_POINT_P(result); } line_construct_pm(&tmp, &lseg->p[yh], invm); /* upper edge of the * "band" */ - if (pt->y > (tmp.A * pt->x + tmp.C)) + if (pt->y > float8_pl(float8_mul(tmp.A, pt->x), tmp.C)) { /* we are below the lower edge */ *result = lseg->p[yh]; /* above the lseg, take higher end pt */ #ifdef GEODEBUG printf("close_ps above: tmp A %f B %f C %f\n", tmp->A, tmp->B, tmp->C); #endif PG_RETURN_POINT_P(result); } /* @@ -2816,32 +2823,32 @@ close_lseg(PG_FUNCTION_ARGS) { LSEG *l1 = PG_GETARG_LSEG_P(0); LSEG *l2 = PG_GETARG_LSEG_P(1); Point *result; double dist, dist0, dist1; dist0 = dist_ps_internal(&l1->p[0], l2); dist1 = dist_ps_internal(&l1->p[1], l2); - dist = Min(dist0, dist1); + dist = float8_min(dist0, dist1); - if (dist_ps_internal(&l2->p[0], l1) < dist) + if (float8_lt(dist_ps_internal(&l2->p[0], l1), dist)) { result = DatumGetPointP(DirectFunctionCall2(close_ps, PointPGetDatum(&l2->p[0]), LsegPGetDatum(l1))); result = DatumGetPointP(DirectFunctionCall2(close_ps, PointPGetDatum(result), LsegPGetDatum(l2))); } - else if (dist_ps_internal(&l2->p[1], l1) < dist) + else if (float8_lt(dist_ps_internal(&l2->p[1], l1), dist)) { result = DatumGetPointP(DirectFunctionCall2(close_ps, PointPGetDatum(&l2->p[1]), LsegPGetDatum(l1))); result = DatumGetPointP(DirectFunctionCall2(close_ps, PointPGetDatum(result), LsegPGetDatum(l2))); } else { @@ -2856,52 +2863,55 @@ close_lseg(PG_FUNCTION_ARGS) * Closest point on or in box to specified point. */ Datum close_pb(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); BOX *box = PG_GETARG_BOX_P(1); LSEG lseg, seg; Point point; - double dist, + float8 dist, d; if (DatumGetBool(DirectFunctionCall2(on_pb, PointPGetDatum(pt), BoxPGetDatum(box)))) PG_RETURN_POINT_P(pt); /* pairwise check lseg distances */ point.x = box->low.x; point.y = box->high.y; statlseg_construct(&lseg, &box->low, &point); dist = dist_ps_internal(pt, &lseg); statlseg_construct(&seg, &box->high, &point); - if ((d = dist_ps_internal(pt, &seg)) < dist) + d = dist_ps_internal(pt, &seg); + if (float8_lt(d, dist)) { dist = d; memcpy(&lseg, &seg, sizeof(lseg)); } point.x = box->high.x; point.y = box->low.y; statlseg_construct(&seg, &box->low, &point); - if ((d = dist_ps_internal(pt, &seg)) < dist) + d = dist_ps_internal(pt, &seg); + if (float8_lt(d, dist)) { dist = d; memcpy(&lseg, &seg, sizeof(lseg)); } statlseg_construct(&seg, &box->high, &point); - if ((d = dist_ps_internal(pt, &seg)) < dist) + d = dist_ps_internal(pt, &seg); + if (float8_lt(d, dist)) { dist = d; memcpy(&lseg, &seg, sizeof(lseg)); } PG_RETURN_DATUM(DirectFunctionCall2(close_ps, PointPGetDatum(pt), LsegPGetDatum(&lseg))); } @@ -2924,21 +2934,21 @@ close_sl(PG_FUNCTION_ARGS) float8 d1, d2; result = (Point *) palloc(sizeof(Point)); if (lseg_interpt_line_internal(result, lseg, line)) PG_RETURN_POINT_P(result); d1 = dist_pl_internal(&lseg->p[0], line); d2 = dist_pl_internal(&lseg->p[1], line); - 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"))); @@ -2958,77 +2968,80 @@ close_ls(PG_FUNCTION_ARGS) float8 d1, d2; result = (Point *) palloc(sizeof(Point)); if (lseg_interpt_line_internal(result, lseg, line)) PG_RETURN_POINT_P(result); d1 = dist_pl_internal(&lseg->p[0], line); d2 = dist_pl_internal(&lseg->p[1], line); - if (d1 < d2) + if (float8_lt(d1, d2)) *result = lseg->p[0]; else *result = lseg->p[1]; PG_RETURN_POINT_P(result); } /* close_sb() * Closest point on or in box to line segment. */ Datum close_sb(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); BOX *box = PG_GETARG_BOX_P(1); Point point; LSEG bseg, seg; - double dist, + float8 dist, d; /* segment intersects box? then just return closest point to center */ if (DatumGetBool(DirectFunctionCall2(inter_sb, LsegPGetDatum(lseg), BoxPGetDatum(box)))) { box_cn(&point, box); PG_RETURN_DATUM(DirectFunctionCall2(close_ps, PointPGetDatum(&point), LsegPGetDatum(lseg))); } /* pairwise check lseg distances */ point.x = box->low.x; point.y = box->high.y; statlseg_construct(&bseg, &box->low, &point); dist = lseg_dt(lseg, &bseg); statlseg_construct(&seg, &box->high, &point); - if ((d = lseg_dt(lseg, &seg)) < dist) + d = lseg_dt(lseg, &seg); + if (float8_lt(d, dist)) { dist = d; memcpy(&bseg, &seg, sizeof(bseg)); } point.x = box->high.x; point.y = box->low.y; statlseg_construct(&seg, &box->low, &point); - if ((d = lseg_dt(lseg, &seg)) < dist) + d = lseg_dt(lseg, &seg); + if (float8_lt(d, dist)) { dist = d; memcpy(&bseg, &seg, sizeof(bseg)); } statlseg_construct(&seg, &box->high, &point); - if ((d = lseg_dt(lseg, &seg)) < dist) + d = lseg_dt(lseg, &seg); + if (float8_lt(d, dist)) { dist = d; memcpy(&bseg, &seg, sizeof(bseg)); } /* OK, we now have the closest line segment on the box boundary */ PG_RETURN_DATUM(DirectFunctionCall2(close_lseg, LsegPGetDatum(lseg), LsegPGetDatum(&bseg))); } @@ -3076,50 +3089,55 @@ on_ps(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); LSEG *lseg = PG_GETARG_LSEG_P(1); PG_RETURN_BOOL(on_ps_internal(pt, lseg)); } static bool on_ps_internal(Point *pt, LSEG *lseg) { - return FPeq(point_dt(pt, &lseg->p[0]) + point_dt(pt, &lseg->p[1]), + return FPeq(float8_pl(point_dt(pt, &lseg->p[0]), + point_dt(pt, &lseg->p[1])), point_dt(&lseg->p[0], &lseg->p[1])); } /* * Check whether the point is in the box or on its border */ Datum on_pb(PG_FUNCTION_ARGS) { Point *pt = PG_GETARG_POINT_P(0); BOX *box = PG_GETARG_BOX_P(1); - PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x && - pt->y <= box->high.y && pt->y >= box->low.y); + PG_RETURN_BOOL(float8_le(pt->x, box->high.x) && + float8_ge(pt->x, box->low.x) && + float8_le(pt->y, box->high.y) && + float8_ge(pt->y, box->low.y)); } /* * Commutator of on_pb() */ Datum box_contain_pt(PG_FUNCTION_ARGS) { BOX *box = PG_GETARG_BOX_P(0); Point *pt = PG_GETARG_POINT_P(1); - PG_RETURN_BOOL(pt->x <= box->high.x && pt->x >= box->low.x && - pt->y <= box->high.y && pt->y >= box->low.y); + PG_RETURN_BOOL(float8_le(pt->x, box->high.x) && + float8_ge(pt->x, box->low.x) && + float8_le(pt->y, box->high.y) && + float8_ge(pt->y, box->low.y)); } /* on_ppath - * Whether a point lies within (on) a polyline. * If open, we have to (groan) check each segment. * (uses same algorithm as for point intersecting segment - tgl 1997-07-09) * If closed, we use the old O(n) ray method for point-in-polygon. * The ray is horizontal, from pt out to the right. * Each segment that crosses the ray counts as an @@ -3127,33 +3145,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); } @@ -3223,24 +3240,24 @@ inter_sl(PG_FUNCTION_ARGS) */ Datum inter_sb(PG_FUNCTION_ARGS) { LSEG *lseg = PG_GETARG_LSEG_P(0); BOX *box = PG_GETARG_BOX_P(1); 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)) PG_RETURN_BOOL(false); /* an endpoint of segment is inside box? then clearly intersects */ if (DatumGetBool(DirectFunctionCall2(on_pb, PointPGetDatum(&lseg->p[0]), BoxPGetDatum(box))) || DatumGetBool(DirectFunctionCall2(on_pb, @@ -3329,27 +3346,27 @@ make_bound_box(POLYGON *poly) { int i; BOX *boundbox = &poly->boundbox; if (poly->npts > 0) { boundbox->low.x = boundbox->high.x = poly->p[0].x; boundbox->low.y = boundbox->high.y = poly->p[0].y; for (i = 1; i < poly->npts; i++) { - if (poly->p[i].x < boundbox->low.x) + if (float8_lt(poly->p[i].x, boundbox->low.x)) boundbox->low.x = poly->p[i].x; - if (poly->p[i].x > boundbox->high.x) + if (float8_gt(poly->p[i].x, boundbox->high.x)) boundbox->high.x = poly->p[i].x; - if (poly->p[i].y < boundbox->low.y) + if (float8_lt(poly->p[i].y, boundbox->low.y)) boundbox->low.y = poly->p[i].y; - if (poly->p[i].y > boundbox->high.y) + if (float8_gt(poly->p[i].y, boundbox->high.y)) boundbox->high.y = poly->p[i].y; } } else ereport(ERROR, (errcode(ERRCODE_INVALID_PARAMETER_VALUE), errmsg("cannot create bounding box for empty polygon"))); } /*------------------------------------------------------------------ @@ -3475,21 +3492,21 @@ poly_send(PG_FUNCTION_ARGS) * the right most point of A left of the left most point * of B? *-------------------------------------------------------*/ Datum poly_left(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.high.x < polyb->boundbox.low.x; + result = float8_lt(polya->boundbox.high.x, polyb->boundbox.low.x); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3498,21 +3515,21 @@ poly_left(PG_FUNCTION_ARGS) * the right most point of A at or left of the right most point * of B? *-------------------------------------------------------*/ Datum poly_overleft(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.high.x <= polyb->boundbox.high.x; + result = float8_le(polya->boundbox.high.x, polyb->boundbox.high.x); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3521,21 +3538,21 @@ poly_overleft(PG_FUNCTION_ARGS) * the left most point of A right of the right most point * of B? *-------------------------------------------------------*/ Datum poly_right(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.low.x > polyb->boundbox.high.x; + result = float8_gt(polya->boundbox.low.x, polyb->boundbox.high.x); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3544,21 +3561,21 @@ poly_right(PG_FUNCTION_ARGS) * the left most point of A at or right of the left most point * of B? *-------------------------------------------------------*/ Datum poly_overright(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.low.x >= polyb->boundbox.low.x; + result = float8_ge(polya->boundbox.low.x, polyb->boundbox.low.x); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3567,21 +3584,21 @@ poly_overright(PG_FUNCTION_ARGS) * the upper most point of A below the lower most point * of B? *-------------------------------------------------------*/ Datum poly_below(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.high.y < polyb->boundbox.low.y; + result = float8_lt(polya->boundbox.high.y, polyb->boundbox.low.y); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3590,21 +3607,21 @@ poly_below(PG_FUNCTION_ARGS) * the upper most point of A at or below the upper most point * of B? *-------------------------------------------------------*/ Datum poly_overbelow(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.high.y <= polyb->boundbox.high.y; + result = float8_le(polya->boundbox.high.y, polyb->boundbox.high.y); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3613,21 +3630,21 @@ poly_overbelow(PG_FUNCTION_ARGS) * the lower most point of A above the upper most point * of B? *-------------------------------------------------------*/ Datum poly_above(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.low.y > polyb->boundbox.high.y; + result = float8_gt(polya->boundbox.low.y, polyb->boundbox.high.y); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3636,21 +3653,21 @@ poly_above(PG_FUNCTION_ARGS) * the lower most point of A at or above the lower most point * of B? *-------------------------------------------------------*/ Datum poly_overabove(PG_FUNCTION_ARGS) { POLYGON *polya = PG_GETARG_POLYGON_P(0); POLYGON *polyb = PG_GETARG_POLYGON_P(1); bool result; - result = polya->boundbox.low.y >= polyb->boundbox.low.y; + result = float8_ge(polya->boundbox.low.y, polyb->boundbox.low.y); /* * Avoid leaking memory for toasted inputs ... needed for rtree indexes */ PG_FREE_IF_COPY(polya, 0); PG_FREE_IF_COPY(polyb, 1); PG_RETURN_BOOL(result); } @@ -3847,22 +3864,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. *-----------------------------------------------------------------*/ @@ -3989,22 +4006,22 @@ point_add(PG_FUNCTION_ARGS) point_add_internal(result, p1, p2); PG_RETURN_POINT_P(result); } static inline void point_add_internal(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_sub(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; @@ -4012,22 +4029,22 @@ point_sub(PG_FUNCTION_ARGS) point_sub_internal(result, p1, p2); PG_RETURN_POINT_P(result); } static inline void point_sub_internal(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_mul(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; @@ -4035,22 +4052,24 @@ point_mul(PG_FUNCTION_ARGS) point_mul_internal(result, p1, p2); PG_RETURN_POINT_P(result); } static inline void point_mul_internal(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_div(PG_FUNCTION_ARGS) { Point *p1 = PG_GETARG_POINT_P(0); Point *p2 = PG_GETARG_POINT_P(1); Point *result; @@ -4059,24 +4078,26 @@ point_div(PG_FUNCTION_ARGS) point_div_internal(result, p1, p2); PG_RETURN_POINT_P(result); } static inline void point_div_internal(Point *result, Point *pt1, Point *pt2) { float8 div; - div = (pt2->x * pt2->x) + (pt2->y * pt2->y); + 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)); } /*********************************************************************** ** ** Routines for 2D boxes. ** ***********************************************************************/ Datum @@ -4185,24 +4206,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. ** ***********************************************************************/ @@ -4511,21 +4532,21 @@ 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) + if (float8_lt(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--; @@ -4578,21 +4599,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 (float8_lt(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 */ @@ -4629,144 +4650,146 @@ circle_same(PG_FUNCTION_ARGS) /* 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) + circle1->radius), circle2->radius)); + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + 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) + circle2->radius), circle1->radius)); + PG_RETURN_BOOL(FPle(point_dt(&circle1->center, &circle2->center), + 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); @@ -4865,36 +4888,38 @@ 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_internal(&result->center, &circle->center, point); - result->radius *= hypot(point->x, point->y); + result->radius = float8_mul(circle->radius, + float8_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_internal(&result->center, &circle->center, point); - result->radius /= hypot(point->x, point->y); + result->radius = float8_div(circle->radius, + float8_hypot(point->x, point->y)); PG_RETURN_CIRCLE_P(result); } /* circle_area - returns the area of the circle. */ Datum circle_area(PG_FUNCTION_ARGS) { @@ -4904,21 +4929,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); @@ -4929,81 +4954,84 @@ 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 (float8_lt(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); + PG_RETURN_BOOL(float8_le(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); + PG_RETURN_BOOL(float8_le(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 (float8_lt(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 (float8_lt(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; @@ -5011,24 +5039,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) { @@ -5043,65 +5071,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"))); @@ -5112,27 +5140,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); } /* poly_circle - convert polygon to circle * * XXX This algorithm should use weighted means of line segments @@ -5151,29 +5182,30 @@ poly_circle(PG_FUNCTION_ARGS) errmsg("cannot convert empty polygon to circle"))); circle = (CIRCLE *) palloc(sizeof(CIRCLE)); circle->center.x = 0; circle->center.y = 0; circle->radius = 0; for (i = 0; i < poly->npts; i++) { - circle->center.x += poly->p[i].x; - circle->center.y += poly->p[i].y; + circle->center.x = float8_pl(circle->center.x, poly->p[i].x); + circle->center.y = float8_pl(circle->center.y, poly->p[i].y); } - circle->center.x /= poly->npts; - circle->center.y /= poly->npts; + circle->center.x = float8_div(circle->center.x, poly->npts); + circle->center.y = float8_div(circle->center.y, poly->npts); for (i = 0; i < poly->npts; i++) - circle->radius += point_dt(&poly->p[i], &circle->center); - circle->radius /= poly->npts; + circle->radius = float8_pl(circle->radius, point_dt(&poly->p[i], + &circle->center)); + circle->radius = float8_div(circle->radius, poly->npts); PG_RETURN_CIRCLE_P(circle); } /*********************************************************************** ** ** Private routines for multiple types. ** ***********************************************************************/ @@ -5187,45 +5219,45 @@ 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; if (npts <= 0) return 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; } @@ -5243,69 +5275,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, diff --git a/src/backend/utils/adt/geo_spgist.c b/src/backend/utils/adt/geo_spgist.c index c800bb1338..f62be1061e 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 { @@ -121,30 +121,30 @@ typedef struct * The quadrant is 8 bit unsigned integer with 4 least bits in use. * This function accepts BOXes as input. They are not casted to * RangeBoxes, yet. All 4 bits are set by comparing a corner of the box. * This makes 16 quadrants in total. */ static uint8 getQuadrant(BOX *centroid, BOX *inBox) { uint8 quadrant = 0; - if (inBox->low.x > centroid->low.x) + if (float8_gt(inBox->low.x, centroid->low.x)) quadrant |= 0x8; - if (inBox->high.x > centroid->high.x) + if (float8_gt(inBox->high.x, centroid->high.x)) quadrant |= 0x4; - if (inBox->low.y > centroid->low.y) + if (float8_gt(inBox->low.y, centroid->low.y)) quadrant |= 0x2; - if (inBox->high.y > centroid->high.y) + if (float8_gt(inBox->high.y, centroid->high.y)) quadrant |= 0x1; return quadrant; } /* * Get RangeBox using BOX * * We are turning the BOX to our structures to emphasize their function * of representing points in 4D space. It also is more convenient to @@ -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 22c68f5b62..eee0ceff11 100644 --- a/src/include/utils/geo_decls.h +++ b/src/include/utils/geo_decls.h @@ -1,23 +1,20 @@ /*------------------------------------------------------------------------- * * geo_decls.h - Declarations for various 2D constructs. * * * Portions Copyright (c) 1996-2017, 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 @@ -28,43 +25,43 @@ /*-------------------------------------------------------------------- * Useful floating point utilities and constants. *-------------------------------------------------------------------*/ #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) -#define FPge(A,B) ((B) - (A) <= EPSILON) +#define FPeq(A, B) (float8_le(fabs(float8_mi(A, B)), EPSILON)) +#define FPne(A, B) (float8_gt(fabs(float8_mi(A, B)), EPSILON)) +#define FPlt(A, B) (float8_gt(float8_mi(B, A), EPSILON)) +#define FPle(A, B) (float8_le(float8_mi(A, B), EPSILON)) +#define FPgt(A, B) (float8_gt(float8_mi(A, B), EPSILON)) +#define FPge(A, B) (float8_le(float8_mi(B, A), EPSILON)) #else -#define FPzero(A) ((A) == 0) -#define FPeq(A,B) ((A) == (B)) -#define FPne(A,B) ((A) != (B)) -#define FPlt(A,B) ((A) < (B)) -#define FPle(A,B) ((A) <= (B)) -#define FPgt(A,B) ((A) > (B)) -#define FPge(A,B) ((A) >= (B)) +#define FPzero(A) ((A) == 0.0) +#define FPeq(A, B) (float8_eq(A, B)) +#define FPne(A, B) (float8_ne(A, B)) +#define FPlt(A, B) (float8_lt(A, B)) +#define FPle(A, B) (float8_le(A, B)) +#define FPgt(A, B) (float8_gt(A, B)) +#define FPge(A, B) (float8_ge(A, B)) #endif /*--------------------------------------------------------------------- * Point - (x,y) *-------------------------------------------------------------------*/ typedef struct { - double x, + float8 x, y; } Point; /*--------------------------------------------------------------------- * LSEG - A straight line, specified by endpoints. *-------------------------------------------------------------------*/ typedef struct { Point p[2]; @@ -72,66 +69,66 @@ typedef struct /*--------------------------------------------------------------------- * PATH - Specified by vertex points. *-------------------------------------------------------------------*/ typedef struct { int32 vl_len_; /* varlena header (do not touch directly!) */ int32 npts; int32 closed; /* is this a closed polygon? */ - int32 dummy; /* padding to make it double align */ + int32 dummy; /* padding to make it float8 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 { Point high, low; /* corner POINTs */ } BOX; /*--------------------------------------------------------------------- - * POLYGON - Specified by an array of doubles defining the points, + * POLYGON - Specified by an array of float8s defining the points, * keeping the number of points and the bounding box for * speed purposes. *-------------------------------------------------------------------*/ typedef struct { int32 vl_len_; /* varlena header (do not touch directly!) */ int32 npts; 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)) @@ -166,11 +163,11 @@ typedef struct #define PolygonPGetDatum(X) PointerGetDatum(X) #define PG_GETARG_POLYGON_P(n) DatumGetPolygonP(PG_GETARG_DATUM(n)) #define PG_GETARG_POLYGON_P_COPY(n) DatumGetPolygonPCopy(PG_GETARG_DATUM(n)) #define PG_RETURN_POLYGON_P(x) return PolygonPGetDatum(x) #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) -#endif /* GEO_DECLS_H */ +#endif /* GEO_DECLS_H */ -- 2.13.5 (Apple Git-94)