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