Back to index

glibc  2.9
e_fmod.c
Go to the documentation of this file.
00001 /* @(#)e_fmod.c 5.1 93/09/24 */
00002 /*
00003  * ====================================================
00004  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00005  *
00006  * Developed at SunPro, a Sun Microsystems, Inc. business.
00007  * Permission to use, copy, modify, and distribute this
00008  * software is freely granted, provided that this notice 
00009  * is preserved.
00010  * ====================================================
00011  */
00012 
00013 #if defined(LIBM_SCCS) && !defined(lint)
00014 static char rcsid[] = "$NetBSD: e_fmod.c,v 1.8 1995/05/10 20:45:07 jtc Exp $";
00015 #endif
00016 
00017 /* 
00018  * __ieee754_fmod(x,y)
00019  * Return x mod y in exact arithmetic
00020  * Method: shift and subtract
00021  */
00022 
00023 #include "math.h"
00024 #include "math_private.h"
00025 
00026 #ifdef __STDC__
00027 static const double one = 1.0, Zero[] = {0.0, -0.0,};
00028 #else
00029 static double one = 1.0, Zero[] = {0.0, -0.0,};
00030 #endif
00031 
00032 #ifdef __STDC__
00033        double __ieee754_fmod(double x, double y)
00034 #else
00035        double __ieee754_fmod(x,y)
00036        double x,y ;
00037 #endif
00038 {
00039        int32_t n,hx,hy,hz,ix,iy,sx,i;
00040        u_int32_t lx,ly,lz;
00041 
00042        EXTRACT_WORDS(hx,lx,x);
00043        EXTRACT_WORDS(hy,ly,y);
00044        sx = hx&0x80000000;         /* sign of x */
00045        hx ^=sx;             /* |x| */
00046        hy &= 0x7fffffff;    /* |y| */
00047 
00048     /* purge off exception values */
00049        if((hy|ly)==0||(hx>=0x7ff00000)||  /* y=0,or x not finite */
00050          ((hy|((ly|-ly)>>31))>0x7ff00000))       /* or y is NaN */
00051            return (x*y)/(x*y);
00052        if(hx<=hy) {
00053            if((hx<hy)||(lx<ly)) return x; /* |x|<|y| return x */
00054            if(lx==ly) 
00055               return Zero[(u_int32_t)sx>>31];    /* |x|=|y| return x*0*/
00056        }
00057 
00058     /* determine ix = ilogb(x) */
00059        if(hx<0x00100000) {  /* subnormal x */
00060            if(hx==0) {
00061               for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
00062            } else {
00063               for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
00064            }
00065        } else ix = (hx>>20)-1023;
00066 
00067     /* determine iy = ilogb(y) */
00068        if(hy<0x00100000) {  /* subnormal y */
00069            if(hy==0) {
00070               for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
00071            } else {
00072               for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
00073            }
00074        } else iy = (hy>>20)-1023;
00075 
00076     /* set up {hx,lx}, {hy,ly} and align y to x */
00077        if(ix >= -1022) 
00078            hx = 0x00100000|(0x000fffff&hx);
00079        else {        /* subnormal x, shift x to normal */
00080            n = -1022-ix;
00081            if(n<=31) {
00082                hx = (hx<<n)|(lx>>(32-n));
00083                lx <<= n;
00084            } else {
00085               hx = lx<<(n-32);
00086               lx = 0;
00087            }
00088        }
00089        if(iy >= -1022) 
00090            hy = 0x00100000|(0x000fffff&hy);
00091        else {        /* subnormal y, shift y to normal */
00092            n = -1022-iy;
00093            if(n<=31) {
00094                hy = (hy<<n)|(ly>>(32-n));
00095                ly <<= n;
00096            } else {
00097               hy = ly<<(n-32);
00098               ly = 0;
00099            }
00100        }
00101 
00102     /* fix point fmod */
00103        n = ix - iy;
00104        while(n--) {
00105            hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
00106            if(hz<0){hx = hx+hx+(lx>>31); lx = lx+lx;}
00107            else {
00108               if((hz|lz)==0)              /* return sign(x)*0 */
00109                   return Zero[(u_int32_t)sx>>31];
00110               hx = hz+hz+(lz>>31); lx = lz+lz;
00111            }
00112        }
00113        hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
00114        if(hz>=0) {hx=hz;lx=lz;}
00115 
00116     /* convert back to floating value and restore the sign */
00117        if((hx|lx)==0)                     /* return sign(x)*0 */
00118            return Zero[(u_int32_t)sx>>31];       
00119        while(hx<0x00100000) {             /* normalize x */
00120            hx = hx+hx+(lx>>31); lx = lx+lx;
00121            iy -= 1;
00122        }
00123        if(iy>= -1022) {     /* normalize output */
00124            hx = ((hx-0x00100000)|((iy+1023)<<20));
00125            INSERT_WORDS(x,hx|sx,lx);
00126        } else {             /* subnormal output */
00127            n = -1022 - iy;
00128            if(n<=20) {
00129               lx = (lx>>n)|((u_int32_t)hx<<(32-n));
00130               hx >>= n;
00131            } else if (n<=31) {
00132               lx = (hx<<(32-n))|(lx>>n); hx = sx;
00133            } else {
00134               lx = hx>>(n-32); hx = sx;
00135            }
00136            INSERT_WORDS(x,hx|sx,lx);
00137            x *= one;        /* create necessary signal */
00138        }
00139        return x;            /* exact output */
00140 }