s_floorl.c revision 140142
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 14140142Sstefanf#ifndef lint 15140142Sstefanfstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/s_floorl.c 140142 2005-01-12 22:10:46Z stefanf $"; 16140142Sstefanf#endif 17140142Sstefanf 18140142Sstefanf/* 19140142Sstefanf * floorl(x) 20140142Sstefanf * Return x rounded toward -inf to integral value 21140142Sstefanf * Method: 22140142Sstefanf * Bit twiddling. 23140142Sstefanf * Exception: 24140142Sstefanf * Inexact flag raised if x not equal to floorl(x). 25140142Sstefanf */ 26140142Sstefanf 27140142Sstefanf#include <float.h> 28140142Sstefanf#include <math.h> 29140142Sstefanf#include <stdint.h> 30140142Sstefanf 31140142Sstefanf#include "fpmath.h" 32140142Sstefanf 33140142Sstefanf#ifdef LDBL_IMPLICIT_NBIT 34140142Sstefanf#define MANH_SIZE (LDBL_MANH_SIZE + 1) 35140142Sstefanf#define INC_MANH(u, c) do { \ 36140142Sstefanf uint64_t o = u.bits.manh; \ 37140142Sstefanf u.bits.manh += (c); \ 38140142Sstefanf if (u.bits.manh < o) \ 39140142Sstefanf u.bits.exp++; \ 40140142Sstefanf} while (0) 41140142Sstefanf#else 42140142Sstefanf#define MANH_SIZE LDBL_MANH_SIZE 43140142Sstefanf#define INC_MANH(u, c) do { \ 44140142Sstefanf uint64_t o = u.bits.manh; \ 45140142Sstefanf u.bits.manh += (c); \ 46140142Sstefanf if (u.bits.manh < o) { \ 47140142Sstefanf u.bits.exp++; \ 48140142Sstefanf u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ 49140142Sstefanf } \ 50140142Sstefanf} while (0) 51140142Sstefanf#endif 52140142Sstefanf 53140142Sstefanflong double 54140142Sstefanffloorl(long double x) 55140142Sstefanf{ 56140142Sstefanf union IEEEl2bits u = { .e = x }; 57140142Sstefanf int e = u.bits.exp - LDBL_MAX_EXP + 1; 58140142Sstefanf 59140142Sstefanf if (e < MANH_SIZE - 1) { 60140142Sstefanf if (e < 0) { /* raise inexact if x != 0 */ 61140142Sstefanf if (u.bits.exp > 0 || (u.bits.manh | u.bits.manl) != 0) 62140142Sstefanf u.e = u.bits.sign ? -1.0 : 0.0; 63140142Sstefanf } else { 64140142Sstefanf uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); 65140142Sstefanf if (((u.bits.manh & m) | u.bits.manl) == 0) 66140142Sstefanf return (x); /* x is integral */ 67140142Sstefanf /* raise inexact flag */ 68140142Sstefanf if (u.bits.sign) { 69140142Sstefanf#ifdef LDBL_IMPLICIT_NBIT 70140142Sstefanf if (e == 0) 71140142Sstefanf u.bits.exp++; 72140142Sstefanf else 73140142Sstefanf#endif 74140142Sstefanf INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); 75140142Sstefanf } 76140142Sstefanf u.bits.manh &= ~m; 77140142Sstefanf u.bits.manl = 0; 78140142Sstefanf } 79140142Sstefanf } else if (e < LDBL_MANT_DIG - 1) { 80140142Sstefanf uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); 81140142Sstefanf if ((u.bits.manl & m) == 0) 82140142Sstefanf return (x); /* x is integral */ 83140142Sstefanf /* raise inexact flag */ 84140142Sstefanf if (u.bits.sign) { 85140142Sstefanf if (e == MANH_SIZE - 1) 86140142Sstefanf INC_MANH(u, 1); 87140142Sstefanf else { 88140142Sstefanf uint64_t o = u.bits.manl; 89140142Sstefanf u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); 90140142Sstefanf if (u.bits.manl < o) /* got a carry */ 91140142Sstefanf INC_MANH(u, 1); 92140142Sstefanf } 93140142Sstefanf } 94140142Sstefanf u.bits.manl &= ~m; 95140142Sstefanf } 96140142Sstefanf return (u.e); 97140142Sstefanf} 98