nsz Git - libm/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         // FIXME: signed shift
  39         if ((hy|ly) == 0 || hx >= 0x7ff00000 ||  /* y = 0, or x not finite */
  40             (hy|((ly|-ly)>>31)) > 0x7ff00000)    /* or y is NaN */
  41                 return (x*y)/(x*y);
  42         if (hx <= hy) {
  43                 if (hx < hy || lx < ly) {  /* |x| < |y| return x or x-y */
  44                         q = 0;
  45                         goto fixup;
  46                 }
  47                 if (lx == ly) {            /* |x| = |y| return x*0 */
  48                         *quo = 1;
  49                         return Zero[(uint32_t)sx>>31];
  50                 }
  51         }
  52
  53         // FIXME: use ilogb?
  54
  55         /* determine ix = ilogb(x) */
  56         if (hx < 0x00100000) {  /* subnormal x */
  57                 if (hx == 0) {
  58                         for (ix = -1043, i=lx; i>0; i<<=1) ix--;
  59                 } else {
  60                         for (ix = -1022, i=hx<<11; i>0; i<<=1) ix--;
  61                 }
  62         } else
  63                 ix = (hx>>20) - 1023;
  64
  65         /* determine iy = ilogb(y) */
  66         if (hy < 0x00100000) {  /* subnormal y */
  67                 if (hy == 0) {
  68                         for (iy = -1043, i=ly; i>0; i<<=1) iy--;
  69                 } else {
  70                         for (iy = -1022, i=hy<<11; i>0; i<<=1) iy--;
  71                 }
  72         } else
  73                 iy = (hy>>20) - 1023;
  74
  75         /* set up {hx,lx}, {hy,ly} and align y to x */
  76         if (ix >= -1022)
  77                 hx = 0x00100000|(0x000fffff&hx);
  78         else {  /* subnormal x, shift x to normal */
  79                 n = -1022 - ix;
  80                 if (n <= 31) {
  81                         hx = (hx<<n)|(lx>>(32-n));
  82                         lx <<= n;
  83                 } else {
  84                         hx = lx<<(n-32);
  85                         lx = 0;
  86                 }
  87         }
  88         if (iy >= -1022)
  89                 hy = 0x00100000|(0x000fffff&hy);
  90         else {  /* subnormal y, shift y to normal */
  91                 n = -1022 - iy;
  92                 if (n <= 31) {
  93                         hy = (hy<<n)|(ly>>(32-n));
  94                         ly <<= n;
  95                 } else {
  96                         hy = ly<<(n-32);
  97                         ly = 0;
  98                 }
  99         }
 100
 101         /* fix point fmod */
 102         n = ix - iy;
 103         q = 0;
 104         while (n--) {
 105                 hz = hx - hy;
 106                 lz = lx - ly;
 107                 if (lx < ly)
 108                         hz--;
 109                 if (hz < 0) {
 110                         hx = hx + hx + (lx>>31);
 111                         lx = lx + lx;
 112                 } else {
 113                         hx = hz + hz + (lz>>31);
 114                         lx = lz + lz;
 115                         q++;
 116                 }
 117                 q <<= 1;
 118         }
 119         hz = hx - hy;
 120         lz = lx - ly;
 121         if (lx < ly)
 122                 hz--;
 123         if (hz >= 0) {
 124                 hx = hz;
 125                 lx = lz;
 126                 q++;
 127         }
 128
 129         /* convert back to floating value and restore the sign */
 130         if ((hx|lx) == 0) {  /* return sign(x)*0 */
 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 = sx;
 149                 } else {
 150                         lx = hx>>(n-32);
 151                         hx = sx;
 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 }