xref: /netbsd-src/external/gpl3/gcc/dist/libquadmath/math/fmodq.c (revision 181254a7b1bdde6873432bffef2d2decc4b5c22f)
1*181254a7Smrg /* e_fmodl.c -- long double version of e_fmod.c.
2*181254a7Smrg  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
3*181254a7Smrg  */
4*181254a7Smrg /*
5*181254a7Smrg  * ====================================================
6*181254a7Smrg  * Copyright (C) 1993, 2011 by Sun Microsystems, Inc. All rights reserved.
7*181254a7Smrg  *
8*181254a7Smrg  * Developed at SunPro, a Sun Microsystems, Inc. business.
9*181254a7Smrg  * Permission to use, copy, modify, and distribute this
10*181254a7Smrg  * software is freely granted, provided that this notice
11*181254a7Smrg  * is preserved.
12*181254a7Smrg  * ====================================================
13*181254a7Smrg  */
14*181254a7Smrg 
15*181254a7Smrg /*
16*181254a7Smrg  * fmodq(x,y)
17*181254a7Smrg  * Return x mod y in exact arithmetic
18*181254a7Smrg  * Method: shift and subtract
19*181254a7Smrg  */
20*181254a7Smrg 
21*181254a7Smrg #include "quadmath-imp.h"
22*181254a7Smrg 
23*181254a7Smrg static const __float128 one = 1.0, Zero[] = {0.0, -0.0,};
24*181254a7Smrg 
25*181254a7Smrg __float128
fmodq(__float128 x,__float128 y)26*181254a7Smrg fmodq (__float128 x, __float128 y)
27*181254a7Smrg {
28*181254a7Smrg 	int64_t n,hx,hy,hz,ix,iy,sx,i;
29*181254a7Smrg 	uint64_t lx,ly,lz;
30*181254a7Smrg 
31*181254a7Smrg 	GET_FLT128_WORDS64(hx,lx,x);
32*181254a7Smrg 	GET_FLT128_WORDS64(hy,ly,y);
33*181254a7Smrg 	sx = hx&0x8000000000000000ULL;		/* sign of x */
34*181254a7Smrg 	hx ^=sx;				/* |x| */
35*181254a7Smrg 	hy &= 0x7fffffffffffffffLL;		/* |y| */
36*181254a7Smrg 
37*181254a7Smrg     /* purge off exception values */
38*181254a7Smrg 	if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
39*181254a7Smrg 	  ((hy|((ly|-ly)>>63))>0x7fff000000000000LL))	/* or y is NaN */
40*181254a7Smrg 	    return (x*y)/(x*y);
41*181254a7Smrg 	if(hx<=hy) {
42*181254a7Smrg 	    if((hx<hy)||(lx<ly)) return x;	/* |x|<|y| return x */
43*181254a7Smrg 	    if(lx==ly)
44*181254a7Smrg 		return Zero[(uint64_t)sx>>63];	/* |x|=|y| return x*0*/
45*181254a7Smrg 	}
46*181254a7Smrg 
47*181254a7Smrg     /* determine ix = ilogb(x) */
48*181254a7Smrg 	if(hx<0x0001000000000000LL) {	/* subnormal x */
49*181254a7Smrg 	    if(hx==0) {
50*181254a7Smrg 		for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
51*181254a7Smrg 	    } else {
52*181254a7Smrg 		for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
53*181254a7Smrg 	    }
54*181254a7Smrg 	} else ix = (hx>>48)-0x3fff;
55*181254a7Smrg 
56*181254a7Smrg     /* determine iy = ilogb(y) */
57*181254a7Smrg 	if(hy<0x0001000000000000LL) {	/* subnormal y */
58*181254a7Smrg 	    if(hy==0) {
59*181254a7Smrg 		for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
60*181254a7Smrg 	    } else {
61*181254a7Smrg 		for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
62*181254a7Smrg 	    }
63*181254a7Smrg 	} else iy = (hy>>48)-0x3fff;
64*181254a7Smrg 
65*181254a7Smrg     /* set up {hx,lx}, {hy,ly} and align y to x */
66*181254a7Smrg 	if(ix >= -16382)
67*181254a7Smrg 	    hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
68*181254a7Smrg 	else {		/* subnormal x, shift x to normal */
69*181254a7Smrg 	    n = -16382-ix;
70*181254a7Smrg 	    if(n<=63) {
71*181254a7Smrg 		hx = (hx<<n)|(lx>>(64-n));
72*181254a7Smrg 		lx <<= n;
73*181254a7Smrg 	    } else {
74*181254a7Smrg 		hx = lx<<(n-64);
75*181254a7Smrg 		lx = 0;
76*181254a7Smrg 	    }
77*181254a7Smrg 	}
78*181254a7Smrg 	if(iy >= -16382)
79*181254a7Smrg 	    hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
80*181254a7Smrg 	else {		/* subnormal y, shift y to normal */
81*181254a7Smrg 	    n = -16382-iy;
82*181254a7Smrg 	    if(n<=63) {
83*181254a7Smrg 		hy = (hy<<n)|(ly>>(64-n));
84*181254a7Smrg 		ly <<= n;
85*181254a7Smrg 	    } else {
86*181254a7Smrg 		hy = ly<<(n-64);
87*181254a7Smrg 		ly = 0;
88*181254a7Smrg 	    }
89*181254a7Smrg 	}
90*181254a7Smrg 
91*181254a7Smrg     /* fix point fmod */
92*181254a7Smrg 	n = ix - iy;
93*181254a7Smrg 	while(n--) {
94*181254a7Smrg 	    hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
95*181254a7Smrg 	    if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
96*181254a7Smrg 	    else {
97*181254a7Smrg 		if((hz|lz)==0)		/* return sign(x)*0 */
98*181254a7Smrg 		    return Zero[(uint64_t)sx>>63];
99*181254a7Smrg 		hx = hz+hz+(lz>>63); lx = lz+lz;
100*181254a7Smrg 	    }
101*181254a7Smrg 	}
102*181254a7Smrg 	hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
103*181254a7Smrg 	if(hz>=0) {hx=hz;lx=lz;}
104*181254a7Smrg 
105*181254a7Smrg     /* convert back to floating value and restore the sign */
106*181254a7Smrg 	if((hx|lx)==0)			/* return sign(x)*0 */
107*181254a7Smrg 	    return Zero[(uint64_t)sx>>63];
108*181254a7Smrg 	while(hx<0x0001000000000000LL) {	/* normalize x */
109*181254a7Smrg 	    hx = hx+hx+(lx>>63); lx = lx+lx;
110*181254a7Smrg 	    iy -= 1;
111*181254a7Smrg 	}
112*181254a7Smrg 	if(iy>= -16382) {	/* normalize output */
113*181254a7Smrg 	    hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
114*181254a7Smrg 	    SET_FLT128_WORDS64(x,hx|sx,lx);
115*181254a7Smrg 	} else {		/* subnormal output */
116*181254a7Smrg 	    n = -16382 - iy;
117*181254a7Smrg 	    if(n<=48) {
118*181254a7Smrg 		lx = (lx>>n)|((uint64_t)hx<<(64-n));
119*181254a7Smrg 		hx >>= n;
120*181254a7Smrg 	    } else if (n<=63) {
121*181254a7Smrg 		lx = (hx<<(64-n))|(lx>>n); hx = sx;
122*181254a7Smrg 	    } else {
123*181254a7Smrg 		lx = hx>>(n-64); hx = sx;
124*181254a7Smrg 	    }
125*181254a7Smrg 	    SET_FLT128_WORDS64(x,hx|sx,lx);
126*181254a7Smrg 	    x *= one;		/* create necessary signal */
127*181254a7Smrg 	}
128*181254a7Smrg 	return x;		/* exact output */
129*181254a7Smrg }
130