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
12/*
13 * floorl(x)
14 * Return x rounded toward -inf to integral value
15 * Method:
16 *	Bit twiddling.
17 * Exception:
18 *	Inexact flag raised if x not equal to floorl(x).
19 */
20
21#include <float.h>
22#include <math.h>
23#include <stdint.h>
24
25#include "fpmath.h"
26
27#ifdef LDBL_IMPLICIT_NBIT
28#define	MANH_SIZE	(LDBL_MANH_SIZE + 1)
29#define	INC_MANH(u, c)	do {					\
30	uint64_t o = u.bits.manh;				\
31	u.bits.manh += (c);					\
32	if (u.bits.manh < o)					\
33		u.bits.exp++;					\
34} while (0)
35#else
36#define	MANH_SIZE	LDBL_MANH_SIZE
37#define	INC_MANH(u, c)	do {					\
38	uint64_t o = u.bits.manh;				\
39	u.bits.manh += (c);					\
40	if (u.bits.manh < o) {					\
41		u.bits.exp++;					\
42		u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1);	\
43	}							\
44} while (0)
45#endif
46
47static const long double huge = 1.0e300;
48
49long double
50floorl(long double x)
51{
52	union IEEEl2bits u = { .e = x };
53	int e = u.bits.exp - LDBL_MAX_EXP + 1;
54
55	if (e < MANH_SIZE - 1) {
56		if (e < 0) {			/* raise inexact if x != 0 */
57			if (huge + x > 0.0)
58				if (u.bits.exp > 0 ||
59				    (u.bits.manh | u.bits.manl) != 0)
60					u.e = u.bits.sign ? -1.0 : 0.0;
61		} else {
62			uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
63			if (((u.bits.manh & m) | u.bits.manl) == 0)
64				return (x);	/* x is integral */
65			if (u.bits.sign) {
66#ifdef LDBL_IMPLICIT_NBIT
67				if (e == 0)
68					u.bits.exp++;
69				else
70#endif
71				INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
72			}
73			if (huge + x > 0.0) {	/* raise inexact flag */
74				u.bits.manh &= ~m;
75				u.bits.manl = 0;
76			}
77		}
78	} else if (e < LDBL_MANT_DIG - 1) {
79		uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
80		if ((u.bits.manl & m) == 0)
81			return (x);	/* x is integral */
82		if (u.bits.sign) {
83			if (e == MANH_SIZE - 1)
84				INC_MANH(u, 1);
85			else {
86				uint64_t o = u.bits.manl;
87				u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
88				if (u.bits.manl < o)	/* got a carry */
89					INC_MANH(u, 1);
90			}
91		}
92		if (huge + x > 0.0)		/* raise inexact flag */
93			u.bits.manl &= ~m;
94	}
95	return (u.e);
96}
97