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 */
65long 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 */
81long 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