/* origin: OpenBSD /usr/src/lib/libm/src/polevll.c */ /* * Copyright (c) 2008 Stephen L. Moshier * * Permission to use, copy, modify, and distribute this software for any * purpose with or without fee is hereby granted, provided that the above * copyright notice and this permission notice appear in all copies. * * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. */ /* * Evaluate polynomial * * * SYNOPSIS: * * int N; * long double x, y, coef[N+1], polevl[]; * * y = polevll( x, coef, N ); * * * DESCRIPTION: * * Evaluates polynomial of degree N: * * 2 N * y = C + C x + C x +...+ C x * 0 1 2 N * * Coefficients are stored in reverse order: * * coef[0] = C , ..., coef[N] = C . * N 0 * * The function p1evll() assumes that coef[N] = 1.0 and is * omitted from the array. Its calling arguments are * otherwise the same as polevll(). * * * SPEED: * * In the interest of speed, there are no checks for out * of bounds arithmetic. This routine is used by most of * the functions in the library. Depending on available * equipment features, the user may wish to rewrite the * program in microcode or assembly language. * */ #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 #else /* * Polynomial evaluator: * P[0] x^n + P[1] x^(n-1) + ... + P[n] */ long double __polevll(long double x, const long double *P, int n) { long double y; y = *P++; do { y = y * x + *P++; } while (--n); return y; } /* * Polynomial evaluator: * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n] */ long double __p1evll(long double x, const long double *P, int n) { long double y; n -= 1; y = x + *P++; do { y = y * x + *P++; } while (--n); return y; } #endif