s_rint.c 2.1 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788
  1. #if !defined(__ppc__)
  2. /* @(#)s_rint.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. #if defined(LIBM_SCCS) && !defined(lint)
  14. static char rcsid[] = "$NetBSD: s_rint.c,v 1.8 1995/05/10 20:48:04 jtc Exp $";
  15. #endif
  16. /*
  17. * rint(x)
  18. * Return x rounded to integral value according to the prevailing
  19. * rounding mode.
  20. * Method:
  21. * Using floating addition.
  22. * Exception:
  23. * Inexact flag raised if x not equal to rint(x).
  24. */
  25. #include "math.h"
  26. #include "math_private.h"
  27. #ifdef __STDC__
  28. static const double
  29. #else
  30. static double
  31. #endif
  32. TWO52[2]={
  33. 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */
  34. -4.50359962737049600000e+15, /* 0xC3300000, 0x00000000 */
  35. };
  36. #ifdef __STDC__
  37. double rint(double x)
  38. #else
  39. double rint(x)
  40. double x;
  41. #endif
  42. {
  43. int32_t i0,j0,sx;
  44. u_int32_t i,i1;
  45. double w,t;
  46. EXTRACT_WORDS(i0,i1,x);
  47. sx = (i0>>31)&1;
  48. j0 = ((i0>>20)&0x7ff)-0x3ff;
  49. if(j0<20) {
  50. if(j0<0) {
  51. if(((i0&0x7fffffff)|i1)==0) return x;
  52. i1 |= (i0&0x0fffff);
  53. i0 &= 0xfffe0000;
  54. i0 |= ((i1|-i1)>>12)&0x80000;
  55. SET_HIGH_WORD(x,i0);
  56. w = TWO52[sx]+x;
  57. t = w-TWO52[sx];
  58. GET_HIGH_WORD(i0,t);
  59. SET_HIGH_WORD(t,(i0&0x7fffffff)|(sx<<31));
  60. return t;
  61. } else {
  62. i = (0x000fffff)>>j0;
  63. if(((i0&i)|i1)==0) return x; /* x is integral */
  64. i>>=1;
  65. if(((i0&i)|i1)!=0) {
  66. if(j0==19) i1 = 0x40000000; else
  67. i0 = (i0&(~i))|((0x20000)>>j0);
  68. }
  69. }
  70. } else if (j0>51) {
  71. if(j0==0x400) return x+x; /* inf or NaN */
  72. else return x; /* x is integral */
  73. } else {
  74. i = ((u_int32_t)(0xffffffff))>>(j0-20);
  75. if((i1&i)==0) return x; /* x is integral */
  76. i>>=1;
  77. if((i1&i)!=0) i1 = (i1&(~i))|((0x40000000)>>(j0-20));
  78. }
  79. INSERT_WORDS(x,i0,i1);
  80. w = TWO52[sx]+x;
  81. return w-TWO52[sx];
  82. }
  83. #endif /* !__ppc__ */