File: | real/kremp2.c |
Location: | line 329, column 19 |
Description: | The left operand of '-' is a garbage value |
1 | /* @(#)k_rem_pio2.c 1.3 95/01/18 */ | |||
2 | ||||
3 | /* | |||
4 | * Copyright (c) 2015-2017 Carsten Sonne Larsen <cs@innolan.dk> | |||
5 | * All rights reserved. | |||
6 | * | |||
7 | * Redistribution and use in source and binary forms, with or without | |||
8 | * modification, are permitted provided that the following conditions | |||
9 | * are met: | |||
10 | * 1. Redistributions of source code must retain the above copyright | |||
11 | * notice, this list of conditions and the following disclaimer. | |||
12 | * 2. Redistributions in binary form must reproduce the above copyright | |||
13 | * notice, this list of conditions and the following disclaimer in the | |||
14 | * documentation and/or other materials provided with the distribution. | |||
15 | * | |||
16 | * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR | |||
17 | * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES | |||
18 | * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. | |||
19 | * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT, | |||
20 | * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT | |||
21 | * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, | |||
22 | * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY | |||
23 | * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT | |||
24 | * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF | |||
25 | * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. | |||
26 | * | |||
27 | * The origin source code can be obtained from: | |||
28 | * http://www.netlib.org/fdlibm/k_rem_pio2.c | |||
29 | * | |||
30 | */ | |||
31 | ||||
32 | /* | |||
33 | * ==================================================== | |||
34 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | |||
35 | * | |||
36 | * Developed at SunSoft, a Sun Microsystems, Inc. business. | |||
37 | * Permission to use, copy, modify, and distribute this | |||
38 | * software is freely granted, provided that this notice | |||
39 | * is preserved. | |||
40 | * ==================================================== | |||
41 | */ | |||
42 | ||||
43 | #include "prim.h" | |||
44 | #include "math.h" | |||
45 | ||||
46 | /* | |||
47 | * Constants: | |||
48 | * The hexadecimal values are the intended ones for the following | |||
49 | * constants. The decimal values may be used, provided that the | |||
50 | * compiler will convert from decimal to binary accurately enough | |||
51 | * to produce the hexadecimal values shown. | |||
52 | */ | |||
53 | ||||
54 | static const int init_jk[] = {2,3,4,6}; /* initial value for jk */ | |||
55 | ||||
56 | static const double PIo2[] = { | |||
57 | 1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */ | |||
58 | 7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */ | |||
59 | 5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */ | |||
60 | 3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */ | |||
61 | 1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */ | |||
62 | 1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */ | |||
63 | 2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */ | |||
64 | 2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */ | |||
65 | }; | |||
66 | ||||
67 | static const double | |||
68 | zero = 0.0, | |||
69 | one = 1.0, | |||
70 | two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */ | |||
71 | twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */ | |||
72 | ||||
73 | /** | |||
74 | * @brief Kernel reduction function. | |||
75 | * @version 1.4 | |||
76 | * @date 96/03/07 | |||
77 | * @details | |||
78 | * <pre> | |||
79 | * __kernel_rem_pio2(x,y,e0,nx,prec,ipio2) | |||
80 | * double x[],y[]; int e0,nx,prec; int ipio2[]; | |||
81 | * | |||
82 | * __kernel_rem_pio2 return the last three digits of N with | |||
83 | * y = x - N*pi/2 | |||
84 | * so that |y| < pi/2. | |||
85 | * | |||
86 | * The method is to compute the integer (mod 8) and fraction parts of | |||
87 | * (2/pi)*x without doing the full multiplication. In general we | |||
88 | * skip the part of the product that are known to be a huge integer ( | |||
89 | * more accurately, = 0 mod 8 ). Thus the number of operations are | |||
90 | * independent of the exponent of the input. | |||
91 | * | |||
92 | * (2/pi) is represented by an array of 24-bit integers in ipio2[]. | |||
93 | * | |||
94 | * Input parameters: | |||
95 | * x[] The input value (must be positive) is broken into nx | |||
96 | * pieces of 24-bit integers in double precision format. | |||
97 | * x[i] will be the i-th 24 bit of x. The scaled exponent | |||
98 | * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0 | |||
99 | * match x's up to 24 bits. | |||
100 | * | |||
101 | * Example of breaking a double positive z into x[0]+x[1]+x[2]: | |||
102 | * e0 = ilogb(z)-23 | |||
103 | * z = scalbn(z,-e0) | |||
104 | * for i = 0,1,2 | |||
105 | * x[i] = floor(z) | |||
106 | * z = (z-x[i])*2**24 | |||
107 | * | |||
108 | * | |||
109 | * y[] ouput result in an array of double precision numbers. | |||
110 | * The dimension of y[] is: | |||
111 | * 24-bit precision 1 | |||
112 | * 53-bit precision 2 | |||
113 | * 64-bit precision 2 | |||
114 | * 113-bit precision 3 | |||
115 | * The actual value is the sum of them. Thus for 113-bit | |||
116 | * precison, one may have to do something like: | |||
117 | * | |||
118 | * long double t,w,r_head, r_tail; | |||
119 | * t = (long double)y[2] + (long double)y[1]; | |||
120 | * w = (long double)y[0]; | |||
121 | * r_head = t+w; | |||
122 | * r_tail = w - (r_head - t); | |||
123 | * | |||
124 | * e0 The exponent of x[0] | |||
125 | * | |||
126 | * nx dimension of x[] | |||
127 | * | |||
128 | * prec an integer indicating the precision: | |||
129 | * 0 24 bits (single) | |||
130 | * 1 53 bits (double) | |||
131 | * 2 64 bits (extended) | |||
132 | * 3 113 bits (quad) | |||
133 | * | |||
134 | * ipio2[] | |||
135 | * integer array, contains the (24*i)-th to (24*i+23)-th | |||
136 | * bit of 2/pi after binary point. The corresponding | |||
137 | * floating value is | |||
138 | * | |||
139 | * ipio2[i] * 2^(-24(i+1)). | |||
140 | * | |||
141 | * External function: | |||
142 | * double scalbn(), floor(); | |||
143 | * | |||
144 | * | |||
145 | * Here is the description of some local variables: | |||
146 | * | |||
147 | * jk jk+1 is the initial number of terms of ipio2[] needed | |||
148 | * in the computation. The recommended value is 2,3,4, | |||
149 | * 6 for single, double, extended,and quad. | |||
150 | * | |||
151 | * jz local integer variable indicating the number of | |||
152 | * terms of ipio2[] used. | |||
153 | * | |||
154 | * jx nx - 1 | |||
155 | * | |||
156 | * jv index for pointing to the suitable ipio2[] for the | |||
157 | * computation. In general, we want | |||
158 | * ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8 | |||
159 | * is an integer. Thus | |||
160 | * e0-3-24*jv >= 0 or (e0-3)/24 >= jv | |||
161 | * Hence jv = max(0,(e0-3)/24). | |||
162 | * | |||
163 | * jp jp+1 is the number of terms in PIo2[] needed, jp = jk. | |||
164 | * | |||
165 | * q[] double array with integral value, representing the | |||
166 | * 24-bits chunk of the product of x and 2/pi. | |||
167 | * | |||
168 | * q0 the corresponding exponent of q[0]. Note that the | |||
169 | * exponent for q[i] would be q0-24*i. | |||
170 | * | |||
171 | * PIo2[] double precision array, obtained by cutting pi/2 | |||
172 | * into 24 bits chunks. | |||
173 | * | |||
174 | * f[] ipio2[] in floating point | |||
175 | * | |||
176 | * iq[] integer array by breaking up q[] in 24-bits chunk. | |||
177 | * | |||
178 | * fq[] final product of x*(2/pi) in fq[0],..,fq[jk] | |||
179 | * | |||
180 | * ih integer. If >0 it indicates q[] is >= 0.5, hence | |||
181 | * it also indicates the *sign* of the result. | |||
182 | * | |||
183 | * </pre> | |||
184 | * @copyright Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. | |||
185 | * @license Developed at SunSoft, a Sun Microsystems, Inc. business. Permission | |||
186 | * to use, copy, modify, and distribute this software is freely granted, | |||
187 | * provided that this notice is preserved. | |||
188 | */ | |||
189 | ||||
190 | int __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec, const int *ipio2) | |||
191 | { | |||
192 | int jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih; | |||
193 | double z,fw,f[20],fq[20],q[20]; | |||
194 | ||||
195 | /* initialize jk*/ | |||
196 | jk = init_jk[prec]; | |||
197 | jp = jk; | |||
198 | ||||
199 | /* determine jx,jv,q0, note that 3>q0 */ | |||
200 | jx = nx-1; | |||
201 | jv = (e0-3)/24; | |||
202 | if(jv<0) jv=0; | |||
| ||||
203 | q0 = e0-24*(jv+1); | |||
204 | ||||
205 | /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */ | |||
206 | j = jv-jx; | |||
207 | m = jx+jk; | |||
208 | for(i=0; i<=m; i++,j++) f[i] = (j<0)? zero : (double) ipio2[j]; | |||
209 | ||||
210 | /* compute q[0],q[1],...q[jk] */ | |||
211 | for (i=0; i<=jk; i++) { | |||
212 | for(j=0,fw=0.0; j<=jx; j++) fw += x[j]*f[jx+i-j]; | |||
213 | q[i] = fw; | |||
214 | } | |||
215 | ||||
216 | jz = jk; | |||
217 | recompute: | |||
218 | /* distill q[] into iq[] reversingly */ | |||
219 | for(i=0,j=jz,z=q[jz]; j>0; i++,j--) { | |||
220 | fw = (double)((int)(twon24* z)); | |||
221 | iq[i] = (int)(z-two24*fw); | |||
222 | z = q[j-1]+fw; | |||
223 | } | |||
224 | ||||
225 | /* compute n */ | |||
226 | z = scalbn(z,q0); /* actual value of z */ | |||
227 | z -= 8.0*floor(z*0.125); /* trim off integer >= 8 */ | |||
228 | n = (int) z; | |||
229 | z -= (double)n; | |||
230 | ih = 0; | |||
231 | if(q0>0) { /* need iq[jz-1] to determine n */ | |||
232 | i = (iq[jz-1]>>(24-q0)); | |||
233 | n += i; | |||
234 | iq[jz-1] -= i<<(24-q0); | |||
235 | ih = iq[jz-1]>>(23-q0); | |||
236 | } | |||
237 | else if(q0==0) ih = iq[jz-1]>>23; | |||
238 | else if(z>=0.5) ih=2; | |||
239 | ||||
240 | if(ih>0) { /* q > 0.5 */ | |||
241 | n += 1; | |||
242 | carry = 0; | |||
243 | for(i=0; i<jz ; i++) { /* compute 1-q */ | |||
244 | j = iq[i]; | |||
245 | if(carry==0) { | |||
246 | if(j!=0) { | |||
247 | carry = 1; | |||
248 | iq[i] = 0x1000000- j; | |||
249 | } | |||
250 | } else iq[i] = 0xffffff - j; | |||
251 | } | |||
252 | if(q0>0) { /* rare case: chance is 1 in 12 */ | |||
253 | switch(q0) { | |||
254 | case 1: | |||
255 | iq[jz-1] &= 0x7fffff; | |||
256 | break; | |||
257 | case 2: | |||
258 | iq[jz-1] &= 0x3fffff; | |||
259 | break; | |||
260 | } | |||
261 | } | |||
262 | if(ih==2) { | |||
263 | z = one - z; | |||
264 | if(carry!=0) z -= scalbn(one,q0); | |||
265 | } | |||
266 | } | |||
267 | ||||
268 | /* check if recomputation is needed */ | |||
269 | if(z==zero) { | |||
270 | j = 0; | |||
271 | for (i=jz-1; i>=jk; i--) j |= iq[i]; | |||
272 | if(j==0) { /* need recomputation */ | |||
273 | for(k=1; iq[jk-k]==0; k++); /* k = no. of terms needed */ | |||
274 | ||||
275 | for(i=jz+1; i<=jz+k; i++) { /* add q[jz+1] to q[jz+k] */ | |||
276 | f[jx+i] = (double) ipio2[jv+i]; | |||
277 | for(j=0,fw=0.0; j<=jx; j++) fw += x[j]*f[jx+i-j]; | |||
278 | q[i] = fw; | |||
279 | } | |||
280 | jz += k; | |||
281 | goto recompute; | |||
282 | } | |||
283 | } | |||
284 | ||||
285 | /* chop off zero terms */ | |||
286 | if(z==0.0) { | |||
287 | jz -= 1; | |||
288 | q0 -= 24; | |||
289 | while(iq[jz]==0) { | |||
290 | jz--; | |||
291 | q0-=24; | |||
292 | } | |||
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; | |||
299 | q0 += 24; | |||
300 | iq[jz] = (int) fw; | |||
301 | } else iq[jz] = (int) z ; | |||
302 | } | |||
303 | ||||
304 | /* convert integer "bit" chunk to floating-point value */ | |||
305 | fw = scalbn(one,q0); | |||
306 | for(i=jz; i>=0; i--) { | |||
307 | q[i] = fw*(double)iq[i]; | |||
308 | fw*=twon24; | |||
309 | } | |||
310 | ||||
311 | /* compute PIo2[0,...,jp]*q[jz,...,0] */ | |||
312 | for(i=jz; i>=0; i--) { | |||
313 | for(fw=0.0,k=0; k<=jp&&k<=jz-i; k++) fw += PIo2[k]*q[i+k]; | |||
314 | fq[jz-i] = fw; | |||
315 | } | |||
316 | ||||
317 | /* compress fq[] into y[] */ | |||
318 | switch(prec) { | |||
319 | case 0: | |||
320 | fw = 0.0; | |||
321 | for (i=jz; i>=0; i--) fw += fq[i]; | |||
322 | y[0] = (ih==0)? fw: -fw; | |||
323 | break; | |||
324 | case 1: | |||
325 | case 2: | |||
326 | fw = 0.0; | |||
327 | for (i=jz; i>=0; i--) fw += fq[i]; | |||
328 | y[0] = (ih==0)? fw: -fw; | |||
329 | fw = fq[0]-fw; | |||
| ||||
330 | for (i=1; i<=jz; i++) fw += fq[i]; | |||
331 | y[1] = (ih==0)? fw: -fw; | |||
332 | break; | |||
333 | case 3: /* painful */ | |||
334 | for (i=jz; i>0; i--) { | |||
335 | fw = fq[i-1]+fq[i]; | |||
336 | fq[i] += fq[i-1]-fw; | |||
337 | fq[i-1] = fw; | |||
338 | } | |||
339 | for (i=jz; i>1; i--) { | |||
340 | fw = fq[i-1]+fq[i]; | |||
341 | fq[i] += fq[i-1]-fw; | |||
342 | fq[i-1] = fw; | |||
343 | } | |||
344 | for (fw=0.0,i=jz; i>=2; i--) fw += fq[i]; | |||
345 | if(ih==0) { | |||
346 | y[0] = fq[0]; | |||
347 | y[1] = fq[1]; | |||
348 | y[2] = fw; | |||
349 | } else { | |||
350 | y[0] = -fq[0]; | |||
351 | y[1] = -fq[1]; | |||
352 | y[2] = -fw; | |||
353 | } | |||
354 | } | |||
355 | return n&7; | |||
356 | } |