Source file src/pkg/math/remainder.go

     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 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	}