Subversion Repositories Kolibri OS

Rev

Rev 4872 | Blame | Compare with Previous | Last modification | View Log | RSS feed

  1.  
  2. /* @(#)s_scalbn.c 5.1 93/09/24 */
  3. /*
  4.  * ====================================================
  5.  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
  6.  *
  7.  * Developed at SunPro, a Sun Microsystems, Inc. business.
  8.  * Permission to use, copy, modify, and distribute this
  9.  * software is freely granted, provided that this notice
  10.  * is preserved.
  11.  * ====================================================
  12.  */
  13.  
  14. /*
  15. FUNCTION
  16. <<scalbn>>, <<scalbnf>>, <<scalbln>>, <<scalblnf>>--scale by power of FLT_RADIX (=2)
  17. INDEX
  18.         scalbn
  19. INDEX
  20.         scalbnf
  21. INDEX
  22.         scalbln
  23. INDEX
  24.         scalblnf
  25.  
  26. ANSI_SYNOPSIS
  27.         #include <math.h>
  28.         double scalbn(double <[x]>, int <[n]>);
  29.         float scalbnf(float <[x]>, int <[n]>);
  30.         double scalbln(double <[x]>, long int <[n]>);
  31.         float scalblnf(float <[x]>, long int <[n]>);
  32.  
  33. DESCRIPTION
  34. The <<scalbn>> and <<scalbln>> functions compute
  35.         @ifnottex
  36.         <[x]> times FLT_RADIX to the power <[n]>.
  37.         @end ifnottex
  38.         @tex
  39.         $x \cdot FLT\_RADIX^n$.
  40.         @end tex
  41. efficiently.  The result is computed by manipulating the exponent, rather than
  42. by actually performing an exponentiation or multiplication.  In this
  43. floating-point implementation FLT_RADIX=2, which makes the <<scalbn>>
  44. functions equivalent to the <<ldexp>> functions.
  45.  
  46. RETURNS
  47. <[x]> times 2 to the power <[n]>.  A range error may occur.
  48.  
  49. PORTABILITY
  50. ANSI C, POSIX
  51.  
  52. SEEALSO
  53. <<ldexp>>
  54.  
  55. */
  56.  
  57. /*
  58.  * scalbn (double x, int n)
  59.  * scalbn(x,n) returns x* 2**n  computed by  exponent  
  60.  * manipulation rather than by actually performing an
  61.  * exponentiation or a multiplication.
  62.  */
  63.  
  64. #include "fdlibm.h"
  65.  
  66. #ifndef _DOUBLE_IS_32BITS
  67.  
  68. #ifdef __STDC__
  69. static const double
  70. #else
  71. static double
  72. #endif
  73. two54   =  1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
  74. twom54  =  5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
  75. huge   = 1.0e+300,
  76. tiny   = 1.0e-300;
  77.  
  78. #ifdef __STDC__
  79.         double scalbn (double x, int n)
  80. #else
  81.         double scalbn (x,n)
  82.         double x; int n;
  83. #endif
  84. {
  85.         __int32_t  k,hx,lx;
  86.         EXTRACT_WORDS(hx,lx,x);
  87.         k = (hx&0x7ff00000)>>20;                /* extract exponent */
  88.         if (k==0) {                             /* 0 or subnormal x */
  89.             if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
  90.             x *= two54;
  91.             GET_HIGH_WORD(hx,x);
  92.             k = ((hx&0x7ff00000)>>20) - 54;
  93.             if (n< -50000) return tiny*x;       /*underflow*/
  94.             }
  95.         if (k==0x7ff) return x+x;               /* NaN or Inf */
  96.         k = k+n;
  97.         if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
  98.         if (k > 0)                              /* normal result */
  99.             {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
  100.         if (k <= -54) {
  101.             if (n > 50000)      /* in case integer overflow in n+k */
  102.                 return huge*copysign(huge,x);   /*overflow*/
  103.             else return tiny*copysign(tiny,x);  /*underflow*/
  104.       }
  105.         k += 54;                                /* subnormal result */
  106.         SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
  107.         return x*twom54;
  108. }
  109.  
  110. #endif /* _DOUBLE_IS_32BITS */
  111.