scalb.c 3.3 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485
  1. /***********************************************************************
  2. ** File: scalb.c
  3. **
  4. ** Contains: C source code for implementations of floating-point
  5. ** scalb functions defined in header <fp.h>. In
  6. ** particular, this file contains implementations of
  7. ** functions scalb and scalbl for double and long double
  8. ** formats on PowerPC platforms.
  9. **
  10. ** Written by: Jon Okada, SANEitation Engineer, ext. 4-4838
  11. **
  12. ** Copyright: © 1992 by Apple Computer, Inc., all rights reserved
  13. **
  14. ** Change History ( most recent first ):
  15. **
  16. ** 28 May 97 ali made an speed improvement for large n,
  17. ** removed scalbl.
  18. ** 12 Dec 92 JPO First created.
  19. **
  20. ***********************************************************************/
  21. typedef union
  22. {
  23. struct {
  24. #if defined(__BIG_ENDIAN__)
  25. unsigned long int hi;
  26. unsigned long int lo;
  27. #else
  28. unsigned long int lo;
  29. unsigned long int hi;
  30. #endif
  31. } words;
  32. double dbl;
  33. } DblInHex;
  34. static const double twoTo1023 = 8.988465674311579539e307; // 0x1p1023
  35. static const double twoToM1022 = 2.225073858507201383e-308; // 0x1p-1022
  36. /***********************************************************************
  37. double scalb( double x, long int n ) returns its argument x scaled
  38. by the factor 2^m. NaNs, signed zeros, and infinities are propagated
  39. by this function regardless of the value of n.
  40. Exceptions: OVERFLOW/INEXACT or UNDERFLOW inexact may occur;
  41. INVALID for signaling NaN inputs ( quiet NaN returned ).
  42. Calls: none.
  43. ***********************************************************************/
  44. double scalb ( double x, int n )
  45. {
  46. DblInHex xInHex;
  47. xInHex.words.lo = 0UL; // init. low half of xInHex
  48. if ( n > 1023 )
  49. { // large positive scaling
  50. if ( n > 2097 ) // huge scaling
  51. return ( ( x * twoTo1023 ) * twoTo1023 ) * twoTo1023;
  52. while ( n > 1023 )
  53. { // scale reduction loop
  54. x *= twoTo1023; // scale x by 2^1023
  55. n -= 1023; // reduce n by 1023
  56. }
  57. }
  58. else if ( n < -1022 )
  59. { // large negative scaling
  60. if ( n < -2098 ) // huge negative scaling
  61. return ( ( x * twoToM1022 ) * twoToM1022 ) * twoToM1022;
  62. while ( n < -1022 )
  63. { // scale reduction loop
  64. x *= twoToM1022; // scale x by 2^( -1022 )
  65. n += 1022; // incr n by 1022
  66. }
  67. }
  68. /*******************************************************************************
  69. * -1022 <= n <= 1023; convert n to double scale factor. *
  70. *******************************************************************************/
  71. xInHex.words.hi = ( ( unsigned long ) ( n + 1023 ) ) << 20;
  72. return ( x * xInHex.dbl );
  73. }