github.com/gopherd/gonum@v0.0.4/mathext/internal/cephes/polevl.go (about)

     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  
    10  package cephes
    11  
    12  import "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.
    18  func 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.
    28  func 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
    37  func 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  }