sysdeps/ieee754/ldbl-128ibm/e_fmodl.c

   1 /* e_fmodl.c -- long double version of e_fmod.c.
   2  * Conversion to IEEE quad long double by Jakub Jelinek, jj@ultra.linux.cz.
   3  */
   4 /*
   5  * ====================================================
   6  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   7  *
   8  * Developed at SunPro, a Sun Microsystems, Inc. business.
   9  * Permission to use, copy, modify, and distribute this
  10  * software is freely granted, provided that this notice
  11  * is preserved.
  12  * ====================================================
  13  */
  14
  15 /*
  16  * __ieee754_fmodl(x,y)
  17  * Return x mod y in exact arithmetic
  18  * Method: shift and subtract
  19  */
  20
  21 #include <math.h>
  22 #include <math_private.h>
  23 #include <ieee754.h>
  24
  25 static const long double one = 1.0, Zero[] = {0.0, -0.0,};
  26
  27 long double
  28 __ieee754_fmodl (long double x, long double y)
  29 {
  30         int64_t n,hx,hy,hz,ix,iy,sx, i;
  31         u_int64_t lx,ly,lz;
  32         int temp;
  33
  34         GET_LDOUBLE_WORDS64(hx,lx,x);
  35         GET_LDOUBLE_WORDS64(hy,ly,y);
  36         sx = hx&0x8000000000000000ULL;          /* sign of x */
  37         hx ^=sx;                                /* |x| */
  38         hy &= 0x7fffffffffffffffLL;             /* |y| */
  39
  40     /* purge off exception values */
  41         if(__builtin_expect((hy|(ly&0x7fffffffffffffff))==0 ||
  42                             (hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
  43                             (hy>0x7ff0000000000000LL),0))       /* or y is NaN */
  44             return (x*y)/(x*y);
  45         if(__builtin_expect(hx<=hy,0)) {
  46             if((hx<hy)||(lx<ly)) return x;      /* |x|<|y| return x */
  47             if(lx==ly)
  48                 return Zero[(u_int64_t)sx>>63]; /* |x|=|y| return x*0*/
  49         }
  50
  51     /* determine ix = ilogb(x) */
  52         if(__builtin_expect(hx<0x0010000000000000LL,0)) {       /* subnormal x */
  53             if(hx==0) {
  54                 for (ix = -1043, i=lx; i>0; i<<=1) ix -=1;
  55             } else {
  56                 for (ix = -1022, i=(hx<<11); i>0; i<<=1) ix -=1;
  57             }
  58         } else ix = (hx>>52)-0x3ff;
  59
  60     /* determine iy = ilogb(y) */
  61         if(__builtin_expect(hy<0x0010000000000000LL,0)) {       /* subnormal y */
  62             if(hy==0) {
  63                 for (iy = -1043, i=ly; i>0; i<<=1) iy -=1;
  64             } else {
  65                 for (iy = -1022, i=(hy<<11); i>0; i<<=1) iy -=1;
  66             }
  67         } else iy = (hy>>52)-0x3ff;
  68
  69     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
  70        bit mantissa so the following operations will give the correct
  71        result.  */
  72         ldbl_extract_mantissa(&hx, &lx, &temp, x);
  73         ldbl_extract_mantissa(&hy, &ly, &temp, y);
  74
  75     /* set up {hx,lx}, {hy,ly} and align y to x */
  76         if(__builtin_expect(ix >= -1022, 1))
  77             hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
  78         else {          /* subnormal x, shift x to normal */
  79             n = -1022-ix;
  80             if(n<=63) {
  81                 hx = (hx<<n)|(lx>>(64-n));
  82                 lx <<= n;
  83             } else {
  84                 hx = lx<<(n-64);
  85                 lx = 0;
  86             }
  87         }
  88         if(__builtin_expect(iy >= -1022, 1))
  89             hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
  90         else {          /* subnormal y, shift y to normal */
  91             n = -1022-iy;
  92             if(n<=63) {
  93                 hy = (hy<<n)|(ly>>(64-n));
  94                 ly <<= n;
  95             } else {
  96                 hy = ly<<(n-64);
  97                 ly = 0;
  98             }
  99         }
 100
 101     /* fix point fmod */
 102         n = ix - iy;
 103         while(n--) {
 104             hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 105             if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
 106             else {
 107                 if((hz|(lz&0x7fffffffffffffff))==0)             /* return sign(x)*0 */
 108                     return Zero[(u_int64_t)sx>>63];
 109                 hx = hz+hz+(lz>>63); lx = lz+lz;
 110             }
 111         }
 112         hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 113         if(hz>=0) {hx=hz;lx=lz;}
 114
 115     /* convert back to floating value and restore the sign */
 116         if((hx|(lx&0x7fffffffffffffff))==0)                     /* return sign(x)*0 */
 117             return Zero[(u_int64_t)sx>>63];
 118         while(hx<0x0001000000000000LL) {        /* normalize x */
 119             hx = hx+hx+(lx>>63); lx = lx+lx;
 120             iy -= 1;
 121         }
 122         if(__builtin_expect(iy>= -1022,0)) {    /* normalize output */
 123             x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
 124         } else {                /* subnormal output */
 125             n = -1022 - iy;
 126             if(n<=48) {
 127                 lx = (lx>>n)|((u_int64_t)hx<<(64-n));
 128                 hx >>= n;
 129             } else if (n<=63) {
 130                 lx = (hx<<(64-n))|(lx>>n); hx = sx;
 131             } else {
 132                 lx = hx>>(n-64); hx = sx;
 133             }
 134             x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
 135             x *= one;           /* create necessary signal */
 136         }
 137         return x;               /* exact output */
 138 }
 139 strong_alias (__ieee754_fmodl, __fmodl_finite)