s_rint.c 7.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157
  1. /*******************************************************************************
  2. ** File: rndint.c
  3. **
  4. ** Contains: C source code for implementations of floating-point
  5. ** functions which round to integral value or format, as
  6. ** defined in header <fp.h>. In particular, this file
  7. ** contains implementations of functions rint, nearbyint,
  8. ** rinttol, round, roundtol, trunc, modf and modfl. This file
  9. ** targets PowerPC or Power platforms.
  10. **
  11. ** Written by: A. Sazegari, Apple AltiVec Group
  12. ** Created originally by Jon Okada, Apple Numerics Group
  13. **
  14. ** Copyright: © 1992-2001 by Apple Computer, Inc., all rights reserved
  15. **
  16. ** Change History (most recent first):
  17. **
  18. ** 13 Jul 01 ram replaced --setflm calls with inline assembly
  19. ** 03 Mar 01 ali first port to os x using gcc, added the crucial __setflm
  20. ** definition.
  21. ** 1. removed double_t, put in double for now.
  22. ** 2. removed iclass from nearbyint.
  23. ** 3. removed wrong comments intrunc.
  24. ** 4.
  25. ** 13 May 97 ali made performance improvements in rint, rinttol, roundtol
  26. ** and trunc by folding some of the taligent ideas into this
  27. ** implementation. nearbyint is faster than the one in taligent,
  28. ** rint is more elegant, but slower by %30 than the taligent one.
  29. ** 09 Apr 97 ali deleted modfl and deferred to AuxiliaryDD.c
  30. ** 15 Sep 94 ali Major overhaul and performance improvements of all functions.
  31. ** 20 Jul 94 PAF New faster version
  32. ** 16 Jul 93 ali Added the modfl function.
  33. ** 18 Feb 93 ali Changed the return value of fenv functions
  34. ** feclearexcept and feraiseexcept to their new
  35. ** NCEG X3J11.1/93-001 definitions.
  36. ** 16 Dec 92 JPO Removed __itrunc implementation to a
  37. ** separate file.
  38. ** 15 Dec 92 JPO Added __itrunc implementation and modified
  39. ** rinttol to include conversion from double
  40. ** to long int format. Modified roundtol to
  41. ** call __itrunc.
  42. ** 10 Dec 92 JPO Added modf (double) implementation.
  43. ** 04 Dec 92 JPO First created.
  44. **
  45. *******************************************************************************/
  46. #include <limits.h>
  47. #include <math.h>
  48. #define SET_INVALID 0x01000000UL
  49. typedef union
  50. {
  51. struct {
  52. #if defined(__BIG_ENDIAN__)
  53. unsigned long int hi;
  54. unsigned long int lo;
  55. #else
  56. unsigned long int lo;
  57. unsigned long int hi;
  58. #endif
  59. } words;
  60. double dbl;
  61. } DblInHex;
  62. static const unsigned long int signMask = 0x80000000ul;
  63. static const double twoTo52 = 4503599627370496.0;
  64. static const double doubleToLong = 4503603922337792.0; // 2^52
  65. static const DblInHex Huge = {{ 0x7FF00000, 0x00000000 }};
  66. static const DblInHex TOWARDZERO = {{ 0x00000000, 0x00000001 }};
  67. /*******************************************************************************
  68. * *
  69. * The function rint rounds its double argument to integral value *
  70. * according to the current rounding direction and returns the result in *
  71. * double format. This function signals inexact if an ordered return *
  72. * value is not equal to the operand. *
  73. * *
  74. ********************************************************************************
  75. * *
  76. * This function calls: fabs. *
  77. * *
  78. *******************************************************************************/
  79. /*******************************************************************************
  80. * First, an elegant implementation. *
  81. ********************************************************************************
  82. *
  83. *double rint ( double x )
  84. * {
  85. * double y;
  86. *
  87. * y = twoTo52.fval;
  88. *
  89. * if ( fabs ( x ) >= y ) // huge case is exact
  90. * return x;
  91. * if ( x < 0 ) y = -y; // negative case
  92. * y = ( x + y ) - y; // force rounding
  93. * if ( y == 0.0 ) // zero results mirror sign of x
  94. * y = copysign ( y, x );
  95. * return ( y );
  96. * }
  97. ********************************************************************************
  98. * Now a bit twidling version that is about %30 faster. *
  99. *******************************************************************************/
  100. double rint ( double x )
  101. {
  102. DblInHex argument;
  103. register double y;
  104. unsigned long int xHead;
  105. register long int target;
  106. argument.dbl = x;
  107. xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
  108. target = ( argument.words.hi < signMask ); // flags positive sign
  109. if ( xHead < 0x43300000ul )
  110. /*******************************************************************************
  111. * Is |x| < 2.0^52? *
  112. *******************************************************************************/
  113. {
  114. if ( xHead < 0x3ff00000ul )
  115. /*******************************************************************************
  116. * Is |x| < 1.0? *
  117. *******************************************************************************/
  118. {
  119. if ( target )
  120. y = ( x + twoTo52 ) - twoTo52; // round at binary point
  121. else
  122. y = ( x - twoTo52 ) + twoTo52; // round at binary point
  123. if ( y == 0.0 )
  124. { // fix sign of zero result
  125. if ( target )
  126. return ( 0.0 );
  127. else
  128. return ( -0.0 );
  129. }
  130. return y;
  131. }
  132. /*******************************************************************************
  133. * Is 1.0 < |x| < 2.0^52? *
  134. *******************************************************************************/
  135. if ( target )
  136. return ( ( x + twoTo52 ) - twoTo52 ); // round at binary pt.
  137. else
  138. return ( ( x - twoTo52 ) + twoTo52 );
  139. }
  140. /*******************************************************************************
  141. * |x| >= 2.0^52 or x is a NaN. *
  142. *******************************************************************************/
  143. return ( x );
  144. }