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