1 | /* |
2 | * Copyright (c) 1998, 2001, Oracle and/or its affiliates. All rights reserved. |
3 | * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER. |
4 | * |
5 | * This code is free software; you can redistribute it and/or modify it |
6 | * under the terms of the GNU General Public License version 2 only, as |
7 | * published by the Free Software Foundation. Oracle designates this |
8 | * particular file as subject to the "Classpath" exception as provided |
9 | * by Oracle in the LICENSE file that accompanied this code. |
10 | * |
11 | * This code is distributed in the hope that it will be useful, but WITHOUT |
12 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
13 | * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
14 | * version 2 for more details (a copy is included in the LICENSE file that |
15 | * accompanied this code). |
16 | * |
17 | * You should have received a copy of the GNU General Public License version |
18 | * 2 along with this work; if not, write to the Free Software Foundation, |
19 | * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA. |
20 | * |
21 | * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA |
22 | * or visit www.oracle.com if you need additional information or have any |
23 | * questions. |
24 | */ |
25 | |
26 | /* |
27 | * scalbn (double x, int n) |
28 | * scalbn(x,n) returns x* 2**n computed by exponent |
29 | * manipulation rather than by actually performing an |
30 | * exponentiation or a multiplication. |
31 | */ |
32 | |
33 | #include "fdlibm.h" |
34 | |
35 | #ifdef __STDC__ |
36 | static const double |
37 | #else |
38 | static double |
39 | #endif |
40 | two54 = 1.80143985094819840000e+16, /* 0x43500000, 0x00000000 */ |
41 | twom54 = 5.55111512312578270212e-17, /* 0x3C900000, 0x00000000 */ |
42 | huge = 1.0e+300, |
43 | tiny = 1.0e-300; |
44 | |
45 | #ifdef __STDC__ |
46 | double scalbn (double x, int n) |
47 | #else |
48 | double scalbn (x,n) |
49 | double x; int n; |
50 | #endif |
51 | { |
52 | int k,hx,lx; |
53 | hx = __HI(x); |
54 | lx = __LO(x); |
55 | k = (hx&0x7ff00000)>>20; /* extract exponent */ |
56 | if (k==0) { /* 0 or subnormal x */ |
57 | if ((lx|(hx&0x7fffffff))==0) return x; /* +-0 */ |
58 | x *= two54; |
59 | hx = __HI(x); |
60 | k = ((hx&0x7ff00000)>>20) - 54; |
61 | if (n< -50000) return tiny*x; /*underflow*/ |
62 | } |
63 | if (k==0x7ff) return x+x; /* NaN or Inf */ |
64 | k = k+n; |
65 | if (k > 0x7fe) return huge*copysign(huge,x); /* overflow */ |
66 | if (k > 0) /* normal result */ |
67 | {__HI(x) = (hx&0x800fffff)|(k<<20); return x;} |
68 | if (k <= -54) { |
69 | if (n > 50000) /* in case integer overflow in n+k */ |
70 | return huge*copysign(huge,x); /*overflow*/ |
71 | else return tiny*copysign(tiny,x); /*underflow*/ |
72 | } |
73 | k += 54; /* subnormal result */ |
74 | __HI(x) = (hx&0x800fffff)|(k<<20); |
75 | return x*twom54; |
76 | } |
77 | |