Bug Summary

File:real/kremp2.c
Location:line 346, column 18
Description:Assigned value is garbage or undefined

Annotated Source Code

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
54static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
55
56static 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
67static const double
68zero = 0.0,
69one = 1.0,
70two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
71twon24 = 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
190int __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;
1
Assuming 'jv' is >= 0
→
2
←
Taking false branch
→
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];
3
←
Assuming 'i' is > 'm'
→
4
←
Loop condition is false. Execution continues on line 211
→
209
210 /* compute q[0],q[1],...q[jk] */
211 for (i=0; i<=jk; i++) {
5
←
Assuming 'i' is > 'jk'
→
6
←
Loop condition is false. Execution continues on line 216
→
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;
217recompute:
218 /* distill q[] into iq[] reversingly */
219 for(i=0,j=jz,z=q[jz]; j>0; i++,j--) {
7
←
Loop condition is false. Execution continues on line 226
→
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 */
8
←
Assuming 'q0' is <= 0
→
9
←
Taking false branch
→
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;
10
←
Assuming 'q0' is not equal to 0
→
11
←
Taking false branch
→
238 else if(z>=0.5) ih=2;
12
←
Taking false branch
→
239
240 if(ih>0) { /* q > 0.5 */
13
←
Taking false branch
→
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) {
14
←
Taking false branch
→
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) {
15
←
Taking true branch
→
287 jz -= 1;
288 q0 -= 24;
289 while(iq[jz]==0) {
16
←
Loop condition is false. Execution continues on line 305
→
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--) {
17
←
Assuming 'i' is < 0
→
18
←
Loop condition is false. Execution continues on line 312
→
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--) {
19
←
Loop condition is false. Execution continues on line 318
→
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) {
20
←
Control jumps to 'case 3:' at line 333
→
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--) {
21
←
Loop condition is false. Execution continues on line 339
→
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--) {
22
←
Loop condition is false. Execution continues on line 344
→
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];
23
←
Loop condition is false. Execution continues on line 345
→
345 if(ih==0) {
24
←
Taking true branch
→
346 y[0] = fq[0];
25
←
Assigned value is garbage or undefined
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}