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