1*05a0b428SJohn Marino /* $OpenBSD: polevll.c,v 1.2 2013/11/12 20:35:09 martynas Exp $ */
2*05a0b428SJohn Marino
3*05a0b428SJohn Marino /*
4*05a0b428SJohn Marino * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
5*05a0b428SJohn Marino *
6*05a0b428SJohn Marino * Permission to use, copy, modify, and distribute this software for any
7*05a0b428SJohn Marino * purpose with or without fee is hereby granted, provided that the above
8*05a0b428SJohn Marino * copyright notice and this permission notice appear in all copies.
9*05a0b428SJohn Marino *
10*05a0b428SJohn Marino * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
11*05a0b428SJohn Marino * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
12*05a0b428SJohn Marino * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR
13*05a0b428SJohn Marino * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
14*05a0b428SJohn Marino * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
15*05a0b428SJohn Marino * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
16*05a0b428SJohn Marino * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
17*05a0b428SJohn Marino */
18*05a0b428SJohn Marino
19*05a0b428SJohn Marino /* polevll.c
20*05a0b428SJohn Marino * p1evll.c
21*05a0b428SJohn Marino *
22*05a0b428SJohn Marino * Evaluate polynomial
23*05a0b428SJohn Marino *
24*05a0b428SJohn Marino *
25*05a0b428SJohn Marino *
26*05a0b428SJohn Marino * SYNOPSIS:
27*05a0b428SJohn Marino *
28*05a0b428SJohn Marino * int N;
29*05a0b428SJohn Marino * long double x, y, coef[N+1], polevl[];
30*05a0b428SJohn Marino *
31*05a0b428SJohn Marino * y = polevll( x, coef, N );
32*05a0b428SJohn Marino *
33*05a0b428SJohn Marino *
34*05a0b428SJohn Marino *
35*05a0b428SJohn Marino * DESCRIPTION:
36*05a0b428SJohn Marino *
37*05a0b428SJohn Marino * Evaluates polynomial of degree N:
38*05a0b428SJohn Marino *
39*05a0b428SJohn Marino * 2 N
40*05a0b428SJohn Marino * y = C + C x + C x +...+ C x
41*05a0b428SJohn Marino * 0 1 2 N
42*05a0b428SJohn Marino *
43*05a0b428SJohn Marino * Coefficients are stored in reverse order:
44*05a0b428SJohn Marino *
45*05a0b428SJohn Marino * coef[0] = C , ..., coef[N] = C .
46*05a0b428SJohn Marino * N 0
47*05a0b428SJohn Marino *
48*05a0b428SJohn Marino * The function p1evll() assumes that coef[N] = 1.0 and is
49*05a0b428SJohn Marino * omitted from the array. Its calling arguments are
50*05a0b428SJohn Marino * otherwise the same as polevll().
51*05a0b428SJohn Marino *
52*05a0b428SJohn Marino *
53*05a0b428SJohn Marino * SPEED:
54*05a0b428SJohn Marino *
55*05a0b428SJohn Marino * In the interest of speed, there are no checks for out
56*05a0b428SJohn Marino * of bounds arithmetic. This routine is used by most of
57*05a0b428SJohn Marino * the functions in the library. Depending on available
58*05a0b428SJohn Marino * equipment features, the user may wish to rewrite the
59*05a0b428SJohn Marino * program in microcode or assembly language.
60*05a0b428SJohn Marino *
61*05a0b428SJohn Marino */
62*05a0b428SJohn Marino
63*05a0b428SJohn Marino #include <math.h>
64*05a0b428SJohn Marino
65*05a0b428SJohn Marino #include "math_private.h"
66*05a0b428SJohn Marino
67*05a0b428SJohn Marino /*
68*05a0b428SJohn Marino * Polynomial evaluator:
69*05a0b428SJohn Marino * P[0] x^n + P[1] x^(n-1) + ... + P[n]
70*05a0b428SJohn Marino */
71*05a0b428SJohn Marino long double
__polevll(long double x,void * PP,int n)72*05a0b428SJohn Marino __polevll(long double x, void *PP, int n)
73*05a0b428SJohn Marino {
74*05a0b428SJohn Marino long double y;
75*05a0b428SJohn Marino long double *P;
76*05a0b428SJohn Marino
77*05a0b428SJohn Marino P = (long double *)PP;
78*05a0b428SJohn Marino y = *P++;
79*05a0b428SJohn Marino do {
80*05a0b428SJohn Marino y = y * x + *P++;
81*05a0b428SJohn Marino } while (--n);
82*05a0b428SJohn Marino
83*05a0b428SJohn Marino return (y);
84*05a0b428SJohn Marino }
85*05a0b428SJohn Marino
86*05a0b428SJohn Marino /*
87*05a0b428SJohn Marino * Polynomial evaluator:
88*05a0b428SJohn Marino * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n]
89*05a0b428SJohn Marino */
90*05a0b428SJohn Marino long double
__p1evll(long double x,void * PP,int n)91*05a0b428SJohn Marino __p1evll(long double x, void *PP, int n)
92*05a0b428SJohn Marino {
93*05a0b428SJohn Marino long double y;
94*05a0b428SJohn Marino long double *P;
95*05a0b428SJohn Marino
96*05a0b428SJohn Marino P = (long double *)PP;
97*05a0b428SJohn Marino n -= 1;
98*05a0b428SJohn Marino y = x + *P++;
99*05a0b428SJohn Marino do {
100*05a0b428SJohn Marino y = y * x + *P++;
101*05a0b428SJohn Marino } while (--n);
102*05a0b428SJohn Marino
103*05a0b428SJohn Marino return (y);
104*05a0b428SJohn Marino }
105