[398] | 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
|
---|