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 | * scalbn (double x, int n) |
14 | * scalbn(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.h> |
20 | #include <math_private.h> |
21 | |
22 | static const double |
23 | two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ |
24 | twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ |
25 | huge = 1.0e+300, |
26 | tiny = 1.0e-300; |
27 | |
28 | double |
29 | __scalbn (double x, int n) |
30 | { |
31 | int64_t ix; |
32 | int64_t k; |
33 | EXTRACT_WORDS64(ix,x); |
34 | k = (ix >> 52) & 0x7ff; /* extract exponent */ |
35 | if (__builtin_expect(k==0, 0)) { /* 0 or subnormal x */ |
36 | if ((ix & UINT64_C(0xfffffffffffff))==0) return x; /* +-0 */ |
37 | x *= two54; |
38 | EXTRACT_WORDS64(ix,x); |
39 | k = ((ix >> 52) & 0x7ff) - 54; |
40 | } |
41 | if (__builtin_expect(k==0x7ff, 0)) return x+x; /* NaN or Inf */ |
42 | if (__builtin_expect(n< -50000, 0)) |
43 | return tiny*copysign(tiny,x); /*underflow*/ |
44 | if (__builtin_expect(n> 50000 || k+n > 0x7fe, 0)) |
45 | return huge*copysign(huge,x); /* overflow */ |
46 | /* Now k and n are bounded we know that k = k+n does not |
47 | overflow. */ |
48 | k = k+n; |
49 | if (__builtin_expect(k > 0, 1)) /* normal result */ |
50 | {INSERT_WORDS64(x,(ix&UINT64_C(0x800fffffffffffff))|(k<<52)); |
51 | return x;} |
52 | if (k <= -54) |
53 | return tiny*copysign(tiny,x); /*underflow*/ |
54 | k += 54; /* subnormal result */ |
55 | INSERT_WORDS64(x,(ix&INT64_C(0x800fffffffffffff))|(k<<52)); |
56 | return x*twom54; |
57 | } |
58 | #ifdef NO_LONG_DOUBLE |
59 | strong_alias (__scalbn, __scalbnl) |
60 | #endif |
61 | |