/fn/fn.go

  1// Copyright 2012 - 2013 The Fn Authors. All rights reserved. See the LICENSE file.
  2
  3package fn
  4
  5// Special math functions.
  6
  7import (
  8	"math"
  9)
 10
 11const π = float64(math.Pi)
 12const ln2 = math.Ln2
 13const lnSqrt2π = 0.918938533204672741780329736406 // log(sqrt(2*pi))
 14const min64 = math.SmallestNonzeroFloat64         //   DBL_MIN
 15const eps64 = 1.1102230246251565e-16              // DBL_EPSILON   
 16const maxExp = 1024.0                             // DBL_MAX_EXP
 17const sqrt2 = math.Sqrt2
 18const lnSqrtπd2 = 0.225791352644727432363097614947 // log(sqrt(pi/2))  M_LN_SQRT_PId2
 19
 20var nan = math.NaN()
 21
 22var fZero float64 = float64(0.0)
 23var fOne float64 = float64(1.0)
 24var iZero int64 = int64(0)
 25var iOne int64 = int64(1)
 26
 27var negInf float64 = math.Inf(-1)
 28var posInf float64 = math.Inf(+1)
 29
 30// Functions imported from "math"
 31var abs func(float64) float64 = math.Abs
 32var floor func(float64) float64 = math.Floor
 33var ceil func(float64) float64 = math.Ceil
 34var log func(float64) float64 = math.Log
 35var log1p func(float64) float64 = math.Log1p
 36var exp func(float64) float64 = math.Exp
 37var sqrt func(float64) float64 = math.Sqrt
 38var pow func(float64, float64) float64 = math.Pow
 39var atan func(float64) float64 = math.Atan
 40var tan func(float64) float64 = math.Tan
 41var sin func(float64) float64 = math.Sin
 42var trunc func(float64) float64 = math.Trunc
 43var erf func(float64) float64 = math.Erf
 44var erfc func(float64) float64 = math.Erfc
 45var isNaN func(float64) bool = math.IsNaN
 46var isInf func(float64, int) bool = math.IsInf
 47var fmod func(float64, float64) float64 = math.Mod
 48
 49var Γ = math.Gamma
 50var GammaF = math.Gamma
 51var sqrt2pi = math.Sqrt(2 * math.Pi)
 52var logsqrt2pi = math.Log(math.Sqrt(2 * math.Pi))
 53var LnΓp = LnGammaP
 54var LnΓpRatio = LnGammaPRatio
 55
 56var lanczos_coef []float64 = []float64{
 57	0.99999999999980993,
 58	676.5203681218851,
 59	-1259.1392167224028,
 60	771.32342877765313,
 61	-176.61502916214059,
 62	12.507343278686905,
 63	-0.13857109526572012,
 64	9.9843695780195716e-6,
 65	1.5056327351493116e-7}
 66
 67func isOdd(k float64) bool {
 68	if k != 2*floor(k/2.0) {
 69		return true
 70	}
 71	return false
 72}
 73
 74func isInt(x float64) bool {
 75	if abs((x)-floor((x)+0.5)) <= 1e-7 {
 76		return true
 77	}
 78	return false
 79}
 80
 81// Round to nearest integer
 82func Round(x float64) float64 {
 83	var i float64
 84	f := math.Floor(x)
 85	c := math.Ceil(x)
 86	if x-f < c-x {
 87		i = f
 88	} else {
 89		i = c
 90	}
 91	return i
 92}
 93
 94// Arithmetic mean
 95func ArithMean(data *Vector) float64 {
 96	n := data.L
 97	sum := 0.0
 98	for i := 0; i < n; i++ {
 99		sum += data.Get(i)
100	}
101	return sum / float64(n)
102}
103
104// Geometric mean
105func GeomMean(data *Vector) float64 {
106	n := data.L
107	sum := 0.0
108	for i := 0; i < n; i++ {
109		sum += math.Log(data.Get(i))
110	}
111	return math.Exp(sum / float64(n))
112}
113
114// Harmonic mean
115func HarmonicMean(data *Vector) float64 {
116	n := data.L
117	sum := 0.0
118	for i := 0; i < n; i++ {
119		sum += 1.0 / data.Get(i)
120	}
121	return float64(n) / sum
122}
123
124// Generalized mean
125func GenMean(data *Vector, p float64) float64 {
126	n := data.L
127	sum := 0.0
128	for i := 0; i < n; i++ {
129		sum += math.Pow(data.Get(i), p)
130	}
131	return math.Pow(sum/float64(n), 1/p)
132}
133
134// Bernoulli number
135// Akiyama–Tanigawa algorithm for Bn
136func Bn(n int64) float64 {
137	var m int64
138	a := make([]float64, n)
139	for m = 0; m <= n; m++ {
140		a[m] = 1 / float64(m+1)
141		for j := m; j >= 1; j-- {
142			a[j-1] = float64(j) * (a[j-1] - a[j])
143		}
144	}
145	return a[0] // (which is Bn)
146}
147
148// H returns the generalized harmonic number of order n of m. 
149func H(n int64, m float64) float64 {
150	var i int64
151	h := 0.0
152	for i = 1; i <= n; i++ {
153		h += math.Pow(float64(i), m)
154	}
155	return h
156}
157
158// Generalized harmonic number
159func H2(n int64, q, s float64) float64 {
160	var i int64
161	h := 0.0
162	for i = 1; i <= n; i++ {
163		h += math.Pow((float64(i) + q), -s)
164	}
165	return h
166}