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
5package 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
37func Remainder(x, y float64) float64 {
38	return remainder(x, y)
39}
40
41func remainder(x, y float64) float64 {
42	const (
43		Tiny    = 4.45014771701440276618e-308 // 0x0020000000000000
44		HalfMax = MaxFloat64 / 2
45	)
46	// special cases
47	switch {
48	case IsNaN(x) || IsNaN(y) || IsInf(x, 0) || y == 0:
49		return NaN()
50	case IsInf(y, 0):
51		return x
52	}
53	sign := false
54	if x < 0 {
55		x = -x
56		sign = true
57	}
58	if y < 0 {
59		y = -y
60	}
61	if x == y {
62		return 0
63	}
64	if y <= HalfMax {
65		x = Mod(x, y+y) // now x < 2y
66	}
67	if y < Tiny {
68		if x+x > y {
69			x -= y
70			if x+x >= y {
71				x -= y
72			}
73		}
74	} else {
75		yHalf := 0.5 * y
76		if x > yHalf {
77			x -= y
78			if x >= yHalf {
79				x -= y
80			}
81		}
82	}
83	if sign {
84		x = -x
85	}
86	return x
87}
88