1 #include "quadmath-imp.h"
2 #include <math.h>
3 
4 
5 /* @(#)k_rem_pio2.c 5.1 93/09/24 */
6 /*
7  * ====================================================
8  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
9  *
10  * Developed at SunPro, a Sun Microsystems, Inc. business.
11  * Permission to use, copy, modify, and distribute this
12  * software is freely granted, provided that this notice
13  * is preserved.
14  * ====================================================
15  */
16 
17 /*
18  * __quadmath_kernel_rem_pio2 (x,y,e0,nx,prec,ipio2)
19  * double x[],y[]; int e0,nx,prec; int ipio2[];
20  *
21  * __quadmath_kernel_rem_pio2  return the last three digits of N with
22  *		y = x - N*pi/2
23  * so that |y| < pi/2.
24  *
25  * The method is to compute the integer (mod 8) and fraction parts of
26  * (2/pi)*x without doing the full multiplication. In general we
27  * skip the part of the product that are known to be a huge integer (
28  * more accurately, = 0 mod 8 ). Thus the number of operations are
29  * independent of the exponent of the input.
30  *
31  * (2/pi) is represented by an array of 24-bit integers in ipio2[].
32  *
33  * Input parameters:
34  * 	x[]	The input value (must be positive) is broken into nx
35  *		pieces of 24-bit integers in double precision format.
36  *		x[i] will be the i-th 24 bit of x. The scaled exponent
37  *		of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
38  *		match x's up to 24 bits.
39  *
40  *		Example of breaking a double positive z into x[0]+x[1]+x[2]:
41  *			e0 = ilogb(z)-23
42  *			z  = scalbn(z,-e0)
43  *		for i = 0,1,2
44  *			x[i] = floor(z)
45  *			z    = (z-x[i])*2**24
46  *
47  *
48  *	y[]	ouput result in an array of double precision numbers.
49  *		The dimension of y[] is:
50  *			24-bit  precision	1
51  *			53-bit  precision	2
52  *			64-bit  precision	2
53  *			113-bit precision	3
54  *		The actual value is the sum of them. Thus for 113-bit
55  *		precision, one may have to do something like:
56  *
57  *		long double t,w,r_head, r_tail;
58  *		t = (long double)y[2] + (long double)y[1];
59  *		w = (long double)y[0];
60  *		r_head = t+w;
61  *		r_tail = w - (r_head - t);
62  *
63  *	e0	The exponent of x[0]
64  *
65  *	nx	dimension of x[]
66  *
67  *  	prec	an integer indicating the precision:
68  *			0	24  bits (single)
69  *			1	53  bits (double)
70  *			2	64  bits (extended)
71  *			3	113 bits (quad)
72  *
73  *	ipio2[]
74  *		integer array, contains the (24*i)-th to (24*i+23)-th
75  *		bit of 2/pi after binary point. The corresponding
76  *		floating value is
77  *
78  *			ipio2[i] * 2^(-24(i+1)).
79  *
80  * External function:
81  *	double scalbn(), floor();
82  *
83  *
84  * Here is the description of some local variables:
85  *
86  * 	jk	jk+1 is the initial number of terms of ipio2[] needed
87  *		in the computation. The recommended value is 2,3,4,
88  *		6 for single, double, extended,and quad.
89  *
90  * 	jz	local integer variable indicating the number of
91  *		terms of ipio2[] used.
92  *
93  *	jx	nx - 1
94  *
95  *	jv	index for pointing to the suitable ipio2[] for the
96  *		computation. In general, we want
97  *			( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
98  *		is an integer. Thus
99  *			e0-3-24*jv >= 0 or (e0-3)/24 >= jv
100  *		Hence jv = max(0,(e0-3)/24).
101  *
102  *	jp	jp+1 is the number of terms in PIo2[] needed, jp = jk.
103  *
104  * 	q[]	double array with integral value, representing the
105  *		24-bits chunk of the product of x and 2/pi.
106  *
107  *	q0	the corresponding exponent of q[0]. Note that the
108  *		exponent for q[i] would be q0-24*i.
109  *
110  *	PIo2[]	double precision array, obtained by cutting pi/2
111  *		into 24 bits chunks.
112  *
113  *	f[]	ipio2[] in floating point
114  *
115  *	iq[]	integer array by breaking up q[] in 24-bits chunk.
116  *
117  *	fq[]	final product of x*(2/pi) in fq[0],..,fq[jk]
118  *
119  *	ih	integer. If >0 it indicates q[] is >= 0.5, hence
120  *		it also indicates the *sign* of the result.
121  *
122  */
123 
124 /*
125  * Constants:
126  * The hexadecimal values are the intended ones for the following
127  * constants. The decimal values may be used, provided that the
128  * compiler will convert from decimal to binary accurately enough
129  * to produce the hexadecimal values shown.
130  */
131 
132 
133 static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
134 
135 static const double PIo2[] = {
136   1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
137   7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
138   5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
139   3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
140   1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
141   1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
142   2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
143   2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
144 };
145 
146 static const double
147   zero   = 0.0,
148   one    = 1.0,
149   two24   =  1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
150   twon24  =  5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
151 
152 
153 static int
__quadmath_kernel_rem_pio2(double * x,double * y,int e0,int nx,int prec,const int32_t * ipio2)154 __quadmath_kernel_rem_pio2 (double *x, double *y, int e0, int nx, int prec, const int32_t *ipio2)
155 {
156 	int32_t jz,jx,jv,jp,jk,carry,n,iq[20],i,j,k,m,q0,ih;
157 	double z,fw,f[20],fq[20],q[20];
158 
159     /* initialize jk*/
160 	jk = init_jk[prec];
161 	jp = jk;
162 
163     /* determine jx,jv,q0, note that 3>q0 */
164 	jx =  nx-1;
165 	jv = (e0-3)/24; if(jv<0) jv=0;
166 	q0 =  e0-24*(jv+1);
167 
168     /* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
169 	j = jv-jx; m = jx+jk;
170 	for(i=0;i<=m;i++,j++) f[i] = (j<0)? zero : (double) ipio2[j];
171 
172     /* compute q[0],q[1],...q[jk] */
173 	for (i=0;i<=jk;i++) {
174 	    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j]; q[i] = fw;
175 	}
176 
177 	jz = jk;
178 recompute:
179     /* distill q[] into iq[] reversingly */
180 	for(i=0,j=jz,z=q[jz];j>0;i++,j--) {
181 	    fw    =  (double)((int32_t)(twon24* z));
182 	    iq[i] =  (int32_t)(z-two24*fw);
183 	    z     =  q[j-1]+fw;
184 	}
185 
186     /* compute n */
187 	z  = scalbn(z,q0);		/* actual value of z */
188 	z -= 8.0*floor(z*0.125);		/* trim off integer >= 8 */
189 	n  = (int32_t) z;
190 	z -= (double)n;
191 	ih = 0;
192 	if(q0>0) {	/* need iq[jz-1] to determine n */
193 	    i  = (iq[jz-1]>>(24-q0)); n += i;
194 	    iq[jz-1] -= i<<(24-q0);
195 	    ih = iq[jz-1]>>(23-q0);
196 	}
197 	else if(q0==0) ih = iq[jz-1]>>23;
198 	else if(z>=0.5) ih=2;
199 
200 	if(ih>0) {	/* q > 0.5 */
201 	    n += 1; carry = 0;
202 	    for(i=0;i<jz ;i++) {	/* compute 1-q */
203 		j = iq[i];
204 		if(carry==0) {
205 		    if(j!=0) {
206 			carry = 1; iq[i] = 0x1000000- j;
207 		    }
208 		} else  iq[i] = 0xffffff - j;
209 	    }
210 	    if(q0>0) {		/* rare case: chance is 1 in 12 */
211 	        switch(q0) {
212 	        case 1:
213 	    	   iq[jz-1] &= 0x7fffff; break;
214 	    	case 2:
215 	    	   iq[jz-1] &= 0x3fffff; break;
216 	        }
217 	    }
218 	    if(ih==2) {
219 		z = one - z;
220 		if(carry!=0) z -= scalbn(one,q0);
221 	    }
222 	}
223 
224     /* check if recomputation is needed */
225 	if(z==zero) {
226 	    j = 0;
227 	    for (i=jz-1;i>=jk;i--) j |= iq[i];
228 	    if(j==0) { /* need recomputation */
229 		for(k=1;iq[jk-k]==0;k++);   /* k = no. of terms needed */
230 
231 		for(i=jz+1;i<=jz+k;i++) {   /* add q[jz+1] to q[jz+k] */
232 		    f[jx+i] = (double) ipio2[jv+i];
233 		    for(j=0,fw=0.0;j<=jx;j++) fw += x[j]*f[jx+i-j];
234 		    q[i] = fw;
235 		}
236 		jz += k;
237 		goto recompute;
238 	    }
239 	}
240 
241     /* chop off zero terms */
242 	if(z==0.0) {
243 	    jz -= 1; q0 -= 24;
244 	    while(iq[jz]==0) { jz--; q0-=24;}
245 	} else { /* break z into 24-bit if necessary */
246 	    z = scalbn(z,-q0);
247 	    if(z>=two24) {
248 		fw = (double)((int32_t)(twon24*z));
249 		iq[jz] = (int32_t)(z-two24*fw);
250 		jz += 1; q0 += 24;
251 		iq[jz] = (int32_t) fw;
252 	    } else iq[jz] = (int32_t) z ;
253 	}
254 
255     /* convert integer "bit" chunk to floating-point value */
256 	fw = scalbn(one,q0);
257 	for(i=jz;i>=0;i--) {
258 	    q[i] = fw*(double)iq[i]; fw*=twon24;
259 	}
260 
261     /* compute PIo2[0,...,jp]*q[jz,...,0] */
262 	for(i=jz;i>=0;i--) {
263 	    for(fw=0.0,k=0;k<=jp&&k<=jz-i;k++) fw += PIo2[k]*q[i+k];
264 	    fq[jz-i] = fw;
265 	}
266 
267     /* compress fq[] into y[] */
268 	switch(prec) {
269 	    case 0:
270 		fw = 0.0;
271 		for (i=jz;i>=0;i--) fw += fq[i];
272 		y[0] = (ih==0)? fw: -fw;
273 		break;
274 	    case 1:
275 	    case 2:
276 		fw = 0.0;
277 		for (i=jz;i>=0;i--) fw += fq[i];
278 		y[0] = (ih==0)? fw: -fw;
279 		fw = fq[0]-fw;
280 		for (i=1;i<=jz;i++) fw += fq[i];
281 		y[1] = (ih==0)? fw: -fw;
282 		break;
283 	    case 3:	/* painful */
284 		for (i=jz;i>0;i--) {
285 #if __FLT_EVAL_METHOD__ != 0
286 		    volatile
287 #endif
288 		    double fv = (double)(fq[i-1]+fq[i]);
289 		    fq[i]  += fq[i-1]-fv;
290 		    fq[i-1] = fv;
291 		}
292 		for (i=jz;i>1;i--) {
293 #if __FLT_EVAL_METHOD__ != 0
294 		    volatile
295 #endif
296 		    double fv = (double)(fq[i-1]+fq[i]);
297 		    fq[i]  += fq[i-1]-fv;
298 		    fq[i-1] = fv;
299 		}
300 		for (fw=0.0,i=jz;i>=2;i--) fw += fq[i];
301 		if(ih==0) {
302 		    y[0] =  fq[0]; y[1] =  fq[1]; y[2] =  fw;
303 		} else {
304 		    y[0] = -fq[0]; y[1] = -fq[1]; y[2] = -fw;
305 		}
306 	}
307 	return n&7;
308 }
309 
310 
311 
312 
313 
314 /* Quad-precision floating point argument reduction.
315    Copyright (C) 1999 Free Software Foundation, Inc.
316    This file is part of the GNU C Library.
317    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
318 
319    The GNU C Library is free software; you can redistribute it and/or
320    modify it under the terms of the GNU Lesser General Public
321    License as published by the Free Software Foundation; either
322    version 2.1 of the License, or (at your option) any later version.
323 
324    The GNU C Library is distributed in the hope that it will be useful,
325    but WITHOUT ANY WARRANTY; without even the implied warranty of
326    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
327    Lesser General Public License for more details.
328 
329    You should have received a copy of the GNU Lesser General Public
330    License along with the GNU C Library; if not, write to the Free
331    Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
332    02111-1307 USA.  */
333 
334 /*
335  * Table of constants for 2/pi, 5628 hexadecimal digits of 2/pi
336  */
337 static const int32_t two_over_pi[] = {
338 0xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
339 0x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
340 0x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
341 0xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
342 0x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
343 0x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
344 0x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
345 0xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
346 0x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
347 0x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
348 0x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
349 0x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
350 0xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
351 0xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
352 0xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
353 0x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
354 0x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
355 0xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
356 0xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
357 0xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
358 0xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
359 0x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
360 0xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
361 0x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
362 0x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
363 0xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
364 0xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
365 0xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
366 0x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
367 0xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
368 0x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
369 0xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
370 0x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
371 0x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
372 0x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
373 0xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
374 0x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
375 0x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
376 0xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
377 0x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
378 0x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
379 0x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
380 0x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
381 0x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
382 0x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
383 0xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
384 0x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
385 0xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
386 0xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
387 0xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
388 0xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
389 0x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
390 0x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
391 0x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
392 0xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
393 0x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
394 0x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
395 0x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
396 0xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
397 0x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
398 0xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
399 0xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
400 0x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
401 0x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
402 0x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
403 0xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
404 0x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
405 0x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
406 0xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
407 0x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
408 0xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
409 0xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
410 0x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
411 0xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
412 0x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
413 0xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
414 0x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
415 0x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
416 0x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
417 0xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
418 0x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
419 0xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
420 0x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
421 0xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
422 0x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
423 0x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
424 0xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
425 0x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
426 0xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
427 0x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
428 0x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
429 0x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
430 0x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
431 0xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
432 0xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
433 0x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
434 0xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
435 0x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
436 0xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
437 0xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
438 0x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
439 0xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
440 0x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
441 0x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
442 0x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
443 0xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
444 0xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
445 0x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
446 0x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
447 0xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
448 0x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
449 0x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
450 0x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
451 0x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
452 0x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
453 0xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
454 0x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
455 0x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
456 0xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
457 0x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
458 0xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
459 0xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
460 0xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
461 0x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
462 0xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
463 0xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
464 0xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
465 0x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
466 0xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
467 0x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
468 0xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
469 0xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
470 0x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
471 0x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
472 0x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
473 0x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
474 0x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
475 0x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
476 0xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
477 0x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
478 0x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
479 0x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
480 0xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
481 0xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
482 0xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
483 0x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
484 0xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
485 0x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
486 0xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
487 0x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
488 0xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
489 0xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
490 0xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
491 0x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
492 0xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
493 0xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
494 0x7b7b89, 0x483d38,
495 };
496 
497 static const __float128 c[] = {
498 /* 93 bits of pi/2 */
499 #define PI_2_1 c[0]
500  1.57079632679489661923132169155131424e+00Q, /* 3fff921fb54442d18469898cc5100000 */
501 
502 /* pi/2 - PI_2_1 */
503 #define PI_2_1t c[1]
504  8.84372056613570112025531863263659260e-29Q, /* 3fa1c06e0e68948127044533e63a0106 */
505 };
506 
507 
508 int32_t
__quadmath_rem_pio2q(__float128 x,__float128 * y)509 __quadmath_rem_pio2q (__float128 x, __float128 *y)
510 {
511   __float128 z, w, t;
512   double tx[8];
513   int64_t exp, n, ix, hx;
514   uint64_t lx;
515 
516   GET_FLT128_WORDS64 (hx, lx, x);
517   ix = hx & 0x7fffffffffffffffLL;
518   if (ix <= 0x3ffe921fb54442d1LL)	/* x in <-pi/4, pi/4> */
519     {
520       y[0] = x;
521       y[1] = 0;
522       return 0;
523     }
524 
525   if (ix < 0x40002d97c7f3321dLL)	/* |x| in <pi/4, 3pi/4) */
526     {
527       if (hx > 0)
528 	{
529 	  /* 113 + 93 bit PI is ok */
530 	  z = x - PI_2_1;
531 	  y[0] = z - PI_2_1t;
532 	  y[1] = (z - y[0]) - PI_2_1t;
533 	  return 1;
534 	}
535       else
536         {
537 	  /* 113 + 93 bit PI is ok */
538 	  z = x + PI_2_1;
539 	  y[0] = z + PI_2_1t;
540 	  y[1] = (z - y[0]) + PI_2_1t;
541 	  return -1;
542 	}
543     }
544 
545   if (ix >= 0x7fff000000000000LL)	/* x is +=oo or NaN */
546     {
547       y[0] = x - x;
548       y[1] = y[0];
549       return 0;
550     }
551 
552   /* Handle large arguments.
553      We split the 113 bits of the mantissa into 5 24bit integers
554      stored in a double array.  */
555   exp = (ix >> 48) - 16383 - 23;
556 
557   /* This is faster than doing this in floating point, because we
558      have to convert it to integers anyway and like this we can keep
559      both integer and floating point units busy.  */
560   tx [0] = (double)(((ix >> 25) & 0x7fffff) | 0x800000);
561   tx [1] = (double)((ix >> 1) & 0xffffff);
562   tx [2] = (double)(((ix << 23) | (lx >> 41)) & 0xffffff);
563   tx [3] = (double)((lx >> 17) & 0xffffff);
564   tx [4] = (double)((lx << 7) & 0xffffff);
565 
566   n = __quadmath_kernel_rem_pio2 (tx, tx + 5, exp,
567 				  ((lx << 7) & 0xffffff) ? 5 : 4,
568 				  3, two_over_pi);
569 
570   /* The result is now stored in 3 double values, we need to convert it into
571      two __float128 values.  */
572   t = (__float128) tx [6] + (__float128) tx [7];
573   w = (__float128) tx [5];
574 
575   if (hx >= 0)
576     {
577       y[0] = w + t;
578       y[1] = t - (y[0] - w);
579       return n;
580     }
581   else
582     {
583       y[0] = -(w + t);
584       y[1] = -t - (y[0] + w);
585       return -n;
586     }
587 }
588