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
13 #if defined(LIBM_SCCS) && !defined(lint)
14 static const char rcsid[] =
15 "$NetBSD: s_scalbn.c,v 1.8 1995/05/10 20:48:08 jtc Exp $";
16 #endif
17
18 /*
19 * scalbn (double x, int n)
20 * scalbn(x,n) returns x* 2**n computed by exponent
21 * manipulation rather than by actually performing an
22 * exponentiation or a multiplication.
23 */
24
25 #include "math_libm.h"
26 #include "math_private.h"
27
28 libm_hidden_proto(copysign)
29 #ifdef __STDC__
30 static const double
31 #else
32 static double
33 #endif
34 two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */
35 twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */
36 huge_val = 1.0e+300, tiny = 1.0e-300;
37
libm_hidden_proto(scalbn)38 libm_hidden_proto(scalbn)
39 #ifdef __STDC__
40 double scalbn(double x, int n)
41 #else
42 double scalbn(x, n)
43 double x;
44 int n;
45 #endif
46 {
47 int32_t k, hx, lx;
48 EXTRACT_WORDS(hx, lx, x);
49 k = (hx & 0x7ff00000) >> 20; /* extract exponent */
50 if (k == 0) { /* 0 or subnormal x */
51 if ((lx | (hx & 0x7fffffff)) == 0)
52 return x; /* +-0 */
53 x *= two54;
54 GET_HIGH_WORD(hx, x);
55 k = ((hx & 0x7ff00000) >> 20) - 54;
56 if (n < -50000)
57 return tiny * x; /* underflow */
58 }
59 if (k == 0x7ff)
60 return x + x; /* NaN or Inf */
61 k = k + n;
62 if (k > 0x7fe)
63 return huge_val * copysign(huge_val, x); /* overflow */
64 if (k > 0) { /* normal result */
65 SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
66 return x;
67 }
68 if (k <= -54) {
69 if (n > 50000) /* in case integer overflow in n+k */
70 return huge_val * copysign(huge_val, x); /* overflow */
71 else
72 return tiny * copysign(tiny, x); /* underflow */
73 }
74 k += 54; /* subnormal result */
75 SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20));
76 return x * twom54;
77 }
78
79 libm_hidden_def(scalbn)
80