1 | /* Compute remainder and a congruent to the quotient. |
---|---|

2 | Copyright (C) 1997-2017 Free Software Foundation, Inc. |

3 | This file is part of the GNU C Library. |

4 | Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997 and |

5 | Jakub Jelinek <jj@ultra.linux.cz>, 1999. |

6 | |

7 | The GNU C Library is free software; you can redistribute it and/or |

8 | modify it under the terms of the GNU Lesser General Public |

9 | License as published by the Free Software Foundation; either |

10 | version 2.1 of the License, or (at your option) any later version. |

11 | |

12 | The GNU C Library is distributed in the hope that it will be useful, |

13 | but WITHOUT ANY WARRANTY; without even the implied warranty of |

14 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |

15 | Lesser General Public License for more details. |

16 | |

17 | You should have received a copy of the GNU Lesser General Public |

18 | License along with the GNU C Library; if not, write to the Free |

19 | Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA |

20 | 02111-1307 USA. */ |

21 | |

22 | #include "quadmath-imp.h" |

23 | |

24 | |

25 | static const __float128 zero = 0.0; |

26 | |

27 | |

28 | __float128 |

29 | remquoq (__float128 x, __float128 y, int *quo) |

30 | { |

31 | int64_t hx,hy; |

32 | uint64_t sx,lx,ly,qs; |

33 | int cquo; |

34 | |

35 | GET_FLT128_WORDS64 (hx, lx, x); |

36 | GET_FLT128_WORDS64 (hy, ly, y); |

37 | sx = hx & 0x8000000000000000ULL; |

38 | qs = sx ^ (hy & 0x8000000000000000ULL); |

39 | hy &= 0x7fffffffffffffffLL; |

40 | hx &= 0x7fffffffffffffffLL; |

41 | |

42 | /* Purge off exception values. */ |

43 | if ((hy | ly) == 0) |

44 | return (x * y) / (x * y); /* y = 0 */ |

45 | if ((hx >= 0x7fff000000000000LL) /* x not finite */ |

46 | || ((hy >= 0x7fff000000000000LL) /* y is NaN */ |

47 | && (((hy - 0x7fff000000000000LL) | ly) != 0))) |

48 | return (x * y) / (x * y); |

49 | |

50 | if (hy <= 0x7ffbffffffffffffLL) |

51 | x = fmodq (x, 8 * y); /* now x < 8y */ |

52 | |

53 | if (((hx - hy) | (lx - ly)) == 0) |

54 | { |

55 | *quo = qs ? -1 : 1; |

56 | return zero * x; |

57 | } |

58 | |

59 | x = fabsq (x); |

60 | y = fabsq (y); |

61 | cquo = 0; |

62 | |

63 | if (hy <= 0x7ffcffffffffffffLL && x >= 4 * y) |

64 | { |

65 | x -= 4 * y; |

66 | cquo += 4; |

67 | } |

68 | if (hy <= 0x7ffdffffffffffffLL && x >= 2 * y) |

69 | { |

70 | x -= 2 * y; |

71 | cquo += 2; |

72 | } |

73 | |

74 | if (hy < 0x0002000000000000LL) |

75 | { |

76 | if (x + x > y) |

77 | { |

78 | x -= y; |

79 | ++cquo; |

80 | if (x + x >= y) |

81 | { |

82 | x -= y; |

83 | ++cquo; |

84 | } |

85 | } |

86 | } |

87 | else |

88 | { |

89 | __float128 y_half = 0.5Q * y; |

90 | if (x > y_half) |

91 | { |

92 | x -= y; |

93 | ++cquo; |

94 | if (x >= y_half) |

95 | { |

96 | x -= y; |

97 | ++cquo; |

98 | } |

99 | } |

100 | } |

101 | |

102 | *quo = qs ? -cquo : cquo; |

103 | |

104 | /* Ensure correct sign of zero result in round-downward mode. */ |

105 | if (x == 0.0Q) |

106 | x = 0.0Q; |

107 | if (sx) |

108 | x = -x; |

109 | return x; |

110 | } |

111 |