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
133static const int init_jk[] = {2,3,4,6}; /* initial value for jk */
134
135static 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
146static 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
153static int
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;
178recompute:
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-2017 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 */
337static const int32_t two_over_pi[] = {
3380xa2f983, 0x6e4e44, 0x1529fc, 0x2757d1, 0xf534dd, 0xc0db62,
3390x95993c, 0x439041, 0xfe5163, 0xabdebb, 0xc561b7, 0x246e3a,
3400x424dd2, 0xe00649, 0x2eea09, 0xd1921c, 0xfe1deb, 0x1cb129,
3410xa73ee8, 0x8235f5, 0x2ebb44, 0x84e99c, 0x7026b4, 0x5f7e41,
3420x3991d6, 0x398353, 0x39f49c, 0x845f8b, 0xbdf928, 0x3b1ff8,
3430x97ffde, 0x05980f, 0xef2f11, 0x8b5a0a, 0x6d1f6d, 0x367ecf,
3440x27cb09, 0xb74f46, 0x3f669e, 0x5fea2d, 0x7527ba, 0xc7ebe5,
3450xf17b3d, 0x0739f7, 0x8a5292, 0xea6bfb, 0x5fb11f, 0x8d5d08,
3460x560330, 0x46fc7b, 0x6babf0, 0xcfbc20, 0x9af436, 0x1da9e3,
3470x91615e, 0xe61b08, 0x659985, 0x5f14a0, 0x68408d, 0xffd880,
3480x4d7327, 0x310606, 0x1556ca, 0x73a8c9, 0x60e27b, 0xc08c6b,
3490x47c419, 0xc367cd, 0xdce809, 0x2a8359, 0xc4768b, 0x961ca6,
3500xddaf44, 0xd15719, 0x053ea5, 0xff0705, 0x3f7e33, 0xe832c2,
3510xde4f98, 0x327dbb, 0xc33d26, 0xef6b1e, 0x5ef89f, 0x3a1f35,
3520xcaf27f, 0x1d87f1, 0x21907c, 0x7c246a, 0xfa6ed5, 0x772d30,
3530x433b15, 0xc614b5, 0x9d19c3, 0xc2c4ad, 0x414d2c, 0x5d000c,
3540x467d86, 0x2d71e3, 0x9ac69b, 0x006233, 0x7cd2b4, 0x97a7b4,
3550xd55537, 0xf63ed7, 0x1810a3, 0xfc764d, 0x2a9d64, 0xabd770,
3560xf87c63, 0x57b07a, 0xe71517, 0x5649c0, 0xd9d63b, 0x3884a7,
3570xcb2324, 0x778ad6, 0x23545a, 0xb91f00, 0x1b0af1, 0xdfce19,
3580xff319f, 0x6a1e66, 0x615799, 0x47fbac, 0xd87f7e, 0xb76522,
3590x89e832, 0x60bfe6, 0xcdc4ef, 0x09366c, 0xd43f5d, 0xd7de16,
3600xde3b58, 0x929bde, 0x2822d2, 0xe88628, 0x4d58e2, 0x32cac6,
3610x16e308, 0xcb7de0, 0x50c017, 0xa71df3, 0x5be018, 0x34132e,
3620x621283, 0x014883, 0x5b8ef5, 0x7fb0ad, 0xf2e91e, 0x434a48,
3630xd36710, 0xd8ddaa, 0x425fae, 0xce616a, 0xa4280a, 0xb499d3,
3640xf2a606, 0x7f775c, 0x83c2a3, 0x883c61, 0x78738a, 0x5a8caf,
3650xbdd76f, 0x63a62d, 0xcbbff4, 0xef818d, 0x67c126, 0x45ca55,
3660x36d9ca, 0xd2a828, 0x8d61c2, 0x77c912, 0x142604, 0x9b4612,
3670xc459c4, 0x44c5c8, 0x91b24d, 0xf31700, 0xad43d4, 0xe54929,
3680x10d5fd, 0xfcbe00, 0xcc941e, 0xeece70, 0xf53e13, 0x80f1ec,
3690xc3e7b3, 0x28f8c7, 0x940593, 0x3e71c1, 0xb3092e, 0xf3450b,
3700x9c1288, 0x7b20ab, 0x9fb52e, 0xc29247, 0x2f327b, 0x6d550c,
3710x90a772, 0x1fe76b, 0x96cb31, 0x4a1679, 0xe27941, 0x89dff4,
3720x9794e8, 0x84e6e2, 0x973199, 0x6bed88, 0x365f5f, 0x0efdbb,
3730xb49a48, 0x6ca467, 0x427271, 0x325d8d, 0xb8159f, 0x09e5bc,
3740x25318d, 0x3974f7, 0x1c0530, 0x010c0d, 0x68084b, 0x58ee2c,
3750x90aa47, 0x02e774, 0x24d6bd, 0xa67df7, 0x72486e, 0xef169f,
3760xa6948e, 0xf691b4, 0x5153d1, 0xf20acf, 0x339820, 0x7e4bf5,
3770x6863b2, 0x5f3edd, 0x035d40, 0x7f8985, 0x295255, 0xc06437,
3780x10d86d, 0x324832, 0x754c5b, 0xd4714e, 0x6e5445, 0xc1090b,
3790x69f52a, 0xd56614, 0x9d0727, 0x50045d, 0xdb3bb4, 0xc576ea,
3800x17f987, 0x7d6b49, 0xba271d, 0x296996, 0xacccc6, 0x5414ad,
3810x6ae290, 0x89d988, 0x50722c, 0xbea404, 0x940777, 0x7030f3,
3820x27fc00, 0xa871ea, 0x49c266, 0x3de064, 0x83dd97, 0x973fa3,
3830xfd9443, 0x8c860d, 0xde4131, 0x9d3992, 0x8c70dd, 0xe7b717,
3840x3bdf08, 0x2b3715, 0xa0805c, 0x93805a, 0x921110, 0xd8e80f,
3850xaf806c, 0x4bffdb, 0x0f9038, 0x761859, 0x15a562, 0xbbcb61,
3860xb989c7, 0xbd4010, 0x04f2d2, 0x277549, 0xf6b6eb, 0xbb22db,
3870xaa140a, 0x2f2689, 0x768364, 0x333b09, 0x1a940e, 0xaa3a51,
3880xc2a31d, 0xaeedaf, 0x12265c, 0x4dc26d, 0x9c7a2d, 0x9756c0,
3890x833f03, 0xf6f009, 0x8c402b, 0x99316d, 0x07b439, 0x15200c,
3900x5bc3d8, 0xc492f5, 0x4badc6, 0xa5ca4e, 0xcd37a7, 0x36a9e6,
3910x9492ab, 0x6842dd, 0xde6319, 0xef8c76, 0x528b68, 0x37dbfc,
3920xaba1ae, 0x3115df, 0xa1ae00, 0xdafb0c, 0x664d64, 0xb705ed,
3930x306529, 0xbf5657, 0x3aff47, 0xb9f96a, 0xf3be75, 0xdf9328,
3940x3080ab, 0xf68c66, 0x15cb04, 0x0622fa, 0x1de4d9, 0xa4b33d,
3950x8f1b57, 0x09cd36, 0xe9424e, 0xa4be13, 0xb52333, 0x1aaaf0,
3960xa8654f, 0xa5c1d2, 0x0f3f0b, 0xcd785b, 0x76f923, 0x048b7b,
3970x721789, 0x53a6c6, 0xe26e6f, 0x00ebef, 0x584a9b, 0xb7dac4,
3980xba66aa, 0xcfcf76, 0x1d02d1, 0x2df1b1, 0xc1998c, 0x77adc3,
3990xda4886, 0xa05df7, 0xf480c6, 0x2ff0ac, 0x9aecdd, 0xbc5c3f,
4000x6dded0, 0x1fc790, 0xb6db2a, 0x3a25a3, 0x9aaf00, 0x9353ad,
4010x0457b6, 0xb42d29, 0x7e804b, 0xa707da, 0x0eaa76, 0xa1597b,
4020x2a1216, 0x2db7dc, 0xfde5fa, 0xfedb89, 0xfdbe89, 0x6c76e4,
4030xfca906, 0x70803e, 0x156e85, 0xff87fd, 0x073e28, 0x336761,
4040x86182a, 0xeabd4d, 0xafe7b3, 0x6e6d8f, 0x396795, 0x5bbf31,
4050x48d784, 0x16df30, 0x432dc7, 0x356125, 0xce70c9, 0xb8cb30,
4060xfd6cbf, 0xa200a4, 0xe46c05, 0xa0dd5a, 0x476f21, 0xd21262,
4070x845cb9, 0x496170, 0xe0566b, 0x015299, 0x375550, 0xb7d51e,
4080xc4f133, 0x5f6e13, 0xe4305d, 0xa92e85, 0xc3b21d, 0x3632a1,
4090xa4b708, 0xd4b1ea, 0x21f716, 0xe4698f, 0x77ff27, 0x80030c,
4100x2d408d, 0xa0cd4f, 0x99a520, 0xd3a2b3, 0x0a5d2f, 0x42f9b4,
4110xcbda11, 0xd0be7d, 0xc1db9b, 0xbd17ab, 0x81a2ca, 0x5c6a08,
4120x17552e, 0x550027, 0xf0147f, 0x8607e1, 0x640b14, 0x8d4196,
4130xdebe87, 0x2afdda, 0xb6256b, 0x34897b, 0xfef305, 0x9ebfb9,
4140x4f6a68, 0xa82a4a, 0x5ac44f, 0xbcf82d, 0x985ad7, 0x95c7f4,
4150x8d4d0d, 0xa63a20, 0x5f57a4, 0xb13f14, 0x953880, 0x0120cc,
4160x86dd71, 0xb6dec9, 0xf560bf, 0x11654d, 0x6b0701, 0xacb08c,
4170xd0c0b2, 0x485551, 0x0efb1e, 0xc37295, 0x3b06a3, 0x3540c0,
4180x7bdc06, 0xcc45e0, 0xfa294e, 0xc8cad6, 0x41f3e8, 0xde647c,
4190xd8649b, 0x31bed9, 0xc397a4, 0xd45877, 0xc5e369, 0x13daf0,
4200x3c3aba, 0x461846, 0x5f7555, 0xf5bdd2, 0xc6926e, 0x5d2eac,
4210xed440e, 0x423e1c, 0x87c461, 0xe9fd29, 0xf3d6e7, 0xca7c22,
4220x35916f, 0xc5e008, 0x8dd7ff, 0xe26a6e, 0xc6fdb0, 0xc10893,
4230x745d7c, 0xb2ad6b, 0x9d6ecd, 0x7b723e, 0x6a11c6, 0xa9cff7,
4240xdf7329, 0xbac9b5, 0x5100b7, 0x0db2e2, 0x24ba74, 0x607de5,
4250x8ad874, 0x2c150d, 0x0c1881, 0x94667e, 0x162901, 0x767a9f,
4260xbefdfd, 0xef4556, 0x367ed9, 0x13d9ec, 0xb9ba8b, 0xfc97c4,
4270x27a831, 0xc36ef1, 0x36c594, 0x56a8d8, 0xb5a8b4, 0x0ecccf,
4280x2d8912, 0x34576f, 0x89562c, 0xe3ce99, 0xb920d6, 0xaa5e6b,
4290x9c2a3e, 0xcc5f11, 0x4a0bfd, 0xfbf4e1, 0x6d3b8e, 0x2c86e2,
4300x84d4e9, 0xa9b4fc, 0xd1eeef, 0xc9352e, 0x61392f, 0x442138,
4310xc8d91b, 0x0afc81, 0x6a4afb, 0xd81c2f, 0x84b453, 0x8c994e,
4320xcc2254, 0xdc552a, 0xd6c6c0, 0x96190b, 0xb8701a, 0x649569,
4330x605a26, 0xee523f, 0x0f117f, 0x11b5f4, 0xf5cbfc, 0x2dbc34,
4340xeebc34, 0xcc5de8, 0x605edd, 0x9b8e67, 0xef3392, 0xb817c9,
4350x9b5861, 0xbc57e1, 0xc68351, 0x103ed8, 0x4871dd, 0xdd1c2d,
4360xa118af, 0x462c21, 0xd7f359, 0x987ad9, 0xc0549e, 0xfa864f,
4370xfc0656, 0xae79e5, 0x362289, 0x22ad38, 0xdc9367, 0xaae855,
4380x382682, 0x9be7ca, 0xa40d51, 0xb13399, 0x0ed7a9, 0x480569,
4390xf0b265, 0xa7887f, 0x974c88, 0x36d1f9, 0xb39221, 0x4a827b,
4400x21cf98, 0xdc9f40, 0x5547dc, 0x3a74e1, 0x42eb67, 0xdf9dfe,
4410x5fd45e, 0xa4677b, 0x7aacba, 0xa2f655, 0x23882b, 0x55ba41,
4420x086e59, 0x862a21, 0x834739, 0xe6e389, 0xd49ee5, 0x40fb49,
4430xe956ff, 0xca0f1c, 0x8a59c5, 0x2bfa94, 0xc5c1d3, 0xcfc50f,
4440xae5adb, 0x86c547, 0x624385, 0x3b8621, 0x94792c, 0x876110,
4450x7b4c2a, 0x1a2c80, 0x12bf43, 0x902688, 0x893c78, 0xe4c4a8,
4460x7bdbe5, 0xc23ac4, 0xeaf426, 0x8a67f7, 0xbf920d, 0x2ba365,
4470xb1933d, 0x0b7cbd, 0xdc51a4, 0x63dd27, 0xdde169, 0x19949a,
4480x9529a8, 0x28ce68, 0xb4ed09, 0x209f44, 0xca984e, 0x638270,
4490x237c7e, 0x32b90f, 0x8ef5a7, 0xe75614, 0x08f121, 0x2a9db5,
4500x4d7e6f, 0x5119a5, 0xabf9b5, 0xd6df82, 0x61dd96, 0x023616,
4510x9f3ac4, 0xa1a283, 0x6ded72, 0x7a8d39, 0xa9b882, 0x5c326b,
4520x5b2746, 0xed3400, 0x7700d2, 0x55f4fc, 0x4d5901, 0x8071e0,
4530xe13f89, 0xb295f3, 0x64a8f1, 0xaea74b, 0x38fc4c, 0xeab2bb,
4540x47270b, 0xabc3a7, 0x34ba60, 0x52dd34, 0xf8563a, 0xeb7e8a,
4550x31bb36, 0x5895b7, 0x47f7a9, 0x94c3aa, 0xd39225, 0x1e7f3e,
4560xd8974e, 0xbba94f, 0xd8ae01, 0xe661b4, 0x393d8e, 0xa523aa,
4570x33068e, 0x1633b5, 0x3bb188, 0x1d3a9d, 0x4013d0, 0xcc1be5,
4580xf862e7, 0x3bf28f, 0x39b5bf, 0x0bc235, 0x22747e, 0xa247c0,
4590xd52d1f, 0x19add3, 0x9094df, 0x9311d0, 0xb42b25, 0x496db2,
4600xe264b2, 0x5ef135, 0x3bc6a4, 0x1a4ad0, 0xaac92e, 0x64e886,
4610x573091, 0x982cfb, 0x311b1a, 0x08728b, 0xbdcee1, 0x60e142,
4620xeb641d, 0xd0bba3, 0xe559d4, 0x597b8c, 0x2a4483, 0xf332ba,
4630xf84867, 0x2c8d1b, 0x2fa9b0, 0x50f3dd, 0xf9f573, 0xdb61b4,
4640xfe233e, 0x6c41a6, 0xeea318, 0x775a26, 0xbc5e5c, 0xcea708,
4650x94dc57, 0xe20196, 0xf1e839, 0xbe4851, 0x5d2d2f, 0x4e9555,
4660xd96ec2, 0xe7d755, 0x6304e0, 0xc02e0e, 0xfc40a0, 0xbbf9b3,
4670x7125a7, 0x222dfb, 0xf619d8, 0x838c1c, 0x6619e6, 0xb20d55,
4680xbb5137, 0x79e809, 0xaf9149, 0x0d73de, 0x0b0da5, 0xce7f58,
4690xac1934, 0x724667, 0x7a1a13, 0x9e26bc, 0x4555e7, 0x585cb5,
4700x711d14, 0x486991, 0x480d60, 0x56adab, 0xd62f64, 0x96ee0c,
4710x212ff3, 0x5d6d88, 0xa67684, 0x95651e, 0xab9e0a, 0x4ddefe,
4720x571010, 0x836a39, 0xf8ea31, 0x9e381d, 0xeac8b1, 0xcac96b,
4730x37f21e, 0xd505e9, 0x984743, 0x9fc56c, 0x0331b7, 0x3b8bf8,
4740x86e56a, 0x8dc343, 0x6230e7, 0x93cfd5, 0x6a8f2d, 0x733005,
4750x1af021, 0xa09fcb, 0x7415a1, 0xd56b23, 0x6ff725, 0x2f4bc7,
4760xb8a591, 0x7fac59, 0x5c55de, 0x212c38, 0xb13296, 0x5cff50,
4770x366262, 0xfa7b16, 0xf4d9a6, 0x2acfe7, 0xf07403, 0xd4d604,
4780x6fd916, 0x31b1bf, 0xcbb450, 0x5bd7c8, 0x0ce194, 0x6bd643,
4790x4fd91c, 0xdf4543, 0x5f3453, 0xe2b5aa, 0xc9aec8, 0x131485,
4800xf9d2bf, 0xbadb9e, 0x76f5b9, 0xaf15cf, 0xca3182, 0x14b56d,
4810xe9fe4d, 0x50fc35, 0xf5aed5, 0xa2d0c1, 0xc96057, 0x192eb6,
4820xe91d92, 0x07d144, 0xaea3c6, 0x343566, 0x26d5b4, 0x3161e2,
4830x37f1a2, 0x209eff, 0x958e23, 0x493798, 0x35f4a6, 0x4bdc02,
4840xc2be13, 0xbe80a0, 0x0b72a3, 0x115c5f, 0x1e1bd1, 0x0db4d3,
4850x869e85, 0x96976b, 0x2ac91f, 0x8a26c2, 0x3070f0, 0x041412,
4860xfc9fa5, 0xf72a38, 0x9c6878, 0xe2aa76, 0x50cfe1, 0x559274,
4870x934e38, 0x0a92f7, 0x5533f0, 0xa63db4, 0x399971, 0xe2b755,
4880xa98a7c, 0x008f19, 0xac54d2, 0x2ea0b4, 0xf5f3e0, 0x60c849,
4890xffd269, 0xae52ce, 0x7a5fdd, 0xe9ce06, 0xfb0ae8, 0xa50cce,
4900xea9d3e, 0x3766dd, 0xb834f5, 0x0da090, 0x846f88, 0x4ae3d5,
4910x099a03, 0x2eae2d, 0xfcb40a, 0xfb9b33, 0xe281dd, 0x1b16ba,
4920xd8c0af, 0xd96b97, 0xb52dc9, 0x9c277f, 0x5951d5, 0x21ccd6,
4930xb6496b, 0x584562, 0xb3baf2, 0xa1a5c4, 0x7ca2cf, 0xa9b93d,
4940x7b7b89, 0x483d38,
495};
496
497static const __float128 c[] = {
498/* 113 bits of pi/2 */
499#define PI_2_1 c[0]
500 0x1.921fb54442d18469898cc51701b8p+0Q,
501
502/* pi/2 - PI_2_1 */
503#define PI_2_1t c[1]
504 0x3.9a252049c1114cf98e804177d4c8p-116Q,
505};
506
507
508int32_t
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 + 113 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 + 113 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