1 /*
2 * Copyright (c) 2002, 2003 Bob Deblier
3 *
4 * This library is free software; you can redistribute it and/or
5 * modify it under the terms of the GNU Lesser General Public
6 * License as published by the Free Software Foundation; either
7 * version 2.1 of the License, or (at your option) any later version.
8 *
9 * This library is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 * Lesser General Public License for more details.
13 *
14 * You should have received a copy of the GNU Lesser General Public
15 * License along with this library; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 *
18 */
19
20 /*!\file mp.c
21 * \brief Multi-precision integer routines.
22 * \author Bob Deblier <bob.deblier@telenet.be>
23 * \ingroup MP_m
24 */
25
26 #define BEECRYPT_DLL_EXPORT
27
28 #if HAVE_CONFIG_H
29 # include "config.h"
30 #endif
31
32 #include "beecrypt/mp.h"
33 #include "beecrypt/mpopt.h"
34
35 #ifndef ASM_MPZERO
mpzero(size_t size,mpw * data)36 void mpzero(size_t size, mpw* data)
37 {
38 while (size--)
39 *(data++) = 0;
40 }
41 #endif
42
43 #ifndef ASM_MPFILL
mpfill(size_t size,mpw * data,mpw fill)44 void mpfill(size_t size, mpw* data, mpw fill)
45 {
46 while (size--)
47 *(data++) = fill;
48 }
49 #endif
50
51 #ifndef ASM_MPODD
mpodd(size_t size,const mpw * data)52 int mpodd(size_t size, const mpw* data)
53 {
54 return (int)(data[size-1] & 0x1);
55 }
56 #endif
57
58 #ifndef ASM_MPEVEN
mpeven(size_t size,const mpw * data)59 int mpeven(size_t size, const mpw* data)
60 {
61 return !(int)(data[size-1] & 0x1);
62 }
63 #endif
64
65 #ifndef ASM_MPZ
mpz(size_t size,const mpw * data)66 int mpz(size_t size, const mpw* data)
67 {
68 while (size--)
69 if (*(data++))
70 return 0;
71 return 1;
72 }
73 #endif
74
75 #ifndef ASM_MPNZ
mpnz(size_t size,const mpw * data)76 int mpnz(size_t size, const mpw* data)
77 {
78 while (size--)
79 if (*(data++))
80 return 1;
81 return 0;
82 }
83 #endif
84
85 #ifndef ASM_MPEQ
mpeq(size_t size,const mpw * xdata,const mpw * ydata)86 int mpeq(size_t size, const mpw* xdata, const mpw* ydata)
87 {
88 while (size--)
89 {
90 if (*xdata == *ydata)
91 {
92 xdata++;
93 ydata++;
94 }
95 else
96 return 0;
97 }
98 return 1;
99 }
100 #endif
101
102 #ifndef ASM_MPEQX
mpeqx(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)103 int mpeqx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
104 {
105 if (xsize > ysize)
106 {
107 register size_t diff = xsize - ysize;
108 return mpeq(ysize, xdata+diff, ydata) && mpz(diff, xdata);
109 }
110 else if (xsize < ysize)
111 {
112 register size_t diff = ysize - xsize;
113 return mpeq(xsize, ydata+diff, xdata) && mpz(diff, ydata);
114 }
115 else
116 return mpeq(xsize, xdata, ydata);
117 }
118 #endif
119
120 #ifndef ASM_MPNE
mpne(size_t size,const mpw * xdata,const mpw * ydata)121 int mpne(size_t size, const mpw* xdata, const mpw* ydata)
122 {
123 while (size--)
124 {
125 if (*xdata == *ydata)
126 {
127 xdata++;
128 ydata++;
129 }
130 else
131 return 1;
132 }
133 return 0;
134 }
135 #endif
136
137 #ifndef ASM_MPNEX
mpnex(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)138 int mpnex(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
139 {
140 if (xsize > ysize)
141 {
142 register size_t diff = xsize - ysize;
143 return mpnz(diff, xdata) || mpne(ysize, xdata+diff, ydata);
144 }
145 else if (xsize < ysize)
146 {
147 register size_t diff = ysize - xsize;
148 return mpnz(diff, ydata) || mpne(xsize, ydata+diff, xdata);
149 }
150 else
151 return mpne(xsize, xdata, ydata);
152 }
153 #endif
154
155 #ifndef ASM_MPGT
mpgt(size_t size,const mpw * xdata,const mpw * ydata)156 int mpgt(size_t size, const mpw* xdata, const mpw* ydata)
157 {
158 while (size--)
159 {
160 if (*xdata < *ydata)
161 return 0;
162 if (*xdata > *ydata)
163 return 1;
164 xdata++; ydata++;
165 }
166 return 0;
167 }
168 #endif
169
170 #ifndef ASM_MPGTX
mpgtx(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)171 int mpgtx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
172 {
173 if (xsize > ysize)
174 {
175 register size_t diff = xsize - ysize;
176 return mpnz(diff, xdata) || mpgt(ysize, xdata + diff, ydata);
177 }
178 else if (xsize < ysize)
179 {
180 register size_t diff = ysize - xsize;
181 return mpz(diff, ydata) && mpgt(xsize, xdata, ydata + diff);
182 }
183 else
184 return mpgt(xsize, xdata, ydata);
185 }
186 #endif
187
188 #ifndef ASM_MPLT
mplt(size_t size,const mpw * xdata,const mpw * ydata)189 int mplt(size_t size, const mpw* xdata, const mpw* ydata)
190 {
191 while (size--)
192 {
193 if (*xdata > *ydata)
194 return 0;
195 if (*xdata < *ydata)
196 return 1;
197 xdata++; ydata++;
198 }
199 return 0;
200 }
201 #endif
202
203 #ifndef ASM_MPLTX
mpltx(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)204 int mpltx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
205 {
206 if (xsize > ysize)
207 {
208 register size_t diff = xsize - ysize;
209 return mpz(diff, xdata) && mplt(ysize, xdata+diff, ydata);
210 }
211 else if (xsize < ysize)
212 {
213 register size_t diff = ysize - xsize;
214 return mpnz(diff, ydata) || mplt(xsize, xdata, ydata+diff);
215 }
216 else
217 return mplt(xsize, xdata, ydata);
218 }
219 #endif
220
221 #ifndef ASM_MPGE
mpge(size_t size,const mpw * xdata,const mpw * ydata)222 int mpge(size_t size, const mpw* xdata, const mpw* ydata)
223 {
224 while (size--)
225 {
226 if (*xdata < *ydata)
227 return 0;
228 if (*xdata > *ydata)
229 return 1;
230 xdata++; ydata++;
231 }
232 return 1;
233 }
234 #endif
235
236 #ifndef ASM_MPGEX
mpgex(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)237 int mpgex(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
238 {
239 if (xsize > ysize)
240 {
241 register size_t diff = xsize - ysize;
242 return mpnz(diff, xdata) || mpge(ysize, xdata+diff, ydata);
243 }
244 else if (xsize < ysize)
245 {
246 register size_t diff = ysize - xsize;
247 return mpz(diff, ydata) && mpge(xsize, xdata, ydata+diff);
248 }
249 else
250 return mpge(xsize, xdata, ydata);
251 }
252 #endif
253
254 #ifndef ASM_MPLE
mple(size_t size,const mpw * xdata,const mpw * ydata)255 int mple(size_t size, const mpw* xdata, const mpw* ydata)
256 {
257 while (size--)
258 {
259 if (*xdata < *ydata)
260 return 1;
261 if (*xdata > *ydata)
262 return 0;
263 xdata++; ydata++;
264 }
265 return 1;
266 }
267 #endif
268
269 #ifndef ASM_MPLEX
mplex(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)270 int mplex(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
271 {
272 if (xsize > ysize)
273 {
274 register size_t diff = xsize - ysize;
275 return mpz(diff, xdata) && mple(ysize, xdata+ diff, ydata);
276 }
277 else if (xsize < ysize)
278 {
279 register size_t diff = ysize - xsize;
280 return mpnz(diff, ydata) || mple(xsize, xdata, ydata+diff);
281 }
282 else
283 return mple(xsize, xdata, ydata);
284 }
285 #endif
286
287 #ifndef ASM_MPCMP
mpcmp(size_t size,const mpw * xdata,const mpw * ydata)288 int mpcmp(size_t size, const mpw* xdata, const mpw* ydata)
289 {
290 while (size--)
291 {
292 if (*xdata < *ydata)
293 return -1;
294 else if (*xdata > *ydata)
295 return 1;
296
297 xdata++;
298 ydata++;
299 }
300 return 0;
301 }
302 #endif
303
304 #ifndef ASM_MPCMPX
mpcmpx(size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)305 int mpcmpx(size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
306 {
307 if (xsize > ysize)
308 {
309 register size_t diff = xsize - ysize;
310 if (mpnz(diff, xdata))
311 return 1;
312
313 xsize -= diff;
314 xdata += diff;
315 }
316 else if (xsize < ysize)
317 {
318 register size_t diff = ysize - xsize;
319 if (mpnz(diff, ydata))
320 return -1;
321
322 ydata += diff;
323 }
324 return mpcmp(xsize, xdata, ydata);
325 }
326 #endif
327
328 #ifndef ASM_MPISONE
mpisone(size_t size,const mpw * data)329 int mpisone(size_t size, const mpw* data)
330 {
331 data += size;
332 if (*(--data) == 1)
333 {
334 while (--size)
335 if (*(--data))
336 return 0;
337 return 1;
338 }
339 return 0;
340 }
341 #endif
342
343 #ifndef ASM_MPISTWO
mpistwo(size_t size,const mpw * data)344 int mpistwo(size_t size, const mpw* data)
345 {
346 data += size;
347 if (*(--data) == 2)
348 {
349 while (--size)
350 if (*(--data))
351 return 0;
352 return 1;
353 }
354 return 0;
355 }
356 #endif
357
358 #ifndef ASM_MPEQMONE
mpeqmone(size_t size,const mpw * xdata,const mpw * ydata)359 int mpeqmone(size_t size, const mpw* xdata, const mpw* ydata)
360 {
361 xdata += size;
362 ydata += size;
363
364 if (*(--xdata)+1 == *(--ydata))
365 {
366 while (--size)
367 if (*(--xdata) != *(--ydata))
368 return 0;
369 return 1;
370 }
371 return 0;
372 }
373 #endif
374
375 #ifndef ASM_MPLEONE
mpleone(size_t size,const mpw * data)376 int mpleone(size_t size, const mpw* data)
377 {
378 data += size;
379 if (*(--data) > 1)
380 return 0;
381 else
382 {
383 while (--size)
384 if (*(--data))
385 return 0;
386 return 1;
387 }
388 }
389 #endif
390
391 #ifndef ASM_MPMSBSET
mpmsbset(size_t size,const mpw * data)392 int mpmsbset(size_t size, const mpw* data)
393 {
394 return (int)((*data) >> (MP_WBITS-1));
395 }
396 #endif
397
398 #ifndef ASM_MPLSBSET
mplsbset(size_t size,const mpw * data)399 int mplsbset(size_t size, const mpw* data)
400 {
401 return (int)(data[size-1] & 0x1);
402 }
403 #endif
404
405 #ifndef ASM_MPSETMSB
mpsetmsb(size_t size,mpw * data)406 void mpsetmsb(size_t size, mpw* data)
407 {
408 *data |= MP_MSBMASK;
409 }
410 #endif
411
412 #ifndef ASM_MPSETLSB
mpsetlsb(size_t size,mpw * data)413 void mpsetlsb(size_t size, mpw* data)
414 {
415 data[size-1] |= MP_LSBMASK;
416 }
417 #endif
418
419 #ifndef ASM_MPCLRMSB
mpclrmsb(size_t size,mpw * data)420 void mpclrmsb(size_t size, mpw* data)
421 {
422 *data &= ~ MP_MSBMASK;
423 }
424 #endif
425
426 #ifndef ASM_MPCLRLSB
mpclrlsb(size_t size,mpw * data)427 void mpclrlsb(size_t size, mpw* data)
428 {
429 data[size-1] &= ~ MP_LSBMASK;
430 }
431 #endif
432
433 #ifndef ASM_MPAND
mpand(size_t size,mpw * xdata,const mpw * ydata)434 void mpand(size_t size, mpw* xdata, const mpw* ydata)
435 {
436 while (size--)
437 xdata[size] &= ydata[size];
438 }
439 #endif
440
441 #ifndef ASM_MPOR
mpor(size_t size,mpw * xdata,const mpw * ydata)442 void mpor(size_t size, mpw* xdata, const mpw* ydata)
443 {
444 while (size--)
445 xdata[size] |= ydata[size];
446 }
447 #endif
448
449 #ifndef ASM_MPXOR
mpxor(size_t size,mpw * xdata,const mpw * ydata)450 void mpxor(size_t size, mpw* xdata, const mpw* ydata)
451 {
452 while (size--)
453 xdata[size] ^= ydata[size];
454 }
455 #endif
456
457 #ifndef ASM_MPNOT
mpnot(size_t size,mpw * data)458 void mpnot(size_t size, mpw* data)
459 {
460 while (size--)
461 data[size] = ~data[size];
462 }
463 #endif
464
465 #ifndef ASM_MPSETW
mpsetw(size_t size,mpw * xdata,mpw y)466 void mpsetw(size_t size, mpw* xdata, mpw y)
467 {
468 while (--size)
469 *(xdata++) = 0;
470 *(xdata++) = y;
471 }
472 #endif
473
474 #ifndef ASM_MPSETWS
mpsetws(size_t size,mpw * xdata,size_t y)475 void mpsetws(size_t size, mpw* xdata, size_t y)
476 {
477 while (size--)
478 {
479 xdata[size] = (mpw) y;
480 #if (MP_WBYTES > SIZEOF_SIZE_T)
481 y >>= MP_WBITS;
482 #else
483 y = 0;
484 #endif
485 }
486 }
487 #endif
488
489 #ifndef ASM_MPSETX
mpsetx(size_t xsize,mpw * xdata,size_t ysize,const mpw * ydata)490 void mpsetx(size_t xsize, mpw* xdata, size_t ysize, const mpw* ydata)
491 {
492 while (xsize > ysize)
493 {
494 xsize--;
495 *(xdata++) = 0;
496 }
497 while (ysize > xsize)
498 {
499 ysize--;
500 ydata++;
501 }
502 while (xsize--)
503 *(xdata++) = *(ydata++);
504 }
505 #endif
506
507 #ifndef ASM_MPADDW
mpaddw(size_t size,mpw * xdata,mpw y)508 int mpaddw(size_t size, mpw* xdata, mpw y)
509 {
510 register mpw load, temp;
511 register int carry = 0;
512
513 xdata += size-1;
514
515 load = *xdata;
516 temp = load + y;
517 *(xdata--) = temp;
518 carry = (load > temp);
519
520 while (--size && carry)
521 {
522 load = *xdata;
523 temp = load + 1;
524 *(xdata--) = temp;
525 carry = (load > temp);
526 }
527 return carry;
528 }
529 #endif
530
531 #ifndef ASM_MPADD
mpadd(size_t size,mpw * xdata,const mpw * ydata)532 int mpadd(size_t size, mpw* xdata, const mpw* ydata)
533 {
534 register mpw load, temp;
535 register int carry = 0;
536
537 xdata += size-1;
538 ydata += size-1;
539
540 while (size--)
541 {
542 temp = *(ydata--);
543 load = *xdata;
544 temp = carry ? (load + temp + 1) : (load + temp);
545 *(xdata--) = temp;
546 carry = carry ? (load >= temp) : (load > temp);
547 }
548 return carry;
549 }
550 #endif
551
552 #ifndef ASM_MPADDX
mpaddx(size_t xsize,mpw * xdata,size_t ysize,const mpw * ydata)553 int mpaddx(size_t xsize, mpw* xdata, size_t ysize, const mpw* ydata)
554 {
555 if (xsize > ysize)
556 {
557 register size_t diff = xsize - ysize;
558 return mpaddw(diff, xdata, (mpw) mpadd(ysize, xdata+diff, ydata));
559 }
560 else
561 {
562 register size_t diff = ysize - xsize;
563 return mpadd(xsize, xdata, ydata+diff);
564 }
565 }
566 #endif
567
568 #ifndef ASM_MPSUBW
mpsubw(size_t size,mpw * xdata,mpw y)569 int mpsubw(size_t size, mpw* xdata, mpw y)
570 {
571 register mpw load, temp;
572 register int carry = 0;
573
574 xdata += size-1;
575
576 load = *xdata;
577 temp = load - y;
578 *(xdata--) = temp;
579 carry = (load < temp);
580
581 while (--size && carry)
582 {
583 load = *xdata;
584 temp = load - 1;
585 *(xdata--) = temp;
586 carry = (load < temp);
587 }
588 return carry;
589 }
590 #endif
591
592 #ifndef ASM_MPSUB
mpsub(size_t size,mpw * xdata,const mpw * ydata)593 int mpsub(size_t size, mpw* xdata, const mpw* ydata)
594 {
595 register mpw load, temp;
596 register int carry = 0;
597
598 xdata += size-1;
599 ydata += size-1;
600
601 while (size--)
602 {
603 temp = *(ydata--);
604 load = *xdata;
605 temp = carry ? (load - temp - 1) : (load - temp);
606 *(xdata--) = temp;
607 carry = carry ? (load <= temp) : (load < temp);
608 }
609 return carry;
610 }
611 #endif
612
613 #ifndef ASM_MPSUBX
mpsubx(size_t xsize,mpw * xdata,size_t ysize,const mpw * ydata)614 int mpsubx(size_t xsize, mpw* xdata, size_t ysize, const mpw* ydata)
615 {
616 if (xsize > ysize)
617 {
618 register size_t diff = xsize - ysize;
619 return mpsubw(diff, xdata, (mpw) mpsub(ysize, xdata+diff, ydata));
620 }
621 else
622 {
623 register size_t diff = ysize - xsize;
624 return mpsub(xsize, xdata, ydata+diff);
625 }
626 }
627 #endif
628
629 #ifndef ASM_MPNEG
mpneg(size_t size,mpw * data)630 void mpneg(size_t size, mpw* data)
631 {
632 mpnot(size, data);
633 mpaddw(size, data, 1);
634 }
635 #endif
636
637 #ifndef ASM_MPSETMUL
mpsetmul(size_t size,mpw * result,const mpw * data,mpw y)638 mpw mpsetmul(size_t size, mpw* result, const mpw* data, mpw y)
639 {
640 #if HAVE_MPDW
641 register mpdw temp;
642 register mpw carry = 0;
643
644 data += size;
645 result += size;
646
647 while (size--)
648 {
649 temp = *(--data);
650 temp *= y;
651 temp += carry;
652 *(--result) = (mpw) temp;
653 carry = (mpw)(temp >> MP_WBITS);
654 }
655 #else
656 register mpw temp, load, carry = 0;
657 register mphw ylo, yhi;
658
659 ylo = (mphw) y;
660 yhi = (mphw) (y >> MP_HWBITS);
661
662 data += size;
663 result += size;
664
665 while (size--)
666 {
667 register mphw xlo, xhi;
668 register mpw rlo, rhi;
669
670 xlo = (mphw) (temp = *(--data));
671 xhi = (mphw) (temp >> MP_HWBITS);
672
673 rlo = (mpw) xlo * ylo;
674 rhi = (mpw) xhi * yhi;
675 load = rlo;
676 temp = (mpw) xhi * ylo;
677 rlo += (temp << MP_HWBITS);
678 rhi += (temp >> MP_HWBITS) + (load > rlo);
679 load = rlo;
680 temp = (mpw) xlo * yhi;
681 rlo += (temp << MP_HWBITS);
682 rhi += (temp >> MP_HWBITS) + (load > rlo);
683 load = rlo;
684 temp = rlo + carry;
685 carry = rhi + (load > temp);
686 *(--result) = temp;
687 }
688 #endif
689 return carry;
690 }
691 #endif
692
693 #ifndef ASM_MPADDMUL
mpaddmul(size_t size,mpw * result,const mpw * data,mpw y)694 mpw mpaddmul(size_t size, mpw* result, const mpw* data, mpw y)
695 {
696 #if HAVE_MPDW
697 register mpdw temp;
698 register mpw carry = 0;
699
700 data += size;
701 result += size;
702
703 while (size--)
704 {
705 temp = *(--data);
706 temp *= y;
707 temp += carry;
708 temp += *(--result);
709 *result = (mpw) temp;
710 carry = (mpw)(temp >> MP_WBITS);
711 }
712 #else
713 register mpw temp, load, carry = 0;
714 register mphw ylo, yhi;
715
716 ylo = (mphw) y;
717 yhi = (mphw) (y >> MP_HWBITS);
718
719 data += size;
720 result += size;
721
722 while (size--)
723 {
724 register mphw xlo, xhi;
725 register mpw rlo, rhi;
726
727 xlo = (mphw) (temp = *(--data));
728 xhi = (mphw) (temp >> MP_HWBITS);
729
730 rlo = (mpw) xlo * ylo;
731 rhi = (mpw) xhi * yhi;
732 load = rlo;
733 temp = (mpw) xhi * ylo;
734 rlo += (temp << MP_HWBITS);
735 rhi += (temp >> MP_HWBITS) + (load > rlo);
736 load = rlo;
737 temp = (mpw) xlo * yhi;
738 rlo += (temp << MP_HWBITS);
739 rhi += (temp >> MP_HWBITS) + (load > rlo);
740 load = rlo;
741 rlo += carry;
742 temp = (load > rlo);
743 load = rhi;
744 rhi += temp;
745 carry = (load > rhi);
746 load = rlo;
747 rlo += *(--result);
748 *result = rlo;
749 carry += rhi + (load > rlo);
750 }
751 #endif
752 return carry;
753 }
754 #endif
755
756 #ifndef ASM_MPMUL
mpmul(mpw * result,size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata)757 void mpmul(mpw* result, size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata)
758 {
759 /* preferred passing of parameters is x the larger of the two numbers */
760 if (xsize >= ysize)
761 {
762 register mpw rc;
763
764 result += ysize;
765 ydata += ysize;
766
767 rc = mpsetmul(xsize, result, xdata, *(--ydata));
768 *(--result) = rc;
769
770 while (--ysize)
771 {
772 rc = mpaddmul(xsize, result, xdata, *(--ydata));
773 *(--result) = rc;
774 }
775 }
776 else
777 {
778 register mpw rc;
779
780 result += xsize;
781 xdata += xsize;
782
783 rc = mpsetmul(ysize, result, ydata, *(--xdata));
784 *(--result) = rc;
785
786 while (--xsize)
787 {
788 rc = mpaddmul(ysize, result, ydata, *(--xdata));
789 *(--result) = rc;
790 }
791 }
792 }
793 #endif
794
795 #ifndef ASM_MPADDSQRTRC
mpaddsqrtrc(size_t size,mpw * result,const mpw * data)796 void mpaddsqrtrc(size_t size, mpw* result, const mpw* data)
797 {
798 #if HAVE_MPDW
799 register mpdw temp;
800 register mpw load, carry = 0;
801
802 result += (size << 1);
803
804 while (size--)
805 {
806 temp = load = data[size];
807 temp *= load;
808 temp += carry;
809 temp += *(--result);
810 *result = (mpw) temp;
811 temp >>= MP_WBITS;
812 temp += *(--result);
813 *result = (mpw) temp;
814 carry = (mpw)(temp >> MP_WBITS);
815 }
816 #else
817 register mpw temp, load, carry = 0;
818
819 result += (size << 1);
820
821 while (size--)
822 {
823 register mphw xlo, xhi;
824 register mpw rlo, rhi;
825
826 xlo = (mphw) (temp = data[size]);
827 xhi = (mphw) (temp >> MP_HWBITS);
828
829 rlo = (mpw) xlo * xlo;
830 rhi = (mpw) xhi * xhi;
831 temp = (mpw) xhi * xlo;
832 load = rlo;
833 rlo += (temp << MP_HWBITS);
834 rhi += (temp >> MP_HWBITS) + (load > rlo);
835 load = rlo;
836 rlo += (temp << MP_HWBITS);
837 rhi += (temp >> MP_HWBITS) + (load > rlo);
838 load = rlo;
839 rlo += carry;
840 rhi += (load > rlo);
841 load = rlo;
842 rlo += *(--result);
843 *result = rlo;
844 temp = (load > rlo);
845 load = rhi;
846 rhi += temp;
847 carry = (load > rhi);
848 load = rhi;
849 rhi += *(--result);
850 *result = rhi;
851 carry += (load > rhi);
852 }
853 #endif
854 }
855 #endif
856
857 #ifndef ASM_MPSQR
mpsqr(mpw * result,size_t size,const mpw * data)858 void mpsqr(mpw* result, size_t size, const mpw* data)
859 {
860 register mpw rc;
861 register size_t n = size-1;
862
863 result += size;
864 result[n] = 0;
865
866 if (n)
867 {
868 rc = mpsetmul(n, result, data, data[n]);
869 *(--result) = rc;
870 while (--n)
871 {
872 rc = mpaddmul(n, result, data, data[n]);
873 *(--result) = rc;
874 }
875 }
876
877 *(--result) = 0;
878
879 mpmultwo(size << 1, result);
880
881 mpaddsqrtrc(size, result, data);
882 }
883 #endif
884
885 #ifndef ASM_MPSIZE
mpsize(size_t size,const mpw * data)886 size_t mpsize(size_t size, const mpw* data)
887 {
888 while (size)
889 {
890 if (*data)
891 return size;
892 data++;
893 size--;
894 }
895 return 0;
896 }
897 #endif
898
899 #ifndef ASM_MPBITS
mpbits(size_t size,const mpw * data)900 size_t mpbits(size_t size, const mpw* data)
901 {
902 return MP_WORDS_TO_BITS(size) - mpmszcnt(size, data);
903 }
904 #endif
905
906 #ifndef ASM_MPNORM
mpnorm(size_t size,mpw * data)907 size_t mpnorm(size_t size, mpw* data)
908 {
909 register size_t shift = mpmszcnt(size, data);
910 mplshift(size, data, shift);
911 return shift;
912 }
913 #endif
914
915 #ifndef ASM_MPDIVTWO
mpdivtwo(size_t size,mpw * data)916 void mpdivtwo(size_t size, mpw* data)
917 {
918 register mpw temp, carry = 0;
919
920 while (size--)
921 {
922 temp = *data;
923 *(data++) = (temp >> 1) | carry;
924 carry = (temp << (MP_WBITS-1));
925 }
926 }
927 #endif
928
929 #ifndef ASM_MPSDIVTWO
mpsdivtwo(size_t size,mpw * data)930 void mpsdivtwo(size_t size, mpw* data)
931 {
932 int carry = mpmsbset(size, data);
933 mpdivtwo(size, data);
934 if (carry)
935 mpsetmsb(size, data);
936 }
937 #endif
938
939 #ifndef ASM_MPMULTWO
mpmultwo(size_t size,mpw * data)940 int mpmultwo(size_t size, mpw* data)
941 {
942 register mpw temp, carry = 0;
943
944 data += size;
945 while (size--)
946 {
947 temp = *(--data);
948 *data = (temp << 1) | carry;
949 carry = (temp >> (MP_WBITS-1));
950 }
951 return (int) carry;
952 }
953 #endif
954
955 #ifndef ASM_MPMSZCNT
mpmszcnt(size_t size,const mpw * data)956 size_t mpmszcnt(size_t size, const mpw* data)
957 {
958 register size_t zbits = 0;
959 register size_t i = 0;
960
961 while (i < size)
962 {
963 register mpw temp = data[i++];
964 if (temp)
965 {
966 while (!(temp & MP_MSBMASK))
967 {
968 zbits++;
969 temp <<= 1;
970 }
971 break;
972 }
973 else
974 zbits += MP_WBITS;
975 }
976 return zbits;
977 }
978 #endif
979
980 #ifndef ASM_MPLSZCNT
mplszcnt(size_t size,const mpw * data)981 size_t mplszcnt(size_t size, const mpw* data)
982 {
983 register size_t zbits = 0;
984
985 while (size--)
986 {
987 register mpw temp = data[size];
988 if (temp)
989 {
990 while (!(temp & MP_LSBMASK))
991 {
992 zbits++;
993 temp >>= 1;
994 }
995 break;
996 }
997 else
998 zbits += MP_WBITS;
999 }
1000 return zbits;
1001 }
1002 #endif
1003
1004 #ifndef ASM_MPLSHIFT
mplshift(size_t size,mpw * data,size_t count)1005 void mplshift(size_t size, mpw* data, size_t count)
1006 {
1007 register size_t words = MP_BITS_TO_WORDS(count);
1008
1009 if (words < size)
1010 {
1011 register short lbits = (short) (count & (MP_WBITS-1));
1012
1013 /* first do the shifting, then do the moving */
1014 if (lbits)
1015 {
1016 register mpw temp, carry = 0;
1017 register short rbits = MP_WBITS - lbits;
1018 register size_t i = size;
1019
1020 while (i > words)
1021 {
1022 temp = data[--i];
1023 data[i] = (temp << lbits) | carry;
1024 carry = (temp >> rbits);
1025 }
1026 }
1027 if (words)
1028 {
1029 mpmove(size-words, data, data+words);
1030 mpzero(words, data+size-words);
1031 }
1032 }
1033 else
1034 mpzero(size, data);
1035 }
1036 #endif
1037
1038 #ifndef ASM_MPRSHIFT
mprshift(size_t size,mpw * data,size_t count)1039 void mprshift(size_t size, mpw* data, size_t count)
1040 {
1041 register size_t words = MP_BITS_TO_WORDS(count);
1042
1043 if (words < size)
1044 {
1045 register short rbits = (short) (count & (MP_WBITS-1));
1046
1047 /* first do the shifting, then do the moving */
1048 if (rbits)
1049 {
1050 register mpw temp, carry = 0;
1051 register short lbits = MP_WBITS - rbits;
1052 register size_t i = 0;
1053
1054 while (i < size-words)
1055 {
1056 temp = data[i];
1057 data[i++] = (temp >> rbits) | carry;
1058 carry = (temp << lbits);
1059 }
1060 }
1061 if (words)
1062 {
1063 mpmove(size-words, data+words, data);
1064 mpzero(words, data);
1065 }
1066 }
1067 else
1068 mpzero(size, data);
1069 }
1070 #endif
1071
1072 #ifndef ASM_MPRSHIFTLSZ
mprshiftlsz(size_t size,mpw * data)1073 size_t mprshiftlsz(size_t size, mpw* data)
1074 {
1075 register mpw* slide = data+size-1;
1076 register size_t zwords = 0; /* counter for 'all zero bit' words */
1077 register short lbits, rbits = 0; /* counter for 'least significant zero' bits */
1078 register mpw temp, carry = 0;
1079
1080 data = slide;
1081
1082 /* count 'all zero' words and move src pointer */
1083 while (size--)
1084 {
1085 /* test if we have a non-zero word */
1086 if ((carry = *(slide--)))
1087 {
1088 /* count 'least signification zero bits and set zbits counter */
1089 while (!(carry & MP_LSBMASK))
1090 {
1091 carry >>= 1;
1092 rbits++;
1093 }
1094 break;
1095 }
1096 zwords++;
1097 }
1098
1099 if ((rbits == 0) && (zwords == 0))
1100 return 0;
1101
1102 /* prepare right-shifting of data */
1103 lbits = MP_WBITS - rbits;
1104
1105 /* shift data */
1106 while (size--)
1107 {
1108 temp = *(slide--);
1109 *(data--) = (temp << lbits) | carry;
1110 carry = (temp >> rbits);
1111 }
1112
1113 /* store the final carry */
1114 *(data--) = carry;
1115
1116 /* store the return value in size */
1117 size = MP_WORDS_TO_BITS(zwords) + rbits;
1118
1119 /* zero the (zwords) most significant words */
1120 while (zwords--)
1121 *(data--) = 0;
1122
1123 return size;
1124 }
1125 #endif
1126
1127 /* try an alternate version here, with descending sizes */
1128 /* also integrate lszcnt and rshift properly into one function */
1129 #ifndef ASM_MPGCD_W
1130 /*
1131 * mpgcd_w
1132 * need workspace of (size) words
1133 */
mpgcd_w(size_t size,const mpw * xdata,const mpw * ydata,mpw * result,mpw * wksp)1134 void mpgcd_w(size_t size, const mpw* xdata, const mpw* ydata, mpw* result, mpw* wksp)
1135 {
1136 register size_t shift, temp;
1137
1138 if (mpge(size, xdata, ydata))
1139 {
1140 mpcopy(size, wksp, xdata);
1141 mpcopy(size, result, ydata);
1142 }
1143 else
1144 {
1145 mpcopy(size, wksp, ydata);
1146 mpcopy(size, result, xdata);
1147 }
1148
1149 /* get the smallest returned values, and set shift to that */
1150
1151 shift = mprshiftlsz(size, wksp);
1152 temp = mprshiftlsz(size, result);
1153
1154 if (shift > temp)
1155 shift = temp;
1156
1157 while (mpnz(size, wksp))
1158 {
1159 mprshiftlsz(size, wksp);
1160 mprshiftlsz(size, result);
1161
1162 if (mpge(size, wksp, result))
1163 mpsub(size, wksp, result);
1164 else
1165 mpsub(size, result, wksp);
1166
1167 /* slide past zero words in both operands by increasing pointers and decreasing size */
1168 if ((*wksp == 0) && (*result == 0))
1169 {
1170 size--;
1171 wksp++;
1172 result++;
1173 }
1174 }
1175
1176 /* figure out if we need to slide the result pointer back */
1177 if ((temp = MP_BITS_TO_WORDS(shift)))
1178 {
1179 size += temp;
1180 result -= temp;
1181 }
1182
1183 mplshift(size, result, shift);
1184 }
1185 #endif
1186
1187 #ifndef ASM_MPEXTGCD_W
1188 /* needs workspace of (6*size+6) words */
1189 /* used to compute the modular inverse */
mpextgcd_w(size_t size,const mpw * xdata,const mpw * ydata,mpw * result,mpw * wksp)1190 int mpextgcd_w(size_t size, const mpw* xdata, const mpw* ydata, mpw* result, mpw* wksp)
1191 {
1192 /*
1193 * For computing a modular inverse, pass the modulus as xdata and the number
1194 * to be inverted as ydata.
1195 *
1196 * Fact: if a element of Zn, then a is invertible if and only if gcd(a,n) = 1
1197 * Hence: if n is even, then a must be odd, otherwise the gcd(a,n) >= 2
1198 *
1199 * The calling routine must guarantee this condition.
1200 */
1201
1202 register size_t sizep = size+1;
1203 register int full;
1204
1205 mpw* udata = wksp;
1206 mpw* vdata = udata+sizep;
1207 mpw* adata = vdata+sizep;
1208 mpw* bdata = adata+sizep;
1209 mpw* cdata = bdata+sizep;
1210 mpw* ddata = cdata+sizep;
1211
1212 mpsetx(sizep, udata, size, xdata);
1213 mpsetx(sizep, vdata, size, ydata);
1214 mpzero(sizep, bdata);
1215 mpsetw(sizep, ddata, 1);
1216
1217 if ((full = mpeven(sizep, udata)))
1218 {
1219 mpsetw(sizep, adata, 1);
1220 mpzero(sizep, cdata);
1221 }
1222
1223 while (1)
1224 {
1225 while (mpeven(sizep, udata))
1226 {
1227 mpdivtwo(sizep, udata);
1228
1229 if (mpodd(sizep, bdata) || (full && mpodd(sizep, adata)))
1230 {
1231 if (full) mpaddx(sizep, adata, size, ydata);
1232 mpsubx(sizep, bdata, size, xdata);
1233 }
1234
1235 if (full) mpsdivtwo(sizep, adata);
1236 mpsdivtwo(sizep, bdata);
1237 }
1238 while (mpeven(sizep, vdata))
1239 {
1240 mpdivtwo(sizep, vdata);
1241
1242 if (mpodd(sizep, ddata) || (full && mpodd(sizep, cdata)))
1243 {
1244 if (full) mpaddx(sizep, cdata, size, ydata);
1245 mpsubx(sizep, ddata, size, xdata);
1246 }
1247
1248 if (full) mpsdivtwo(sizep, cdata);
1249 mpsdivtwo(sizep, ddata);
1250 }
1251 if (mpge(sizep, udata, vdata))
1252 {
1253 mpsub(sizep, udata, vdata);
1254 if (full) mpsub(sizep, adata, cdata);
1255 mpsub(sizep, bdata, ddata);
1256 }
1257 else
1258 {
1259 mpsub(sizep, vdata, udata);
1260 if (full) mpsub(sizep, cdata, adata);
1261 mpsub(sizep, ddata, bdata);
1262 }
1263 if (mpz(sizep, udata))
1264 {
1265 if (mpisone(sizep, vdata))
1266 {
1267 if (result)
1268 {
1269 if (*ddata & MP_MSBMASK)
1270 {
1271 /* keep adding the modulus until we get a carry */
1272 while (!mpaddx(sizep, ddata, size, xdata));
1273 }
1274 else
1275 {
1276 /* in some computations, d ends up > x, hence:
1277 * keep subtracting n from d until d < x
1278 */
1279 while (mpgtx(sizep, ddata, size, xdata))
1280 mpsubx(sizep, ddata, size, xdata);
1281 }
1282 mpsetx(size, result, sizep, ddata);
1283 }
1284 return 1;
1285 }
1286 return 0;
1287 }
1288 }
1289 }
1290 #endif
1291
1292 #ifndef ASM_MPPNDIV
mppndiv(mpw xhi,mpw xlo,mpw y)1293 mpw mppndiv(mpw xhi, mpw xlo, mpw y)
1294 {
1295 register mpw result = 0;
1296 register short count = MP_WBITS;
1297 register int carry = 0;
1298
1299 while (count--)
1300 {
1301 if (carry | (xhi >= y))
1302 {
1303 xhi -= y;
1304 result++;
1305 }
1306 carry = (xhi >> (MP_WBITS-1));
1307 xhi <<= 1;
1308 xhi |= (xlo >> (MP_WBITS-1));
1309 xlo <<= 1;
1310 result <<= 1;
1311 }
1312 if (carry | (xhi >= y))
1313 {
1314 xhi -= y;
1315 result++;
1316 }
1317 return result;
1318 }
1319 #endif
1320
1321 #ifndef ASM_MPMOD
mpmod(mpw * result,size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata,mpw * workspace)1322 void mpmod(mpw* result, size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata, mpw* workspace)
1323 {
1324 /* result size ysize, workspace size 2*ysize+1 */
1325 mpw q, msw;
1326 mpw* rdata = result;
1327 mpw* ynorm = workspace+ysize+1;
1328 size_t shift, qsize = xsize-ysize;
1329
1330 mpcopy(ysize, ynorm, ydata);
1331 shift = mpnorm(ysize, ynorm);
1332 msw = *ynorm;
1333 mpcopy(xsize, rdata, xdata);
1334 if (mpge(ysize, rdata, ynorm))
1335 mpsub(ysize, rdata, ynorm);
1336
1337 while (qsize--)
1338 {
1339 q = mppndiv(rdata[0], rdata[1], msw);
1340
1341 *workspace = mpsetmul(ysize, workspace+1, ynorm, q);
1342
1343 while (mplt(ysize+1, rdata, workspace))
1344 {
1345 mpsubx(ysize+1, workspace, ysize, ynorm);
1346 q--;
1347 }
1348 mpsub(ysize+1, rdata, workspace);
1349 rdata++;
1350 }
1351 /* de-normalization steps */
1352 while (shift--)
1353 {
1354 mpdivtwo(ysize, ynorm);
1355 if (mpge(ysize, rdata, ynorm))
1356 mpsub(ysize, rdata, ynorm);
1357 }
1358 }
1359 #endif
1360
1361 #ifndef ASM_MPNDIVMOD
mpndivmod(mpw * result,size_t xsize,const mpw * xdata,size_t ysize,const mpw * ydata,register mpw * workspace)1362 void mpndivmod(mpw* result, size_t xsize, const mpw* xdata, size_t ysize, const mpw* ydata, register mpw* workspace)
1363 {
1364 /* result must be xsize+1 in length */
1365 /* workspace must be ysize+1 in length */
1366 /* expect ydata to be normalized */
1367 mpw q;
1368 mpw msw = *ydata;
1369 size_t qsize = xsize-ysize;
1370
1371 *result = (mpge(ysize, xdata, ydata) ? 1 : 0);
1372 mpcopy(xsize, result+1, xdata);
1373
1374 if (*result)
1375 (void) mpsub(ysize, result+1, ydata);
1376
1377 result++;
1378
1379 while (qsize--)
1380 {
1381 q = mppndiv(result[0], result[1], msw);
1382
1383 *workspace = mpsetmul(ysize, workspace+1, ydata, q);
1384
1385 while (mplt(ysize+1, result, workspace))
1386 {
1387 mpsubx(ysize+1, workspace, ysize, ydata);
1388 q--;
1389 }
1390 mpsub(ysize+1, result, workspace);
1391 *(result++) = q;
1392 }
1393 }
1394 #endif
1395
mpprint(size_t size,const mpw * data)1396 void mpprint(size_t size, const mpw* data)
1397 {
1398 mpfprint(stdout, size, data);
1399 }
1400
mpprintln(size_t size,const mpw * data)1401 void mpprintln(size_t size, const mpw* data)
1402 {
1403 mpfprintln(stdout, size, data);
1404 }
1405
mpfprint(FILE * f,size_t size,const mpw * data)1406 void mpfprint(FILE* f, size_t size, const mpw* data)
1407 {
1408 if (data == (mpw*) 0)
1409 return;
1410
1411 if (f == (FILE*) 0)
1412 return;
1413
1414 while (size--)
1415 {
1416 #if (MP_WBITS == 32)
1417 fprintf(f, "%08x", (unsigned) *(data++));
1418 #elif (MP_WBITS == 64)
1419 # if WIN32
1420 fprintf(f, "%016I64x", *(data++));
1421 # elif SIZEOF_UNSIGNED_LONG == 8
1422 fprintf(f, "%016lx", *(data++));
1423 # else
1424 fprintf(f, "%016llx", *(data++));
1425 # endif
1426 #else
1427 # error
1428 #endif
1429 }
1430 fflush(f);
1431 }
1432
mpfprintln(FILE * f,size_t size,const mpw * data)1433 void mpfprintln(FILE* f, size_t size, const mpw* data)
1434 {
1435 if (data == (mpw*) 0)
1436 return;
1437
1438 if (f == (FILE*) 0)
1439 return;
1440
1441 while (size--)
1442 {
1443 #if (MP_WBITS == 32)
1444 fprintf(f, "%08x", *(data++));
1445 #elif (MP_WBITS == 64)
1446 # if WIN32
1447 fprintf(f, "%016I64x", *(data++));
1448 # elif SIZEOF_UNSIGNED_LONG == 8
1449 fprintf(f, "%016lx", *(data++));
1450 # else
1451 fprintf(f, "%016llx", *(data++));
1452 # endif
1453 #else
1454 # error
1455 #endif
1456 }
1457 fprintf(f, "\n");
1458 fflush(f);
1459 }
1460
i2osp(byte * osdata,size_t ossize,const mpw * idata,size_t isize)1461 int i2osp(byte *osdata, size_t ossize, const mpw* idata, size_t isize)
1462 {
1463 #if WORDS_BIGENDIAN
1464 size_t max_bytes = MP_WORDS_TO_BYTES(isize);
1465 #endif
1466 size_t significant_bytes = (mpbits(isize, idata) + 7) >> 3;
1467
1468 /* verify that ossize is large enough to contain the significant bytes */
1469 if (ossize >= significant_bytes)
1470 {
1471 /* looking good; check if we have more space than significant bytes */
1472 if (ossize > significant_bytes)
1473 { /* fill most significant bytes with zero */
1474 memset(osdata, 0, ossize - significant_bytes);
1475 osdata += ossize - significant_bytes;
1476 }
1477 if (significant_bytes)
1478 { /* fill remaining bytes with endian-adjusted data */
1479 #if !WORDS_BIGENDIAN
1480 mpw w = idata[--isize];
1481 byte shift = 0;
1482
1483 /* fill right-to-left; much easier than left-to-right */
1484 do
1485 {
1486 osdata[--significant_bytes] = (byte)(w >> shift);
1487 shift += 8;
1488 if ((shift == MP_WBITS) && isize)
1489 {
1490 shift = 0;
1491 w = idata[--isize];
1492 }
1493 } while (significant_bytes);
1494 #else
1495 /* just copy data past zero bytes */
1496 memcpy(osdata, ((byte*) idata) + (max_bytes - significant_bytes), significant_bytes);
1497 #endif
1498 }
1499 return 0;
1500 }
1501 return -1;
1502 }
1503
os2ip(mpw * idata,size_t isize,const byte * osdata,size_t ossize)1504 int os2ip(mpw* idata, size_t isize, const byte* osdata, size_t ossize)
1505 {
1506 size_t required;
1507
1508 /* skip non-significant leading zero bytes */
1509 while (!(*osdata) && ossize)
1510 {
1511 osdata++;
1512 ossize--;
1513 }
1514
1515 required = MP_BYTES_TO_WORDS(ossize + MP_WBYTES - 1);
1516
1517 if (isize >= required)
1518 {
1519 /* yes, we have enough space and can proceed */
1520 mpw w = 0;
1521 /* adjust counter so that the loop will start by skipping the proper
1522 * amount of leading bytes in the first significant word
1523 */
1524 byte b = (byte) (ossize % MP_WBYTES);
1525
1526 if (isize > required)
1527 { /* fill initials words with zero */
1528 mpzero(isize-required, idata);
1529 idata += isize-required;
1530 }
1531
1532 if (b == 0)
1533 b = MP_WBYTES;
1534
1535 while (ossize--)
1536 {
1537 w <<= 8;
1538 w |= *(osdata++);
1539 b--;
1540
1541 if (b == 0)
1542 {
1543 *(idata++) = w;
1544 w = 0;
1545 b = MP_WBYTES;
1546 }
1547 }
1548
1549 return 0;
1550 }
1551 return -1;
1552 }
1553
hs2ip(mpw * idata,size_t isize,const char * hsdata,size_t hssize)1554 int hs2ip(mpw* idata, size_t isize, const char* hsdata, size_t hssize)
1555 {
1556 size_t required = MP_NIBBLES_TO_WORDS(hssize + MP_WNIBBLES - 1);
1557
1558 if (isize >= required)
1559 {
1560 register size_t i;
1561
1562
1563 if (isize > required)
1564 { /* fill initial words with zero */
1565 for (i = required; i < isize; i++)
1566 *(idata++) = 0;
1567 }
1568 while (hssize)
1569 {
1570 register mpw w = 0;
1571 register size_t chunk = hssize & (MP_WNIBBLES - 1);
1572 register char ch;
1573
1574 if (chunk == 0) chunk = MP_WNIBBLES;
1575
1576 for (i = 0; i < chunk; i++)
1577 {
1578 ch = *(hsdata++);
1579 w <<= 4;
1580 if (ch >= '0' && ch <= '9')
1581 w += (ch - '0');
1582 else if (ch >= 'A' && ch <= 'F')
1583 w += (ch - 'A') + 10;
1584 else if (ch >= 'a' && ch <= 'f')
1585 w += (ch - 'a') + 10;
1586 }
1587 *(idata++) = w;
1588 hssize -= chunk;
1589 }
1590 return 0;
1591 }
1592 return -1;
1593 }
1594