1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
// Derived from SciPy's special/cephes/polevl.h
// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/polevl.h
// Made freely available by Stephen L. Moshier without support or guarantee.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
// Copyright ©1984, ©1987, ©1988 by Stephen L. Moshier
// Portions Copyright ©2016 The Gonum Authors. All rights reserved.
package cephes
import "math"
// polevl evaluates a polynomial of degree N
//
// y = c_0 + c_1 x_1 + c_2 x_2^2 ...
//
// where the coefficients are stored in reverse order, i.e. coef[0] = c_n and
// coef[n] = c_0.
func polevl(x float64, coef []float64, n int) float64 {
ans := coef[0]
for i := 1; i <= n; i++ {
ans = ans*x + coef[i]
}
return ans
}
// p1evl is the same as polevl, except c_n is assumed to be 1 and is not included
// in the slice.
func p1evl(x float64, coef []float64, n int) float64 {
ans := x + coef[0]
for i := 1; i <= n-1; i++ {
ans = ans*x + coef[i]
}
return ans
}
// ratevl evaluates a rational function
func ratevl(x float64, num []float64, m int, denom []float64, n int) float64 {
// Source: Holin et. al., "Polynomial and Rational Function Evaluation",
// http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/roots/rational.html
absx := math.Abs(x)
var dir, idx int
var y float64
if absx > 1 {
// Evaluate as a polynomial in 1/x
dir = -1
idx = m
y = 1 / x
} else {
dir = 1
idx = 0
y = x
}
// Evaluate the numerator
numAns := num[idx]
idx += dir
for i := 0; i < m; i++ {
numAns = numAns*y + num[idx]
idx += dir
}
// Evaluate the denominator
if absx > 1 {
idx = n
} else {
idx = 0
}
denomAns := denom[idx]
idx += dir
for i := 0; i < n; i++ {
denomAns = denomAns*y + denom[idx]
idx += dir
}
if absx > 1 {
pow := float64(n - m)
return math.Pow(x, pow) * numAns / denomAns
}
return numAns / denomAns
}
|
