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