sysdeps/ieee754/dbl-64/e_fmod.c

   1 /*
   2  * ====================================================
   3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   4  *
   5  * Developed at SunPro, a Sun Microsystems, Inc. business.
   6  * Permission to use, copy, modify, and distribute this
   7  * software is freely granted, provided that this notice
   8  * is preserved.
   9  * ====================================================
  10  */
  11
  12 /*
  13  * __ieee754_fmod(x,y)
  14  * Return x mod y in exact arithmetic
  15  * Method: shift and subtract
  16  */
  17
  18 #include <math.h>
  19 #include <math_private.h>
  20
  21 static const double one = 1.0, Zero[] = { 0.0, -0.0, };
  22
  23 double
  24 __ieee754_fmod (double x, double y)
  25 {
  26   int32_t n, hx, hy, hz, ix, iy, sx, i;
  27   u_int32_t lx, ly, lz;
  28
  29   EXTRACT_WORDS (hx, lx, x);
  30   EXTRACT_WORDS (hy, ly, y);
  31   sx = hx & 0x80000000;                 /* sign of x */
  32   hx ^= sx;                     /* |x| */
  33   hy &= 0x7fffffff;             /* |y| */
  34
  35   /* purge off exception values */
  36   if ((hy | ly) == 0 || (hx >= 0x7ff00000) ||   /* y=0,or x not finite */
  37       ((hy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* or y is NaN */
  38     return (x * y) / (x * y);
  39   if (hx <= hy)
  40     {
  41       if ((hx < hy) || (lx < ly))
  42         return x;                               /* |x|<|y| return x */
  43       if (lx == ly)
  44         return Zero[(u_int32_t) sx >> 31];      /* |x|=|y| return x*0*/
  45     }
  46
  47   /* determine ix = ilogb(x) */
  48   if (__glibc_unlikely (hx < 0x00100000))                  /* subnormal x */
  49     {
  50       if (hx == 0)
  51         {
  52           for (ix = -1043, i = lx; i > 0; i <<= 1)
  53             ix -= 1;
  54         }
  55       else
  56         {
  57           for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
  58             ix -= 1;
  59         }
  60     }
  61   else
  62     ix = (hx >> 20) - 1023;
  63
  64   /* determine iy = ilogb(y) */
  65   if (__glibc_unlikely (hy < 0x00100000))                  /* subnormal y */
  66     {
  67       if (hy == 0)
  68         {
  69           for (iy = -1043, i = ly; i > 0; i <<= 1)
  70             iy -= 1;
  71         }
  72       else
  73         {
  74           for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
  75             iy -= 1;
  76         }
  77     }
  78   else
  79     iy = (hy >> 20) - 1023;
  80
  81   /* set up {hx,lx}, {hy,ly} and align y to x */
  82   if (__glibc_likely (ix >= -1022))
  83     hx = 0x00100000 | (0x000fffff & hx);
  84   else                  /* subnormal x, shift x to normal */
  85     {
  86       n = -1022 - ix;
  87       if (n <= 31)
  88         {
  89           hx = (hx << n) | (lx >> (32 - n));
  90           lx <<= n;
  91         }
  92       else
  93         {
  94           hx = lx << (n - 32);
  95           lx = 0;
  96         }
  97     }
  98   if (__glibc_likely (iy >= -1022))
  99     hy = 0x00100000 | (0x000fffff & hy);
 100   else                  /* subnormal y, shift y to normal */
 101     {
 102       n = -1022 - iy;
 103       if (n <= 31)
 104         {
 105           hy = (hy << n) | (ly >> (32 - n));
 106           ly <<= n;
 107         }
 108       else
 109         {
 110           hy = ly << (n - 32);
 111           ly = 0;
 112         }
 113     }
 114
 115   /* fix point fmod */
 116   n = ix - iy;
 117   while (n--)
 118     {
 119       hz = hx - hy; lz = lx - ly; if (lx < ly)
 120         hz -= 1;
 121       if (hz < 0)
 122         {
 123           hx = hx + hx + (lx >> 31); lx = lx + lx;
 124         }
 125       else
 126         {
 127           if ((hz | lz) == 0)           /* return sign(x)*0 */
 128             return Zero[(u_int32_t) sx >> 31];
 129           hx = hz + hz + (lz >> 31); lx = lz + lz;
 130         }
 131     }
 132   hz = hx - hy; lz = lx - ly; if (lx < ly)
 133     hz -= 1;
 134   if (hz >= 0)
 135     {
 136       hx = hz; lx = lz;
 137     }
 138
 139   /* convert back to floating value and restore the sign */
 140   if ((hx | lx) == 0)                   /* return sign(x)*0 */
 141     return Zero[(u_int32_t) sx >> 31];
 142   while (hx < 0x00100000)               /* normalize x */
 143     {
 144       hx = hx + hx + (lx >> 31); lx = lx + lx;
 145       iy -= 1;
 146     }
 147   if (__glibc_likely (iy >= -1022))              /* normalize output */
 148     {
 149       hx = ((hx - 0x00100000) | ((iy + 1023) << 20));
 150       INSERT_WORDS (x, hx | sx, lx);
 151     }
 152   else                          /* subnormal output */
 153     {
 154       n = -1022 - iy;
 155       if (n <= 20)
 156         {
 157           lx = (lx >> n) | ((u_int32_t) hx << (32 - n));
 158           hx >>= n;
 159         }
 160       else if (n <= 31)
 161         {
 162           lx = (hx << (32 - n)) | (lx >> n); hx = sx;
 163         }
 164       else
 165         {
 166           lx = hx >> (n - 32); hx = sx;
 167         }
 168       INSERT_WORDS (x, hx | sx, lx);
 169       x *= one;                 /* create necessary signal */
 170     }
 171   return x;                     /* exact output */
 172 }
 173 strong_alias (__ieee754_fmod, __fmod_finite)