vendor/gonum.org/v1/gonum/floats/floats.go

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