Subversion Repositories openrisc

// Copyright 2009 The Go Authors. All rights reserved.

// Use of this source code is governed by a BSD-style

// license that can be found in the LICENSE file.

package math

// Sqrt returns the square root of x.

//

// Special cases are:

//      Sqrt(+Inf) = +Inf

//      Sqrt(±0) = ±0

//      Sqrt(x < 0) = NaN

//      Sqrt(NaN) = NaN

//extern sqrt

func libc_sqrt(float64) float64

func Sqrt(x float64) float64 {

        return libc_sqrt(x)

}

// The original C code and the long comment below are

// from FreeBSD's /usr/src/lib/msun/src/e_sqrt.c and

// came with this notice.  The go code is a simplified

// version of the original C.

//

// ====================================================

// Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.

//

// Developed at SunPro, a Sun Microsystems, Inc. business.

// Permission to use, copy, modify, and distribute this

// software is freely granted, provided that this notice

// is preserved.

// ====================================================

//

// __ieee754_sqrt(x)

// Return correctly rounded sqrt.

//           -----------------------------------------

//           | Use the hardware sqrt if you have one |

//           -----------------------------------------

// Method:

//   Bit by bit method using integer arithmetic. (Slow, but portable)

//   1. Normalization

//      Scale x to y in [1,4) with even powers of 2:

//      find an integer k such that  1 <= (y=x*2**(2k)) < 4, then

//              sqrt(x) = 2**k * sqrt(y)

//   2. Bit by bit computation

//      Let q  = sqrt(y) truncated to i bit after binary point (q = 1),

//           i                                                   0

//                                     i+1         2

//          s  = 2*q , and      y  =  2   * ( y - q  ).          (1)

//           i      i            i                 i

//

//      To compute q    from q , one checks whether

//                  i+1       i

//

//                            -(i+1) 2

//                      (q + 2      )  <= y.                     (2)

//                        i

//                                                            -(i+1)

//      If (2) is false, then q   = q ; otherwise q   = q  + 2      .

//                             i+1   i             i+1   i

//

//      With some algebraic manipulation, it is not difficult to see

//      that (2) is equivalent to

//                             -(i+1)

//                      s  +  2       <= y                       (3)

//                       i                i

//

//      The advantage of (3) is that s  and y  can be computed by

//                                    i      i

//      the following recurrence formula:

//          if (3) is false

//

//          s     =  s  ,       y    = y   ;                     (4)

//           i+1      i          i+1    i

//

//      otherwise,

//                         -i                      -(i+1)

//          s     =  s  + 2  ,  y    = y  -  s  - 2              (5)

//           i+1      i          i+1    i     i

//

//      One may easily use induction to prove (4) and (5).

//      Note. Since the left hand side of (3) contain only i+2 bits,

//            it does not necessary to do a full (53-bit) comparison

//            in (3).

//   3. Final rounding

//      After generating the 53 bits result, we compute one more bit.

//      Together with the remainder, we can decide whether the

//      result is exact, bigger than 1/2ulp, or less than 1/2ulp

//      (it will never equal to 1/2ulp).

//      The rounding mode can be detected by checking whether

//      huge + tiny is equal to huge, and whether huge - tiny is

//      equal to huge for some floating point number "huge" and "tiny".

//

//

// Notes:  Rounding mode detection omitted.  The constants "mask", "shift",

// and "bias" are found in src/pkg/math/bits.go

// Sqrt returns the square root of x.

//

// Special cases are:

//      Sqrt(+Inf) = +Inf

//      Sqrt(±0) = ±0

//      Sqrt(x < 0) = NaN

//      Sqrt(NaN) = NaN

func sqrt(x float64) float64 {

        // special cases

        switch {

        case x == 0 || IsNaN(x) || IsInf(x, 1):

                return x

        case x < 0:

                return NaN()

        }

        ix := Float64bits(x)

        // normalize x

        exp := int((ix >> shift) & mask)

        if exp == 0 { // subnormal x

                for ix&1<

                        ix <<= 1

                        exp--

                }

                exp++

        }

        exp -= bias // unbias exponent

        ix &^= mask << shift

        ix |= 1 << shift

        if exp&1 == 1 { // odd exp, double x to make it even

                ix <<= 1

        }

        exp >>= 1 // exp = exp/2, exponent of square root

        // generate sqrt(x) bit by bit

        ix <<= 1

        var q, s uint64               // q = sqrt(x)

        r := uint64(1 << (shift + 1)) // r = moving bit from MSB to LSB

        for r != 0 {

                t := s + r

                if t <= ix {

                        s = t + r

                        ix -= t

                        q += r

                }

                ix <<= 1

                r >>= 1

        }

        // final rounding

        if ix != 0 { // remainder, result not exact

                q += q & 1 // round according to extra bit

        }

        ix = q>>1 + uint64(exp-1+bias)<

        return Float64frombits(ix)

}

func sqrtC(f float64, r *float64) {

        *r = sqrt(f)

}