modular arithmetics
[openssl.git] / crypto / bn / bn_exp.c
1 /* crypto/bn/bn_exp.c */
2 /* Copyright (C) 1995-1998 Eric Young (eay@cryptsoft.com)
3  * All rights reserved.
4  *
5  * This package is an SSL implementation written
6  * by Eric Young (eay@cryptsoft.com).
7  * The implementation was written so as to conform with Netscapes SSL.
8  * 
9  * This library is free for commercial and non-commercial use as long as
10  * the following conditions are aheared to.  The following conditions
11  * apply to all code found in this distribution, be it the RC4, RSA,
12  * lhash, DES, etc., code; not just the SSL code.  The SSL documentation
13  * included with this distribution is covered by the same copyright terms
14  * except that the holder is Tim Hudson (tjh@cryptsoft.com).
15  * 
16  * Copyright remains Eric Young's, and as such any Copyright notices in
17  * the code are not to be removed.
18  * If this package is used in a product, Eric Young should be given attribution
19  * as the author of the parts of the library used.
20  * This can be in the form of a textual message at program startup or
21  * in documentation (online or textual) provided with the package.
22  * 
23  * Redistribution and use in source and binary forms, with or without
24  * modification, are permitted provided that the following conditions
25  * are met:
26  * 1. Redistributions of source code must retain the copyright
27  *    notice, this list of conditions and the following disclaimer.
28  * 2. Redistributions in binary form must reproduce the above copyright
29  *    notice, this list of conditions and the following disclaimer in the
30  *    documentation and/or other materials provided with the distribution.
31  * 3. All advertising materials mentioning features or use of this software
32  *    must display the following acknowledgement:
33  *    "This product includes cryptographic software written by
34  *     Eric Young (eay@cryptsoft.com)"
35  *    The word 'cryptographic' can be left out if the rouines from the library
36  *    being used are not cryptographic related :-).
37  * 4. If you include any Windows specific code (or a derivative thereof) from 
38  *    the apps directory (application code) you must include an acknowledgement:
39  *    "This product includes software written by Tim Hudson (tjh@cryptsoft.com)"
40  * 
41  * THIS SOFTWARE IS PROVIDED BY ERIC YOUNG ``AS IS'' AND
42  * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
43  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
44  * ARE DISCLAIMED.  IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
45  * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
46  * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
47  * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
48  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
49  * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
50  * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
51  * SUCH DAMAGE.
52  * 
53  * The licence and distribution terms for any publically available version or
54  * derivative of this code cannot be changed.  i.e. this code cannot simply be
55  * copied and put under another distribution licence
56  * [including the GNU Public Licence.]
57  */
58 /* ====================================================================
59  * Copyright (c) 1998-2000 The OpenSSL Project.  All rights reserved.
60  *
61  * Redistribution and use in source and binary forms, with or without
62  * modification, are permitted provided that the following conditions
63  * are met:
64  *
65  * 1. Redistributions of source code must retain the above copyright
66  *    notice, this list of conditions and the following disclaimer. 
67  *
68  * 2. Redistributions in binary form must reproduce the above copyright
69  *    notice, this list of conditions and the following disclaimer in
70  *    the documentation and/or other materials provided with the
71  *    distribution.
72  *
73  * 3. All advertising materials mentioning features or use of this
74  *    software must display the following acknowledgment:
75  *    "This product includes software developed by the OpenSSL Project
76  *    for use in the OpenSSL Toolkit. (http://www.openssl.org/)"
77  *
78  * 4. The names "OpenSSL Toolkit" and "OpenSSL Project" must not be used to
79  *    endorse or promote products derived from this software without
80  *    prior written permission. For written permission, please contact
81  *    openssl-core@openssl.org.
82  *
83  * 5. Products derived from this software may not be called "OpenSSL"
84  *    nor may "OpenSSL" appear in their names without prior written
85  *    permission of the OpenSSL Project.
86  *
87  * 6. Redistributions of any form whatsoever must retain the following
88  *    acknowledgment:
89  *    "This product includes software developed by the OpenSSL Project
90  *    for use in the OpenSSL Toolkit (http://www.openssl.org/)"
91  *
92  * THIS SOFTWARE IS PROVIDED BY THE OpenSSL PROJECT ``AS IS'' AND ANY
93  * EXPRESSED OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
94  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
95  * PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE OpenSSL PROJECT OR
96  * ITS CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
97  * SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
98  * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
99  * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
100  * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
101  * STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
102  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
103  * OF THE POSSIBILITY OF SUCH DAMAGE.
104  * ====================================================================
105  *
106  * This product includes cryptographic software written by Eric Young
107  * (eay@cryptsoft.com).  This product includes software written by Tim
108  * Hudson (tjh@cryptsoft.com).
109  *
110  */
111
112
113 #include "cryptlib.h"
114 #include "bn_lcl.h"
115
116 #define TABLE_SIZE      32
117
118 /* this one works - simple but works */
119 int BN_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, BN_CTX *ctx)
120         {
121         int i,bits,ret=0;
122         BIGNUM *v,*rr;
123
124         BN_CTX_start(ctx);
125         if ((r == a) || (r == p))
126                 rr = BN_CTX_get(ctx);
127         else
128                 rr = r;
129         if ((v = BN_CTX_get(ctx)) == NULL) goto err;
130
131         if (BN_copy(v,a) == NULL) goto err;
132         bits=BN_num_bits(p);
133
134         if (BN_is_odd(p))
135                 { if (BN_copy(rr,a) == NULL) goto err; }
136         else    { if (!BN_one(rr)) goto err; }
137
138         for (i=1; i<bits; i++)
139                 {
140                 if (!BN_sqr(v,v,ctx)) goto err;
141                 if (BN_is_bit_set(p,i))
142                         {
143                         if (!BN_mul(rr,rr,v,ctx)) goto err;
144                         }
145                 }
146         ret=1;
147 err:
148         if (r != rr) BN_copy(r,rr);
149         BN_CTX_end(ctx);
150         return(ret);
151         }
152
153
154 int BN_mod_exp(BIGNUM *r, const BIGNUM *a, const BIGNUM *p, const BIGNUM *m,
155                BN_CTX *ctx)
156         {
157         int ret;
158
159         bn_check_top(a);
160         bn_check_top(p);
161         bn_check_top(m);
162
163         /* For even modulus  m = 2^k*m_odd,  it might make sense to compute
164          * a^p mod m_odd  and  a^p mod 2^k  separately (with Montgomery
165          * exponentiation for the odd part), using appropriate exponent
166          * reductions, and combine the results using the CRT.
167          *
168          * For now, we use Montgomery only if the modulus is odd; otherwise,
169          * exponentiation using the reciprocal-based quick remaindering
170          * algorithm is used.
171          *
172          * (For computations  a^p mod m  where  a, p, m  are of the same
173          * length, BN_mod_exp_recp takes roughly 50 .. 70 % the time
174          * required by the standard algorithm, and BN_mod_exp takes
175          * about 33 .. 40 % of it.
176          * [Timings obtained with expspeed.c on a AMD K6-2 platform under Linux,
177          * with various OpenSSL debugging macros defined.  YMMV.])
178          */
179
180 #define MONT_MUL_MOD
181 #define RECP_MUL_MOD
182
183 #ifdef MONT_MUL_MOD
184         /* I have finally been able to take out this pre-condition of
185          * the top bit being set.  It was caused by an error in BN_div
186          * with negatives.  There was also another problem when for a^b%m
187          * a >= m.  eay 07-May-97 */
188 /*      if ((m->d[m->top-1]&BN_TBIT) && BN_is_odd(m)) */
189
190         if (BN_is_odd(m))
191                 {
192                 if (a->top == 1 && !a->neg)
193                         {
194                         BN_ULONG A = a->d[0];
195                         ret=BN_mod_exp_mont_word(r,A,p,m,ctx,NULL);
196                         }
197                 else
198                         ret=BN_mod_exp_mont(r,a,p,m,ctx,NULL);
199                 }
200         else
201 #endif
202 #ifdef RECP_MUL_MOD
203                 { ret=BN_mod_exp_recp(r,a,p,m,ctx); }
204 #else
205                 { ret=BN_mod_exp_simple(r,a,p,m,ctx); }
206 #endif
207
208         return(ret);
209         }
210
211
212 int BN_mod_exp_recp(BIGNUM *r, const BIGNUM *a, const BIGNUM *p,
213                     const BIGNUM *m, BN_CTX *ctx)
214         {
215         int i,j,bits,ret=0,wstart,wend,window,wvalue;
216         int start=1,ts=0;
217         BIGNUM *aa;
218         BIGNUM val[TABLE_SIZE];
219         BN_RECP_CTX recp;
220
221         bits=BN_num_bits(p);
222
223         if (bits == 0)
224                 {
225                 BN_one(r);
226                 return(1);
227                 }
228
229         BN_CTX_start(ctx);
230         if ((aa = BN_CTX_get(ctx)) == NULL) goto err;
231
232         BN_RECP_CTX_init(&recp);
233         if (BN_RECP_CTX_set(&recp,m,ctx) <= 0) goto err;
234
235         BN_init(&(val[0]));
236         ts=1;
237
238         if (!BN_nnmod(&(val[0]),a,m,ctx)) goto err;             /* 1 */
239
240         window = BN_window_bits_for_exponent_size(bits);
241         if (window > 1)
242                 {
243                 if (!BN_mod_mul_reciprocal(aa,&(val[0]),&(val[0]),&recp,ctx))
244                         goto err;                               /* 2 */
245                 j=1<<(window-1);
246                 for (i=1; i<j; i++)
247                         {
248                         BN_init(&val[i]);
249                         if (!BN_mod_mul_reciprocal(&(val[i]),&(val[i-1]),aa,&recp,ctx))
250                                 goto err;
251                         }
252                 ts=i;
253                 }
254                 
255         start=1;        /* This is used to avoid multiplication etc
256                          * when there is only the value '1' in the
257                          * buffer. */
258         wvalue=0;       /* The 'value' of the window */
259         wstart=bits-1;  /* The top bit of the window */
260         wend=0;         /* The bottom bit of the window */
261
262         if (!BN_one(r)) goto err;
263
264         for (;;)
265                 {
266                 if (BN_is_bit_set(p,wstart) == 0)
267                         {
268                         if (!start)
269                                 if (!BN_mod_mul_reciprocal(r,r,r,&recp,ctx))
270                                 goto err;
271                         if (wstart == 0) break;
272                         wstart--;
273                         continue;
274                         }
275                 /* We now have wstart on a 'set' bit, we now need to work out
276                  * how bit a window to do.  To do this we need to scan
277                  * forward until the last set bit before the end of the
278                  * window */
279                 j=wstart;
280                 wvalue=1;
281                 wend=0;
282                 for (i=1; i<window; i++)
283                         {
284                         if (wstart-i < 0) break;
285                         if (BN_is_bit_set(p,wstart-i))
286                                 {
287                                 wvalue<<=(i-wend);
288                                 wvalue|=1;
289                                 wend=i;
290                                 }
291                         }
292
293                 /* wend is the size of the current window */
294                 j=wend+1;
295                 /* add the 'bytes above' */
296                 if (!start)
297                         for (i=0; i<j; i++)
298                                 {
299                                 if (!BN_mod_mul_reciprocal(r,r,r,&recp,ctx))
300                                         goto err;
301                                 }
302                 
303                 /* wvalue will be an odd number < 2^window */
304                 if (!BN_mod_mul_reciprocal(r,r,&(val[wvalue>>1]),&recp,ctx))
305                         goto err;
306
307                 /* move the 'window' down further */
308                 wstart-=wend+1;
309                 wvalue=0;
310                 start=0;
311                 if (wstart < 0) break;
312                 }
313         ret=1;
314 err:
315         BN_CTX_end(ctx);
316         for (i=0; i<ts; i++)
317                 BN_clear_free(&(val[i]));
318         BN_RECP_CTX_free(&recp);
319         return(ret);
320         }
321
322
323 int BN_mod_exp_mont(BIGNUM *rr, const BIGNUM *a, const BIGNUM *p,
324                     const BIGNUM *m, BN_CTX *ctx, BN_MONT_CTX *in_mont)
325         {
326         int i,j,bits,ret=0,wstart,wend,window,wvalue;
327         int start=1,ts=0;
328         BIGNUM *d,*r;
329         const BIGNUM *aa;
330         BIGNUM val[TABLE_SIZE];
331         BN_MONT_CTX *mont=NULL;
332
333         bn_check_top(a);
334         bn_check_top(p);
335         bn_check_top(m);
336
337         if (!(m->d[0] & 1))
338                 {
339                 BNerr(BN_F_BN_MOD_EXP_MONT,BN_R_CALLED_WITH_EVEN_MODULUS);
340                 return(0);
341                 }
342         bits=BN_num_bits(p);
343         if (bits == 0)
344                 {
345                 BN_one(rr);
346                 return(1);
347                 }
348         BN_CTX_start(ctx);
349         d = BN_CTX_get(ctx);
350         r = BN_CTX_get(ctx);
351         if (d == NULL || r == NULL) goto err;
352
353         /* If this is not done, things will break in the montgomery
354          * part */
355
356         if (in_mont != NULL)
357                 mont=in_mont;
358         else
359                 {
360                 if ((mont=BN_MONT_CTX_new()) == NULL) goto err;
361                 if (!BN_MONT_CTX_set(mont,m,ctx)) goto err;
362                 }
363
364         BN_init(&val[0]);
365         ts=1;
366         if (!a->neg && BN_ucmp(a,m) >= 0)
367                 {
368                 if (!BN_nnmod(&(val[0]),a,m,ctx))
369                         goto err;
370                 aa= &(val[0]);
371                 }
372         else
373                 aa=a;
374         if (!BN_to_montgomery(&(val[0]),aa,mont,ctx)) goto err; /* 1 */
375
376         window = BN_window_bits_for_exponent_size(bits);
377         if (window > 1)
378                 {
379                 if (!BN_mod_mul_montgomery(d,&(val[0]),&(val[0]),mont,ctx)) goto err; /* 2 */
380                 j=1<<(window-1);
381                 for (i=1; i<j; i++)
382                         {
383                         BN_init(&(val[i]));
384                         if (!BN_mod_mul_montgomery(&(val[i]),&(val[i-1]),d,mont,ctx))
385                                 goto err;
386                         }
387                 ts=i;
388                 }
389
390         start=1;        /* This is used to avoid multiplication etc
391                          * when there is only the value '1' in the
392                          * buffer. */
393         wvalue=0;       /* The 'value' of the window */
394         wstart=bits-1;  /* The top bit of the window */
395         wend=0;         /* The bottom bit of the window */
396
397         if (!BN_to_montgomery(r,BN_value_one(),mont,ctx)) goto err;
398         for (;;)
399                 {
400                 if (BN_is_bit_set(p,wstart) == 0)
401                         {
402                         if (!start)
403                                 {
404                                 if (!BN_mod_mul_montgomery(r,r,r,mont,ctx))
405                                 goto err;
406                                 }
407                         if (wstart == 0) break;
408                         wstart--;
409                         continue;
410                         }
411                 /* We now have wstart on a 'set' bit, we now need to work out
412                  * how bit a window to do.  To do this we need to scan
413                  * forward until the last set bit before the end of the
414                  * window */
415                 j=wstart;
416                 wvalue=1;
417                 wend=0;
418                 for (i=1; i<window; i++)
419                         {
420                         if (wstart-i < 0) break;
421                         if (BN_is_bit_set(p,wstart-i))
422                                 {
423                                 wvalue<<=(i-wend);
424                                 wvalue|=1;
425                                 wend=i;
426                                 }
427                         }
428
429                 /* wend is the size of the current window */
430                 j=wend+1;
431                 /* add the 'bytes above' */
432                 if (!start)
433                         for (i=0; i<j; i++)
434                                 {
435                                 if (!BN_mod_mul_montgomery(r,r,r,mont,ctx))
436                                         goto err;
437                                 }
438                 
439                 /* wvalue will be an odd number < 2^window */
440                 if (!BN_mod_mul_montgomery(r,r,&(val[wvalue>>1]),mont,ctx))
441                         goto err;
442
443                 /* move the 'window' down further */
444                 wstart-=wend+1;
445                 wvalue=0;
446                 start=0;
447                 if (wstart < 0) break;
448                 }
449         if (!BN_from_montgomery(rr,r,mont,ctx)) goto err;
450         ret=1;
451 err:
452         if ((in_mont == NULL) && (mont != NULL)) BN_MONT_CTX_free(mont);
453         BN_CTX_end(ctx);
454         for (i=0; i<ts; i++)
455                 BN_clear_free(&(val[i]));
456         return(ret);
457         }
458
459 int BN_mod_exp_mont_word(BIGNUM *rr, BN_ULONG a, const BIGNUM *p,
460                          const BIGNUM *m, BN_CTX *ctx, BN_MONT_CTX *in_mont)
461         {
462         BN_MONT_CTX *mont = NULL;
463         int b, bits, ret=0;
464         int r_is_one;
465         BN_ULONG w, next_w;
466         BIGNUM *d, *r, *t;
467         BIGNUM *swap_tmp;
468 #define BN_MOD_MUL_WORD(r, w, m) \
469                 (BN_mul_word(r, (w)) && \
470                 (/* BN_ucmp(r, (m)) < 0 ? 1 :*/  \
471                         (BN_mod(t, r, m, ctx) && (swap_tmp = r, r = t, t = swap_tmp, 1))))
472                 /* BN_MOD_MUL_WORD is only used with 'w' large,
473                   * so the BN_ucmp test is probably more overhead
474                   * than always using BN_mod (which uses BN_copy if
475                   * a similar test returns true). */
476 #define BN_TO_MONTGOMERY_WORD(r, w, mont) \
477                 (BN_set_word(r, (w)) && BN_to_montgomery(r, r, (mont), ctx))
478
479         bn_check_top(p);
480         bn_check_top(m);
481
482         if (!(m->d[0] & 1))
483                 {
484                 BNerr(BN_F_BN_MOD_EXP_MONT_WORD,BN_R_CALLED_WITH_EVEN_MODULUS);
485                 return(0);
486                 }
487         bits = BN_num_bits(p);
488         if (bits == 0)
489                 {
490                 BN_one(rr);
491                 return(1);
492                 }
493         BN_CTX_start(ctx);
494         d = BN_CTX_get(ctx);
495         r = BN_CTX_get(ctx);
496         t = BN_CTX_get(ctx);
497         if (d == NULL || r == NULL || t == NULL) goto err;
498
499         if (in_mont != NULL)
500                 mont=in_mont;
501         else
502                 {
503                 if ((mont = BN_MONT_CTX_new()) == NULL) goto err;
504                 if (!BN_MONT_CTX_set(mont, m, ctx)) goto err;
505                 }
506
507         r_is_one = 1; /* except for Montgomery factor */
508
509         /* bits-1 >= 0 */
510
511         /* The result is accumulated in the product r*w. */
512         w = a; /* bit 'bits-1' of 'p' is always set */
513         for (b = bits-2; b >= 0; b--)
514                 {
515                 /* First, square r*w. */
516                 next_w = w*w;
517                 if ((next_w/w) != w) /* overflow */
518                         {
519                         if (r_is_one)
520                                 {
521                                 if (!BN_TO_MONTGOMERY_WORD(r, w, mont)) goto err;
522                                 r_is_one = 0;
523                                 }
524                         else
525                                 {
526                                 if (!BN_MOD_MUL_WORD(r, w, m)) goto err;
527                                 }
528                         next_w = 1;
529                         }
530                 w = next_w;
531                 if (!r_is_one)
532                         {
533                         if (!BN_mod_mul_montgomery(r, r, r, mont, ctx)) goto err;
534                         }
535
536                 /* Second, multiply r*w by 'a' if exponent bit is set. */
537                 if (BN_is_bit_set(p, b))
538                         {
539                         next_w = w*a;
540                         if ((next_w/a) != w) /* overflow */
541                                 {
542                                 if (r_is_one)
543                                         {
544                                         if (!BN_TO_MONTGOMERY_WORD(r, w, mont)) goto err;
545                                         r_is_one = 0;
546                                         }
547                                 else
548                                         {
549                                         if (!BN_MOD_MUL_WORD(r, w, m)) goto err;
550                                         }
551                                 next_w = a;
552                                 }
553                         w = next_w;
554                         }
555                 }
556
557         /* Finally, set r:=r*w. */
558         if (w != 1)
559                 {
560                 if (r_is_one)
561                         {
562                         if (!BN_TO_MONTGOMERY_WORD(r, w, mont)) goto err;
563                         r_is_one = 0;
564                         }
565                 else
566                         {
567                         if (!BN_MOD_MUL_WORD(r, w, m)) goto err;
568                         }
569                 }
570
571         if (r_is_one) /* can happen only if a == 1*/
572                 {
573                 if (!BN_one(rr)) goto err;
574                 }
575         else
576                 {
577                 if (!BN_from_montgomery(rr, r, mont, ctx)) goto err;
578                 }
579         ret = 1;
580 err:
581         if ((in_mont == NULL) && (mont != NULL)) BN_MONT_CTX_free(mont);
582         BN_CTX_end(ctx);
583         return(ret);
584         }
585
586
587 /* The old fallback, simple version :-) */
588 int BN_mod_exp_simple(BIGNUM *r,
589         const BIGNUM *a, const BIGNUM *p, const BIGNUM *m,
590         BN_CTX *ctx)
591         {
592         int i,j,bits,ret=0,wstart,wend,window,wvalue,ts=0;
593         int start=1;
594         BIGNUM *d;
595         BIGNUM val[TABLE_SIZE];
596
597         bits=BN_num_bits(p);
598
599         if (bits == 0)
600                 {
601                 BN_one(r);
602                 return(1);
603                 }
604
605         BN_CTX_start(ctx);
606         if ((d = BN_CTX_get(ctx)) == NULL) goto err;
607
608         BN_init(&(val[0]));
609         ts=1;
610         if (!BN_nnmod(&(val[0]),a,m,ctx)) goto err;             /* 1 */
611
612         window = BN_window_bits_for_exponent_size(bits);
613         if (window > 1)
614                 {
615                 if (!BN_mod_mul(d,&(val[0]),&(val[0]),m,ctx))
616                         goto err;                               /* 2 */
617                 j=1<<(window-1);
618                 for (i=1; i<j; i++)
619                         {
620                         BN_init(&(val[i]));
621                         if (!BN_mod_mul(&(val[i]),&(val[i-1]),d,m,ctx))
622                                 goto err;
623                         }
624                 ts=i;
625                 }
626
627         start=1;        /* This is used to avoid multiplication etc
628                          * when there is only the value '1' in the
629                          * buffer. */
630         wvalue=0;       /* The 'value' of the window */
631         wstart=bits-1;  /* The top bit of the window */
632         wend=0;         /* The bottom bit of the window */
633
634         if (!BN_one(r)) goto err;
635
636         for (;;)
637                 {
638                 if (BN_is_bit_set(p,wstart) == 0)
639                         {
640                         if (!start)
641                                 if (!BN_mod_mul(r,r,r,m,ctx))
642                                 goto err;
643                         if (wstart == 0) break;
644                         wstart--;
645                         continue;
646                         }
647                 /* We now have wstart on a 'set' bit, we now need to work out
648                  * how bit a window to do.  To do this we need to scan
649                  * forward until the last set bit before the end of the
650                  * window */
651                 j=wstart;
652                 wvalue=1;
653                 wend=0;
654                 for (i=1; i<window; i++)
655                         {
656                         if (wstart-i < 0) break;
657                         if (BN_is_bit_set(p,wstart-i))
658                                 {
659                                 wvalue<<=(i-wend);
660                                 wvalue|=1;
661                                 wend=i;
662                                 }
663                         }
664
665                 /* wend is the size of the current window */
666                 j=wend+1;
667                 /* add the 'bytes above' */
668                 if (!start)
669                         for (i=0; i<j; i++)
670                                 {
671                                 if (!BN_mod_mul(r,r,r,m,ctx))
672                                         goto err;
673                                 }
674                 
675                 /* wvalue will be an odd number < 2^window */
676                 if (!BN_mod_mul(r,r,&(val[wvalue>>1]),m,ctx))
677                         goto err;
678
679                 /* move the 'window' down further */
680                 wstart-=wend+1;
681                 wvalue=0;
682                 start=0;
683                 if (wstart < 0) break;
684                 }
685         ret=1;
686 err:
687         BN_CTX_end(ctx);
688         for (i=0; i<ts; i++)
689                 BN_clear_free(&(val[i]));
690         return(ret);
691         }
692