Home | History | Annotate | Line # | Download | only in src 
      1 /*
      2  * ====================================================
      3  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
      4  *
      5  * Developed at SunPro, a Sun Microsystems, Inc. business.
      6  * Permission to use, copy, modify, and distribute this
      7  * software is freely granted, provided that this notice
      8  * is preserved.
      9  * ====================================================
     10  *
     11  * From: @(#)s_ceil.c 5.1 93/09/24
     12  */
     13 
     14 #include <sys/cdefs.h>
     15 __RCSID("$NetBSD: s_ceill.c,v 1.1 2013/11/11 23:57:34 joerg Exp $");
     16 #if 0
     17 __FBSDID("$FreeBSD: head/lib/msun/src/s_ceill.c 176280 2008-02-14 15:10:34Z bde $");
     18 #endif
     19 
     20 /*
     21  * ceill(x)
     22  * Return x rounded toward -inf to integral value
     23  * Method:
     24  *	Bit twiddling.
     25  * Exception:
     26  *	Inexact flag raised if x not equal to ceill(x).
     27  */
     28 #include "namespace.h"
     29 
     30 #include <float.h>
     31 #include <math.h>
     32 #include <stdint.h>
     33 #include <machine/ieee.h>
     34 
     35 #ifdef __HAVE_LONG_DOUBLE
     36 
     37 #ifdef __weak_alias
     38 __weak_alias(ceill, _ceill)
     39 #endif
     40 
     41 #ifdef LDBL_IMPLICIT_NBIT
     42 #define	MANH_SIZE	(EXT_FRACHBITS + 1)
     43 #define	INC_MANH(ux, c)	do {					\
     44 	uint64_t oi = ux.extu_frach;				\
     45 	ux.extu_frach += (c);					\
     46 	if (ux.extu_frach < oi)					\
     47 		ux.extu_exp++;					\
     48 } while (0)
     49 #else
     50 #define	MANH_SIZE	EXT_FRACHBITS
     51 #define	INC_MANH(ux, c)	do {					\
     52 	uint64_t oi = ux.extu_frach;				\
     53 	ux.extu_frach += (c);					\
     54 	if (ux.extu_frach < oi) {				\
     55 		ux.extu_exp++;					\
     56 		ux.extu_frach |= 1llu << (EXT_FRACHBITS - 1);	\
     57 	}							\
     58 } while (0)
     59 #endif
     60 
     61 static const long double huge = 1.0e300;
     62 
     63 long double
     64 ceill(long double x)
     65 {
     66 	union ieee_ext_u ux = { .extu_ld = x, };
     67 	int e = ux.extu_exp - LDBL_MAX_EXP + 1;
     68 
     69 	if (e < MANH_SIZE - 1) {
     70 		if (e < 0) {			/* raise inexact if x != 0 */
     71 			if (huge + x > 0.0)
     72 				if (ux.extu_exp > 0 ||
     73 				    (ux.extu_frach | ux.extu_fracl) != 0)
     74 					ux.extu_ld = ux.extu_sign ? -0.0 : 1.0;
     75 		} else {
     76 			uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
     77 			if (((ux.extu_frach & m) | ux.extu_fracl) == 0)
     78 				return (x);	/* x is integral */
     79 			if (!ux.extu_sign) {
     80 #ifdef LDBL_IMPLICIT_NBIT
     81 				if (e == 0)
     82 					ux.extu_exp++;
     83 				else
     84 #endif
     85 				INC_MANH(ux, 1llu << (MANH_SIZE - e - 1));
     86 			}
     87 			if (huge + x > 0.0) {	/* raise inexact flag */
     88 				ux.extu_frach &= ~m;
     89 				ux.extu_fracl = 0;
     90 			}
     91 		}
     92 	} else if (e < LDBL_MANT_DIG - 1) {
     93 		uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
     94 		if ((ux.extu_fracl & m) == 0)
     95 			return (x);	/* x is integral */
     96 		if (!ux.extu_sign) {
     97 			if (e == MANH_SIZE - 1)
     98 				INC_MANH(ux, 1);
     99 			else {
    100 				uint64_t o = ux.extu_fracl;
    101 				ux.extu_fracl += 1llu << (LDBL_MANT_DIG - e - 1);
    102 				if (ux.extu_fracl < o)	/* got a carry */
    103 					INC_MANH(ux, 1);
    104 			}
    105 		}
    106 		if (huge + x > 0.0)		/* raise inexact flag */
    107 			ux.extu_fracl &= ~m;
    108 	}
    109 	return (ux.extu_ld);
    110 }
    111 
    112 #endif /* __HAVE_LONG_DOUBLE */
    113