Coverage Report

Created: 2024-11-21 07:03

/src/libgmp/mpn/powlo.c
Line
Count
Source
1
/* mpn_powlo -- Compute R = U^E mod B^n, where B is the limb base.
2
3
Copyright 2007-2009, 2012, 2015, 2016, 2018, 2020 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 either:
10
11
  * the GNU Lesser General Public License as published by the Free
12
    Software Foundation; either version 3 of the License, or (at your
13
    option) any later version.
14
15
or
16
17
  * the GNU General Public License as published by the Free Software
18
    Foundation; either version 2 of the License, or (at your option) any
19
    later version.
20
21
or both in parallel, as here.
22
23
The GNU MP Library is distributed in the hope that it will be useful, but
24
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
25
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
26
for more details.
27
28
You should have received copies of the GNU General Public License and the
29
GNU Lesser General Public License along with the GNU MP Library.  If not,
30
see https://www.gnu.org/licenses/.  */
31
32
33
#include "gmp-impl.h"
34
#include "longlong.h"
35
36
37
#define getbit(p,bi) \
38
333k
  ((p[(bi - 1) / GMP_LIMB_BITS] >> (bi - 1) % GMP_LIMB_BITS) & 1)
39
40
static inline mp_limb_t
41
getbits (const mp_limb_t *p, mp_bitcnt_t bi, unsigned nbits)
42
112k
{
43
112k
  unsigned nbits_in_r;
44
112k
  mp_limb_t r;
45
112k
  mp_size_t i;
46
47
112k
  if (bi <= nbits)
48
2.11k
    {
49
2.11k
      return p[0] & (((mp_limb_t) 1 << bi) - 1);
50
2.11k
    }
51
110k
  else
52
110k
    {
53
110k
      bi -= nbits;      /* bit index of low bit to extract */
54
110k
      i = bi / GMP_NUMB_BITS;   /* word index of low bit to extract */
55
110k
      bi %= GMP_NUMB_BITS;   /* bit index in low word */
56
110k
      r = p[i] >> bi;     /* extract (low) bits */
57
110k
      nbits_in_r = GMP_NUMB_BITS - bi;  /* number of bits now in r */
58
110k
      if (nbits_in_r < nbits)    /* did we get enough bits? */
59
9.82k
  r += p[i + 1] << nbits_in_r; /* prepend bits from higher word */
60
110k
      return r & (((mp_limb_t) 1 << nbits) - 1);
61
110k
    }
62
112k
}
63
64
static inline unsigned
65
win_size (mp_bitcnt_t eb)
66
3.04k
{
67
3.04k
  unsigned k;
68
3.04k
  static mp_bitcnt_t x[] = {7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0};
69
3.04k
  ASSERT (eb > 1);
70
6.19k
  for (k = 0; eb > x[k++];)
71
3.14k
    ;
72
3.04k
  return k;
73
3.04k
}
74
75
/* rp[n-1..0] = bp[n-1..0] ^ ep[en-1..0] mod B^n, B is the limb base.
76
   Requires that ep[en-1] is non-zero.
77
   Uses scratch space tp[3n-1..0], i.e., 3n words.  */
78
/* We only use n words in the scratch space, we should pass tp + n to
79
   mullo/sqrlo as a temporary area, it is needed. */
80
void
81
mpn_powlo (mp_ptr rp, mp_srcptr bp,
82
     mp_srcptr ep, mp_size_t en,
83
     mp_size_t n, mp_ptr tp)
84
3.04k
{
85
3.04k
  unsigned cnt;
86
3.04k
  mp_bitcnt_t ebi;
87
3.04k
  unsigned windowsize, this_windowsize;
88
3.04k
  mp_limb_t expbits;
89
3.04k
  mp_limb_t *pp;
90
3.04k
  long i;
91
3.04k
  int flipflop;
92
3.04k
  TMP_DECL;
93
94
3.04k
  ASSERT (en > 1 || (en == 1 && ep[0] > 1));
95
96
3.04k
  TMP_MARK;
97
98
3.04k
  MPN_SIZEINBASE_2EXP(ebi, ep, en, 1);
99
100
3.04k
  windowsize = win_size (ebi);
101
3.04k
  if (windowsize > 1)
102
1.18k
    {
103
1.18k
      mp_limb_t *this_pp, *last_pp;
104
1.18k
      ASSERT (windowsize < ebi);
105
106
1.18k
      pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)));
107
108
1.18k
      this_pp = pp;
109
110
1.18k
      MPN_COPY (this_pp, bp, n);
111
112
      /* Store b^2 in tp.  */
113
1.18k
      mpn_sqrlo (tp, bp, n);
114
115
      /* Precompute odd powers of b and put them in the temporary area at pp.  */
116
1.18k
      i = (1 << (windowsize - 1)) - 1;
117
1.18k
      do
118
18.1k
  {
119
18.1k
    last_pp = this_pp;
120
18.1k
    this_pp += n;
121
18.1k
    mpn_mullo_n (this_pp, last_pp, tp, n);
122
18.1k
  } while (--i != 0);
123
124
1.18k
      expbits = getbits (ep, ebi, windowsize);
125
1.18k
      ebi -= windowsize;
126
127
      /* THINK: Should we initialise the case expbits % 4 == 0 with a mullo? */
128
1.18k
      count_trailing_zeros (cnt, expbits);
129
1.18k
      ebi += cnt;
130
1.18k
      expbits >>= cnt;
131
132
1.18k
      MPN_COPY (rp, pp + n * (expbits >> 1), n);
133
1.18k
    }
134
1.86k
  else
135
1.86k
    {
136
1.86k
      pp = tp + n;
137
1.86k
      MPN_COPY (pp, bp, n);
138
1.86k
      MPN_COPY (rp, bp, n);
139
1.86k
      --ebi;
140
1.86k
    }
141
142
3.04k
  flipflop = 0;
143
144
3.04k
  do
145
112k
    {
146
333k
      while (getbit (ep, ebi) == 0)
147
222k
  {
148
222k
    mpn_sqrlo (tp, rp, n);
149
222k
    MP_PTR_SWAP (rp, tp);
150
222k
    flipflop = ! flipflop;
151
222k
    if (--ebi == 0)
152
1.18k
      goto done;
153
222k
  }
154
155
      /* The next bit of the exponent is 1.  Now extract the largest block of
156
   bits <= windowsize, and such that the least significant bit is 1.  */
157
158
111k
      expbits = getbits (ep, ebi, windowsize);
159
111k
      this_windowsize = MIN (windowsize, ebi);
160
161
111k
      count_trailing_zeros (cnt, expbits);
162
111k
      this_windowsize -= cnt;
163
111k
      ebi -= this_windowsize;
164
111k
      expbits >>= cnt;
165
166
405k
      while (this_windowsize > 1)
167
294k
  {
168
294k
    mpn_sqrlo (tp, rp, n);
169
294k
    mpn_sqrlo (rp, tp, n);
170
294k
    this_windowsize -= 2;
171
294k
  }
172
173
111k
      if (this_windowsize != 0)
174
55.6k
  mpn_sqrlo (tp, rp, n);
175
55.3k
      else
176
55.3k
  {
177
55.3k
    MP_PTR_SWAP (rp, tp);
178
55.3k
    flipflop = ! flipflop;
179
55.3k
  }
180
181
111k
      mpn_mullo_n (rp, tp, pp + n * (expbits >> 1), n);
182
111k
    } while (ebi != 0);
183
184
3.04k
 done:
185
3.04k
  if (flipflop)
186
1.61k
    MPN_COPY (tp, rp, n);
187
3.04k
  TMP_FREE;
188
3.04k
}