Back to index

glibc  2.9
e_fmodl.c
Go to the documentation of this file.
00001 /* e_fmodl.c -- long double version of e_fmod.c.
00002  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
00003  */
00004 /*
00005  * ====================================================
00006  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
00007  *
00008  * Developed at SunPro, a Sun Microsystems, Inc. business.
00009  * Permission to use, copy, modify, and distribute this
00010  * software is freely granted, provided that this notice
00011  * is preserved.
00012  * ====================================================
00013  */
00014 
00015 /*
00016  * __ieee754_fmodl(x,y)
00017  * Return x mod y in exact arithmetic
00018  * Method: shift and subtract
00019  */
00020 
00021 #include "math.h"
00022 #include "math_private.h"
00023 #include <ieee754.h>
00024 
00025 #ifdef __STDC__
00026 static const long double one = 1.0, Zero[] = {0.0, -0.0,};
00027 #else
00028 static long double one = 1.0, Zero[] = {0.0, -0.0,};
00029 #endif
00030 
00031 #ifdef __STDC__
00032        long double __ieee754_fmodl(long double x, long double y)
00033 #else
00034        long double __ieee754_fmodl(x,y)
00035        long double x,y;
00036 #endif
00037 {
00038        int64_t n,hx,hy,hz,ix,iy,sx,i;
00039        u_int64_t lx,ly,lz;
00040        int temp;
00041 
00042        GET_LDOUBLE_WORDS64(hx,lx,x);
00043        GET_LDOUBLE_WORDS64(hy,ly,y);
00044        sx = hx&0x8000000000000000ULL;            /* sign of x */
00045        hx ^=sx;                           /* |x| */
00046        hy &= 0x7fffffffffffffffLL;        /* |y| */
00047 
00048     /* purge off exception values */
00049        if((hy|(ly&0x7fffffffffffffff))==0||(hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
00050          (hy>0x7ff0000000000000LL))       /* 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_int64_t)sx>>63];    /* |x|=|y| return x*0*/
00056        }
00057 
00058     /* determine ix = ilogb(x) */
00059        if(hx<0x0010000000000000LL) {      /* 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<<19; i>0; i<<=1) ix -=1;
00064            }
00065        } else ix = (hx>>52)-0x3ff;
00066 
00067     /* determine iy = ilogb(y) */
00068        if(hy<0x0010000000000000LL) {      /* 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<<19; i>0; i<<=1) iy -=1;
00073            }
00074        } else iy = (hy>>52)-0x3ff;
00075 
00076     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
00077        bit mantissa so the following operatations will give the correct
00078        result.  */
00079         ldbl_extract_mantissa(&hx, &lx, &temp, x);
00080         ldbl_extract_mantissa(&hy, &ly, &temp, y);
00081 
00082     /* set up {hx,lx}, {hy,ly} and align y to x */
00083        if(ix >= -1022)
00084            hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
00085        else {        /* subnormal x, shift x to normal */
00086            n = -1022-ix;
00087            if(n<=63) {
00088                hx = (hx<<n)|(lx>>(64-n));
00089                lx <<= n;
00090            } else {
00091               hx = lx<<(n-64);
00092               lx = 0;
00093            }
00094        }
00095        if(iy >= -1022)
00096            hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
00097        else {        /* subnormal y, shift y to normal */
00098            n = -1022-iy;
00099            if(n<=63) {
00100                hy = (hy<<n)|(ly>>(64-n));
00101                ly <<= n;
00102            } else {
00103               hy = ly<<(n-64);
00104               ly = 0;
00105            }
00106        }
00107 
00108     /* fix point fmod */
00109        n = ix - iy;
00110        while(n--) {
00111            hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
00112            if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
00113            else {
00114               if((hz|(lz&0x7fffffffffffffff))==0)              /* return sign(x)*0 */
00115                   return Zero[(u_int64_t)sx>>63];
00116               hx = hz+hz+(lz>>63); lx = lz+lz;
00117            }
00118        }
00119        hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
00120        if(hz>=0) {hx=hz;lx=lz;}
00121 
00122     /* convert back to floating value and restore the sign */
00123        if((hx|(lx&0x7fffffffffffffff))==0)                     /* return sign(x)*0 */
00124            return Zero[(u_int64_t)sx>>63];
00125        while(hx<0x0001000000000000LL) {   /* normalize x */
00126            hx = hx+hx+(lx>>63); lx = lx+lx;
00127            iy -= 1;
00128        }
00129        if(iy>= -1022) {     /* normalize output */
00130            x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
00131        } else {             /* subnormal output */
00132            n = -1022 - iy;
00133            if(n<=48) {
00134               lx = (lx>>n)|((u_int64_t)hx<<(64-n));
00135               hx >>= n;
00136            } else if (n<=63) {
00137               lx = (hx<<(64-n))|(lx>>n); hx = sx;
00138            } else {
00139               lx = hx>>(n-64); hx = sx;
00140            }
00141            x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
00142            x *= one;        /* create necessary signal */
00143        }
00144        return x;            /* exact output */
00145 }