Line data Source code
1 : #include <tommath.h>
2 : #ifdef BN_S_MP_EXPTMOD_C
3 : /* LibTomMath, multiple-precision integer library -- Tom St Denis
4 : *
5 : * LibTomMath is a library that provides multiple-precision
6 : * integer arithmetic as well as number theoretic functionality.
7 : *
8 : * The library was designed directly after the MPI library by
9 : * Michael Fromberger but has been written from scratch with
10 : * additional optimizations in place.
11 : *
12 : * The library is free for all purposes without any express
13 : * guarantee it works.
14 : *
15 : * Tom St Denis, tomstdenis@gmail.com, http://libtom.org
16 : */
17 : #ifdef MP_LOW_MEM
18 : #define TAB_SIZE 32
19 : #else
20 : #define TAB_SIZE 256
21 : #endif
22 :
23 0 : int s_mp_exptmod (mp_int * G, mp_int * X, mp_int * P, mp_int * Y, int redmode)
24 : {
25 : mp_int M[TAB_SIZE], res, mu;
26 : mp_digit buf;
27 : int err, bitbuf, bitcpy, bitcnt, mode, digidx, x, y, winsize;
28 : int (*redux)(mp_int*,mp_int*,mp_int*);
29 :
30 : /* find window size */
31 0 : x = mp_count_bits (X);
32 0 : if (x <= 7) {
33 0 : winsize = 2;
34 0 : } else if (x <= 36) {
35 0 : winsize = 3;
36 0 : } else if (x <= 140) {
37 0 : winsize = 4;
38 0 : } else if (x <= 450) {
39 0 : winsize = 5;
40 0 : } else if (x <= 1303) {
41 0 : winsize = 6;
42 0 : } else if (x <= 3529) {
43 0 : winsize = 7;
44 : } else {
45 0 : winsize = 8;
46 : }
47 :
48 : #ifdef MP_LOW_MEM
49 : if (winsize > 5) {
50 : winsize = 5;
51 : }
52 : #endif
53 :
54 : /* init M array */
55 : /* init first cell */
56 0 : if ((err = mp_init(&M[1])) != MP_OKAY) {
57 0 : return err;
58 : }
59 :
60 : /* now init the second half of the array */
61 0 : for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
62 0 : if ((err = mp_init(&M[x])) != MP_OKAY) {
63 0 : for (y = 1<<(winsize-1); y < x; y++) {
64 0 : mp_clear (&M[y]);
65 : }
66 0 : mp_clear(&M[1]);
67 0 : return err;
68 : }
69 : }
70 :
71 : /* create mu, used for Barrett reduction */
72 0 : if ((err = mp_init (&mu)) != MP_OKAY) {
73 0 : goto LBL_M;
74 : }
75 :
76 0 : if (redmode == 0) {
77 0 : if ((err = mp_reduce_setup (&mu, P)) != MP_OKAY) {
78 0 : goto LBL_MU;
79 : }
80 0 : redux = mp_reduce;
81 : } else {
82 0 : if ((err = mp_reduce_2k_setup_l (P, &mu)) != MP_OKAY) {
83 0 : goto LBL_MU;
84 : }
85 0 : redux = mp_reduce_2k_l;
86 : }
87 :
88 : /* create M table
89 : *
90 : * The M table contains powers of the base,
91 : * e.g. M[x] = G**x mod P
92 : *
93 : * The first half of the table is not
94 : * computed though accept for M[0] and M[1]
95 : */
96 0 : if ((err = mp_mod (G, P, &M[1])) != MP_OKAY) {
97 0 : goto LBL_MU;
98 : }
99 :
100 : /* compute the value at M[1<<(winsize-1)] by squaring
101 : * M[1] (winsize-1) times
102 : */
103 0 : if ((err = mp_copy (&M[1], &M[1 << (winsize - 1)])) != MP_OKAY) {
104 0 : goto LBL_MU;
105 : }
106 :
107 0 : for (x = 0; x < (winsize - 1); x++) {
108 : /* square it */
109 0 : if ((err = mp_sqr (&M[1 << (winsize - 1)],
110 0 : &M[1 << (winsize - 1)])) != MP_OKAY) {
111 0 : goto LBL_MU;
112 : }
113 :
114 : /* reduce modulo P */
115 0 : if ((err = redux (&M[1 << (winsize - 1)], P, &mu)) != MP_OKAY) {
116 0 : goto LBL_MU;
117 : }
118 : }
119 :
120 : /* create upper table, that is M[x] = M[x-1] * M[1] (mod P)
121 : * for x = (2**(winsize - 1) + 1) to (2**winsize - 1)
122 : */
123 0 : for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
124 0 : if ((err = mp_mul (&M[x - 1], &M[1], &M[x])) != MP_OKAY) {
125 0 : goto LBL_MU;
126 : }
127 0 : if ((err = redux (&M[x], P, &mu)) != MP_OKAY) {
128 0 : goto LBL_MU;
129 : }
130 : }
131 :
132 : /* setup result */
133 0 : if ((err = mp_init (&res)) != MP_OKAY) {
134 0 : goto LBL_MU;
135 : }
136 0 : mp_set (&res, 1);
137 :
138 : /* set initial mode and bit cnt */
139 0 : mode = 0;
140 0 : bitcnt = 1;
141 0 : buf = 0;
142 0 : digidx = X->used - 1;
143 0 : bitcpy = 0;
144 0 : bitbuf = 0;
145 :
146 : for (;;) {
147 : /* grab next digit as required */
148 0 : if (--bitcnt == 0) {
149 : /* if digidx == -1 we are out of digits */
150 0 : if (digidx == -1) {
151 0 : break;
152 : }
153 : /* read next digit and reset the bitcnt */
154 0 : buf = X->dp[digidx--];
155 0 : bitcnt = (int) DIGIT_BIT;
156 : }
157 :
158 : /* grab the next msb from the exponent */
159 0 : y = (buf >> (mp_digit)(DIGIT_BIT - 1)) & 1;
160 0 : buf <<= (mp_digit)1;
161 :
162 : /* if the bit is zero and mode == 0 then we ignore it
163 : * These represent the leading zero bits before the first 1 bit
164 : * in the exponent. Technically this opt is not required but it
165 : * does lower the # of trivial squaring/reductions used
166 : */
167 0 : if (mode == 0 && y == 0) {
168 0 : continue;
169 : }
170 :
171 : /* if the bit is zero and mode == 1 then we square */
172 0 : if (mode == 1 && y == 0) {
173 0 : if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
174 0 : goto LBL_RES;
175 : }
176 0 : if ((err = redux (&res, P, &mu)) != MP_OKAY) {
177 0 : goto LBL_RES;
178 : }
179 0 : continue;
180 : }
181 :
182 : /* else we add it to the window */
183 0 : bitbuf |= (y << (winsize - ++bitcpy));
184 0 : mode = 2;
185 :
186 0 : if (bitcpy == winsize) {
187 : /* ok window is filled so square as required and multiply */
188 : /* square first */
189 0 : for (x = 0; x < winsize; x++) {
190 0 : if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
191 0 : goto LBL_RES;
192 : }
193 0 : if ((err = redux (&res, P, &mu)) != MP_OKAY) {
194 0 : goto LBL_RES;
195 : }
196 : }
197 :
198 : /* then multiply */
199 0 : if ((err = mp_mul (&res, &M[bitbuf], &res)) != MP_OKAY) {
200 0 : goto LBL_RES;
201 : }
202 0 : if ((err = redux (&res, P, &mu)) != MP_OKAY) {
203 0 : goto LBL_RES;
204 : }
205 :
206 : /* empty window and reset */
207 0 : bitcpy = 0;
208 0 : bitbuf = 0;
209 0 : mode = 1;
210 : }
211 : }
212 :
213 : /* if bits remain then square/multiply */
214 0 : if (mode == 2 && bitcpy > 0) {
215 : /* square then multiply if the bit is set */
216 0 : for (x = 0; x < bitcpy; x++) {
217 0 : if ((err = mp_sqr (&res, &res)) != MP_OKAY) {
218 0 : goto LBL_RES;
219 : }
220 0 : if ((err = redux (&res, P, &mu)) != MP_OKAY) {
221 0 : goto LBL_RES;
222 : }
223 :
224 0 : bitbuf <<= 1;
225 0 : if ((bitbuf & (1 << winsize)) != 0) {
226 : /* then multiply */
227 0 : if ((err = mp_mul (&res, &M[1], &res)) != MP_OKAY) {
228 0 : goto LBL_RES;
229 : }
230 0 : if ((err = redux (&res, P, &mu)) != MP_OKAY) {
231 0 : goto LBL_RES;
232 : }
233 : }
234 : }
235 : }
236 :
237 0 : mp_exch (&res, Y);
238 0 : err = MP_OKAY;
239 0 : LBL_RES:mp_clear (&res);
240 0 : LBL_MU:mp_clear (&mu);
241 0 : LBL_M:
242 0 : mp_clear(&M[1]);
243 0 : for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
244 0 : mp_clear (&M[x]);
245 : }
246 0 : return err;
247 : }
248 : #endif
249 :
250 : /* $Source: /cvs/libtom/libtommath/bn_s_mp_exptmod.c,v $ */
251 : /* $Revision: 1.5 $ */
252 : /* $Date: 2006/12/28 01:25:13 $ */
|
