s_ceill.c revision 140172
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 140172 2005-01-13 09:11:41Z stefanf $"; 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 53140172Sstefanflong double 54140172Sstefanfceill(long double x) 55140172Sstefanf{ 56140172Sstefanf union IEEEl2bits u = { .e = x }; 57140172Sstefanf int e = u.bits.exp - LDBL_MAX_EXP + 1; 58140172Sstefanf 59140172Sstefanf if (e < MANH_SIZE - 1) { 60140172Sstefanf if (e < 0) { /* raise inexact if x != 0 */ 61140172Sstefanf if (u.bits.exp > 0 || (u.bits.manh | u.bits.manl) != 0) 62140172Sstefanf u.e = u.bits.sign ? 0.0 : 1.0; 63140172Sstefanf } else { 64140172Sstefanf uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); 65140172Sstefanf if (((u.bits.manh & m) | u.bits.manl) == 0) 66140172Sstefanf return (x); /* x is integral */ 67140172Sstefanf /* raise inexact flag */ 68140172Sstefanf if (!u.bits.sign) { 69140172Sstefanf#ifdef LDBL_IMPLICIT_NBIT 70140172Sstefanf if (e == 0) 71140172Sstefanf u.bits.exp++; 72140172Sstefanf else 73140172Sstefanf#endif 74140172Sstefanf INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); 75140172Sstefanf } 76140172Sstefanf u.bits.manh &= ~m; 77140172Sstefanf u.bits.manl = 0; 78140172Sstefanf } 79140172Sstefanf } else if (e < LDBL_MANT_DIG - 1) { 80140172Sstefanf uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); 81140172Sstefanf if ((u.bits.manl & m) == 0) 82140172Sstefanf return (x); /* x is integral */ 83140172Sstefanf /* raise inexact flag */ 84140172Sstefanf if (!u.bits.sign) { 85140172Sstefanf if (e == MANH_SIZE - 1) 86140172Sstefanf INC_MANH(u, 1); 87140172Sstefanf else { 88140172Sstefanf uint64_t o = u.bits.manl; 89140172Sstefanf u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); 90140172Sstefanf if (u.bits.manl < o) /* got a carry */ 91140172Sstefanf INC_MANH(u, 1); 92140172Sstefanf } 93140172Sstefanf } 94140172Sstefanf u.bits.manl &= ~m; 95140172Sstefanf } 96140172Sstefanf return (u.e); 97140172Sstefanf} 98