File: ./doc.go
   1 /*
   2 # mathplus
   3 
   4 This is an add-on to the stdlib package `math`, with statistics-gathering
   5 functionality and plenty more math-related functions.
   6 
   7 This package also wraps types/methods from math/rand with nil-safe ones, which
   8 act on the default value-generator when nil pointers are used.
   9 */
  10 package mathplus

     File: ./functions.go
   1 package mathplus
   2 
   3 import "math"
   4 
   5 func Decade(x float64) float64  { return 10 * math.Floor(0.1*x) }
   6 func Century(x float64) float64 { return 100 * math.Floor(0.01*x) }
   7 
   8 func Round1(x float64) float64 { return 1e-1 * math.Round(1e+1*x) }
   9 func Round2(x float64) float64 { return 1e-2 * math.Round(1e+2*x) }
  10 func Round3(x float64) float64 { return 1e-3 * math.Round(1e+3*x) }
  11 func Round4(x float64) float64 { return 1e-4 * math.Round(1e+4*x) }
  12 func Round5(x float64) float64 { return 1e-5 * math.Round(1e+5*x) }
  13 func Round6(x float64) float64 { return 1e-6 * math.Round(1e+6*x) }
  14 func Round7(x float64) float64 { return 1e-7 * math.Round(1e+7*x) }
  15 func Round8(x float64) float64 { return 1e-8 * math.Round(1e+8*x) }
  16 func Round9(x float64) float64 { return 1e-9 * math.Round(1e+9*x) }
  17 
  18 // Ln calculates the natural logarithm.
  19 func Ln(x float64) float64 {
  20     return math.Log(x)
  21 }
  22 
  23 // Log calculates logarithms using the base given.
  24 // func Log(base float64, x float64) float64 {
  25 //  return math.Log(x) / math.Log(base)
  26 // }
  27 
  28 // Round rounds the number given to the number of decimal places given.
  29 func Round(x float64, decimals int) float64 {
  30     k := math.Pow(10, float64(decimals))
  31     return math.Round(k*x) / k
  32 }
  33 
  34 // RoundBy rounds the number given to the unit-size given
  35 // func RoundBy(x float64, by float64) float64 {
  36 //  return math.Round(by*x) / by
  37 // }
  38 
  39 // FloorBy quantizes a number downward by the unit-size given
  40 // func FloorBy(x float64, by float64) float64 {
  41 //  mod := math.Mod(x, by)
  42 //  if x < 0 {
  43 //      if mod == 0 {
  44 //          return x - mod
  45 //      }
  46 //      return x - mod - by
  47 //  }
  48 //  return x - mod
  49 // }
  50 
  51 // Wrap linearly interpolates the number given to the range [0...1],
  52 // according to its continuous domain, specified by the limits given
  53 func Wrap(x float64, min, max float64) float64 {
  54     return (x - min) / (max - min)
  55 }
  56 
  57 // AnchoredWrap works like Wrap, except when both domain boundaries are positive,
  58 // the minimum becomes 0, and when both domain boundaries are negative, the
  59 // maximum becomes 0.
  60 func AnchoredWrap(x float64, min, max float64) float64 {
  61     if min > 0 && max > 0 {
  62         min = 0
  63     } else if max < 0 && min < 0 {
  64         max = 0
  65     }
  66     return (x - min) / (max - min)
  67 }
  68 
  69 // Unwrap undoes Wrap, by turning the [0...1] number given into its equivalent
  70 // in the new domain given.
  71 func Unwrap(x float64, min, max float64) float64 {
  72     return (max-min)*x + min
  73 }
  74 
  75 // Clamp constrains the domain of the value given
  76 func Clamp(x float64, min, max float64) float64 {
  77     return math.Min(math.Max(x, min), max)
  78 }
  79 
  80 // Scale transforms a number according to its position in [xMin, xMax] into its
  81 // correspondingly-positioned counterpart in [yMin, yMax]: if value isn't in its
  82 // assumed domain, its result will be extrapolated accordingly
  83 func Scale(x float64, xmin, xmax, ymin, ymax float64) float64 {
  84     k := (x - xmin) / (xmax - xmin)
  85     return (ymax-ymin)*k + ymin
  86 }
  87 
  88 // AnchoredScale works like Scale, except when both domain boundaries are positive,
  89 // the minimum becomes 0, and when both domain boundaries are negative, the maximum
  90 // becomes 0. This allows for proportionally-correct scaling of quantities, such
  91 // as when showing data visually.
  92 func AnchoredScale(x float64, xmin, xmax, ymin, ymax float64) float64 {
  93     if xmin > 0 && xmax > 0 {
  94         xmin = 0
  95     } else if xmax < 0 && xmin < 0 {
  96         xmax = 0
  97     }
  98     k := (x - xmin) / (xmax - xmin)
  99     return (ymax-ymin)*k + ymin
 100 }
 101 
 102 // ForceRange is handy when trying to mold floating-point values into numbers
 103 // valid for JSON, since NaN and Infinity are replaced by the values given;
 104 // the infinity-replacement value is negated for negative infinity.
 105 func ForceRange(x float64, nan, inf float64) float64 {
 106     if math.IsNaN(x) {
 107         return nan
 108     }
 109     if math.IsInf(x, -1) {
 110         return -inf
 111     }
 112     if math.IsInf(x, +1) {
 113         return inf
 114     }
 115     return x
 116 }
 117 
 118 // Sign returns the standardized sign of a value, either as -1, 0, or +1: NaN
 119 // values stay as NaN, as is expected when using floating-point values.
 120 func Sign(x float64) float64 {
 121     if x > 0 {
 122         return +1
 123     }
 124     if x < 0 {
 125         return -1
 126     }
 127     if x == 0 {
 128         return 0
 129     }
 130     return math.NaN()
 131 }
 132 
 133 // IsInteger checks if a floating-point value is already rounded to an integer
 134 // value.
 135 func IsInteger(x float64) bool {
 136     _, frac := math.Modf(x)
 137     return frac == 0
 138 }
 139 
 140 func Square(x float64) float64 {
 141     return x * x
 142 }
 143 
 144 func Cube(x float64) float64 {
 145     return x * x * x
 146 }
 147 
 148 func RoundTo(x float64, decimals float64) float64 {
 149     k := math.Pow10(int(decimals))
 150     return math.Round(k*x) / k
 151 }
 152 
 153 func RelDiff(x, y float64) float64 {
 154     return (x - y) / y
 155 }
 156 
 157 // Bool booleanizes a number: 0 stays 0, and anything else becomes 1.
 158 func Bool(x float64) float64 {
 159     if x == 0 {
 160         return 0
 161     }
 162     return 1
 163 }
 164 
 165 // DeNaN replaces NaN with the alternative value given: regular values are
 166 // returned as given.
 167 func DeNaN(x float64, instead float64) float64 {
 168     if !math.IsNaN(x) {
 169         return x
 170     }
 171     return instead
 172 }
 173 
 174 // DeInf replaces either infinity with the alternative value given: regular
 175 // values are returned as given.
 176 func DeInf(x float64, instead float64) float64 {
 177     if !math.IsInf(x, 0) {
 178         return x
 179     }
 180     return instead
 181 }
 182 
 183 // Revalue replaces NaN and either infinity with the alternative value given:
 184 // regular values are returned as given.
 185 func Revalue(x float64, instead float64) float64 {
 186     if !math.IsNaN(x) && !math.IsInf(x, 0) {
 187         return x
 188     }
 189     return instead
 190 }
 191 
 192 // Linear evaluates a linear polynomial, the first arg being the main input,
 193 // followed by the polynomial coefficients in decreasing-power order.
 194 func Linear(x float64, a, b float64) float64 {
 195     return a*x + b
 196 }
 197 
 198 // Quadratic evaluates a 2nd-degree polynomial, the first arg being the main
 199 // input, followed by the polynomial coefficients in decreasing-power order.
 200 func Quadratic(x float64, a, b, c float64) float64 {
 201     return (a*x+b)*x + c
 202 }
 203 
 204 // Cubic evaluates a cubic polynomial, the first arg being the main input,
 205 // followed by the polynomial coefficients in decreasing-power order.
 206 func Cubic(x float64, a, b, c, d float64) float64 {
 207     return ((a*x+b)*x+c)*x + d
 208 }
 209 
 210 func LinearFMA(x float64, a, b float64) float64 {
 211     return math.FMA(x, a, b)
 212 }
 213 
 214 func QuadraticFMA(x float64, a, b, c float64) float64 {
 215     lin := math.FMA(x, a, b)
 216     return math.FMA(lin, x, c)
 217 }
 218 
 219 func CubicFMA(x float64, a, b, c, d float64) float64 {
 220     lin := math.FMA(x, a, b)
 221     quad := math.FMA(lin, x, c)
 222     return math.FMA(quad, x, d)
 223 }
 224 
 225 // Radians converts angular degrees into angular radians: 180 degrees are pi
 226 // pi radians.
 227 func Radians(deg float64) float64 {
 228     const k = math.Pi / 180
 229     return k * deg
 230 }
 231 
 232 // Degrees converts angular radians into angular degrees: pi radians are 180
 233 // degrees.
 234 func Degrees(rad float64) float64 {
 235     const k = 180 / math.Pi
 236     return k * rad
 237 }
 238 
 239 // Fract calculates the non-integer/fractional part of a number.
 240 func Fract(x float64) float64 {
 241     return x - math.Floor(x)
 242 }
 243 
 244 // Mix interpolates 2 numbers using a third number, used as an interpolation
 245 // coefficient. This parameter naturally falls in the range [0, 1], but doesn't
 246 // have to be: when given outside that range, the parameter can extrapolate in
 247 // either direction instead.
 248 func Mix(x, y float64, k float64) float64 {
 249     return (1-k)*(y-x) + x
 250 }
 251 
 252 // Step implements a step function with a parametric threshold.
 253 func Step(edge, x float64) float64 {
 254     if x < edge {
 255         return 0
 256     }
 257     return 1
 258 }
 259 
 260 // SmoothStep is like the `smoothstep` func found in GLSL, using a cubic
 261 // interpolator in the transition region.
 262 func SmoothStep(edge0, edge1, x float64) float64 {
 263     if x <= edge0 {
 264         return 0
 265     }
 266     if x >= edge1 {
 267         return 1
 268     }
 269 
 270     // use the cubic hermite interpolator 3x^2 - 2x^3 in the transition band
 271     return x * x * (3 - 2*x)
 272 }
 273 
 274 // Logistic approximates the math func of the same name.
 275 func Logistic(x float64) float64 {
 276     return 1 / (1 + math.Exp(-x))
 277 }
 278 
 279 // Sinc approximates the math func of the same name.
 280 func Sinc(x float64) float64 {
 281     if x != 0 {
 282         return math.Sin(x) / x
 283     }
 284     return 1
 285 }
 286 
 287 // Sum adds all the numbers in an array.
 288 func Sum(v ...float64) float64 {
 289     s := 0.0
 290     for _, f := range v {
 291         s += f
 292     }
 293     return s
 294 }
 295 
 296 // Product multiplies all the numbers in an array.
 297 func Product(v ...float64) float64 {
 298     p := 1.0
 299     for _, f := range v {
 300         p *= f
 301     }
 302     return p
 303 }
 304 
 305 // Length calculates the Euclidean length of the vector given.
 306 func Length(v ...float64) float64 {
 307     ss := 0.0
 308     for _, f := range v {
 309         ss += f * f
 310     }
 311     return math.Sqrt(ss)
 312 }
 313 
 314 // Dot calculates the dot product of 2 vectors
 315 func Dot(x []float64, y []float64) float64 {
 316     l := len(x)
 317     if len(y) < l {
 318         l = len(y)
 319     }
 320 
 321     dot := 0.0
 322     for i := 0; i < l; i++ {
 323         dot += x[i] * y[i]
 324     }
 325     return dot
 326 }
 327 
 328 // Min finds the minimum value from the numbers given.
 329 func Min(v ...float64) float64 {
 330     min := +math.Inf(+1)
 331     for _, f := range v {
 332         min = math.Min(min, f)
 333     }
 334     return min
 335 }
 336 
 337 // Max finds the maximum value from the numbers given.
 338 func Max(v ...float64) float64 {
 339     max := +math.Inf(-1)
 340     for _, f := range v {
 341         max = math.Max(max, f)
 342     }
 343     return max
 344 }
 345 
 346 // Hypot calculates the Euclidean n-dimensional hypothenuse from the numbers
 347 // given: all numbers can be lengths, or simply positional coordinates.
 348 func Hypot(v ...float64) float64 {
 349     sumsq := 0.0
 350     for _, f := range v {
 351         sumsq += f * f
 352     }
 353     return math.Sqrt(sumsq)
 354 }
 355 
 356 // Polyval evaluates a polynomial using Horner's algorithm. The array has all
 357 // the coefficients in textbook order, from the highest power down to the final
 358 // constant.
 359 func Polyval(x float64, v ...float64) float64 {
 360     if len(v) == 0 {
 361         return 0
 362     }
 363 
 364     x0 := x
 365     x = 1.0
 366     y := 0.0
 367     for i := len(v) - 1; i >= 0; i-- {
 368         y += v[i] * x
 369         x *= x0
 370     }
 371     return y
 372 }
 373 
 374 // LnGamma is a 1-input 1-output version of math.Lgamma from the stdlib.
 375 func LnGamma(x float64) float64 {
 376     y, s := math.Lgamma(x)
 377     if s < 0 {
 378         return math.NaN()
 379     }
 380     return y
 381 }
 382 
 383 // LnBeta calculates the natural-logarithm of the beta function.
 384 func LnBeta(x float64, y float64) float64 {
 385     return LnGamma(x) + LnGamma(y) - LnGamma(x+y)
 386 }
 387 
 388 // Beta calculates the beta function.
 389 func Beta(x float64, y float64) float64 {
 390     return math.Exp(LnBeta(x, y))
 391 }
 392 
 393 // Factorial calculates the product of all integers in [1, n]
 394 func Factorial(n int) int64 {
 395     return int64(math.Round(math.Gamma(float64(n + 1))))
 396 }
 397 
 398 // IsPrime checks whether an integer is bigger than 1 and can only be fully
 399 // divided by 1 and itself, which is the definition of a prime number.
 400 func IsPrime(n int64) bool {
 401     // prime numbers start at 2
 402     if n < 2 {
 403         return false
 404     }
 405 
 406     // 2 is the only even prime
 407     if n%2 == 0 {
 408         return n == 2
 409     }
 410 
 411     // no divisor can be more than the square root of the target number:
 412     // this limit makes the loop an O(sqrt(n)) one, instead of O(n); this
 413     // is a major algorithmic speedup both in theory and in practice
 414     max := int64(math.Floor(math.Sqrt(float64(n))))
 415 
 416     // the only possible full-divisors are odd integers 3..sqrt(n),
 417     // since reaching this point guarantees n is odd and n > 2
 418     for d := int64(3); d <= max; d += 2 {
 419         if n%d == 0 {
 420             return false
 421         }
 422     }
 423     return true
 424 }
 425 
 426 // LCM finds the least common-multiple of 2 positive integers; when one or
 427 // both inputs aren't positive, this func returns 0.
 428 func LCM(x, y int64) int64 {
 429     if gcd := GCD(x, y); gcd > 0 {
 430         return x * y / gcd
 431     }
 432     return 0
 433 }
 434 
 435 // GCD finds the greatest common-divisor of 2 positive integers; when one or
 436 // both inputs aren't positive, this func returns 0.
 437 func GCD(x, y int64) int64 {
 438     if x < 1 || y < 1 {
 439         return 0
 440     }
 441 
 442     // the loop below requires a >= b
 443     a, b := x, y
 444     if a < b {
 445         a, b = y, x
 446     }
 447 
 448     for b > 0 {
 449         a, b = b, a%b
 450     }
 451     return a
 452 }
 453 
 454 // Perm counts the number of all possible permutations from n objects when
 455 // picking k times. When one or both inputs aren't positive, the result is 0.
 456 func Perm(n, k int) int64 {
 457     if n < k || n < 0 || k < 0 {
 458         return 0
 459     }
 460 
 461     perm := int64(1)
 462     for i := n - k + 1; i <= n; i++ {
 463         perm *= int64(i)
 464     }
 465     return perm
 466 }
 467 
 468 // Choose counts the number of all possible combinations from n objects when
 469 // picking k times. When one or both inputs aren't positive, the result is 0.
 470 func Choose(n, k int) int64 {
 471     if n < k || n < 0 || k < 0 {
 472         return 0
 473     }
 474 
 475     // the log trick isn't always more accurate when there's no overflow:
 476     // for those cases calculate using the textbook definition
 477     f := math.Round(float64(Perm(n, k) / Factorial(k)))
 478     if !math.IsInf(f, 0) {
 479         return int64(f)
 480     }
 481 
 482     // calculate using the log-factorial of n, k, and n - k
 483     a, _ := math.Lgamma(float64(n + 1))
 484     b, _ := math.Lgamma(float64(k + 1))
 485     c, _ := math.Lgamma(float64(n - k + 1))
 486     return int64(math.Round(math.Exp(a - b - c)))
 487 }
 488 
 489 // BinomialMass calculates the probability mass of the binomial random process
 490 // given. When the probability given isn't between 0 and 1, the result is NaN.
 491 func BinomialMass(x, n int, p float64) float64 {
 492     // invalid probability input
 493     if p < 0 || p > 1 {
 494         return math.NaN()
 495     }
 496     // events outside the support are impossible by definition
 497     if x < 0 || x > n {
 498         return 0
 499     }
 500 
 501     q := 1 - p
 502     m := n - x
 503     ncx := float64(Choose(n, x))
 504     return ncx * math.Pow(p, float64(x)) * math.Pow(q, float64(m))
 505 }
 506 
 507 // CumulativeBinomialDensity calculates cumulative probabilities/masses up to
 508 // the value given for the binomial random process given. When the probability
 509 // given isn't between 0 and 1, the result is NaN.
 510 func CumulativeBinomialDensity(x, n int, p float64) float64 {
 511     // invalid probability input
 512     if p < 0 || p > 1 {
 513         return math.NaN()
 514     }
 515     if x < 0 {
 516         return 0
 517     }
 518     if x >= n {
 519         return 1
 520     }
 521 
 522     p0 := p
 523     q0 := 1 - p0
 524     q := math.Pow(q0, float64(n))
 525 
 526     pbinom := 0.0
 527     np1 := float64(n + 1)
 528     for k := 0; k < x; k++ {
 529         a, _ := math.Lgamma(np1)
 530         b, _ := math.Lgamma(float64(k + 1))
 531         c, _ := math.Lgamma(float64(n - k + 1))
 532         // count all possible combinations for this event
 533         ncomb := math.Round(math.Exp(a - b - c))
 534         pbinom += ncomb * p * q
 535         p *= p0
 536         q /= q0
 537     }
 538     return pbinom
 539 }
 540 
 541 // NormalDensity calculates the density at a point along the normal distribution
 542 // given.
 543 func NormalDensity(x float64, mu, sigma float64) float64 {
 544     z := (x - mu) / sigma
 545     return math.Sqrt(0.5/sigma) * math.Exp(-(z*z)/sigma)
 546 }
 547 
 548 // CumulativeNormalDensity calculates the probability of a normal variate of
 549 // being up to the value given.
 550 func CumulativeNormalDensity(x float64, mu, sigma float64) float64 {
 551     z := (x - mu) / sigma
 552     return 0.5 + 0.5*math.Erf(z/math.Sqrt2)
 553 }
 554 
 555 // Epanechnikov is a commonly-used kernel function.
 556 func Epanechnikov(x float64) float64 {
 557     if math.Abs(x) > 1 {
 558         // func is 0 ouside -1..+1
 559         return 0
 560     }
 561     return 0.75 * (1 - x*x)
 562 }
 563 
 564 // Gauss is the commonly-used Gaussian kernel function.
 565 func Gauss(x float64) float64 {
 566     return math.Exp(-(x * x))
 567 }
 568 
 569 // Tricube is a commonly-used kernel function.
 570 func Tricube(x float64) float64 {
 571     a := math.Abs(x)
 572     if a > 1 {
 573         // func is 0 ouside -1..+1
 574         return 0
 575     }
 576 
 577     b := a * a * a
 578     c := 1 - b
 579     return 70.0 / 81.0 * c * c * c
 580 }
 581 
 582 // SolveQuad finds the solutions of a 2nd-degree polynomial, using a formula
 583 // which is more accurate than the textbook one.
 584 func SolveQuad(a, b, c float64) (x1 float64, x2 float64) {
 585     div := 2 * c
 586     disc := math.Sqrt(b*b - 4*a*c)
 587     x1 = div / (-b - disc)
 588     x2 = div / (-b + disc)
 589     return x1, x2
 590 }

     File: ./functions_test.go
   1 package mathplus
   2 
   3 import (
   4     "math"
   5     "testing"
   6 )
   7 
   8 func TestRounding(t *testing.T) {
   9     var decadeTests = []struct {
  10         Number   float64
  11         Expected float64
  12     }{
  13         {0, 0},
  14         {-1, -10},
  15         {-0.75, -10},
  16         {-0.725, -10},
  17         {123.532123143, 120},
  18         {1932.532123143, 1930},
  19         {2023.4, 2020},
  20     }
  21 
  22     for _, tc := range decadeTests {
  23         y := Decade(tc.Number)
  24         if y != tc.Expected {
  25             const fs = `decade(%f): expected %f, got %f`
  26             t.Fatalf(fs, tc.Number, tc.Expected, y)
  27         }
  28     }
  29 
  30     var centuryTests = []struct {
  31         Number   float64
  32         Expected float64
  33     }{
  34         {0, 0},
  35         {-1, -100},
  36         {-0.75, -100},
  37         {-0.725, -100},
  38         {123.532123143, 100},
  39         {1932.532123143, 1900},
  40         {2023.4, 2000},
  41     }
  42 
  43     for _, tc := range centuryTests {
  44         y := Century(tc.Number)
  45         if y != tc.Expected {
  46             const fs = `century(%f): expected %f, got %f`
  47             t.Fatalf(fs, tc.Number, tc.Expected, y)
  48         }
  49     }
  50 
  51     var roundingTests = []struct {
  52         Number   float64
  53         Expected [6]float64
  54     }{
  55         {0, [6]float64{0, 0, 0, 0, 0, 0}},
  56         {-1, [6]float64{-1, -1, -1, -1, -1, -1}},
  57         {-0.75, [6]float64{-0.8, -0.75, -0.75, -0.75, -0.75, -0.75}},
  58         {-0.725, [6]float64{-0.7, -0.73, -0.725, -0.725, -0.725, -0.725}},
  59         {123.532123143, [6]float64{123.5, 123.53, 123.532, 123.5321, 123.53212, 123.532123}},
  60         {1932.532123143, [6]float64{1932.5, 1932.53, 1932.532, 1932.5321, 1932.53212, 1932.532123}},
  61         {2023.4, [6]float64{2023.4, 2023.4, 2023.4, 2023.4, 2023.4, 2023.4}},
  62     }
  63 
  64     for _, tc := range roundingTests {
  65         x := tc.Number
  66         y := []float64{
  67             Round1(x), Round2(x), Round3(x), Round4(x), Round5(x), Round6(x),
  68         }
  69 
  70         for i, f := range y {
  71             exp := tc.Expected[i]
  72             if math.Abs(exp-f) > 1e-12 {
  73                 const fs = `r%d(%f): expected %f, got %f`
  74                 t.Fatalf(fs, i+1, tc.Number, exp, f)
  75             }
  76         }
  77     }
  78 }
  79 
  80 var scaleTests = []struct {
  81     Input    float64
  82     InMin    float64
  83     InMax    float64
  84     OutMin   float64
  85     OutMax   float64
  86     Expected float64
  87 }{
  88     {-2, -5, 4, 0, 1, 1.0 / 3},
  89     {0.1, 0, 0.5, -3, 5, -1.4},
  90 }
  91 
  92 func TestScale(t *testing.T) {
  93     for _, tc := range scaleTests {
  94         in := tc.Input
  95         exp := tc.Expected
  96         got := Scale(in, tc.InMin, tc.InMax, tc.OutMin, tc.OutMax)
  97         if got != exp {
  98             const fs = `Scale(%f%f%f%f%f): expected %f, got %f`
  99             t.Fatalf(fs, in, tc.InMin, tc.InMax, tc.OutMin, tc.OutMax, exp, got)
 100         }
 101     }
 102 }
 103 
 104 func TestIsPrime(t *testing.T) {
 105     var tests = []struct {
 106         Input    int64
 107         Expected bool
 108     }{
 109         {-3, false},
 110         {0, false},
 111         {1, false},
 112         {4, false},
 113         {9, false},
 114         {21, false},
 115 
 116         {2, true},
 117         {3, true},
 118         {5, true},
 119         {19, true},
 120         // 15,485,863 is the millionth prime
 121         {15_485_863, true},
 122     }
 123 
 124     for _, tc := range tests {
 125         if v := IsPrime(tc.Input); v != tc.Expected {
 126             const fs = `isprime(%d) wrongly returned %v`
 127             t.Fatalf(fs, tc.Input, v)
 128         }
 129     }
 130 }
 131 
 132 func TestHorner(t *testing.T) {
 133     var tests = []struct {
 134         X        float64
 135         C        []float64
 136         Expected float64
 137     }{
 138         {2, []float64{1, 2, 3}, 11},
 139         {3, []float64{3, 5, -1}, 41},
 140     }
 141 
 142     for _, tc := range tests {
 143         got := Polyval(tc.X, tc.C...)
 144         if got != tc.Expected {
 145             const fs = `horner(%f%#v) gave %f, instead of %f`
 146             t.Fatalf(fs, tc.X, tc.C, got, tc.Expected)
 147             return
 148         }
 149     }
 150 }
 151 
 152 func TestGCD(t *testing.T) {
 153     var tests = []struct {
 154         X        int64
 155         Y        int64
 156         Expected int64
 157     }{
 158         {0, 0, 0},
 159         {-1, 10, 0},
 160         {1, -10, 0},
 161         {1, 1, 1},
 162         {1, 7, 1},
 163         {3 * 12, 12, 12},
 164         {1280, 1920, 640},
 165     }
 166 
 167     for _, tc := range tests {
 168         got := GCD(tc.X, tc.Y)
 169         if got != tc.Expected {
 170             const fs = `gcd(%d%d) gave %d, instead of %d`
 171             t.Fatalf(fs, tc.X, tc.Y, got, tc.Expected)
 172             return
 173         }
 174     }
 175 }
 176 
 177 func TestPerm(t *testing.T) {
 178     var tests = []struct {
 179         X        int
 180         Y        int
 181         Expected int64
 182     }{
 183         {10, 4, 5_040},
 184         {5, 0, 1},
 185         {5, 5, 120},
 186     }
 187 
 188     for _, tc := range tests {
 189         got := Perm(tc.X, tc.Y)
 190         if got != tc.Expected {
 191             const fs = `perm(%d%d) gave %d, instead of %d`
 192             t.Fatalf(fs, tc.X, tc.Y, got, tc.Expected)
 193             return
 194         }
 195     }
 196 }
 197 
 198 func TestChoose(t *testing.T) {
 199     var tests = []struct {
 200         X        int
 201         Y        int
 202         Expected int64
 203     }{
 204         {10, 4, 210},
 205         {10, 0, 1},
 206         {10, 10, 1},
 207     }
 208 
 209     for _, tc := range tests {
 210         got := Choose(tc.X, tc.Y)
 211         if got != tc.Expected {
 212             const fs = `comb(%d%d) gave %d, instead of %d`
 213             t.Fatalf(fs, tc.X, tc.Y, got, tc.Expected)
 214             return
 215         }
 216     }
 217 }

     File: ./go.mod
   1 module mathplus
   2 
   3 go 1.16

     File: ./integers.go
   1 package mathplus
   2 
   3 import "math"
   4 
   5 const (
   6     MinInt16 = -(1 << 15)
   7     MinInt32 = -(1 << 31)
   8     MinInt64 = -(1 << 63)
   9 
  10     MaxInt16 = 1<<15 - 1
  11     MaxInt32 = 1<<31 - 1
  12     MaxInt64 = 1<<63 - 1
  13 
  14     MaxUint16 = 1<<16 - 1
  15     MaxUint32 = 1<<32 - 1
  16     MaxUint64 = 1<<64 - 1
  17 )
  18 
  19 func CountIntegerDigits(n int64) int {
  20     if n < 0 {
  21         n = -n
  22     }
  23 
  24     // 0 doesn't have a log10
  25     if n == 0 {
  26         return 1
  27     }
  28     // add 1 to the floored logarithm, since digits are always full and
  29     // with ceiling-like behavior, to use an analogy for this formula
  30     return int(math.Floor(math.Log10(float64(n)))) + 1
  31 }
  32 
  33 func LoopThousandsGroups(n int64, fn func(i, n int)) {
  34     // 0 doesn't have a log10
  35     if n == 0 {
  36         fn(0, 0)
  37         return
  38     }
  39 
  40     sign := +1
  41     if n < 0 {
  42         n = -n
  43         sign = -1
  44     }
  45 
  46     intLog1000 := int(math.Log10(float64(n)) / 3)
  47     remBase := int64(math.Pow10(3 * intLog1000))
  48 
  49     for i := 0; remBase > 0; i++ {
  50         group := (1000 * n) / remBase / 1000
  51         fn(i, sign*int(group))
  52         // if original number was negative, ensure only first
  53         // group gives a negative input to the callback
  54         sign = +1
  55 
  56         n %= remBase
  57         remBase /= 1000
  58     }
  59 }
  60 
  61 var pow10 = []int64{
  62     1,
  63     10,
  64     100,
  65     1000,
  66     10000,
  67     100000,
  68     1000000,
  69     10000000,
  70     100000000,
  71     1000000000, // last entry for int32
  72     10000000000,
  73     100000000000,
  74     1000000000000,
  75     10000000000000,
  76     100000000000000,
  77     1000000000000000,
  78     10000000000000000,
  79     100000000000000000,
  80     1000000000000000000,
  81     // 10000000000000000000,
  82 }
  83 
  84 var pow2 = []int64{
  85     1,
  86     2,
  87     4,
  88     8,
  89     16,
  90     32,
  91     64,
  92     128,
  93     256,
  94     512,
  95     1024,
  96     2048,
  97     4096,
  98     8192,
  99     16384,
 100     32768,
 101     65536,
 102     131072,
 103     262144,
 104     524288,
 105     1048576,
 106     2097152,
 107     4194304,
 108     8388608,
 109     16777216,
 110     33554432,
 111     67108864,
 112     134217728,
 113     268435456,
 114     536870912,
 115     1073741824,
 116     2147483648, // last entry for int32
 117     4294967296,
 118     8589934592,
 119     17179869184,
 120     34359738368,
 121     68719476736,
 122     137438953472,
 123     274877906944,
 124     549755813888,
 125     1099511627776,
 126     2199023255552,
 127     4398046511104,
 128     8796093022208,
 129     17592186044416,
 130     35184372088832,
 131     70368744177664,
 132     140737488355328,
 133     281474976710656,
 134     562949953421312,
 135     1125899906842624,
 136     2251799813685248,
 137     4503599627370496,
 138     9007199254740992,
 139     18014398509481984,
 140     36028797018963968,
 141     72057594037927936,
 142     144115188075855872,
 143     288230376151711744,
 144     576460752303423488,
 145     1152921504606846976,
 146     2305843009213693952,
 147     4611686018427387904,
 148     // 9223372036854775808,
 149 }
 150 
 151 // Log2Int gives you the floor of the base-2 logarithm, or negative results
 152 // for invalid inputs. Values less than 1 aren't supported, since they don't
 153 // have logarithms of any valid base, let alone base-2 logarithms.
 154 func Log2Int(n int64) (log2 int, ok bool) {
 155     if n < 1 {
 156         return -1, false
 157     }
 158 
 159     v := uint64(n)
 160     mask := uint64(1 << 63) // 2**63
 161     for i := int64(63); i > 0; i-- {
 162         if v&mask != 0 {
 163             return int(i), true
 164         }
 165         mask >>= 1
 166     }
 167     return 0, true
 168 }
 169 
 170 // Log10Int gives you the floor of the base-10 logarithm, or negative results
 171 // for invalid inputs. Values less than 1 aren't supported, since they don't
 172 // have logarithms of any valid base, let alone base-10 logarithms.
 173 func Log10Int(n int64) (log10 int, ok bool) {
 174     if n < 1 {
 175         return -1, false
 176     }
 177 
 178     for i := int64(0); i < 19; i++ {
 179         n /= 10
 180         if n == 0 {
 181             return int(i), true
 182         }
 183     }
 184     return 0, true
 185 }
 186 
 187 // Pow10Int gives you the integers powers of 10 as far as int64 allows: negative
 188 // inputs aren't supported, since negative exponents don't give integer results.
 189 func Pow10Int(n int) (power10 int64, ok bool) {
 190     if 0 <= n && n < len(pow10) {
 191         return pow10[n], true
 192     }
 193     return -1, false
 194 }
 195 
 196 // Pow2Int gives you the integers powers of 2 as far as int64 allows: negative
 197 // inputs aren't supported, since negative exponents don't give integer results.
 198 func Pow2Int(n int) (power2 int64, ok bool) {
 199     if 0 <= n && n < len(pow2) {
 200         return pow2[n], true
 201     }
 202     return -1, false
 203 }

     File: ./integers_test.go
   1 package mathplus
   2 
   3 import "testing"
   4 
   5 func TestCountIntegerDigits(t *testing.T) {
   6     var tests = []struct {
   7         Input    int64
   8         Expected int
   9     }{
  10         {0, 1},
  11         {-32, 2},
  12         {999, 3},
  13         {3_490, 4},
  14         {12_332, 5},
  15         {999_999, 6},
  16         {1_000_000, 7},
  17         {1_000_001, 7},
  18         {12_345_678, 8},
  19     }
  20 
  21     for _, tc := range tests {
  22         if n := CountIntegerDigits(tc.Input); n != tc.Expected {
  23             const fs = `integer digits in %d: got %d instead of %d`
  24             t.Errorf(fs, tc.Input, n, tc.Expected)
  25         }
  26     }
  27 }
  28 
  29 func TestLoopThousandsGroups(t *testing.T) {
  30     var tests = []struct {
  31         Input    int64
  32         Expected []int
  33     }{
  34         {0, []int{0}},
  35         {-32, []int{-32}}, // negatives not supported yet
  36         {999, []int{999}},
  37         {1_670, []int{1, 670}},
  38         {3_490, []int{3, 490}},
  39         {12_332, []int{12, 332}},
  40         {999_999, []int{999, 999}},
  41         {1_000_000, []int{1, 0, 0}},
  42         {1_000_001, []int{1, 0, 1}},
  43         {1_234_567, []int{1, 234, 567}},
  44     }
  45 
  46     // return
  47     for _, tc := range tests {
  48         count := 0
  49         LoopThousandsGroups(tc.Input, func(i, n int) {
  50             // t.Log(tc.Input, i, n)
  51             if n != tc.Expected[i] {
  52                 const fs = `group %d in %d: got %d instead of %d`
  53                 t.Errorf(fs, i, tc.Input, n, tc.Expected[i])
  54             }
  55             count++
  56         })
  57 
  58         if count != len(tc.Expected) {
  59             const fs = `thousands-groups from %d: got %d instead of %d`
  60             t.Errorf(fs, tc.Input, count, len(tc.Expected))
  61         }
  62     }
  63 }
  64 
  65 func TestLog2Int(t *testing.T) {
  66     var tests = []struct {
  67         Value    int64
  68         Expected int
  69     }{
  70         {-3, -1},
  71         {1, 0},
  72         {2, 1},
  73         {3, 1},
  74         {4, 2},
  75         {1024, 10},
  76         {1_025, 10},
  77         {2*1024 - 1, 10},
  78     }
  79 
  80     for _, tc := range tests {
  81         got, ok := Log2Int(tc.Value)
  82         if got != tc.Expected || (ok && tc.Value < 1) {
  83             const fs = `log2int(%d) = %d, but got %d instead`
  84             t.Fatalf(fs, tc.Value, tc.Expected, got)
  85         }
  86     }
  87 }
  88 
  89 func TestLog10Int(t *testing.T) {
  90     var tests = []struct {
  91         Value    int64
  92         Expected int
  93     }{
  94         {-3, -1},
  95         {1, 0},
  96         {10, 1},
  97         {100, 2},
  98         {101, 2},
  99         {199, 2},
 100         {1_000_000, 6},
 101     }
 102 
 103     for _, tc := range tests {
 104         got, ok := Log10Int(tc.Value)
 105         if got != tc.Expected || (ok && tc.Value < 1) {
 106             const fs = `log10int(%d) = %d, but got %d instead`
 107             t.Fatalf(fs, tc.Value, tc.Expected, got)
 108         }
 109     }
 110 }

     File: ./numbers.go
   1 package mathplus
   2 
   3 import (
   4     "math"
   5 )
   6 
   7 const (
   8     // the maximum integer a float64 can represent exactly; since float64
   9     // uses a sign bit, instead of a 2-complement mantissa, -maxflint is
  10     // the minimum integer
  11     maxflint = 2 << 52
  12 )
  13 
  14 // StringWidth counts how many runes it takes to represent the value given as a
  15 // string both as a plain no-commas string and as a string with digit-grouping
  16 // commas, if needed. Each group is 3 digits. Showing numbers with commas makes
  17 // long numbers easier to read.
  18 func StringWidth(f float64, decimals int) (plain, nice int) {
  19     if math.IsNaN(f) || math.IsInf(f, 0) {
  20         return 0, 0
  21     }
  22 
  23     // avoid wrong results, since decimals will be added at the end
  24     if decimals < 0 {
  25         decimals = 0
  26     }
  27 
  28     extras := 0 // count non-digits, such as negatives and dots
  29     if f < 0 {
  30         f = -f   // fix value for uintWidth
  31         extras++ // count the leading negative sign
  32     }
  33     if decimals > 0 {
  34         extras++ // count the decimal dot
  35     }
  36 
  37     // at this point, f >= 0 for sure
  38     plain, nice = uintWidth(f)
  39     plain += extras + decimals
  40     nice += extras + decimals
  41     return plain, nice
  42 }
  43 
  44 // uintWidth can't handle negatives correctly, as its name suggests
  45 func uintWidth(f float64) (plain, nice int) {
  46     // only 1 digit and 0 commas for 0 <= x < 10
  47     if f < 10 {
  48         return 1, 1
  49     }
  50 
  51     mag := math.Log10(math.Floor(f)) // order of magnitude
  52     digits := int(mag) + 1           // integer digits
  53     commas := int(mag) / 3           // commas separating digits in groups of 3
  54     return digits, digits + commas
  55 }
  56 
  57 // Equals allows easily comparing numbers, including NaN values, which otherwise
  58 // never equal anything they're compared to, including themselves.
  59 func Equals(x, y float64) bool {
  60     if x == y {
  61         return true
  62     }
  63     return math.IsNaN(x) && math.IsNaN(y)
  64 }
  65 
  66 // ApproxEqual allows approximate comparisons, as well as checking if both values
  67 // are NaN, which otherwise never equal themselves when compared directly. The
  68 // last parameter must be positive for approx. comparisons, or zero for exact
  69 // comparisons.
  70 func ApproxEqual(x, y float64, maxdiff float64) bool {
  71     if math.Abs(x-y) <= maxdiff {
  72         return true
  73     }
  74     return math.IsNaN(x) && math.IsNaN(y)
  75 }
  76 
  77 // GuessDecimalsCount tries to guess the number of decimal digits of the number
  78 // given.
  79 func GuessDecimalsCount(x float64, max int) int {
  80     if x < 0 {
  81         x = -x
  82     }
  83 
  84     // only up to 16, because log10(2**-53) ~= -15.9546
  85     if !(0 <= max && max <= 16) {
  86         max = 16
  87     }
  88 
  89     const tol = 5e-13 // 1e-11, 1e-12, 5e-13
  90     for digits := 0; digits <= max; digits++ {
  91         _, frac := math.Modf(x)
  92         frac = math.Abs(frac)
  93         // when it's time to stop, the absolute value of the fraction part is
  94         // either extremely close to 0 or extremely close to 1
  95         if frac < tol || (1-frac) < tol {
  96             return digits
  97         }
  98         x *= 10
  99     }
 100     return max
 101 }
 102 
 103 // Default returns the first non-NaN value among those given: failing that,
 104 // the result will be NaN.
 105 func Default(args ...float64) float64 {
 106     for _, x := range args {
 107         if !math.IsNaN(x) {
 108             return x
 109         }
 110     }
 111     return math.NaN()
 112 }
 113 
 114 // TrimSlice ignores all leading/trailing NaN values from the slice given: its
 115 // main use-case is after sorting via sort.Float64s, since all NaN values are
 116 // moved in place to the start/end of the now-sorted slice.
 117 //
 118 // # Example
 119 //
 120 // sort.Float64(values)
 121 // res := mathplus.TrimSlice(values)
 122 func TrimSlice(x []float64) []float64 {
 123     // ignore all leading NaNs
 124     for len(x) > 0 && math.IsNaN(x[0]) {
 125         x = x[1:]
 126     }
 127     // ignore all trailing NaNs
 128     for len(x) > 0 && math.IsNaN(x[len(x)-1]) {
 129         x = x[:len(x)-1]
 130     }
 131     return x
 132 }
 133 
 134 // Linspace works like the Matlab function of the same name, except it takes a
 135 // callback, generalizing its behavior.
 136 func Linspace(a, incr, b float64, f func(x float64)) {
 137     if incr <= 0 {
 138         return
 139     }
 140 
 141     if a < b {
 142         forwardLinspace(a, incr, b, f)
 143     } else if a > b {
 144         backwardLinspace(a, incr, b, f)
 145     } else {
 146         f(a)
 147     }
 148 }
 149 
 150 func forwardLinspace(a, incr, b float64, f func(x float64)) {
 151     for i := 0; true; i++ {
 152         x := float64(i)*incr + a
 153         if x <= b {
 154             f(x)
 155         } else {
 156             return
 157         }
 158     }
 159 }
 160 
 161 func backwardLinspace(a, incr, b float64, f func(x float64)) {
 162     for i := 0; true; i++ {
 163         x := b - float64(i)*incr
 164         if x >= a {
 165             f(x)
 166         } else {
 167             return
 168         }
 169     }
 170 }
 171 
 172 // Seq is a special case of Linspace, where the increment is +/-1.
 173 func Seq(a, b float64, f func(x float64)) {
 174     Linspace(a, 1, b, f)
 175 }
 176 
 177 // Increment increments the number given by adding 1, when it's in the int-safe
 178 // range of float64s; when outside that range, the smallest available delta is
 179 // used instead.
 180 func Increment(x float64) float64 {
 181     if math.IsNaN(x) || math.IsInf(x, 0) {
 182         return x
 183     }
 184     if incr := x + 1; incr != x {
 185         return incr
 186     }
 187     return math.Float64frombits(math.Float64bits(x) + 1)
 188 }
 189 
 190 // Decrement decrements the number given by subtracting 1, when it's in the
 191 // int-safe range of float64s: when outside that range, the smallest available
 192 // delta is used instead.
 193 func Decrement(x float64) float64 {
 194     if math.IsNaN(x) || math.IsInf(x, 0) {
 195         return x
 196     }
 197     if decr := x - 1; decr != x {
 198         return decr
 199     }
 200     // subtracting 1 from the value's uint64 counterpart doesn't work for values
 201     // less than the minimum exact-integer: uint64s use 2-complement arithmetic
 202     return -math.Float64frombits(math.Float64bits(-x) + 1)
 203 }
 204 
 205 func Deapproximate(x float64) (min float64, max float64) {
 206     if math.IsNaN(x) || math.IsInf(x, 0) {
 207         return x, x
 208     }
 209 
 210     if x == 0 {
 211         return -0.5, +0.5
 212     }
 213 
 214     if math.Remainder(x, 1) == 0 {
 215         if math.Remainder(x, 10) != 0 {
 216             return x - 0.5, x + 0.5
 217         }
 218 
 219         sign := Sign(x)
 220         abs := math.Abs(x)
 221         delta := float64(maxExactPow10(int(abs))) / 2
 222         return sign*abs - delta, sign*abs + delta
 223     }
 224 
 225     // return surrounding integers when given non-integers
 226     return math.Floor(x), math.Ceil(x)
 227 }
 228 
 229 func maxExactPow10(x int) int {
 230     pow10 := 1
 231     for {
 232         if x%pow10 != 0 {
 233             return pow10 / 10
 234         }
 235         pow10 *= 10
 236     }
 237 }

     File: ./numbers_test.go
   1 package mathplus
   2 
   3 import (
   4     "math"
   5     "strconv"
   6     "testing"
   7 )
   8 
   9 func TestNumberWidth(t *testing.T) {
  10     var tests = []struct {
  11         Input       float64
  12         Decimals    int
  13         PlainWidth  int
  14         PrettyWidth int
  15     }{
  16         {0, 0, len(`0`), len(`0`)},
  17         {0.923123, 0, len(`0`), len(`0`)},
  18         {0.923123, 2, len(`0.92`), len(`0.92`)},
  19         {-3, 0, len(`-3`), len(`-3`)},
  20         {-300, 0, len(`-300`), len(`-300`)},
  21         {+2_300, 0, len(`2300`), len(`2,300`)},
  22         {-1_000_000, 0, len(`-1000000`), len(`-1,000,000`)},
  23         {+1_000_000, 0, len(`1000000`), len(`1,000,000`)},
  24         {+1_610_058.23423, 4, len(`1610058.2342`), len(`1_610_058.2342`)},
  25         {-317_289, 4, len(`-317289.0000`), len(`-317_289.0000`)},
  26         {5_457_013.0000, 4, len(`5457013.0000`), len(`5_457_013.0000`)},
  27         {-118.406800000000004, 15, 20, 20},
  28     }
  29 
  30     for _, tc := range tests {
  31         pl, pr := StringWidth(tc.Input, tc.Decimals)
  32         if pl != tc.PlainWidth {
  33             const fs = `PlainWidth(%f%d): expected %d, got %d`
  34             t.Fatalf(fs, tc.Input, tc.Decimals, tc.PlainWidth, pl)
  35         }
  36         if pr != tc.PrettyWidth {
  37             const fs = `PrettyWidth(%f%d): expected %d, got %d`
  38             t.Fatalf(fs, tc.Input, tc.Decimals, tc.PrettyWidth, pr)
  39         }
  40     }
  41 }
  42 
  43 func TestNumberWidthInternal(t *testing.T) {
  44     var tests = []struct {
  45         Number    float64
  46         Decimals  int
  47         IntDigits int
  48         Commas    int
  49     }{
  50         {0, 0, 1, 0},
  51         {-3, 0, 1, 0},
  52         {-300, 0, 3, 0},
  53         {+2_300, 0, 4, 1},
  54         {-1_000_000, 0, 7, 2},
  55         {+1_000_000, 0, 7, 2},
  56         {+1_610_058.23423, 4, 7, 2},
  57         {-317_289, 4, 6, 1},
  58         {5_457_013.0000, 4, 7, 2},
  59     }
  60 
  61     for _, tc := range tests {
  62         // remember that uintWidth can't handle negatives
  63         digits, pretty := uintWidth(math.Abs(tc.Number))
  64 
  65         commas := pretty - digits
  66         if digits != tc.IntDigits || commas != tc.Commas {
  67             const fs = `%f: expected %d and %d, but got %d and %d instead`
  68             t.Fatalf(fs, tc.Number, tc.IntDigits, tc.Commas, digits, commas)
  69         }
  70     }
  71 }
  72 
  73 func TestGuessDecimalsCount(t *testing.T) {
  74     var tests = []struct {
  75         Input    float64
  76         Expected int
  77     }{
  78         {-1, 0},
  79         {-1.3, 1},
  80         {-1.1, 1},
  81         {10.103, 3},
  82         {3.9500, 2},
  83         {1.789000, 3},
  84         {5.0409999000000000, 7},
  85     }
  86 
  87     for _, tc := range tests {
  88         got := GuessDecimalsCount(tc.Input, -1)
  89         if got != tc.Expected {
  90             const fs = `guess(%f%d) gave %d, but should have been %d`
  91             t.Fatalf(fs, tc.Input, -1, got, tc.Expected)
  92         }
  93     }
  94 
  95     var v = []float64{
  96         1.773, 1.789, 1.773, 1.779, 1.817, 1.797, 1.780, 1.737,
  97         1.723, 1.725, 1.733, 1.742, 1.746, 1.728, 1.726, 1.701,
  98         1.690, 1.675, 1.680, 1.672, 1.685, 1.679, 1.674, 1.653,
  99         1.668, 1.669, 1.680, 1.693, 1.713, 1.712, 1.733, 1.744,
 100         1.739, 1.708, 1.695, 1.691, 1.715, 1.733, 1.730, 1.722,
 101         1.737, 1.749, 1.774, 1.773, 1.789, 1.804, 1.795, 1.752,
 102         1.770, 1.736, 1.717, 1.690, 1.696, 1.703, 1.679, 1.673,
 103         1.694, 1.710, 1.711, 1.739, 1.725, 1.748, 1.764, 1.781,
 104         1.766, 1.741, 1.739, 1.750, 1.752, 1.731, 1.717, 1.699,
 105         1.689, 1.701, 1.713, 1.715, 1.713, 1.780, 1.795, 1.822,
 106         1.832, 1.835, 1.846, 1.857, 1.869, 1.847, 1.831, 1.874,
 107         1.928, 1.900, 1.925, 1.888, 1.842, 1.836, 1.825, 1.804,
 108         1.745, 1.753, 1.779, 1.797, 1.812, 1.820, 1.829, 1.818,
 109         1.799, 1.795, 1.795, 1.792, 1.792, 1.760, 1.741, 1.721,
 110         1.718, 1.734, 1.751, 1.776, 1.768, 1.774, 1.790, 1.813,
 111         1.840, 1.861, 1.853, 1.834, 1.848, 1.859, 1.856, 1.839,
 112         1.826, 1.833, 1.850, 1.861, 1.874, 1.905, 1.931, 1.951,
 113         1.971, 1.985, 1.997, 2.023, 2.047, 2.080, 2.103, 2.139,
 114         2.155, 2.151, 2.142, 2.147, 2.147, 2.165, 2.171, 2.176,
 115         2.174, 2.190, 2.192, 2.184, 2.170, 2.155, 2.155, 2.140,
 116         2.129, 2.124, 2.128, 2.146, 2.135, 2.151, 2.141, 2.134,
 117         2.086, 2.059, 2.014, 1.980, 1.969, 1.972, 1.953, 1.927,
 118         1.915, 1.901, 1.905, 1.941, 1.951, 1.940, 1.949, 1.959,
 119         1.985, 1.995, 2.001, 2.002, 2.017, 2.018, 2.006, 1.997,
 120         1.994, 1.987, 1.978, 1.977, 1.940, 1.928, 1.909, 1.813,
 121         1.736, 1.710, 1.712, 1.726, 1.740, 1.747, 1.748, 1.745,
 122         1.728, 1.714, 1.656, 1.649, 1.642, 1.630, 1.612, 1.588,
 123         1.583, 1.570, 1.568, 1.561, 1.537, 1.542, 1.548, 1.538,
 124         1.520, 1.521, 1.513, 1.494, 1.469, 1.458, 1.449, 1.441,
 125         1.434, 1.428, 1.431, 1.440, 1.447, 1.457, 1.456, 1.457,
 126         1.450, 1.451, 1.439, 1.423, 1.374, 1.100, 1.147, 1.134,
 127         1.121, 1.119, 1.115,
 128     }
 129 
 130     for i, x := range v {
 131         guess := gdc(t, x)
 132         // guess := GuessDecimalsCount(v, -1)
 133         if guess > 3 {
 134             _, diff := math.Modf(1000 * x)
 135             t.Logf("len: %d\n", len(v))
 136             const fs = `(item %2d%f doesn't have %d decimals; diff: %.20f`
 137             t.Fatalf(fs, i+1, x, guess, diff)
 138         }
 139     }
 140 
 141     v = []float64{
 142         6000.00, 2300.00, 2000.00, 1700.00, 1500.00, 1500.00, 1400.00, 1300.00,
 143         900.00, 800.00, 800.00, 700.00, 600.00, 550.00, 550.00, 400.00,
 144         320.39, 300.00, 200.00, 198.56, 155.94, 155.73, 98.15, 80.00,
 145         80.00, 74.50, 70.00, 70.00, 70.00, 70.00, 70.00, 65.00,
 146         60.00, 55.45, 49.00, 35.00, 35.00, 25.00, 25.00, 25.00,
 147         20.00, 15.00, 15.00, 12.00, 10.00, 10.00, 10.00, 10.00,
 148         7.00, 6.00, 6.00, 5.00, 3.00, 2.00, 1.00,
 149     }
 150     for i, x := range v {
 151         guess := gdc(t, x)
 152         // guess := GuessDecimalsCount(v, -1)
 153         if guess > 2 {
 154             _, diff := math.Modf(100 * x)
 155             t.Logf("len: %d\n", len(v))
 156             const fs = `(item %2d%f doesn't have %d decimals; diff: %.20f`
 157             t.Fatalf(fs, i+1, x, guess, diff)
 158         }
 159     }
 160 }
 161 
 162 // gdc logs each step of the loop inside GuessDecimalsCount
 163 func gdc(t *testing.T, x float64) int {
 164     const max = 16
 165     const tol = 5e-11 // 5e-13
 166     for digits := 0; digits <= max; digits++ {
 167         _, frac := math.Modf(x)
 168         const fs = "x: %f\tdigits: %d\tdiff: %.20f\tcdiff: %.20f\n"
 169         t.Logf(fs, x, digits, frac, 1-frac)
 170         ar := math.Abs(frac)
 171         if ar < tol || (1-ar) < tol {
 172             return digits
 173         }
 174         x *= 10
 175     }
 176     return max
 177 }
 178 
 179 func TestIncrement(t *testing.T) {
 180     var tests = []struct {
 181         Input    float64
 182         Expected float64
 183     }{
 184         {math.NaN(), math.NaN()},
 185         {math.Inf(-1), math.Inf(-1)},
 186         {math.Inf(+1), math.Inf(+1)},
 187 
 188         {-maxflint - 2, -maxflint},
 189         {-maxflint, -maxflint + 1},
 190         {-1, 0},
 191         {0, 1},
 192         {10.25, 11.25},
 193         {maxflint - 1, maxflint},
 194         {maxflint, maxflint + 2},
 195     }
 196 
 197     for _, tc := range tests {
 198         name := strconv.FormatFloat(tc.Input, 'f', 2, 64)
 199         t.Run(name, func(t *testing.T) {
 200             got := Increment(tc.Input)
 201             if !Equals(got, tc.Expected) {
 202                 const fs = `got %v, instead of %v`
 203                 t.Fatalf(fs, got, tc.Expected)
 204             }
 205         })
 206     }
 207 }
 208 
 209 func TestDecrement(t *testing.T) {
 210     var tests = []struct {
 211         Input    float64
 212         Expected float64
 213     }{
 214         {math.NaN(), math.NaN()},
 215         {math.Inf(-1), math.Inf(-1)},
 216         {math.Inf(+1), math.Inf(+1)},
 217 
 218         {-maxflint, -maxflint - 2},
 219         {-maxflint + 1, -maxflint},
 220         {-1, -2},
 221         {0, -1},
 222         {11.25, 10.25},
 223         {maxflint, maxflint - 1},
 224         {maxflint + 2, maxflint},
 225     }
 226 
 227     for _, tc := range tests {
 228         name := strconv.FormatFloat(tc.Input, 'f', 2, 64)
 229         t.Run(name, func(t *testing.T) {
 230             got := Decrement(tc.Input)
 231             if !Equals(got, tc.Expected) {
 232                 const fs = `got %v, instead of %v`
 233                 t.Fatalf(fs, got, tc.Expected)
 234             }
 235         })
 236     }
 237 }
 238 
 239 func TestDeapproximate(t *testing.T) {
 240     var tests = []struct {
 241         Input       float64
 242         ExpectedMin float64
 243         ExpectedMax float64
 244     }{
 245         {0, -0.5, +0.5},
 246         {1, 0.5, 1.5},
 247         {10, 5, 15},
 248         {100, 50, 150},
 249         {101, 100.5, 101.5},
 250         {102, 101.5, 102.5},
 251         {120, 115, 125},
 252     }
 253 
 254     for _, tc := range tests {
 255         name := strconv.FormatFloat(tc.Input, 'f', 6, 64)
 256         t.Run(name, func(t *testing.T) {
 257             min, max := Deapproximate(tc.Input)
 258             if !Equals(min, tc.ExpectedMin) || !Equals(max, tc.ExpectedMax) {
 259                 const fs = `got %f and %f, instead of %f and %f`
 260                 t.Fatalf(fs, min, max, tc.ExpectedMin, tc.ExpectedMax)
 261             }
 262         })
 263     }
 264 }
 265 
 266 func TestMaxExactPow10(t *testing.T) {
 267     var tests = []struct {
 268         Input    int
 269         Expected int
 270     }{
 271         // {0, 1},
 272         {1, 1},
 273         {4, 1},
 274         {10, 10},
 275         {100, 100},
 276         {101, 1},
 277         {102, 1},
 278         {120, 10},
 279     }
 280 
 281     for _, tc := range tests {
 282         name := strconv.Itoa(tc.Input)
 283         t.Run(name, func(t *testing.T) {
 284             got := maxExactPow10(int(tc.Input))
 285             if got != int(tc.Expected) {
 286                 const fs = `got %d, instead of %d`
 287                 t.Fatalf(fs, got, tc.Expected)
 288             }
 289         })
 290     }
 291 }

     File: ./statistics.go
   1 package mathplus
   2 
   3 import (
   4     "math"
   5 )
   6 
   7 // Quantile calculates a sorted array's quantile, parameterized by a number in
   8 // [0, 1]; the sorted array must be in increasing order and can't contain any
   9 // NaN value.
  10 func Quantile(x []float64, q float64) float64 {
  11     l := len(x)
  12     // quantiles aren't defined for empty arrays/samples
  13     if l == 0 {
  14         return math.NaN()
  15     }
  16 
  17     // calculate indices of surrounding values, which match when index isn't
  18     // fractional
  19     mid := float64(l-1) * q
  20     low := math.Floor(mid)
  21     high := math.Ceil(mid)
  22 
  23     // for fractional indices, interpolate their 2 surrounding values
  24     a := x[int(low)]
  25     b := x[int(high)]
  26     return a + (mid-low)*(b-a)
  27 }
  28 
  29 // welford has almost everything needed to implement Welford's running-stats
  30 // algorithm for the arithmetic mean and the standard deviation: the only
  31 // thing missing is the count of values so far, which you must provide every
  32 // time you call its methods.
  33 type welford struct {
  34     Mean          float64
  35     meanOfSquares float64
  36 }
  37 
  38 // Update advances Welford's algorithm, using the value and item-count given.
  39 func (w *welford) Update(x float64, n int) {
  40     d1 := x - w.Mean
  41     w.Mean += d1 / float64(n)
  42     d2 := x - w.Mean
  43     w.meanOfSquares += d1 * d2
  44 }
  45 
  46 // SD calculates the current standard-deviation. The first parameter is the
  47 // bias-correction to use: 0 gives you the `population` SD, 1 gives you the
  48 // `sample` SD; 1.5 is very-rarely used and only in special circumstances.
  49 func (w welford) SD(bias float64, n int) float64 {
  50     denom := float64(n) - bias
  51     return math.Sqrt(w.meanOfSquares / denom)
  52 }
  53 
  54 // RMS calculates the current root-mean-square statistic.
  55 func (w welford) RMS() float64 {
  56     return math.Sqrt(w.meanOfSquares)
  57 }
  58 
  59 // NumberSummary has/updates numeric constant-space running stats.
  60 type NumberSummary struct {
  61     // struct welford has the Mean public field
  62     welford
  63 
  64     // Count is how many values there are so far, including NaNs
  65     Count int
  66 
  67     // NaN counts NaN values so far
  68     NaN int
  69 
  70     // Integers counts all integers so far
  71     Integers int
  72 
  73     // Negatives counts all negative number so far
  74     Negatives int
  75 
  76     // Zeros counts all zeros so far
  77     Zeros int
  78 
  79     // Positives counts all positive numbers so far
  80     Positives int
  81 
  82     // Min is the least number so far, ignoring NaNs
  83     Min float64
  84 
  85     // Max is the highest number so far, ignoring NaNs
  86     Max float64
  87 
  88     // Sum is the sum of all numbers so far, ignoring NaNs
  89     Sum float64
  90 
  91     // sumOfLogs is used to calculate the geometric mean
  92     sumOfLogs float64
  93 }
  94 
  95 // Update does exactly what it says.
  96 func (ns *NumberSummary) Update(f float64) {
  97     if ns.Count == 0 {
  98         ns.Min = math.Inf(+1)
  99         ns.Max = math.Inf(-1)
 100     }
 101 
 102     ns.Count++
 103     if math.IsNaN(f) {
 104         ns.NaN++
 105         return
 106     }
 107 
 108     if _, r := math.Modf(f); r == 0 {
 109         ns.Integers++
 110     }
 111 
 112     if f > 0 {
 113         ns.Positives++
 114     } else if f == 0 {
 115         ns.Zeros++
 116     } else if f < 0 {
 117         ns.Negatives++
 118     }
 119 
 120     ns.Sum += f
 121     ns.sumOfLogs += math.Log(f)
 122     ns.Min = math.Min(ns.Min, f)
 123     ns.Max = math.Max(ns.Max, f)
 124     ns.welford.Update(f, ns.Valid())
 125 }
 126 
 127 // Valid finds how many numbers are valid so far
 128 func (ns NumberSummary) Valid() int {
 129     return ns.Count - ns.NaN
 130 }
 131 
 132 // Invalid finds how many numbers are invalid so far
 133 func (ns NumberSummary) Invalid() int {
 134     return ns.NaN
 135 }
 136 
 137 // Geomean calculates the current geometric mean
 138 func (ns NumberSummary) Geomean() float64 {
 139     if ns.Negatives > 0 || ns.Zeros > 0 {
 140         return math.NaN()
 141     }
 142     return math.Exp(ns.sumOfLogs / float64(ns.Valid()))
 143 }
 144 
 145 // SD calculates the current standard-deviation. The only parameter is the
 146 // bias-correction to use:
 147 //
 148 //  0 means calculate the current population standard-deviation
 149 //  1 means calculate the current sample standard-deviation
 150 //  1.5 is very-rarely used and only in special circumstances
 151 func (ns NumberSummary) SD(bias float64) float64 {
 152     return ns.welford.SD(bias, ns.Valid())
 153 }
 154 
 155 // CommonQuantiles groups all the most commonly-used quantiles in practice.
 156 //
 157 // Funcs AppendAllDeciles and AppendAllPercentiles give you other sets of
 158 // commonly-used ranking stats.
 159 type CommonQuantiles struct {
 160     Min float64
 161     P01 float64 // 1st percentile
 162     P05 float64 // 5th percentile
 163     P10 float64 // 10th percentile, also the 1st decile
 164     P25 float64 // 1st quartile, also the 25th percentile
 165     P50 float64 // 2nd quartile, also the 50th percentile
 166     P75 float64 // 3rd quartile, also the 75th percentile
 167     P90 float64
 168     P95 float64
 169     P99 float64
 170     Max float64
 171 }
 172 
 173 // NewCommonQuantiles is a convenience constructor for struct CommonQuantiles.
 174 func NewCommonQuantiles(x []float64) CommonQuantiles {
 175     return CommonQuantiles{
 176         Min: Quantile(x, 0.00),
 177         P01: Quantile(x, 0.01),
 178         P05: Quantile(x, 0.05),
 179         P10: Quantile(x, 0.10),
 180         P25: Quantile(x, 0.25),
 181         P50: Quantile(x, 0.50),
 182         P75: Quantile(x, 0.75),
 183         P90: Quantile(x, 0.90),
 184         P95: Quantile(x, 0.95),
 185         P99: Quantile(x, 0.99),
 186         Max: Quantile(x, 1.00),
 187     }
 188 }
 189 
 190 // AppendAllDeciles appends 11 items to the slice given, which ultimately
 191 // lets you reuse slices for multiple such calculations, thus avoiding extra
 192 // allocations.
 193 //
 194 // The 1st item returned is the minimum, and the last one is the maximum.
 195 func AppendAllDeciles(dest []float64, x []float64) []float64 {
 196     for i := 0; i <= 10; i++ {
 197         dest = append(dest, Quantile(x, float64(i)/10))
 198     }
 199     return dest
 200 }
 201 
 202 // AppendAllPercentiles appends 101 items to the slice given, which ultimately
 203 // lets you reuse slices for multiple such calculations, thus avoiding extra
 204 // allocations.
 205 //
 206 // The 1st item returned is the minimum, while the last one is the maximum.
 207 func AppendAllPercentiles(dest []float64, x []float64) []float64 {
 208     for i := 0; i <= 100; i++ {
 209         dest = append(dest, Quantile(x, float64(i)/100))
 210     }
 211     return dest
 212 }

     File: ./statistics_test.go
   1 package mathplus
   2 
   3 import (
   4     "math"
   5     "testing"
   6 )
   7 
   8 type NumericTestResult struct {
   9     Count     int
  10     Valid     int
  11     Integers  int
  12     Negatives int
  13     Zeros     int
  14     Positives int
  15 
  16     Min     float64
  17     Max     float64
  18     Sum     float64
  19     Mean    float64
  20     Geomean float64
  21     SD      float64
  22     RMS     float64
  23 }
  24 
  25 func (tr NumericTestResult) Match(ns NumberSummary) bool {
  26     return true &&
  27         tr.Count == ns.Count &&
  28         tr.Valid == ns.Valid() &&
  29         tr.Integers == ns.Integers &&
  30         tr.Negatives == ns.Negatives &&
  31         tr.Zeros == ns.Zeros &&
  32         tr.Positives == ns.Positives &&
  33         equals(tr.Min, ns.Min) &&
  34         equals(tr.Max, ns.Max) &&
  35         equals(tr.Sum, ns.Sum) &&
  36         equals(tr.Mean, ns.Mean) &&
  37         equals(tr.Geomean, ns.Geomean()) &&
  38         equals(tr.SD, ns.SD(0)) &&
  39         equals(tr.RMS, ns.RMS()) &&
  40         true
  41 }
  42 
  43 func (tr NumericTestResult) Test(t *testing.T, ns NumberSummary) {
  44     var checks = []struct {
  45         X       float64
  46         Y       float64
  47         Message string
  48     }{
  49         {float64(tr.Count), float64(ns.Count), `count`},
  50         {float64(tr.Valid), float64(ns.Valid()), `valid`},
  51         {float64(tr.Integers), float64(ns.Integers), `integers`},
  52         {float64(tr.Negatives), float64(ns.Negatives), `negatives`},
  53         {float64(tr.Zeros), float64(ns.Zeros), `zeros`},
  54         {float64(tr.Positives), float64(ns.Positives), `positives`},
  55         {float64(tr.Min), float64(ns.Min), `min`},
  56         {float64(tr.Max), float64(ns.Max), `max`},
  57         {float64(tr.Sum), float64(ns.Sum), `sum`},
  58         {float64(tr.Mean), float64(ns.Mean), `mean`},
  59         {float64(tr.Geomean), float64(ns.Geomean()), `geomean`},
  60         {float64(tr.SD), float64(ns.SD(0)), `sd`},
  61         {float64(tr.RMS), float64(ns.RMS()), `rms`},
  62     }
  63 
  64     const fs = "field %q failed: %f and %f differ\nexpected %#v\ngot      %#v"
  65     for _, tc := range checks {
  66         if !equals(tc.X, tc.Y) {
  67             t.Fatalf(fs, tc.Message, tc.X, tc.Y, tr, ns)
  68             return
  69         }
  70     }
  71 }
  72 
  73 var numericTests = []struct {
  74     Description string
  75     Input       []float64
  76     Expected    NumericTestResult
  77 }{
  78     {
  79         Description: `no values`,
  80         Input:       []float64{},
  81         Expected: NumericTestResult{
  82             Count:    0,
  83             Valid:    0,
  84             Integers: 0,
  85             // Min:      math.Inf(+1),
  86             // Max:      math.Inf(-1),
  87             Min:     0.0,
  88             Max:     0.0,
  89             Sum:     0.0,
  90             Mean:    0.0,
  91             Geomean: math.NaN(),
  92             SD:      math.NaN(),
  93             RMS:     0.0,
  94         },
  95     },
  96     {
  97         Description: `just a value`,
  98         Input:       []float64{-3.5},
  99         Expected: NumericTestResult{
 100             Count:     1,
 101             Valid:     1,
 102             Integers:  0,
 103             Negatives: 1,
 104             Zeros:     0,
 105             Positives: 0,
 106             Min:       -3.5,
 107             Max:       -3.5,
 108             Sum:       -3.5,
 109             Mean:      -3.5,
 110             Geomean:   math.NaN(),
 111             SD:        0.0,
 112             RMS:       0.0,
 113         },
 114     },
 115     {
 116         Description: `integers 1..10`,
 117         Input:       []float64{math.NaN(), 1, 2, 3, 4, 5, 6, 7, 8, 9, 10},
 118         Expected: NumericTestResult{
 119             Count:     11,
 120             Valid:     10,
 121             Integers:  10,
 122             Negatives: 0,
 123             Zeros:     0,
 124             Positives: 10,
 125             Min:       1,
 126             Max:       10,
 127             Sum:       55,
 128             Mean:      5.5,
 129             Geomean:   4.528728688116766,
 130             RMS:       9.082951062292475,
 131             SD:        2.8722813232690143,
 132         },
 133     },
 134     {
 135         Description: `nothing valid`,
 136         Input: []float64{
 137             math.NaN(), math.NaN(), math.NaN(), math.NaN(), math.NaN(),
 138         },
 139         Expected: NumericTestResult{
 140             Count:    5,
 141             Valid:    0,
 142             Integers: 0,
 143             Min:      math.Inf(+1),
 144             Max:      math.Inf(-1),
 145             Sum:      0.0,
 146             Mean:     0.0,
 147             Geomean:  math.NaN(),
 148             SD:       math.NaN(),
 149             RMS:      0.0,
 150         },
 151     },
 152 }
 153 
 154 func TestNumberSummary(t *testing.T) {
 155     for _, tc := range numericTests {
 156         var s NumberSummary
 157         for _, x := range tc.Input {
 158             s.Update(x)
 159         }
 160         tc.Expected.Test(t, s)
 161     }
 162 }
 163 
 164 func equals(x, y float64) bool {
 165     if x == y {
 166         return true
 167     }
 168     return math.IsNaN(x) && math.IsNaN(y)
 169 }