sqrt.go 1.2 KB

1234567891011121314151617181920212223242526272829303132333435363738394041424344454647484950515253545556575859606162636465
  1. package math32
  2. // Sqrt returns the square root of x.
  3. // Special cases are:
  4. // Sqrt(+Inf) = +Inf
  5. // Sqrt(±0) = ±0
  6. // Sqrt(x < 0) = NaN
  7. // Sqrt(NaN) = NaN
  8. func Sqrt(x float32) float32 {
  9. if haveArchSqrt {
  10. return archSqrt(x)
  11. }
  12. return sqrt(x)
  13. }
  14. // TODO: add assembly for !build noasm
  15. func sqrt(x float32) float32 {
  16. // special cases
  17. switch {
  18. case x == 0 || IsNaN(x) || IsInf(x, 1):
  19. return x
  20. case x < 0:
  21. return NaN()
  22. }
  23. ix := Float32bits(x)
  24. // normalize x
  25. exp := int((ix >> shift) & mask)
  26. if exp == 0 { // subnormal x
  27. for ix&(1<<shift) == 0 {
  28. ix <<= 1
  29. exp--
  30. }
  31. exp++
  32. }
  33. exp -= bias // unbias exponent
  34. ix &^= mask << shift
  35. ix |= 1 << shift
  36. if exp&1 == 1 { // odd exp, double x to make it even
  37. ix <<= 1
  38. }
  39. exp >>= 1 // exp = exp/2, exponent of square root
  40. // generate sqrt(x) bit by bit
  41. ix <<= 1
  42. var q, s uint32 // q = sqrt(x)
  43. r := uint32(1 << (shift + 1)) // r = moving bit from MSB to LSB
  44. for r != 0 {
  45. t := s + r
  46. if t <= ix {
  47. s = t + r
  48. ix -= t
  49. q += r
  50. }
  51. ix <<= 1
  52. r >>= 1
  53. }
  54. // final rounding
  55. if ix != 0 { // remainder, result not exact
  56. q += q & 1 // round according to extra bit
  57. }
  58. ix = q>>1 + uint32(exp-1+bias)<<shift // significand + biased exponent
  59. return Float32frombits(ix)
  60. }