s_ceil.c 4.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113
  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. #include <math.h>
  24. #include <endian.h>
  25. static const double twoTo52 = 4503599627370496.0;
  26. static const unsigned long signMask = 0x80000000ul;
  27. typedef union
  28. {
  29. struct {
  30. #if (__BYTE_ORDER == __BIG_ENDIAN)
  31. unsigned long int hi;
  32. unsigned long int lo;
  33. #else
  34. unsigned long int lo;
  35. unsigned long int hi;
  36. #endif
  37. } words;
  38. double dbl;
  39. } DblInHex;
  40. /*******************************************************************************
  41. * Functions needed for the computation. *
  42. *******************************************************************************/
  43. /*******************************************************************************
  44. * Ceil(x) returns the smallest integer not less than x. *
  45. *******************************************************************************/
  46. libm_hidden_proto(ceil)
  47. double ceil ( double x )
  48. {
  49. DblInHex xInHex,OldEnvironment;
  50. register double y;
  51. register unsigned long int xhi;
  52. register int target;
  53. xInHex.dbl = x;
  54. xhi = xInHex.words.hi & 0x7fffffffUL; // xhi is the high half of |x|
  55. target = ( xInHex.words.hi < signMask );
  56. if ( xhi < 0x43300000ul )
  57. /*******************************************************************************
  58. * Is |x| < 2.0^52? *
  59. *******************************************************************************/
  60. {
  61. if ( xhi < 0x3ff00000ul )
  62. /*******************************************************************************
  63. * Is |x| < 1.0? *
  64. *******************************************************************************/
  65. {
  66. if ( ( xhi | xInHex.words.lo ) == 0ul ) // zero x is exact case
  67. return ( x );
  68. else
  69. { // inexact case
  70. asm ("mffs %0" : "=f" (OldEnvironment.dbl));
  71. OldEnvironment.words.lo |= 0x02000000ul;
  72. asm ("mtfsf 255,%0" : /*NULLOUT*/ : /*IN*/ "f" ( OldEnvironment.dbl ));
  73. if ( target )
  74. return ( 1.0 );
  75. else
  76. return ( -0.0 );
  77. }
  78. }
  79. /*******************************************************************************
  80. * Is 1.0 < |x| < 2.0^52? *
  81. *******************************************************************************/
  82. if ( target )
  83. {
  84. y = ( x + twoTo52 ) - twoTo52; // round at binary pt.
  85. if ( y < x )
  86. return ( y + 1.0 );
  87. else
  88. return ( y );
  89. }
  90. else
  91. {
  92. y = ( x - twoTo52 ) + twoTo52; // round at binary pt.
  93. if ( y < x )
  94. return ( y + 1.0 );
  95. else
  96. return ( y );
  97. }
  98. }
  99. /*******************************************************************************
  100. * |x| >= 2.0^52 or x is a NaN. *
  101. *******************************************************************************/
  102. return ( x );
  103. }
  104. libm_hidden_def(ceil)