xref: /netbsd-src/lib/libm/src/s_truncl.c (revision a4ddc2c8fb9af816efe3b1c375a5530aef0e89e9)
1 /*-
2  * Copyright (c) 2012 The NetBSD Foundation, Inc.
3  * All rights reserved.
4  *
5  * This code is derived from software contributed to The NetBSD Foundation
6  * by Matt Thomas of 3am Software Foundry.
7  *
8  * Redistribution and use in source and binary forms, with or without
9  * modification, are permitted provided that the following conditions
10  * are met:
11  * 1. Redistributions of source code must retain the above copyright
12  *    notice, this list of conditions and the following disclaimer.
13  * 2. Redistributions in binary form must reproduce the above copyright
14  *    notice, this list of conditions and the following disclaimer in the
15  *    documentation and/or other materials provided with the distribution.
16  *
17  * THIS SOFTWARE IS PROVIDED BY THE NETBSD FOUNDATION, INC. AND CONTRIBUTORS
18  * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
19  * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
20  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE FOUNDATION OR CONTRIBUTORS
21  * BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
22  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
23  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
24  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
25  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
26  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
27  * POSSIBILITY OF SUCH DAMAGE.
28  */
29 
30 #include <sys/cdefs.h>
31 #if defined(LIBM_SCCS) && !defined(lint)
32 __RCSID("$NetBSD: s_truncl.c,v 1.2 2012/08/08 16:58:28 matt Exp $");
33 #endif
34 
35 /*
36  * trunc(x)
37  * Return x rounded toward 0 to integral value
38  * Method:
39  *	Bit twiddling.
40  * Exception:
41  *	Inexact flag raised if x was not equal to trunc(x).
42  */
43 
44 #include <machine/ieee.h>
45 #include <float.h>
46 #include "math.h"
47 #include "math_private.h"
48 
49 #ifdef __HAVE_LONG_DOUBLE
50 
51 static const long double huge = LDBL_MAX;
52 
53 long double
54 truncl(long double x)
55 {
56 	struct ieee_ext_u ux = { .extu_ld = x, };
57 	int32_t exponent = ux.extu_exp - EXT_EXP_BIAS;
58 #ifdef LDBL_MANT_DIG == EXT_FRACBITS
59 	/*
60 	 * If there is no hidden bit, don't count it
61 	 */
62 	const u_int frach_bits = EXT_FRACHBITS - 1;
63 	const u_int frac_bits = EXT_FRACBITS - 1;
64 #else
65 	const u_int frach_bits = EXT_FRACHBITS;
66 	const u_int frac_bits = EXT_FRACBITS;
67 #endif
68 
69 	/*
70 	 * If this number is big enough to have no fractional digits...
71 	 * (or is an inf or nan).
72 	 */
73 	if (exponent >= frac_bits) {
74 		if (exponent == EXT_EXP_NAN - EXT_EXP_BIAS)
75 			return x+x;	/* x is inf or nan */
76 		return x;		/* x is integral */
77 	}
78 
79 	/*
80 	 * If this number is too small enough to have any integral digits...
81 	 */
82 	if (exponent < 0 && (huge - x > 0.0 || true)) {
83 		/* set inexact if x != 0 */
84 		/* |x|<1, so return 0*sign(x) */
85 		return ux.extu_sign ? -0.0 : 0.0;
86 	}
87 
88 	uint32_t frach_mask = __BIT(frach_bits) - 1;
89 #ifdef EXT_FRACHMBITS
90 	uint32_t frachm_mask = __BIT(EXT_FRACHMBITS) - 1;
91 #endif
92 #ifdef EXT_FRACHMBITS
93 	uint32_t frachl_mask = __BIT(EXT_FRACLMBITS) - 1;
94 #endif
95 	uint32_t fracl_mask = __BIT(EXT_FRACLBITS) - 1;
96 
97 	if (exponent < frach_bits) {
98 		frach_mask >>= exponent;
99 #ifdef EXT_FRACHMBITS
100 	} else if (exponent < frach_bits + EXT_FRACHM_BITS) {
101 		frach_mask = 0;		exponent -= frach_bits;
102 		frachm_mask >>= exponent;
103 #endif
104 #ifdef EXT_FRACLMBITS
105 	} else if (exponent < frach_bits + EXT_FRACHM_BITS + EXT_FRACLMBITS) {
106 		frach_mask = 0;		exponent -= frach_bits;
107 		frachm_mask = 0;	exponent -= EXT_FRACHMBITS;
108 		fraclm_mask >>= exponent;
109 #endif
110 	} else {
111 		frach_mask = 0;		exponent -= frach_bits;
112 #ifdef EXT_FRACHMBITS
113 		frachm_mask = 0;	exponent -= EXT_FRACHMBITS;
114 #endif
115 #ifdef EXT_FRACLMBITS
116 		fraclm_mask = 0;	exponent -= EXT_FRACLMBITS;
117 #endif
118 		fraclm_mask >>= exponent;
119 	}
120 
121 	if ((ux.extu_frach & frach_mask) == 0
122 #ifdef EXT_FRACHMBITS
123 	    && (ux.extu_frachm & frachm_mask) == 0
124 #endif
125 #ifdef EXT_FRACLMBITS
126 	    && (ux.extu_fraclm & frachm_mask) == 0
127 #endif
128 	    && (ux.extu_fracl & fracl_mask) == 0)
129 		return x; /* x is integral */
130 
131 	if (huge - x > 0.0 || true) {		/* set inexact flag */
132 		/*
133 		 * Clear any fractional bits...
134 		 */
135 		ux.extu_frach &= ~frach_mask;
136 #ifdef EXT_FRACHMBITS
137 		ux.extu_frachm &= ~frachm_mask;
138 #endif
139 #ifdef EXT_FRACLMBITS
140 		ux.extu_fraclm &= ~fraclm_mask;
141 #endif
142 		ux.extu_fracl &= ~fracl_mask;
143 		return ux.extu_ld;
144 	}
145 }
146 
147 #endif	/* __HAVE_LONG_DOUBLE */
148