143295facSStefan Farfeleder /*
243295facSStefan Farfeleder * ====================================================
343295facSStefan Farfeleder * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
443295facSStefan Farfeleder *
543295facSStefan Farfeleder * Developed at SunPro, a Sun Microsystems, Inc. business.
643295facSStefan Farfeleder * Permission to use, copy, modify, and distribute this
743295facSStefan Farfeleder * software is freely granted, provided that this notice
843295facSStefan Farfeleder * is preserved.
943295facSStefan Farfeleder * ====================================================
1043295facSStefan Farfeleder */
1143295facSStefan Farfeleder
1243295facSStefan Farfeleder /*
1343295facSStefan Farfeleder * ceill(x)
1443295facSStefan Farfeleder * Return x rounded toward -inf to integral value
1543295facSStefan Farfeleder * Method:
1643295facSStefan Farfeleder * Bit twiddling.
1743295facSStefan Farfeleder * Exception:
1843295facSStefan Farfeleder * Inexact flag raised if x not equal to ceill(x).
1943295facSStefan Farfeleder */
2043295facSStefan Farfeleder
2143295facSStefan Farfeleder #include <float.h>
2243295facSStefan Farfeleder #include <math.h>
2343295facSStefan Farfeleder #include <stdint.h>
2443295facSStefan Farfeleder
2543295facSStefan Farfeleder #include "fpmath.h"
2643295facSStefan Farfeleder
2743295facSStefan Farfeleder #ifdef LDBL_IMPLICIT_NBIT
2843295facSStefan Farfeleder #define MANH_SIZE (LDBL_MANH_SIZE + 1)
2943295facSStefan Farfeleder #define INC_MANH(u, c) do { \
3043295facSStefan Farfeleder uint64_t o = u.bits.manh; \
3143295facSStefan Farfeleder u.bits.manh += (c); \
3243295facSStefan Farfeleder if (u.bits.manh < o) \
3343295facSStefan Farfeleder u.bits.exp++; \
3443295facSStefan Farfeleder } while (0)
3543295facSStefan Farfeleder #else
3643295facSStefan Farfeleder #define MANH_SIZE LDBL_MANH_SIZE
3743295facSStefan Farfeleder #define INC_MANH(u, c) do { \
3843295facSStefan Farfeleder uint64_t o = u.bits.manh; \
3943295facSStefan Farfeleder u.bits.manh += (c); \
4043295facSStefan Farfeleder if (u.bits.manh < o) { \
4143295facSStefan Farfeleder u.bits.exp++; \
4243295facSStefan Farfeleder u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \
4343295facSStefan Farfeleder } \
4443295facSStefan Farfeleder } while (0)
4543295facSStefan Farfeleder #endif
4643295facSStefan Farfeleder
4766116c07SStefan Farfeleder static const long double huge = 1.0e300;
489eb30792SStefan Farfeleder
4943295facSStefan Farfeleder long double
ceill(long double x)5043295facSStefan Farfeleder ceill(long double x)
5143295facSStefan Farfeleder {
5243295facSStefan Farfeleder union IEEEl2bits u = { .e = x };
5343295facSStefan Farfeleder int e = u.bits.exp - LDBL_MAX_EXP + 1;
5443295facSStefan Farfeleder
5543295facSStefan Farfeleder if (e < MANH_SIZE - 1) {
5643295facSStefan Farfeleder if (e < 0) { /* raise inexact if x != 0 */
5766116c07SStefan Farfeleder if (huge + x > 0.0)
589eb30792SStefan Farfeleder if (u.bits.exp > 0 ||
599eb30792SStefan Farfeleder (u.bits.manh | u.bits.manl) != 0)
6065971872SBruce Evans u.e = u.bits.sign ? -0.0 : 1.0;
6143295facSStefan Farfeleder } else {
6243295facSStefan Farfeleder uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1);
6343295facSStefan Farfeleder if (((u.bits.manh & m) | u.bits.manl) == 0)
6443295facSStefan Farfeleder return (x); /* x is integral */
6543295facSStefan Farfeleder if (!u.bits.sign) {
6643295facSStefan Farfeleder #ifdef LDBL_IMPLICIT_NBIT
6743295facSStefan Farfeleder if (e == 0)
6843295facSStefan Farfeleder u.bits.exp++;
6943295facSStefan Farfeleder else
7043295facSStefan Farfeleder #endif
7143295facSStefan Farfeleder INC_MANH(u, 1llu << (MANH_SIZE - e - 1));
7243295facSStefan Farfeleder }
7366116c07SStefan Farfeleder if (huge + x > 0.0) { /* raise inexact flag */
7443295facSStefan Farfeleder u.bits.manh &= ~m;
7543295facSStefan Farfeleder u.bits.manl = 0;
7643295facSStefan Farfeleder }
779eb30792SStefan Farfeleder }
7843295facSStefan Farfeleder } else if (e < LDBL_MANT_DIG - 1) {
7943295facSStefan Farfeleder uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1);
8043295facSStefan Farfeleder if ((u.bits.manl & m) == 0)
8143295facSStefan Farfeleder return (x); /* x is integral */
8243295facSStefan Farfeleder if (!u.bits.sign) {
8343295facSStefan Farfeleder if (e == MANH_SIZE - 1)
8443295facSStefan Farfeleder INC_MANH(u, 1);
8543295facSStefan Farfeleder else {
8643295facSStefan Farfeleder uint64_t o = u.bits.manl;
8743295facSStefan Farfeleder u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1);
8843295facSStefan Farfeleder if (u.bits.manl < o) /* got a carry */
8943295facSStefan Farfeleder INC_MANH(u, 1);
9043295facSStefan Farfeleder }
9143295facSStefan Farfeleder }
9266116c07SStefan Farfeleder if (huge + x > 0.0) /* raise inexact flag */
9343295facSStefan Farfeleder u.bits.manl &= ~m;
9443295facSStefan Farfeleder }
9543295facSStefan Farfeleder return (u.e);
9643295facSStefan Farfeleder }
97