Coverage Report

Created: 2024-11-21 06:47

/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
672k
  ((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
192k
{
43
192k
  unsigned nbits_in_r;
44
192k
  mp_limb_t r;
45
192k
  mp_size_t i;
46
47
192k
  if (bi <= nbits)
48
318
    {
49
318
      return p[0] & (((mp_limb_t) 1 << bi) - 1);
50
318
    }
51
192k
  else
52
192k
    {
53
192k
      bi -= nbits;      /* bit index of low bit to extract */
54
192k
      i = bi / GMP_NUMB_BITS;   /* word index of low bit to extract */
55
192k
      bi %= GMP_NUMB_BITS;   /* bit index in low word */
56
192k
      r = p[i] >> bi;     /* extract (low) bits */
57
192k
      nbits_in_r = GMP_NUMB_BITS - bi;  /* number of bits now in r */
58
192k
      if (nbits_in_r < nbits)    /* did we get enough bits? */
59
17.7k
  r += p[i + 1] << nbits_in_r; /* prepend bits from higher word */
60
192k
      return r & (((mp_limb_t) 1 << nbits) - 1);
61
192k
    }
62
192k
}
63
64
static inline unsigned
65
win_size (mp_bitcnt_t eb)
66
622
{
67
622
  unsigned k;
68
622
  static mp_bitcnt_t x[] = {7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0};
69
622
  ASSERT (eb > 1);
70
3.79k
  for (k = 0; eb > x[k++];)
71
3.17k
    ;
72
622
  return k;
73
622
}
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
622
{
85
622
  unsigned cnt;
86
622
  mp_bitcnt_t ebi;
87
622
  unsigned windowsize, this_windowsize;
88
622
  mp_limb_t expbits;
89
622
  mp_limb_t *pp;
90
622
  long i;
91
622
  int flipflop;
92
622
  TMP_DECL;
93
94
622
  ASSERT (en > 1 || (en == 1 && ep[0] > 1));
95
96
622
  TMP_MARK;
97
98
622
  MPN_SIZEINBASE_2EXP(ebi, ep, en, 1);
99
100
622
  windowsize = win_size (ebi);
101
622
  if (windowsize > 1)
102
607
    {
103
607
      mp_limb_t *this_pp, *last_pp;
104
607
      ASSERT (windowsize < ebi);
105
106
607
      pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)));
107
108
607
      this_pp = pp;
109
110
607
      MPN_COPY (this_pp, bp, n);
111
112
      /* Store b^2 in tp.  */
113
607
      mpn_sqrlo (tp, bp, n);
114
115
      /* Precompute odd powers of b and put them in the temporary area at pp.  */
116
607
      i = (1 << (windowsize - 1)) - 1;
117
607
      do
118
29.6k
  {
119
29.6k
    last_pp = this_pp;
120
29.6k
    this_pp += n;
121
29.6k
    mpn_mullo_n (this_pp, last_pp, tp, n);
122
29.6k
  } while (--i != 0);
123
124
607
      expbits = getbits (ep, ebi, windowsize);
125
607
      ebi -= windowsize;
126
127
      /* THINK: Should we initialise the case expbits % 4 == 0 with a mullo? */
128
607
      count_trailing_zeros (cnt, expbits);
129
607
      ebi += cnt;
130
607
      expbits >>= cnt;
131
132
607
      MPN_COPY (rp, pp + n * (expbits >> 1), n);
133
607
    }
134
15
  else
135
15
    {
136
15
      pp = tp + n;
137
15
      MPN_COPY (pp, bp, n);
138
15
      MPN_COPY (rp, bp, n);
139
15
      --ebi;
140
15
    }
141
142
622
  flipflop = 0;
143
144
622
  do
145
192k
    {
146
672k
      while (getbit (ep, ebi) == 0)
147
479k
  {
148
479k
    mpn_sqrlo (tp, rp, n);
149
479k
    MP_PTR_SWAP (rp, tp);
150
479k
    flipflop = ! flipflop;
151
479k
    if (--ebi == 0)
152
373
      goto done;
153
479k
  }
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
192k
      expbits = getbits (ep, ebi, windowsize);
159
192k
      this_windowsize = MIN (windowsize, ebi);
160
161
192k
      count_trailing_zeros (cnt, expbits);
162
192k
      this_windowsize -= cnt;
163
192k
      ebi -= this_windowsize;
164
192k
      expbits >>= cnt;
165
166
699k
      while (this_windowsize > 1)
167
507k
  {
168
507k
    mpn_sqrlo (tp, rp, n);
169
507k
    mpn_sqrlo (rp, tp, n);
170
507k
    this_windowsize -= 2;
171
507k
  }
172
173
192k
      if (this_windowsize != 0)
174
126k
  mpn_sqrlo (tp, rp, n);
175
65.7k
      else
176
65.7k
  {
177
65.7k
    MP_PTR_SWAP (rp, tp);
178
65.7k
    flipflop = ! flipflop;
179
65.7k
  }
180
181
192k
      mpn_mullo_n (rp, tp, pp + n * (expbits >> 1), n);
182
192k
    } while (ebi != 0);
183
184
622
 done:
185
622
  if (flipflop)
186
309
    MPN_COPY (tp, rp, n);
187
622
  TMP_FREE;
188
622
}