Home | History | Annotate | Download | only in math 
      1 // Copyright 2010 The Go Authors. All rights reserved.
      2 // Use of this source code is governed by a BSD-style
      3 // license that can be found in the LICENSE file.
      4 
      5 package math
      6 
      7 // The original C code and the comment below are from
      8 // FreeBSD's /usr/src/lib/msun/src/e_remainder.c and came
      9 // with this notice.  The go code is a simplified version of
     10 // the original C.
     11 //
     12 // ====================================================
     13 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
     14 //
     15 // Developed at SunPro, a Sun Microsystems, Inc. business.
     16 // Permission to use, copy, modify, and distribute this
     17 // software is freely granted, provided that this notice
     18 // is preserved.
     19 // ====================================================
     20 //
     21 // __ieee754_remainder(x,y)
     22 // Return :
     23 //      returns  x REM y  =  x - [x/y]*y  as if in infinite
     24 //      precision arithmetic, where [x/y] is the (infinite bit)
     25 //      integer nearest x/y (in half way cases, choose the even one).
     26 // Method :
     27 //      Based on Mod() returning  x - [x/y]chopped * y  exactly.
     28 
     29 // Remainder returns the IEEE 754 floating-point remainder of x/y.
     30 //
     31 // Special cases are:
     32 //	Remainder(Inf, y) = NaN
     33 //	Remainder(NaN, y) = NaN
     34 //	Remainder(x, 0) = NaN
     35 //	Remainder(x, Inf) = x
     36 //	Remainder(x, NaN) = NaN
     37 func Remainder(x, y float64) float64
     38 
     39 func remainder(x, y float64) float64 {
     40 	const (
     41 		Tiny    = 4.45014771701440276618e-308 // 0x0020000000000000
     42 		HalfMax = MaxFloat64 / 2
     43 	)
     44 	// special cases
     45 	switch {
     46 	case IsNaN(x) || IsNaN(y) || IsInf(x, 0) || y == 0:
     47 		return NaN()
     48 	case IsInf(y, 0):
     49 		return x
     50 	}
     51 	sign := false
     52 	if x < 0 {
     53 		x = -x
     54 		sign = true
     55 	}
     56 	if y < 0 {
     57 		y = -y
     58 	}
     59 	if x == y {
     60 		return 0
     61 	}
     62 	if y <= HalfMax {
     63 		x = Mod(x, y+y) // now x < 2y
     64 	}
     65 	if y < Tiny {
     66 		if x+x > y {
     67 			x -= y
     68 			if x+x >= y {
     69 				x -= y
     70 			}
     71 		}
     72 	} else {
     73 		yHalf := 0.5 * y
     74 		if x > yHalf {
     75 			x -= y
     76 			if x >= yHalf {
     77 				x -= y
     78 			}
     79 		}
     80 	}
     81 	if sign {
     82 		x = -x
     83 	}
     84 	return x
     85 }
     86