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

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

Line | |
---|---|

1 | /* mpz_powm_ui(res,base,exp,mod) -- Set RES to (base**exp) mod MOD. |

2 | |

3 | Copyright 1991, 1993, 1994, 1996, 1997, 2000, 2001, 2002 Free Software |

4 | Foundation, Inc. |

5 | |

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

7 | |

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

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

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

11 | option) any later version. |

12 | |

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

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

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

16 | License for more details. |

17 | |

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

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

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

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

22 | |

23 | |

24 | #include "gmp.h" |

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

26 | #include "longlong.h" |

27 | |

28 | /* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and |

29 | t is defined by (tp,mn). */ |

30 | static void |

31 | reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn) |

32 | { |

33 | mp_ptr qp; |

34 | TMP_DECL (marker); |

35 | |

36 | TMP_MARK (marker); |

37 | qp = TMP_ALLOC_LIMBS (an - mn + 1); |

38 | |

39 | mpn_tdiv_qr (qp, tp, 0L, ap, an, mp, mn); |

40 | |

41 | TMP_FREE (marker); |

42 | } |

43 | |

44 | void |

45 | mpz_powm_ui (mpz_ptr r, mpz_srcptr b, unsigned long int el, mpz_srcptr m) |

46 | { |

47 | mp_ptr xp, tp, qp, mp, bp; |

48 | mp_size_t xn, tn, mn, bn; |

49 | int m_zero_cnt; |

50 | int c; |

51 | mp_limb_t e; |

52 | TMP_DECL (marker); |

53 | |

54 | mp = PTR(m); |

55 | mn = ABSIZ(m); |

56 | if (mn == 0) |

57 | DIVIDE_BY_ZERO; |

58 | |

59 | if (el == 0) |

60 | { |

61 | /* Exponent is zero, result is 1 mod MOD, i.e., 1 or 0 |

62 | depending on if MOD equals 1. */ |

63 | SIZ(r) = (mn == 1 && mp[0] == 1) ? 0 : 1; |

64 | PTR(r)[0] = 1; |

65 | return; |

66 | } |

67 | |

68 | TMP_MARK (marker); |

69 | |

70 | /* Normalize m (i.e. make its most significant bit set) as required by |

71 | division functions below. */ |

72 | count_leading_zeros (m_zero_cnt, mp[mn - 1]); |

73 | m_zero_cnt -= GMP_NAIL_BITS; |

74 | if (m_zero_cnt != 0) |

75 | { |

76 | mp_ptr new_mp = TMP_ALLOC_LIMBS (mn); |

77 | mpn_lshift (new_mp, mp, mn, m_zero_cnt); |

78 | mp = new_mp; |

79 | } |

80 | |

81 | bn = ABSIZ(b); |

82 | bp = PTR(b); |

83 | if (bn > mn) |

84 | { |

85 | /* Reduce possibly huge base. Use a function call to reduce, since we |

86 | don't want the quotient allocation to live until function return. */ |

87 | mp_ptr new_bp = TMP_ALLOC_LIMBS (mn); |

88 | reduce (new_bp, bp, bn, mp, mn); |

89 | bp = new_bp; |

90 | bn = mn; |

91 | /* Canonicalize the base, since we are potentially going to multiply with |

92 | it quite a few times. */ |

93 | MPN_NORMALIZE (bp, bn); |

94 | } |

95 | |

96 | if (bn == 0) |

97 | { |

98 | SIZ(r) = 0; |

99 | TMP_FREE (marker); |

100 | return; |

101 | } |

102 | |

103 | tp = TMP_ALLOC_LIMBS (2 * mn + 1); |

104 | xp = TMP_ALLOC_LIMBS (mn); |

105 | |

106 | qp = TMP_ALLOC_LIMBS (mn + 1); |

107 | |

108 | MPN_COPY (xp, bp, bn); |

109 | xn = bn; |

110 | |

111 | e = el; |

112 | count_leading_zeros (c, e); |

113 | e = (e << c) << 1; /* shift the exp bits to the left, lose msb */ |

114 | c = BITS_PER_MP_LIMB - 1 - c; |

115 | |

116 | /* Main loop. */ |

117 | |

118 | /* If m is already normalized (high bit of high limb set), and b is the |

119 | same size, but a bigger value, and e==1, then there's no modular |

120 | reductions done and we can end up with a result out of range at the |

121 | end. */ |

122 | if (c == 0) |

123 | { |

124 | if (xn == mn && mpn_cmp (xp, mp, mn) >= 0) |

125 | mpn_sub_n (xp, xp, mp, mn); |

126 | goto finishup; |

127 | } |

128 | |

129 | while (c != 0) |

130 | { |

131 | mpn_sqr_n (tp, xp, xn); |

132 | tn = 2 * xn; tn -= tp[tn - 1] == 0; |

133 | if (tn < mn) |

134 | { |

135 | MPN_COPY (xp, tp, tn); |

136 | xn = tn; |

137 | } |

138 | else |

139 | { |

140 | mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn); |

141 | xn = mn; |

142 | } |

143 | |

144 | if ((mp_limb_signed_t) e < 0) |

145 | { |

146 | mpn_mul (tp, xp, xn, bp, bn); |

147 | tn = xn + bn; tn -= tp[tn - 1] == 0; |

148 | if (tn < mn) |

149 | { |

150 | MPN_COPY (xp, tp, tn); |

151 | xn = tn; |

152 | } |

153 | else |

154 | { |

155 | mpn_tdiv_qr (qp, xp, 0L, tp, tn, mp, mn); |

156 | xn = mn; |

157 | } |

158 | } |

159 | e <<= 1; |

160 | c--; |

161 | } |

162 | |

163 | finishup: |

164 | /* We shifted m left m_zero_cnt steps. Adjust the result by reducing |

165 | it with the original MOD. */ |

166 | if (m_zero_cnt != 0) |

167 | { |

168 | mp_limb_t cy; |

169 | cy = mpn_lshift (tp, xp, xn, m_zero_cnt); |

170 | tp[xn] = cy; xn += cy != 0; |

171 | |

172 | if (xn < mn) |

173 | { |

174 | MPN_COPY (xp, tp, xn); |

175 | } |

176 | else |

177 | { |

178 | mpn_tdiv_qr (qp, xp, 0L, tp, xn, mp, mn); |

179 | xn = mn; |

180 | } |

181 | mpn_rshift (xp, xp, xn, m_zero_cnt); |

182 | } |

183 | MPN_NORMALIZE (xp, xn); |

184 | |

185 | if ((el & 1) != 0 && SIZ(b) < 0 && xn != 0) |

186 | { |

187 | mp = PTR(m); /* want original, unnormalized m */ |

188 | mpn_sub (xp, mp, mn, xp, xn); |

189 | xn = mn; |

190 | MPN_NORMALIZE (xp, xn); |

191 | } |

192 | MPZ_REALLOC (r, xn); |

193 | SIZ (r) = xn; |

194 | MPN_COPY (PTR(r), xp, xn); |

195 | |

196 | TMP_FREE (marker); |

197 | } |

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