1 /* $OpenBSD: polevll.c,v 1.2 2011/06/02 21:47:40 martynas Exp $ */ 2 3 /* 4 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 5 * 6 * Permission to use, copy, modify, and distribute this software for any 7 * purpose with or without fee is hereby granted, provided that the above 8 * copyright notice and this permission notice appear in all copies. 9 * 10 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 11 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 12 * MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR 13 * ANY SPECIAL, DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES 14 * WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN 15 * ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF 16 * OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. 17 */ 18 19 /* polevll.c 20 * p1evll.c 21 * 22 * Evaluate polynomial 23 * 24 * 25 * 26 * SYNOPSIS: 27 * 28 * int N; 29 * long double x, y, coef[N+1], polevl[]; 30 * 31 * y = polevll( x, coef, N ); 32 * 33 * 34 * 35 * DESCRIPTION: 36 * 37 * Evaluates polynomial of degree N: 38 * 39 * 2 N 40 * y = C + C x + C x +...+ C x 41 * 0 1 2 N 42 * 43 * Coefficients are stored in reverse order: 44 * 45 * coef[0] = C , ..., coef[N] = C . 46 * N 0 47 * 48 * The function p1evll() assumes that coef[N] = 1.0 and is 49 * omitted from the array. Its calling arguments are 50 * otherwise the same as polevll(). 51 * 52 * This module also contains the following globally declared constants: 53 * MAXNUML = 1.189731495357231765021263853E4932L; 54 * MACHEPL = 5.42101086242752217003726400434970855712890625E-20L; 55 * MAXLOGL = 1.1356523406294143949492E4L; 56 * MINLOGL = -1.1355137111933024058873E4L; 57 * LOGE2L = 6.9314718055994530941723E-1L; 58 * LOG2EL = 1.4426950408889634073599E0L; 59 * PIL = 3.1415926535897932384626L; 60 * PIO2L = 1.5707963267948966192313L; 61 * PIO4L = 7.8539816339744830961566E-1L; 62 * 63 * SPEED: 64 * 65 * In the interest of speed, there are no checks for out 66 * of bounds arithmetic. This routine is used by most of 67 * the functions in the library. Depending on available 68 * equipment features, the user may wish to rewrite the 69 * program in microcode or assembly language. 70 * 71 */ 72 73 #include <float.h> 74 75 #if LDBL_MANT_DIG == 64 76 #include "mconf.h" 77 78 #if UNK 79 /* almost 2^16384 */ 80 long double MAXNUML = 1.189731495357231765021263853E4932L; 81 /* 2^-64 */ 82 long double MACHEPL = 5.42101086242752217003726400434970855712890625E-20L; 83 /* log( MAXNUML ) */ 84 long double MAXLOGL = 1.1356523406294143949492E4L; 85 #ifdef DENORMAL 86 /* log(smallest denormal number = 2^-16446) */ 87 long double MINLOGL = -1.13994985314888605586758E4L; 88 #else 89 /* log( underflow threshold = 2^(-16382) ) */ 90 long double MINLOGL = -1.1355137111933024058873E4L; 91 #endif 92 long double LOGE2L = 6.9314718055994530941723E-1L; 93 long double LOG2EL = 1.4426950408889634073599E0L; 94 long double PIL = 3.1415926535897932384626L; 95 long double PIO2L = 1.5707963267948966192313L; 96 long double PIO4L = 7.8539816339744830961566E-1L; 97 #ifdef INFINITIES 98 long double NANL = 0.0L / 0.0L; 99 long double INFINITYL = 1.0L / 0.0L; 100 #else 101 long double INFINITYL = 1.189731495357231765021263853E4932L; 102 long double NANL = 0.0L; 103 #endif 104 #endif 105 #if IBMPC 106 short MAXNUML[] = {0xffff,0xffff,0xffff,0xffff,0x7ffe, XPD}; 107 short MAXLOGL[] = {0x79ab,0xd1cf,0x17f7,0xb172,0x400c, XPD}; 108 #ifdef INFINITIES 109 short INFINITYL[] = {0,0,0,0x8000,0x7fff, XPD}; 110 short NANL[] = {0,0,0,0xc000,0x7fff, XPD}; 111 #else 112 short INFINITYL[] = {0xffff,0xffff,0xffff,0xffff,0x7ffe, XPD}; 113 long double NANL = 0.0L; 114 #endif 115 #ifdef DENORMAL 116 short MINLOGL[] = {0xbaaa,0x09e2,0xfe7f,0xb21d,0xc00c, XPD}; 117 #else 118 short MINLOGL[] = {0xeb2f,0x1210,0x8c67,0xb16c,0xc00c, XPD}; 119 #endif 120 short MACHEPL[] = {0x0000,0x0000,0x0000,0x8000,0x3fbf, XPD}; 121 short LOGE2L[] = {0x79ac,0xd1cf,0x17f7,0xb172,0x3ffe, XPD}; 122 short LOG2EL[] = {0xf0bc,0x5c17,0x3b29,0xb8aa,0x3fff, XPD}; 123 short PIL[] = {0xc235,0x2168,0xdaa2,0xc90f,0x4000, XPD}; 124 short PIO2L[] = {0xc235,0x2168,0xdaa2,0xc90f,0x3fff, XPD}; 125 short PIO4L[] = {0xc235,0x2168,0xdaa2,0xc90f,0x3ffe, XPD}; 126 #endif 127 #if MIEEE 128 long MAXNUML[] = {0x7ffe0000,0xffffffff,0xffffffff}; 129 long MAXLOGL[] = {0x400c0000,0xb17217f7,0xd1cf79ab}; 130 #ifdef INFINITIES 131 long INFINITY[] = {0x7fff0000,0x80000000,0x00000000}; 132 long NANL[] = {0x7fff0000,0xffffffff,0xffffffff}; 133 #else 134 long INFINITYL[] = {0x7ffe0000,0xffffffff,0xffffffff}; 135 long double NANL = 0.0L; 136 #endif 137 #ifdef DENORMAL 138 long MINLOGL[] = {0xc00c0000,0xb21dfe7f,0x09e2baaa}; 139 #else 140 long MINLOGL[] = {0xc00c0000,0xb16c8c67,0x1210eb2f}; 141 #endif 142 long MACHEPL[] = {0x3fbf0000,0x80000000,0x00000000}; 143 long LOGE2L[] = {0x3ffe0000,0xb17217f7,0xd1cf79ac}; 144 long LOG2EL[] = {0x3fff0000,0xb8aa3b29,0x5c17f0bc}; 145 long PIL[] = {0x40000000,0xc90fdaa2,0x2168c235}; 146 long PIO2L[] = {0x3fff0000,0xc90fdaa2,0x2168c235}; 147 long PIO4L[] = {0x3ffe0000,0xc90fdaa2,0x2168c235}; 148 #endif 149 150 #ifdef MINUSZERO 151 long double NEGZEROL = -0.0L; 152 #else 153 long double NEGZEROL = 0.0L; 154 #endif 155 156 /* Polynomial evaluator: 157 * P[0] x^n + P[1] x^(n-1) + ... + P[n] 158 */ 159 long double polevll( x, p, n ) 160 long double x; 161 void *p; 162 int n; 163 { 164 register long double y; 165 register long double *P = (long double *)p; 166 167 y = *P++; 168 do 169 { 170 y = y * x + *P++; 171 } 172 while( --n ); 173 return(y); 174 } 175 176 177 178 /* Polynomial evaluator: 179 * x^n + P[0] x^(n-1) + P[1] x^(n-2) + ... + P[n] 180 */ 181 long double p1evll( x, p, n ) 182 long double x; 183 void *p; 184 int n; 185 { 186 register long double y; 187 register long double *P = (long double *)p; 188 189 n -= 1; 190 y = x + *P++; 191 do 192 { 193 y = y * x + *P++; 194 } 195 while( --n ); 196 return( y ); 197 } 198 #endif /* LDBL_MANT_DIG == 64 */ 199