1/* mpn_powlo -- Compute R = U^E mod B^n, where B is the limb base.
2
3Copyright 2007, 2008, 2009 Free Software Foundation, Inc.
4
5This file is part of the GNU MP Library.
6
7The GNU MP Library is free software; you can redistribute it and/or modify
8it under the terms of the GNU Lesser General Public License as published by
9the Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12The GNU MP Library is distributed in the hope that it will be useful, but
13WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
15License for more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along with the GNU MP Library.  If not, see http://www.gnu.org/licenses/.  */
19
20
21#include "gmp.h"
22#include "gmp-impl.h"
23#include "longlong.h"
24
25
26#define getbit(p,bi) \
27  ((p[(bi - 1) / GMP_LIMB_BITS] >> (bi - 1) % GMP_LIMB_BITS) & 1)
28
29static inline mp_limb_t
30getbits (const mp_limb_t *p, mp_bitcnt_t bi, int nbits)
31{
32  int nbits_in_r;
33  mp_limb_t r;
34  mp_size_t i;
35
36  if (bi < nbits)
37    {
38      return p[0] & (((mp_limb_t) 1 << bi) - 1);
39    }
40  else
41    {
42      bi -= nbits;			/* bit index of low bit to extract */
43      i = bi / GMP_NUMB_BITS;		/* word index of low bit to extract */
44      bi %= GMP_NUMB_BITS;		/* bit index in low word */
45      r = p[i] >> bi;			/* extract (low) bits */
46      nbits_in_r = GMP_NUMB_BITS - bi;	/* number of bits now in r */
47      if (nbits_in_r < nbits)		/* did we get enough bits? */
48	r += p[i + 1] << nbits_in_r;	/* prepend bits from higher word */
49      return r & (((mp_limb_t ) 1 << nbits) - 1);
50    }
51}
52
53static inline int
54win_size (mp_bitcnt_t eb)
55{
56  int k;
57  static mp_bitcnt_t x[] = {1,7,25,81,241,673,1793,4609,11521,28161,~(mp_bitcnt_t)0};
58  for (k = 0; eb > x[k]; k++)
59    ;
60  return k;
61}
62
63/* rp[n-1..0] = bp[n-1..0] ^ ep[en-1..0] mod B^n, B is the limb base.
64   Requires that ep[en-1] is non-zero.
65   Uses scratch space tp[3n-1..0], i.e., 3n words.  */
66void
67mpn_powlo (mp_ptr rp, mp_srcptr bp,
68	   mp_srcptr ep, mp_size_t en,
69	   mp_size_t n, mp_ptr tp)
70{
71  int cnt;
72  mp_bitcnt_t ebi;
73  int windowsize, this_windowsize;
74  mp_limb_t expbits;
75  mp_limb_t *pp, *this_pp, *last_pp;
76  mp_limb_t *b2p;
77  long i;
78  TMP_DECL;
79
80  ASSERT (en > 1 || (en == 1 && ep[0] > 1));
81
82  TMP_MARK;
83
84  count_leading_zeros (cnt, ep[en - 1]);
85  ebi = (mp_bitcnt_t) en * GMP_LIMB_BITS - cnt;
86
87  windowsize = win_size (ebi);
88
89  pp = TMP_ALLOC_LIMBS ((n << (windowsize - 1)) + n); /* + n is for mullo ign part */
90
91  this_pp = pp;
92
93  MPN_COPY (this_pp, bp, n);
94
95  b2p = tp + 2*n;
96
97  /* Store b^2 in b2.  */
98  mpn_sqr (tp, bp, n);	/* FIXME: Use "mpn_sqrlo" */
99  MPN_COPY (b2p, tp, n);
100
101  /* Precompute odd powers of b and put them in the temporary area at pp.  */
102  for (i = (1 << (windowsize - 1)) - 1; i > 0; i--)
103    {
104      last_pp = this_pp;
105      this_pp += n;
106      mpn_mullo_n (this_pp, last_pp, b2p, n);
107    }
108
109  expbits = getbits (ep, ebi, windowsize);
110  if (ebi < windowsize)
111    ebi = 0;
112  else
113    ebi -= windowsize;
114
115  count_trailing_zeros (cnt, expbits);
116  ebi += cnt;
117  expbits >>= cnt;
118
119  MPN_COPY (rp, pp + n * (expbits >> 1), n);
120
121  while (ebi != 0)
122    {
123      while (getbit (ep, ebi) == 0)
124	{
125	  mpn_sqr (tp, rp, n);	/* FIXME: Use "mpn_sqrlo" */
126	  MPN_COPY (rp, tp, n);
127	  ebi--;
128	  if (ebi == 0)
129	    goto done;
130	}
131
132      /* The next bit of the exponent is 1.  Now extract the largest block of
133	 bits <= windowsize, and such that the least significant bit is 1.  */
134
135      expbits = getbits (ep, ebi, windowsize);
136      this_windowsize = windowsize;
137      if (ebi < windowsize)
138	{
139	  this_windowsize -= windowsize - ebi;
140	  ebi = 0;
141	}
142      else
143	ebi -= windowsize;
144
145      count_trailing_zeros (cnt, expbits);
146      this_windowsize -= cnt;
147      ebi += cnt;
148      expbits >>= cnt;
149
150      do
151	{
152	  mpn_sqr (tp, rp, n);
153	  MPN_COPY (rp, tp, n);
154	  this_windowsize--;
155	}
156      while (this_windowsize != 0);
157
158      mpn_mullo_n (tp, rp, pp + n * (expbits >> 1), n);
159      MPN_COPY (rp, tp, n);
160    }
161
162 done:
163  TMP_FREE;
164}
165