1 /*
2 * Copyright (C) 2017 Denys Vlasenko
3 *
4 * Licensed under GPLv2, see file LICENSE in this source tree.
5 */
6 #include "tls.h"
7
8 /* The file is taken almost verbatim from matrixssl-3-7-2b-open/crypto/math/.
9 * Changes are flagged with //bbox
10 */
11
12 /**
13 * @file pstm.c
14 * @version 33ef80f (HEAD, tag: MATRIXSSL-3-7-2-OPEN, tag: MATRIXSSL-3-7-2-COMM, origin/master, origin/HEAD, master)
15 *
16 * Multiprecision number implementation.
17 */
18 /*
19 * Copyright (c) 2013-2015 INSIDE Secure Corporation
20 * Copyright (c) PeerSec Networks, 2002-2011
21 * All Rights Reserved
22 *
23 * The latest version of this code is available at http://www.matrixssl.org
24 *
25 * This software is open source; you can redistribute it and/or modify
26 * it under the terms of the GNU General Public License as published by
27 * the Free Software Foundation; either version 2 of the License, or
28 * (at your option) any later version.
29 *
30 * This General Public License does NOT permit incorporating this software
31 * into proprietary programs. If you are unable to comply with the GPL, a
32 * commercial license for this software may be purchased from INSIDE at
33 * http://www.insidesecure.com/eng/Company/Locations
34 *
35 * This program is distributed in WITHOUT ANY WARRANTY; without even the
36 * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
37 * See the GNU General Public License for more details.
38 *
39 * You should have received a copy of the GNU General Public License
40 * along with this program; if not, write to the Free Software
41 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
42 * http://www.gnu.org/copyleft/gpl.html
43 */
44 /******************************************************************************/
45
46 //bbox
47 //#include "../cryptoApi.h"
48 #ifndef DISABLE_PSTM
49
50 #undef pstm_mul_2d
51 static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c); //bbox: was int16 b
52 #define pstm_mul_2d(a, b, c) (pstm_mul_2d(a, b, c), PSTM_OKAY)
53
54 /******************************************************************************/
55 /*
56 init an pstm_int for a given size
57 */
58 #undef pstm_init_size
59 #define pstm_init_size(pool, a, size) \
60 pstm_init_size( a, size)
pstm_init_size(psPool_t * pool,pstm_int * a,uint32 size)61 int32 FAST_FUNC pstm_init_size(psPool_t *pool, pstm_int * a, uint32 size)
62 {
63 //bbox
64 // uint16 x;
65
66 /*
67 alloc mem
68 */
69 a->dp = xzalloc(sizeof (pstm_digit) * size);//bbox
70 //bbox a->pool = pool;
71 a->used = 0;
72 a->alloc = size;
73 a->sign = PSTM_ZPOS;
74 /*
75 zero the digits
76 */
77 //bbox
78 // for (x = 0; x < size; x++) {
79 // a->dp[x] = 0;
80 // }
81 return PSTM_OKAY;
82 }
83 #undef pstm_init_size
84 #define pstm_init_size(pool, a, size) (pstm_init_size(a, size), PSTM_OKAY)
85
86 /******************************************************************************/
87 /*
88 Init a new pstm_int.
89 */
90 #undef pstm_init
91 #define pstm_init(pool, a) \
92 pstm_init( a)
pstm_init(psPool_t * pool,pstm_int * a)93 static int32 pstm_init(psPool_t *pool, pstm_int * a)
94 {
95 //bbox
96 // int32 i;
97 /*
98 allocate memory required and clear it
99 */
100 a->dp = xzalloc(sizeof (pstm_digit) * PSTM_DEFAULT_INIT);//bbox
101 /*
102 set the digits to zero
103 */
104 //bbox
105 // for (i = 0; i < PSTM_DEFAULT_INIT; i++) {
106 // a->dp[i] = 0;
107 // }
108 /*
109 set the used to zero, allocated digits to the default precision and sign
110 to positive
111 */
112 //bbox a->pool = pool;
113 a->used = 0;
114 a->alloc = PSTM_DEFAULT_INIT;
115 a->sign = PSTM_ZPOS;
116
117 return PSTM_OKAY;
118 }
119 #undef pstm_init
120 #define pstm_init(pool, a) (pstm_init(a), PSTM_OKAY)
121
122 /******************************************************************************/
123 /*
124 Grow as required
125 */
126 #undef pstm_grow
pstm_grow(pstm_int * a,int size)127 int32 FAST_FUNC pstm_grow(pstm_int * a, int size)
128 {
129 int i; //bbox: was int16
130 pstm_digit *tmp;
131
132 /*
133 If the alloc size is smaller alloc more ram.
134 */
135 if (a->alloc < size) {
136 /*
137 Reallocate the array a->dp
138
139 We store the return in a temporary variable in case the operation
140 failed we don't want to overwrite the dp member of a.
141 */
142 tmp = xrealloc(a->dp, sizeof (pstm_digit) * size);//bbox
143 /*
144 reallocation succeeded so set a->dp
145 */
146 a->dp = tmp;
147 /*
148 zero excess digits
149 */
150 i = a->alloc;
151 a->alloc = size;
152 for (; i < a->alloc; i++) {
153 a->dp[i] = 0;
154 }
155 }
156 return PSTM_OKAY;
157 }
158 #define pstm_grow(a, size) (pstm_grow(a, size), PSTM_OKAY)
159
160 /******************************************************************************/
161 /*
162 copy, b = a (b must be pre-allocated)
163 */
164 #undef pstm_copy
pstm_copy(pstm_int * a,pstm_int * b)165 int32 pstm_copy(pstm_int * a, pstm_int * b)
166 {
167 int32 res, n;
168
169 /*
170 If dst == src do nothing
171 */
172 if (a == b) {
173 return PSTM_OKAY;
174 }
175 /*
176 Grow dest
177 */
178 if (b->alloc < a->used) {
179 if ((res = pstm_grow (b, a->used)) != PSTM_OKAY) {
180 return res;
181 }
182 }
183 /*
184 Zero b and copy the parameters over
185 */
186 {
187 register pstm_digit *tmpa, *tmpb;
188
189 /* pointer aliases */
190 /* source */
191 tmpa = a->dp;
192
193 /* destination */
194 tmpb = b->dp;
195
196 /* copy all the digits */
197 for (n = 0; n < a->used; n++) {
198 *tmpb++ = *tmpa++;
199 }
200
201 /* clear high digits */
202 for (; n < b->used; n++) {
203 *tmpb++ = 0;
204 }
205 }
206 /*
207 copy used count and sign
208 */
209 b->used = a->used;
210 b->sign = a->sign;
211 return PSTM_OKAY;
212 }
213 #define pstm_copy(a, b) (pstm_copy(a, b), PSTM_OKAY)
214
215 /******************************************************************************/
216 /*
217 Trim unused digits
218
219 This is used to ensure that leading zero digits are trimed and the
220 leading "used" digit will be non-zero. Typically very fast. Also fixes
221 the sign if there are no more leading digits
222 */
pstm_clamp(pstm_int * a)223 void FAST_FUNC pstm_clamp(pstm_int * a)
224 {
225 /* decrease used while the most significant digit is zero. */
226 while (a->used > 0 && a->dp[a->used - 1] == 0) {
227 --(a->used);
228 }
229 /* reset the sign flag if used == 0 */
230 if (a->used == 0) {
231 a->sign = PSTM_ZPOS;
232 }
233 }
234
235 /******************************************************************************/
236 /*
237 clear one (frees).
238 */
pstm_clear(pstm_int * a)239 void FAST_FUNC pstm_clear(pstm_int * a)
240 {
241 int32 i;
242 /*
243 only do anything if a hasn't been freed previously
244 */
245 if (a != NULL && a->dp != NULL) {
246 /*
247 first zero the digits
248 */
249 for (i = 0; i < a->used; i++) {
250 a->dp[i] = 0;
251 }
252
253 psFree (a->dp, a->pool);
254 /*
255 reset members to make debugging easier
256 */
257 a->dp = NULL;
258 a->alloc = a->used = 0;
259 a->sign = PSTM_ZPOS;
260 }
261 }
262
263 /******************************************************************************/
264 /*
265 clear many (frees).
266 */
267 #if 0 //UNUSED
268 void pstm_clear_multi(pstm_int *mp0, pstm_int *mp1, pstm_int *mp2,
269 pstm_int *mp3, pstm_int *mp4, pstm_int *mp5,
270 pstm_int *mp6, pstm_int *mp7)
271 {
272 int32 n; /* Number of ok inits */
273
274 pstm_int *tempArray[9];
275
276 tempArray[0] = mp0;
277 tempArray[1] = mp1;
278 tempArray[2] = mp2;
279 tempArray[3] = mp3;
280 tempArray[4] = mp4;
281 tempArray[5] = mp5;
282 tempArray[6] = mp6;
283 tempArray[7] = mp7;
284 tempArray[8] = NULL;
285
286 for (n = 0; tempArray[n] != NULL; n++) {
287 if ((tempArray[n] != NULL) && (tempArray[n]->dp != NULL)) {
288 pstm_clear(tempArray[n]);
289 }
290 }
291 }
292 #endif
293
294 /******************************************************************************/
295 /*
296 Set to zero.
297 */
pstm_zero(pstm_int * a)298 static void pstm_zero(pstm_int * a)
299 {
300 int32 n;
301 pstm_digit *tmp;
302
303 a->sign = PSTM_ZPOS;
304 a->used = 0;
305
306 tmp = a->dp;
307 for (n = 0; n < a->alloc; n++) {
308 *tmp++ = 0;
309 }
310 }
311
312
313 /******************************************************************************/
314 /*
315 Compare maginitude of two ints (unsigned).
316 */
pstm_cmp_mag(pstm_int * a,pstm_int * b)317 int32 FAST_FUNC pstm_cmp_mag(pstm_int * a, pstm_int * b)
318 {
319 int n; //bbox: was int16
320 pstm_digit *tmpa, *tmpb;
321
322 /*
323 compare based on # of non-zero digits
324 */
325 if (a->used > b->used) {
326 return PSTM_GT;
327 }
328
329 if (a->used < b->used) {
330 return PSTM_LT;
331 }
332
333 /* alias for a */
334 tmpa = a->dp + (a->used - 1);
335
336 /* alias for b */
337 tmpb = b->dp + (a->used - 1);
338
339 /*
340 compare based on digits
341 */
342 for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
343 if (*tmpa > *tmpb) {
344 return PSTM_GT;
345 }
346 if (*tmpa < *tmpb) {
347 return PSTM_LT;
348 }
349 }
350 return PSTM_EQ;
351 }
352
353 /******************************************************************************/
354 /*
355 Compare two ints (signed)
356 */
pstm_cmp(pstm_int * a,pstm_int * b)357 int32 FAST_FUNC pstm_cmp(pstm_int * a, pstm_int * b)
358 {
359 /*
360 compare based on sign
361 */
362 if (a->sign != b->sign) {
363 if (a->sign == PSTM_NEG) {
364 return PSTM_LT;
365 } else {
366 return PSTM_GT;
367 }
368 }
369 /*
370 compare digits
371 */
372 if (a->sign == PSTM_NEG) {
373 /* if negative compare opposite direction */
374 return pstm_cmp_mag(b, a);
375 } else {
376 return pstm_cmp_mag(a, b);
377 }
378 }
379
380 /******************************************************************************/
381 /*
382 pstm_ints can be initialized more precisely when they will populated
383 using pstm_read_unsigned_bin since the length of the byte stream is known
384 */
pstm_init_for_read_unsigned_bin(psPool_t * pool,pstm_int * a,uint32 len)385 int32 FAST_FUNC pstm_init_for_read_unsigned_bin(psPool_t *pool, pstm_int *a, uint32 len)
386 {
387 int32 size;
388 /*
389 Need to set this based on how many words max it will take to store the bin.
390 The magic + 2:
391 1 to round up for the remainder of this integer math
392 1 for the initial carry of '1' bits that fall between DIGIT_BIT and 8
393 */
394 size = (((len / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
395 / DIGIT_BIT) + 2;
396 return pstm_init_size(pool, a, size);
397 }
398
399
400 /******************************************************************************/
401 /*
402 Reads a unsigned char array into pstm_int format. User should have
403 called pstm_init_for_read_unsigned_bin first. There is some grow logic
404 here if the default pstm_init was used but we don't really want to hit it.
405 */
pstm_read_unsigned_bin(pstm_int * a,unsigned char * b,int32 c)406 int32 FAST_FUNC pstm_read_unsigned_bin(pstm_int *a, unsigned char *b, int32 c)
407 {
408 /* zero the int */
409 pstm_zero (a);
410
411 /*
412 If we know the endianness of this architecture, and we're using
413 32-bit pstm_digits, we can optimize this
414 */
415 #if (defined(ENDIAN_LITTLE) || defined(ENDIAN_BIG)) && !defined(PSTM_64BIT)
416 /* But not for both simultaneously */
417 #if defined(ENDIAN_LITTLE) && defined(ENDIAN_BIG)
418 #error Both ENDIAN_LITTLE and ENDIAN_BIG defined.
419 #endif
420 {
421 unsigned char *pd;
422 if ((unsigned)c > (PSTM_MAX_SIZE * sizeof(pstm_digit))) {
423 uint32 excess = c - (PSTM_MAX_SIZE * sizeof(pstm_digit));
424 c -= excess;
425 b += excess;
426 }
427 a->used = ((c + sizeof(pstm_digit) - 1)/sizeof(pstm_digit));
428 if (a->alloc < a->used) {
429 if (pstm_grow(a, a->used) != PSTM_OKAY) {
430 return PSTM_MEM;
431 }
432 }
433 pd = (unsigned char *)a->dp;
434 /* read the bytes in */
435 #ifdef ENDIAN_BIG
436 {
437 /* Use Duff's device to unroll the loop. */
438 int32 idx = (c - 1) & ~3;
439 switch (c % 4) {
440 case 0: do { pd[idx+0] = *b++;
441 case 3: pd[idx+1] = *b++;
442 case 2: pd[idx+2] = *b++;
443 case 1: pd[idx+3] = *b++;
444 idx -= 4;
445 } while ((c -= 4) > 0);
446 }
447 }
448 #else
449 for (c -= 1; c >= 0; c -= 1) {
450 pd[c] = *b++;
451 }
452 #endif
453 }
454 #else
455 /* Big enough based on the len? */
456 a->used = (((c / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
457 / DIGIT_BIT) + 2;
458
459 if (a->alloc < a->used) {
460 if (pstm_grow(a, a->used) != PSTM_OKAY) {
461 return PSTM_MEM;
462 }
463 }
464 /* read the bytes in */
465 for (; c > 0; c--) {
466 if (pstm_mul_2d (a, 8, a) != PSTM_OKAY) {
467 return PS_MEM_FAIL;
468 }
469 a->dp[0] |= *b++;
470 a->used += 1;
471 }
472 #endif
473
474 pstm_clamp (a);
475 return PS_SUCCESS;
476 }
477
478 /******************************************************************************/
479 /*
480 */
pstm_count_bits(pstm_int * a)481 static int pstm_count_bits(pstm_int * a)
482 {
483 int r; //bbox: was int16
484 pstm_digit q;
485
486 if (a->used == 0) {
487 return 0;
488 }
489
490 /* get number of digits and add that */
491 r = (a->used - 1) * DIGIT_BIT;
492
493 /* take the last digit and count the bits in it */
494 q = a->dp[a->used - 1];
495 while (q > ((pstm_digit) 0)) {
496 ++r;
497 q >>= ((pstm_digit) 1);
498 }
499 return r;
500 }
501
502 /******************************************************************************/
pstm_unsigned_bin_size(pstm_int * a)503 int32 FAST_FUNC pstm_unsigned_bin_size(pstm_int *a)
504 {
505 int32 size = pstm_count_bits (a);
506 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
507 }
508
509 /******************************************************************************/
pstm_set(pstm_int * a,pstm_digit b)510 static void pstm_set(pstm_int *a, pstm_digit b)
511 {
512 pstm_zero(a);
513 a->dp[0] = b;
514 a->used = a->dp[0] ? 1 : 0;
515 }
516
517 /******************************************************************************/
518 /*
519 Right shift
520 */
pstm_rshd(pstm_int * a,int x)521 static void pstm_rshd(pstm_int *a, int x)
522 {
523 int y; //bbox: was int16
524
525 /* too many digits just zero and return */
526 if (x >= a->used) {
527 pstm_zero(a);
528 return;
529 }
530
531 /* shift */
532 for (y = 0; y < a->used - x; y++) {
533 a->dp[y] = a->dp[y+x];
534 }
535
536 /* zero rest */
537 for (; y < a->used; y++) {
538 a->dp[y] = 0;
539 }
540
541 /* decrement count */
542 a->used -= x;
543 pstm_clamp(a);
544 }
545
546 /******************************************************************************/
547 /*
548 Shift left a certain amount of digits.
549 */
550 #undef pstm_lshd
pstm_lshd(pstm_int * a,int b)551 static int32 pstm_lshd(pstm_int * a, int b)
552 {
553 int x; //bbox: was int16
554 int32 res;
555
556 /*
557 If its less than zero return.
558 */
559 if (b <= 0) {
560 return PSTM_OKAY;
561 }
562 /*
563 Grow to fit the new digits.
564 */
565 if (a->alloc < a->used + b) {
566 if ((res = pstm_grow (a, a->used + b)) != PSTM_OKAY) {
567 return res;
568 }
569 }
570
571 {
572 register pstm_digit *top, *bottom;
573 /*
574 Increment the used by the shift amount then copy upwards.
575 */
576 a->used += b;
577
578 /* top */
579 top = a->dp + a->used - 1;
580
581 /* base */
582 bottom = a->dp + a->used - 1 - b;
583 /*
584 This is implemented using a sliding window except the window goes the
585 other way around. Copying from the bottom to the top.
586 */
587 for (x = a->used - 1; x >= b; x--) {
588 *top-- = *bottom--;
589 }
590
591 /* zero the lower digits */
592 top = a->dp;
593 for (x = 0; x < b; x++) {
594 *top++ = 0;
595 }
596 }
597 return PSTM_OKAY;
598 }
599 #define pstm_lshd(a, b) (pstm_lshd(a, b), PSTM_OKAY)
600
601 /******************************************************************************/
602 /*
603 computes a = 2**b
604 */
pstm_2expt(pstm_int * a,int b)605 static int32 pstm_2expt(pstm_int *a, int b)
606 {
607 int z; //bbox: was int16
608
609 /* zero a as per default */
610 pstm_zero (a);
611
612 if (b < 0) {
613 return PSTM_OKAY;
614 }
615
616 z = b / DIGIT_BIT;
617 if (z >= PSTM_MAX_SIZE) {
618 return PS_LIMIT_FAIL;
619 }
620
621 /* set the used count of where the bit will go */
622 a->used = z + 1;
623
624 if (a->used > a->alloc) {
625 if (pstm_grow(a, a->used) != PSTM_OKAY) {
626 return PS_MEM_FAIL;
627 }
628 }
629
630 /* put the single bit in its place */
631 a->dp[z] = ((pstm_digit)1) << (b % DIGIT_BIT);
632 return PSTM_OKAY;
633 }
634
635 /******************************************************************************/
636 /*
637
638 */
pstm_mul_2(pstm_int * a,pstm_int * b)639 int32 FAST_FUNC pstm_mul_2(pstm_int * a, pstm_int * b)
640 {
641 int32 res;
642 int x, oldused; //bbox: was int16
643
644 /*
645 grow to accomodate result
646 */
647 if (b->alloc < a->used + 1) {
648 if ((res = pstm_grow (b, a->used + 1)) != PSTM_OKAY) {
649 return res;
650 }
651 }
652 oldused = b->used;
653 b->used = a->used;
654
655 {
656 register pstm_digit r, rr, *tmpa, *tmpb;
657
658 /* alias for source */
659 tmpa = a->dp;
660
661 /* alias for dest */
662 tmpb = b->dp;
663
664 /* carry */
665 r = 0;
666 for (x = 0; x < a->used; x++) {
667 /*
668 get what will be the *next* carry bit from the
669 MSB of the current digit
670 */
671 rr = *tmpa >> ((pstm_digit)(DIGIT_BIT - 1));
672 /*
673 now shift up this digit, add in the carry [from the previous]
674 */
675 *tmpb++ = ((*tmpa++ << ((pstm_digit)1)) | r);
676 /*
677 copy the carry that would be from the source
678 digit into the next iteration
679 */
680 r = rr;
681 }
682
683 /* new leading digit? */
684 if (r != 0 && b->used != (PSTM_MAX_SIZE-1)) {
685 /* add a MSB which is always 1 at this point */
686 *tmpb = 1;
687 ++(b->used);
688 }
689 /*
690 now zero any excess digits on the destination that we didn't write to
691 */
692 tmpb = b->dp + b->used;
693 for (x = b->used; x < oldused; x++) {
694 *tmpb++ = 0;
695 }
696 }
697 b->sign = a->sign;
698 return PSTM_OKAY;
699 }
700
701 /******************************************************************************/
702 /*
703 unsigned subtraction ||a|| >= ||b|| ALWAYS!
704 */
s_pstm_sub(pstm_int * a,pstm_int * b,pstm_int * c)705 int32 FAST_FUNC s_pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
706 {
707 int oldbused, oldused; //bbox: was int16
708 int32 x;
709 pstm_word t;
710
711 if (b->used > a->used) {
712 return PS_LIMIT_FAIL;
713 }
714 if (c->alloc < a->used) {
715 if ((x = pstm_grow (c, a->used)) != PSTM_OKAY) {
716 return x;
717 }
718 }
719 oldused = c->used;
720 oldbused = b->used;
721 c->used = a->used;
722 t = 0;
723
724 for (x = 0; x < oldbused; x++) {
725 t = ((pstm_word)a->dp[x]) - (((pstm_word)b->dp[x]) + t);
726 c->dp[x] = (pstm_digit)t;
727 t = (t >> DIGIT_BIT)&1;
728 }
729 for (; x < a->used; x++) {
730 t = ((pstm_word)a->dp[x]) - t;
731 c->dp[x] = (pstm_digit)t;
732 t = (t >> DIGIT_BIT);
733 }
734 for (; x < oldused; x++) {
735 c->dp[x] = 0;
736 }
737 pstm_clamp(c);
738 return PSTM_OKAY;
739 }
740
741 /******************************************************************************/
742 /*
743 unsigned addition
744 */
s_pstm_add(pstm_int * a,pstm_int * b,pstm_int * c)745 static int32 s_pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
746 {
747 int x, y, oldused; //bbox: was int16
748 register pstm_word t, adp, bdp;
749
750 y = a->used;
751 if (b->used > y) {
752 y = b->used;
753 }
754 oldused = c->used;
755 c->used = y;
756
757 if (c->used > c->alloc) {
758 if (pstm_grow(c, c->used) != PSTM_OKAY) {
759 return PS_MEM_FAIL;
760 }
761 }
762
763 t = 0;
764 for (x = 0; x < y; x++) {
765 if (a->used < x) {
766 adp = 0;
767 } else {
768 adp = (pstm_word)a->dp[x];
769 }
770 if (b->used < x) {
771 bdp = 0;
772 } else {
773 bdp = (pstm_word)b->dp[x];
774 }
775 t += (adp) + (bdp);
776 c->dp[x] = (pstm_digit)t;
777 t >>= DIGIT_BIT;
778 }
779 if (t != 0 && x < PSTM_MAX_SIZE) {
780 if (c->used == c->alloc) {
781 if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
782 return PS_MEM_FAIL;
783 }
784 }
785 c->dp[c->used++] = (pstm_digit)t;
786 ++x;
787 }
788
789 c->used = x;
790 for (; x < oldused; x++) {
791 c->dp[x] = 0;
792 }
793 pstm_clamp(c);
794 return PSTM_OKAY;
795 }
796
797
798 /******************************************************************************/
799 /*
800
801 */
pstm_sub(pstm_int * a,pstm_int * b,pstm_int * c)802 int32 FAST_FUNC pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
803 {
804 int32 res;
805 int sa, sb; //bbox: was int16
806
807 sa = a->sign;
808 sb = b->sign;
809
810 if (sa != sb) {
811 /*
812 subtract a negative from a positive, OR a positive from a negative.
813 For both, ADD their magnitudes, and use the sign of the first number.
814 */
815 c->sign = sa;
816 if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
817 return res;
818 }
819 } else {
820 /*
821 subtract a positive from a positive, OR a negative from a negative.
822 First, take the difference between their magnitudes, then...
823 */
824 if (pstm_cmp_mag (a, b) != PSTM_LT) {
825 /* Copy the sign from the first */
826 c->sign = sa;
827 /* The first has a larger or equal magnitude */
828 if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
829 return res;
830 }
831 } else {
832 /* The result has the _opposite_ sign from the first number. */
833 c->sign = (sa == PSTM_ZPOS) ? PSTM_NEG : PSTM_ZPOS;
834 /* The second has a larger magnitude */
835 if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
836 return res;
837 }
838 }
839 }
840 return PS_SUCCESS;
841 }
842
843 /******************************************************************************/
844 /*
845 c = a - b
846 */
847 #if 0 //UNUSED
848 int32 pstm_sub_d(psPool_t *pool, pstm_int *a, pstm_digit b, pstm_int *c)
849 {
850 pstm_int tmp;
851 int32 res;
852
853 if (pstm_init_size(pool, &tmp, sizeof(pstm_digit)) != PSTM_OKAY) {
854 return PS_MEM_FAIL;
855 }
856 pstm_set(&tmp, b);
857 res = pstm_sub(a, &tmp, c);
858 pstm_clear(&tmp);
859 return res;
860 }
861 #endif
862
863 /******************************************************************************/
864 /*
865 setups the montgomery reduction
866 */
pstm_montgomery_setup(pstm_int * a,pstm_digit * rho)867 static int32 pstm_montgomery_setup(pstm_int *a, pstm_digit *rho)
868 {
869 pstm_digit x, b;
870
871 /*
872 fast inversion mod 2**k
873 Based on the fact that
874 XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
875 => 2*X*A - X*X*A*A = 1
876 => 2*(1) - (1) = 1
877 */
878 b = a->dp[0];
879
880 if ((b & 1) == 0) {
881 psTraceCrypto("pstm_montogomery_setup failure\n");
882 return PS_ARG_FAIL;
883 }
884
885 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
886 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
887 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
888 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
889 #ifdef PSTM_64BIT
890 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
891 #endif
892 /* rho = -1/m mod b */
893 *rho = (pstm_digit)(((pstm_word) 1 << ((pstm_word) DIGIT_BIT)) -
894 ((pstm_word)x));
895 return PSTM_OKAY;
896 }
897
898 /******************************************************************************/
899 /*
900 * computes a = B**n mod b without division or multiplication useful for
901 * normalizing numbers in a Montgomery system.
902 */
pstm_montgomery_calc_normalization(pstm_int * a,pstm_int * b)903 static int32 pstm_montgomery_calc_normalization(pstm_int *a, pstm_int *b)
904 {
905 int32 x;
906 int bits; //bbox: was int16
907
908 /* how many bits of last digit does b use */
909 bits = pstm_count_bits (b) % DIGIT_BIT;
910 if (!bits) bits = DIGIT_BIT;
911
912 /* compute A = B^(n-1) * 2^(bits-1) */
913 if (b->used > 1) {
914 if ((x = pstm_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) !=
915 PSTM_OKAY) {
916 return x;
917 }
918 } else {
919 pstm_set(a, 1);
920 bits = 1;
921 }
922
923 /* now compute C = A * B mod b */
924 for (x = bits - 1; x < (int32)DIGIT_BIT; x++) {
925 if (pstm_mul_2 (a, a) != PSTM_OKAY) {
926 return PS_MEM_FAIL;
927 }
928 if (pstm_cmp_mag (a, b) != PSTM_LT) {
929 if (s_pstm_sub (a, b, a) != PSTM_OKAY) {
930 return PS_MEM_FAIL;
931 }
932 }
933 }
934 return PSTM_OKAY;
935 }
936
937 /******************************************************************************/
938 /*
939 c = a * 2**d
940 */
941 #undef pstm_mul_2d
pstm_mul_2d(pstm_int * a,int b,pstm_int * c)942 static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c)
943 {
944 pstm_digit carry, carrytmp, shift;
945 int x; //bbox: was int16
946
947 /* copy it */
948 if (pstm_copy(a, c) != PSTM_OKAY) {
949 return PS_MEM_FAIL;
950 }
951
952 /* handle whole digits */
953 if (b >= DIGIT_BIT) {
954 if (pstm_lshd(c, b/DIGIT_BIT) != PSTM_OKAY) {
955 return PS_MEM_FAIL;
956 }
957 }
958 b %= DIGIT_BIT;
959
960 /* shift the digits */
961 if (b != 0) {
962 carry = 0;
963 shift = DIGIT_BIT - b;
964 for (x = 0; x < c->used; x++) {
965 carrytmp = c->dp[x] >> shift;
966 c->dp[x] = (c->dp[x] << b) + carry;
967 carry = carrytmp;
968 }
969 /* store last carry if room */
970 if (carry && x < PSTM_MAX_SIZE) {
971 if (c->used == c->alloc) {
972 if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
973 return PS_MEM_FAIL;
974 }
975 }
976 c->dp[c->used++] = carry;
977 }
978 }
979 pstm_clamp(c);
980 return PSTM_OKAY;
981 }
982 #define pstm_mul_2d(a, b, c) (pstm_mul_2d(a, b, c), PSTM_OKAY)
983
984 /******************************************************************************/
985 /*
986 c = a mod 2**d
987 */
988 #undef pstm_mod_2d
pstm_mod_2d(pstm_int * a,int b,pstm_int * c)989 static int32 pstm_mod_2d(pstm_int *a, int b, pstm_int *c) //bbox: was int16 b
990 {
991 int x; //bbox: was int16
992
993 /* zero if count less than or equal to zero */
994 if (b <= 0) {
995 pstm_zero(c);
996 return PSTM_OKAY;
997 }
998
999 /* get copy of input */
1000 if (pstm_copy(a, c) != PSTM_OKAY) {
1001 return PS_MEM_FAIL;
1002 }
1003
1004 /* if 2**d is larger than we just return */
1005 if (b >= (DIGIT_BIT * a->used)) {
1006 return PSTM_OKAY;
1007 }
1008
1009 /* zero digits above the last digit of the modulus */
1010 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++)
1011 {
1012 c->dp[x] = 0;
1013 }
1014 /* clear the digit that is not completely outside/inside the modulus */
1015 c->dp[b / DIGIT_BIT] &= ~((pstm_digit)0) >> (DIGIT_BIT - b);
1016 pstm_clamp (c);
1017 return PSTM_OKAY;
1018 }
1019 #define pstm_mod_2d(a, b, c) (pstm_mod_2d(a, b, c), PSTM_OKAY)
1020
1021
1022 /******************************************************************************/
1023 /*
1024 c = a * b
1025 */
1026 #undef pstm_mul_d
pstm_mul_d(pstm_int * a,pstm_digit b,pstm_int * c)1027 static int32 pstm_mul_d(pstm_int *a, pstm_digit b, pstm_int *c)
1028 {
1029 pstm_word w;
1030 int32 res;
1031 int x, oldused; //bbox: was int16
1032
1033 if (c->alloc < a->used + 1) {
1034 if ((res = pstm_grow (c, a->used + 1)) != PSTM_OKAY) {
1035 return res;
1036 }
1037 }
1038 oldused = c->used;
1039 c->used = a->used;
1040 c->sign = a->sign;
1041 w = 0;
1042 for (x = 0; x < a->used; x++) {
1043 w = ((pstm_word)a->dp[x]) * ((pstm_word)b) + w;
1044 c->dp[x] = (pstm_digit)w;
1045 w = w >> DIGIT_BIT;
1046 }
1047 if (w != 0 && (a->used != PSTM_MAX_SIZE)) {
1048 c->dp[c->used++] = (pstm_digit)w;
1049 ++x;
1050 }
1051 for (; x < oldused; x++) {
1052 c->dp[x] = 0;
1053 }
1054 pstm_clamp(c);
1055 return PSTM_OKAY;
1056 }
1057 #define pstm_mul_d(a, b, c) (pstm_mul_d(a, b, c), PSTM_OKAY)
1058
1059 /******************************************************************************/
1060 /*
1061 c = a / 2**b
1062 */
1063 #undef pstm_div_2d
1064 #define pstm_div_2d(pool, a, b, c, d) \
1065 pstm_div_2d( a, b, c, d)
pstm_div_2d(psPool_t * pool,pstm_int * a,int b,pstm_int * c,pstm_int * d)1066 static int32 pstm_div_2d(psPool_t *pool, pstm_int *a, int b, pstm_int *c,
1067 pstm_int *d)
1068 {
1069 pstm_digit D, r, rr;
1070 int32 res;
1071 int x; //bbox: was int16
1072 pstm_int t;
1073
1074 /* if the shift count is <= 0 then we do no work */
1075 if (b <= 0) {
1076 if (pstm_copy (a, c) != PSTM_OKAY) {
1077 return PS_MEM_FAIL;
1078 }
1079 if (d != NULL) {
1080 pstm_zero (d);
1081 }
1082 return PSTM_OKAY;
1083 }
1084
1085 /* get the remainder */
1086 if (d != NULL) {
1087 if (pstm_init(pool, &t) != PSTM_OKAY) {
1088 return PS_MEM_FAIL;
1089 }
1090 if (pstm_mod_2d (a, b, &t) != PSTM_OKAY) {
1091 res = PS_MEM_FAIL;
1092 goto LBL_DONE;
1093 }
1094 }
1095
1096 /* copy */
1097 if (pstm_copy(a, c) != PSTM_OKAY) {
1098 res = PS_MEM_FAIL;
1099 goto LBL_DONE;
1100 }
1101
1102 /* shift by as many digits in the bit count */
1103 if (b >= (int32)DIGIT_BIT) {
1104 pstm_rshd (c, b / DIGIT_BIT);
1105 }
1106
1107 /* shift any bit count < DIGIT_BIT */
1108 D = (pstm_digit) (b % DIGIT_BIT);
1109 if (D != 0) {
1110 register pstm_digit *tmpc, mask, shift;
1111
1112 /* mask */
1113 mask = (((pstm_digit)1) << D) - 1;
1114
1115 /* shift for lsb */
1116 shift = DIGIT_BIT - D;
1117
1118 /* alias */
1119 tmpc = c->dp + (c->used - 1);
1120
1121 /* carry */
1122 r = 0;
1123 for (x = c->used - 1; x >= 0; x--) {
1124 /* get the lower bits of this word in a temp */
1125 rr = *tmpc & mask;
1126
1127 /* shift the current word and mix in the carry bits from previous */
1128 *tmpc = (*tmpc >> D) | (r << shift);
1129 --tmpc;
1130
1131 /* set the carry to the carry bits of the current word above */
1132 r = rr;
1133 }
1134 }
1135 pstm_clamp (c);
1136
1137 res = PSTM_OKAY;
1138 LBL_DONE:
1139 if (d != NULL) {
1140 if (pstm_copy(&t, d) != PSTM_OKAY) {
1141 res = PS_MEM_FAIL;
1142 }
1143 pstm_clear(&t);
1144 }
1145 return res;
1146 }
1147 #undef pstm_div_2d
1148 #define pstm_div_2d(pool, a, b, c, d) (pstm_div_2d(a, b, c, d), PSTM_OKAY)
1149
1150 /******************************************************************************/
1151 /*
1152 b = a/2
1153 */
1154 #if 0 //UNUSED
1155 int32 pstm_div_2(pstm_int * a, pstm_int * b)
1156 {
1157 int x, oldused; //bbox: was int16
1158
1159 if (b->alloc < a->used) {
1160 if (pstm_grow(b, a->used) != PSTM_OKAY) {
1161 return PS_MEM_FAIL;
1162 }
1163 }
1164 oldused = b->used;
1165 b->used = a->used;
1166 {
1167 register pstm_digit r, rr, *tmpa, *tmpb;
1168
1169 /* source alias */
1170 tmpa = a->dp + b->used - 1;
1171
1172 /* dest alias */
1173 tmpb = b->dp + b->used - 1;
1174
1175 /* carry */
1176 r = 0;
1177 for (x = b->used - 1; x >= 0; x--) {
1178 /* get the carry for the next iteration */
1179 rr = *tmpa & 1;
1180
1181 /* shift the current digit, add in carry and store */
1182 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1183
1184 /* forward carry to next iteration */
1185 r = rr;
1186 }
1187
1188 /* zero excess digits */
1189 tmpb = b->dp + b->used;
1190 for (x = b->used; x < oldused; x++) {
1191 *tmpb++ = 0;
1192 }
1193 }
1194 b->sign = a->sign;
1195 pstm_clamp (b);
1196 return PSTM_OKAY;
1197 }
1198 #endif
1199
1200 /******************************************************************************/
1201 /*
1202 Creates "a" then copies b into it
1203 */
1204 #undef pstm_init_copy
1205 #define pstm_init_copy(pool, a, b, toSqr) \
1206 pstm_init_copy( a, b, toSqr)
pstm_init_copy(psPool_t * pool,pstm_int * a,pstm_int * b,int toSqr)1207 static int32 pstm_init_copy(psPool_t *pool, pstm_int * a, pstm_int * b, int toSqr)
1208 {
1209 int x; //bbox: was int16
1210 int32 res;
1211
1212 if (a == b) {
1213 return PSTM_OKAY;
1214 }
1215 x = b->alloc;
1216
1217 if (toSqr) {
1218 /*
1219 Smart-size: Increasing size of a if b->used is roughly half
1220 of b->alloc because usage has shown that a lot of these copies
1221 go on to be squared and need these extra digits
1222 */
1223 if ((b->used * 2) + 2 >= x) {
1224 x = (b->used * 2) + 3;
1225 }
1226 }
1227 if ((res = pstm_init_size(pool, a, x)) != PSTM_OKAY) {
1228 return res;
1229 }
1230 return pstm_copy(b, a);
1231 }
1232 #undef pstm_init_copy
1233 #define pstm_init_copy(pool, a, b, toSqr) (pstm_init_copy(a, b, toSqr), PSTM_OKAY)
1234
1235 /******************************************************************************/
1236 /*
1237 With some compilers, we have seen issues linking with the builtin
1238 64 bit division routine. The issues with either manifest in a failure
1239 to find 'udivdi3' at link time, or a runtime invalid instruction fault
1240 during an RSA operation.
1241 The routine below divides a 64 bit unsigned int by a 32 bit unsigned int
1242 explicitly, rather than using the division operation
1243 The 64 bit result is placed in the 'numerator' parameter
1244 The 32 bit mod (remainder) of the division is the return parameter
1245 Based on implementations by:
1246 Copyright (C) 2003 Bernardo Innocenti <bernie@develer.com>
1247 Copyright (C) 1999 Hewlett-Packard Co
1248 Copyright (C) 1999 David Mosberger-Tang <davidm@hpl.hp.com>
1249 */
1250 #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
psDiv64(uint64 * numerator,uint32 denominator)1251 static uint32 psDiv64(uint64 *numerator, uint32 denominator)
1252 {
1253 uint64 rem = *numerator;
1254 uint64 b = denominator;
1255 uint64 res = 0;
1256 uint64 d = 1;
1257 uint32 high = rem >> 32;
1258
1259 if (high >= denominator) {
1260 high /= denominator;
1261 res = (uint64) high << 32;
1262 rem -= (uint64) (high * denominator) << 32;
1263 }
1264 while ((int64)b > 0 && b < rem) {
1265 b = b+b;
1266 d = d+d;
1267 }
1268 do {
1269 if (rem >= b) {
1270 rem -= b;
1271 res += d;
1272 }
1273 b >>= 1;
1274 d >>= 1;
1275 } while (d);
1276 *numerator = res;
1277 return rem;
1278 }
1279 #endif /* USE_MATRIX_DIV64 */
1280
1281 #if defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1282 typedef unsigned long uint128 __attribute__ ((mode(TI)));
psDiv128(uint128 * numerator,uint64 denominator)1283 static uint64 psDiv128(uint128 *numerator, uint64 denominator)
1284 {
1285 uint128 rem = *numerator;
1286 uint128 b = denominator;
1287 uint128 res = 0;
1288 uint128 d = 1;
1289 uint64 high = rem >> 64;
1290
1291 if (high >= denominator) {
1292 high /= denominator;
1293 res = (uint128) high << 64;
1294 rem -= (uint128) (high * denominator) << 64;
1295 }
1296 while ((uint128)b > 0 && b < rem) {
1297 b = b+b;
1298 d = d+d;
1299 }
1300 do {
1301 if (rem >= b) {
1302 rem -= b;
1303 res += d;
1304 }
1305 b >>= 1;
1306 d >>= 1;
1307 } while (d);
1308 *numerator = res;
1309 return rem;
1310 }
1311 #endif /* USE_MATRIX_DIV128 */
1312
1313 /******************************************************************************/
1314 /*
1315 a/b => cb + d == a
1316 */
pstm_div(psPool_t * pool,pstm_int * a,pstm_int * b,pstm_int * c,pstm_int * d)1317 static int32 pstm_div(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
1318 pstm_int *d)
1319 {
1320 pstm_int q, x, y, t1, t2;
1321 int32 res;
1322 int n, t, i, norm, neg; //bbox: was int16
1323
1324 /* is divisor zero ? */
1325 if (pstm_iszero (b) == 1) {
1326 return PS_LIMIT_FAIL;
1327 }
1328
1329 /* if a < b then q=0, r = a */
1330 if (pstm_cmp_mag (a, b) == PSTM_LT) {
1331 if (d != NULL) {
1332 if (pstm_copy(a, d) != PSTM_OKAY) {
1333 return PS_MEM_FAIL;
1334 }
1335 }
1336 if (c != NULL) {
1337 pstm_zero (c);
1338 }
1339 return PSTM_OKAY;
1340 }
1341 /*
1342 Smart-size inits
1343 */
1344 if ((res = pstm_init_size(pool, &t1, a->alloc)) != PSTM_OKAY) {
1345 return res;
1346 }
1347 if ((res = pstm_init_size(pool, &t2, 3)) != PSTM_OKAY) {
1348 goto LBL_T1;
1349 }
1350 if ((res = pstm_init_copy(pool, &x, a, 0)) != PSTM_OKAY) {
1351 goto LBL_T2;
1352 }
1353 /*
1354 Used to be an init_copy on b but pstm_grow was always hit with triple size
1355 */
1356 if ((res = pstm_init_size(pool, &y, b->used * 3)) != PSTM_OKAY) {
1357 goto LBL_X;
1358 }
1359 if ((res = pstm_copy(b, &y)) != PSTM_OKAY) {
1360 goto LBL_Y;
1361 }
1362
1363 /* fix the sign */
1364 neg = (a->sign == b->sign) ? PSTM_ZPOS : PSTM_NEG;
1365 x.sign = y.sign = PSTM_ZPOS;
1366
1367 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1368 norm = pstm_count_bits(&y) % DIGIT_BIT;
1369 if (norm < (int32)(DIGIT_BIT-1)) {
1370 norm = (DIGIT_BIT-1) - norm;
1371 if ((res = pstm_mul_2d(&x, norm, &x)) != PSTM_OKAY) {
1372 goto LBL_Y;
1373 }
1374 if ((res = pstm_mul_2d(&y, norm, &y)) != PSTM_OKAY) {
1375 goto LBL_Y;
1376 }
1377 } else {
1378 norm = 0;
1379 }
1380
1381 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1382 n = x.used - 1;
1383 t = y.used - 1;
1384
1385 if ((res = pstm_init_size(pool, &q, n - t + 1)) != PSTM_OKAY) {
1386 goto LBL_Y;
1387 }
1388 q.used = n - t + 1;
1389
1390 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1391 if ((res = pstm_lshd(&y, n - t)) != PSTM_OKAY) { /* y = y*b**{n-t} */
1392 goto LBL_Q;
1393 }
1394
1395 while (pstm_cmp (&x, &y) != PSTM_LT) {
1396 ++(q.dp[n - t]);
1397 if ((res = pstm_sub(&x, &y, &x)) != PSTM_OKAY) {
1398 goto LBL_Q;
1399 }
1400 }
1401
1402 /* reset y by shifting it back down */
1403 pstm_rshd (&y, n - t);
1404
1405 /* step 3. for i from n down to (t + 1) */
1406 for (i = n; i >= (t + 1); i--) {
1407 if (i > x.used) {
1408 continue;
1409 }
1410
1411 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1412 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1413 if (x.dp[i] == y.dp[t]) {
1414 q.dp[i - t - 1] = (pstm_digit)((((pstm_word)1) << DIGIT_BIT) - 1);
1415 } else {
1416 pstm_word tmp;
1417 tmp = ((pstm_word) x.dp[i]) << ((pstm_word) DIGIT_BIT);
1418 tmp |= ((pstm_word) x.dp[i - 1]);
1419 #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1420 psDiv64(&tmp, y.dp[t]);
1421 #elif defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1422 psDiv128(&tmp, y.dp[t]);
1423 #else
1424 tmp /= ((pstm_word) y.dp[t]);
1425 #endif /* USE_MATRIX_DIV64 */
1426 q.dp[i - t - 1] = (pstm_digit) (tmp);
1427 }
1428
1429 /* while (q{i-t-1} * (yt * b + y{t-1})) >
1430 xi * b**2 + xi-1 * b + xi-2
1431
1432 do q{i-t-1} -= 1;
1433 */
1434 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
1435 do {
1436 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
1437
1438 /* find left hand */
1439 pstm_zero (&t1);
1440 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1441 t1.dp[1] = y.dp[t];
1442 t1.used = 2;
1443 if ((res = pstm_mul_d (&t1, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
1444 goto LBL_Q;
1445 }
1446
1447 /* find right hand */
1448 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1449 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1450 t2.dp[2] = x.dp[i];
1451 t2.used = 3;
1452 } while (pstm_cmp_mag(&t1, &t2) == PSTM_GT);
1453
1454 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1455 if ((res = pstm_mul_d(&y, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
1456 goto LBL_Q;
1457 }
1458
1459 if ((res = pstm_lshd(&t1, i - t - 1)) != PSTM_OKAY) {
1460 goto LBL_Q;
1461 }
1462
1463 if ((res = pstm_sub(&x, &t1, &x)) != PSTM_OKAY) {
1464 goto LBL_Q;
1465 }
1466
1467 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1468 if (x.sign == PSTM_NEG) {
1469 if ((res = pstm_copy(&y, &t1)) != PSTM_OKAY) {
1470 goto LBL_Q;
1471 }
1472 if ((res = pstm_lshd (&t1, i - t - 1)) != PSTM_OKAY) {
1473 goto LBL_Q;
1474 }
1475 if ((res = pstm_add (&x, &t1, &x)) != PSTM_OKAY) {
1476 goto LBL_Q;
1477 }
1478 q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
1479 }
1480 }
1481 /*
1482 now q is the quotient and x is the remainder (which we have to normalize)
1483 */
1484 /* get sign before writing to c */
1485 x.sign = x.used == 0 ? PSTM_ZPOS : a->sign;
1486
1487 if (c != NULL) {
1488 pstm_clamp (&q);
1489 if (pstm_copy (&q, c) != PSTM_OKAY) {
1490 res = PS_MEM_FAIL;
1491 goto LBL_Q;
1492 }
1493 c->sign = neg;
1494 }
1495
1496 if (d != NULL) {
1497 if ((res = pstm_div_2d (pool, &x, norm, &x, NULL)) != PSTM_OKAY) {
1498 goto LBL_Q;
1499 }
1500 /*
1501 the following is a kludge, essentially we were seeing the right
1502 remainder but with excess digits that should have been zero
1503 */
1504 for (i = b->used; i < x.used; i++) {
1505 x.dp[i] = 0;
1506 }
1507 pstm_clamp(&x);
1508 if (pstm_copy (&x, d) != PSTM_OKAY) {
1509 res = PS_MEM_FAIL;
1510 goto LBL_Q;
1511 }
1512 }
1513
1514 res = PSTM_OKAY;
1515
1516 LBL_Q:pstm_clear (&q);
1517 LBL_Y:pstm_clear (&y);
1518 LBL_X:pstm_clear (&x);
1519 LBL_T2:pstm_clear (&t2);
1520 LBL_T1:pstm_clear (&t1);
1521
1522 return res;
1523 }
1524
1525 /******************************************************************************/
1526 /*
1527 Swap the elements of two integers, for cases where you can't simply swap
1528 the pstm_int pointers around
1529 */
pstm_exch(pstm_int * a,pstm_int * b)1530 static void pstm_exch(pstm_int * a, pstm_int * b)
1531 {
1532 pstm_int t;
1533
1534 t = *a;
1535 *a = *b;
1536 *b = t;
1537 }
1538
1539 /******************************************************************************/
1540 /*
1541 c = a mod b, 0 <= c < b
1542 */
pstm_mod(psPool_t * pool,pstm_int * a,pstm_int * b,pstm_int * c)1543 static int32 pstm_mod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
1544 {
1545 pstm_int t;
1546 int32 err;
1547 /*
1548 Smart-size
1549 */
1550 if ((err = pstm_init_size(pool, &t, b->alloc)) != PSTM_OKAY) {
1551 return err;
1552 }
1553 if ((err = pstm_div(pool, a, b, NULL, &t)) != PSTM_OKAY) {
1554 pstm_clear (&t);
1555 return err;
1556 }
1557 if (t.sign != b->sign) {
1558 err = pstm_add(&t, b, c);
1559 } else {
1560 pstm_exch (&t, c);
1561 }
1562 pstm_clear (&t);
1563 return err;
1564 }
1565
1566 /******************************************************************************/
1567 /*
1568 d = a * b (mod c)
1569 */
pstm_mulmod(psPool_t * pool,pstm_int * a,pstm_int * b,pstm_int * c,pstm_int * d)1570 int32 FAST_FUNC pstm_mulmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
1571 pstm_int *d)
1572 {
1573 int32 res;
1574 int size; //bbox: was int16
1575 pstm_int tmp;
1576
1577 /*
1578 Smart-size pstm_inits. d is an output that is influenced by this local 't'
1579 so don't shrink 'd' if it wants to becuase this will lead to an pstm_grow
1580 in RSA operations
1581 */
1582 size = a->used + b->used + 1;
1583 if ((a == d) && (size < a->alloc)) {
1584 size = a->alloc;
1585 }
1586 if ((res = pstm_init_size(pool, &tmp, size)) != PSTM_OKAY) {
1587 return res;
1588 }
1589 if ((res = pstm_mul_comba(pool, a, b, &tmp, NULL, 0)) != PSTM_OKAY) {
1590 pstm_clear(&tmp);
1591 return res;
1592 }
1593 res = pstm_mod(pool, &tmp, c, d);
1594 pstm_clear(&tmp);
1595 return res;
1596 }
1597
1598 /******************************************************************************/
1599 /*
1600 * y = g**x (mod b)
1601 * Some restrictions... x must be positive and < b
1602 */
pstm_exptmod(psPool_t * pool,pstm_int * G,pstm_int * X,pstm_int * P,pstm_int * Y)1603 int32 FAST_FUNC pstm_exptmod(psPool_t *pool, pstm_int *G, pstm_int *X, pstm_int *P,
1604 pstm_int *Y)
1605 {
1606 pstm_int M[32], res; /* Keep this winsize based: (1 << max_winsize) */
1607 pstm_digit buf, mp;
1608 pstm_digit *paD;
1609 int32 err, bitbuf;
1610 int bitcpy, bitcnt, mode, digidx, x, y, winsize; //bbox: was int16
1611 uint32 paDlen;
1612
1613 /* set window size from what user set as optimization */
1614 x = pstm_count_bits(X);
1615 if (x < 50) {
1616 winsize = 2;
1617 } else {
1618 winsize = PS_EXPTMOD_WINSIZE;
1619 }
1620
1621 /* now setup montgomery */
1622 if ((err = pstm_montgomery_setup (P, &mp)) != PSTM_OKAY) {
1623 return err;
1624 }
1625
1626 /* setup result */
1627 if ((err = pstm_init_size(pool, &res, (P->used * 2) + 1)) != PSTM_OKAY) {
1628 return err;
1629 }
1630 /*
1631 create M table
1632 The M table contains powers of the input base, e.g. M[x] = G^x mod P
1633 The first half of the table is not computed though except for M[0] and M[1]
1634 */
1635 /* now we need R mod m */
1636 if ((err = pstm_montgomery_calc_normalization (&res, P)) != PSTM_OKAY) {
1637 goto LBL_RES;
1638 }
1639 /*
1640 init M array
1641 init first cell
1642 */
1643 if ((err = pstm_init_size(pool, &M[1], res.used)) != PSTM_OKAY) {
1644 goto LBL_RES;
1645 }
1646
1647 /* now set M[1] to G * R mod m */
1648 if (pstm_cmp_mag(P, G) != PSTM_GT) {
1649 /* G > P so we reduce it first */
1650 if ((err = pstm_mod(pool, G, P, &M[1])) != PSTM_OKAY) {
1651 goto LBL_M;
1652 }
1653 } else {
1654 if ((err = pstm_copy(G, &M[1])) != PSTM_OKAY) {
1655 goto LBL_M;
1656 }
1657 }
1658 if ((err = pstm_mulmod (pool, &M[1], &res, P, &M[1])) != PSTM_OKAY) {
1659 goto LBL_M;
1660 }
1661 /*
1662 Pre-allocated digit. Used for mul, sqr, AND reduce
1663 */
1664 paDlen = ((M[1].used + 3) * 2) * sizeof(pstm_digit);
1665 paD = xzalloc(paDlen);//bbox
1666 /*
1667 compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times
1668 */
1669 if (pstm_init_copy(pool, &M[1 << (winsize - 1)], &M[1], 1) != PSTM_OKAY) {
1670 err = PS_MEM_FAIL;
1671 goto LBL_PAD;
1672 }
1673 for (x = 0; x < (winsize - 1); x++) {
1674 if ((err = pstm_sqr_comba (pool, &M[1 << (winsize - 1)],
1675 &M[1 << (winsize - 1)], paD, paDlen)) != PSTM_OKAY) {
1676 goto LBL_PAD;
1677 }
1678 if ((err = pstm_montgomery_reduce(pool, &M[1 << (winsize - 1)], P, mp,
1679 paD, paDlen)) != PSTM_OKAY) {
1680 goto LBL_PAD;
1681 }
1682 }
1683 /*
1684 now init the second half of the array
1685 */
1686 for (x = (1<<(winsize-1)) + 1; x < (1 << winsize); x++) {
1687 if ((err = pstm_init_size(pool, &M[x], M[1<<(winsize-1)].alloc + 1))
1688 != PSTM_OKAY) {
1689 for (y = 1<<(winsize-1); y < x; y++) {
1690 pstm_clear(&M[y]);
1691 }
1692 goto LBL_PAD;
1693 }
1694 }
1695
1696 /* create upper table */
1697 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
1698 if ((err = pstm_mul_comba(pool, &M[x - 1], &M[1], &M[x], paD, paDlen))
1699 != PSTM_OKAY) {
1700 goto LBL_MARRAY;
1701 }
1702 if ((err = pstm_montgomery_reduce(pool, &M[x], P, mp, paD, paDlen)) !=
1703 PSTM_OKAY) {
1704 goto LBL_MARRAY;
1705 }
1706 }
1707
1708 /* set initial mode and bit cnt */
1709 mode = 0;
1710 bitcnt = 1;
1711 buf = 0;
1712 digidx = X->used - 1;
1713 bitcpy = 0;
1714 bitbuf = 0;
1715
1716 for (;;) {
1717 /* grab next digit as required */
1718 if (--bitcnt == 0) {
1719 /* if digidx == -1 we are out of digits so break */
1720 if (digidx == -1) {
1721 break;
1722 }
1723 /* read next digit and reset bitcnt */
1724 buf = X->dp[digidx--];
1725 bitcnt = (int32)DIGIT_BIT;
1726 }
1727
1728 /* grab the next msb from the exponent */
1729 y = (pstm_digit)(buf >> (DIGIT_BIT - 1)) & 1;
1730 buf <<= (pstm_digit)1;
1731 /*
1732 If the bit is zero and mode == 0 then we ignore it.
1733 These represent the leading zero bits before the first 1 bit
1734 in the exponent. Technically this opt is not required but it
1735 does lower the # of trivial squaring/reductions used
1736 */
1737 if (mode == 0 && y == 0) {
1738 continue;
1739 }
1740
1741 /* if the bit is zero and mode == 1 then we square */
1742 if (mode == 1 && y == 0) {
1743 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1744 PSTM_OKAY) {
1745 goto LBL_MARRAY;
1746 }
1747 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1748 != PSTM_OKAY) {
1749 goto LBL_MARRAY;
1750 }
1751 continue;
1752 }
1753
1754 /* else we add it to the window */
1755 bitbuf |= (y << (winsize - ++bitcpy));
1756 mode = 2;
1757
1758 if (bitcpy == winsize) {
1759 /* ok window is filled so square as required and mul square first */
1760 for (x = 0; x < winsize; x++) {
1761 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1762 PSTM_OKAY) {
1763 goto LBL_MARRAY;
1764 }
1765 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
1766 paDlen)) != PSTM_OKAY) {
1767 goto LBL_MARRAY;
1768 }
1769 }
1770
1771 /* then multiply */
1772 if ((err = pstm_mul_comba(pool, &res, &M[bitbuf], &res, paD,
1773 paDlen)) != PSTM_OKAY) {
1774 goto LBL_MARRAY;
1775 }
1776 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1777 != PSTM_OKAY) {
1778 goto LBL_MARRAY;
1779 }
1780
1781 /* empty window and reset */
1782 bitcpy = 0;
1783 bitbuf = 0;
1784 mode = 1;
1785 }
1786 }
1787
1788 /* if bits remain then square/multiply */
1789 if (mode == 2 && bitcpy > 0) {
1790 /* square then multiply if the bit is set */
1791 for (x = 0; x < bitcpy; x++) {
1792 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1793 PSTM_OKAY) {
1794 goto LBL_MARRAY;
1795 }
1796 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1797 != PSTM_OKAY) {
1798 goto LBL_MARRAY;
1799 }
1800
1801 /* get next bit of the window */
1802 bitbuf <<= 1;
1803 if ((bitbuf & (1 << winsize)) != 0) {
1804 /* then multiply */
1805 if ((err = pstm_mul_comba(pool, &res, &M[1], &res, paD, paDlen))
1806 != PSTM_OKAY) {
1807 goto LBL_MARRAY;
1808 }
1809 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
1810 paDlen)) != PSTM_OKAY) {
1811 goto LBL_MARRAY;
1812 }
1813 }
1814 }
1815 }
1816 /*
1817 Fix up result if Montgomery reduction is used recall that any value in a
1818 Montgomery system is actually multiplied by R mod n. So we have to reduce
1819 one more time to cancel out the factor of R.
1820 */
1821 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen)) !=
1822 PSTM_OKAY) {
1823 goto LBL_MARRAY;
1824 }
1825 /* swap res with Y */
1826 if ((err = pstm_copy (&res, Y)) != PSTM_OKAY) {
1827 goto LBL_MARRAY;
1828 }
1829 err = PSTM_OKAY;
1830 LBL_MARRAY:
1831 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1832 pstm_clear(&M[x]);
1833 }
1834 LBL_PAD:psFree(paD, pool);
1835 LBL_M: pstm_clear(&M[1]);
1836 LBL_RES:pstm_clear(&res);
1837 return err;
1838 }
1839
1840 /******************************************************************************/
1841 /*
1842
1843 */
pstm_add(pstm_int * a,pstm_int * b,pstm_int * c)1844 int32 FAST_FUNC pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
1845 {
1846 int32 res;
1847 int sa, sb; //bbox: was int16
1848
1849 /* get sign of both inputs */
1850 sa = a->sign;
1851 sb = b->sign;
1852
1853 /* handle two cases, not four */
1854 if (sa == sb) {
1855 /* both positive or both negative, add their mags, copy the sign */
1856 c->sign = sa;
1857 if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
1858 return res;
1859 }
1860 } else {
1861 /*
1862 one positive, the other negative
1863 subtract the one with the greater magnitude from the one of the lesser
1864 magnitude. The result gets the sign of the one with the greater mag.
1865 */
1866 if (pstm_cmp_mag (a, b) == PSTM_LT) {
1867 c->sign = sb;
1868 if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
1869 return res;
1870 }
1871 } else {
1872 c->sign = sa;
1873 if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
1874 return res;
1875 }
1876 }
1877 }
1878 return PS_SUCCESS;
1879 }
1880
1881 /******************************************************************************/
1882 /*
1883 reverse an array, used for radix code
1884 */
pstm_reverse(unsigned char * s,int len)1885 static void pstm_reverse (unsigned char *s, int len) //bbox: was int16 len
1886 {
1887 int32 ix, iy;
1888 unsigned char t;
1889
1890 ix = 0;
1891 iy = len - 1;
1892 while (ix < iy) {
1893 t = s[ix];
1894 s[ix] = s[iy];
1895 s[iy] = t;
1896 ++ix;
1897 --iy;
1898 }
1899 }
1900 /******************************************************************************/
1901 /*
1902 No reverse. Useful in some of the EIP-154 PKA stuff where special byte
1903 order seems to come into play more often
1904 */
1905 #if 0 //UNUSED
1906 int32 pstm_to_unsigned_bin_nr(psPool_t *pool, pstm_int *a, unsigned char *b)
1907 {
1908 int32 res;
1909 int x; //bbox: was int16
1910 pstm_int t = { 0 };
1911
1912 if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
1913 return res;
1914 }
1915
1916 x = 0;
1917 while (pstm_iszero (&t) == 0) {
1918 b[x++] = (unsigned char) (t.dp[0] & 255);
1919 if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
1920 pstm_clear(&t);
1921 return res;
1922 }
1923 }
1924 pstm_clear(&t);
1925 return PS_SUCCESS;
1926 }
1927 #endif
1928 /******************************************************************************/
1929 /*
1930
1931 */
pstm_to_unsigned_bin(psPool_t * pool,pstm_int * a,unsigned char * b)1932 int32 FAST_FUNC pstm_to_unsigned_bin(psPool_t *pool, pstm_int *a, unsigned char *b)
1933 {
1934 int32 res;
1935 int x; //bbox: was int16
1936 pstm_int t = { 0 };
1937
1938 if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
1939 return res;
1940 }
1941
1942 x = 0;
1943 while (pstm_iszero (&t) == 0) {
1944 b[x++] = (unsigned char) (t.dp[0] & 255);
1945 if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
1946 pstm_clear(&t);
1947 return res;
1948 }
1949 }
1950 pstm_reverse (b, x);
1951 pstm_clear(&t);
1952 return PS_SUCCESS;
1953 }
1954
1955 #if 0 //UNUSED
1956 /******************************************************************************/
1957 /*
1958 compare against a single digit
1959 */
1960 static int32 pstm_cmp_d(pstm_int *a, pstm_digit b)
1961 {
1962 /* compare based on sign */
1963 if ((b && a->used == 0) || a->sign == PSTM_NEG) {
1964 return PSTM_LT;
1965 }
1966
1967 /* compare based on magnitude */
1968 if (a->used > 1) {
1969 return PSTM_GT;
1970 }
1971
1972 /* compare the only digit of a to b */
1973 if (a->dp[0] > b) {
1974 return PSTM_GT;
1975 } else if (a->dp[0] < b) {
1976 return PSTM_LT;
1977 } else {
1978 return PSTM_EQ;
1979 }
1980 }
1981
1982 /*
1983 Need invmod for ECC and also private key loading for hardware crypto
1984 in cases where dQ > dP. The values must be switched and a new qP must be
1985 calculated using this function
1986 */
1987 //bbox: pool unused
1988 #define pstm_invmod_slow(pool, a, b, c) \
1989 pstm_invmod_slow( a, b, c)
1990 static int32 pstm_invmod_slow(psPool_t *pool, pstm_int * a, pstm_int * b,
1991 pstm_int * c)
1992 {
1993 pstm_int x, y, u, v, A, B, C, D;
1994 int32 res;
1995
1996 /* b cannot be negative */
1997 if (b->sign == PSTM_NEG || pstm_iszero(b) == 1) {
1998 return PS_LIMIT_FAIL;
1999 }
2000
2001 /* init temps */
2002 if (pstm_init_size(pool, &x, b->used) != PSTM_OKAY) {
2003 return PS_MEM_FAIL;
2004 }
2005
2006 /* x = a, y = b */
2007 if ((res = pstm_mod(pool, a, b, &x)) != PSTM_OKAY) {
2008 goto LBL_X;
2009 }
2010
2011 if (pstm_init_copy(pool, &y, b, 0) != PSTM_OKAY) {
2012 goto LBL_X;
2013 }
2014
2015 /* 2. [modified] if x,y are both even then return an error! */
2016 if (pstm_iseven (&x) == 1 && pstm_iseven (&y) == 1) {
2017 res = PS_FAILURE;
2018 goto LBL_Y;
2019 }
2020
2021 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2022 if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
2023 goto LBL_Y;
2024 }
2025 if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
2026 goto LBL_U;
2027 }
2028
2029 if ((res = pstm_init_size(pool, &A, sizeof(pstm_digit))) != PSTM_OKAY) {
2030 goto LBL_V;
2031 }
2032
2033 if ((res = pstm_init_size(pool, &D, sizeof(pstm_digit))) != PSTM_OKAY) {
2034 goto LBL_A;
2035 }
2036 pstm_set (&A, 1);
2037 pstm_set (&D, 1);
2038
2039 if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
2040 goto LBL_D;
2041 }
2042 if ((res = pstm_init(pool, &C)) != PSTM_OKAY) {
2043 goto LBL_B;
2044 }
2045
2046 top:
2047 /* 4. while u is even do */
2048 while (pstm_iseven (&u) == 1) {
2049 /* 4.1 u = u/2 */
2050 if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
2051 goto LBL_C;
2052 }
2053
2054 /* 4.2 if A or B is odd then */
2055 if (pstm_isodd (&A) == 1 || pstm_isodd (&B) == 1) {
2056 /* A = (A+y)/2, B = (B-x)/2 */
2057 if ((res = pstm_add (&A, &y, &A)) != PSTM_OKAY) {
2058 goto LBL_C;
2059 }
2060 if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
2061 goto LBL_C;
2062 }
2063 }
2064 /* A = A/2, B = B/2 */
2065 if ((res = pstm_div_2 (&A, &A)) != PSTM_OKAY) {
2066 goto LBL_C;
2067 }
2068 if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
2069 goto LBL_C;
2070 }
2071 }
2072
2073 /* 5. while v is even do */
2074 while (pstm_iseven (&v) == 1) {
2075 /* 5.1 v = v/2 */
2076 if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
2077 goto LBL_C;
2078 }
2079
2080 /* 5.2 if C or D is odd then */
2081 if (pstm_isodd (&C) == 1 || pstm_isodd (&D) == 1) {
2082 /* C = (C+y)/2, D = (D-x)/2 */
2083 if ((res = pstm_add (&C, &y, &C)) != PSTM_OKAY) {
2084 goto LBL_C;
2085 }
2086 if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
2087 goto LBL_C;
2088 }
2089 }
2090 /* C = C/2, D = D/2 */
2091 if ((res = pstm_div_2 (&C, &C)) != PSTM_OKAY) {
2092 goto LBL_C;
2093 }
2094 if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
2095 goto LBL_C;
2096 }
2097 }
2098
2099 /* 6. if u >= v then */
2100 if (pstm_cmp (&u, &v) != PSTM_LT) {
2101 /* u = u - v, A = A - C, B = B - D */
2102 if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
2103 goto LBL_C;
2104 }
2105 if ((res = pstm_sub (&A, &C, &A)) != PSTM_OKAY) {
2106 goto LBL_C;
2107 }
2108 if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
2109 goto LBL_C;
2110 }
2111 } else {
2112 /* v - v - u, C = C - A, D = D - B */
2113 if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
2114 goto LBL_C;
2115 }
2116 if ((res = pstm_sub (&C, &A, &C)) != PSTM_OKAY) {
2117 goto LBL_C;
2118 }
2119 if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
2120 goto LBL_C;
2121 }
2122 }
2123
2124 /* if not zero goto step 4 */
2125 if (pstm_iszero (&u) == 0)
2126 goto top;
2127
2128 /* now a = C, b = D, gcd == g*v */
2129
2130 /* if v != 1 then there is no inverse */
2131 if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
2132 res = PS_FAILURE;
2133 goto LBL_C;
2134 }
2135
2136 /* if its too low */
2137 while (pstm_cmp_d(&C, 0) == PSTM_LT) {
2138 if ((res = pstm_add(&C, b, &C)) != PSTM_OKAY) {
2139 goto LBL_C;
2140 }
2141 }
2142
2143 /* too big */
2144 while (pstm_cmp_mag(&C, b) != PSTM_LT) {
2145 if ((res = pstm_sub(&C, b, &C)) != PSTM_OKAY) {
2146 goto LBL_C;
2147 }
2148 }
2149
2150 /* C is now the inverse */
2151 if ((res = pstm_copy(&C, c)) != PSTM_OKAY) {
2152 goto LBL_C;
2153 }
2154 res = PSTM_OKAY;
2155
2156 LBL_C: pstm_clear(&C);
2157 LBL_D: pstm_clear(&D);
2158 LBL_B: pstm_clear(&B);
2159 LBL_A: pstm_clear(&A);
2160 LBL_V: pstm_clear(&v);
2161 LBL_U: pstm_clear(&u);
2162 LBL_Y: pstm_clear(&y);
2163 LBL_X: pstm_clear(&x);
2164
2165 return res;
2166 }
2167
2168 /* c = 1/a (mod b) for odd b only */
2169 int32 pstm_invmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
2170 {
2171 pstm_int x, y, u, v, B, D;
2172 int32 res;
2173 int neg, sanity; //bbox: was uint16
2174
2175 /* 2. [modified] b must be odd */
2176 if (pstm_iseven (b) == 1) {
2177 return pstm_invmod_slow(pool, a,b,c);
2178 }
2179
2180 /* x == modulus, y == value to invert */
2181 if ((res = pstm_init_copy(pool, &x, b, 0)) != PSTM_OKAY) {
2182 return res;
2183 }
2184
2185 if ((res = pstm_init_size(pool, &y, a->alloc)) != PSTM_OKAY) {
2186 goto LBL_X;
2187 }
2188
2189 /* we need y = |a| */
2190 pstm_abs(a, &y);
2191
2192 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2193 if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
2194 goto LBL_Y;
2195 }
2196 if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
2197 goto LBL_U;
2198 }
2199 if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
2200 goto LBL_V;
2201 }
2202 if ((res = pstm_init(pool, &D)) != PSTM_OKAY) {
2203 goto LBL_B;
2204 }
2205
2206 pstm_set (&D, 1);
2207
2208 sanity = 0;
2209 top:
2210 /* 4. while u is even do */
2211 while (pstm_iseven (&u) == 1) {
2212 /* 4.1 u = u/2 */
2213 if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
2214 goto LBL_D;
2215 }
2216
2217 /* 4.2 if B is odd then */
2218 if (pstm_isodd (&B) == 1) {
2219 if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
2220 goto LBL_D;
2221 }
2222 }
2223 /* B = B/2 */
2224 if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
2225 goto LBL_D;
2226 }
2227 }
2228
2229 /* 5. while v is even do */
2230 while (pstm_iseven (&v) == 1) {
2231 /* 5.1 v = v/2 */
2232 if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
2233 goto LBL_D;
2234 }
2235 /* 5.2 if D is odd then */
2236 if (pstm_isodd (&D) == 1) {
2237 /* D = (D-x)/2 */
2238 if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
2239 goto LBL_D;
2240 }
2241 }
2242 /* D = D/2 */
2243 if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
2244 goto LBL_D;
2245 }
2246 }
2247
2248 /* 6. if u >= v then */
2249 if (pstm_cmp (&u, &v) != PSTM_LT) {
2250 /* u = u - v, B = B - D */
2251 if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
2252 goto LBL_D;
2253 }
2254 if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
2255 goto LBL_D;
2256 }
2257 } else {
2258 /* v - v - u, D = D - B */
2259 if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
2260 goto LBL_D;
2261 }
2262 if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
2263 goto LBL_D;
2264 }
2265 }
2266
2267 /* if not zero goto step 4 */
2268 if (sanity++ > 1000) {
2269 res = PS_LIMIT_FAIL;
2270 goto LBL_D;
2271 }
2272 if (pstm_iszero (&u) == 0) {
2273 goto top;
2274 }
2275
2276 /* now a = C, b = D, gcd == g*v */
2277
2278 /* if v != 1 then there is no inverse */
2279 if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
2280 res = PS_FAILURE;
2281 goto LBL_D;
2282 }
2283
2284 /* b is now the inverse */
2285 neg = a->sign;
2286 while (D.sign == PSTM_NEG) {
2287 if ((res = pstm_add (&D, b, &D)) != PSTM_OKAY) {
2288 goto LBL_D;
2289 }
2290 }
2291 if ((res = pstm_copy (&D, c)) != PSTM_OKAY) {
2292 goto LBL_D;
2293 }
2294 c->sign = neg;
2295 res = PSTM_OKAY;
2296
2297 LBL_D: pstm_clear(&D);
2298 LBL_B: pstm_clear(&B);
2299 LBL_V: pstm_clear(&v);
2300 LBL_U: pstm_clear(&u);
2301 LBL_Y: pstm_clear(&y);
2302 LBL_X: pstm_clear(&x);
2303 return res;
2304 }
2305 #endif //UNUSED
2306
2307 #endif /* !DISABLE_PSTM */
2308 /******************************************************************************/
2309