libgo/go/math/sqrt.go

   1 // Copyright 2009 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 // Sqrt returns the square root of x.
   8 //
   9 // Special cases are:
  10 //      Sqrt(+Inf) = +Inf
  11 //      Sqrt(±0) = ±0
  12 //      Sqrt(x < 0) = NaN
  13 //      Sqrt(NaN) = NaN
  14 func libc_sqrt(float64) float64 __asm__("sqrt")
  15 func Sqrt(x float64) float64 {
  16         return libc_sqrt(x)
  17 }
  18
  19 // The original C code and the long comment below are
  20 // from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and
  21 // came with this notice.  The go code is a simplified
  22 // version of the original C.
  23 //
  24 // ====================================================
  25 // Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  26 //
  27 // Developed at SunPro, a Sun Microsystems, Inc. business.
  28 // Permission to use, copy, modify, and distribute this
  29 // software is freely granted, provided that this notice
  30 // is preserved.
  31 // ====================================================
  32 //
  33 // __ieee754_sqrt(x)
  34 // Return correctly rounded sqrt.
  35 //           -----------------------------------------
  36 //           | Use the hardware sqrt if you have one |
  37 //           -----------------------------------------
  38 // Method:
  39 //   Bit by bit method using integer arithmetic. (Slow, but portable)
  40 //   1. Normalization
  41 //      Scale x to y in [1,4) with even powers of 2:
  42 //      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then
  43 //              sqrt(x) = 2**k * sqrt(y)
  44 //   2. Bit by bit computation
  45 //      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),
  46 //           i                                                   0
  47 //                                     i+1         2
  48 //          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)
  49 //           i      i            i                 i
  50 //
  51 //      To compute q    from q , one checks whether
  52 //                  i+1       i
  53 //
  54 //                            -(i+1) 2
  55 //                      (q + 2      )  <= y.                     (2)
  56 //                        i
  57 //                                                            -(i+1)
  58 //      If (2) is false, then q   = q ; otherwise q   = q  + 2      .
  59 //                             i+1   i             i+1   i
  60 //
  61 //      With some algebraic manipulation, it is not difficult to see
  62 //      that (2) is equivalent to
  63 //                             -(i+1)
  64 //                      s  +  2       <= y                       (3)
  65 //                       i                i
  66 //
  67 //      The advantage of (3) is that s  and y  can be computed by
  68 //                                    i      i
  69 //      the following recurrence formula:
  70 //          if (3) is false
  71 //
  72 //          s     =  s  ,       y    = y   ;                     (4)
  73 //           i+1      i          i+1    i
  74 //
  75 //      otherwise,
  76 //                         -i                      -(i+1)
  77 //          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)
  78 //           i+1      i          i+1    i     i
  79 //
  80 //      One may easily use induction to prove (4) and (5).
  81 //      Note. Since the left hand side of (3) contain only i+2 bits,
  82 //            it does not necessary to do a full (53-bit) comparison
  83 //            in (3).
  84 //   3. Final rounding
  85 //      After generating the 53 bits result, we compute one more bit.
  86 //      Together with the remainder, we can decide whether the
  87 //      result is exact, bigger than 1/2ulp, or less than 1/2ulp
  88 //      (it will never equal to 1/2ulp).
  89 //      The rounding mode can be detected by checking whether
  90 //      huge + tiny is equal to huge, and whether huge - tiny is
  91 //      equal to huge for some floating point number "huge" and "tiny".
  92 //
  93 //
  94 // Notes:  Rounding mode detection omitted.  The constants "mask", "shift",
  95 // and "bias" are found in src/pkg/math/bits.go
  96
  97 // Sqrt returns the square root of x.
  98 //
  99 // Special cases are:
 100 //      Sqrt(+Inf) = +Inf
 101 //      Sqrt(±0) = ±0
 102 //      Sqrt(x < 0) = NaN
 103 //      Sqrt(NaN) = NaN
 104 func sqrt(x float64) float64 {
 105         // special cases
 106         // TODO(rsc): Remove manual inlining of IsNaN, IsInf
 107         // when compiler does it for us
 108         switch {
 109         case x == 0 || x != x || x > MaxFloat64: // x == 0 || IsNaN(x) || IsInf(x, 1):
 110                 return x
 111         case x < 0:
 112                 return NaN()
 113         }
 114         ix := Float64bits(x)
 115         // normalize x
 116         exp := int((ix >> shift) & mask)
 117         if exp == 0 { // subnormal x
 118                 for ix&1<<shift == 0 {
 119                         ix <<= 1
 120                         exp--
 121                 }
 122                 exp++
 123         }
 124         exp -= bias // unbias exponent
 125         ix &^= mask << shift
 126         ix |= 1 << shift
 127         if exp&1 == 1 { // odd exp, double x to make it even
 128                 ix <<= 1
 129         }
 130         exp >>= 1 // exp = exp/2, exponent of square root
 131         // generate sqrt(x) bit by bit
 132         ix <<= 1
 133         var q, s uint64               // q = sqrt(x)
 134         r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB
 135         for r != 0 {
 136                 t := s + r
 137                 if t <= ix {
 138                         s = t + r
 139                         ix -= t
 140                         q += r
 141                 }
 142                 ix <<= 1
 143                 r >>= 1
 144         }
 145         // final rounding
 146         if ix != 0 { // remainder, result not exact
 147                 q += q & 1 // round according to extra bit
 148         }
 149         ix = q>>1 + uint64(exp-1+bias)<<shift // significand + biased exponent
 150         return Float64frombits(ix)
 151 }
 152
 153 func sqrtC(f float64, r *float64) {
 154         *r = sqrt(f)
 155 }