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 hx, hy, hz, sx, sy;
  31         uint64_t lx, ly, lz;
  32         int n, ix, iy;
  33         double xhi, xlo, yhi, ylo;
  34
  35         ldbl_unpack (x, &xhi, &xlo);
  36         EXTRACT_WORDS64 (hx, xhi);
  37         EXTRACT_WORDS64 (lx, xlo);
  38         ldbl_unpack (y, &yhi, &ylo);
  39         EXTRACT_WORDS64 (hy, yhi);
  40         EXTRACT_WORDS64 (ly, ylo);
  41         sx = hx&0x8000000000000000ULL;          /* sign of x */
  42         hx ^= sx;                               /* |x| */
  43         sy = hy&0x8000000000000000ULL;          /* sign of y */
  44         hy ^= sy;                               /* |y| */
  45
  46     /* purge off exception values */
  47         if(__builtin_expect(hy==0 ||
  48                             (hx>=0x7ff0000000000000LL)|| /* y=0,or x not finite */
  49                             (hy>0x7ff0000000000000LL),0))       /* or y is NaN */
  50             return (x*y)/(x*y);
  51         if (__builtin_expect (hx <= hy, 0))
  52           {
  53             /* If |x| < |y| return x.  */
  54             if (hx < hy)
  55               return x;
  56             /* At this point the absolute value of the high doubles of
  57                x and y must be equal.  */
  58             /* If the low double of y is the same sign as the high
  59                double of y (ie. the low double increases |y|)...  */
  60             if (((ly ^ sy) & 0x8000000000000000LL) == 0
  61                 /* ... then a different sign low double to high double
  62                    for x or same sign but lower magnitude...  */
  63                 && (int64_t) (lx ^ sx) < (int64_t) (ly ^ sy))
  64               /* ... means |x| < |y|.  */
  65               return x;
  66             /* If the low double of x differs in sign to the high
  67                double of x (ie. the low double decreases |x|)...  */
  68             if (((lx ^ sx) & 0x8000000000000000LL) != 0
  69                 /* ... then a different sign low double to high double
  70                    for y with lower magnitude (we've already caught
  71                    the same sign for y case above)...  */
  72                 && (int64_t) (lx ^ sx) > (int64_t) (ly ^ sy))
  73               /* ... means |x| < |y|.  */
  74               return x;
  75             /* If |x| == |y| return x*0.  */
  76             if ((lx ^ sx) == (ly ^ sy))
  77               return Zero[(uint64_t) sx >> 63];
  78         }
  79
  80     /* Make the IBM extended format 105 bit mantissa look like the ieee854 112
  81        bit mantissa so the following operations will give the correct
  82        result.  */
  83         ldbl_extract_mantissa(&hx, &lx, &ix, x);
  84         ldbl_extract_mantissa(&hy, &ly, &iy, y);
  85
  86         if (__builtin_expect (ix == -IEEE754_DOUBLE_BIAS, 0))
  87           {
  88             /* subnormal x, shift x to normal.  */
  89             while ((hx & (1LL << 48)) == 0)
  90               {
  91                 hx = (hx << 1) | (lx >> 63);
  92                 lx = lx << 1;
  93                 ix -= 1;
  94               }
  95           }
  96
  97         if (__builtin_expect (iy == -IEEE754_DOUBLE_BIAS, 0))
  98           {
  99             /* subnormal y, shift y to normal.  */
 100             while ((hy & (1LL << 48)) == 0)
 101               {
 102                 hy = (hy << 1) | (ly >> 63);
 103                 ly = ly << 1;
 104                 iy -= 1;
 105               }
 106           }
 107
 108     /* fix point fmod */
 109         n = ix - iy;
 110         while(n--) {
 111             hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 112             if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
 113             else {
 114                 if((hz|lz)==0)          /* return sign(x)*0 */
 115                     return Zero[(u_int64_t)sx>>63];
 116                 hx = hz+hz+(lz>>63); lx = lz+lz;
 117             }
 118         }
 119         hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 120         if(hz>=0) {hx=hz;lx=lz;}
 121
 122     /* convert back to floating value and restore the sign */
 123         if((hx|lx)==0)                  /* return sign(x)*0 */
 124             return Zero[(u_int64_t)sx>>63];
 125         while(hx<0x0001000000000000LL) {        /* normalize x */
 126             hx = hx+hx+(lx>>63); lx = lx+lx;
 127             iy -= 1;
 128         }
 129         if(__builtin_expect(iy>= -1022,0)) {    /* normalize output */
 130             x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
 131         } else {                /* subnormal output */
 132             n = -1022 - iy;
 133             if(n<=48) {
 134                 lx = (lx>>n)|((u_int64_t)hx<<(64-n));
 135                 hx >>= n;
 136             } else if (n<=63) {
 137                 lx = (hx<<(64-n))|(lx>>n); hx = sx;
 138             } else {
 139                 lx = hx>>(n-64); hx = sx;
 140             }
 141             x = ldbl_insert_mantissa((sx>>63), iy, hx, lx);
 142             x *= one;           /* create necessary signal */
 143         }
 144         return x;               /* exact output */
 145 }
 146 strong_alias (__ieee754_fmodl, __fmodl_finite)