Merge pull request #213 from prashantkhoje/ppc64le-fma-failures

Fix floating point precision failures caused by FMA on ppc64le

xy/area_centroid.go

@@ -155,9 +155,9 @@ func (calc *AreaCentroidCalculator) addTriangle(p0, p1, p2 geom.Coord, isPositiv
 	}
 	centroid3(p0, p1, p2, calc.triangleCent3)
 	area2 := area2(p0, p1, p2)
-	calc.cg3[0] += sign * area2 * calc.triangleCent3[0]
-	calc.cg3[1] += sign * area2 * calc.triangleCent3[1]
-	calc.areasum2 += sign * area2
+	calc.cg3[0] += float64(sign * area2 * calc.triangleCent3[0]) // nolint:unconvert
+	calc.cg3[1] += float64(sign * area2 * calc.triangleCent3[1]) // nolint:unconvert
+	calc.areasum2 += float64(sign * area2)                       // nolint:unconvert
 }
 
 // Returns three times the centroid of the triangle p1-p2-p3.
@@ -170,7 +170,7 @@ func centroid3(p1, p2, p3, c geom.Coord) {
 // Returns twice the signed area of the triangle p1-p2-p3,
 // positive if a,b,c are oriented ccw, and negative if cw.
 func area2(p1, p2, p3 geom.Coord) float64 {
-	return (p2[0]-p1[0])*(p3[1]-p1[1]) - (p3[0]-p1[0])*(p2[1]-p1[1])
+	return float64((p2[0]-p1[0])*(p3[1]-p1[1])) - float64((p3[0]-p1[0])*(p2[1]-p1[1])) // nolint:unconvert
 }
 
 // Adds the linear segments defined by an array of coordinates
@@ -186,8 +186,8 @@ func (calc *AreaCentroidCalculator) addLinearSegments(pts []float64) {
 		calc.totalLength += segmentLen
 
 		midx := (pts[i] + pts[i+stride]) / 2
-		calc.centSum[0] += segmentLen * midx
+		calc.centSum[0] += float64(segmentLen * midx) // nolint:unconvert
 		midy := (pts[i+1] + pts[i+stride+1]) / 2
-		calc.centSum[1] += segmentLen * midy
+		calc.centSum[1] += float64(segmentLen * midy) // nolint:unconvert
 	}
 }

xy/cga.go

@@ -338,7 +338,7 @@ func SignedArea(layout geom.Layout, ring []float64) float64 {
 		x := ring[i] - x0
 		y1 := ring[i+stride+1]
 		y2 := ring[i-stride+1]
-		sum += x * (y2 - y1)
+		sum += float64(x * (y2 - y1)) // nolint:unconvert
 	}
 	return sum / 2.0
 }

xy/internal/cga.go

@@ -60,5 +60,5 @@ func Distance2D(c1, c2 geom.Coord) float64 {
 	dx := c1[0] - c2[0]
 	dy := c1[1] - c2[1]
 
-	return math.Sqrt(dx*dx + dy*dy)
+	return math.Sqrt(float64(dx*dx) + float64(dy*dy)) // nolint:unconvert
 }

xy/internal/hcoords/hcoords.go

@@ -17,15 +17,15 @@ func GetIntersection(line1End1, line1End2, line2End1, line2End2 geom.Coord) (geo
 	// unrolled computation
 	line1Xdiff := line1End1[1] - line1End2[1]
 	line1Ydiff := line1End2[0] - line1End1[0]
-	line1W := line1End1[0]*line1End2[1] - line1End2[0]*line1End1[1]
+	line1W := float64(line1End1[0]*line1End2[1]) - float64(line1End2[0]*line1End1[1]) // nolint:unconvert
 
 	line2X := line2End1[1] - line2End2[1]
 	line2Y := line2End2[0] - line2End1[0]
-	line2W := line2End1[0]*line2End2[1] - line2End2[0]*line2End1[1]
+	line2W := float64(line2End1[0]*line2End2[1]) - float64(line2End2[0]*line2End1[1]) // nolint:unconvert
 
-	x := line1Ydiff*line2W - line2Y*line1W
-	y := line2X*line1W - line1Xdiff*line2W
-	w := line1Xdiff*line2Y - line2X*line1Ydiff
+	x := float64(line1Ydiff*line2W) - float64(line2Y*line1W)     // nolint:unconvert
+	y := float64(line2X*line1W) - float64(line1Xdiff*line2W)     // nolint:unconvert
+	w := float64(line1Xdiff*line2Y) - float64(line2X*line1Ydiff) // nolint:unconvert
 
 	xIntersection := x / w
 	yIntersection := y / w

xy/line_centroid.go

@@ -107,8 +107,8 @@ func (calc *LineCentroidCalculator) addLine(line []float64, startLine, endLine i
 		calc.totalLength += segmentLen
 
 		midx := (line[i] + line[i+calc.stride]) / 2
-		calc.centSum[0] += segmentLen * midx
+		calc.centSum[0] += float64(segmentLen * midx) // nolint:unconvert
 		midy := (line[i+1] + line[i+calc.stride+1]) / 2
-		calc.centSum[1] += segmentLen * midy
+		calc.centSum[1] += float64(segmentLen * midy) // nolint:unconvert
 	}
 }

xyz/vector.go

@@ -14,7 +14,7 @@ func VectorDot(v1Start, v1End, v2Start, v2End geom.Coord) float64 {
 	v2Startv2Endx := v2End[0] - v2Start[0]
 	v2Startv2Endy := v2End[1] - v2Start[1]
 	v2Startv2Endz := v2End[2] - v2Start[2]
-	return v1Startv2Endx*v2Startv2Endx + v1Startv2Endy*v2Startv2Endy + v1Startv2Endz*v2Startv2Endz
+	return float64(v1Startv2Endx*v2Startv2Endx) + float64(v1Startv2Endy*v2Startv2Endy) + float64(v1Startv2Endz*v2Startv2Endz) // nolint:unconvert
 }
 
 // VectorNormalize creates a coordinate that is the normalized vector from 0,0,0 to vector
@@ -25,5 +25,5 @@ func VectorNormalize(vector geom.Coord) geom.Coord {
 
 // VectorLength calculates the length of the vector from 0,0,0 to vector
 func VectorLength(vector geom.Coord) float64 {
-	return math.Sqrt(vector[0]*vector[0] + vector[1]*vector[1] + vector[2]*vector[2])
+	return math.Sqrt(float64(vector[0]*vector[0]) + float64(vector[1]*vector[1]) + float64(vector[2]*vector[2])) // nolint:unconvert
 }

xyz/xyz.go

@@ -20,7 +20,7 @@ func Distance(point1, point2 geom.Coord) float64 {
 	dx := point1[0] - point2[0]
 	dy := point1[1] - point2[1]
 	dz := point1[2] - point2[2]
-	return math.Sqrt(dx*dx + dy*dy + dz*dz)
+	return math.Sqrt(float64(dx*dx) + float64(dy*dy) + float64(dz*dz)) // nolint:unconvert
 }
 
 // DistancePointToLine calculates the distance from point to a point on the line
@@ -44,11 +44,15 @@ func DistancePointToLine(point, lineStart, lineEnd geom.Coord) float64 {
 	 *   0<r<1 P is interior to AB
 	 */
 
-	len2 := (lineEnd[0]-lineStart[0])*(lineEnd[0]-lineStart[0]) + (lineEnd[1]-lineStart[1])*(lineEnd[1]-lineStart[1]) + (lineEnd[2]-lineStart[2])*(lineEnd[2]-lineStart[2])
+	len2 := float64((lineEnd[0]-lineStart[0])*(lineEnd[0]-lineStart[0])) + // nolint:unconvert
+		float64((lineEnd[1]-lineStart[1])*(lineEnd[1]-lineStart[1])) +     // nolint:unconvert
+		float64((lineEnd[2]-lineStart[2])*(lineEnd[2]-lineStart[2]))       // nolint:unconvert
 	if math.IsNaN(len2) {
 		panic("Ordinates must not be NaN")
 	}
-	r := ((point[0]-lineStart[0])*(lineEnd[0]-lineStart[0]) + (point[1]-lineStart[1])*(lineEnd[1]-lineStart[1]) + (point[2]-lineStart[2])*(lineEnd[2]-lineStart[2])) / len2
+	r := (float64((point[0]-lineStart[0])*(lineEnd[0]-lineStart[0])) +     // nolint:unconvert
+		float64((point[1]-lineStart[1])*(lineEnd[1]-lineStart[1])) +       // nolint:unconvert
+		float64((point[2]-lineStart[2])*(lineEnd[2]-lineStart[2]))) / len2 // nolint:unconvert
 
 	if r <= 0.0 {
 		return Distance(point, lineStart)
@@ -58,14 +62,14 @@ func DistancePointToLine(point, lineStart, lineEnd geom.Coord) float64 {
 	}
 
 	// compute closest point q on line segment
-	qx := lineStart[0] + r*(lineEnd[0]-lineStart[0])
-	qy := lineStart[1] + r*(lineEnd[1]-lineStart[1])
-	qz := lineStart[2] + r*(lineEnd[2]-lineStart[2])
+	qx := lineStart[0] + float64(r*(lineEnd[0]-lineStart[0])) // nolint:unconvert
+	qy := lineStart[1] + float64(r*(lineEnd[1]-lineStart[1])) // nolint:unconvert
+	qz := lineStart[2] + float64(r*(lineEnd[2]-lineStart[2])) // nolint:unconvert
 	// result is distance from p to q
 	dx := point[0] - qx
 	dy := point[1] - qy
 	dz := point[2] - qz
-	return math.Sqrt(dx*dx + dy*dy + dz*dz)
+	return math.Sqrt(float64(dx*dx) + float64(dy*dy) + float64(dz*dz)) // nolint:unconvert
 }
 
 // Equals determines if the two coordinates have equal in 3d space
@@ -98,7 +102,7 @@ func DistanceLineToLine(line1Start, line1End, line2Start, line2End geom.Coord) f
 	d := VectorDot(line1Start, line1End, line2Start, line1Start)
 	e := VectorDot(line2Start, line2End, line2Start, line1Start)
 
-	denom := a*c - b*b
+	denom := float64(a*c) - float64(b*b) // nolint:unconvert
 	if math.IsNaN(denom) {
 		panic("Ordinates must not be NaN")
 	}
@@ -117,8 +121,8 @@ func DistanceLineToLine(line1Start, line1End, line2Start, line2End geom.Coord) f
 			t = e / c
 		}
 	} else {
-		s = (b*e - c*d) / denom
-		t = (a*e - b*d) / denom
+		s = (float64(b*e) - float64(c*d)) / denom // nolint:unconvert
+		t = (float64(a*e) - float64(b*d)) / denom // nolint:unconvert
 	}
 	switch {
 	case s < 0:
@@ -134,13 +138,13 @@ func DistanceLineToLine(line1Start, line1End, line2Start, line2End geom.Coord) f
 	 * The closest points are in interiors of segments,
 	 * so compute them directly
 	 */
-	x1 := line1Start[0] + s*(line1End[0]-line1Start[0])
-	y1 := line1Start[1] + s*(line1End[1]-line1Start[1])
-	z1 := line1Start[2] + s*(line1End[2]-line1Start[2])
+	x1 := line1Start[0] + float64(s*(line1End[0]-line1Start[0])) // nolint:unconvert
+	y1 := line1Start[1] + float64(s*(line1End[1]-line1Start[1])) // nolint:unconvert
+	z1 := line1Start[2] + float64(s*(line1End[2]-line1Start[2])) // nolint:unconvert
 
-	x2 := line2Start[0] + t*(line2End[0]-line2Start[0])
-	y2 := line2Start[1] + t*(line2End[1]-line2Start[1])
-	z2 := line2Start[2] + t*(line2End[2]-line2Start[2])
+	x2 := line2Start[0] + float64(t*(line2End[0]-line2Start[0])) // nolint:unconvert
+	y2 := line2Start[1] + float64(t*(line2End[1]-line2Start[1])) // nolint:unconvert
+	z2 := line2Start[2] + float64(t*(line2End[2]-line2Start[2])) // nolint:unconvert
 
 	// length (p1-p2)
 	return Distance(geom.Coord{x1, y1, z1}, geom.Coord{x2, y2, z2})