1 files changed, 94 insertions, 0 deletions
diff --git a/src/math/remainder.go b/src/math/remainder.go
new file mode 100644
index 0000000..bf8bfd5
--- /dev/null
+++ b/src/math/remainder.go
@@ -0,0 +1,94 @@
+// Copyright 2010 The Go Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package math
+
+// The original C code and the comment below are from
+// FreeBSD's /usr/src/lib/msun/src/e_remainder.c and came
+// with this notice. The go code is a simplified version of
+// the original C.
+//
+// ====================================================
+// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
+//
+// Developed at SunPro, a Sun Microsystems, Inc. business.
+// Permission to use, copy, modify, and distribute this
+// software is freely granted, provided that this notice
+// is preserved.
+// ====================================================
+//
+// __ieee754_remainder(x,y)
+// Return :
+//      returns  x REM y  =  x - [x/y]*y  as if in infinite
+//      precision arithmetic, where [x/y] is the (infinite bit)
+//      integer nearest x/y (in half way cases, choose the even one).
+// Method :
+//      Based on Mod() returning  x - [x/y]chopped * y  exactly.
+
+// Remainder returns the IEEE 754 floating-point remainder of x/y.
+//
+// Special cases are:
+//	Remainder(±Inf, y) = NaN
+//	Remainder(NaN, y) = NaN
+//	Remainder(x, 0) = NaN
+//	Remainder(x, ±Inf) = x
+//	Remainder(x, NaN) = NaN
+func Remainder(x, y float64) float64 {
+	if haveArchRemainder {
+		return archRemainder(x, y)
+	}
+	return remainder(x, y)
+}
+
+func remainder(x, y float64) float64 {
+	const (
+		Tiny    = 4.45014771701440276618e-308 // 0x0020000000000000
+		HalfMax = MaxFloat64 / 2
+	)
+	// special cases
+	switch {
+	case IsNaN(x) || IsNaN(y) || IsInf(x, 0) || y == 0:
+		return NaN()
+	case IsInf(y, 0):
+		return x
+	}
+	sign := false
+	if x < 0 {
+		x = -x
+		sign = true
+	}
+	if y < 0 {
+		y = -y
+	}
+	if x == y {
+		if sign {
+			zero := 0.0
+			return -zero
+		}
+		return 0
+	}
+	if y <= HalfMax {
+		x = Mod(x, y+y) // now x < 2y
+	}
+	if y < Tiny {
+		if x+x > y {
+			x -= y
+			if x+x >= y {
+				x -= y
+			}
+		}
+	} else {
+		yHalf := 0.5 * y
+		if x > yHalf {
+			x -= y
+			if x >= yHalf {
+				x -= y
+			}
+		}
+	}
+	if sign {
+		x = -x
+	}
+	return x
+}