scalar.go 4.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192
  1. // Copyright ©2013 The Gonum Authors. All rights reserved.
  2. // Use of this code is governed by a BSD-style
  3. // license that can be found in the LICENSE file.
  4. package scalar
  5. import (
  6. "math"
  7. "strconv"
  8. )
  9. // EqualWithinAbs returns true when a and b have an absolute difference
  10. // not greater than tol.
  11. func EqualWithinAbs(a, b, tol float64) bool {
  12. return a == b || math.Abs(a-b) <= tol
  13. }
  14. // minNormalFloat64 is the smallest normal number. For 64 bit IEEE-754
  15. // floats this is 2^{-1022}.
  16. const minNormalFloat64 = 0x1p-1022
  17. // EqualWithinRel returns true when the difference between a and b
  18. // is not greater than tol times the greater absolute value of a and b,
  19. // abs(a-b) <= tol * max(abs(a), abs(b)).
  20. func EqualWithinRel(a, b, tol float64) bool {
  21. if a == b {
  22. return true
  23. }
  24. delta := math.Abs(a - b)
  25. if delta <= minNormalFloat64 {
  26. return delta <= tol*minNormalFloat64
  27. }
  28. // We depend on the division in this relationship to identify
  29. // infinities (we rely on the NaN to fail the test) otherwise
  30. // we compare Infs of the same sign and evaluate Infs as equal
  31. // independent of sign.
  32. return delta/math.Max(math.Abs(a), math.Abs(b)) <= tol
  33. }
  34. // EqualWithinAbsOrRel returns true when a and b are equal to within
  35. // the absolute or relative tolerances. See EqualWithinAbs and
  36. // EqualWithinRel for details.
  37. func EqualWithinAbsOrRel(a, b, absTol, relTol float64) bool {
  38. return EqualWithinAbs(a, b, absTol) || EqualWithinRel(a, b, relTol)
  39. }
  40. // EqualWithinULP returns true when a and b are equal to within
  41. // the specified number of floating point units in the last place.
  42. func EqualWithinULP(a, b float64, ulp uint) bool {
  43. if a == b {
  44. return true
  45. }
  46. if math.IsNaN(a) || math.IsNaN(b) {
  47. return false
  48. }
  49. if math.Signbit(a) != math.Signbit(b) {
  50. return math.Float64bits(math.Abs(a))+math.Float64bits(math.Abs(b)) <= uint64(ulp)
  51. }
  52. return ulpDiff(math.Float64bits(a), math.Float64bits(b)) <= uint64(ulp)
  53. }
  54. func ulpDiff(a, b uint64) uint64 {
  55. if a > b {
  56. return a - b
  57. }
  58. return b - a
  59. }
  60. const (
  61. nanBits = 0x7ff8000000000000
  62. nanMask = 0xfff8000000000000
  63. )
  64. // NaNWith returns an IEEE 754 "quiet not-a-number" value with the
  65. // payload specified in the low 51 bits of payload.
  66. // The NaN returned by math.NaN has a bit pattern equal to NaNWith(1).
  67. func NaNWith(payload uint64) float64 {
  68. return math.Float64frombits(nanBits | (payload &^ nanMask))
  69. }
  70. // NaNPayload returns the lowest 51 bits payload of an IEEE 754 "quiet
  71. // not-a-number". For values of f other than quiet-NaN, NaNPayload
  72. // returns zero and false.
  73. func NaNPayload(f float64) (payload uint64, ok bool) {
  74. b := math.Float64bits(f)
  75. if b&nanBits != nanBits {
  76. return 0, false
  77. }
  78. return b &^ nanMask, true
  79. }
  80. // ParseWithNA converts the string s to a float64 in value.
  81. // If s equals missing, weight is returned as 0, otherwise 1.
  82. func ParseWithNA(s, missing string) (value, weight float64, err error) {
  83. if s == missing {
  84. return 0, 0, nil
  85. }
  86. value, err = strconv.ParseFloat(s, 64)
  87. if err == nil {
  88. weight = 1
  89. }
  90. return value, weight, err
  91. }
  92. // Round returns the half away from zero rounded value of x with prec precision.
  93. //
  94. // Special cases are:
  95. // Round(±0) = +0
  96. // Round(±Inf) = ±Inf
  97. // Round(NaN) = NaN
  98. func Round(x float64, prec int) float64 {
  99. if x == 0 {
  100. // Make sure zero is returned
  101. // without the negative bit set.
  102. return 0
  103. }
  104. // Fast path for positive precision on integers.
  105. if prec >= 0 && x == math.Trunc(x) {
  106. return x
  107. }
  108. pow := math.Pow10(prec)
  109. intermed := x * pow
  110. if math.IsInf(intermed, 0) {
  111. return x
  112. }
  113. if x < 0 {
  114. x = math.Ceil(intermed - 0.5)
  115. } else {
  116. x = math.Floor(intermed + 0.5)
  117. }
  118. if x == 0 {
  119. return 0
  120. }
  121. return x / pow
  122. }
  123. // RoundEven returns the half even rounded value of x with prec precision.
  124. //
  125. // Special cases are:
  126. // RoundEven(±0) = +0
  127. // RoundEven(±Inf) = ±Inf
  128. // RoundEven(NaN) = NaN
  129. func RoundEven(x float64, prec int) float64 {
  130. if x == 0 {
  131. // Make sure zero is returned
  132. // without the negative bit set.
  133. return 0
  134. }
  135. // Fast path for positive precision on integers.
  136. if prec >= 0 && x == math.Trunc(x) {
  137. return x
  138. }
  139. pow := math.Pow10(prec)
  140. intermed := x * pow
  141. if math.IsInf(intermed, 0) {
  142. return x
  143. }
  144. if isHalfway(intermed) {
  145. correction, _ := math.Modf(math.Mod(intermed, 2))
  146. intermed += correction
  147. if intermed > 0 {
  148. x = math.Floor(intermed)
  149. } else {
  150. x = math.Ceil(intermed)
  151. }
  152. } else {
  153. if x < 0 {
  154. x = math.Ceil(intermed - 0.5)
  155. } else {
  156. x = math.Floor(intermed + 0.5)
  157. }
  158. }
  159. if x == 0 {
  160. return 0
  161. }
  162. return x / pow
  163. }
  164. func isHalfway(x float64) bool {
  165. _, frac := math.Modf(x)
  166. frac = math.Abs(frac)
  167. return frac == 0.5 || (math.Nextafter(frac, math.Inf(-1)) < 0.5 && math.Nextafter(frac, math.Inf(1)) > 0.5)
  168. }
  169. // Same returns true when the inputs have the same value, allowing NaN equality.
  170. func Same(a, b float64) bool {
  171. return a == b || (math.IsNaN(a) && math.IsNaN(b))
  172. }