newlib/libm/common/sf_lrint.c

   1
   2 /* @(#)sf_lrint.c 5.1 93/09/24 */
   3 /*
   4  * ====================================================
   5  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
   6  *
   7  * Developed at SunPro, a Sun Microsystems, Inc. business.
   8  * Permission to use, copy, modify, and distribute this
   9  * software is freely granted, provided that this notice
  10  * is preserved.
  11  * ====================================================
  12  */
  13
  14 /*
  15  * lrintf(x)
  16  * Return x rounded to integral value according to the prevailing
  17  * rounding mode.
  18  * Method:
  19  *      Using floating addition.
  20  * Exception:
  21  *      Inexact flag raised if x not equal to lrintf(x).
  22  */
  23
  24 #include "fdlibm.h"
  25
  26 #ifdef __STDC__
  27 static const float
  28 #else
  29 static float
  30 #endif
  31 /* Adding a float, x, to 2^23 will cause the result to be rounded based on
  32    the fractional part of x, according to the implementation's current rounding
  33    mode.  2^23 is the smallest float that can be represented using all 23 significant
  34    digits. */
  35 TWO23[2]={
  36   8.3886080000e+06, /* 0x4b000000 */
  37  -8.3886080000e+06, /* 0xcb000000 */
  38 };
  39
  40 #ifdef __STDC__
  41         long int lrintf(float x)
  42 #else
  43         long int lrintf(x)
  44         float x;
  45 #endif
  46 {
  47   __int32_t j0,sx;
  48   __uint32_t i0;
  49   float t;
  50   volatile float w;
  51   long int result;
  52
  53   GET_FLOAT_WORD(i0,x);
  54
  55   /* Extract sign bit. */
  56   sx = (i0 >> 31);
  57
  58   /* Extract exponent field. */
  59   j0 = ((i0 & 0x7f800000) >> 23) - 127;
  60
  61   if (j0 < (int)(sizeof (long int) * 8) - 1)
  62     {
  63       if (j0 >= 23)
  64         result = (long int) ((i0 & 0x7fffff) | 0x800000) << (j0 - 23);
  65       else
  66         {
  67           w = TWO23[sx] + x;
  68           t = w - TWO23[sx];
  69           GET_FLOAT_WORD (i0, t);
  70           /* Detect the all-zeros representation of plus and
  71              minus zero, which fails the calculation below. */
  72           if ((i0 & ~(1L << 31)) == 0)
  73               return 0;
  74           j0 = ((i0 >> 23) & 0xff) - 0x7f;
  75           i0 &= 0x7fffff;
  76           i0 |= 0x800000;
  77           result = (j0 < 0 ? 0 : i0 >> (23 - j0));
  78         }
  79     }
  80   else
  81     {
  82       return (long int) x;
  83     }
  84   return sx ? -result : result;
  85 }
  86
  87 #ifdef _DOUBLE_IS_32BITS
  88
  89 #ifdef __STDC__
  90         long int lrint(double x)
  91 #else
  92         long int lrint(x)
  93         double x;
  94 #endif
  95 {
  96   return lrintf((float) x);
  97 }
  98
  99 #endif /* defined(_DOUBLE_IS_32BITS) */