Print this page
11210 libm should be cstyle(1ONBLD) clean


   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 }


   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 }