src/math/remquol.c

   1 #include "libm.h"
   2
   3 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
   4 long double remquol(long double x, long double y, int *quo)
   5 {
   6         return remquo(x, y, quo);
   7 }
   8 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
   9 long double remquol(long double x, long double y, int *quo)
  10 {
  11         union ldshape ux = {x}, uy = {y};
  12         int ex = ux.i.se & 0x7fff;
  13         int ey = uy.i.se & 0x7fff;
  14         int sx = ux.i.se >> 15;
  15         int sy = uy.i.se >> 15;
  16         uint32_t q;
  17
  18         *quo = 0;
  19         if (y == 0 || isnan(y) || ex == 0x7fff)
  20                 return (x*y)/(x*y);
  21         if (x == 0)
  22                 return x;
  23
  24         /* normalize x and y */
  25         if (!ex) {
  26                 ux.i.se = ex;
  27                 ux.f *= 0x1p120f;
  28                 ex = ux.i.se - 120;
  29         }
  30         if (!ey) {
  31                 uy.i.se = ey;
  32                 uy.f *= 0x1p120f;
  33                 ey = uy.i.se - 120;
  34         }
  35
  36         q = 0;
  37         if (ex >= ey) {
  38                 /* x mod y */
  39 #if LDBL_MANT_DIG == 64
  40                 uint64_t i, mx, my;
  41                 mx = ux.i.m;
  42                 my = uy.i.m;
  43                 for (; ex > ey; ex--) {
  44                         i = mx - my;
  45                         if (mx >= my) {
  46                                 mx = 2*i;
  47                                 q++;
  48                                 q <<= 1;
  49                         } else if (2*mx < mx) {
  50                                 mx = 2*mx - my;
  51                                 q <<= 1;
  52                                 q++;
  53                         } else {
  54                                 mx = 2*mx;
  55                                 q <<= 1;
  56                         }
  57                 }
  58                 i = mx - my;
  59                 if (mx >= my) {
  60                         mx = i;
  61                         q++;
  62                 }
  63                 if (mx == 0)
  64                         ex = -120;
  65                 else
  66                         for (; mx >> 63 == 0; mx *= 2, ex--);
  67                 ux.i.m = mx;
  68 #elif LDBL_MANT_DIG == 113
  69                 uint64_t hi, lo, xhi, xlo, yhi, ylo;
  70                 xhi = (ux.i2.hi & -1ULL>>16) | 1ULL<<48;
  71                 yhi = (uy.i2.hi & -1ULL>>16) | 1ULL<<48;
  72                 xlo = ux.i2.lo;
  73                 ylo = ux.i2.lo;
  74                 for (; ex > ey; ex--) {
  75                         hi = xhi - yhi;
  76                         lo = xlo - ylo;
  77                         if (xlo < ylo)
  78                                 hi -= 1;
  79                         if (hi >> 63 == 0) {
  80                                 xhi = 2*hi + (lo>>63);
  81                                 xlo = 2*lo;
  82                                 q++;
  83                         } else {
  84                                 xhi = 2*xhi + (xlo>>63);
  85                                 xlo = 2*xlo;
  86                         }
  87                         q <<= 1;
  88                 }
  89                 hi = xhi - yhi;
  90                 lo = xlo - ylo;
  91                 if (xlo < ylo)
  92                         hi -= 1;
  93                 if (hi >> 63 == 0) {
  94                         xhi = hi;
  95                         xlo = lo;
  96                         q++;
  97                 }
  98                 if ((xhi|xlo) == 0)
  99                         ex = -120;
 100                 else
 101                         for (; xhi >> 48 == 0; xhi = 2*xhi + (xlo>>63), xlo = 2*xlo, ex--);
 102                 ux.i2.hi = xhi;
 103                 ux.i2.lo = xlo;
 104 #endif
 105         }
 106
 107         /* scale result and decide between |x| and |x|-|y| */
 108         if (ex <= 0) {
 109                 ux.i.se = ex + 120;
 110                 ux.f *= 0x1p-120f;
 111         } else
 112                 ux.i.se = ex;
 113         x = ux.f;
 114         if (sy)
 115                 y = -y;
 116         if (ex == ey || (ex+1 == ey && (2*x > y || (2*x == y && q%2)))) {
 117                 x -= y;
 118                 q++;
 119         }
 120         q &= 0x7fffffff;
 121         *quo = sx^sy ? -(int)q : (int)q;
 122         return sx ? -x : x;
 123 }
 124 #endif