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 fmod() 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(x, NaN) = NaN
  33 //      Remainder(NaN, y) = NaN
  34 //      Remainder(Inf, y) = NaN
  35 //      Remainder(x, 0) = NaN
  36 //      Remainder(x, Inf) = x
  37 func Remainder(x, y float64) float64 {
  38         const (
  39                 Tiny    = 4.45014771701440276618e-308 // 0x0020000000000000
  40                 HalfMax = MaxFloat64 / 2
  41         )
  42         // TODO(rsc): Remove manual inlining of IsNaN, IsInf
  43         // when compiler does it for us
  44         // special cases
  45         switch {
  46         case x != x || y != y || x < -MaxFloat64 || x > MaxFloat64 || y == 0: // IsNaN(x) || IsNaN(y) || IsInf(x, 0) || y == 0:
  47                 return NaN()
  48         case y < -MaxFloat64 || y > MaxFloat64: // IsInf(y):
  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 = Fmod(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 }