1/*	$NetBSD: s_scalbnl.c,v 1.8 2013/05/19 20:50:02 martin 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.8 2013/05/19 20:50:02 martin 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
44long double
45scalbnl(long double x, int n)
46{
47	return scalblnl(x, n);
48}
49#else
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
67long double
68scalblnl(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
109underflow:
110	return LDBL_MIN * copysignl(LDBL_MIN, x);
111
112overflow:
113	return LDBL_MAX * copysignl(LDBL_MAX, x);
114}
115
116#endif
117