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 | |
---|

12 | /* |
---|

13 | * Modified for ALMOS-MKH OS at UPMC, France, August 2018. (Alain Greiner) |
---|

14 | */ |
---|

15 | |
---|

16 | /* |
---|

17 | * ceil(x) |
---|

18 | * Return x rounded toward -inf to integral value |
---|

19 | * Method: |
---|

20 | * Bit twiddling. |
---|

21 | * Exception: |
---|

22 | * Inexact flag raised if x not equal to ceil(x). |
---|

23 | */ |
---|

24 | |
---|

25 | #include "math.h" |
---|

26 | #include "math_private.h" |
---|

27 | |
---|

28 | static const double huge = 1.0e300; |
---|

29 | |
---|

30 | double ceil(double x) |
---|

31 | { |
---|

32 | int32_t i0, i1, j0; |
---|

33 | uint32_t i, j; |
---|

34 | EXTRACT_WORDS(i0, i1, x); |
---|

35 | j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; |
---|

36 | if (j0 < 20) { |
---|

37 | if (j0 < 0) { /* raise inexact if x != 0 */ |
---|

38 | if (huge + x > 0.0) {/* return 0*sign(x) if |x|<1 */ |
---|

39 | if (i0 < 0) { |
---|

40 | i0 = 0x80000000; |
---|

41 | i1 = 0; |
---|

42 | } |
---|

43 | else if ((i0 | i1) != 0) { |
---|

44 | i0 = 0x3ff00000; |
---|

45 | i1 = 0; |
---|

46 | } |
---|

47 | } |
---|

48 | } |
---|

49 | else { |
---|

50 | i = (0x000fffff) >> j0; |
---|

51 | if (((i0 & i) | i1) == 0) { |
---|

52 | return x; /* x is integral */ |
---|

53 | } |
---|

54 | if (huge + x > 0.0) { /* raise inexact flag */ |
---|

55 | if (i0 > 0) { |
---|

56 | i0 += (0x00100000) >> j0; |
---|

57 | } |
---|

58 | i0 &= (~i); |
---|

59 | i1 = 0; |
---|

60 | } |
---|

61 | } |
---|

62 | } |
---|

63 | else if (j0 > 51) { |
---|

64 | if (j0 == 0x400) { |
---|

65 | return x + x; /* inf or NaN */ |
---|

66 | } |
---|

67 | else { |
---|

68 | return x; /* x is integral */ |
---|

69 | } |
---|

70 | } |
---|

71 | else { |
---|

72 | i = ((uint32_t) (0xffffffff)) >> (j0 - 20); |
---|

73 | if ((i1 & i) == 0) { |
---|

74 | return x; /* x is integral */ |
---|

75 | } |
---|

76 | if (huge + x > 0.0) { /* raise inexact flag */ |
---|

77 | if (i0 > 0) { |
---|

78 | if (j0 == 20) { |
---|

79 | i0 += 1; |
---|

80 | } |
---|

81 | else { |
---|

82 | j = i1 + (1 << (52 - j0)); |
---|

83 | if (j < i1) { |
---|

84 | i0 += 1; /* got a carry */ |
---|

85 | } |
---|

86 | i1 = j; |
---|

87 | } |
---|

88 | } |
---|

89 | i1 &= (~i); |
---|

90 | } |
---|

91 | } |
---|

92 | INSERT_WORDS(x, i0, i1); |
---|

93 | return x; |
---|

94 | } |
---|

95 | |
---|