1/* mpz_powm(res,base,exp,mod) -- Set R to (U^E) mod M.
2
3   Contributed to the GNU project by Torbjorn Granlund.
4
5Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2008, 2009,
62011, 2012, 2015, 2019 Free Software Foundation, Inc.
7
8This file is part of the GNU MP Library.
9
10The GNU MP Library is free software; you can redistribute it and/or modify
11it under the terms of either:
12
13  * the GNU Lesser General Public License as published by the Free
14    Software Foundation; either version 3 of the License, or (at your
15    option) any later version.
16
17or
18
19  * the GNU General Public License as published by the Free Software
20    Foundation; either version 2 of the License, or (at your option) any
21    later version.
22
23or both in parallel, as here.
24
25The GNU MP Library is distributed in the hope that it will be useful, but
26WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
28for more details.
29
30You should have received copies of the GNU General Public License and the
31GNU Lesser General Public License along with the GNU MP Library.  If not,
32see https://www.gnu.org/licenses/.  */
33
34
35#include "gmp-impl.h"
36#include "longlong.h"
37
38
39/* TODO
40
41 * Improve handling of buffers.  It is pretty ugly now.
42
43 * For even moduli, we compute a binvert of its odd part both here and in
44   mpn_powm.  How can we avoid this recomputation?
45*/
46
47/*
48  b ^ e mod m   res
49  0   0     0    ?
50  0   e     0    ?
51  0   0     m    ?
52  0   e     m    0
53  b   0     0    ?
54  b   e     0    ?
55  b   0     m    1 mod m
56  b   e     m    b^e mod m
57*/
58
59#define HANDLE_NEGATIVE_EXPONENT 1
60
61void
62mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m)
63{
64  mp_size_t n, nodd, ncnt;
65  int cnt;
66  mp_ptr rp, tp;
67  mp_srcptr bp, ep, mp;
68  mp_size_t rn, bn, es, en, itch;
69  mpz_t new_b;			/* note: value lives long via 'b' */
70  TMP_DECL;
71
72  n = ABSIZ(m);
73  if (UNLIKELY (n == 0))
74    DIVIDE_BY_ZERO;
75
76  mp = PTR(m);
77
78  TMP_MARK;
79
80  es = SIZ(e);
81  if (UNLIKELY (es <= 0))
82    {
83      if (es == 0)
84	{
85	  /* b^0 mod m,  b is anything and m is non-zero.
86	     Result is 1 mod m, i.e., 1 or 0 depending on if m = 1.  */
87	  SIZ(r) = n != 1 || mp[0] != 1;
88	  MPZ_NEWALLOC (r, 1)[0] = 1;
89	  TMP_FREE;	/* we haven't really allocated anything here */
90	  return;
91	}
92#if HANDLE_NEGATIVE_EXPONENT
93      MPZ_TMP_INIT (new_b, n + 1);
94
95      if (UNLIKELY (! mpz_invert (new_b, b, m)))
96	DIVIDE_BY_ZERO;
97      b = new_b;
98      es = -es;
99#else
100      DIVIDE_BY_ZERO;
101#endif
102    }
103  en = es;
104
105  bn = ABSIZ(b);
106
107  if (UNLIKELY (bn == 0))
108    {
109      SIZ(r) = 0;
110      TMP_FREE;
111      return;
112    }
113
114  ep = PTR(e);
115
116  /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case.  */
117  if (UNLIKELY (en == 1 && ep[0] == 1))
118    {
119      rp = TMP_ALLOC_LIMBS (n);
120      bp = PTR(b);
121      if (bn >= n)
122	{
123	  mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1);
124	  mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n);
125	  rn = n;
126	  MPN_NORMALIZE (rp, rn);
127
128	  if (rn != 0 && SIZ(b) < 0)
129	    {
130	      mpn_sub (rp, mp, n, rp, rn);
131	      rn = n;
132	      MPN_NORMALIZE_NOT_ZERO (rp, rn);
133	    }
134	}
135      else
136	{
137	  if (SIZ(b) < 0)
138	    {
139	      mpn_sub (rp, mp, n, bp, bn);
140	      rn = n;
141	      MPN_NORMALIZE_NOT_ZERO (rp, rn);
142	    }
143	  else
144	    {
145	      MPN_COPY (rp, bp, bn);
146	      rn = bn;
147	    }
148	}
149      goto ret;
150    }
151
152  /* Remove low zero limbs from M.  This loop will terminate for correctly
153     represented mpz numbers.  */
154  ncnt = 0;
155  while (UNLIKELY (mp[0] == 0))
156    {
157      mp++;
158      ncnt++;
159    }
160  nodd = n - ncnt;
161  cnt = 0;
162  if (mp[0] % 2 == 0)
163    {
164      mp_ptr newmp = TMP_ALLOC_LIMBS (nodd);
165      count_trailing_zeros (cnt, mp[0]);
166      mpn_rshift (newmp, mp, nodd, cnt);
167      nodd -= newmp[nodd - 1] == 0;
168      mp = newmp;
169      ncnt++;
170    }
171
172  if (ncnt != 0)
173    {
174      /* We will call both mpn_powm and mpn_powlo.  */
175      /* rp needs n, mpn_powlo needs 4n, the 2 mpn_binvert might need more */
176      mp_size_t n_largest_binvert = MAX (ncnt, nodd);
177      mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert);
178      itch = 3 * n + MAX (itch_binvert, 2 * n);
179    }
180  else
181    {
182      /* We will call just mpn_powm.  */
183      mp_size_t itch_binvert = mpn_binvert_itch (nodd);
184      itch = n + MAX (itch_binvert, 2 * n);
185    }
186  tp = TMP_ALLOC_LIMBS (itch);
187
188  rp = tp;  tp += n;
189
190  bp = PTR(b);
191  mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp);
192
193  rn = n;
194
195  if (ncnt != 0)
196    {
197      mp_ptr r2, xp, yp, odd_inv_2exp;
198      unsigned long t;
199      int bcnt;
200
201      if (bn < ncnt)
202	{
203	  mp_ptr newbp = TMP_ALLOC_LIMBS (ncnt);
204	  MPN_COPY (newbp, bp, bn);
205	  MPN_ZERO (newbp + bn, ncnt - bn);
206	  bp = newbp;
207	}
208
209      r2 = tp;
210
211      if (bp[0] % 2 == 0)
212	{
213	  if (en > 1)
214	    {
215	      MPN_ZERO (r2, ncnt);
216	      goto zero;
217	    }
218
219	  ASSERT (en == 1);
220	  t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt;
221
222	  /* Count number of low zero bits in B, up to 3.  */
223	  bcnt = (0x1213 >> ((bp[0] & 7) << 1)) & 0x3;
224	  /* Note that ep[0] * bcnt might overflow, but that just results
225	     in a missed optimization.  */
226	  if (ep[0] * bcnt >= t)
227	    {
228	      MPN_ZERO (r2, ncnt);
229	      goto zero;
230	    }
231	}
232
233      mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt);
234
235    zero:
236      if (nodd < ncnt)
237	{
238	  mp_ptr newmp = TMP_ALLOC_LIMBS (ncnt);
239	  MPN_COPY (newmp, mp, nodd);
240	  MPN_ZERO (newmp + nodd, ncnt - nodd);
241	  mp = newmp;
242	}
243
244      odd_inv_2exp = tp + n;
245      mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n);
246
247      mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd);
248
249      xp = tp + 2 * n;
250      mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt);
251
252      if (cnt != 0)
253	xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1;
254
255      yp = tp;
256      if (ncnt > nodd)
257	mpn_mul (yp, xp, ncnt, mp, nodd);
258      else
259	mpn_mul (yp, mp, nodd, xp, ncnt);
260
261      mpn_add (rp, yp, n, rp, nodd);
262
263      ASSERT (nodd + ncnt >= n);
264      ASSERT (nodd + ncnt <= n + 1);
265    }
266
267  MPN_NORMALIZE (rp, rn);
268
269  if ((ep[0] & 1) && SIZ(b) < 0 && rn != 0)
270    {
271      mpn_sub (rp, PTR(m), n, rp, rn);
272      rn = n;
273      MPN_NORMALIZE (rp, rn);
274    }
275
276 ret:
277  MPZ_NEWALLOC (r, rn);
278  SIZ(r) = rn;
279  MPN_COPY (PTR(r), rp, rn);
280
281  TMP_FREE;
282}
283