1 | /* s_scalblnl.c -- long double version of s_scalbln.c. |
2 | */ |
3 | |
4 | /* @(#)s_scalbln.c 5.1 93/09/24 */ |
5 | /* |
6 | * ==================================================== |
7 | * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. |
8 | * |
9 | * Developed at SunPro, a Sun Microsystems, Inc. business. |
10 | * Permission to use, copy, modify, and distribute this |
11 | * software is freely granted, provided that this notice |
12 | * is preserved. |
13 | * ==================================================== |
14 | */ |
15 | |
16 | #if defined(LIBM_SCCS) && !defined(lint) |
17 | static char rcsid[] = "$NetBSD: $" ; |
18 | #endif |
19 | |
20 | /* |
21 | * scalblnl (long double x, long int n) |
22 | * scalblnl(x,n) returns x* 2**n computed by exponent |
23 | * manipulation rather than by actually performing an |
24 | * exponentiation or a multiplication. |
25 | */ |
26 | |
27 | #include <math.h> |
28 | #include <math_private.h> |
29 | #include <math_ldbl_opt.h> |
30 | |
31 | static const long double |
32 | twolm54 = 5.55111512312578270212e-17, /* 0x3C90000000000000, 0 */ |
33 | huge = 1.0E+300L, |
34 | tiny = 1.0E-300L; |
35 | static const double |
36 | two54 = 1.80143985094819840000e+16, /* 0x4350000000000000 */ |
37 | twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */ |
38 | |
39 | long double __scalblnl (long double x, long int n) |
40 | { |
41 | int64_t k,l,hx,lx; |
42 | union { int64_t i; double d; } u; |
43 | double xhi, xlo; |
44 | |
45 | ldbl_unpack (x, &xhi, &xlo); |
46 | EXTRACT_WORDS64 (hx, xhi); |
47 | EXTRACT_WORDS64 (lx, xlo); |
48 | k = (hx>>52)&0x7ff; /* extract exponent */ |
49 | l = (lx>>52)&0x7ff; |
50 | if (k==0) { /* 0 or subnormal x */ |
51 | if ((hx&0x7fffffffffffffffULL)==0) return x; /* +-0 */ |
52 | u.i = hx; |
53 | u.d *= two54; |
54 | hx = u.i; |
55 | k = ((hx>>52)&0x7ff) - 54; |
56 | } |
57 | else if (k==0x7ff) return x+x; /* NaN or Inf */ |
58 | if (n< -50000) return tiny*copysignl(tiny,x); /*underflow */ |
59 | if (n> 50000 || k+n > 0x7fe) |
60 | return huge*copysignl(huge,x); /* overflow */ |
61 | /* Now k and n are bounded we know that k = k+n does not |
62 | overflow. */ |
63 | k = k+n; |
64 | if (k > 0) { /* normal result */ |
65 | hx = (hx&0x800fffffffffffffULL)|(k<<52); |
66 | if ((lx & 0x7fffffffffffffffULL) == 0) { /* low part +-0 */ |
67 | INSERT_WORDS64 (xhi, hx); |
68 | INSERT_WORDS64 (xlo, lx); |
69 | x = ldbl_pack (xhi, xlo); |
70 | return x; |
71 | } |
72 | if (l == 0) { /* low part subnormal */ |
73 | u.i = lx; |
74 | u.d *= two54; |
75 | lx = u.i; |
76 | l = ((lx>>52)&0x7ff) - 54; |
77 | } |
78 | l = l + n; |
79 | if (l > 0) |
80 | lx = (lx&0x800fffffffffffffULL)|(l<<52); |
81 | else if (l <= -54) |
82 | lx = (lx&0x8000000000000000ULL); |
83 | else { |
84 | l += 54; |
85 | u.i = (lx&0x800fffffffffffffULL)|(l<<52); |
86 | u.d *= twom54; |
87 | lx = u.i; |
88 | } |
89 | INSERT_WORDS64 (xhi, hx); |
90 | INSERT_WORDS64 (xlo, lx); |
91 | x = ldbl_pack (xhi, xlo); |
92 | return x; |
93 | } |
94 | if (k <= -54) |
95 | return tiny*copysignl(tiny,x); /*underflow*/ |
96 | k += 54; /* subnormal result */ |
97 | lx &= 0x8000000000000000ULL; |
98 | hx &= 0x800fffffffffffffULL; |
99 | INSERT_WORDS64 (xhi, hx|(k<<52)); |
100 | INSERT_WORDS64 (xlo, lx); |
101 | x = ldbl_pack (xhi, xlo); |
102 | return x*twolm54; |
103 | } |
104 | |