1
2 /***************************************************************************
3 *
4 Copyright 2012 CertiVox IOM Ltd. *
5 *
6 This file is part of CertiVox MIRACL Crypto SDK. *
7 *
8 The CertiVox MIRACL Crypto SDK provides developers with an *
9 extensive and efficient set of cryptographic functions. *
10 For further information about its features and functionalities please *
11 refer to http://www.certivox.com *
12 *
13 * The CertiVox MIRACL Crypto SDK is free software: you can *
14 redistribute it and/or modify it under the terms of the *
15 GNU Affero General Public License as published by the *
16 Free Software Foundation, either version 3 of the License, *
17 or (at your option) any later version. *
18 *
19 * The CertiVox MIRACL Crypto SDK is distributed in the hope *
20 that it will be useful, but WITHOUT ANY WARRANTY; without even the *
21 implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. *
22 See the GNU Affero General Public License for more details. *
23 *
24 * You should have received a copy of the GNU Affero General Public *
25 License along with CertiVox MIRACL Crypto SDK. *
26 If not, see <http://www.gnu.org/licenses/>. *
27 *
28 You can be released from the requirements of the license by purchasing *
29 a commercial license. Buying such a license is mandatory as soon as you *
30 develop commercial activities involving the CertiVox MIRACL Crypto SDK *
31 without disclosing the source code of your own applications, or shipping *
32 the CertiVox MIRACL Crypto SDK with a closed source product. *
33 *
34 ***************************************************************************/
35 /*
36 * MIRACL methods for modular exponentiation
37 * mrpower.c
38 */
39
40 #include <stdlib.h>
41 #include "miracl.h"
42
nres_powltr(_MIPD_ int x,big y,big w)43 void nres_powltr(_MIPD_ int x,big y,big w)
44 { /* calculates w=x^y mod z using Left to Right Method */
45 /* uses only n^2 modular squarings, because x is small */
46 /* Note: x is NOT an nresidue */
47 int i,nb;
48 #ifdef MR_OS_THREADS
49 miracl *mr_mip=get_mip();
50 #endif
51
52 if (mr_mip->ERNUM) return;
53 copy(y,mr_mip->w1);
54
55 MR_IN(86)
56
57 zero(w);
58 if (x==0)
59 {
60 if (size(mr_mip->w1)==0)
61 { /* 0^0 = 1 */
62 copy(mr_mip->one,w);
63 }
64 MR_OUT
65 return;
66 }
67
68 copy(mr_mip->one,w);
69 if (size(mr_mip->w1)==0)
70 {
71 MR_OUT
72 return;
73 }
74 if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
75 if (mr_mip->ERNUM)
76 {
77 MR_OUT
78 return;
79 }
80
81 #ifndef MR_ALWAYS_BINARY
82 if (mr_mip->base==mr_mip->base2)
83 {
84 #endif
85 nb=logb2(_MIPP_ mr_mip->w1);
86 convert(_MIPP_ x,w);
87 nres(_MIPP_ w,w);
88 if (nb>1) for (i=nb-2;i>=0;i--)
89 { /* Left to Right binary method */
90
91 if (mr_mip->user!=NULL) (*mr_mip->user)();
92
93 nres_modmult(_MIPP_ w,w,w);
94 if (mr_testbit(_MIPP_ mr_mip->w1,i))
95 { /* this is quick */
96 premult(_MIPP_ w,x,w);
97 divide(_MIPP_ w,mr_mip->modulus,mr_mip->modulus);
98 }
99 }
100 #ifndef MR_ALWAYS_BINARY
101 }
102 else
103 {
104 expb2(_MIPP_ logb2(_MIPP_ mr_mip->w1)-1,mr_mip->w2);
105 while (size(mr_mip->w2)!=0)
106 { /* Left to Right binary method */
107
108 if (mr_mip->user!=NULL) (*mr_mip->user)();
109 if (mr_mip->ERNUM) break;
110 nres_modmult(_MIPP_ w,w,w);
111 if (mr_compare(mr_mip->w1,mr_mip->w2)>=0)
112 {
113 premult(_MIPP_ w,x,w);
114 divide(_MIPP_ w,mr_mip->modulus,mr_mip->modulus);
115 subtract(_MIPP_ mr_mip->w1,mr_mip->w2,mr_mip->w1);
116 }
117 subdiv(_MIPP_ mr_mip->w2,2,mr_mip->w2);
118 }
119 }
120 #endif
121 if (size(w)<0) add(_MIPP_ w,mr_mip->modulus,w);
122 MR_OUT
123 return;
124 }
125
126 #ifndef MR_STATIC
127
nres_powmodn(_MIPD_ int n,big * x,big * y,big w)128 void nres_powmodn(_MIPD_ int n,big *x,big *y,big w)
129 {
130 int i,j,k,m,nb,ea;
131 big *G;
132 #ifdef MR_OS_THREADS
133 miracl *mr_mip=get_mip();
134 #endif
135 if (mr_mip->ERNUM) return;
136
137 MR_IN(112)
138
139 m=1<<n;
140 G=(big *)mr_alloc(_MIPP_ m,sizeof(big));
141
142 /* 2^n - n - 1 modmults */
143 /* 4(n=3) 11(n=4) etc */
144
145 for (i=0,k=1;i<n;i++)
146 {
147 for (j=0; j < (1<<i) ;j++)
148 {
149 G[k]=mirvar(_MIPP_ 0);
150 if (j==0) copy(x[i],G[k]);
151 else nres_modmult(_MIPP_ G[j],x[i],G[k]);
152 k++;
153 }
154 }
155
156 nb=0;
157 for (j=0;j<n;j++)
158 if ((k=logb2(_MIPP_ y[j]))>nb) nb=k;
159
160 copy(mr_mip->one,w);
161
162 #ifndef MR_ALWAYS_BINARY
163
164 if (mr_mip->base==mr_mip->base2)
165 {
166 #endif
167 for (i=nb-1;i>=0;i--)
168 {
169 if (mr_mip->user!=NULL) (*mr_mip->user)();
170 ea=0;
171 k=1;
172 for (j=0;j<n;j++)
173 {
174 if (mr_testbit(_MIPP_ y[j],i)) ea+=k;
175 k<<=1;
176 }
177 nres_modmult(_MIPP_ w,w,w);
178 if (ea!=0) nres_modmult(_MIPP_ w,G[ea],w);
179 }
180
181 #ifndef MR_ALWAYS_BINARY
182 }
183 else mr_berror(_MIPP_ MR_ERR_NOT_SUPPORTED);
184 #endif
185
186 for (i=1;i<m;i++) mirkill(G[i]);
187 mr_free(G);
188
189 MR_OUT
190 }
191
powmodn(_MIPD_ int n,big * x,big * y,big p,big w)192 void powmodn(_MIPD_ int n,big *x,big *y,big p,big w)
193 {/* w=x[0]^y[0].x[1]^y[1] .... x[n-1]^y[n-1] mod n */
194 int j;
195 #ifdef MR_OS_THREADS
196 miracl *mr_mip=get_mip();
197 #endif
198 if (mr_mip->ERNUM) return;
199
200 MR_IN(113)
201
202 prepare_monty(_MIPP_ p);
203 for (j=0;j<n;j++) nres(_MIPP_ x[j],x[j]);
204 nres_powmodn(_MIPP_ n,x,y,w);
205 for (j=0;j<n;j++) redc(_MIPP_ x[j],x[j]);
206
207 redc(_MIPP_ w,w);
208
209 MR_OUT
210 }
211
212 #endif
213
nres_powmod2(_MIPD_ big x,big y,big a,big b,big w)214 void nres_powmod2(_MIPD_ big x,big y,big a,big b,big w)
215 { /* finds w = x^y.a^b mod n. Fast for some cryptosystems */
216 int i,j,nb,nb2,nbw,nzs,n;
217 big table[16];
218 #ifdef MR_OS_THREADS
219 miracl *mr_mip=get_mip();
220 #endif
221 if (mr_mip->ERNUM) return;
222 copy(y,mr_mip->w1);
223 copy(x,mr_mip->w2);
224 copy(b,mr_mip->w3);
225 copy(a,mr_mip->w4);
226 zero(w);
227 if (size(mr_mip->w2)==0 || size(mr_mip->w4)==0) return;
228
229 MR_IN(99)
230
231 copy(mr_mip->one,w);
232 if (size(mr_mip->w1)==0 && size(mr_mip->w3)==0)
233 {
234 MR_OUT
235 return;
236 }
237 if (size(mr_mip->w1)<0 || size(mr_mip->w3)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
238 if (mr_mip->ERNUM)
239 {
240 MR_OUT
241 return;
242 }
243
244 #ifndef MR_ALWAYS_BINARY
245
246 if (mr_mip->base==mr_mip->base2)
247 { /* uses 2-bit sliding window. This is 25% faster! */
248 #endif
249 nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w4,mr_mip->w5);
250 nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w2,mr_mip->w12);
251 nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w4,mr_mip->w13);
252 nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w13,mr_mip->w14);
253 nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w13,mr_mip->w6);
254 nres_modmult(_MIPP_ mr_mip->w6,mr_mip->w4,mr_mip->w15);
255 nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w12,mr_mip->w8);
256 nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w12,mr_mip->w9);
257 nres_modmult(_MIPP_ mr_mip->w4,mr_mip->w9,mr_mip->w10);
258 nres_modmult(_MIPP_ mr_mip->w14,mr_mip->w12,mr_mip->w11);
259 nres_modmult(_MIPP_ mr_mip->w9,mr_mip->w13,mr_mip->w12);
260 nres_modmult(_MIPP_ mr_mip->w10,mr_mip->w13,mr_mip->w13);
261
262 table[0]=NULL; table[1]=mr_mip->w4; table[2]=mr_mip->w2; table[3]=mr_mip->w5;
263 table[4]=NULL; table[5]=mr_mip->w14; table[6]=mr_mip->w6; table[7]=mr_mip->w15;
264 table[8]=NULL; table[9]=mr_mip->w8; table[10]=mr_mip->w9; table[11]=mr_mip->w10;
265 table[12]=NULL;table[13]=mr_mip->w11;table[14]=mr_mip->w12; table[15]=mr_mip->w13;
266 nb=logb2(_MIPP_ mr_mip->w1);
267 if ((nb2=logb2(_MIPP_ mr_mip->w3))>nb) nb=nb2;
268
269 for (i=nb-1;i>=0;)
270 {
271 if (mr_mip->user!=NULL) (*mr_mip->user)();
272
273 n=mr_window2(_MIPP_ mr_mip->w1,mr_mip->w3,i,&nbw,&nzs);
274 for (j=0;j<nbw;j++)
275 nres_modmult(_MIPP_ w,w,w);
276 if (n>0) nres_modmult(_MIPP_ w,table[n],w);
277 i-=nbw;
278 if (nzs)
279 {
280 nres_modmult(_MIPP_ w,w,w);
281 i--;
282 }
283 }
284 #ifndef MR_ALWAYS_BINARY
285 }
286 else
287 {
288 nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w4,mr_mip->w5);
289
290 if (mr_compare(mr_mip->w1,mr_mip->w3)>=0)
291 expb2(_MIPP_ logb2(_MIPP_ mr_mip->w1)-1,mr_mip->w6);
292 else expb2(_MIPP_ logb2(_MIPP_ mr_mip->w3)-1,mr_mip->w6);
293 while (size(mr_mip->w6)!=0)
294 {
295 if (mr_mip->user!=NULL) (*mr_mip->user)();
296 if (mr_mip->ERNUM) break;
297 nres_modmult(_MIPP_ w,w,w);
298 if (mr_compare(mr_mip->w1,mr_mip->w6)>=0)
299 {
300 if (mr_compare(mr_mip->w3,mr_mip->w6)>=0)
301 {
302 nres_modmult(_MIPP_ w,mr_mip->w5,w);
303 subtract(_MIPP_ mr_mip->w3,mr_mip->w6,mr_mip->w3);
304 }
305
306 else nres_modmult(_MIPP_ w,mr_mip->w2,w);
307 subtract(_MIPP_ mr_mip->w1,mr_mip->w6,mr_mip->w1);
308 }
309 else
310 {
311 if (mr_compare(mr_mip->w3,mr_mip->w6)>=0)
312 {
313 nres_modmult(_MIPP_ w,mr_mip->w4,w);
314 subtract(_MIPP_ mr_mip->w3,mr_mip->w6,mr_mip->w3);
315 }
316 }
317 subdiv(_MIPP_ mr_mip->w6,2,mr_mip->w6);
318 }
319 }
320 #endif
321 MR_OUT
322 }
323
powmod2(_MIPD_ big x,big y,big a,big b,big n,big w)324 void powmod2(_MIPD_ big x,big y,big a,big b,big n,big w)
325 {/* w=x^y.a^b mod n */
326 #ifdef MR_OS_THREADS
327 miracl *mr_mip=get_mip();
328 #endif
329 if (mr_mip->ERNUM) return;
330
331 MR_IN(79)
332
333 prepare_monty(_MIPP_ n);
334 nres(_MIPP_ x,mr_mip->w2);
335 nres(_MIPP_ a,mr_mip->w4);
336 nres_powmod2(_MIPP_ mr_mip->w2,y,mr_mip->w4,b,w);
337 redc(_MIPP_ w,w);
338
339 MR_OUT
340 }
341
342
powmod(_MIPD_ big x,big y,big n,big w)343 void powmod(_MIPD_ big x,big y,big n,big w)
344 { /* fast powmod, using Montgomery's method internally */
345
346 mr_small norm;
347 BOOL mty;
348 #ifdef MR_OS_THREADS
349 miracl *mr_mip=get_mip();
350 #endif
351 if (mr_mip->ERNUM) return;
352
353 MR_IN(18)
354
355 mty=TRUE;
356
357 if (mr_mip->base!=mr_mip->base2)
358 {
359 if (size(n)<2 || sgcd(n->w[0],mr_mip->base)!=1) mty=FALSE;
360 }
361 else
362 if (subdivisible(_MIPP_ n,2)) mty=FALSE;
363
364 if (!mty)
365 { /* can't use Montgomery */
366 copy(y,mr_mip->w1);
367 copy(x,mr_mip->w3);
368 zero(w);
369 if (size(mr_mip->w3)==0)
370 {
371 MR_OUT
372 return;
373 }
374 convert(_MIPP_ 1,w);
375 if (size(mr_mip->w1)==0)
376 {
377 MR_OUT
378 return;
379 }
380 if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
381 if (w==n) mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS) ;
382 if (mr_mip->ERNUM)
383 {
384 MR_OUT
385 return;
386 }
387
388 norm=normalise(_MIPP_ n,n);
389 divide(_MIPP_ mr_mip->w3,n,n);
390 forever
391 {
392 if (mr_mip->user!=NULL) (*mr_mip->user)();
393
394 if (subdiv(_MIPP_ mr_mip->w1,2,mr_mip->w1)!=0)
395 mad(_MIPP_ w,mr_mip->w3,mr_mip->w3,n,n,w);
396 if (mr_mip->ERNUM || size(mr_mip->w1)==0) break;
397 mad(_MIPP_ mr_mip->w3,mr_mip->w3,mr_mip->w3,n,n,mr_mip->w3);
398 }
399 if (norm!=1)
400 {
401 #ifdef MR_FP_ROUNDING
402 mr_sdiv(_MIPP_ n,norm,mr_invert(norm),n);
403 #else
404 mr_sdiv(_MIPP_ n,norm,n);
405 #endif
406 divide(_MIPP_ w,n,n);
407 }
408 }
409 else
410 { /* optimized code for odd moduli */
411 prepare_monty(_MIPP_ n);
412 nres(_MIPP_ x,mr_mip->w3);
413 nres_powmod(_MIPP_ mr_mip->w3,y,w);
414 redc(_MIPP_ w,w);
415 }
416
417 MR_OUT
418 }
419
powltr(_MIPD_ int x,big y,big n,big w)420 int powltr(_MIPD_ int x, big y, big n, big w)
421 {
422 mr_small norm;
423 BOOL clean_up,mty;
424 #ifdef MR_OS_THREADS
425 miracl *mr_mip=get_mip();
426 #endif
427 if (mr_mip->ERNUM) return 0;
428
429 MR_IN(71)
430 mty=TRUE;
431
432 if (mr_mip->base!=mr_mip->base2)
433 {
434 if (size(n)<2 || sgcd(n->w[0],mr_mip->base)!=1) mty=FALSE;
435 }
436 else
437 {
438 if (subdivisible(_MIPP_ n,2)) mty=FALSE;
439 }
440
441 /* Is Montgomery modulus in use... */
442
443 clean_up=FALSE;
444 if (mty)
445 { /* still a chance to use monty... */
446 if (size(mr_mip->modulus)!=0)
447 { /* somebody else is using it */
448 if (mr_compare(n,mr_mip->modulus)!=0) mty=FALSE;
449 }
450 else
451 { /* its unused, so use it, but clean up after */
452 clean_up=TRUE;
453 }
454 }
455 if (!mty)
456 { /* can't use Montgomery! */
457 copy(y,mr_mip->w1);
458 zero(w);
459 if (x==0)
460 {
461 MR_OUT
462 return 0;
463 }
464 convert(_MIPP_ 1,w);
465 if (size(mr_mip->w1)==0)
466 {
467 MR_OUT
468 return 1;
469 }
470
471 if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
472 if (w==n) mr_berror(_MIPP_ MR_ERR_BAD_PARAMETERS) ;
473 if (mr_mip->ERNUM)
474 {
475 MR_OUT
476 return 0;
477 }
478
479 norm=normalise(_MIPP_ n,n);
480
481 expb2(_MIPP_ logb2(_MIPP_ mr_mip->w1)-1,mr_mip->w2);
482 while (size(mr_mip->w2)!=0)
483 { /* Left to Right binary method */
484
485 if (mr_mip->user!=NULL) (*mr_mip->user)();
486 if (mr_mip->ERNUM) break;
487 mad(_MIPP_ w,w,w,n,n,w);
488 if (mr_compare(mr_mip->w1,mr_mip->w2)>=0)
489 {
490 premult(_MIPP_ w,x,w);
491 divide(_MIPP_ w,n,n);
492 subtract(_MIPP_ mr_mip->w1,mr_mip->w2,mr_mip->w1);
493 }
494 subdiv(_MIPP_ mr_mip->w2,2,mr_mip->w2);
495 }
496 if (norm!=1)
497 {
498 #ifdef MR_FP_ROUNDING
499 mr_sdiv(_MIPP_ n,norm,mr_invert(norm),n);
500 #else
501 mr_sdiv(_MIPP_ n,norm,n);
502 #endif
503 divide(_MIPP_ w,n,n);
504 }
505 }
506 else
507 { /* optimized code for odd moduli */
508 prepare_monty(_MIPP_ n);
509 nres_powltr(_MIPP_ x,y,w);
510 redc(_MIPP_ w,w);
511 if (clean_up) kill_monty(_MIPPO_ );
512 }
513 MR_OUT
514 return (size(w));
515 }
516
nres_powmod(_MIPD_ big x,big y,big w)517 void nres_powmod(_MIPD_ big x,big y,big w)
518 { /* calculates w=x^y mod z, using m-residues *
519 * See "Analysis of Sliding Window Techniques for *
520 * Exponentiation, C.K. Koc, Computers. Math. & *
521 * Applic. Vol. 30 pp17-24 1995. Uses work-space *
522 * variables for pre-computed table. */
523 int i,j,k,t,nb,nbw,nzs,n;
524 big table[16];
525 #ifdef MR_OS_THREADS
526 miracl *mr_mip=get_mip();
527 #endif
528 if (mr_mip->ERNUM) return;
529
530 copy(y,mr_mip->w1);
531 copy(x,mr_mip->w3);
532
533 MR_IN(84)
534 zero(w);
535 if (size(x)==0)
536 {
537 if (size(mr_mip->w1)==0)
538 { /* 0^0 = 1 */
539 copy(mr_mip->one,w);
540 }
541 MR_OUT
542 return;
543 }
544
545 copy(mr_mip->one,w);
546 if (size(mr_mip->w1)==0)
547 {
548 MR_OUT
549 return;
550 }
551
552 if (size(mr_mip->w1)<0) mr_berror(_MIPP_ MR_ERR_NEG_POWER);
553
554 if (mr_mip->ERNUM)
555 {
556 MR_OUT
557 return;
558 }
559
560 #ifndef MR_ALWAYS_BINARY
561 if (mr_mip->base==mr_mip->base2)
562 { /* build a table. Up to 5-bit sliding windows. Windows with
563 * two adjacent 0 bits simply won't happen */
564 #endif
565 table[0]=mr_mip->w3; table[1]=mr_mip->w4; table[2]=mr_mip->w5; table[3]=mr_mip->w14;
566 table[4]=NULL; table[5]=mr_mip->w6; table[6]=mr_mip->w15; table[7]=mr_mip->w8;
567 table[8]=NULL; table[9]=NULL; table[10]=mr_mip->w9; table[11]=mr_mip->w10;
568 table[12]=NULL; table[13]=mr_mip->w11; table[14]=mr_mip->w12; table[15]=mr_mip->w13;
569
570 nres_modmult(_MIPP_ mr_mip->w3,mr_mip->w3,mr_mip->w2); /* x^2 */
571 n=15;
572 j=0;
573 do
574 { /* pre-computations */
575 t=1; k=j+1;
576 while (table[k]==NULL) {k++; t++;}
577 copy(table[j],table[k]);
578 for (i=0;i<t;i++) nres_modmult(_MIPP_ table[k],mr_mip->w2,table[k]);
579 j=k;
580 } while (j<n);
581
582 nb=logb2(_MIPP_ mr_mip->w1);
583 copy(mr_mip->w3,w);
584 if (nb>1) for (i=nb-2;i>=0;)
585 { /* Left to Right method */
586
587 if (mr_mip->user!=NULL) (*mr_mip->user)();
588
589 n=mr_window(_MIPP_ mr_mip->w1,i,&nbw,&nzs,5);
590
591 for (j=0;j<nbw;j++)
592 nres_modmult(_MIPP_ w,w,w);
593 if (n>0)
594 nres_modmult(_MIPP_ w,table[n/2],w);
595
596 i-=nbw;
597 if (nzs)
598 {
599 for (j=0;j<nzs;j++)
600 nres_modmult(_MIPP_ w,w,w);
601 i-=nzs;
602 }
603 }
604
605 #ifndef MR_ALWAYS_BINARY
606 }
607 else
608 {
609 copy(mr_mip->w3,mr_mip->w2);
610 forever
611 { /* "Russian peasant" Right-to-Left exponentiation */
612
613 if (mr_mip->user!=NULL) (*mr_mip->user)();
614
615 if (subdiv(_MIPP_ mr_mip->w1,2,mr_mip->w1)!=0)
616 nres_modmult(_MIPP_ w,mr_mip->w2,w);
617 if (mr_mip->ERNUM || size(mr_mip->w1)==0) break;
618 nres_modmult(_MIPP_ mr_mip->w2,mr_mip->w2,mr_mip->w2);
619 }
620 }
621 #endif
622 MR_OUT
623 }
624
625