6
6
*
7
7
* Developed at SunSoft, a Sun Microsystems, Inc. business.
8
8
* Permission to use, copy, modify, and distribute this
9
- * software is freely granted, provided that this notice
9
+ * software is freely granted, provided that this notice
10
10
* is preserved.
11
11
* ====================================================
12
12
*/
16
16
/*
17
17
* __kernel_rem_pio2(x,y,e0,nx,prec)
18
18
* double x[],y[]; int e0,nx,prec;
19
- *
20
- * __kernel_rem_pio2 return the last three digits of N with
19
+ *
20
+ * __kernel_rem_pio2 return the last three digits of N with
21
21
* y = x - N*pi/2
22
22
* so that |y| < pi/2.
23
23
*
24
- * The method is to compute the integer (mod 8) and fraction parts of
24
+ * The method is to compute the integer (mod 8) and fraction parts of
25
25
* (2/pi)*x without doing the full multiplication. In general we
26
26
* skip the part of the product that are known to be a huge integer (
27
27
* more accurately, = 0 mod 8 ). Thus the number of operations are
30
30
* (2/pi) is represented by an array of 24-bit integers in ipio2[].
31
31
*
32
32
* Input parameters:
33
- * x[] The input value (must be positive) is broken into nx
33
+ * x[] The input value (must be positive) is broken into nx
34
34
* pieces of 24-bit integers in double precision format.
35
- * x[i] will be the i-th 24 bit of x. The scaled exponent
36
- * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
35
+ * x[i] will be the i-th 24 bit of x. The scaled exponent
36
+ * of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
37
37
* match x's up to 24 bits.
38
38
*
39
39
* Example of breaking a double positive z into x[0]+x[1]+x[2]:
85
85
* the fraction part may be lost to cancelation before we
86
86
* recompute.)
87
87
*
88
- * jz local integer variable indicating the number of
89
- * terms of ipio2[] used.
88
+ * jz local integer variable indicating the number of
89
+ * terms of ipio2[] used.
90
90
*
91
91
* jx nx - 1
92
92
*
106
106
* exponent for q[i] would be q0-24*i.
107
107
*
108
108
* PIo2[] double precision array, obtained by cutting pi/2
109
- * into 24 bits chunks.
109
+ * into 24 bits chunks.
110
110
*
111
- * f[] ipio2[] in floating point
111
+ * f[] ipio2[] in floating point
112
112
*
113
113
* iq[] integer array by breaking up q[] in 24-bits chunk.
114
114
*
122
122
123
123
/*
124
124
* Constants:
125
- * The hexadecimal values are the intended ones for the following
126
- * constants. The decimal values may be used, provided that the
127
- * compiler will convert from decimal to binary accurately enough
125
+ * The hexadecimal values are the intended ones for the following
126
+ * constants. The decimal values may be used, provided that the
127
+ * compiler will convert from decimal to binary accurately enough
128
128
* to produce the hexadecimal values shown.
129
129
*/
130
130
@@ -142,8 +142,8 @@ static const int init_jk[] = {3,4,4,6}; /* initial value for jk */
142
142
/*
143
143
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
144
144
*
145
- * integer array, contains the (24*i)-th to (24*i+23)-th
146
- * bit of 2/pi after binary point. The corresponding
145
+ * integer array, contains the (24*i)-th to (24*i+23)-th
146
+ * bit of 2/pi after binary point. The corresponding
147
147
* floating value is
148
148
*
149
149
* ipio2[i] * 2^(-24(i+1)).
@@ -152,17 +152,17 @@ static const int init_jk[] = {3,4,4,6}; /* initial value for jk */
152
152
* For quad precision (e0 <= 16360, jk = 6), this is 686.
153
153
*/
154
154
static const int32_t ipio2 [] = {
155
- 0xA2F983 , 0x6E4E44 , 0x1529FC , 0x2757D1 , 0xF534DD , 0xC0DB62 ,
156
- 0x95993C , 0x439041 , 0xFE5163 , 0xABDEBB , 0xC561B7 , 0x246E3A ,
157
- 0x424DD2 , 0xE00649 , 0x2EEA09 , 0xD1921C , 0xFE1DEB , 0x1CB129 ,
158
- 0xA73EE8 , 0x8235F5 , 0x2EBB44 , 0x84E99C , 0x7026B4 , 0x5F7E41 ,
159
- 0x3991D6 , 0x398353 , 0x39F49C , 0x845F8B , 0xBDF928 , 0x3B1FF8 ,
160
- 0x97FFDE , 0x05980F , 0xEF2F11 , 0x8B5A0A , 0x6D1F6D , 0x367ECF ,
161
- 0x27CB09 , 0xB74F46 , 0x3F669E , 0x5FEA2D , 0x7527BA , 0xC7EBE5 ,
162
- 0xF17B3D , 0x0739F7 , 0x8A5292 , 0xEA6BFB , 0x5FB11F , 0x8D5D08 ,
163
- 0x560330 , 0x46FC7B , 0x6BABF0 , 0xCFBC20 , 0x9AF436 , 0x1DA9E3 ,
164
- 0x91615E , 0xE61B08 , 0x659985 , 0x5F14A0 , 0x68408D , 0xFFD880 ,
165
- 0x4D7327 , 0x310606 , 0x1556CA , 0x73A8C9 , 0x60E27B , 0xC08C6B ,
155
+ 0xA2F983 , 0x6E4E44 , 0x1529FC , 0x2757D1 , 0xF534DD , 0xC0DB62 ,
156
+ 0x95993C , 0x439041 , 0xFE5163 , 0xABDEBB , 0xC561B7 , 0x246E3A ,
157
+ 0x424DD2 , 0xE00649 , 0x2EEA09 , 0xD1921C , 0xFE1DEB , 0x1CB129 ,
158
+ 0xA73EE8 , 0x8235F5 , 0x2EBB44 , 0x84E99C , 0x7026B4 , 0x5F7E41 ,
159
+ 0x3991D6 , 0x398353 , 0x39F49C , 0x845F8B , 0xBDF928 , 0x3B1FF8 ,
160
+ 0x97FFDE , 0x05980F , 0xEF2F11 , 0x8B5A0A , 0x6D1F6D , 0x367ECF ,
161
+ 0x27CB09 , 0xB74F46 , 0x3F669E , 0x5FEA2D , 0x7527BA , 0xC7EBE5 ,
162
+ 0xF17B3D , 0x0739F7 , 0x8A5292 , 0xEA6BFB , 0x5FB11F , 0x8D5D08 ,
163
+ 0x560330 , 0x46FC7B , 0x6BABF0 , 0xCFBC20 , 0x9AF436 , 0x1DA9E3 ,
164
+ 0x91615E , 0xE61B08 , 0x659985 , 0x5F14A0 , 0x68408D , 0xFFD880 ,
165
+ 0x4D7327 , 0x310606 , 0x1556CA , 0x73A8C9 , 0x60E27B , 0xC08C6B ,
166
166
167
167
#if LDBL_MAX_EXP > 1024
168
168
#if LDBL_MAX_EXP > 16384
@@ -287,7 +287,7 @@ static const double PIo2[] = {
287
287
2.16741683877804819444e-51 , /* 0x3569F31D, 0x00000000 */
288
288
};
289
289
290
- static const double
290
+ static const double
291
291
zero = 0.0 ,
292
292
one = 1.0 ,
293
293
two24 = 1.67772160000000000000e+07 , /* 0x41700000, 0x00000000 */
@@ -336,7 +336,7 @@ __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
336
336
i = (iq [jz - 1 ]>>(24 - q0 )); n += i ;
337
337
iq [jz - 1 ] -= i <<(24 - q0 );
338
338
ih = iq [jz - 1 ]>>(23 - q0 );
339
- }
339
+ }
340
340
else if (q0 == 0 ) ih = iq [jz - 1 ]>>23 ;
341
341
else if (z >=0.5 ) ih = 2 ;
342
342
@@ -387,7 +387,7 @@ __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
387
387
while (iq [jz ]== 0 ) { jz -- ; q0 -= 24 ;}
388
388
} else { /* break z into 24-bit if necessary */
389
389
z = scalbn (z ,- q0 );
390
- if (z >=two24 ) {
390
+ if (z >=two24 ) {
391
391
fw = (double )((int32_t )(twon24 * z ));
392
392
iq [jz ] = (int32_t )(z - two24 * fw );
393
393
jz += 1 ; q0 += 24 ;
@@ -412,30 +412,30 @@ __kernel_rem_pio2(double *x, double *y, int e0, int nx, int prec)
412
412
case 0 :
413
413
fw = 0.0 ;
414
414
for (i = jz ;i >=0 ;i -- ) fw += fq [i ];
415
- y [0 ] = (ih == 0 )? fw : - fw ;
415
+ y [0 ] = (ih == 0 )? fw : - fw ;
416
416
break ;
417
417
case 1 :
418
418
case 2 :
419
419
fw = 0.0 ;
420
- for (i = jz ;i >=0 ;i -- ) fw += fq [i ];
420
+ for (i = jz ;i >=0 ;i -- ) fw += fq [i ];
421
421
STRICT_ASSIGN (double ,fw ,fw );
422
- y [0 ] = (ih == 0 )? fw : - fw ;
422
+ y [0 ] = (ih == 0 )? fw : - fw ;
423
423
fw = fq [0 ]- fw ;
424
424
for (i = 1 ;i <=jz ;i ++ ) fw += fq [i ];
425
- y [1 ] = (ih == 0 )? fw : - fw ;
425
+ y [1 ] = (ih == 0 )? fw : - fw ;
426
426
break ;
427
427
case 3 : /* painful */
428
428
for (i = jz ;i > 0 ;i -- ) {
429
- fw = fq [i - 1 ]+ fq [i ];
429
+ fw = fq [i - 1 ]+ fq [i ];
430
430
fq [i ] += fq [i - 1 ]- fw ;
431
431
fq [i - 1 ] = fw ;
432
432
}
433
433
for (i = jz ;i > 1 ;i -- ) {
434
- fw = fq [i - 1 ]+ fq [i ];
434
+ fw = fq [i - 1 ]+ fq [i ];
435
435
fq [i ] += fq [i - 1 ]- fw ;
436
436
fq [i - 1 ] = fw ;
437
437
}
438
- for (fw = 0.0 ,i = jz ;i >=2 ;i -- ) fw += fq [i ];
438
+ for (fw = 0.0 ,i = jz ;i >=2 ;i -- ) fw += fq [i ];
439
439
if (ih == 0 ) {
440
440
y [0 ] = fq [0 ]; y [1 ] = fq [1 ]; y [2 ] = fw ;
441
441
} else {
0 commit comments