Coverage Report

Created: 2024-11-21 07:03

/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
27.8k
{
64
27.8k
  mp_size_t n, nodd, ncnt;
65
27.8k
  int cnt;
66
27.8k
  mp_ptr rp, tp;
67
27.8k
  mp_srcptr bp, ep, mp;
68
27.8k
  mp_size_t rn, bn, es, en, itch;
69
27.8k
  mpz_t new_b;      /* note: value lives long via 'b' */
70
27.8k
  TMP_DECL;
71
72
27.8k
  n = ABSIZ(m);
73
27.8k
  if (UNLIKELY (n == 0))
74
0
    DIVIDE_BY_ZERO;
75
76
27.8k
  mp = PTR(m);
77
78
27.8k
  TMP_MARK;
79
80
27.8k
  es = SIZ(e);
81
27.8k
  if (UNLIKELY (es <= 0))
82
65
    {
83
65
      if (es == 0)
84
65
  {
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
65
    SIZ(r) = n != 1 || mp[0] != 1;
88
65
    MPZ_NEWALLOC (r, 1)[0] = 1;
89
65
    TMP_FREE; /* we haven't really allocated anything here */
90
65
    return;
91
65
  }
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
27.7k
  en = es;
104
105
27.7k
  bn = ABSIZ(b);
106
107
27.7k
  if (UNLIKELY (bn == 0))
108
12
    {
109
12
      SIZ(r) = 0;
110
12
      TMP_FREE;
111
12
      return;
112
12
    }
113
114
27.7k
  ep = PTR(e);
115
116
  /* Handle (b^1 mod m) early, since mpn_pow* do not handle that case.  */
117
27.7k
  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
1
  {
123
1
    mp_ptr qp = TMP_ALLOC_LIMBS (bn - n + 1);
124
1
    mpn_tdiv_qr (qp, rp, 0L, bp, bn, mp, n);
125
1
    rn = n;
126
1
    MPN_NORMALIZE (rp, rn);
127
128
1
    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
1
  }
135
18
      else
136
18
  {
137
18
    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
18
    else
144
18
      {
145
18
        MPN_COPY (rp, bp, bn);
146
18
        rn = bn;
147
18
      }
148
18
  }
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
27.7k
  ncnt = 0;
155
27.7k
  while (UNLIKELY (mp[0] == 0))
156
823
    {
157
823
      mp++;
158
823
      ncnt++;
159
823
    }
160
27.7k
  nodd = n - ncnt;
161
27.7k
  cnt = 0;
162
27.7k
  if (mp[0] % 2 == 0)
163
1.98k
    {
164
1.98k
      mp_ptr newmp = TMP_ALLOC_LIMBS (nodd);
165
1.98k
      count_trailing_zeros (cnt, mp[0]);
166
1.98k
      mpn_rshift (newmp, mp, nodd, cnt);
167
1.98k
      nodd -= newmp[nodd - 1] == 0;
168
1.98k
      mp = newmp;
169
1.98k
      ncnt++;
170
1.98k
    }
171
172
27.7k
  if (ncnt != 0)
173
2.00k
    {
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.00k
      mp_size_t n_largest_binvert = MAX (ncnt, nodd);
177
2.00k
      mp_size_t itch_binvert = mpn_binvert_itch (n_largest_binvert);
178
2.00k
      itch = 3 * n + MAX (itch_binvert, 2 * n);
179
2.00k
    }
180
25.7k
  else
181
25.7k
    {
182
      /* We will call just mpn_powm.  */
183
25.7k
      mp_size_t itch_binvert = mpn_binvert_itch (nodd);
184
25.7k
      itch = n + MAX (itch_binvert, 2 * n);
185
25.7k
    }
186
27.7k
  tp = TMP_ALLOC_LIMBS (itch);
187
188
27.7k
  rp = tp;  tp += n;
189
190
27.7k
  bp = PTR(b);
191
27.7k
  mpn_powm (rp, bp, bn, ep, en, mp, nodd, tp);
192
193
27.7k
  rn = n;
194
195
27.7k
  if (ncnt != 0)
196
2.00k
    {
197
2.00k
      mp_ptr r2, xp, yp, odd_inv_2exp;
198
2.00k
      unsigned long t;
199
2.00k
      int bcnt;
200
201
2.00k
      if (bn < ncnt)
202
83
  {
203
83
    mp_ptr newbp = TMP_ALLOC_LIMBS (ncnt);
204
83
    MPN_COPY (newbp, bp, bn);
205
83
    MPN_ZERO (newbp + bn, ncnt - bn);
206
83
    bp = newbp;
207
83
  }
208
209
2.00k
      r2 = tp;
210
211
2.00k
      if (bp[0] % 2 == 0)
212
1.26k
  {
213
1.26k
    if (en > 1)
214
913
      {
215
913
        MPN_ZERO (r2, ncnt);
216
913
        goto zero;
217
913
      }
218
219
356
    ASSERT (en == 1);
220
356
    t = (ncnt - (cnt != 0)) * GMP_NUMB_BITS + cnt;
221
222
    /* Count number of low zero bits in B, up to 3.  */
223
356
    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
356
    if (ep[0] * bcnt >= t)
227
276
      {
228
276
        MPN_ZERO (r2, ncnt);
229
276
        goto zero;
230
276
      }
231
356
  }
232
233
817
      mpn_powlo (r2, bp, ep, en, ncnt, tp + ncnt);
234
235
2.00k
    zero:
236
2.00k
      if (nodd < ncnt)
237
0
  {
238
0
    mp_ptr newmp = TMP_ALLOC_LIMBS (ncnt);
239
0
    MPN_COPY (newmp, mp, nodd);
240
0
    MPN_ZERO (newmp + nodd, ncnt - nodd);
241
0
    mp = newmp;
242
0
  }
243
244
2.00k
      odd_inv_2exp = tp + n;
245
2.00k
      mpn_binvert (odd_inv_2exp, mp, ncnt, tp + 2 * n);
246
247
2.00k
      mpn_sub (r2, r2, ncnt, rp, nodd > ncnt ? ncnt : nodd);
248
249
2.00k
      xp = tp + 2 * n;
250
2.00k
      mpn_mullo_n (xp, odd_inv_2exp, r2, ncnt);
251
252
2.00k
      if (cnt != 0)
253
1.98k
  xp[ncnt - 1] &= (CNST_LIMB(1) << cnt) - 1;
254
255
2.00k
      yp = tp;
256
2.00k
      if (ncnt > nodd)
257
0
  mpn_mul (yp, xp, ncnt, mp, nodd);
258
2.00k
      else
259
2.00k
  mpn_mul (yp, mp, nodd, xp, ncnt);
260
261
2.00k
      mpn_add (rp, yp, n, rp, nodd);
262
263
2.00k
      ASSERT (nodd + ncnt >= n);
264
2.00k
      ASSERT (nodd + ncnt <= n + 1);
265
2.00k
    }
266
267
27.7k
  MPN_NORMALIZE (rp, rn);
268
269
27.7k
  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
27.7k
 ret:
277
27.7k
  MPZ_NEWALLOC (r, rn);
278
27.7k
  SIZ(r) = rn;
279
27.7k
  MPN_COPY (PTR(r), rp, rn);
280
281
27.7k
  TMP_FREE;
282
27.7k
}