s_ceil.c 4.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109
  1. /*******************************************************************************
  2. * *
  3. * File ceilfloor.c, *
  4. * Function ceil(x) and floor(x), *
  5. * Implementation of ceil and floor for the PowerPC. *
  6. * *
  7. * Copyright © 1991 Apple Computer, Inc. All rights reserved. *
  8. * *
  9. * Written by Ali Sazegari, started on November 1991, *
  10. * *
  11. * based on math.h, library code for Macintoshes with a 68881/68882 *
  12. * by Jim Thomas. *
  13. * *
  14. * W A R N I N G: This routine expects a 64 bit double model. *
  15. * *
  16. * December 03 1992: first rs6000 port. *
  17. * July 14 1993: comment changes and addition of #pragma fenv_access. *
  18. * May 06 1997: port of the ibm/taligent ceil and floor routines. *
  19. * April 11 2001: first port to os x using gcc. *
  20. * June 13 2001: replaced __setflm with in-line assembly *
  21. * *
  22. *******************************************************************************/
  23. static const double twoTo52 = 4503599627370496.0;
  24. static const unsigned long signMask = 0x80000000ul;
  25. typedef union
  26. {
  27. struct {
  28. #if defined(__BIG_ENDIAN__)
  29. unsigned long int hi;
  30. unsigned long int lo;
  31. #else
  32. unsigned long int lo;
  33. unsigned long int hi;
  34. #endif
  35. } words;
  36. double dbl;
  37. } DblInHex;
  38. /*******************************************************************************
  39. * Functions needed for the computation. *
  40. *******************************************************************************/
  41. /*******************************************************************************
  42. * Ceil(x) returns the smallest integer not less than x. *
  43. *******************************************************************************/
  44. double ceil ( double x )
  45. {
  46. DblInHex xInHex,OldEnvironment;
  47. register double y;
  48. register unsigned long int xhi;
  49. register int target;
  50. xInHex.dbl = x;
  51. xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x|
  52. target = ( xInHex.words.hi < signMask );
  53. if ( xhi < 0x43300000ul )
  54. /*******************************************************************************
  55. * Is |x| < 2.0^52? *
  56. *******************************************************************************/
  57. {
  58. if ( xhi < 0x3ff00000ul )
  59. /*******************************************************************************
  60. * Is |x| < 1.0? *
  61. *******************************************************************************/
  62. {
  63. if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case
  64. return ( x );
  65. else
  66. { // inexact case
  67. asm ("mffs %0" : "=f" (OldEnvironment.dbl));
  68. OldEnvironment.words.lo |= 0x02000000ul;
  69. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  70. if ( target )
  71. return ( 1.0 );
  72. else
  73. return ( -0.0 );
  74. }
  75. }
  76. /*******************************************************************************
  77. * Is 1.0 < |x| < 2.0^52? *
  78. *******************************************************************************/
  79. if ( target )
  80. {
  81. y = ( x + twoTo52 ) - twoTo52; // round at binary pt.
  82. if ( y < x )
  83. return ( y + 1.0 );
  84. else
  85. return ( y );
  86. }
  87. else
  88. {
  89. y = ( x - twoTo52 ) + twoTo52; // round at binary pt.
  90. if ( y < x )
  91. return ( y + 1.0 );
  92. else
  93. return ( y );
  94. }
  95. }
  96. /*******************************************************************************
  97. * |x| >= 2.0^52 or x is a NaN. *
  98. *******************************************************************************/
  99. return ( x );
  100. }