nsz Git - musl/blob - src/math/remquo.c

   1 /* origin: FreeBSD /usr/src/lib/msun/src/s_remquo.c */
   2 /*-
   3  * ====================================================
   4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   5  *
   6  * Developed at SunSoft, 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  * Return the IEEE remainder and set *quo to the last n bits of the
  14  * quotient, rounded to the nearest integer.  We choose n=31 because
  15  * we wind up computing all the integer bits of the quotient anyway as
  16  * a side-effect of computing the remainder by the shift and subtract
  17  * method.  In practice, this is far more bits than are needed to use
  18  * remquo in reduction algorithms.
  19  */
  20
  21 #include "libm.h"
  22
  23 static const double Zero[] = {0.0, -0.0,};
  24
  25 double remquo(double x, double y, int *quo)
  26 {
  27         int32_t n,hx,hy,hz,ix,iy,sx,i;
  28         uint32_t lx,ly,lz,q,sxy;
  29
  30         EXTRACT_WORDS(hx, lx, x);
  31         EXTRACT_WORDS(hy, ly, y);
  32         sxy = (hx ^ hy) & 0x80000000;
  33         sx = hx & 0x80000000;   /* sign of x */
  34         hx ^= sx;               /* |x| */
  35         hy &= 0x7fffffff;       /* |y| */
  36
  37         /* purge off exception values */
  38         if ((hy|ly) == 0 || hx >= 0x7ff00000 ||  /* y = 0, or x not finite */
  39             (hy|((ly|-ly)>>31)) > 0x7ff00000)    /* or y is NaN */
  40                 return (x*y)/(x*y);
  41         if (hx <= hy) {
  42                 if (hx < hy || lx < ly) {  /* |x| < |y| return x or x-y */
  43                         q = 0;
  44                         goto fixup;
  45                 }
  46                 if (lx == ly) {            /* |x| = |y| return x*0 */
  47                         *quo = sxy ? -1 : 1;
  48                         return Zero[(uint32_t)sx>>31];
  49                 }
  50         }
  51
  52         // FIXME: use ilogb?
  53
  54         /* determine ix = ilogb(x) */
  55         if (hx < 0x00100000) {  /* subnormal x */
  56                 if (hx == 0) {
  57                         for (ix = -1043, i=lx; i>0; i<<=1) ix--;
  58                 } else {
  59                         for (ix = -1022, i=hx<<11; i>0; i<<=1) ix--;
  60                 }
  61         } else
  62                 ix = (hx>>20) - 1023;
  63
  64         /* determine iy = ilogb(y) */
  65         if (hy < 0x00100000) {  /* subnormal y */
  66                 if (hy == 0) {
  67                         for (iy = -1043, i=ly; i>0; i<<=1) iy--;
  68                 } else {
  69                         for (iy = -1022, i=hy<<11; i>0; i<<=1) iy--;
  70                 }
  71         } else
  72                 iy = (hy>>20) - 1023;
  73
  74         /* set up {hx,lx}, {hy,ly} and align y to x */
  75         if (ix >= -1022)
  76                 hx = 0x00100000|(0x000fffff&hx);
  77         else {  /* subnormal x, shift x to normal */
  78                 n = -1022 - ix;
  79                 if (n <= 31) {
  80                         hx = (hx<<n)|(lx>>(32-n));
  81                         lx <<= n;
  82                 } else {
  83                         hx = lx<<(n-32);
  84                         lx = 0;
  85                 }
  86         }
  87         if (iy >= -1022)
  88                 hy = 0x00100000|(0x000fffff&hy);
  89         else {  /* subnormal y, shift y to normal */
  90                 n = -1022 - iy;
  91                 if (n <= 31) {
  92                         hy = (hy<<n)|(ly>>(32-n));
  93                         ly <<= n;
  94                 } else {
  95                         hy = ly<<(n-32);
  96                         ly = 0;
  97                 }
  98         }
  99
 100         /* fix point fmod */
 101         n = ix - iy;
 102         q = 0;
 103         while (n--) {
 104                 hz = hx - hy;
 105                 lz = lx - ly;
 106                 if (lx < ly)
 107                         hz--;
 108                 if (hz < 0) {
 109                         hx = hx + hx + (lx>>31);
 110                         lx = lx + lx;
 111                 } else {
 112                         hx = hz + hz + (lz>>31);
 113                         lx = lz + lz;
 114                         q++;
 115                 }
 116                 q <<= 1;
 117         }
 118         hz = hx - hy;
 119         lz = lx - ly;
 120         if (lx < ly)
 121                 hz--;
 122         if (hz >= 0) {
 123                 hx = hz;
 124                 lx = lz;
 125                 q++;
 126         }
 127
 128         /* convert back to floating value and restore the sign */
 129         if ((hx|lx) == 0) {  /* return sign(x)*0 */
 130                 q &= 0x7fffffff;
 131                 *quo = sxy ? -q : q;
 132                 return Zero[(uint32_t)sx>>31];
 133         }
 134         while (hx < 0x00100000) {  /* normalize x */
 135                 hx = hx + hx + (lx>>31);
 136                 lx = lx + lx;
 137                 iy--;
 138         }
 139         if (iy >= -1022) {         /* normalize output */
 140                 hx = (hx-0x00100000)|((iy+1023)<<20);
 141         } else {                   /* subnormal output */
 142                 n = -1022 - iy;
 143                 if (n <= 20) {
 144                         lx = (lx>>n)|((uint32_t)hx<<(32-n));
 145                         hx >>= n;
 146                 } else if (n <= 31) {
 147                         lx = (hx<<(32-n))|(lx>>n);
 148                         hx = 0;
 149                 } else {
 150                         lx = hx>>(n-32);
 151                         hx = 0;
 152                 }
 153         }
 154 fixup:
 155         INSERT_WORDS(x, hx, lx);
 156         y = fabs(y);
 157         if (y < 0x1p-1021) {
 158                 if (x + x > y || (x + x == y && (q & 1))) {
 159                         q++;
 160                         x -= y;
 161                 }
 162         } else if (x > 0.5*y || (x == 0.5*y && (q & 1))) {
 163                 q++;
 164                 x -= y;
 165         }
 166         GET_HIGH_WORD(hx, x);
 167         SET_HIGH_WORD(x, hx ^ sx);
 168         q &= 0x7fffffff;
 169         *quo = sxy ? -q : q;
 170         return x;
 171 }