Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpz/powm_ui.c
Line
Count
Source (jump to first uncovered line)
1
/* mpz_powm_ui(res,base,exp,mod) -- Set R to (B^E) mod M.
2
3
   Contributed to the GNU project by Torbjörn Granlund.
4
5
Copyright 1991, 1993, 1994, 1996, 1997, 2000-2002, 2005, 2008, 2009,
6
2011-2013, 2015 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
/* This code is very old, and should be rewritten to current GMP standard.  It
40
   is slower than mpz_powm for large exponents, but also for small exponents
41
   when the mod argument is small.
42
43
   As an intermediate solution, we now deflect to mpz_powm for exponents >= 20.
44
*/
45
46
/*
47
  b ^ e mod m   res
48
  0   0     0    ?
49
  0   e     0    ?
50
  0   0     m    ?
51
  0   e     m    0
52
  b   0     0    ?
53
  b   e     0    ?
54
  b   0     m    1 mod m
55
  b   e     m    b^e mod m
56
*/
57
58
static void
59
mod (mp_ptr np, mp_size_t nn, mp_srcptr dp, mp_size_t dn, gmp_pi1_t *dinv, mp_ptr tp)
60
45.2k
{
61
45.2k
  mp_ptr qp = tp;
62
63
45.2k
  if (dn == 1)
64
2.56k
    {
65
2.56k
      np[0] = mpn_divrem_1 (qp, (mp_size_t) 0, np, nn, dp[0]);
66
2.56k
    }
67
42.6k
  else if (dn == 2)
68
3.22k
    {
69
3.22k
      mpn_div_qr_2n_pi1 (qp, np, np, nn, dp[1], dp[0], dinv->inv32);
70
3.22k
    }
71
39.4k
  else if (BELOW_THRESHOLD (dn, DC_DIV_QR_THRESHOLD) ||
72
39.4k
     BELOW_THRESHOLD (nn - dn, DC_DIV_QR_THRESHOLD))
73
39.0k
    {
74
39.0k
      mpn_sbpi1_div_qr (qp, np, nn, dp, dn, dinv->inv32);
75
39.0k
    }
76
427
  else if (BELOW_THRESHOLD (dn, MUPI_DIV_QR_THRESHOLD) ||   /* fast condition */
77
427
     BELOW_THRESHOLD (nn, 2 * MU_DIV_QR_THRESHOLD) || /* fast condition */
78
427
     (double) (2 * (MU_DIV_QR_THRESHOLD - MUPI_DIV_QR_THRESHOLD)) * dn /* slow... */
79
0
     + (double) MUPI_DIV_QR_THRESHOLD * nn > (double) dn * nn)    /* ...condition */
80
427
    {
81
427
      mpn_dcpi1_div_qr (qp, np, nn, dp, dn, dinv);
82
427
    }
83
0
  else
84
0
    {
85
      /* We need to allocate separate remainder area, since mpn_mu_div_qr does
86
   not handle overlap between the numerator and remainder areas.
87
   FIXME: Make it handle such overlap.  */
88
0
      mp_ptr rp, scratch;
89
0
      mp_size_t itch;
90
0
      TMP_DECL;
91
0
      TMP_MARK;
92
93
0
      itch = mpn_mu_div_qr_itch (nn, dn, 0);
94
0
      rp = TMP_BALLOC_LIMBS (dn);
95
0
      scratch = TMP_BALLOC_LIMBS (itch);
96
97
0
      mpn_mu_div_qr (qp, rp, np, nn, dp, dn, scratch);
98
0
      MPN_COPY (np, rp, dn);
99
100
0
      TMP_FREE;
101
0
    }
102
45.2k
}
103
104
/* Compute t = a mod m, a is defined by (ap,an), m is defined by (mp,mn), and
105
   t is defined by (tp,mn).  */
106
static void
107
reduce (mp_ptr tp, mp_srcptr ap, mp_size_t an, mp_srcptr mp, mp_size_t mn, gmp_pi1_t *dinv)
108
0
{
109
0
  mp_ptr rp, scratch;
110
0
  TMP_DECL;
111
0
  TMP_MARK;
112
113
0
  TMP_ALLOC_LIMBS_2 (rp, an, scratch, an - mn + 1);
114
0
  MPN_COPY (rp, ap, an);
115
0
  mod (rp, an, mp, mn, dinv, scratch);
116
0
  MPN_COPY (tp, rp, mn);
117
118
0
  TMP_FREE;
119
0
}
120
121
void
122
mpz_powm_ui (mpz_ptr r, mpz_srcptr b, unsigned long int el, mpz_srcptr m)
123
23.5k
{
124
23.5k
  if (el < 20)
125
23.4k
    {
126
23.4k
      mp_ptr xp, tp, mp, bp, scratch;
127
23.4k
      mp_size_t xn, tn, mn, bn;
128
23.4k
      int m_zero_cnt;
129
23.4k
      int c;
130
23.4k
      mp_limb_t e, m2;
131
23.4k
      gmp_pi1_t dinv;
132
23.4k
      TMP_DECL;
133
134
23.4k
      mp = PTR(m);
135
23.4k
      mn = ABSIZ(m);
136
23.4k
      if (UNLIKELY (mn == 0))
137
0
  DIVIDE_BY_ZERO;
138
139
23.4k
      if (el <= 1)
140
2
  {
141
2
    if (el == 1)
142
0
      {
143
0
        mpz_mod (r, b, m);
144
0
        return;
145
0
      }
146
    /* Exponent is zero, result is 1 mod M, i.e., 1 or 0 depending on if
147
       M equals 1.  */
148
2
    SIZ(r) = mn != 1 || mp[0] != 1;
149
2
    MPZ_NEWALLOC (r, 1)[0] = 1;
150
2
    return;
151
2
  }
152
153
23.4k
      TMP_MARK;
154
155
      /* Normalize m (i.e. make its most significant bit set) as required by
156
   division functions below.  */
157
23.4k
      count_leading_zeros (m_zero_cnt, mp[mn - 1]);
158
23.4k
      m_zero_cnt -= GMP_NAIL_BITS;
159
23.4k
      if (m_zero_cnt != 0)
160
21.7k
  {
161
21.7k
    mp_ptr new_mp = TMP_ALLOC_LIMBS (mn);
162
21.7k
    mpn_lshift (new_mp, mp, mn, m_zero_cnt);
163
21.7k
    mp = new_mp;
164
21.7k
  }
165
166
23.4k
      m2 = mn == 1 ? 0 : mp[mn - 2];
167
23.4k
      invert_pi1 (dinv, mp[mn - 1], m2);
168
169
23.4k
      bn = ABSIZ(b);
170
23.4k
      bp = PTR(b);
171
23.4k
      if (bn > mn)
172
0
  {
173
    /* Reduce possibly huge base.  Use a function call to reduce, since we
174
       don't want the quotient allocation to live until function return.  */
175
0
    mp_ptr new_bp = TMP_ALLOC_LIMBS (mn);
176
0
    reduce (new_bp, bp, bn, mp, mn, &dinv);
177
0
    bp = new_bp;
178
0
    bn = mn;
179
    /* Canonicalize the base, since we are potentially going to multiply with
180
       it quite a few times.  */
181
0
    MPN_NORMALIZE (bp, bn);
182
0
  }
183
184
23.4k
      if (bn == 0)
185
0
  {
186
0
    SIZ(r) = 0;
187
0
    TMP_FREE;
188
0
    return;
189
0
  }
190
191
23.4k
      TMP_ALLOC_LIMBS_3 (xp, mn, scratch, mn + 1, tp, 2 * mn + 1);
192
193
23.4k
      MPN_COPY (xp, bp, bn);
194
23.4k
      xn = bn;
195
196
23.4k
      e = el;
197
23.4k
      count_leading_zeros (c, e);
198
23.4k
      e = (e << c) << 1;    /* shift the exp bits to the left, lose msb */
199
23.4k
      c = GMP_LIMB_BITS - 1 - c;
200
201
23.4k
      ASSERT (c != 0); /* el > 1 */
202
23.4k
  {
203
    /* Main loop. */
204
23.4k
    do
205
23.4k
      {
206
23.4k
        mpn_sqr (tp, xp, xn);
207
23.4k
        tn = 2 * xn; tn -= tp[tn - 1] == 0;
208
23.4k
        if (tn < mn)
209
0
    {
210
0
      MPN_COPY (xp, tp, tn);
211
0
      xn = tn;
212
0
    }
213
23.4k
        else
214
23.4k
    {
215
23.4k
      mod (tp, tn, mp, mn, &dinv, scratch);
216
23.4k
      MPN_COPY (xp, tp, mn);
217
23.4k
      xn = mn;
218
23.4k
    }
219
220
23.4k
        if ((mp_limb_signed_t) e < 0)
221
0
    {
222
0
      mpn_mul (tp, xp, xn, bp, bn);
223
0
      tn = xn + bn; tn -= tp[tn - 1] == 0;
224
0
      if (tn < mn)
225
0
        {
226
0
          MPN_COPY (xp, tp, tn);
227
0
          xn = tn;
228
0
        }
229
0
      else
230
0
        {
231
0
          mod (tp, tn, mp, mn, &dinv, scratch);
232
0
          MPN_COPY (xp, tp, mn);
233
0
          xn = mn;
234
0
        }
235
0
    }
236
23.4k
        e <<= 1;
237
23.4k
        c--;
238
23.4k
      }
239
23.4k
    while (c != 0);
240
23.4k
  }
241
242
      /* We shifted m left m_zero_cnt steps.  Adjust the result by reducing it
243
   with the original M.  */
244
23.4k
      if (m_zero_cnt != 0)
245
21.7k
  {
246
21.7k
    mp_limb_t cy;
247
21.7k
    cy = mpn_lshift (tp, xp, xn, m_zero_cnt);
248
21.7k
    tp[xn] = cy; xn += cy != 0;
249
250
21.7k
    if (xn >= mn)
251
21.7k
      {
252
21.7k
        mod (tp, xn, mp, mn, &dinv, scratch);
253
21.7k
        xn = mn;
254
21.7k
      }
255
21.7k
    mpn_rshift (xp, tp, xn, m_zero_cnt);
256
21.7k
  }
257
23.4k
      MPN_NORMALIZE (xp, xn);
258
259
23.4k
      if ((el & 1) != 0 && SIZ(b) < 0 && xn != 0)
260
0
  {
261
0
    mp = PTR(m);      /* want original, unnormalized m */
262
0
    mpn_sub (xp, mp, mn, xp, xn);
263
0
    xn = mn;
264
0
    MPN_NORMALIZE (xp, xn);
265
0
  }
266
23.4k
      MPZ_NEWALLOC (r, xn);
267
23.4k
      SIZ (r) = xn;
268
23.4k
      MPN_COPY (PTR(r), xp, xn);
269
270
23.4k
      TMP_FREE;
271
23.4k
    }
272
29
  else
273
29
    {
274
      /* For large exponents, fake an mpz_t exponent and deflect to the more
275
   sophisticated mpz_powm.  */
276
29
      mpz_t e;
277
29
      mp_limb_t ep[LIMBS_PER_ULONG];
278
29
      MPZ_FAKE_UI (e, ep, el);
279
29
      mpz_powm (r, b, e, m);
280
29
    }
281
23.5k
}