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 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 #ifdef _NEED_FLOAT64
67
68 static const __float64
69 two54 = _F_64(1.80143985094819840000e+16), /* 0x43500000, 0x00000000 */
70 twom54 = _F_64(5.55111512312578270212e-17); /* 0x3C900000, 0x00000000 */
71
72 __float64
scalbn64(__float64 x,int n)73 scalbn64(__float64 x, int n)
74 {
75 __int32_t k,hx,lx;
76 EXTRACT_WORDS(hx,lx,x);
77 k = (hx&0x7ff00000)>>20; /* extract exponent */
78 if (k==0) { /* 0 or subnormal x */
79 if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */
80 x *= two54;
81 GET_HIGH_WORD(hx,x);
82 k = ((hx&0x7ff00000)>>20) - 54;
83 #if __SIZEOF_INT__ > 2
84 if (n< -50000) return __math_uflow(hx<0); /*underflow*/
85 #endif
86 }
87 if (k==0x7ff) return x+x; /* NaN or Inf */
88 #if __SIZEOF_INT__ > 2
89 if (n > 50000) /* in case integer overflow in n+k */
90 return __math_oflow(hx<0); /*overflow*/
91 #endif
92 k = k+n;
93 if (k > 0x7fe) return __math_oflow(hx<0); /* overflow */
94 if (k > 0) /* normal result */
95 {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
96 if (k <= -54)
97 return __math_uflow(hx<0); /*underflow*/
98 k += 54; /* subnormal result */
99 SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
100 return check_uflow(x*twom54);
101 }
102
103 #if defined(_HAVE_ALIAS_ATTRIBUTE)
104 #ifndef __clang__
105 #pragma GCC diagnostic ignored "-Wmissing-attributes"
106 #endif
107 __strong_reference(scalbn64, ldexp64);
108 #else
109
110 __float64
ldexp64(__float64 value,int exp)111 ldexp64(__float64 value, int exp)
112 {
113 return scalbn64(value, exp);
114 }
115
116 #endif
117
118 _MATH_ALIAS_d_di(scalbn)
119 _MATH_ALIAS_d_di(ldexp)
120
121 #endif /* _NEED_FLOAT64 */
122
