1 /* $NetBSD: s_scalbnl.c,v 1.9 2013/05/20 19:40:09 joerg Exp $ */
2
3 /*-
4 * Copyright (c) 2011 The NetBSD Foundation, Inc.
5 * All rights reserved.
6 *
7 * This code is derived from software contributed to The NetBSD Foundation
8 * by Joerg Sonnenberger.
9 *
10 * Redistribution and use in source and binary forms, with or without
11 * modification, are permitted provided that the following conditions
12 * are met:
13 * 1. Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in the
17 * documentation and/or other materials provided with the distribution.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
20 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
21 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
22 * PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
23 * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 */
31
32 #include <sys/cdefs.h>
33 __RCSID("$NetBSD: s_scalbnl.c,v 1.9 2013/05/20 19:40:09 joerg Exp $");
34
35 #include "namespace.h"
36
37 #include <float.h>
38 #include <math.h>
39 #include <machine/ieee.h>
40
41 #ifdef __HAVE_LONG_DOUBLE
42
43 #ifdef _LP64
44 long double
scalbnl(long double x,int n)45 scalbnl(long double x, int n)
46 {
47 return scalblnl(x, n);
48 }
49 #else
__strong_alias(_scalbnl,_scalblnl)50 __strong_alias(_scalbnl, _scalblnl)
51 #endif
52
53 __weak_alias(scalbnl, _scalbnl)
54 __weak_alias(scalblnl, _scalblnl)
55 __weak_alias(ldexpl, _scalbnl)
56
57 #if LDBL_MANT_DIG == 64
58 #define FROM_UNDERFLOW 0x1p65L
59 #define TO_UNDERFLOW 0x1p-65L
60 #elif LDBL_MANT_DIG == 113
61 #define FROM_UNDERFLOW 0x1p114L
62 #define TO_UNDERFLOW 0x1p-114L
63 #else
64 #error Unsupported long double format
65 #endif
66
67 long double
68 scalblnl(long double x, long n)
69 {
70 union ieee_ext_u u;
71
72 /* Trivial cases first */
73 if (n == 0 || x == 0.0L)
74 return x;
75
76 u.extu_ld = x;
77
78 /* NaN and infinite don't change either, but trigger exception */
79 if (u.extu_ext.ext_exp == EXT_EXP_INFNAN)
80 return x + x;
81
82 /* Protect against integer overflow in calculation of new exponent */
83 if (n > LDBL_MAX_EXP - LDBL_MIN_EXP + LDBL_MANT_DIG)
84 goto overflow;
85 if (n < LDBL_MIN_EXP - LDBL_MAX_EXP - LDBL_MANT_DIG)
86 goto underflow;
87
88 /* Scale denormalized numbers slightly, so that they are normal */
89 if (u.extu_ext.ext_exp == 0) {
90 u.extu_ld *= FROM_UNDERFLOW;
91 n -= LDBL_MANT_DIG + 1;
92 }
93
94 n += u.extu_ext.ext_exp;
95 if (n >= LDBL_MAX_EXP + EXT_EXP_BIAS)
96 goto overflow;
97 /* Positive exponent (incl. bias) means normal result */
98 if (n > 0) {
99 u.extu_ext.ext_exp = n;
100 return u.extu_ld;
101 }
102 /* Shift the exponent and let the multiply below handle subnormal */
103 n += LDBL_MANT_DIG + 1;
104 if (n <= 0)
105 goto underflow;
106 u.extu_ext.ext_exp = n;
107 return u.extu_ld * TO_UNDERFLOW;
108
109 underflow:
110 return LDBL_MIN * copysignl(LDBL_MIN, x);
111
112 overflow:
113 return LDBL_MAX * copysignl(LDBL_MAX, x);
114 }
115
116 #endif
117