gonum.org/v1/gonum@v0.14.0/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  //
    16  //	y = c_0 + c_1 x_1 + c_2 x_2^2 ...
    17  //
    18  // where the coefficients are stored in reverse order, i.e. coef[0] = c_n and
    19  // coef[n] = c_0.
    20  func polevl(x float64, coef []float64, n int) float64 {
    21  	ans := coef[0]
    22  	for i := 1; i <= n; i++ {
    23  		ans = ans*x + coef[i]
    24  	}
    25  	return ans
    26  }
    27  
    28  // p1evl is the same as polevl, except c_n is assumed to be 1 and is not included
    29  // in the slice.
    30  func p1evl(x float64, coef []float64, n int) float64 {
    31  	ans := x + coef[0]
    32  	for i := 1; i <= n-1; i++ {
    33  		ans = ans*x + coef[i]
    34  	}
    35  	return ans
    36  }
    37  
    38  // ratevl evaluates a rational function
    39  func ratevl(x float64, num []float64, m int, denom []float64, n int) float64 {
    40  	// Source: Holin et. al., "Polynomial and Rational Function Evaluation",
    41  	// http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/roots/rational.html
    42  	absx := math.Abs(x)
    43  
    44  	var dir, idx int
    45  	var y float64
    46  	if absx > 1 {
    47  		// Evaluate as a polynomial in 1/x
    48  		dir = -1
    49  		idx = m
    50  		y = 1 / x
    51  	} else {
    52  		dir = 1
    53  		idx = 0
    54  		y = x
    55  	}
    56  
    57  	// Evaluate the numerator
    58  	numAns := num[idx]
    59  	idx += dir
    60  	for i := 0; i < m; i++ {
    61  		numAns = numAns*y + num[idx]
    62  		idx += dir
    63  	}
    64  
    65  	// Evaluate the denominator
    66  	if absx > 1 {
    67  		idx = n
    68  	} else {
    69  		idx = 0
    70  	}
    71  
    72  	denomAns := denom[idx]
    73  	idx += dir
    74  	for i := 0; i < n; i++ {
    75  		denomAns = denomAns*y + denom[idx]
    76  		idx += dir
    77  	}
    78  
    79  	if absx > 1 {
    80  		pow := float64(n - m)
    81  		return math.Pow(x, pow) * numAns / denomAns
    82  	}
    83  	return numAns / denomAns
    84  }