s_round.c 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115
  1. #include <limits.h>
  2. #include <math.h>
  3. #include <endian.h>
  4. typedef union
  5. {
  6. struct {
  7. #if (__BYTE_ORDER == __BIG_ENDIAN)
  8. unsigned long int hi;
  9. unsigned long int lo;
  10. #else
  11. unsigned long int lo;
  12. unsigned long int hi;
  13. #endif
  14. } words;
  15. double dbl;
  16. } DblInHex;
  17. static const unsigned long int signMask = 0x80000000ul;
  18. static const double twoTo52 = 4503599627370496.0;
  19. /*******************************************************************************
  20. * *
  21. * The function round rounds its double argument to integral value *
  22. * according to the "add half to the magnitude and truncate" rounding of *
  23. * Pascal's Round function and FORTRAN's ANINT function and returns the *
  24. * result in double format. This function signals inexact if an ordered *
  25. * return value is not equal to the operand. *
  26. * *
  27. *******************************************************************************/
  28. libm_hidden_proto(round)
  29. double round ( double x )
  30. {
  31. DblInHex argument, OldEnvironment;
  32. register double y, z;
  33. register unsigned long int xHead;
  34. register long int target;
  35. argument.dbl = x;
  36. xHead = argument.words.hi & 0x7fffffffUL; // xHead <- high half of |x|
  37. target = ( argument.words.hi < signMask ); // flag positive sign
  38. if ( xHead < 0x43300000ul )
  39. /*******************************************************************************
  40. * Is |x| < 2.0^52? *
  41. *******************************************************************************/
  42. {
  43. if ( xHead < 0x3ff00000ul )
  44. /*******************************************************************************
  45. * Is |x| < 1.0? *
  46. *******************************************************************************/
  47. {
  48. asm ("mffs %0" : "=f" (OldEnvironment.dbl)); // get environment
  49. if ( xHead < 0x3fe00000ul )
  50. /*******************************************************************************
  51. * Is |x| < 0.5? *
  52. *******************************************************************************/
  53. {
  54. if ( ( xHead | argument.words.lo ) != 0ul )
  55. OldEnvironment.words.lo |= 0x02000000ul;
  56. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  57. if ( target )
  58. return ( 0.0 );
  59. else
  60. return ( -0.0 );
  61. }
  62. /*******************************************************************************
  63. * Is 0.5 ² |x| < 1.0? *
  64. *******************************************************************************/
  65. OldEnvironment.words.lo |= 0x02000000ul;
  66. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  67. if ( target )
  68. return ( 1.0 );
  69. else
  70. return ( -1.0 );
  71. }
  72. /*******************************************************************************
  73. * Is 1.0 < |x| < 2.0^52? *
  74. *******************************************************************************/
  75. if ( target )
  76. { // positive x
  77. y = ( x + twoTo52 ) - twoTo52; // round at binary point
  78. if ( y == x ) // exact case
  79. return ( x );
  80. z = x + 0.5; // inexact case
  81. y = ( z + twoTo52 ) - twoTo52; // round at binary point
  82. if ( y > z )
  83. return ( y - 1.0 );
  84. else
  85. return ( y );
  86. }
  87. /*******************************************************************************
  88. * Is x < 0? *
  89. *******************************************************************************/
  90. else
  91. {
  92. y = ( x - twoTo52 ) + twoTo52; // round at binary point
  93. if ( y == x )
  94. return ( x );
  95. z = x - 0.5;
  96. y = ( z - twoTo52 ) + twoTo52; // round at binary point
  97. if ( y < z )
  98. return ( y + 1.0 );
  99. else
  100. return ( y );
  101. }
  102. }
  103. /*******************************************************************************
  104. * |x| >= 2.0^52 or x is a NaN. *
  105. *******************************************************************************/
  106. return ( x );
  107. }
  108. libm_hidden_def(round)