12116Sjkh/* e_fmodf.c -- float version of e_fmod.c. 22116Sjkh * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 32116Sjkh */ 42116Sjkh 52116Sjkh/* 62116Sjkh * ==================================================== 72116Sjkh * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 82116Sjkh * 92116Sjkh * Developed at SunPro, a Sun Microsystems, Inc. business. 102116Sjkh * Permission to use, copy, modify, and distribute this 118870Srgrimes * software is freely granted, provided that this notice 122116Sjkh * is preserved. 132116Sjkh * ==================================================== 142116Sjkh */ 152116Sjkh 16176451Sdas#include <sys/cdefs.h> 17176451Sdas__FBSDID("$FreeBSD$"); 182116Sjkh 198870Srgrimes/* 202116Sjkh * __ieee754_fmodf(x,y) 212116Sjkh * Return x mod y in exact arithmetic 222116Sjkh * Method: shift and subtract 232116Sjkh */ 242116Sjkh 252116Sjkh#include "math.h" 262116Sjkh#include "math_private.h" 272116Sjkh 282116Sjkhstatic const float one = 1.0, Zero[] = {0.0, -0.0,}; 292116Sjkh 3097407Salfredfloat 3197407Salfred__ieee754_fmodf(float x, float y) 322116Sjkh{ 332116Sjkh int32_t n,hx,hy,hz,ix,iy,sx,i; 342116Sjkh 352116Sjkh GET_FLOAT_WORD(hx,x); 362116Sjkh GET_FLOAT_WORD(hy,y); 372116Sjkh sx = hx&0x80000000; /* sign of x */ 382116Sjkh hx ^=sx; /* |x| */ 392116Sjkh hy &= 0x7fffffff; /* |y| */ 402116Sjkh 412116Sjkh /* purge off exception values */ 422116Sjkh if(hy==0||(hx>=0x7f800000)|| /* y=0,or x not finite */ 432116Sjkh (hy>0x7f800000)) /* or y is NaN */ 442116Sjkh return (x*y)/(x*y); 452116Sjkh if(hx<hy) return x; /* |x|<|y| return x */ 462116Sjkh if(hx==hy) 472116Sjkh return Zero[(u_int32_t)sx>>31]; /* |x|=|y| return x*0*/ 482116Sjkh 492116Sjkh /* determine ix = ilogb(x) */ 502116Sjkh if(hx<0x00800000) { /* subnormal x */ 512116Sjkh for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1; 522116Sjkh } else ix = (hx>>23)-127; 532116Sjkh 542116Sjkh /* determine iy = ilogb(y) */ 552116Sjkh if(hy<0x00800000) { /* subnormal y */ 562116Sjkh for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1; 572116Sjkh } else iy = (hy>>23)-127; 582116Sjkh 592116Sjkh /* set up {hx,lx}, {hy,ly} and align y to x */ 608870Srgrimes if(ix >= -126) 612116Sjkh hx = 0x00800000|(0x007fffff&hx); 622116Sjkh else { /* subnormal x, shift x to normal */ 632116Sjkh n = -126-ix; 642116Sjkh hx = hx<<n; 652116Sjkh } 668870Srgrimes if(iy >= -126) 672116Sjkh hy = 0x00800000|(0x007fffff&hy); 682116Sjkh else { /* subnormal y, shift y to normal */ 692116Sjkh n = -126-iy; 702116Sjkh hy = hy<<n; 712116Sjkh } 722116Sjkh 732116Sjkh /* fix point fmod */ 742116Sjkh n = ix - iy; 752116Sjkh while(n--) { 762116Sjkh hz=hx-hy; 772116Sjkh if(hz<0){hx = hx+hx;} 782116Sjkh else { 792116Sjkh if(hz==0) /* return sign(x)*0 */ 802116Sjkh return Zero[(u_int32_t)sx>>31]; 812116Sjkh hx = hz+hz; 822116Sjkh } 832116Sjkh } 842116Sjkh hz=hx-hy; 852116Sjkh if(hz>=0) {hx=hz;} 862116Sjkh 872116Sjkh /* convert back to floating value and restore the sign */ 882116Sjkh if(hx==0) /* return sign(x)*0 */ 898870Srgrimes return Zero[(u_int32_t)sx>>31]; 902116Sjkh while(hx<0x00800000) { /* normalize x */ 912116Sjkh hx = hx+hx; 922116Sjkh iy -= 1; 932116Sjkh } 942116Sjkh if(iy>= -126) { /* normalize output */ 952116Sjkh hx = ((hx-0x00800000)|((iy+127)<<23)); 962116Sjkh SET_FLOAT_WORD(x,hx|sx); 972116Sjkh } else { /* subnormal output */ 982116Sjkh n = -126 - iy; 992116Sjkh hx >>= n; 1002116Sjkh SET_FLOAT_WORD(x,hx|sx); 1012116Sjkh x *= one; /* create necessary signal */ 1022116Sjkh } 1032116Sjkh return x; /* exact output */ 1042116Sjkh} 105