libgo/go/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         return remainder(x, y)
  39 }
  40
  41 func 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 }