gitee.com/quant1x/num@v0.3.2/math32/remainder.go (about)

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