s_ceill.c revision 176243
18478Swollman/* 28478Swollman * ==================================================== 38478Swollman * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 48478Swollman * 58478Swollman * Developed at SunPro, a Sun Microsystems, Inc. business. 68478Swollman * Permission to use, copy, modify, and distribute this 78478Swollman * software is freely granted, provided that this notice 88478Swollman * is preserved. 98478Swollman * ==================================================== 108478Swollman * 118478Swollman * From: @(#)s_ceil.c 5.1 93/09/24 128478Swollman */ 138478Swollman 148478Swollman#ifndef lint 158478Swollmanstatic char rcsid[] = "$FreeBSD: head/lib/msun/src/s_ceill.c 176243 2008-02-13 16:56:52Z bde $"; 168478Swollman#endif 178478Swollman 188478Swollman/* 198478Swollman * ceill(x) 208478Swollman * Return x rounded toward -inf to integral value 218478Swollman * Method: 228478Swollman * Bit twiddling. 238478Swollman * Exception: 248478Swollman * Inexact flag raised if x not equal to ceill(x). 258478Swollman */ 268478Swollman 278478Swollman#include <float.h> 288478Swollman#include <math.h> 298478Swollman#include <stdint.h> 30114589Sobrien 318478Swollman#include "fpmath.h" 3236999Scharnier 338478Swollman#ifdef LDBL_IMPLICIT_NBIT 348478Swollman#define MANH_SIZE (LDBL_MANH_SIZE + 1) 358478Swollman#define INC_MANH(u, c) do { \ 368478Swollman uint64_t o = u.bits.manh; \ 378478Swollman u.bits.manh += (c); \ 3836999Scharnier if (u.bits.manh < o) \ 39114589Sobrien u.bits.exp++; \ 4036999Scharnier} while (0) 41114589Sobrien#else 42114589Sobrien#define MANH_SIZE LDBL_MANH_SIZE 438478Swollman#define INC_MANH(u, c) do { \ 44136104Sdes uint64_t o = u.bits.manh; \ 45136104Sdes u.bits.manh += (c); \ 46136104Sdes if (u.bits.manh < o) { \ 47136104Sdes u.bits.exp++; \ 4836999Scharnier u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ 49136104Sdes } \ 50136104Sdes} while (0) 51136104Sdes#endif 528478Swollman 5378732Sddstatic const long double huge = 1.0e300; 5496381Salfred 55136104Sdeslong double 568478Swollmanceill(long double x) 578478Swollman{ 58136104Sdes union IEEEl2bits u = { .e = x }; 598478Swollman int e = u.bits.exp - LDBL_MAX_EXP + 1; 60136104Sdes 61136104Sdes if (e < MANH_SIZE - 1) { 62136104Sdes if (e < 0) { /* raise inexact if x != 0 */ 63136104Sdes if (huge + x > 0.0) 64136104Sdes if (u.bits.exp > 0 || 65136104Sdes (u.bits.manh | u.bits.manl) != 0) 66136104Sdes u.e = u.bits.sign ? -0.0 : 1.0; 67136104Sdes } else { 68136104Sdes uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); 69136104Sdes if (((u.bits.manh & m) | u.bits.manl) == 0) 70136104Sdes return (x); /* x is integral */ 71136104Sdes if (!u.bits.sign) { 72136104Sdes#ifdef LDBL_IMPLICIT_NBIT 73163852Sjhb if (e == 0) 74136104Sdes u.bits.exp++; 75158083Sps else 76136104Sdes#endif 77158083Sps INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); 78136104Sdes } 79163852Sjhb if (huge + x > 0.0) { /* raise inexact flag */ 80158083Sps u.bits.manh &= ~m; 81158083Sps u.bits.manl = 0; 82158083Sps } 83163852Sjhb } 84136104Sdes } else if (e < LDBL_MANT_DIG - 1) { 85136104Sdes uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); 86136104Sdes if ((u.bits.manl & m) == 0) 87136104Sdes return (x); /* x is integral */ 88136110Sdes if (!u.bits.sign) { 89136104Sdes if (e == MANH_SIZE - 1) 90136104Sdes INC_MANH(u, 1); 91136104Sdes else { 92136104Sdes uint64_t o = u.bits.manl; 93136104Sdes u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); 94136104Sdes if (u.bits.manl < o) /* got a carry */ 958478Swollman INC_MANH(u, 1); 9692542Simp } 978478Swollman } 98136104Sdes if (huge + x > 0.0) /* raise inexact flag */ 9993491Sphk u.bits.manl &= ~m; 10093491Sphk } 1018478Swollman return (u.e); 10224359Simp} 1038478Swollman 1048478Swollman#if LDBL_MANT_DIG == 53 1058478Swollman__weak_reference(ceil, ceill); 1068478Swollman#endif 1078478Swollman