1 | /* |
2 | * ==================================================== |
3 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
4 | * |
5 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
6 | * Permission to use, copy, modify, and distribute this |
7 | * software is freely granted, provided that this notice |
8 | * is preserved. |
9 | * ==================================================== |
10 | */ |
11 | |
12 | /* |
13 | * scalbln(double x, long n) |
14 | * scalbln(x,n) returns x * 2**n computed by exponent |
15 | * manipulation rather than by actually performing an |
16 | * exponentiation or a multiplication. |
17 | */ |
18 | |
19 | #include "math_libm.h" |
20 | #include "math_private.h" |
21 | #include <limits.h> |
22 | |
23 | #ifdef __WATCOMC__ /* Watcom defines huge=__huge */ |
24 | #undef huge |
25 | #endif |
26 | |
27 | static const double |
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 | |
33 | double scalbln(double x, long n) |
34 | { |
35 | int32_t k, hx, lx; |
36 | |
37 | EXTRACT_WORDS(hx, lx, x); |
38 | k = (hx & 0x7ff00000) >> 20; /* extract exponent */ |
39 | if (k == 0) { /* 0 or subnormal x */ |
40 | if ((lx | (hx & 0x7fffffff)) == 0) |
41 | return x; /* +-0 */ |
42 | x *= two54; |
43 | GET_HIGH_WORD(hx, x); |
44 | k = ((hx & 0x7ff00000) >> 20) - 54; |
45 | } |
46 | if (k == 0x7ff) |
47 | return x + x; /* NaN or Inf */ |
48 | k = (int32_t)(k + n); |
49 | if (k > 0x7fe) |
50 | return huge * copysign(huge, x); /* overflow */ |
51 | if (n < -50000) |
52 | return tiny * copysign(tiny, x); /* underflow */ |
53 | if (k > 0) { /* normal result */ |
54 | SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20)); |
55 | return x; |
56 | } |
57 | if (k <= -54) { |
58 | if (n > 50000) /* in case integer overflow in n+k */ |
59 | return huge * copysign(huge, x); /* overflow */ |
60 | return tiny * copysign(tiny, x); /* underflow */ |
61 | } |
62 | k += 54; /* subnormal result */ |
63 | SET_HIGH_WORD(x, (hx & 0x800fffff) | (k << 20)); |
64 | return x * twom54; |
65 | } |
66 | libm_hidden_def(scalbln) |
67 | |
68 | |
69 | double scalbn(double x, int n) |
70 | { |
71 | return scalbln(x, n); |
72 | } |
73 | libm_hidden_def(scalbn) |
74 | |