summaryrefslogtreecommitdiff log msg author committer range
path: root/src/math/floorl.c
blob: 3901b060b2823613a1757900683577b9810c325e (plain)
 ```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 ``` ``````/* origin: FreeBSD /usr/src/lib/msun/src/s_floorl.c */ /* * ==================================================== * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. * * Developed at SunPro, a Sun Microsystems, Inc. business. * Permission to use, copy, modify, and distribute this * software is freely granted, provided that this notice * is preserved. * ==================================================== */ /* * floorl(x) * Return x rounded toward -inf to integral value * Method: * Bit twiddling. * Exception: * Inexact flag raised if x not equal to floorl(x). */ #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 long double floorl(long double x) { return floor(x); } #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #ifdef LDBL_IMPLICIT_NBIT #define MANH_SIZE (LDBL_MANH_SIZE + 1) #define INC_MANH(u, c) do { \ uint64_t o = u.bits.manh; \ u.bits.manh += (c); \ if (u.bits.manh < o) \ u.bits.exp++; \ } while (0) #else #define MANH_SIZE LDBL_MANH_SIZE #define INC_MANH(u, c) do { \ uint64_t o = u.bits.manh; \ u.bits.manh += (c); \ if (u.bits.manh < o) { \ u.bits.exp++; \ u.bits.manh |= 1llu << (LDBL_MANH_SIZE - 1); \ } \ } while (0) #endif static const long double huge = 1.0e300; long double floorl(long double x) { union IEEEl2bits u = { .e = x }; int e = u.bits.exp - LDBL_MAX_EXP + 1; if (e < MANH_SIZE - 1) { if (e < 0) { /* raise inexact if x != 0 */ if (huge + x > 0.0) if (u.bits.exp > 0 || (u.bits.manh | u.bits.manl) != 0) u.e = u.bits.sign ? -1.0 : 0.0; } else { uint64_t m = ((1llu << MANH_SIZE) - 1) >> (e + 1); if (((u.bits.manh & m) | u.bits.manl) == 0) return x; /* x is integral */ if (u.bits.sign) { #ifdef LDBL_IMPLICIT_NBIT if (e == 0) u.bits.exp++; else #endif INC_MANH(u, 1llu << (MANH_SIZE - e - 1)); } /* raise inexact flag */ if (huge + x > 0.0) { u.bits.manh &= ~m; u.bits.manl = 0; } } } else if (e < LDBL_MANT_DIG - 1) { uint64_t m = (uint64_t)-1 >> (64 - LDBL_MANT_DIG + e + 1); if ((u.bits.manl & m) == 0) return x; /* x is integral */ if (u.bits.sign) { if (e == MANH_SIZE - 1) INC_MANH(u, 1); else { uint64_t o = u.bits.manl; u.bits.manl += 1llu << (LDBL_MANT_DIG - e - 1); if (u.bits.manl < o) /* got a carry */ INC_MANH(u, 1); } } /* raise inexact flag */ if (huge + x > 0.0) u.bits.manl &= ~m; } return u.e; } #endif ``````