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 #if defined(LIBM_SCCS) && !defined(lint)
  19 static char rcsid[] = "$NetBSD: $";
  20 #endif
  21
  22 /*
  23  * rintl(x)
  24  * Return x rounded to integral value according to the prevailing
  25  * rounding mode.
  26  * Method:
  27  *      Using floating addition.
  28  * Exception:
  29  *      Inexact flag raised if x not equal to rintl(x).
  30  */
  31
  32 #include <fenv.h>
  33 #include "math.h"
  34 #include "math_private.h"
  35
  36 #ifdef __STDC__
  37 static const long double
  38 #else
  39 static long double
  40 #endif
  41 TWO63[2]={
  42   9.223372036854775808000000e+18, /* 0x403E, 0x00000000, 0x00000000 */
  43  -9.223372036854775808000000e+18  /* 0xC03E, 0x00000000, 0x00000000 */
  44 };
  45
  46 #ifdef __STDC__
  47         long double __rintl(long double x)
  48 #else
  49         long double __rintl(x)
  50         long double x;
  51 #endif
  52 {
  53         fenv_t env;
  54         int32_t se,j0,sx;
  55         u_int32_t i,i0,i1;
  56         long double w,t;
  57         GET_LDOUBLE_WORDS(se,i0,i1,x);
  58         sx = (se>>15)&1;
  59         j0 = (se&0x7fff)-0x3fff;
  60         if(j0<31) {
  61             if(j0<0) {
  62                 if(((se&0x7fff)|i0|i1)==0) return x;
  63                 i1 |= i0;
  64                 i0 &= 0xe0000000;
  65                 i0 |= (i1|-i1)&0x80000000;
  66                 SET_LDOUBLE_MSW(x,i0);
  67                 feholdexcept (&env);
  68                 w = TWO63[sx]+x;
  69                 t = w-TWO63[sx];
  70                 fesetenv (&env);
  71                 GET_LDOUBLE_EXP(i0,t);
  72                 SET_LDOUBLE_EXP(t,(i0&0x7fff)|(sx<<15));
  73                 return t;
  74             } else {
  75                 i = (0x7fffffff)>>j0;
  76                 if(((i0&i)|i1)==0) return x; /* x is integral */
  77                 i>>=1;
  78                 if(((i0&i)|i1)!=0) {
  79                     if(j0==31) i1 = 0x40000000; else
  80                     i0 = (i0&(~i))|((0x20000000)>>j0);
  81                     /* Shouldn't this be
  82                          if (j0 >= 30) i1 = 0x80000000 >> (j0 - 30);
  83                          i0 = (i0&(~i))|((0x20000000)>>j0);
  84                        If yes, this should be correct in s_rint and
  85                        s_rintf, too.  -- drepper@cygnus.com */
  86                 }
  87             }
  88         } else if (j0>62) {
  89             if(j0==0x4000) return x+x;  /* inf or NaN */
  90             else return x;              /* x is integral */
  91         } else {
  92             i = ((u_int32_t)(0xffffffff))>>(j0-31);
  93             if((i1&i)==0) return x;     /* x is integral */
  94             i>>=1;
  95             if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-31));
  96         }
  97         SET_LDOUBLE_WORDS(x,se,i0,i1);
  98         feholdexcept (&env);
  99         w = TWO63[sx]+x;
 100         t = w-TWO63[sx];
 101         fesetenv (&env);
 102         return t;
 103 }
 104 weak_alias (__rintl, rintl)