#
source:
trunk/third/gmp/mpz/cong.c
@
18191

Revision 18191, 5.2 KB checked in by ghudson, 22 years ago (diff) |
---|

Line | |
---|---|

1 | /* mpz_congruent_p -- test congruence of two mpz's. |

2 | |

3 | Copyright 2001, 2002 Free Software Foundation, Inc. |

4 | |

5 | This file is part of the GNU MP Library. |

6 | |

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

8 | it under the terms of the GNU Lesser General Public License as published by |

9 | the Free Software Foundation; either version 2.1 of the License, or (at your |

10 | option) any later version. |

11 | |

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

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

14 | or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public |

15 | License for more details. |

16 | |

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

18 | along with the GNU MP Library; see the file COPYING.LIB. If not, write to |

19 | the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, |

20 | MA 02111-1307, USA. */ |

21 | |

22 | #include "gmp.h" |

23 | #include "gmp-impl.h" |

24 | #include "longlong.h" |

25 | |

26 | |

27 | /* For big divisors this code is only very slightly better than the user |

28 | doing a combination of mpz_sub and mpz_tdiv_r, but it's quite convenient, |

29 | and perhaps in the future can be improved, in similar ways to |

30 | mpn_divisible_p perhaps. |

31 | |

32 | The csize==1 / dsize==1 special case makes mpz_congruent_p as good as |

33 | mpz_congruent_ui_p on relevant operands, though such a combination |

34 | probably doesn't occur often. |

35 | |

36 | Alternatives: |

37 | |

38 | If c<d then it'd work to just form a%d and compare a and c (either as |

39 | a==c or a+c==d depending on the signs), but the saving from avoiding the |

40 | abs(a-c) calculation would be small compared to the division. |

41 | |

42 | Similarly if both a<d and c<d then it would work to just compare a and c |

43 | (a==c or a+c==d), but this isn't considered a particularly important case |

44 | and so isn't done for the moment. |

45 | |

46 | Low zero limbs on d could be stripped and the corresponding limbs of a |

47 | and c tested and skipped, but doing so would introduce a borrow when a |

48 | and c differ in sign and have non-zero skipped limbs. It doesn't seem |

49 | worth the complications to do this, since low zero limbs on d should |

50 | occur only rarely. */ |

51 | |

52 | int |

53 | mpz_congruent_p (mpz_srcptr a, mpz_srcptr c, mpz_srcptr d) |

54 | { |

55 | mp_size_t asize, csize, dsize, sign; |

56 | mp_srcptr ap, cp, dp; |

57 | mp_ptr xp; |

58 | mp_limb_t alow, clow, dlow, dmask, r; |

59 | int result; |

60 | TMP_DECL (marker); |

61 | |

62 | dsize = SIZ(d); |

63 | if (dsize == 0) |

64 | DIVIDE_BY_ZERO; |

65 | dsize = ABS(dsize); |

66 | dp = PTR(d); |

67 | |

68 | if (ABSIZ(a) < ABSIZ(c)) |

69 | MPZ_SRCPTR_SWAP (a, c); |

70 | |

71 | asize = SIZ(a); |

72 | csize = SIZ(c); |

73 | sign = (asize ^ csize); |

74 | |

75 | asize = ABS(asize); |

76 | ap = PTR(a); |

77 | |

78 | if (csize == 0) |

79 | return mpn_divisible_p (ap, asize, dp, dsize); |

80 | |

81 | csize = ABS(csize); |

82 | cp = PTR(c); |

83 | |

84 | alow = ap[0]; |

85 | clow = cp[0]; |

86 | dlow = dp[0]; |

87 | |

88 | /* Check a==c mod low zero bits of dlow. This might catch a few cases of |

89 | a!=c quickly, and it helps the csize==1 special cases below. */ |

90 | dmask = LOW_ZEROS_MASK (dlow) & GMP_NUMB_MASK; |

91 | alow = (sign >= 0 ? alow : -alow); |

92 | if (((alow-clow) & dmask) != 0) |

93 | return 0; |

94 | |

95 | if (csize == 1) |

96 | { |

97 | if (dsize == 1) |

98 | { |

99 | cong_1: |

100 | if (sign < 0) |

101 | NEG_MOD (clow, clow, dlow); |

102 | |

103 | if (BELOW_THRESHOLD (asize, MODEXACT_1_ODD_THRESHOLD)) |

104 | { |

105 | r = mpn_mod_1 (ap, asize, dlow); |

106 | if (clow < dlow) |

107 | return r == clow; |

108 | else |

109 | return r == (clow % dlow); |

110 | } |

111 | |

112 | if ((dlow & 1) == 0) |

113 | { |

114 | /* Strip low zero bits to get odd d required by modexact. If |

115 | d==e*2^n then a==c mod d if and only if both a==c mod e and |

116 | a==c mod 2^n, the latter having been done above. */ |

117 | unsigned twos; |

118 | count_trailing_zeros (twos, dlow); |

119 | dlow >>= twos; |

120 | } |

121 | |

122 | r = mpn_modexact_1c_odd (ap, asize, dlow, clow); |

123 | return r == 0 || r == dlow; |

124 | } |

125 | |

126 | /* dlow==0 is avoided since we don't want to bother handling extra low |

127 | zero bits if dsecond is even (would involve borrow if a,c differ in |

128 | sign and alow,clow!=0). */ |

129 | if (dsize == 2 && dlow != 0) |

130 | { |

131 | mp_limb_t dsecond = dp[1]; |

132 | |

133 | if (dsecond <= dmask) |

134 | { |

135 | unsigned twos; |

136 | count_trailing_zeros (twos, dlow); |

137 | dlow = (dlow >> twos) | (dsecond << (GMP_NUMB_BITS-twos)); |

138 | ASSERT_LIMB (dlow); |

139 | |

140 | /* dlow will be odd here, so the test for it even under cong_1 |

141 | is unnecessary, but the rest of that code is wanted. */ |

142 | goto cong_1; |

143 | } |

144 | } |

145 | } |

146 | |

147 | TMP_MARK (marker); |

148 | xp = TMP_ALLOC_LIMBS (asize+1); |

149 | |

150 | /* calculate abs(a-c) */ |

151 | if (sign >= 0) |

152 | { |

153 | /* same signs, subtract */ |

154 | if (asize > csize || mpn_cmp (ap, cp, asize) >= 0) |

155 | ASSERT_NOCARRY (mpn_sub (xp, ap, asize, cp, csize)); |

156 | else |

157 | ASSERT_NOCARRY (mpn_sub_n (xp, cp, ap, asize)); |

158 | MPN_NORMALIZE (xp, asize); |

159 | } |

160 | else |

161 | { |

162 | /* different signs, add */ |

163 | mp_limb_t carry; |

164 | carry = mpn_add (xp, ap, asize, cp, csize); |

165 | xp[asize] = carry; |

166 | asize += (carry != 0); |

167 | } |

168 | |

169 | result = mpn_divisible_p (xp, asize, dp, dsize); |

170 | |

171 | TMP_FREE (marker); |

172 | return result; |

173 | } |

**Note:**See TracBrowser for help on using the repository browser.