sysdeps/libm-ieee754/s_nearbyintl.c

   1 /* s_rintl.c -- long double version of s_rint.c.
   2  * Conversion to long double by Ulrich Drepper,
   3  * Cygnus Support, drepper@cygnus.com.
   4  */
   5 /* Adapted for use as nearbyint by Ulrich Drepper <drepper@cygnus.com>.  */
   6
   7 /*
   8  * ====================================================
   9  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  10  *
  11  * Developed at SunPro, a Sun Microsystems, Inc. business.
  12  * Permission to use, copy, modify, and distribute this
  13  * software is freely granted, provided that this notice
  14  * is preserved.
  15  * ====================================================
  16  */
  17
  18 /*
  19  * rintl(x)
  20  * Return x rounded to integral value according to the prevailing
  21  * rounding mode.
  22  * Method:
  23  *      Using floating addition.
  24  * Exception:
  25  *      Inexact flag raised if x not equal to rintl(x).
  26  */
  27
  28 #include <fenv.h>
  29 #include "math.h"
  30 #include "math_private.h"
  31
  32 #ifdef __STDC__
  33 static const long double
  34 #else
  35 static long double
  36 #endif
  37 TWO63[2]={
  38   9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */
  39  -9.223372036854775808000000e+18  /* 0xC03E, 0x00000000, 0x00000000 */
  40 };
  41
  42 #ifdef __STDC__
  43         long double __nearbyintl(long double x)
  44 #else
  45         long double __nearbyintl(x)
  46         long double x;
  47 #endif
  48 {
  49         fenv_t env;
  50         int32_t se,j0,sx;
  51         u_int32_t i,i0,i1;
  52         long double w,t;
  53         GET_LDOUBLE_WORDS(se,i0,i1,x);
  54         sx = (se>>15)&1;
  55         j0 = (se&0x7fff)-0x3fff;
  56         if(j0<31) {
  57             if(j0<0) {
  58                 if(((se&0x7fff)|i0|i1)==0) return x;
  59                 i1 |= i0;
  60                 i0 &= 0xe0000000;
  61                 i0 |= (i1|-i1)&0x80000000;
  62                 SET_LDOUBLE_MSW(x,i0);
  63                 feholdexcept (&env);
  64                 w = TWO63[sx]+x;
  65                 t = w-TWO63[sx];
  66                 fesetenv (&env);
  67                 GET_LDOUBLE_EXP(i0,t);
  68                 SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15));
  69                 return t;
  70             } else {
  71                 i = (0x7fffffff)>>j0;
  72                 if(((i0&i)|i1)==0) return x; /* x is integral */
  73                 i>>=1;
  74                 if(((i0&i)|i1)!=0) {
  75                     if(j0==31) i1 = 0x40000000; else
  76                     i0 = (i0&(~i))|((0x20000000)>>j0);
  77                     /* Shouldn't this be
  78                          if (j0 >= 30) i1 = 0x80000000 >> (j0 - 30);
  79                          i0 = (i0&(~i))|((0x20000000)>>j0);
  80                        If yes, this should be correct in s_rint and
  81                        s_rintf, too.  -- drepper@cygnus.com */
  82                 }
  83             }
  84         } else if (j0>62) {
  85             if(j0==0x4000) return x+x;  /* inf or NaN */
  86             else return x;              /* x is integral */
  87         } else {
  88             i = ((u_int32_t)(0xffffffff))>>(j0-31);
  89             if((i1&i)==0) return x;     /* x is integral */
  90             i>>=1;
  91             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31));
  92         }
  93         SET_LDOUBLE_WORDS(x,se,i0,i1);
  94         feholdexcept (&env);
  95         w = TWO63[sx]+x;
  96         t = w-TWO63[sx];
  97         fesetenv (&env);
  98         return t;
  99 }
 100 weak_alias (__nearbyintl, nearbyintl)