]>
Commit | Line | Data |
---|---|---|
f964490f | 1 | /* s_scalblnl.c -- long double version of s_scalbln.c. |
f964490f RM |
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 | ||
1ed0291c RH |
27 | #include <math.h> |
28 | #include <math_private.h> | |
f964490f RM |
29 | #include <math_ldbl_opt.h> |
30 | ||
f964490f | 31 | static const long double |
f964490f RM |
32 | twolm54 = 5.55111512312578270212e-17, /* 0x3C90000000000000, 0 */ |
33 | huge = 1.0E+300L, | |
34 | tiny = 1.0E-300L; | |
f964490f | 35 | static const double |
f964490f RM |
36 | two54 = 1.80143985094819840000e+16, /* 0x4350000000000000 */ |
37 | twom54 = 5.55111512312578270212e-17; /* 0x3C90000000000000 */ | |
38 | ||
8db21882 | 39 | long double __scalblnl (long double x, long int n) |
f964490f RM |
40 | { |
41 | int64_t k,l,hx,lx; | |
42 | union { int64_t i; double d; } u; | |
765714ca AM |
43 | double xhi, xlo; |
44 | ||
45 | ldbl_unpack (x, &xhi, &xlo); | |
46 | EXTRACT_WORDS64 (hx, xhi); | |
47 | EXTRACT_WORDS64 (lx, xlo); | |
f964490f RM |
48 | k = (hx>>52)&0x7ff; /* extract exponent */ |
49 | l = (lx>>52)&0x7ff; | |
50 | if (k==0) { /* 0 or subnormal x */ | |
765714ca | 51 | if ((hx&0x7fffffffffffffffULL)==0) return x; /* +-0 */ |
f964490f RM |
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 */ | |
81dca813 | 58 | if (n< -50000) return tiny*copysignl(tiny,x); /*underflow */ |
a6d06d7b | 59 | if (n> 50000 || k+n > 0x7fe) |
81dca813 | 60 | return huge*copysignl(huge,x); /* overflow */ |
a6d06d7b JM |
61 | /* Now k and n are bounded we know that k = k+n does not |
62 | overflow. */ | |
63 | k = k+n; | |
f964490f RM |
64 | if (k > 0) { /* normal result */ |
65 | hx = (hx&0x800fffffffffffffULL)|(k<<52); | |
66 | if ((lx & 0x7fffffffffffffffULL) == 0) { /* low part +-0 */ | |
765714ca AM |
67 | INSERT_WORDS64 (xhi, hx); |
68 | INSERT_WORDS64 (xlo, lx); | |
69 | x = ldbl_pack (xhi, xlo); | |
350635a5 | 70 | return x; |
f964490f RM |
71 | } |
72 | if (l == 0) { /* low part subnormal */ | |
350635a5 OB |
73 | u.i = lx; |
74 | u.d *= two54; | |
75 | lx = u.i; | |
76 | l = ((lx>>52)&0x7ff) - 54; | |
f964490f RM |
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 { | |
350635a5 OB |
84 | l += 54; |
85 | u.i = (lx&0x800fffffffffffffULL)|(l<<52); | |
86 | u.d *= twom54; | |
87 | lx = u.i; | |
f964490f | 88 | } |
765714ca AM |
89 | INSERT_WORDS64 (xhi, hx); |
90 | INSERT_WORDS64 (xlo, lx); | |
91 | x = ldbl_pack (xhi, xlo); | |
f964490f RM |
92 | return x; |
93 | } | |
94 | if (k <= -54) | |
81dca813 | 95 | return tiny*copysignl(tiny,x); /*underflow*/ |
f964490f RM |
96 | k += 54; /* subnormal result */ |
97 | lx &= 0x8000000000000000ULL; | |
765714ca AM |
98 | hx &= 0x800fffffffffffffULL; |
99 | INSERT_WORDS64 (xhi, hx|(k<<52)); | |
100 | INSERT_WORDS64 (xlo, lx); | |
101 | x = ldbl_pack (xhi, xlo); | |
f964490f RM |
102 | return x*twolm54; |
103 | } |