usr/src/lib/libm/common/m9x/remquof.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 __remquof = remquof
  31
  32 /* INDENT OFF */
  33 /*
  34  * float remquof(float x, float y, int *quo) return remainderf(x,y) and an
  35  * integer pointer quo such that *quo = N mod (2**31),  where N is the
  36  * exact integeral part of x/y rounded to nearest even.
  37  *
  38  * remquof call internal fmodquof
  39  */
  40
  41 #include "libm.h"
  42 #include "libm_protos.h"
  43 #include <math.h>
  44 extern float fabsf(float);
  45
  46 static const int
  47         is = (int) 0x80000000,
  48         im = 0x007fffff,
  49         ii = 0x7f800000,
  50         iu = 0x00800000;
  51
  52 static const float zero = 0.0F, half = 0.5F;
  53 /* INDENT ON */
  54
  55 static float
  56 fmodquof(float x, float y, int *quo) {
  57         float w;
  58         int hx, ix, iy, iz, k, ny, nd, m, sq;
  59
  60         hx = *(int *) &x;
  61         ix = hx & 0x7fffffff;
  62         iy = *(int *) &y;
  63         sq = (iy ^ hx) & is;    /* sign of x/y */
  64         iy &= 0x7fffffff;
  65
  66         /* purge off exception values */
  67         *quo = 0;
  68         if (ix >= ii || iy > ii || iy == 0) {
  69                 w = x * y;
  70                 w = w / w;
  71         } else if (ix <= iy) {
  72                 if (ix < iy)
  73                         w = x;  /* return x if |x|<|y| */
  74                 else {
  75                         *quo = 1 + (sq >> 30);
  76                         w = zero * x;   /* return sign(x)*0.0  */
  77                 }
  78         } else {
  79                 /* INDENT OFF */
  80                 /*
  81                  * scale x,y to "normal" with
  82                  *      ny = exponent of y
  83                  *      nd = exponent of x minus exponent of y
  84                  */
  85                 /* INDENT ON */
  86                 ny = iy >> 23;
  87                 k = ix >> 23;
  88
  89                 /* special case for subnormal y or x */
  90                 if (ny == 0) {
  91                         ny = 1;
  92                         while (iy < iu) {
  93                                 ny -= 1;
  94                                 iy += iy;
  95                         }
  96                         nd = k - ny;
  97                         if (k == 0) {
  98                                 nd += 1;
  99                                 while (ix < iu) {
 100                                         nd -= 1;
 101                                         ix += ix;
 102                                 }
 103                         } else
 104                                 ix = iu | (ix & im);
 105                 } else {
 106                         nd = k - ny;
 107                         ix = iu | (ix & im);
 108                         iy = iu | (iy & im);
 109                 }
 110                 /* INDENT OFF */
 111                 /* fix point fmod for normalized ix and iy */
 112                 /*
 113                  * while (nd--) {
 114                  *      iz = ix - iy;
 115                  *      if (iz < 0)
 116                  *              ix = ix + ix;
 117                  *      else if (iz == 0) {
 118                  *              *(int *) &w = is & hx;
 119                  *              return w;
 120                  *      } else
 121                  *              ix = iz + iz;
 122                  * }
 123                  */
 124                 /* INDENT ON */
 125                 /* unroll the above loop 4 times to gain performance */
 126                 m = 0;
 127                 k = nd >> 2;
 128                 nd -= (k << 2);
 129                 while (k--) {
 130                         iz = ix - iy;
 131                         if (iz >= 0) {
 132                                 m += 1;
 133                                 ix = iz + iz;
 134                         } else
 135                                 ix += ix;
 136                         m += m;
 137                         iz = ix - iy;
 138                         if (iz >= 0) {
 139                                 m += 1;
 140                                 ix = iz + iz;
 141                         } else
 142                                 ix += ix;
 143                         m += m;
 144                         iz = ix - iy;
 145                         if (iz >= 0) {
 146                                 m += 1;
 147                                 ix = iz + iz;
 148                         } else
 149                                 ix += ix;
 150                         m += m;
 151                         iz = ix - iy;
 152                         if (iz >= 0) {
 153                                 m += 1;
 154                                 ix = iz + iz;
 155                         } else
 156                                 ix += ix;
 157                         m += m;
 158                         if (iz == 0) {
 159                                 iz = (k << 2) + nd;
 160                                 if (iz < 32)
 161                                         m <<= iz;
 162                                 else
 163                                         m = 0;
 164                                 m &= 0x7fffffff;
 165                                 *quo = sq >= 0 ? m : -m;
 166                                 *(int *) &w = is & hx;
 167                                 return (w);
 168                         }
 169                 }
 170                 while (nd--) {
 171                         iz = ix - iy;
 172                         if (iz >= 0) {
 173                                 m += 1;
 174                                 ix = iz + iz;
 175                         } else
 176                                 ix += ix;
 177                         m += m;
 178                 }
 179                 /* end of unrolling */
 180
 181                 iz = ix - iy;
 182                 if (iz >= 0) {
 183                         m += 1;
 184                         ix = iz;
 185                 }
 186                 m &= 0x7fffffff;
 187                 *quo = sq >= 0 ? m : -m;
 188
 189                 /* convert back to floating value and restore the sign */
 190                 if (ix == 0) {
 191                         *(int *) &w = is & hx;
 192                         return (w);
 193                 }
 194                 while (ix < iu) {
 195                         ix += ix;
 196                         ny -= 1;
 197                 }
 198                 while (ix > (iu + iu)) {
 199                         ny += 1;
 200                         ix >>= 1;
 201                 }
 202                 if (ny > 0)
 203                         *(int *) &w = (is & hx) | (ix & im) | (ny << 23);
 204                 else {          /* subnormal output */
 205                         k = -ny + 1;
 206                         ix >>= k;
 207                         *(int *) &w = (is & hx) | ix;
 208                 }
 209         }
 210         return (w);
 211 }
 212
 213 float
 214 remquof(float x, float y, int *quo) {
 215         int hx, hy, sx, sq;
 216         float v;
 217
 218         hx = *(int *) &x;       /* high word of x */
 219         hy = *(int *) &y;       /* high word of y */
 220         sx = hx & is;           /* sign of x */
 221         sq = (hx ^ hy) & is;    /* sign of x/y */
 222         hx ^= sx;               /* |x| */
 223         hy &= 0x7fffffff;       /* |y| */
 224
 225         /* purge off exception values: y is 0 or NaN, x is Inf or NaN */
 226         *quo = 0;
 227         if (hx >= ii || hy > ii || hy == 0) {
 228                 v = x * y;
 229                 return (v / v);
 230         }
 231
 232         y = fabsf(y);
 233         x = fabsf(x);
 234         if (hy <= 0x7f7fffff) {
 235                 x = fmodquof(x, y + y, quo);
 236                 *quo = ((*quo) & 0x3fffffff) << 1;
 237         }
 238         if (hy < 0x01000000) {
 239                 if (x + x > y) {
 240                         *quo += 1;
 241                         if (x == y)
 242                                 x = zero;
 243                         else
 244                                 x -= y;
 245                         if (x + x >= y) {
 246                                 x -= y;
 247                                 *quo += 1;
 248                         }
 249                 }
 250         } else {
 251                 v = half * y;
 252                 if (x > v) {
 253                         *quo += 1;
 254                         if (x == y)
 255                                 x = zero;
 256                         else
 257                                 x -= y;
 258                         if (x >= v) {
 259                                 x -= y;
 260                                 *quo += 1;
 261                         }
 262                 }
 263         }
 264         if (sq != 0)
 265                 *quo = -(*quo);
 266         return (sx == 0 ? x : -x);
 267 }