libquadmath/math/fmodq.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  * fmodq(x,y)
  17  * Return x mod y in exact arithmetic
  18  * Method: shift and subtract
  19  */
  20
  21 #include "quadmath-imp.h"
  22
  23 static const __float128 one = 1.0, Zero[] = {0.0, -0.0,};
  24
  25 __float128
  26 fmodq (__float128 x, __float128 y)
  27 {
  28   int64_t n,hx,hy,hz,ix,iy,sx,i;
  29   uint64_t lx,ly,lz;
  30
  31   GET_FLT128_WORDS64(hx,lx,x);
  32   GET_FLT128_WORDS64(hy,ly,y);
  33   sx = hx&0x8000000000000000ULL;        /* sign of x */
  34   hx ^=sx;                              /* |x| */
  35   hy &= 0x7fffffffffffffffLL;           /* |y| */
  36
  37   /* purge off exception values */
  38   if((hy|ly)==0||(hx>=0x7fff000000000000LL)|| /* y=0,or x not finite */
  39     ((hy|((ly|-ly)>>63))>0x7fff000000000000LL)) /* or y is NaN */
  40       return (x*y)/(x*y);
  41   if(hx<=hy) {
  42       if((hx<hy)||(lx<ly)) return x;    /* |x|<|y| return x */
  43       if(lx==ly)
  44           return Zero[(uint64_t)sx>>63];        /* |x|=|y| return x*0*/
  45   }
  46
  47   /* determine ix = ilogb(x) */
  48   if(hx<0x0001000000000000LL) { /* subnormal x */
  49       if(hx==0) {
  50           for (ix = -16431, i=lx; i>0; i<<=1) ix -=1;
  51       } else {
  52           for (ix = -16382, i=hx<<15; i>0; i<<=1) ix -=1;
  53       }
  54   } else ix = (hx>>48)-0x3fff;
  55
  56   /* determine iy = ilogb(y) */
  57       if(hy<0x0001000000000000LL) {     /* subnormal y */
  58           if(hy==0) {
  59               for (iy = -16431, i=ly; i>0; i<<=1) iy -=1;
  60           } else {
  61               for (iy = -16382, i=hy<<15; i>0; i<<=1) iy -=1;
  62           }
  63       } else iy = (hy>>48)-0x3fff;
  64
  65   /* set up {hx,lx}, {hy,ly} and align y to x */
  66       if(ix >= -16382)
  67           hx = 0x0001000000000000LL|(0x0000ffffffffffffLL&hx);
  68       else {            /* subnormal x, shift x to normal */
  69           n = -16382-ix;
  70           if(n<=63) {
  71               hx = (hx<<n)|(lx>>(64-n));
  72               lx <<= n;
  73           } else {
  74               hx = lx<<(n-64);
  75               lx = 0;
  76           }
  77       }
  78       if(iy >= -16382)
  79           hy = 0x0001000000000000LL|(0x0000ffffffffffffLL&hy);
  80       else {            /* subnormal y, shift y to normal */
  81           n = -16382-iy;
  82           if(n<=63) {
  83               hy = (hy<<n)|(ly>>(64-n));
  84               ly <<= n;
  85           } else {
  86               hy = ly<<(n-64);
  87               ly = 0;
  88           }
  89       }
  90
  91   /* fix point fmod */
  92       n = ix - iy;
  93       while(n--) {
  94           hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
  95           if(hz<0){hx = hx+hx+(lx>>63); lx = lx+lx;}
  96           else {
  97               if((hz|lz)==0)            /* return sign(x)*0 */
  98                   return Zero[(uint64_t)sx>>63];
  99               hx = hz+hz+(lz>>63); lx = lz+lz;
 100           }
 101       }
 102       hz=hx-hy;lz=lx-ly; if(lx<ly) hz -= 1;
 103       if(hz>=0) {hx=hz;lx=lz;}
 104
 105   /* convert back to floating value and restore the sign */
 106       if((hx|lx)==0)                    /* return sign(x)*0 */
 107           return Zero[(uint64_t)sx>>63];
 108       while(hx<0x0001000000000000LL) {  /* normalize x */
 109           hx = hx+hx+(lx>>63); lx = lx+lx;
 110           iy -= 1;
 111       }
 112       if(iy>= -16382) { /* normalize output */
 113           hx = ((hx-0x0001000000000000LL)|((iy+16383)<<48));
 114           SET_FLT128_WORDS64(x,hx|sx,lx);
 115       } else {          /* subnormal output */
 116           n = -16382 - iy;
 117           if(n<=48) {
 118               lx = (lx>>n)|((uint64_t)hx<<(64-n));
 119               hx >>= n;
 120           } else if (n<=63) {
 121               lx = (hx<<(64-n))|(lx>>n); hx = sx;
 122           } else {
 123               lx = hx>>(n-64); hx = sx;
 124           }
 125           SET_FLT128_WORDS64(x,hx|sx,lx);
 126           x *= one;             /* create necessary signal */
 127       }
 128       return x;         /* exact output */
 129 }