floats.go 22 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934
  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 floats
  5. import (
  6. "errors"
  7. "math"
  8. "sort"
  9. "strconv"
  10. "gonum.org/v1/gonum/internal/asm/f64"
  11. )
  12. // Add adds, element-wise, the elements of s and dst, and stores in dst.
  13. // Panics if the lengths of dst and s do not match.
  14. func Add(dst, s []float64) {
  15. if len(dst) != len(s) {
  16. panic("floats: length of the slices do not match")
  17. }
  18. f64.AxpyUnitaryTo(dst, 1, s, dst)
  19. }
  20. // AddTo adds, element-wise, the elements of s and t and
  21. // stores the result in dst. Panics if the lengths of s, t and dst do not match.
  22. func AddTo(dst, s, t []float64) []float64 {
  23. if len(s) != len(t) {
  24. panic("floats: length of adders do not match")
  25. }
  26. if len(dst) != len(s) {
  27. panic("floats: length of destination does not match length of adder")
  28. }
  29. f64.AxpyUnitaryTo(dst, 1, s, t)
  30. return dst
  31. }
  32. // AddConst adds the scalar c to all of the values in dst.
  33. func AddConst(c float64, dst []float64) {
  34. f64.AddConst(c, dst)
  35. }
  36. // AddScaled performs dst = dst + alpha * s.
  37. // It panics if the lengths of dst and s are not equal.
  38. func AddScaled(dst []float64, alpha float64, s []float64) {
  39. if len(dst) != len(s) {
  40. panic("floats: length of destination and source to not match")
  41. }
  42. f64.AxpyUnitaryTo(dst, alpha, s, dst)
  43. }
  44. // AddScaledTo performs dst = y + alpha * s, where alpha is a scalar,
  45. // and dst, y and s are all slices.
  46. // It panics if the lengths of dst, y, and s are not equal.
  47. //
  48. // At the return of the function, dst[i] = y[i] + alpha * s[i]
  49. func AddScaledTo(dst, y []float64, alpha float64, s []float64) []float64 {
  50. if len(dst) != len(s) || len(dst) != len(y) {
  51. panic("floats: lengths of slices do not match")
  52. }
  53. f64.AxpyUnitaryTo(dst, alpha, s, y)
  54. return dst
  55. }
  56. // argsort is a helper that implements sort.Interface, as used by
  57. // Argsort.
  58. type argsort struct {
  59. s []float64
  60. inds []int
  61. }
  62. func (a argsort) Len() int {
  63. return len(a.s)
  64. }
  65. func (a argsort) Less(i, j int) bool {
  66. return a.s[i] < a.s[j]
  67. }
  68. func (a argsort) Swap(i, j int) {
  69. a.s[i], a.s[j] = a.s[j], a.s[i]
  70. a.inds[i], a.inds[j] = a.inds[j], a.inds[i]
  71. }
  72. // Argsort sorts the elements of dst while tracking their original order.
  73. // At the conclusion of Argsort, dst will contain the original elements of dst
  74. // but sorted in increasing order, and inds will contain the original position
  75. // of the elements in the slice such that dst[i] = origDst[inds[i]].
  76. // It panics if the lengths of dst and inds do not match.
  77. func Argsort(dst []float64, inds []int) {
  78. if len(dst) != len(inds) {
  79. panic("floats: length of inds does not match length of slice")
  80. }
  81. for i := range dst {
  82. inds[i] = i
  83. }
  84. a := argsort{s: dst, inds: inds}
  85. sort.Sort(a)
  86. }
  87. // Count applies the function f to every element of s and returns the number
  88. // of times the function returned true.
  89. func Count(f func(float64) bool, s []float64) int {
  90. var n int
  91. for _, val := range s {
  92. if f(val) {
  93. n++
  94. }
  95. }
  96. return n
  97. }
  98. // CumProd finds the cumulative product of the first i elements in
  99. // s and puts them in place into the ith element of the
  100. // destination dst. A panic will occur if the lengths of arguments
  101. // do not match.
  102. //
  103. // At the return of the function, dst[i] = s[i] * s[i-1] * s[i-2] * ...
  104. func CumProd(dst, s []float64) []float64 {
  105. if len(dst) != len(s) {
  106. panic("floats: length of destination does not match length of the source")
  107. }
  108. if len(dst) == 0 {
  109. return dst
  110. }
  111. return f64.CumProd(dst, s)
  112. }
  113. // CumSum finds the cumulative sum of the first i elements in
  114. // s and puts them in place into the ith element of the
  115. // destination dst. A panic will occur if the lengths of arguments
  116. // do not match.
  117. //
  118. // At the return of the function, dst[i] = s[i] + s[i-1] + s[i-2] + ...
  119. func CumSum(dst, s []float64) []float64 {
  120. if len(dst) != len(s) {
  121. panic("floats: length of destination does not match length of the source")
  122. }
  123. if len(dst) == 0 {
  124. return dst
  125. }
  126. return f64.CumSum(dst, s)
  127. }
  128. // Distance computes the L-norm of s - t. See Norm for special cases.
  129. // A panic will occur if the lengths of s and t do not match.
  130. func Distance(s, t []float64, L float64) float64 {
  131. if len(s) != len(t) {
  132. panic("floats: slice lengths do not match")
  133. }
  134. if len(s) == 0 {
  135. return 0
  136. }
  137. var norm float64
  138. if L == 2 {
  139. for i, v := range s {
  140. diff := t[i] - v
  141. norm = math.Hypot(norm, diff)
  142. }
  143. return norm
  144. }
  145. if L == 1 {
  146. for i, v := range s {
  147. norm += math.Abs(t[i] - v)
  148. }
  149. return norm
  150. }
  151. if math.IsInf(L, 1) {
  152. for i, v := range s {
  153. absDiff := math.Abs(t[i] - v)
  154. if absDiff > norm {
  155. norm = absDiff
  156. }
  157. }
  158. return norm
  159. }
  160. for i, v := range s {
  161. norm += math.Pow(math.Abs(t[i]-v), L)
  162. }
  163. return math.Pow(norm, 1/L)
  164. }
  165. // Div performs element-wise division dst / s
  166. // and stores the value in dst. It panics if the
  167. // lengths of s and t are not equal.
  168. func Div(dst, s []float64) {
  169. if len(dst) != len(s) {
  170. panic("floats: slice lengths do not match")
  171. }
  172. f64.Div(dst, s)
  173. }
  174. // DivTo performs element-wise division s / t
  175. // and stores the value in dst. It panics if the
  176. // lengths of s, t, and dst are not equal.
  177. func DivTo(dst, s, t []float64) []float64 {
  178. if len(s) != len(t) || len(dst) != len(t) {
  179. panic("floats: slice lengths do not match")
  180. }
  181. return f64.DivTo(dst, s, t)
  182. }
  183. // Dot computes the dot product of s1 and s2, i.e.
  184. // sum_{i = 1}^N s1[i]*s2[i].
  185. // A panic will occur if lengths of arguments do not match.
  186. func Dot(s1, s2 []float64) float64 {
  187. if len(s1) != len(s2) {
  188. panic("floats: lengths of the slices do not match")
  189. }
  190. return f64.DotUnitary(s1, s2)
  191. }
  192. // Equal returns true if the slices have equal lengths and
  193. // all elements are numerically identical.
  194. func Equal(s1, s2 []float64) bool {
  195. if len(s1) != len(s2) {
  196. return false
  197. }
  198. for i, val := range s1 {
  199. if s2[i] != val {
  200. return false
  201. }
  202. }
  203. return true
  204. }
  205. // EqualApprox returns true if the slices have equal lengths and
  206. // all element pairs have an absolute tolerance less than tol or a
  207. // relative tolerance less than tol.
  208. func EqualApprox(s1, s2 []float64, tol float64) bool {
  209. if len(s1) != len(s2) {
  210. return false
  211. }
  212. for i, a := range s1 {
  213. if !EqualWithinAbsOrRel(a, s2[i], tol, tol) {
  214. return false
  215. }
  216. }
  217. return true
  218. }
  219. // EqualFunc returns true if the slices have the same lengths
  220. // and the function returns true for all element pairs.
  221. func EqualFunc(s1, s2 []float64, f func(float64, float64) bool) bool {
  222. if len(s1) != len(s2) {
  223. return false
  224. }
  225. for i, val := range s1 {
  226. if !f(val, s2[i]) {
  227. return false
  228. }
  229. }
  230. return true
  231. }
  232. // EqualWithinAbs returns true if a and b have an absolute
  233. // difference of less than tol.
  234. func EqualWithinAbs(a, b, tol float64) bool {
  235. return a == b || math.Abs(a-b) <= tol
  236. }
  237. const minNormalFloat64 = 2.2250738585072014e-308
  238. // EqualWithinRel returns true if the difference between a and b
  239. // is not greater than tol times the greater value.
  240. func EqualWithinRel(a, b, tol float64) bool {
  241. if a == b {
  242. return true
  243. }
  244. delta := math.Abs(a - b)
  245. if delta <= minNormalFloat64 {
  246. return delta <= tol*minNormalFloat64
  247. }
  248. // We depend on the division in this relationship to identify
  249. // infinities (we rely on the NaN to fail the test) otherwise
  250. // we compare Infs of the same sign and evaluate Infs as equal
  251. // independent of sign.
  252. return delta/math.Max(math.Abs(a), math.Abs(b)) <= tol
  253. }
  254. // EqualWithinAbsOrRel returns true if a and b are equal to within
  255. // the absolute tolerance.
  256. func EqualWithinAbsOrRel(a, b, absTol, relTol float64) bool {
  257. if EqualWithinAbs(a, b, absTol) {
  258. return true
  259. }
  260. return EqualWithinRel(a, b, relTol)
  261. }
  262. // EqualWithinULP returns true if a and b are equal to within
  263. // the specified number of floating point units in the last place.
  264. func EqualWithinULP(a, b float64, ulp uint) bool {
  265. if a == b {
  266. return true
  267. }
  268. if math.IsNaN(a) || math.IsNaN(b) {
  269. return false
  270. }
  271. if math.Signbit(a) != math.Signbit(b) {
  272. return math.Float64bits(math.Abs(a))+math.Float64bits(math.Abs(b)) <= uint64(ulp)
  273. }
  274. return ulpDiff(math.Float64bits(a), math.Float64bits(b)) <= uint64(ulp)
  275. }
  276. func ulpDiff(a, b uint64) uint64 {
  277. if a > b {
  278. return a - b
  279. }
  280. return b - a
  281. }
  282. // EqualLengths returns true if all of the slices have equal length,
  283. // and false otherwise. Returns true if there are no input slices.
  284. func EqualLengths(slices ...[]float64) bool {
  285. // This length check is needed: http://play.golang.org/p/sdty6YiLhM
  286. if len(slices) == 0 {
  287. return true
  288. }
  289. l := len(slices[0])
  290. for i := 1; i < len(slices); i++ {
  291. if len(slices[i]) != l {
  292. return false
  293. }
  294. }
  295. return true
  296. }
  297. // Find applies f to every element of s and returns the indices of the first
  298. // k elements for which the f returns true, or all such elements
  299. // if k < 0.
  300. // Find will reslice inds to have 0 length, and will append
  301. // found indices to inds.
  302. // If k > 0 and there are fewer than k elements in s satisfying f,
  303. // all of the found elements will be returned along with an error.
  304. // At the return of the function, the input inds will be in an undetermined state.
  305. func Find(inds []int, f func(float64) bool, s []float64, k int) ([]int, error) {
  306. // inds is also returned to allow for calling with nil
  307. // Reslice inds to have zero length
  308. inds = inds[:0]
  309. // If zero elements requested, can just return
  310. if k == 0 {
  311. return inds, nil
  312. }
  313. // If k < 0, return all of the found indices
  314. if k < 0 {
  315. for i, val := range s {
  316. if f(val) {
  317. inds = append(inds, i)
  318. }
  319. }
  320. return inds, nil
  321. }
  322. // Otherwise, find the first k elements
  323. nFound := 0
  324. for i, val := range s {
  325. if f(val) {
  326. inds = append(inds, i)
  327. nFound++
  328. if nFound == k {
  329. return inds, nil
  330. }
  331. }
  332. }
  333. // Finished iterating over the loop, which means k elements were not found
  334. return inds, errors.New("floats: insufficient elements found")
  335. }
  336. // HasNaN returns true if the slice s has any values that are NaN and false
  337. // otherwise.
  338. func HasNaN(s []float64) bool {
  339. for _, v := range s {
  340. if math.IsNaN(v) {
  341. return true
  342. }
  343. }
  344. return false
  345. }
  346. // LogSpan returns a set of n equally spaced points in log space between,
  347. // l and u where N is equal to len(dst). The first element of the
  348. // resulting dst will be l and the final element of dst will be u.
  349. // Panics if len(dst) < 2
  350. // Note that this call will return NaNs if either l or u are negative, and
  351. // will return all zeros if l or u is zero.
  352. // Also returns the mutated slice dst, so that it can be used in range, like:
  353. //
  354. // for i, x := range LogSpan(dst, l, u) { ... }
  355. func LogSpan(dst []float64, l, u float64) []float64 {
  356. Span(dst, math.Log(l), math.Log(u))
  357. for i := range dst {
  358. dst[i] = math.Exp(dst[i])
  359. }
  360. return dst
  361. }
  362. // LogSumExp returns the log of the sum of the exponentials of the values in s.
  363. // Panics if s is an empty slice.
  364. func LogSumExp(s []float64) float64 {
  365. // Want to do this in a numerically stable way which avoids
  366. // overflow and underflow
  367. // First, find the maximum value in the slice.
  368. maxval := Max(s)
  369. if math.IsInf(maxval, 0) {
  370. // If it's infinity either way, the logsumexp will be infinity as well
  371. // returning now avoids NaNs
  372. return maxval
  373. }
  374. var lse float64
  375. // Compute the sumexp part
  376. for _, val := range s {
  377. lse += math.Exp(val - maxval)
  378. }
  379. // Take the log and add back on the constant taken out
  380. return math.Log(lse) + maxval
  381. }
  382. // Max returns the maximum value in the input slice. If the slice is empty, Max will panic.
  383. func Max(s []float64) float64 {
  384. return s[MaxIdx(s)]
  385. }
  386. // MaxIdx returns the index of the maximum value in the input slice. If several
  387. // entries have the maximum value, the first such index is returned. If the slice
  388. // is empty, MaxIdx will panic.
  389. func MaxIdx(s []float64) int {
  390. if len(s) == 0 {
  391. panic("floats: zero slice length")
  392. }
  393. max := math.NaN()
  394. var ind int
  395. for i, v := range s {
  396. if math.IsNaN(v) {
  397. continue
  398. }
  399. if v > max || math.IsNaN(max) {
  400. max = v
  401. ind = i
  402. }
  403. }
  404. return ind
  405. }
  406. // Min returns the maximum value in the input slice. If the slice is empty, Min will panic.
  407. func Min(s []float64) float64 {
  408. return s[MinIdx(s)]
  409. }
  410. // MinIdx returns the index of the minimum value in the input slice. If several
  411. // entries have the maximum value, the first such index is returned. If the slice
  412. // is empty, MinIdx will panic.
  413. func MinIdx(s []float64) int {
  414. if len(s) == 0 {
  415. panic("floats: zero slice length")
  416. }
  417. min := math.NaN()
  418. var ind int
  419. for i, v := range s {
  420. if math.IsNaN(v) {
  421. continue
  422. }
  423. if v < min || math.IsNaN(min) {
  424. min = v
  425. ind = i
  426. }
  427. }
  428. return ind
  429. }
  430. // Mul performs element-wise multiplication between dst
  431. // and s and stores the value in dst. Panics if the
  432. // lengths of s and t are not equal.
  433. func Mul(dst, s []float64) {
  434. if len(dst) != len(s) {
  435. panic("floats: slice lengths do not match")
  436. }
  437. for i, val := range s {
  438. dst[i] *= val
  439. }
  440. }
  441. // MulTo performs element-wise multiplication between s
  442. // and t and stores the value in dst. Panics if the
  443. // lengths of s, t, and dst are not equal.
  444. func MulTo(dst, s, t []float64) []float64 {
  445. if len(s) != len(t) || len(dst) != len(t) {
  446. panic("floats: slice lengths do not match")
  447. }
  448. for i, val := range t {
  449. dst[i] = val * s[i]
  450. }
  451. return dst
  452. }
  453. const (
  454. nanBits = 0x7ff8000000000000
  455. nanMask = 0xfff8000000000000
  456. )
  457. // NaNWith returns an IEEE 754 "quiet not-a-number" value with the
  458. // payload specified in the low 51 bits of payload.
  459. // The NaN returned by math.NaN has a bit pattern equal to NaNWith(1).
  460. func NaNWith(payload uint64) float64 {
  461. return math.Float64frombits(nanBits | (payload &^ nanMask))
  462. }
  463. // NaNPayload returns the lowest 51 bits payload of an IEEE 754 "quiet
  464. // not-a-number". For values of f other than quiet-NaN, NaNPayload
  465. // returns zero and false.
  466. func NaNPayload(f float64) (payload uint64, ok bool) {
  467. b := math.Float64bits(f)
  468. if b&nanBits != nanBits {
  469. return 0, false
  470. }
  471. return b &^ nanMask, true
  472. }
  473. // NearestIdx returns the index of the element in s
  474. // whose value is nearest to v. If several such
  475. // elements exist, the lowest index is returned.
  476. // NearestIdx panics if len(s) == 0.
  477. func NearestIdx(s []float64, v float64) int {
  478. if len(s) == 0 {
  479. panic("floats: zero length slice")
  480. }
  481. switch {
  482. case math.IsNaN(v):
  483. return 0
  484. case math.IsInf(v, 1):
  485. return MaxIdx(s)
  486. case math.IsInf(v, -1):
  487. return MinIdx(s)
  488. }
  489. var ind int
  490. dist := math.NaN()
  491. for i, val := range s {
  492. newDist := math.Abs(v - val)
  493. // A NaN distance will not be closer.
  494. if math.IsNaN(newDist) {
  495. continue
  496. }
  497. if newDist < dist || math.IsNaN(dist) {
  498. dist = newDist
  499. ind = i
  500. }
  501. }
  502. return ind
  503. }
  504. // NearestIdxForSpan return the index of a hypothetical vector created
  505. // by Span with length n and bounds l and u whose value is closest
  506. // to v. That is, NearestIdxForSpan(n, l, u, v) is equivalent to
  507. // Nearest(Span(make([]float64, n),l,u),v) without an allocation.
  508. // NearestIdxForSpan panics if n is less than two.
  509. func NearestIdxForSpan(n int, l, u float64, v float64) int {
  510. if n <= 1 {
  511. panic("floats: span must have length >1")
  512. }
  513. if math.IsNaN(v) {
  514. return 0
  515. }
  516. // Special cases for Inf and NaN.
  517. switch {
  518. case math.IsNaN(l) && !math.IsNaN(u):
  519. return n - 1
  520. case math.IsNaN(u):
  521. return 0
  522. case math.IsInf(l, 0) && math.IsInf(u, 0):
  523. if l == u {
  524. return 0
  525. }
  526. if n%2 == 1 {
  527. if !math.IsInf(v, 0) {
  528. return n / 2
  529. }
  530. if math.Copysign(1, v) == math.Copysign(1, l) {
  531. return 0
  532. }
  533. return n/2 + 1
  534. }
  535. if math.Copysign(1, v) == math.Copysign(1, l) {
  536. return 0
  537. }
  538. return n / 2
  539. case math.IsInf(l, 0):
  540. if v == l {
  541. return 0
  542. }
  543. return n - 1
  544. case math.IsInf(u, 0):
  545. if v == u {
  546. return n - 1
  547. }
  548. return 0
  549. case math.IsInf(v, -1):
  550. if l <= u {
  551. return 0
  552. }
  553. return n - 1
  554. case math.IsInf(v, 1):
  555. if u <= l {
  556. return 0
  557. }
  558. return n - 1
  559. }
  560. // Special cases for v outside (l, u) and (u, l).
  561. switch {
  562. case l < u:
  563. if v <= l {
  564. return 0
  565. }
  566. if v >= u {
  567. return n - 1
  568. }
  569. case l > u:
  570. if v >= l {
  571. return 0
  572. }
  573. if v <= u {
  574. return n - 1
  575. }
  576. default:
  577. return 0
  578. }
  579. // Can't guarantee anything about exactly halfway between
  580. // because of floating point weirdness.
  581. return int((float64(n)-1)/(u-l)*(v-l) + 0.5)
  582. }
  583. // Norm returns the L norm of the slice S, defined as
  584. // (sum_{i=1}^N s[i]^L)^{1/L}
  585. // Special cases:
  586. // L = math.Inf(1) gives the maximum absolute value.
  587. // Does not correctly compute the zero norm (use Count).
  588. func Norm(s []float64, L float64) float64 {
  589. // Should this complain if L is not positive?
  590. // Should this be done in log space for better numerical stability?
  591. // would be more cost
  592. // maybe only if L is high?
  593. if len(s) == 0 {
  594. return 0
  595. }
  596. if L == 2 {
  597. twoNorm := math.Abs(s[0])
  598. for i := 1; i < len(s); i++ {
  599. twoNorm = math.Hypot(twoNorm, s[i])
  600. }
  601. return twoNorm
  602. }
  603. var norm float64
  604. if L == 1 {
  605. for _, val := range s {
  606. norm += math.Abs(val)
  607. }
  608. return norm
  609. }
  610. if math.IsInf(L, 1) {
  611. for _, val := range s {
  612. norm = math.Max(norm, math.Abs(val))
  613. }
  614. return norm
  615. }
  616. for _, val := range s {
  617. norm += math.Pow(math.Abs(val), L)
  618. }
  619. return math.Pow(norm, 1/L)
  620. }
  621. // ParseWithNA converts the string s to a float64 in v.
  622. // If s equals missing, w is returned as 0, otherwise 1.
  623. func ParseWithNA(s, missing string) (v, w float64, err error) {
  624. if s == missing {
  625. return 0, 0, nil
  626. }
  627. v, err = strconv.ParseFloat(s, 64)
  628. if err == nil {
  629. w = 1
  630. }
  631. return v, w, err
  632. }
  633. // Prod returns the product of the elements of the slice.
  634. // Returns 1 if len(s) = 0.
  635. func Prod(s []float64) float64 {
  636. prod := 1.0
  637. for _, val := range s {
  638. prod *= val
  639. }
  640. return prod
  641. }
  642. // Reverse reverses the order of elements in the slice.
  643. func Reverse(s []float64) {
  644. for i, j := 0, len(s)-1; i < j; i, j = i+1, j-1 {
  645. s[i], s[j] = s[j], s[i]
  646. }
  647. }
  648. // Round returns the half away from zero rounded value of x with prec precision.
  649. //
  650. // Special cases are:
  651. // Round(±0) = +0
  652. // Round(±Inf) = ±Inf
  653. // Round(NaN) = NaN
  654. func Round(x float64, prec int) float64 {
  655. if x == 0 {
  656. // Make sure zero is returned
  657. // without the negative bit set.
  658. return 0
  659. }
  660. // Fast path for positive precision on integers.
  661. if prec >= 0 && x == math.Trunc(x) {
  662. return x
  663. }
  664. pow := math.Pow10(prec)
  665. intermed := x * pow
  666. if math.IsInf(intermed, 0) {
  667. return x
  668. }
  669. if x < 0 {
  670. x = math.Ceil(intermed - 0.5)
  671. } else {
  672. x = math.Floor(intermed + 0.5)
  673. }
  674. if x == 0 {
  675. return 0
  676. }
  677. return x / pow
  678. }
  679. // RoundEven returns the half even rounded value of x with prec precision.
  680. //
  681. // Special cases are:
  682. // RoundEven(±0) = +0
  683. // RoundEven(±Inf) = ±Inf
  684. // RoundEven(NaN) = NaN
  685. func RoundEven(x float64, prec int) float64 {
  686. if x == 0 {
  687. // Make sure zero is returned
  688. // without the negative bit set.
  689. return 0
  690. }
  691. // Fast path for positive precision on integers.
  692. if prec >= 0 && x == math.Trunc(x) {
  693. return x
  694. }
  695. pow := math.Pow10(prec)
  696. intermed := x * pow
  697. if math.IsInf(intermed, 0) {
  698. return x
  699. }
  700. if isHalfway(intermed) {
  701. correction, _ := math.Modf(math.Mod(intermed, 2))
  702. intermed += correction
  703. if intermed > 0 {
  704. x = math.Floor(intermed)
  705. } else {
  706. x = math.Ceil(intermed)
  707. }
  708. } else {
  709. if x < 0 {
  710. x = math.Ceil(intermed - 0.5)
  711. } else {
  712. x = math.Floor(intermed + 0.5)
  713. }
  714. }
  715. if x == 0 {
  716. return 0
  717. }
  718. return x / pow
  719. }
  720. func isHalfway(x float64) bool {
  721. _, frac := math.Modf(x)
  722. frac = math.Abs(frac)
  723. return frac == 0.5 || (math.Nextafter(frac, math.Inf(-1)) < 0.5 && math.Nextafter(frac, math.Inf(1)) > 0.5)
  724. }
  725. // Same returns true if the input slices have the same length and the all elements
  726. // have the same value with NaN treated as the same.
  727. func Same(s, t []float64) bool {
  728. if len(s) != len(t) {
  729. return false
  730. }
  731. for i, v := range s {
  732. w := t[i]
  733. if v != w && !(math.IsNaN(v) && math.IsNaN(w)) {
  734. return false
  735. }
  736. }
  737. return true
  738. }
  739. // Scale multiplies every element in dst by the scalar c.
  740. func Scale(c float64, dst []float64) {
  741. if len(dst) > 0 {
  742. f64.ScalUnitary(c, dst)
  743. }
  744. }
  745. // ScaleTo multiplies the elements in s by c and stores the result in dst.
  746. func ScaleTo(dst []float64, c float64, s []float64) []float64 {
  747. if len(dst) != len(s) {
  748. panic("floats: lengths of slices do not match")
  749. }
  750. if len(dst) > 0 {
  751. f64.ScalUnitaryTo(dst, c, s)
  752. }
  753. return dst
  754. }
  755. // Span returns a set of N equally spaced points between l and u, where N
  756. // is equal to the length of the destination. The first element of the destination
  757. // is l, the final element of the destination is u.
  758. //
  759. // Panics if len(dst) < 2.
  760. //
  761. // Span also returns the mutated slice dst, so that it can be used in range expressions,
  762. // like:
  763. //
  764. // for i, x := range Span(dst, l, u) { ... }
  765. func Span(dst []float64, l, u float64) []float64 {
  766. n := len(dst)
  767. if n < 2 {
  768. panic("floats: destination must have length >1")
  769. }
  770. // Special cases for Inf and NaN.
  771. switch {
  772. case math.IsNaN(l):
  773. for i := range dst[:len(dst)-1] {
  774. dst[i] = math.NaN()
  775. }
  776. dst[len(dst)-1] = u
  777. return dst
  778. case math.IsNaN(u):
  779. for i := range dst[1:] {
  780. dst[i+1] = math.NaN()
  781. }
  782. dst[0] = l
  783. return dst
  784. case math.IsInf(l, 0) && math.IsInf(u, 0):
  785. for i := range dst[:len(dst)/2] {
  786. dst[i] = l
  787. dst[len(dst)-i-1] = u
  788. }
  789. if len(dst)%2 == 1 {
  790. if l != u {
  791. dst[len(dst)/2] = 0
  792. } else {
  793. dst[len(dst)/2] = l
  794. }
  795. }
  796. return dst
  797. case math.IsInf(l, 0):
  798. for i := range dst[:len(dst)-1] {
  799. dst[i] = l
  800. }
  801. dst[len(dst)-1] = u
  802. return dst
  803. case math.IsInf(u, 0):
  804. for i := range dst[1:] {
  805. dst[i+1] = u
  806. }
  807. dst[0] = l
  808. return dst
  809. }
  810. step := (u - l) / float64(n-1)
  811. for i := range dst {
  812. dst[i] = l + step*float64(i)
  813. }
  814. return dst
  815. }
  816. // Sub subtracts, element-wise, the elements of s from dst. Panics if
  817. // the lengths of dst and s do not match.
  818. func Sub(dst, s []float64) {
  819. if len(dst) != len(s) {
  820. panic("floats: length of the slices do not match")
  821. }
  822. f64.AxpyUnitaryTo(dst, -1, s, dst)
  823. }
  824. // SubTo subtracts, element-wise, the elements of t from s and
  825. // stores the result in dst. Panics if the lengths of s, t and dst do not match.
  826. func SubTo(dst, s, t []float64) []float64 {
  827. if len(s) != len(t) {
  828. panic("floats: length of subtractor and subtractee do not match")
  829. }
  830. if len(dst) != len(s) {
  831. panic("floats: length of destination does not match length of subtractor")
  832. }
  833. f64.AxpyUnitaryTo(dst, -1, t, s)
  834. return dst
  835. }
  836. // Sum returns the sum of the elements of the slice.
  837. func Sum(s []float64) float64 {
  838. return f64.Sum(s)
  839. }
  840. // Within returns the first index i where s[i] <= v < s[i+1]. Within panics if:
  841. // - len(s) < 2
  842. // - s is not sorted
  843. func Within(s []float64, v float64) int {
  844. if len(s) < 2 {
  845. panic("floats: slice length less than 2")
  846. }
  847. if !sort.Float64sAreSorted(s) {
  848. panic("floats: input slice not sorted")
  849. }
  850. if v < s[0] || v >= s[len(s)-1] || math.IsNaN(v) {
  851. return -1
  852. }
  853. for i, f := range s[1:] {
  854. if v < f {
  855. return i
  856. }
  857. }
  858. return -1
  859. }