sysdeps/ieee754/dbl-64/wordsize-64/e_fmod.c

   1 /* Rewritten for 64-bit machines by Ulrich Drepper <drepper@gmail.com>.  */
   2 /*
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunPro, a Sun Microsystems, Inc. business.
   7  * Permission to use, copy, modify, and distribute this
   8  * software is freely granted, provided that this notice
   9  * is preserved.
  10  * ====================================================
  11  */
  12
  13 /*
  14  * __ieee754_fmod(x,y)
  15  * Return x mod y in exact arithmetic
  16  * Method: shift and subtract
  17  */
  18
  19 #include <math.h>
  20 #include <math_private.h>
  21
  22 static const double one = 1.0, Zero[] = {0.0, -0.0,};
  23
  24 double
  25 __ieee754_fmod (double x, double y)
  26 {
  27         int32_t n,ix,iy;
  28         int64_t hx,hy,hz,sx,i;
  29
  30         EXTRACT_WORDS64(hx,x);
  31         EXTRACT_WORDS64(hy,y);
  32         sx = hx&UINT64_C(0x8000000000000000);   /* sign of x */
  33         hx ^=sx;                                /* |x| */
  34         hy &= UINT64_C(0x7fffffffffffffff);     /* |y| */
  35
  36     /* purge off exception values */
  37         if(__builtin_expect(hy==0
  38                             || hx >= UINT64_C(0x7ff0000000000000)
  39                             || hy > UINT64_C(0x7ff0000000000000), 0))
  40           /* y=0,or x not finite or y is NaN */
  41             return (x*y)/(x*y);
  42         if(__builtin_expect(hx<=hy, 0)) {
  43             if(hx<hy) return x; /* |x|<|y| return x */
  44             return Zero[(uint64_t)sx>>63];      /* |x|=|y| return x*0*/
  45         }
  46
  47     /* determine ix = ilogb(x) */
  48         if(__builtin_expect(hx<UINT64_C(0x0010000000000000), 0)) {
  49           /* subnormal x */
  50           for (ix = -1022,i=(hx<<11); i>0; i<<=1) ix -=1;
  51         } else ix = (hx>>52)-1023;
  52
  53     /* determine iy = ilogb(y) */
  54         if(__builtin_expect(hy<UINT64_C(0x0010000000000000), 0)) {      /* subnormal y */
  55           for (iy = -1022,i=(hy<<11); i>0; i<<=1) iy -=1;
  56         } else iy = (hy>>52)-1023;
  57
  58     /* set up hx, hy and align y to x */
  59         if(__builtin_expect(ix >= -1022, 1))
  60             hx = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hx);
  61         else {          /* subnormal x, shift x to normal */
  62             n = -1022-ix;
  63             hx<<=n;
  64         }
  65         if(__builtin_expect(iy >= -1022, 1))
  66             hy = UINT64_C(0x0010000000000000)|(UINT64_C(0x000fffffffffffff)&hy);
  67         else {          /* subnormal y, shift y to normal */
  68             n = -1022-iy;
  69             hy<<=n;
  70         }
  71
  72     /* fix point fmod */
  73         n = ix - iy;
  74         while(n--) {
  75             hz=hx-hy;
  76             if(hz<0){hx = hx+hx;}
  77             else {
  78                 if(hz==0)               /* return sign(x)*0 */
  79                     return Zero[(uint64_t)sx>>63];
  80                 hx = hz+hz;
  81             }
  82         }
  83         hz=hx-hy;
  84         if(hz>=0) {hx=hz;}
  85
  86     /* convert back to floating value and restore the sign */
  87         if(hx==0)                       /* return sign(x)*0 */
  88             return Zero[(uint64_t)sx>>63];
  89         while(hx<UINT64_C(0x0010000000000000)) {        /* normalize x */
  90             hx = hx+hx;
  91             iy -= 1;
  92         }
  93         if(__builtin_expect(iy>= -1022, 1)) {   /* normalize output */
  94           hx = ((hx-UINT64_C(0x0010000000000000))|((uint64_t)(iy+1023)<<52));
  95             INSERT_WORDS64(x,hx|sx);
  96         } else {                /* subnormal output */
  97             n = -1022 - iy;
  98             hx>>=n;
  99             INSERT_WORDS64(x,hx|sx);
 100             x *= one;           /* create necessary signal */
 101         }
 102         return x;               /* exact output */
 103 }
 104 strong_alias (__ieee754_fmod, __fmod_finite)