s_logb.c 5.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103
  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. typedef union
  35. {
  36. struct {
  37. #if defined(__BIG_ENDIAN__)
  38. unsigned long int hi;
  39. unsigned long int lo;
  40. #else
  41. unsigned long int lo;
  42. unsigned long int hi;
  43. #endif
  44. } words;
  45. double dbl;
  46. } DblInHex;
  47. static const double twoTo52 = 4.50359962737049600e15; // 0x1p52
  48. static const double klTod = 4503601774854144.0; // 0x1.000008p52
  49. static const unsigned long int signMask = 0x80000000ul;
  50. static const DblInHex minusInf = {{ 0xFFF00000, 0x00000000 }};
  51. /*******************************************************************************
  52. ********************************************************************************
  53. * L O G B *
  54. ********************************************************************************
  55. *******************************************************************************/
  56. double logb ( double x )
  57. {
  58. DblInHex xInHex;
  59. long int shiftedExp;
  60. xInHex.dbl = x;
  61. shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
  62. if ( shiftedExp == 2047 )
  63. { // NaN or INF
  64. if ( ( ( xInHex.words.hi & signMask ) == 0 ) || ( x != x ) )
  65. return x; // NaN or +INF return x
  66. else
  67. return -x; // -INF returns +INF
  68. }
  69. if ( shiftedExp != 0 ) // normal number
  70. shiftedExp -= 1023; // unbias exponent
  71. else if ( x == 0.0 )
  72. { // zero
  73. xInHex.words.hi = 0x0UL; // return -infinity
  74. return ( minusInf.dbl );
  75. }
  76. else
  77. { // subnormal number
  78. xInHex.dbl *= twoTo52; // scale up
  79. shiftedExp = ( xInHex.words.hi & 0x7ff00000UL ) >> 20;
  80. shiftedExp -= 1075; // unbias exponent
  81. }
  82. if ( shiftedExp == 0 ) // zero result
  83. return ( 0.0 );
  84. else
  85. { // nonzero result
  86. xInHex.dbl = klTod;
  87. xInHex.words.lo += shiftedExp;
  88. return ( xInHex.dbl - klTod );
  89. }
  90. }