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 */
80long double MAXNUML = 1.189731495357231765021263853E4932L;
81/* 2^-64 */
82long double MACHEPL = 5.42101086242752217003726400434970855712890625E-20L;
83/* log( MAXNUML ) */
84long double MAXLOGL =  1.1356523406294143949492E4L;
85#ifdef DENORMAL
86/* log(smallest denormal number = 2^-16446) */
87long double MINLOGL = -1.13994985314888605586758E4L;
88#else
89/* log( underflow threshold = 2^(-16382) ) */
90long double MINLOGL = -1.1355137111933024058873E4L;
91#endif
92long double LOGE2L  = 6.9314718055994530941723E-1L;
93long double LOG2EL  = 1.4426950408889634073599E0L;
94long double PIL     = 3.1415926535897932384626L;
95long double PIO2L   = 1.5707963267948966192313L;
96long double PIO4L   = 7.8539816339744830961566E-1L;
97#ifdef INFINITIES
98long double NANL = 0.0L / 0.0L;
99long double INFINITYL = 1.0L / 0.0L;
100#else
101long double INFINITYL = 1.189731495357231765021263853E4932L;
102long double NANL = 0.0L;
103#endif
104#endif
105#if IBMPC
106short MAXNUML[] = {0xffff,0xffff,0xffff,0xffff,0x7ffe, XPD};
107short MAXLOGL[] = {0x79ab,0xd1cf,0x17f7,0xb172,0x400c, XPD};
108#ifdef INFINITIES
109short INFINITYL[] = {0,0,0,0x8000,0x7fff, XPD};
110short NANL[] = {0,0,0,0xc000,0x7fff, XPD};
111#else
112short INFINITYL[] = {0xffff,0xffff,0xffff,0xffff,0x7ffe, XPD};
113long double NANL = 0.0L;
114#endif
115#ifdef DENORMAL
116short MINLOGL[] = {0xbaaa,0x09e2,0xfe7f,0xb21d,0xc00c, XPD};
117#else
118short MINLOGL[] = {0xeb2f,0x1210,0x8c67,0xb16c,0xc00c, XPD};
119#endif
120short MACHEPL[] = {0x0000,0x0000,0x0000,0x8000,0x3fbf, XPD};
121short LOGE2L[]  = {0x79ac,0xd1cf,0x17f7,0xb172,0x3ffe, XPD};
122short LOG2EL[]  = {0xf0bc,0x5c17,0x3b29,0xb8aa,0x3fff, XPD};
123short PIL[]     = {0xc235,0x2168,0xdaa2,0xc90f,0x4000, XPD};
124short PIO2L[]   = {0xc235,0x2168,0xdaa2,0xc90f,0x3fff, XPD};
125short PIO4L[]   = {0xc235,0x2168,0xdaa2,0xc90f,0x3ffe, XPD};
126#endif
127#if MIEEE
128long MAXNUML[] = {0x7ffe0000,0xffffffff,0xffffffff};
129long MAXLOGL[] = {0x400c0000,0xb17217f7,0xd1cf79ab};
130#ifdef INFINITIES
131long INFINITY[] = {0x7fff0000,0x80000000,0x00000000};
132long NANL[] = {0x7fff0000,0xffffffff,0xffffffff};
133#else
134long INFINITYL[] = {0x7ffe0000,0xffffffff,0xffffffff};
135long double NANL = 0.0L;
136#endif
137#ifdef DENORMAL
138long MINLOGL[] = {0xc00c0000,0xb21dfe7f,0x09e2baaa};
139#else
140long MINLOGL[] = {0xc00c0000,0xb16c8c67,0x1210eb2f};
141#endif
142long MACHEPL[] = {0x3fbf0000,0x80000000,0x00000000};
143long LOGE2L[]  = {0x3ffe0000,0xb17217f7,0xd1cf79ac};
144long LOG2EL[]  = {0x3fff0000,0xb8aa3b29,0x5c17f0bc};
145long PIL[]     = {0x40000000,0xc90fdaa2,0x2168c235};
146long PIO2L[]   = {0x3fff0000,0xc90fdaa2,0x2168c235};
147long PIO4L[]   = {0x3ffe0000,0xc90fdaa2,0x2168c235};
148#endif
149
150#ifdef MINUSZERO
151long double NEGZEROL = -0.0L;
152#else
153long double NEGZEROL = 0.0L;
154#endif
155
156/* Polynomial evaluator:
157 *  P[0] x^n  +  P[1] x^(n-1)  +  ...  +  P[n]
158 */
159long double polevll( x, p, n )
160long double x;
161void *p;
162int n;
163{
164register long double y;
165register long double *P = (long double *)p;
166
167y = *P++;
168do
169	{
170	y = y * x + *P++;
171	}
172while( --n );
173return(y);
174}
175
176
177
178/* Polynomial evaluator:
179 *  x^n  +  P[0] x^(n-1)  +  P[1] x^(n-2)  +  ...  +  P[n]
180 */
181long double p1evll( x, p, n )
182long double x;
183void *p;
184int n;
185{
186register long double y;
187register long double *P = (long double *)p;
188
189n -= 1;
190y = x + *P++;
191do
192	{
193	y = y * x + *P++;
194	}
195while( --n );
196return( y );
197}
198#endif	/* LDBL_MANT_DIG == 64 */
199