1141296Sdas 2141296Sdas/* @(#)e_fmod.c 1.3 95/01/18 */ 32116Sjkh/* 42116Sjkh * ==================================================== 52116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 62116Sjkh * 7141296Sdas * Developed at SunSoft, a Sun Microsystems, Inc. business. 82116Sjkh * Permission to use, copy, modify, and distribute this 9141296Sdas * software is freely granted, provided that this notice 102116Sjkh * is preserved. 112116Sjkh * ==================================================== 122116Sjkh */ 132116Sjkh 14176451Sdas#include <sys/cdefs.h> 15176451Sdas__FBSDID("$FreeBSD$"); 162116Sjkh 17141296Sdas/* 182116Sjkh * __ieee754_fmod(x,y) 192116Sjkh * Return x mod y in exact arithmetic 202116Sjkh * Method: shift and subtract 212116Sjkh */ 222116Sjkh 232116Sjkh#include "math.h" 242116Sjkh#include "math_private.h" 252116Sjkh 262116Sjkhstatic const double one = 1.0, Zero[] = {0.0, -0.0,}; 272116Sjkh 2897407Salfreddouble 29117912Speter__ieee754_fmod(double x, double y) 302116Sjkh{ 312116Sjkh int32_t n,hx,hy,hz,ix,iy,sx,i; 322116Sjkh u_int32_t lx,ly,lz; 332116Sjkh 342116Sjkh EXTRACT_WORDS(hx,lx,x); 352116Sjkh EXTRACT_WORDS(hy,ly,y); 362116Sjkh sx = hx&0x80000000; /* sign of x */ 372116Sjkh hx ^=sx; /* |x| */ 382116Sjkh hy &= 0x7fffffff; /* |y| */ 392116Sjkh 402116Sjkh /* purge off exception values */ 412116Sjkh if((hy|ly)==0||(hx>=0x7ff00000)|| /* y=0,or x not finite */ 422116Sjkh ((hy|((ly|-ly)>>31))>0x7ff00000)) /* or y is NaN */ 432116Sjkh return (x*y)/(x*y); 442116Sjkh if(hx<=hy) { 452116Sjkh if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */ 46141296Sdas if(lx==ly) 472116Sjkh return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 482116Sjkh } 492116Sjkh 502116Sjkh /* determine ix = ilogb(x) */ 512116Sjkh if(hx<0x00100000) { /* subnormal x */ 522116Sjkh if(hx==0) { 532116Sjkh for (ix = -1043, i=lx; i>0; i<<=1) ix -=1; 542116Sjkh } else { 552116Sjkh for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1; 562116Sjkh } 572116Sjkh } else ix = (hx>>20)-1023; 582116Sjkh 592116Sjkh /* determine iy = ilogb(y) */ 602116Sjkh if(hy<0x00100000) { /* subnormal y */ 612116Sjkh if(hy==0) { 622116Sjkh for (iy = -1043, i=ly; i>0; i<<=1) iy -=1; 632116Sjkh } else { 642116Sjkh for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1; 652116Sjkh } 662116Sjkh } else iy = (hy>>20)-1023; 672116Sjkh 682116Sjkh /* set up {hx,lx}, {hy,ly} and align y to x */ 69141296Sdas if(ix >= -1022) 702116Sjkh hx = 0x00100000|(0x000fffff&hx); 712116Sjkh else { /* subnormal x, shift x to normal */ 722116Sjkh n = -1022-ix; 732116Sjkh if(n<=31) { 742116Sjkh hx = (hx<<n)|(lx>>(32-n)); 752116Sjkh lx <<= n; 762116Sjkh } else { 772116Sjkh hx = lx<<(n-32); 782116Sjkh lx = 0; 792116Sjkh } 802116Sjkh } 81141296Sdas if(iy >= -1022) 822116Sjkh hy = 0x00100000|(0x000fffff&hy); 832116Sjkh else { /* subnormal y, shift y to normal */ 842116Sjkh n = -1022-iy; 852116Sjkh if(n<=31) { 862116Sjkh hy = (hy<<n)|(ly>>(32-n)); 872116Sjkh ly <<= n; 882116Sjkh } else { 892116Sjkh hy = ly<<(n-32); 902116Sjkh ly = 0; 912116Sjkh } 922116Sjkh } 932116Sjkh 942116Sjkh /* fix point fmod */ 952116Sjkh n = ix - iy; 962116Sjkh while(n--) { 972116Sjkh hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 982116Sjkh if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;} 992116Sjkh else { 1002116Sjkh if((hz|lz)==0) /* return sign(x)*0 */ 1012116Sjkh return Zero[(u_int32_t)sx>>31]; 1022116Sjkh hx = hz+hz+(lz>>31); lx = lz+lz; 1032116Sjkh } 1042116Sjkh } 1052116Sjkh hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1; 1062116Sjkh if(hz>=0) {hx=hz;lx=lz;} 1072116Sjkh 1082116Sjkh /* convert back to floating value and restore the sign */ 1092116Sjkh if((hx|lx)==0) /* return sign(x)*0 */ 1108870Srgrimes return Zero[(u_int32_t)sx>>31]; 1112116Sjkh while(hx<0x00100000) { /* normalize x */ 1122116Sjkh hx = hx+hx+(lx>>31); lx = lx+lx; 1132116Sjkh iy -= 1; 1142116Sjkh } 1152116Sjkh if(iy>= -1022) { /* normalize output */ 1162116Sjkh hx = ((hx-0x00100000)|((iy+1023)<<20)); 1172116Sjkh INSERT_WORDS(x,hx|sx,lx); 1182116Sjkh } else { /* subnormal output */ 1192116Sjkh n = -1022 - iy; 1202116Sjkh if(n<=20) { 1212116Sjkh lx = (lx>>n)|((u_int32_t)hx<<(32-n)); 1222116Sjkh hx >>= n; 1232116Sjkh } else if (n<=31) { 1242116Sjkh lx = (hx<<(32-n))|(lx>>n); hx = sx; 1252116Sjkh } else { 1262116Sjkh lx = hx>>(n-32); hx = sx; 1272116Sjkh } 1282116Sjkh INSERT_WORDS(x,hx|sx,lx); 1292116Sjkh x *= one; /* create necessary signal */ 1302116Sjkh } 1312116Sjkh return x; /* exact output */ 1322116Sjkh} 133