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 94