s_logb.c 5.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107
  1. /*******************************************************************************
  2. * *
  3. * File logb.c, *
  4. * Functions logb. *
  5. * Implementation of logb for the PowerPC. *
  6. * *
  7. * Copyright © 1991 Apple Computer, Inc. All rights reserved. *
  8. * *
  9. * Written by Ali Sazegari, started on June 1991, *
  10. * *
  11. * August 26 1991: removed CFront Version 1.1d17 warnings. *
  12. * August 27 1991: no errors reported by the test suite. *
  13. * November 11 1991: changed CLASSEXTENDED to the macro CLASSIFY and *
  14. * + or - infinity to constants. *
  15. * November 18 1991: changed the macro CLASSIFY to CLASSEXTENDEDint to *
  16. * improve performance. *
  17. * February 07 1992: changed bit operations to macros ( object size is *
  18. * unchanged ). *
  19. * September24 1992: took the "#include support.h" out. *
  20. * December 03 1992: first rs/6000 port. *
  21. * August 30 1992: set the divide by zero for the zero argument case. *
  22. * October 05 1993: corrected the environment. *
  23. * October 17 1994: replaced all environmental functions with __setflm. *
  24. * May 28 1997: made speed improvements. *
  25. * April 30 2001: forst mac os x port using gcc. *
  26. * *
  27. ********************************************************************************
  28. * The C math library offers a similar function called "frexp". It is *
  29. * different in details from logb, but similar in spirit. This current *
  30. * implementation of logb follows the recommendation in IEEE Standard 854 *
  31. * which is different in its handling of denormalized numbers from the IEEE *
  32. * Standard 754. *
  33. *******************************************************************************/
  34. #include <math.h>
  35. #include <endian.h>
  36. typedef union
  37. {
  38. struct {
  39. #if (__BYTE_ORDER == __BIG_ENDIAN)
  40. unsigned long int hi;
  41. unsigned long int lo;
  42. #else
  43. unsigned long int lo;
  44. unsigned long int hi;
  45. #endif
  46. } words;
  47. double dbl;
  48. } DblInHex;
  49. static const double twoTo52 = 4.50359962737049600e15; // 0x1p52
  50. static const double klTod = 4503601774854144.0; // 0x1.000008p52
  51. static const unsigned long int signMask = 0x80000000ul;
  52. static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }};
  53. /*******************************************************************************
  54. ********************************************************************************
  55. * L O G B *
  56. ********************************************************************************
  57. *******************************************************************************/
  58. libm_hidden_proto(logb)
  59. double logb ( double x )
  60. {
  61. DblInHex xInHex;
  62. long int shiftedExp;
  63. xInHex.dbl = x;
  64. shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
  65. if ( shiftedExp == 2047 )
  66. { // NaN or INF
  67. if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) )
  68. return x; // NaN or +INF return x
  69. else
  70. return -x; // -INF returns +INF
  71. }
  72. if ( shiftedExp != 0 ) // normal number
  73. shiftedExp -= 1023; // unbias exponent
  74. else if ( x == 0.0 )
  75. { // zero
  76. xInHex.words.hi = 0x0UL; // return -infinity
  77. return ( minusInf.dbl );
  78. }
  79. else
  80. { // subnormal number
  81. xInHex.dbl *= twoTo52; // scale up
  82. shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
  83. shiftedExp -= 1075; // unbias exponent
  84. }
  85. if ( shiftedExp == 0 ) // zero result
  86. return ( 0.0 );
  87. else
  88. { // nonzero result
  89. xInHex.dbl = klTod;
  90. xInHex.words.lo += shiftedExp;
  91. return ( xInHex.dbl - klTod );
  92. }
  93. }
  94. libm_hidden_def(logb)