Coverage Report

Created: 2024-11-21 06:47

/src/libgmp/mpz/powm.c
Line
Count
Source (jump to first uncovered line)
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
5
Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2008, 2009,
6
2011, 2012, 2015, 2019 Free Software Foundation, Inc.
7
8
This file is part of the GNU MP Library.
9
10
The GNU MP Library is free software; you can redistribute it and/or modify
11
it 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
17
or
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
23
or both in parallel, as here.
24
25
The GNU MP Library is distributed in the hope that it will be useful, but
26
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
27
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
28
for more details.
29
30
You should have received copies of the GNU General Public License and the
31
GNU Lesser General Public License along with the GNU MP Library.  If not,
32
see 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
61
void
62
mpz_powm (mpz_ptr r, mpz_srcptr b, mpz_srcptr e, mpz_srcptr m)
63
2.99k
{
64
2.99k
  mp_size_t n, nodd, ncnt;
65
2.99k
  int cnt;
66
2.99k
  mp_ptr rp, tp;
67
2.99k
  mp_srcptr bp, ep, mp;
68
2.99k
  mp_size_t rn, bn, es, en, itch;
69
2.99k
  mpz_t new_b;      /* note: value lives long via 'b' */
70
2.99k
  TMP_DECL;
71
72
2.99k
  n = ABSIZ(m);
73
2.99k
  if (UNLIKELY (n == 0))
74
0
    DIVIDE_BY_ZERO;
75
76
2.99k
  mp = PTR(m);
77
78
2.99k
  TMP_MARK;
79
80
2.99k
  es = SIZ(e);
81
2.99k
  if (UNLIKELY (es <= 0))
82
0
    {
83
0
      if (es == 0)
84
0
  {
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
0
    SIZ(r) = n != 1 || mp[0] != 1;
88
0
    MPZ_NEWALLOC (r, 1)[0] = 1;
89
0
    TMP_FREE; /* we haven't really allocated anything here */
90
0
    return;
91
0
  }
92
0
#if HANDLE_NEGATIVE_EXPONENT
93
0
      MPZ_TMP_INIT (new_b, n + 1);
94
95
0
      if (UNLIKELY (! mpz_invert (new_b, b, m)))
96
0
  DIVIDE_BY_ZERO;
97
0
      b = new_b;
98
0
      es = -es;
99
#else
100
      DIVIDE_BY_ZERO;
101
#endif
102
0
    }
103
2.99k
  en = es;
104
105
2.99k
  bn = ABSIZ(b);
106
107
2.99k
  if (UNLIKELY (bn == 0))
108
0
    {
109
0
      SIZ(r) = 0;
110
0
      TMP_FREE;
111
0
      return;
112
0
    }
113
114
2.99k
  ep = PTR(e);
115
116
  /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case.  */
117
2.99k
  if (UNLIKELY (en == 1 && ep[0] == 1))
118
19
    {
119
19
      rp = TMP_ALLOC_LIMBS (n);
120
19
      bp = PTR(b);
121
19
      if (bn >= n)
122
12
  {
123
12
    mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1);
124
12
    mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n);
125
12
    rn = n;
126
12
    MPN_NORMALIZE (rp, rn);
127
128
12
    if (rn != 0 && SIZ(b) < 0)
129
0
      {
130
0
        mpn_sub (rp, mp, n, rp, rn);
131
0
        rn = n;
132
0
        MPN_NORMALIZE_NOT_ZERO (rp, rn);
133
0
      }
134
12
  }
135
7
      else
136
7
  {
137
7
    if (SIZ(b) < 0)
138
0
      {
139
0
        mpn_sub (rp, mp, n, bp, bn);
140
0
        rn = n;
141
0
        MPN_NORMALIZE_NOT_ZERO (rp, rn);
142
0
      }
143
7
    else
144
7
      {
145
7
        MPN_COPY (rp, bp, bn);
146
7
        rn = bn;
147
7
      }
148
7
  }
149
19
      goto ret;
150
19
    }
151
152
  /* Remove low zero limbs from M.  This loop will terminate for correctly
153
     represented mpz numbers.  */
154
2.97k
  ncnt = 0;
155
2.97k
  while (UNLIKELY (mp[0] == 0))
156
2.54k
    {
157
2.54k
      mp++;
158
2.54k
      ncnt++;
159
2.54k
    }
160
2.97k
  nodd = n - ncnt;
161
2.97k
  cnt = 0;
162
2.97k
  if (mp[0] % 2 == 0)
163
2.05k
    {
164
2.05k
      mp_ptr newmp = TMP_ALLOC_LIMBS (nodd);
165
2.05k
      count_trailing_zeros (cnt, mp[0]);
166
2.05k
      mpn_rshift (newmp, mp, nodd, cnt);
167
2.05k
      nodd -= newmp[nodd - 1] == 0;
168
2.05k
      mp = newmp;
169
2.05k
      ncnt++;
170
2.05k
    }
171
172
2.97k
  if (ncnt != 0)
173
2.06k
    {
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
2.06k
      mp_size_t n_largest_binvert = MAX (ncnt, nodd);
177
2.06k
      mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert);
178
2.06k
      itch = 3 * n + MAX (itch_binvert, 2 * n);
179
2.06k
    }
180
911
  else
181
911
    {
182
      /* We will call just mpn_powm.  */
183
911
      mp_size_t itch_binvert = mpn_binvert_itch (nodd);
184
911
      itch = n + MAX (itch_binvert, 2 * n);
185
911
    }
186
2.97k
  tp = TMP_ALLOC_LIMBS (itch);
187
188
2.97k
  rp = tp;  tp += n;
189
190
2.97k
  bp = PTR(b);
191
2.97k
  mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp);
192
193
2.97k
  rn = n;
194
195
2.97k
  if (ncnt != 0)
196
2.06k
    {
197
2.06k
      mp_ptr r2, xp, yp, odd_inv_2exp;
198
2.06k
      unsigned long t;
199
2.06k
      int bcnt;
200
201
2.06k
      if (bn < ncnt)
202
99
  {
203
99
    mp_ptr newbp = TMP_ALLOC_LIMBS (ncnt);
204
99
    MPN_COPY (newbp, bp, bn);
205
99
    MPN_ZERO (newbp + bn, ncnt - bn);
206
99
    bp = newbp;
207
99
  }
208
209
2.06k
      r2 = tp;
210
211
2.06k
      if (bp[0] % 2 == 0)
212
1.45k
  {
213
1.45k
    if (en > 1)
214
1.31k
      {
215
1.31k
        MPN_ZERO (r2, ncnt);
216
1.31k
        goto zero;
217
1.31k
      }
218
219
137
    ASSERT (en == 1);
220
137
    t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt;
221
222
    /* Count number of low zero bits in B, up to 3.  */
223
137
    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
137
    if (ep[0] * bcnt >= t)
227
130
      {
228
130
        MPN_ZERO (r2, ncnt);
229
130
        goto zero;
230
130
      }
231
137
  }
232
233
622
      mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt);
234
235
2.06k
    zero:
236
2.06k
      if (nodd < ncnt)
237
25
  {
238
25
    mp_ptr newmp = TMP_ALLOC_LIMBS (ncnt);
239
25
    MPN_COPY (newmp, mp, nodd);
240
25
    MPN_ZERO (newmp + nodd, ncnt - nodd);
241
25
    mp = newmp;
242
25
  }
243
244
2.06k
      odd_inv_2exp = tp + n;
245
2.06k
      mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n);
246
247
2.06k
      mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd);
248
249
2.06k
      xp = tp + 2 * n;
250
2.06k
      mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt);
251
252
2.06k
      if (cnt != 0)
253
2.05k
  xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1;
254
255
2.06k
      yp = tp;
256
2.06k
      if (ncnt > nodd)
257
25
  mpn_mul (yp, xp, ncnt, mp, nodd);
258
2.04k
      else
259
2.04k
  mpn_mul (yp, mp, nodd, xp, ncnt);
260
261
2.06k
      mpn_add (rp, yp, n, rp, nodd);
262
263
2.06k
      ASSERT (nodd + ncnt >= n);
264
2.06k
      ASSERT (nodd + ncnt <= n + 1);
265
2.06k
    }
266
267
2.97k
  MPN_NORMALIZE (rp, rn);
268
269
2.97k
  if ((ep[0] & 1) && SIZ(b) < 0 && rn != 0)
270
0
    {
271
0
      mpn_sub (rp, PTR(m), n, rp, rn);
272
0
      rn = n;
273
0
      MPN_NORMALIZE (rp, rn);
274
0
    }
275
276
2.99k
 ret:
277
2.99k
  MPZ_NEWALLOC (r, rn);
278
2.99k
  SIZ(r) = rn;
279
2.99k
  MPN_COPY (PTR(r), rp, rn);
280
281
2.99k
  TMP_FREE;
282
2.99k
}