1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 *
11 * From: @(#)s_floor.c 5.1 93/09/24
12 */
13
14/*
15 * truncl(x)
16 * Return x rounded toward 0 to integral value
17 * Method:
18 *	Bit twiddling.
19 * Exception:
20 *	Inexact flag raised if x not equal to truncl(x).
21 */
22
23#include <sys/types.h>
24#include <machine/ieee.h>
25
26#include <float.h>
27#include <math.h>
28#include <stdint.h>
29
30#include "math_private.h"
31
32#ifdef LDBL_IMPLICIT_NBIT
33#define	MANH_SIZE	(EXT_FRACHBITS + EXT_FRACHMBITS + 1)
34#else
35#define	MANH_SIZE	(EXT_FRACHBITS + EXT_FRACHMBITS)
36#endif
37
38static const long double huge = 1.0e300;
39static const float zero[] = { 0.0, -0.0 };
40
41long double
42truncl(long double x)
43{
44	int e;
45	int64_t ix0, ix1;
46
47	GET_LDOUBLE_WORDS64(ix0,ix1,x);
48	e = ((ix0>>48)&0x7fff) - LDBL_MAX_EXP + 1;
49
50	if (e < MANH_SIZE - 1) {
51		if (e < 0) {			/* raise inexact if x != 0 */
52			if (huge + x > 0.0)
53				return (zero[((ix0>>48)&0x8000)!=0]);
54		} else {
55			uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
56			if (((ix0 & m) | ix1) == 0)
57				return (x);	/* x is integral */
58			if (huge + x > 0.0) {	/* raise inexact flag */
59				ix0 &= ~m;
60				ix1 = 0;
61			}
62		}
63	} else if (e < LDBL_MANT_DIG - 1) {
64		uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
65		if ((ix1 & m) == 0)
66			return (x);	/* x is integral */
67		if (huge + x > 0.0)		/* raise inexact flag */
68			ix1 &= ~m;
69	}
70	SET_LDOUBLE_WORDS64(x,ix0,ix1);
71	return (x);
72}
73