From 3b948e0f20751df4d5b2d4081b80cae35a4b0702 Mon Sep 17 00:00:00 2001 From: Erica Fischer Date: Mon, 18 Mar 2024 11:37:47 -0700 Subject: [PATCH] Use 128-bit arithmetic to avoid overflow --- clip.cpp | 12 +++++------- geometry.cpp | 44 +++++++++++++++++--------------------------- geometry.hpp | 2 +- main.cpp | 2 +- 4 files changed, 24 insertions(+), 36 deletions(-) diff --git a/clip.cpp b/clip.cpp index d1af5deb..627c6038 100644 --- a/clip.cpp +++ b/clip.cpp @@ -177,28 +177,26 @@ int clip(long long *x0, long long *y0, long long *x1, long long *y1, long long x } else { // failed both tests, so calculate the line segment to clip // from an outside point to an intersection with clip edge - long long x = *x0, y = *y0; + __int128 x = *x0, y = *y0; // At least one endpoint is outside the clip rectangle; pick it. int outcodeOut = outcode0 ? outcode0 : outcode1; // XXX truncating division - long long shift = 1LL << (GLOBAL_DETAIL - 32); - // Now find the intersection point; // use formulas y = y0 + slope * (x - x0), x = x0 + (1 / slope) * (y - y0) if (outcodeOut & TOP) { // point is above the clip rectangle - x = *x0 + ((*x1 - *x0) / shift) * ((ymax - *y0) / shift) / (*y1 - *y0) * shift * shift; + x = *x0 + ((__int128) (*x1 - *x0)) * ((ymax - *y0)) / (*y1 - *y0); y = ymax; } else if (outcodeOut & BOTTOM) { // point is below the clip rectangle - x = *x0 + ((*x1 - *x0) / shift) * ((ymin - *y0) / shift) / (*y1 - *y0) * shift * shift; + x = *x0 + ((__int128) (*x1 - *x0)) * ((ymin - *y0)) / (*y1 - *y0); y = ymin; } else if (outcodeOut & RIGHT) { // point is to the right of clip rectangle - y = *y0 + ((*y1 - *y0) / shift) * ((xmax - *x0) / shift) / (*x1 - *x0) * shift * shift; + y = *y0 + ((__int128) (*y1 - *y0)) * ((xmax - *x0)) / (*x1 - *x0); x = xmax; } else if (outcodeOut & LEFT) { // point is to the left of clip rectangle - y = *y0 + ((*y1 - *y0) / shift) * ((xmin - *x0) / shift) / (*x1 - *x0) * shift * shift; + y = *y0 + ((__int128) (*y1 - *y0)) * ((xmin - *x0)) / (*x1 - *x0); x = xmin; } diff --git a/geometry.cpp b/geometry.cpp index 3d155628..5a107b6f 100644 --- a/geometry.cpp +++ b/geometry.cpp @@ -300,21 +300,15 @@ bool point_within_tile(long long x, long long y, int z) { return x >= 0 && y >= 0 && x < area && y < area; } -double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) { - long long p2x = segB_x - segA_x; - long long p2y = segB_y - segA_y; +double distance_from_line(__int128 point_x, __int128 point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y) { + __int128 p2x = segB_x - segA_x; + __int128 p2y = segB_y - segA_y; // These calculations must be made in integers instead of floating point // to make them consistent between x86 and arm floating point implementations. - // - // In a 32-bit world, coordinates may be up to 34 bits, so their product is up to 68 bits, - // making their sum up to 69 bits. Downshift before multiplying to keep them in range. - // - // If the world is bigger than 32 bits, scale down to 32 bits. - long long shift = 1LL << (GLOBAL_DETAIL - 32); - double something = ((p2x / 4 / shift) * (p2x / 8 / shift) + (p2y / 4 / shift) * (p2y / 8 / shift)) * 32.0 * shift * shift; + double something = (p2x) * (p2x) + (p2y) * (p2y); // likewise - double u = (0 == something) ? 0 : ((point_x - segA_x) / 4 / shift * (p2x / 8 / shift) + (point_y - segA_y) / 4 / shift * (p2y / 8 / shift)) * 32.0 * shift * shift / (something); + double u = (0 == something) ? 0 : ((point_x - segA_x) * (p2x) + (point_y - segA_y) * (p2y)) / (something); if (u >= 1) { u = 1; @@ -674,8 +668,8 @@ drawvec fix_polygon(const drawvec &geom) { // calculate centroid // a + 1 < size() because point 0 is duplicated at the end - long long xtotal = 0; - long long ytotal = 0; + __int128 xtotal = 0; + __int128 ytotal = 0; long long count = 0; for (size_t a = 0; a + 1 < ring.size(); a++) { xtotal += ring[a].x; @@ -685,16 +679,13 @@ drawvec fix_polygon(const drawvec &geom) { xtotal /= count; ytotal /= count; - long long shift = 1LL << (GLOBAL_DETAIL - 32); - // figure out which point is furthest from the centroid - long long dist2 = 0; - long long furthest = 0; + __int128 dist2 = 0; + size_t furthest = 0; for (size_t a = 0; a + 1 < ring.size(); a++) { - // division by 16 because these are z0 coordinates and we need to avoid overflow - long long xd = (ring[a].x - xtotal) / 16 / shift; - long long yd = (ring[a].y - ytotal) / 16 / shift; - long long d2 = xd * xd + yd * yd; + __int128 xd = (ring[a].x - xtotal); + __int128 yd = (ring[a].y - ytotal); + __int128 d2 = xd * xd + yd * yd; if (d2 > dist2 || (d2 == dist2 && ring[a] < ring[furthest])) { dist2 = d2; furthest = a; @@ -704,13 +695,12 @@ drawvec fix_polygon(const drawvec &geom) { // then figure out which point is furthest from *that*, // which will hopefully be a good origin point since it should be // at a far edge of the shape. - long long dist2b = 0; - long long furthestb = 0; + __int128 dist2b = 0; + size_t furthestb = 0; for (size_t a = 0; a + 1 < ring.size(); a++) { - // division by 16 because these are z0 coordinates and we need to avoid overflow - long long xd = (ring[a].x - ring[furthest].x) / 16 / shift; - long long yd = (ring[a].y - ring[furthest].y) / 16 / shift; - long long d2 = xd * xd + yd * yd; + __int128 xd = (ring[a].x - ring[furthest].x); + __int128 yd = (ring[a].y - ring[furthest].y); + __int128 d2 = xd * xd + yd * yd; if (d2 > dist2b || (d2 == dist2b && ring[a] < ring[furthestb])) { dist2b = d2; furthestb = a; diff --git a/geometry.hpp b/geometry.hpp index 48c70646..4a864c41 100644 --- a/geometry.hpp +++ b/geometry.hpp @@ -98,7 +98,7 @@ drawvec clip_lines(drawvec &geom, long long x1, long long y1, long long x2, long drawvec clip_point(drawvec &geom, long long x1, long long y1, long long x2, long long y2); void visvalingam(drawvec &ls, size_t start, size_t end, double threshold, size_t retain); int pnpoly(const drawvec &vert, size_t start, size_t nvert, long long testx, long long testy); -double distance_from_line(long long point_x, long long point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y); +double distance_from_line(__int128 point_x, __int128 point_y, long long segA_x, long long segA_y, long long segB_x, long long segB_y); std::string overzoom(const mvt_tile &tile, int oz, int ox, int oy, int nz, int nx, int ny, int detail, int buffer, std::set const &keep, bool do_compress, diff --git a/main.cpp b/main.cpp index f06f4620..a7e329c9 100644 --- a/main.cpp +++ b/main.cpp @@ -2431,7 +2431,7 @@ std::pair read_input(std::vector &sources, char *fname, i double total_tile_count = 0; for (int i = 1; i <= maxzoom; i++) { - double tile_count = ceil(area_sum / ((1LL << (GLOBAL_DETAIL - i)) * (1LL << (GLOBAL_DETAIL - i)))); + double tile_count = ceil(area_sum / ((__int128) (1LL << (GLOBAL_DETAIL - i)) * (1LL << (GLOBAL_DETAIL - i)))); total_tile_count += tile_count; // 2M tiles is an arbitrary limit, chosen to make tiling jobs