polevl.h
Go to the documentation of this file.
1 /* polevl.c
2  * p1evl.c
3  *
4  * Evaluate polynomial
5  *
6  *
7  *
8  * SYNOPSIS:
9  *
10  * int N;
11  * double x, y, coef[N+1], polevl[];
12  *
13  * y = polevl( x, coef, N );
14  *
15  *
16  *
17  * DESCRIPTION:
18  *
19  * Evaluates polynomial of degree N:
20  *
21  * 2 N
22  * y = C + C x + C x +...+ C x
23  * 0 1 2 N
24  *
25  * Coefficients are stored in reverse order:
26  *
27  * coef[0] = C , ..., coef[N] = C .
28  * N 0
29  *
30  * The function p1evl() assumes that c_N = 1.0 so that coefficent
31  * is omitted from the array. Its calling arguments are
32  * otherwise the same as polevl().
33  *
34  *
35  * SPEED:
36  *
37  * In the interest of speed, there are no checks for out
38  * of bounds arithmetic. This routine is used by most of
39  * the functions in the library. Depending on available
40  * equipment features, the user may wish to rewrite the
41  * program in microcode or assembly language.
42  *
43  */
44 
45 /*
46  * Cephes Math Library Release 2.1: December, 1988
47  * Copyright 1984, 1987, 1988 by Stephen L. Moshier
48  * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
49  */
50 
51 /* Sources:
52  * [1] Holin et. al., "Polynomial and Rational Function Evaluation",
53  * https://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/roots/rational.html
54  */
55 
56 /* Scipy changes:
57  * - 06-23-2016: add code for evaluating rational functions
58  */
59 
60 #ifndef CEPHES_POLEV
61 #define CEPHES_POLEV
62 
63 #include <math.h>
64 
65 static inline double polevl(double x, const double coef[], int N)
66 {
67  double ans;
68  int i;
69  const double *p;
70 
71  p = coef;
72  ans = *p++;
73  i = N;
74 
75  do
76  ans = ans * x + *p++;
77  while (--i);
78 
79  return (ans);
80 }
81 
82 /* p1evl() */
83 /* N
84  * Evaluate polynomial when coefficient of x is 1.0.
85  * That is, C_{N} is assumed to be 1, and that coefficient
86  * is not included in the input array coef.
87  * coef must have length N and contain the polynomial coefficients
88  * stored as
89  * coef[0] = C_{N-1}
90  * coef[1] = C_{N-2}
91  * ...
92  * coef[N-2] = C_1
93  * coef[N-1] = C_0
94  * Otherwise same as polevl.
95  */
96 
97 static inline double p1evl(double x, const double coef[], int N)
98 {
99  double ans;
100  const double *p;
101  int i;
102 
103  p = coef;
104  ans = x + *p++;
105  i = N - 1;
106 
107  do
108  ans = ans * x + *p++;
109  while (--i);
110 
111  return (ans);
112 }
113 
114 /* Evaluate a rational function. See [1]. */
115 
116 static inline double ratevl(double x, const double num[], int M,
117  const double denom[], int N)
118 {
119  int i, dir;
120  double y, num_ans, denom_ans;
121  double absx = fabs(x);
122  const double *p;
123 
124  if (absx > 1) {
125  /* Evaluate as a polynomial in 1/x. */
126  dir = -1;
127  p = num + M;
128  y = 1 / x;
129  } else {
130  dir = 1;
131  p = num;
132  y = x;
133  }
134 
135  /* Evaluate the numerator */
136  num_ans = *p;
137  p += dir;
138  for (i = 1; i <= M; i++) {
139  num_ans = num_ans * y + *p;
140  p += dir;
141  }
142 
143  /* Evaluate the denominator */
144  if (absx > 1) {
145  p = denom + N;
146  } else {
147  p = denom;
148  }
149 
150  denom_ans = *p;
151  p += dir;
152  for (i = 1; i <= N; i++) {
153  denom_ans = denom_ans * y + *p;
154  p += dir;
155  }
156 
157  if (absx > 1) {
158  i = N - M;
159  return pow(x, i) * num_ans / denom_ans;
160  } else {
161  return num_ans / denom_ans;
162  }
163 }
164 
165 #endif
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
p1evl
static double p1evl(double x, const double coef[], int N)
Definition: polevl.h:97
ratevl
static double ratevl(double x, const double num[], int M, const double denom[], int N)
Definition: polevl.h:116
polevl
static double polevl(double x, const double coef[], int N)
Definition: polevl.h:65
y
Scalar * y
Definition: level1_cplx_impl.h:124
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
p
float * p
Definition: Tutorial_Map_using.cpp:9
N
#define N
Definition: igam.h:9
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
M
Matrix< RealScalar, Dynamic, Dynamic > M
Definition: bench_gemm.cpp:51


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:34:23