1140172Sstefanf/* 2140172Sstefanf * ==================================================== 3140172Sstefanf * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 4140172Sstefanf * 5140172Sstefanf * Developed at SunPro, a Sun Microsystems, Inc. business. 6140172Sstefanf * Permission to use, copy, modify, and distribute this 7140172Sstefanf * software is freely granted, provided that this notice 8140172Sstefanf * is preserved. 9140172Sstefanf * ==================================================== 10140172Sstefanf * 11140172Sstefanf * From: @(#)s_ceil.c 5.1 93/09/24 12140172Sstefanf */ 13140172Sstefanf 14176245Sbde#include <sys/cdefs.h> 15176245Sbde__FBSDID("$FreeBSD: releng/10.2/lib/msun/src/s_ceill.c 176280 2008-02-14 15:10:34Z bde $"); 16140172Sstefanf 17140172Sstefanf/* 18140172Sstefanf * ceill(x) 19140172Sstefanf * Return x rounded toward -inf to integral value 20140172Sstefanf * Method: 21140172Sstefanf * Bit twiddling. 22140172Sstefanf * Exception: 23140172Sstefanf * Inexact flag raised if x not equal to ceill(x). 24140172Sstefanf */ 25140172Sstefanf 26140172Sstefanf#include <float.h> 27140172Sstefanf#include <math.h> 28140172Sstefanf#include <stdint.h> 29140172Sstefanf 30140172Sstefanf#include "fpmath.h" 31140172Sstefanf 32140172Sstefanf#ifdef LDBL_IMPLICIT_NBIT 33140172Sstefanf#define MANH_SIZE (LDBL_MANH_SIZE + 1) 34140172Sstefanf#define INC_MANH(u, c) do { \ 35140172Sstefanf uint64_t o = u.bits.manh; \ 36140172Sstefanf u.bits.manh += (c); \ 37140172Sstefanf if (u.bits.manh < o) \ 38140172Sstefanf u.bits.exp++; \ 39140172Sstefanf} while (0) 40140172Sstefanf#else 41140172Sstefanf#define MANH_SIZE LDBL_MANH_SIZE 42140172Sstefanf#define INC_MANH(u, c) do { \ 43140172Sstefanf uint64_t o = u.bits.manh; \ 44140172Sstefanf u.bits.manh += (c); \ 45140172Sstefanf if (u.bits.manh < o) { \ 46140172Sstefanf u.bits.exp++; \ 47140172Sstefanf u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ 48140172Sstefanf } \ 49140172Sstefanf} while (0) 50140172Sstefanf#endif 51140172Sstefanf 52145637Sstefanfstatic const long double huge = 1.0e300; 53145394Sstefanf 54140172Sstefanflong double 55140172Sstefanfceill(long double x) 56140172Sstefanf{ 57140172Sstefanf union IEEEl2bits u = { .e = x }; 58140172Sstefanf int e = u.bits.exp - LDBL_MAX_EXP + 1; 59140172Sstefanf 60140172Sstefanf if (e < MANH_SIZE - 1) { 61140172Sstefanf 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) 65176236Sbde u.e = u.bits.sign ? -0.0 : 1.0; 66140172Sstefanf } else { 67140172Sstefanf uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); 68140172Sstefanf if (((u.bits.manh & m) | u.bits.manl) == 0) 69140172Sstefanf return (x); /* x is integral */ 70140172Sstefanf if (!u.bits.sign) { 71140172Sstefanf#ifdef LDBL_IMPLICIT_NBIT 72140172Sstefanf if (e == 0) 73140172Sstefanf u.bits.exp++; 74140172Sstefanf else 75140172Sstefanf#endif 76140172Sstefanf INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); 77140172Sstefanf } 78145637Sstefanf if (huge + x > 0.0) { /* raise inexact flag */ 79145394Sstefanf u.bits.manh &= ~m; 80145394Sstefanf u.bits.manl = 0; 81145394Sstefanf } 82140172Sstefanf } 83140172Sstefanf } else if (e < LDBL_MANT_DIG - 1) { 84140172Sstefanf uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); 85140172Sstefanf if ((u.bits.manl & m) == 0) 86140172Sstefanf return (x); /* x is integral */ 87140172Sstefanf if (!u.bits.sign) { 88140172Sstefanf if (e == MANH_SIZE - 1) 89140172Sstefanf INC_MANH(u, 1); 90140172Sstefanf else { 91140172Sstefanf uint64_t o = u.bits.manl; 92140172Sstefanf u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); 93140172Sstefanf if (u.bits.manl < o) /* got a carry */ 94140172Sstefanf INC_MANH(u, 1); 95140172Sstefanf } 96140172Sstefanf } 97145637Sstefanf if (huge + x > 0.0) /* raise inexact flag */ 98145394Sstefanf u.bits.manl &= ~m; 99140172Sstefanf } 100140172Sstefanf return (u.e); 101140172Sstefanf} 102