Mercurial > nand-ecoscentric

comparison packages/language/c/libm/current/src/double/internal/k_rem_pio2.c @ 0:3111d98ba7b3 ecos-v1_1-release

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial commit of eCos version 1.1
author jlarmour
date Tue, 11 May 1999 11:16:07 +0000
parents
children 443894e2e912
comparison
equal deleted inserted replaced
-1:000000000000 0:3111d98ba7b3
1 //===========================================================================
2 //
3 // k_rem_pio2.c
4 //
5 // Part of the standard mathematical function library
6 //
7 //===========================================================================
8 //####COPYRIGHTBEGIN####
9 //
10 // -------------------------------------------
11 // The contents of this file are subject to the Cygnus eCos Public License
12 // Version 1.0 (the "License"); you may not use this file except in
13 // compliance with the License. You may obtain a copy of the License at
14 // http://sourceware.cygnus.com/ecos
15 //
16 // Software distributed under the License is distributed on an "AS IS"
17 // basis, WITHOUT WARRANTY OF ANY KIND, either express or implied. See the
18 // License for the specific language governing rights and limitations under
19 // the License.
20 //
21 // The Original Code is eCos - Embedded Cygnus Operating System, released
22 // September 30, 1998.
23 //
24 // The Initial Developer of the Original Code is Cygnus. Portions created
25 // by Cygnus are Copyright (C) 1998 Cygnus Solutions. All Rights Reserved.
26 // -------------------------------------------
27 //
28 //####COPYRIGHTEND####
29 //===========================================================================
30 //#####DESCRIPTIONBEGIN####
31 //
32 // Author(s): jlarmour@cygnus.co.uk
33 // Contributors: jlarmour@cygnus.co.uk
34 // Date: 1998-02-13
35 // Purpose:
36 // Description:
37 // Usage:
38 //
39 //####DESCRIPTIONEND####
40 //
41 //===========================================================================
42
43 // CONFIGURATION
44
45 #include <pkgconf/libm.h> // Configuration header
46
47 // Include the Math library?
48 #ifdef CYGPKG_LIBM
49
50 // Derived from code with the following copyright
51
52
53 /* @(#)k_rem_pio2.c 1.3 95/01/18 */
54 /*
55 * ====================================================
56 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
57 *
58 * Developed at SunSoft, a Sun Microsystems, Inc. business.
59 * Permission to use, copy, modify, and distribute this
60 * software is freely granted, provided that this notice
61 * is preserved.
62 * ====================================================
63 */
64
65 /*
66 * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
67 * double x[],y[]; int e0,nx,prec; int ipio2[];
68 *
69 * __kernel_rem_pio2 return the last three digits of N with
70 * y = x - N*pi/2
71 * so that |y| < pi/2.
72 *
73 * The method is to compute the integer (mod 8) and fraction parts of
74 * (2/pi)*x without doing the full multiplication. In general we
75 * skip the part of the product that are known to be a huge integer (
76 * more accurately, = 0 mod 8 ). Thus the number of operations are
77 * independent of the exponent of the input.
78 *
79 * (2/pi) is represented by an array of 24-bit integers in ipio2[].
80 *
81 * Input parameters:
82 * x[] The input value (must be positive) is broken into nx
83 * pieces of 24-bit integers in double precision format.
84 * x[i] will be the i-th 24 bit of x. The scaled exponent
85 * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
86 * match x's up to 24 bits.
87 *
88 * Example of breaking a double positive z into x[0]+x[1]+x[2]:
89 * e0 = ilogb(z)-23
90 * z = scalbn(z,-e0)
91 * for i = 0,1,2
92 * x[i] = floor(z)
93 * z = (z-x[i])*2**24
94 *
95 *
96 * y[] ouput result in an array of double precision numbers.
97 * The dimension of y[] is:
98 * 24-bit precision 1
99 * 53-bit precision 2
100 * 64-bit precision 2
101 * 113-bit precision 3
102 * The actual value is the sum of them. Thus for 113-bit
103 * precison, one may have to do something like:
104 *
105 * long double t,w,r_head, r_tail;
106 * t = (long double)y[2] + (long double)y[1];
107 * w = (long double)y[0];
108 * r_head = t+w;
109 * r_tail = w - (r_head - t);
110 *
111 * e0 The exponent of x[0]
112 *
113 * nx dimension of x[]
114 *
115 * prec an integer indicating the precision:
116 * 0 24 bits (single)
117 * 1 53 bits (double)
118 * 2 64 bits (extended)
119 * 3 113 bits (quad)
120 *
121 * ipio2[]
122 * integer array, contains the (24*i)-th to (24*i+23)-th
123 * bit of 2/pi after binary point. The corresponding
124 * floating value is
125 *
126 * ipio2[i] * 2^(-24(i+1)).
127 *
128 * External function:
129 * double scalbn(), floor();
130 *
131 *
132 * Here is the description of some local variables:
133 *
134 * jk jk+1 is the initial number of terms of ipio2[] needed
135 * in the computation. The recommended value is 2,3,4,
136 * 6 for single, double, extended,and quad.
137 *
138 * jz local integer variable indicating the number of
139 * terms of ipio2[] used.
140 *
141 * jx nx - 1
142 *
143 * jv index for pointing to the suitable ipio2[] for the
144 * computation. In general, we want
145 * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
146 * is an integer. Thus
147 * e0-3-24*jv >= 0 or (e0-3)/24 >= jv
148 * Hence jv = max(0,(e0-3)/24).
149 *
150 * jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
151 *
152 * q[] double array with integral value, representing the
153 * 24-bits chunk of the product of x and 2/pi.
154 *
155 * q0 the corresponding exponent of q[0]. Note that the
156 * exponent for q[i] would be q0-24*i.
157 *
158 * PIo2[] double precision array, obtained by cutting pi/2
159 * into 24 bits chunks.
160 *
161 * f[] ipio2[] in floating point
162 *
163 * iq[] integer array by breaking up q[] in 24-bits chunk.
164 *
165 * fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
166 *
167 * ih integer. If >0 it indicates q[] is >= 0.5, hence
168 * it also indicates the *sign* of the result.
169 *
170 */
171
172
173 /*
174 * Constants:
175 * The hexadecimal values are the intended ones for the following
176 * constants. The decimal values may be used, provided that the
177 * compiler will convert from decimal to binary accurately enough
178 * to produce the hexadecimal values shown.
179 */
180
181 #include "mathincl/fdlibm.h"
182
183 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
184
185 static const double PIo2[] = {
186 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
187 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
188 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
189 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
190 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
191 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
192 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
193 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
194 };
195
196 static const double
197 zero = 0.0,
198 one = 1.0,
199 two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
200 twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
201
202 int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2)
203 {
204 int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
205 double z,fw,f[20],fq[20],q[20];
206
207 /* initialize jk*/
208 jk = init_jk[prec];
209 jp = jk;
210
211 /* determine jx,jv,q0, note that 3>q0 */
212 jx = nx-1;
213 jv = (e0-3)/24; if(jv<0) jv=0;
214 q0 = e0-24*(jv+1);
215
216 /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
217 j = jv-jx; m = jx+jk;
218 for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
219
220 /* compute q[0],q[1],...q[jk] */
221 for (i=0;i<=jk;i++) {
222 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
223 }
224
225 jz = jk;
226 recompute:
227 /* distill q[] into iq[] reversingly */
228 for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
229 fw = (double)((int)(twon24* z));
230 iq[i] = (int)(z-two24*fw);
231 z = q[j-1]+fw;
232 }
233
234 /* compute n */
235 z = scalbn(z,q0); /* actual value of z */
236 z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */
237 n = (int) z;
238 z -= (double)n;
239 ih = 0;
240 if(q0>0) { /* need iq[jz-1] to determine n */
241 i = (iq[jz-1]>>(24-q0)); n += i;
242 iq[jz-1] -= i<<(24-q0);
243 ih = iq[jz-1]>>(23-q0);
244 }
245 else if(q0==0) ih = iq[jz-1]>>23;
246 else if(z>=0.5) ih=2;
247
248 if(ih>0) { /* q > 0.5 */
249 n += 1; carry = 0;
250 for(i=0;i<jz ;i++) { /* compute 1-q */
251 j = iq[i];
252 if(carry==0) {
253 if(j!=0) {
254 carry = 1; iq[i] = 0x1000000- j;
255 }
256 } else iq[i] = 0xffffff - j;
257 }
258 if(q0>0) { /* rare case: chance is 1 in 12 */
259 switch(q0) {
260 case 1:
261 iq[jz-1] &= 0x7fffff; break;
262 case 2:
263 iq[jz-1] &= 0x3fffff; break;
264 }
265 }
266 if(ih==2) {
267 z = one - z;
268 if(carry!=0) z -= scalbn(one,q0);
269 }
270 }
271
272 /* check if recomputation is needed */
273 if(z==zero) {
274 j = 0;
275 for (i=jz-1;i>=jk;i--) j |= iq[i];
276 if(j==0) { /* need recomputation */
277 for(k=1;iq[jk-k]==0;k++); /* k = no. of terms needed */
278
279 for(i=jz+1;i<=jz+k;i++) { /* add q[jz+1] to q[jz+k] */
280 f[jx+i] = (double) ipio2[jv+i];
281 for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
282 q[i] = fw;
283 }
284 jz += k;
285 goto recompute;
286 }
287 }
288
289 /* chop off zero terms */
290 if(z==0.0) {
291 jz -= 1; q0 -= 24;
292 while(iq[jz]==0) { jz--; q0-=24;}
293 } else { /* break z into 24-bit if necessary */
294 z = scalbn(z,-q0);
295 if(z>=two24) {
296 fw = (double)((int)(twon24*z));
297 iq[jz] = (int)(z-two24*fw);
298 jz += 1; q0 += 24;
299 iq[jz] = (int) fw;
300 } else iq[jz] = (int) z ;
301 }
302
303 /* convert integer "bit" chunk to floating-point value */
304 fw = scalbn(one,q0);
305 for(i=jz;i>=0;i--) {
306 q[i] = fw*(double)iq[i]; fw*=twon24;
307 }
308
309 /* compute PIo2[0,...,jp]*q[jz,...,0] */
310 for(i=jz;i>=0;i--) {
311 for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
312 fq[jz-i] = fw;
313 }
314
315 /* compress fq[] into y[] */
316 switch(prec) {
317 case 0:
318 fw = 0.0;
319 for (i=jz;i>=0;i--) fw += fq[i];
320 y[0] = (ih==0)? fw: -fw;
321 break;
322 case 1:
323 case 2:
324 fw = 0.0;
325 for (i=jz;i>=0;i--) fw += fq[i];
326 y[0] = (ih==0)? fw: -fw;
327 fw = fq[0]-fw;
328 for (i=1;i<=jz;i++) fw += fq[i];
329 y[1] = (ih==0)? fw: -fw;
330 break;
331 case 3: /* painful */
332 for (i=jz;i>0;i--) {
333 fw = fq[i-1]+fq[i];
334 fq[i] += fq[i-1]-fw;
335 fq[i-1] = fw;
336 }
337 for (i=jz;i>1;i--) {
338 fw = fq[i-1]+fq[i];
339 fq[i] += fq[i-1]-fw;
340 fq[i-1] = fw;
341 }
342 for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
343 if(ih==0) {
344 y[0] = fq[0]; y[1] = fq[1]; y[2] = fw;
345 } else {
346 y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
347 }
348 }
349 return n&7;
350 }
351
352 #endif // ifdef CYGPKG_LIBM
353
354 // EOF k_rem_pio2.c