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