### Bug Summary

File: | real/kremp2.c
Location: | line 351, 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 |
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;

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

21 |

Assuming 'j' is <= 0

22 |

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

23 |

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

24 |

Taking false branch
238 else if(z>=0.5) ih=2;

12 |

Taking false branch

25 |

Taking true branch
239
240 if(ih>0) { /* q > 0.5 */

13 |

Taking false branch

26 |

Taking true branch
241 n += 1;
242 carry = 0;
243 for(i=0; i<jz ; i++) { /* compute 1-q */

27 |

Loop condition is false. Execution continues on line 252
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 */

28 |

Taking false branch
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) {

29 |

Taking true branch
263 z = one - z;
264 if(carry!=0) z -= scalbn(one,q0);

30 |

Taking false branch
265 }
266 }
267
268 /* check if recomputation is needed */
269 if(z==zero) {

14 |

Taking true branch

31 |

Taking false branch
270 j = 0;
271 for (i=jz-1; i>=jk; i--) j |= iq[i];

15 |

Loop condition is false. Execution continues on line 272
272 if(j==0) { /* need recomputation */

16 |

Taking true branch
273 for(k=1; iq[jk-k]==0; k++); /* k = no. of terms needed */

17 |

Loop condition is true. Entering loop body

18 |

Loop condition is false. Execution continues on line 275
274
275 for(i=jz+1; i<=jz+k; i++) { /* add q[jz+1] to q[jz+k] */

19 |

Loop condition is false. Execution continues on line 280
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;

20 |

Control jumps to line 219
282 }
283 }
284
285 /* chop off zero terms */
286 if(z==0.0) {

32 |

Taking false branch
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) {

33 |

Taking false branch
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--) {

34 |

Assuming 'i' is >= 0

35 |

Loop condition is true. Entering loop body

36 |

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--) {

37 |

Loop condition is true. Entering loop body

38 |

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) {

39 |

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--) {

40 |

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--) {

41 |

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];

42 |

Loop condition is false. Execution continues on line 345
345 if(ih==0) {

43 |

Taking false branch
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];

44 |

Assigned value is garbage or undefined
352 y[2] = -fw;
353 }
354 }
355 return n&7;
356}