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