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