usr/src/lib/libm/common/m9x/remquo.c

   1 /*
   2  * CDDL HEADER START
   3  *
   4  * The contents of this file are subject to the terms of the
   5  * Common Development and Distribution License (the "License").
   6  * You may not use this file except in compliance with the License.
   7  *
   8  * You can obtain a copy of the license at usr/src/OPENSOLARIS.LICENSE
   9  * or http://www.opensolaris.org/os/licensing.
  10  * See the License for the specific language governing permissions
  11  * and limitations under the License.
  12  *
  13  * When distributing Covered Code, include this CDDL HEADER in each
  14  * file and include the License file at usr/src/OPENSOLARIS.LICENSE.
  15  * If applicable, add the following below this CDDL HEADER, with the
  16  * fields enclosed by brackets "[]" replaced with your own identifying
  17  * information: Portions Copyright [yyyy] [name of copyright owner]
  18  *
  19  * CDDL HEADER END
  20  */
  21
  22 /*
  23  * Copyright 2011 Nexenta Systems, Inc.  All rights reserved.
  24  */
  25 /*
  26  * Copyright 2006 Sun Microsystems, Inc.  All rights reserved.
  27  * Use is subject to license terms.
  28  */
  29
  30 #pragma weak __remquo = remquo
  31
  32 /* INDENT OFF */
  33 /*
  34  * double remquo(double x, double y, int *quo) return remainder(x,y) and an
  35  * integer pointer quo such that *quo = N mod {2**31}, where N is the
  36  * exact integral part of x/y rounded to nearest even.
  37  *
  38  * remquo call internal fmodquo
  39  */
  40 /* INDENT ON */
  41
  42 #include "libm.h"
  43 #include "libm_protos.h"
  44 #include <math.h>               /* fabs() */
  45 #include <sys/isa_defs.h>
  46
  47 #if defined(_BIG_ENDIAN)
  48 #define HIWORD  0
  49 #define LOWORD  1
  50 #else
  51 #define HIWORD  1
  52 #define LOWORD  0
  53 #endif
  54 #define __HI(x) ((int *) &x)[HIWORD]
  55 #define __LO(x) ((int *) &x)[LOWORD]
  56
  57 static const double one = 1.0, Zero[] = {0.0, -0.0};
  58
  59 static double
  60 fmodquo(double x, double y, int *quo) {
  61         int n, hx, hy, hz, ix, iy, sx, sq, i, m;
  62         unsigned lx, ly, lz;
  63
  64         hx = __HI(x);           /* high word of x */
  65         lx = __LO(x);           /* low  word of x */
  66         hy = __HI(y);           /* high word of y */
  67         ly = __LO(y);           /* low  word of y */
  68         sx = hx & 0x80000000;   /* sign of x */
  69         sq = (hx ^ hy) & 0x80000000;    /* sign of x/y */
  70         hx ^= sx;               /* |x| */
  71         hy &= 0x7fffffff;       /* |y| */
  72
  73         /* purge off exception values */
  74         *quo = 0;
  75         if ((hy | ly) == 0 || hx >= 0x7ff00000 ||       /* y=0, or x !finite */
  76             (hy | ((ly | -ly) >> 31)) > 0x7ff00000)     /* or y is NaN */
  77                 return ((x * y) / (x * y));
  78         if (hx <= hy) {
  79                 if (hx < hy || lx < ly)
  80                         return (x);     /* |x|<|y| return x */
  81                 if (lx == ly) {
  82                         *quo = 1 + (sq >> 30);
  83                         /* |x|=|y| return x*0 */
  84                         return (Zero[(unsigned) sx >> 31]);
  85                 }
  86         }
  87
  88         /* determine ix = ilogb(x) */
  89         if (hx < 0x00100000) {  /* subnormal x */
  90                 if (hx == 0) {
  91                         for (ix = -1043, i = lx; i > 0; i <<= 1)
  92                                 ix -= 1;
  93                 } else {
  94                         for (ix = -1022, i = (hx << 11); i > 0; i <<= 1)
  95                                 ix -= 1;
  96                 }
  97         } else
  98                 ix = (hx >> 20) - 1023;
  99
 100         /* determine iy = ilogb(y) */
 101         if (hy < 0x00100000) {  /* subnormal y */
 102                 if (hy == 0) {
 103                         for (iy = -1043, i = ly; i > 0; i <<= 1)
 104                                 iy -= 1;
 105                 } else {
 106                         for (iy = -1022, i = (hy << 11); i > 0; i <<= 1)
 107                                 iy -= 1;
 108                 }
 109         } else
 110                 iy = (hy >> 20) - 1023;
 111
 112         /* set up {hx,lx}, {hy,ly} and align y to x */
 113         if (ix >= -1022)
 114                 hx = 0x00100000 | (0x000fffff & hx);
 115         else {                  /* subnormal x, shift x to normal */
 116                 n = -1022 - ix;
 117                 if (n <= 31) {
 118                         hx = (hx << n) | (lx >> (32 - n));
 119                         lx <<= n;
 120                 } else {
 121                         hx = lx << (n - 32);
 122                         lx = 0;
 123                 }
 124         }
 125         if (iy >= -1022)
 126                 hy = 0x00100000 | (0x000fffff & hy);
 127         else {                  /* subnormal y, shift y to normal */
 128                 n = -1022 - iy;
 129                 if (n <= 31) {
 130                         hy = (hy << n) | (ly >> (32 - n));
 131                         ly <<= n;
 132                 } else {
 133                         hy = ly << (n - 32);
 134                         ly = 0;
 135                 }
 136         }
 137
 138         /* fix point fmod */
 139         n = ix - iy;
 140         m = 0;
 141         while (n--) {
 142                 hz = hx - hy;
 143                 lz = lx - ly;
 144                 if (lx < ly)
 145                         hz -= 1;
 146                 if (hz < 0) {
 147                         hx = hx + hx + (lx >> 31);
 148                         lx = lx + lx;
 149                 } else {
 150                         m += 1;
 151                         if ((hz | lz) == 0) {   /* return sign(x)*0 */
 152                                 if (n < 31)
 153                                         m <<= 1 + n;
 154                                 else
 155                                         m = 0;
 156                                 m &= 0x7fffffff;
 157                                 *quo = sq >= 0 ? m : -m;
 158                                 return (Zero[(unsigned) sx >> 31]);
 159                         }
 160                         hx = hz + hz + (lz >> 31);
 161                         lx = lz + lz;
 162                 }
 163                 m += m;
 164         }
 165         hz = hx - hy;
 166         lz = lx - ly;
 167         if (lx < ly)
 168                 hz -= 1;
 169         if (hz >= 0) {
 170                 hx = hz;
 171                 lx = lz;
 172                 m += 1;
 173         }
 174         m &= 0x7fffffff;
 175         *quo = sq >= 0 ? m : -m;
 176
 177         /* convert back to floating value and restore the sign */
 178         if ((hx | lx) == 0) {   /* return sign(x)*0 */
 179                 return (Zero[(unsigned) sx >> 31]);
 180         }
 181         while (hx < 0x00100000) {       /* normalize x */
 182                 hx = hx + hx + (lx >> 31);
 183                 lx = lx + lx;
 184                 iy -= 1;
 185         }
 186         if (iy >= -1022) {      /* normalize output */
 187                 hx = (hx - 0x00100000) | ((iy + 1023) << 20);
 188                 __HI(x) = hx | sx;
 189                 __LO(x) = lx;
 190         } else {                        /* subnormal output */
 191                 n = -1022 - iy;
 192                 if (n <= 20) {
 193                         lx = (lx >> n) | ((unsigned) hx << (32 - n));
 194                         hx >>= n;
 195                 } else if (n <= 31) {
 196                         lx = (hx << (32 - n)) | (lx >> n);
 197                         hx = sx;
 198                 } else {
 199                         lx = hx >> (n - 32);
 200                         hx = sx;
 201                 }
 202                 __HI(x) = hx | sx;
 203                 __LO(x) = lx;
 204                 x *= one;       /* create necessary signal */
 205         }
 206         return (x);             /* exact output */
 207 }
 208
 209 #define zero    Zero[0]
 210
 211 double
 212 remquo(double x, double y, int *quo) {
 213         int hx, hy, sx, sq;
 214         double v;
 215         unsigned ly;
 216
 217         hx = __HI(x);           /* high word of x */
 218         hy = __HI(y);           /* high word of y */
 219         ly = __LO(y);           /* low  word of y */
 220         sx = hx & 0x80000000;   /* sign of x */
 221         sq = (hx ^ hy) & 0x80000000;    /* sign of x/y */
 222         hx ^= sx;               /* |x| */
 223         hy &= 0x7fffffff;       /* |y| */
 224
 225         /* purge off exception values */
 226         *quo = 0;
 227         if ((hy | ly) == 0 || hx >= 0x7ff00000 ||       /* y=0, or x !finite */
 228             (hy | ((ly | -ly) >> 31)) > 0x7ff00000)     /* or y is NaN */
 229                 return ((x * y) / (x * y));
 230
 231         y = fabs(y);
 232         x = fabs(x);
 233         if (hy <= 0x7fdfffff) {
 234                 x = fmodquo(x, y + y, quo);
 235                 *quo = ((*quo) & 0x3fffffff) << 1;
 236         }
 237         if (hy < 0x00200000) {
 238                 if (x + x > y) {
 239                         *quo += 1;
 240                         if (x == y)
 241                                 x = zero;
 242                         else
 243                                 x -= y;
 244                         if (x + x >= y) {
 245                                 x -= y;
 246                                 *quo += 1;
 247                         }
 248                 }
 249         } else {
 250                 v = 0.5 * y;
 251                 if (x > v) {
 252                         *quo += 1;
 253                         if (x == y)
 254                                 x = zero;
 255                         else
 256                                 x -= y;
 257                         if (x >= v) {
 258                                 x -= y;
 259                                 *quo += 1;
 260                         }
 261                 }
 262         }
 263         if (sq != 0)
 264                 *quo = -(*quo);
 265         return (sx == 0 ? x : -x);
 266 }