s_scalbn.c 2.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869
  1. /* @(#)s_scalbn.c 5.1 93/09/24 */
  2. /*
  3. * ====================================================
  4. * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  5. *
  6. * Developed at SunPro, a Sun Microsystems, Inc. business.
  7. * Permission to use, copy, modify, and distribute this
  8. * software is freely granted, provided that this notice
  9. * is preserved.
  10. * ====================================================
  11. */
  12. #if defined(LIBM_SCCS) && !defined(lint)
  13. static char rcsid[] = "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
  14. #endif
  15. /*
  16. * scalbn (double x, int n)
  17. * scalbn(x,n) returns x* 2**n computed by exponent
  18. * manipulation rather than by actually performing an
  19. * exponentiation or a multiplication.
  20. */
  21. #include "math.h"
  22. #include "math_private.h"
  23. #ifdef __STDC__
  24. static const double
  25. #else
  26. static double
  27. #endif
  28. two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  29. twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
  30. huge = 1.0e+300,
  31. tiny = 1.0e-300;
  32. #ifdef __STDC__
  33. double scalbn (double x, int n)
  34. #else
  35. double scalbn (x,n)
  36. double x; int n;
  37. #endif
  38. {
  39. int32_t k,hx,lx;
  40. EXTRACT_WORDS(hx,lx,x);
  41. k = (hx&0x7ff00000)>>20; /* extract exponent */
  42. if (k==0) { /* 0 or subnormal x */
  43. if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
  44. x *= two54;
  45. GET_HIGH_WORD(hx,x);
  46. k = ((hx&0x7ff00000)>>20) - 54;
  47. if (n< -50000) return tiny*x; /*underflow*/
  48. }
  49. if (k==0x7ff) return x+x; /* NaN or Inf */
  50. k = k+n;
  51. if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */
  52. if (k > 0) /* normal result */
  53. {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
  54. if (k <= -54) {
  55. if (n > 50000) /* in case integer overflow in n+k */
  56. return huge*copysign(huge,x); /*overflow*/
  57. else return tiny*copysign(tiny,x); /*underflow*/
  58. }
  59. k += 54; /* subnormal result */
  60. SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
  61. return x*twom54;
  62. }
  63. libm_hidden_def(scalbn)