1 /* @(#)e_fmod.c 5.1 93/09/24 */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 13 /* 14 * fmodl(x,y) 15 * Return x mod y in exact arithmetic 16 * Method: shift and subtract 17 */ 18 19 #include <math.h> 20 21 #include "math_private.h" 22 23 static const long double one = 1.0, Zero[] = {0.0, -0.0,}; 24 25 long double 26 fmodl(long double x, long double y) 27 { 28 int64_t n,hx,hy,hz,ix,iy,sx,i; 29 u_int64_t lx,ly,lz; 30 31 GET_LDOUBLE_WORDS64(hx,lx,x); 32 GET_LDOUBLE_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[(u_int64_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[(u_int64_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[(u_int64_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_LDOUBLE_WORDS64(x,hx|sx,lx); 115 } else { /* subnormal output */ 116 n = -16382 - iy; 117 if(n<=48) { 118 lx = (lx>>n)|((u_int64_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_LDOUBLE_WORDS64(x,hx|sx,lx); 126 x *= one; /* create necessary signal */ 127 } 128 return x; /* exact output */ 129 } 130