1140142Sstefanf/*
2140142Sstefanf * ====================================================
3140142Sstefanf * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4140142Sstefanf *
5140142Sstefanf * Developed at SunPro, a Sun Microsystems, Inc. business.
6140142Sstefanf * Permission to use, copy, modify, and distribute this
7140142Sstefanf * software is freely granted, provided that this notice
8140142Sstefanf * is preserved.
9140142Sstefanf * ====================================================
10140142Sstefanf *
11140142Sstefanf * From: @(#)s_floor.c 5.1 93/09/24
12140142Sstefanf */
13140142Sstefanf
14176245Sbde#include <sys/cdefs.h>
15176245Sbde__FBSDID("$FreeBSD$");
16140142Sstefanf
17140142Sstefanf/*
18140142Sstefanf * floorl(x)
19140142Sstefanf * Return x rounded toward -inf to integral value
20140142Sstefanf * Method:
21140142Sstefanf *	Bit twiddling.
22140142Sstefanf * Exception:
23140142Sstefanf *	Inexact flag raised if x not equal to floorl(x).
24140142Sstefanf */
25140142Sstefanf
26140142Sstefanf#include <float.h>
27140142Sstefanf#include <math.h>
28140142Sstefanf#include <stdint.h>
29140142Sstefanf
30140142Sstefanf#include "fpmath.h"
31140142Sstefanf
32140142Sstefanf#ifdef LDBL_IMPLICIT_NBIT
33140142Sstefanf#define	MANH_SIZE	(LDBL_MANH_SIZE + 1)
34140142Sstefanf#define	INC_MANH(u, c)	do {					\
35140142Sstefanf	uint64_t o = u.bits.manh;				\
36140142Sstefanf	u.bits.manh += (c);					\
37140142Sstefanf	if (u.bits.manh < o)					\
38140142Sstefanf		u.bits.exp++;					\
39140142Sstefanf} while (0)
40140142Sstefanf#else
41140142Sstefanf#define	MANH_SIZE	LDBL_MANH_SIZE
42140142Sstefanf#define	INC_MANH(u, c)	do {					\
43140142Sstefanf	uint64_t o = u.bits.manh;				\
44140142Sstefanf	u.bits.manh += (c);					\
45140142Sstefanf	if (u.bits.manh < o) {					\
46140142Sstefanf		u.bits.exp++;					\
47140142Sstefanf		u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1);	\
48140142Sstefanf	}							\
49140142Sstefanf} while (0)
50140142Sstefanf#endif
51140142Sstefanf
52145637Sstefanfstatic const long double huge = 1.0e300;
53145394Sstefanf
54140142Sstefanflong double
55140142Sstefanffloorl(long double x)
56140142Sstefanf{
57140142Sstefanf	union IEEEl2bits u = { .e = x };
58140142Sstefanf	int e = u.bits.exp - LDBL_MAX_EXP + 1;
59140142Sstefanf
60140142Sstefanf	if (e < MANH_SIZE - 1) {
61140142Sstefanf		if (e < 0) {			/* raise inexact if x != 0 */
62145637Sstefanf			if (huge + x > 0.0)
63145394Sstefanf				if (u.bits.exp > 0 ||
64145394Sstefanf				    (u.bits.manh | u.bits.manl) != 0)
65145394Sstefanf					u.e = u.bits.sign ? -1.0 : 0.0;
66140142Sstefanf		} else {
67140142Sstefanf			uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
68140142Sstefanf			if (((u.bits.manh & m) | u.bits.manl) == 0)
69140142Sstefanf				return (x);	/* x is integral */
70140142Sstefanf			if (u.bits.sign) {
71140142Sstefanf#ifdef LDBL_IMPLICIT_NBIT
72140142Sstefanf				if (e == 0)
73140142Sstefanf					u.bits.exp++;
74140142Sstefanf				else
75140142Sstefanf#endif
76140142Sstefanf				INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
77140142Sstefanf			}
78145637Sstefanf			if (huge + x > 0.0) {	/* raise inexact flag */
79145394Sstefanf				u.bits.manh &= ~m;
80145394Sstefanf				u.bits.manl = 0;
81145394Sstefanf			}
82140142Sstefanf		}
83140142Sstefanf	} else if (e < LDBL_MANT_DIG - 1) {
84140142Sstefanf		uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
85140142Sstefanf		if ((u.bits.manl & m) == 0)
86140142Sstefanf			return (x);	/* x is integral */
87140142Sstefanf		if (u.bits.sign) {
88140142Sstefanf			if (e == MANH_SIZE - 1)
89140142Sstefanf				INC_MANH(u, 1);
90140142Sstefanf			else {
91140142Sstefanf				uint64_t o = u.bits.manl;
92140142Sstefanf				u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
93140142Sstefanf				if (u.bits.manl < o)	/* got a carry */
94140142Sstefanf					INC_MANH(u, 1);
95140142Sstefanf			}
96140142Sstefanf		}
97145637Sstefanf		if (huge + x > 0.0)		/* raise inexact flag */
98145394Sstefanf			u.bits.manl &= ~m;
99140142Sstefanf	}
100140142Sstefanf	return (u.e);
101140142Sstefanf}
102