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 */
polevll(x,p,n)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 */
p1evll(x,p,n)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