s_ceill.c revision 176236
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 14140172Sstefanf#ifndef lint 15140172Sstefanfstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/s_ceill.c 176236 2008-02-13 15:22:53Z bde $"; 16140172Sstefanf#endif 17140172Sstefanf 18140172Sstefanf/* 19140172Sstefanf * ceill(x) 20140172Sstefanf * Return x rounded toward -inf to integral value 21140172Sstefanf * Method: 22140172Sstefanf * Bit twiddling. 23140172Sstefanf * Exception: 24140172Sstefanf * Inexact flag raised if x not equal to ceill(x). 25140172Sstefanf */ 26140172Sstefanf 27140172Sstefanf#include <float.h> 28140172Sstefanf#include <math.h> 29140172Sstefanf#include <stdint.h> 30140172Sstefanf 31140172Sstefanf#include "fpmath.h" 32140172Sstefanf 33140172Sstefanf#ifdef LDBL_IMPLICIT_NBIT 34140172Sstefanf#define MANH_SIZE (LDBL_MANH_SIZE + 1) 35140172Sstefanf#define INC_MANH(u, c) do { \ 36140172Sstefanf uint64_t o = u.bits.manh; \ 37140172Sstefanf u.bits.manh += (c); \ 38140172Sstefanf if (u.bits.manh < o) \ 39140172Sstefanf u.bits.exp++; \ 40140172Sstefanf} while (0) 41140172Sstefanf#else 42140172Sstefanf#define MANH_SIZE LDBL_MANH_SIZE 43140172Sstefanf#define INC_MANH(u, c) do { \ 44140172Sstefanf uint64_t o = u.bits.manh; \ 45140172Sstefanf u.bits.manh += (c); \ 46140172Sstefanf if (u.bits.manh < o) { \ 47140172Sstefanf u.bits.exp++; \ 48140172Sstefanf u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ 49140172Sstefanf } \ 50140172Sstefanf} while (0) 51140172Sstefanf#endif 52140172Sstefanf 53145637Sstefanfstatic const long double huge = 1.0e300; 54145394Sstefanf 55140172Sstefanflong double 56140172Sstefanfceill(long double x) 57140172Sstefanf{ 58140172Sstefanf union IEEEl2bits u = { .e = x }; 59140172Sstefanf int e = u.bits.exp - LDBL_MAX_EXP + 1; 60140172Sstefanf 61140172Sstefanf if (e < MANH_SIZE - 1) { 62140172Sstefanf if (e < 0) { /* raise inexact if x != 0 */ 63145637Sstefanf if (huge + x > 0.0) 64145394Sstefanf if (u.bits.exp > 0 || 65145394Sstefanf (u.bits.manh | u.bits.manl) != 0) 66176236Sbde u.e = u.bits.sign ? -0.0 : 1.0; 67140172Sstefanf } else { 68140172Sstefanf uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); 69140172Sstefanf if (((u.bits.manh & m) | u.bits.manl) == 0) 70140172Sstefanf return (x); /* x is integral */ 71140172Sstefanf if (!u.bits.sign) { 72140172Sstefanf#ifdef LDBL_IMPLICIT_NBIT 73140172Sstefanf if (e == 0) 74140172Sstefanf u.bits.exp++; 75140172Sstefanf else 76140172Sstefanf#endif 77140172Sstefanf INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); 78140172Sstefanf } 79145637Sstefanf if (huge + x > 0.0) { /* raise inexact flag */ 80145394Sstefanf u.bits.manh &= ~m; 81145394Sstefanf u.bits.manl = 0; 82145394Sstefanf } 83140172Sstefanf } 84140172Sstefanf } else if (e < LDBL_MANT_DIG - 1) { 85140172Sstefanf uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); 86140172Sstefanf if ((u.bits.manl & m) == 0) 87140172Sstefanf return (x); /* x is integral */ 88140172Sstefanf if (!u.bits.sign) { 89140172Sstefanf if (e == MANH_SIZE - 1) 90140172Sstefanf INC_MANH(u, 1); 91140172Sstefanf else { 92140172Sstefanf uint64_t o = u.bits.manl; 93140172Sstefanf u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); 94140172Sstefanf if (u.bits.manl < o) /* got a carry */ 95140172Sstefanf INC_MANH(u, 1); 96140172Sstefanf } 97140172Sstefanf } 98145637Sstefanf if (huge + x > 0.0) /* raise inexact flag */ 99145394Sstefanf u.bits.manl &= ~m; 100140172Sstefanf } 101140172Sstefanf return (u.e); 102140172Sstefanf} 103