Back to index

glibc  2.9
e_fmodf.c
Go to the documentation of this file.
00001 /* e_fmodf.c -- float version of e_fmod.c.
00002  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
00003  */
00004 
00005 /*
00006  * ====================================================
00007  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00008  *
00009  * Developed at SunPro, a Sun Microsystems, Inc. business.
00010  * Permission to use, copy, modify, and distribute this
00011  * software is freely granted, provided that this notice 
00012  * is preserved.
00013  * ====================================================
00014  */
00015 
00016 #if defined(LIBM_SCCS) && !defined(lint)
00017 static char rcsid[] = "$NetBSD: e_fmodf.c,v 1.4 1995/05/10 20:45:10 jtc Exp $";
00018 #endif
00019 
00020 /* 
00021  * __ieee754_fmodf(x,y)
00022  * Return x mod y in exact arithmetic
00023  * Method: shift and subtract
00024  */
00025 
00026 #include "math.h"
00027 #include "math_private.h"
00028 
00029 #ifdef __STDC__
00030 static const float one = 1.0, Zero[] = {0.0, -0.0,};
00031 #else
00032 static float one = 1.0, Zero[] = {0.0, -0.0,};
00033 #endif
00034 
00035 #ifdef __STDC__
00036        float __ieee754_fmodf(float x, float y)
00037 #else
00038        float __ieee754_fmodf(x,y)
00039        float x,y ;
00040 #endif
00041 {
00042        int32_t n,hx,hy,hz,ix,iy,sx,i;
00043 
00044        GET_FLOAT_WORD(hx,x);
00045        GET_FLOAT_WORD(hy,y);
00046        sx = hx&0x80000000;         /* sign of x */
00047        hx ^=sx;             /* |x| */
00048        hy &= 0x7fffffff;    /* |y| */
00049 
00050     /* purge off exception values */
00051        if(hy==0||(hx>=0x7f800000)||              /* y=0,or x not finite */
00052           (hy>0x7f800000))                /* or y is NaN */
00053            return (x*y)/(x*y);
00054        if(hx<hy) return x;                /* |x|<|y| return x */
00055        if(hx==hy)
00056            return Zero[(u_int32_t)sx>>31];       /* |x|=|y| return x*0*/
00057 
00058     /* determine ix = ilogb(x) */
00059        if(hx<0x00800000) {  /* subnormal x */
00060            for (ix = -126,i=(hx<<8); i>0; i<<=1) ix -=1;
00061        } else ix = (hx>>23)-127;
00062 
00063     /* determine iy = ilogb(y) */
00064        if(hy<0x00800000) {  /* subnormal y */
00065            for (iy = -126,i=(hy<<8); i>=0; i<<=1) iy -=1;
00066        } else iy = (hy>>23)-127;
00067 
00068     /* set up {hx,lx}, {hy,ly} and align y to x */
00069        if(ix >= -126) 
00070            hx = 0x00800000|(0x007fffff&hx);
00071        else {        /* subnormal x, shift x to normal */
00072            n = -126-ix;
00073            hx = hx<<n;
00074        }
00075        if(iy >= -126) 
00076            hy = 0x00800000|(0x007fffff&hy);
00077        else {        /* subnormal y, shift y to normal */
00078            n = -126-iy;
00079            hy = hy<<n;
00080        }
00081 
00082     /* fix point fmod */
00083        n = ix - iy;
00084        while(n--) {
00085            hz=hx-hy;
00086            if(hz<0){hx = hx+hx;}
00087            else {
00088               if(hz==0)            /* return sign(x)*0 */
00089                   return Zero[(u_int32_t)sx>>31];
00090               hx = hz+hz;
00091            }
00092        }
00093        hz=hx-hy;
00094        if(hz>=0) {hx=hz;}
00095 
00096     /* convert back to floating value and restore the sign */
00097        if(hx==0)                   /* return sign(x)*0 */
00098            return Zero[(u_int32_t)sx>>31];       
00099        while(hx<0x00800000) {             /* normalize x */
00100            hx = hx+hx;
00101            iy -= 1;
00102        }
00103        if(iy>= -126) {             /* normalize output */
00104            hx = ((hx-0x00800000)|((iy+127)<<23));
00105            SET_FLOAT_WORD(x,hx|sx);
00106        } else {             /* subnormal output */
00107            n = -126 - iy;
00108            hx >>= n;
00109            SET_FLOAT_WORD(x,hx|sx);
00110            x *= one;        /* create necessary signal */
00111        }
00112        return x;            /* exact output */
00113 }