File Coverage

inc/matrixssl-3-9-3-open/crypto/math/pstm.c
Criterion Covered Total %
statement 517 886 58.3
branch 275 678 40.5
condition n/a
subroutine n/a
pod n/a
total 792 1564 50.6


line stmt bran cond sub pod time code
1             /**
2             * @file pstm.c
3             * @version 950bba4 (HEAD -> master)
4             *
5             * Multiprecision number implementation.
6             */
7             /*
8             * Copyright (c) 2013-2017 INSIDE Secure Corporation
9             * Copyright (c) PeerSec Networks, 2002-2011
10             * All Rights Reserved
11             *
12             * The latest version of this code is available at http://www.matrixssl.org
13             *
14             * This software is open source; you can redistribute it and/or modify
15             * it under the terms of the GNU General Public License as published by
16             * the Free Software Foundation; either version 2 of the License, or
17             * (at your option) any later version.
18             *
19             * This General Public License does NOT permit incorporating this software
20             * into proprietary programs. If you are unable to comply with the GPL, a
21             * commercial license for this software may be purchased from INSIDE at
22             * http://www.insidesecure.com/
23             *
24             * This program is distributed in WITHOUT ANY WARRANTY; without even the
25             * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
26             * See the GNU General Public License for more details.
27             *
28             * You should have received a copy of the GNU General Public License
29             * along with this program; if not, write to the Free Software
30             * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
31             * http://www.gnu.org/copyleft/gpl.html
32             */
33             /******************************************************************************/
34              
35             /* This pstm mathematics library is
36             based on libraries by Tom St Denis. */
37              
38             #include "../cryptoImpl.h"
39             #include "pstmnt.h"
40              
41             #include /* toupper() */
42              
43             #if defined(USE_MATRIX_RSA) || defined(USE_MATRIX_ECC) || defined(USE_MATRIX_DH) || defined(USE_CL_RSA) || defined(USE_CL_DH) || defined(USE_QUICK_ASSIST_RSA) || defined(USE_QUICK_ASSIST_ECC)
44              
45             static int32_t pstm_mul_2d(const pstm_int *a, int16_t b, pstm_int *c);
46              
47             /******************************************************************************/
48             /**
49             Initialize a pstm_int and allocate working memory for a given initial size.
50              
51             @param[in] pool Memory pool to use for allocation.
52             @param[in,out] a Allocated pstm_int to initialize.
53             @param[in] size Number of digits to pre-allocate for integer. Typically
54             a digit is 32 or 64 bits.
55              
56             @return < 0 on failure, >=0 on success.
57             */
58 4545796           int32_t pstm_init_size(psPool_t *pool, pstm_int *a, psSize_t size)
59             {
60             uint16_t x;
61              
62 4545796 50         if (size > PSTM_MAX_SIZE)
63             {
64 0           return PSTM_MEM;
65             }
66 4545796           a->dp = psMalloc(pool, sizeof(pstm_digit) * size);
67 4545796 50         if (a->dp == NULL)
68             {
69 0           return PSTM_MEM;
70             }
71 4545796           a->pool = pool; /* Pool to use when growing or shrinking digits */
72 4545796           a->used = 0; /* Zero of the digits are currently used */
73 4545796           a->alloc = size; /* How many digits are pre-allocated */
74 4545796           a->sign = PSTM_ZPOS; /* Number is positive */
75             /* zero the digits */
76 83222287 100         for (x = 0; x < size; x++)
77             {
78 78676491           a->dp[x] = 0;
79             }
80 4545796           return PSTM_OKAY;
81             }
82              
83             /******************************************************************************/
84             /*
85             Init a new pstm_int with a default size.
86             */
87 74500           int32_t pstm_init(psPool_t *pool, pstm_int *a)
88             {
89 74500           return pstm_init_size(pool, a, (MIN_RSA_BITS / DIGIT_BIT) * 3);
90             }
91              
92             /******************************************************************************/
93             /**
94             Grow a pstm_int to the give size in digits.
95              
96             @param[in,out] a Allocated and initialized pstm_int to grow.
97             @param[in] size Number of digits to grow to. This is not the the number
98             of digits to grow _by_. 'size' <= current size is ignored, so this
99             api cannot be used to shrink the integer.
100             */
101 27           int32_t pstm_grow(pstm_int *a, psSize_t size)
102             {
103             uint16_t i;
104             pstm_digit *tmp;
105              
106 27 50         if (size > PSTM_MAX_SIZE)
107             {
108 0           return PSTM_MEM;
109             }
110             /* If the alloc size is smaller alloc more ram. */
111 27 50         if (a->alloc < size)
112             {
113             /*
114             Reallocate the array a->dp
115             We store the return in a temporary variable in case the operation
116             failed we don't want to overwrite the dp member of a.
117             */
118 27           tmp = psRealloc(a->dp, sizeof(pstm_digit) * size, a->pool);
119 27 50         if (tmp == NULL)
120             {
121             /* reallocation failed but "a" is still valid [can be freed] */
122 0           return PSTM_MEM;
123             }
124             /* reallocation succeeded so set a->dp */
125 27           a->dp = tmp;
126 27           i = a->alloc;
127 27           a->alloc = size;
128             /* zero excess digits */
129 54 100         for (; i < a->alloc; i++)
130             {
131 27           a->dp[i] = 0;
132             }
133             }
134 27           return PSTM_OKAY;
135             }
136              
137             /******************************************************************************/
138             /*
139             copy, b = a (b must be pre-allocated)
140             */
141 22355403           int32_t pstm_copy(const pstm_int *a, pstm_int *b)
142             {
143             int32_t res, n;
144              
145             /* If dst == src do nothing */
146 22355403 100         if (a == b)
147             {
148 20801529           return PSTM_OKAY;
149             }
150             /* Grow dest */
151 1553874 50         if (b->alloc < a->used)
152             {
153 0 0         if ((res = pstm_grow(b, a->used)) != PSTM_OKAY)
154             {
155 0           return res;
156             }
157             }
158             /* Zero b and copy the parameters over */
159             {
160             /* pointer aliases */
161             register pstm_digit *tmpa, *tmpb;
162              
163             /* source */
164 1553874           tmpa = a->dp;
165             /* destination */
166 1553874           tmpb = b->dp;
167              
168             /* copy all the digits */
169 15765710 100         for (n = 0; n < a->used; n++)
170             {
171 14211836           *tmpb++ = *tmpa++;
172             }
173             /* clear high digits */
174 1553981 100         for (; n < b->used; n++)
175             {
176 107           *tmpb++ = 0;
177             }
178             }
179             /* copy used count and sign */
180 1553874           b->used = a->used;
181 1553874           b->sign = a->sign;
182 1553874           return PSTM_OKAY;
183             }
184              
185             /******************************************************************************/
186             /**
187             b = |a|.
188             Copy 'a' to 'b' and make positive.
189             */
190 35032           int32_t pstm_abs(const pstm_int *a, pstm_int *b)
191             {
192 35032 50         if (pstm_copy(a, b) != PSTM_OKAY)
193             {
194 0           return PSTM_MEM;
195             }
196 35032           b->sign = 0;
197 35032           return PSTM_OKAY;
198             }
199              
200             /******************************************************************************/
201             /**
202             Trim unused digits.
203              
204             This is used to ensure that leading zero digits are trimed and the
205             leading "used" digit will be non-zero. Typically very fast. Also fixes
206             the sign if there are no more leading digits
207             */
208 84916264           void pstm_clamp(pstm_int *a)
209             {
210             /* decrease used while the most significant digit is zero. */
211 114742355 100         while (a->used > 0 && a->dp[a->used - 1] == 0)
    100          
212             {
213 29826091           --(a->used);
214             }
215             /* reset the sign flag if used == 0 */
216 84916264 100         if (a->used == 0)
217             {
218 63767           a->sign = PSTM_ZPOS;
219             }
220 84916264           }
221              
222             /******************************************************************************/
223             /**
224             Clear a big integer, and free associated working memory.
225             */
226 4769630           void pstm_clear(pstm_int *a)
227             {
228             int32 i;
229              
230             /* only do anything if a hasn't been freed previously */
231 4769630 50         if (a != NULL && a->dp != NULL)
    100          
232             {
233             /* first zero the digits */
234 38251275 100         for (i = 0; i < a->used; i++)
235             {
236 33705479           a->dp[i] = 0;
237             }
238 4545796           psFree(a->dp, a->pool);
239             /* reset members to make debugging easier */
240 4545796           a->dp = NULL;
241 4545796           a->alloc = a->used = 0;
242 4545796           a->sign = PSTM_ZPOS;
243             }
244 4769630           }
245              
246             /******************************************************************************/
247             /**
248             Clears mp0 - mp7, stopping at the first NULL mp value.
249             @example pstm_clear_multi(a, b, c, d, NULL, NULL, NULL, NULL);
250             */
251 1117973           void pstm_clear_multi(pstm_int *mp0, pstm_int *mp1, pstm_int *mp2,
252             pstm_int *mp3, pstm_int *mp4, pstm_int *mp5,
253             pstm_int *mp6, pstm_int *mp7)
254             {
255 1117973 50         if (mp0)
256             {
257 1117973           pstm_clear(mp0);
258 1117973 50         if (mp1)
259             {
260 1117973           pstm_clear(mp1);
261 1117973 100         if (mp2)
262             {
263 4477           pstm_clear(mp2);
264 4477 50         if (mp3)
265             {
266 0           pstm_clear(mp3);
267 0 0         if (mp4)
268             {
269 0           pstm_clear(mp4);
270 0 0         if (mp5)
271             {
272 0           pstm_clear(mp5);
273 0 0         if (mp6)
274             {
275 0           pstm_clear(mp6);
276 0 0         if (mp7)
277             {
278 0           pstm_clear(mp7);
279             }
280             }
281             }
282             }
283             }
284             }
285             }
286             }
287 1117973           }
288              
289             /******************************************************************************/
290             /**
291             Set to zero.
292             */
293 914205           void pstm_zero(pstm_int *a)
294             {
295             uint16_t n;
296             pstm_digit *tmp;
297              
298 914205           a->sign = PSTM_ZPOS;
299 914205           a->used = 0;
300              
301 914205           tmp = a->dp;
302 9220918 100         for (n = 0; n < a->alloc; n++)
303             {
304 8306713           *tmp++ = 0;
305             }
306 914205           }
307              
308             /******************************************************************************/
309             /*
310             Compare maginitude of two ints (unsigned).
311             */
312 48827887           int32_t pstm_cmp_mag(const pstm_int *a, const pstm_int *b)
313             {
314             uint16_t n;
315             const pstm_digit *tmpa, *tmpb;
316              
317             /* compare based on # of non-zero digits */
318 48827887 100         if (a->used > b->used)
319             {
320 469170           return PSTM_GT;
321             }
322 48358717 100         else if (a->used < b->used)
323             {
324 563774           return PSTM_LT;
325             }
326             /* alias for a */
327 47794943           tmpa = a->dp + (a->used - 1);
328             /* alias for b */
329 47794943           tmpb = b->dp + (a->used - 1);
330             /* compare based on digits */
331 48382323 100         for (n = 0; n < a->used; ++n, --tmpa, --tmpb)
332             {
333 48375109 100         if (*tmpa > *tmpb)
334             {
335 19643635           return PSTM_GT;
336             }
337 28731474 100         else if (*tmpa < *tmpb)
338             {
339 28144094           return PSTM_LT;
340             }
341             }
342 7214           return PSTM_EQ;
343             }
344              
345             /******************************************************************************/
346             /*
347             Compare two ints (signed)
348             */
349 15811863           int32_t pstm_cmp(const pstm_int *a, const pstm_int *b)
350             {
351             /* compare based on sign */
352 15811863 50         if (a->sign != b->sign)
353             {
354 0 0         if (a->sign == PSTM_NEG)
355             {
356 0           return PSTM_LT;
357             }
358             else
359             {
360 0           return PSTM_GT;
361             }
362             }
363             /* compare digits */
364 15811863 50         if (a->sign == PSTM_NEG)
365             {
366             /* if negative compare opposite direction */
367 0           return pstm_cmp_mag(b, a);
368             }
369             else
370             {
371 15811863           return pstm_cmp_mag(a, b);
372             }
373             }
374              
375             /******************************************************************************/
376             /*
377             compare against a single digit
378             */
379 7523409           int32_t pstm_cmp_d(const pstm_int *a, pstm_digit b)
380             {
381             /* compare based on sign */
382 7523409 100         if ((b && a->used == 0) || a->sign == PSTM_NEG)
    50          
    100          
383             {
384 3497506           return PSTM_LT;
385             }
386             /* compare based on magnitude */
387 4025903 100         if (a->used > 1)
388             {
389 3985823           return PSTM_GT;
390             }
391             /* compare the only digit of a to b */
392 40080 100         if (a->dp[0] > b)
393             {
394 37938           return PSTM_GT;
395             }
396 2142 50         else if (a->dp[0] < b)
397             {
398 0           return PSTM_LT;
399             }
400             else
401             {
402 2142           return PSTM_EQ;
403             }
404             }
405              
406             /******************************************************************************/
407             /*
408             pstm_ints can be initialized more precisely when they will populated
409             using pstm_read_unsigned_bin since the length of the byte stream is known
410             */
411 26513           int32_t pstm_init_for_read_unsigned_bin(psPool_t *pool, pstm_int *a, psSize_t len)
412             {
413             psSize_t size;
414              
415             /*
416             Need to set this based on how many words max it will take to store the bin.
417             The magic + 2:
418             1 to round up for the remainder of this integer math
419             1 for the initial carry of '1' bits that fall between DIGIT_BIT and 8
420             */
421 53026           size = (((len / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
422 26513           / DIGIT_BIT) + 2;
423 26513           return pstm_init_size(pool, a, size);
424             }
425              
426              
427             /******************************************************************************/
428             /**
429             Reads a unsigned char array into pstm_int format.
430              
431             @param[in,out] a The allocated and initialized pstm_int
432             @param[in] b Pointer to a byte array of length 'c'.
433             @param[in] c Length in bytes of 'b'.
434              
435             @pre User should have called pstm_init_for_read_unsigned_bin first.
436             @note There is some grow logic here if the default pstm_init was used but we
437             don't really want to hit it.
438             */
439 19896           int32_t pstm_read_unsigned_bin(pstm_int *a, const unsigned char *buf, psSize_t len)
440             {
441             /* zero the int */
442 19896           pstm_zero(a);
443              
444             /*
445             If we know the endianness of this architecture, and we're using
446             32-bit pstm_digits, we can optimize this
447             TODO Can optimize 64 bit case as well.
448             */
449             # if (defined(ENDIAN_LITTLE) || defined(ENDIAN_BIG)) && !defined(PSTM_64BIT)
450             /* But not for both simultaneously */
451             # if defined(ENDIAN_LITTLE) && defined(ENDIAN_BIG)
452             # error Both ENDIAN_LITTLE and ENDIAN_BIG defined.
453             # endif
454             {
455             unsigned char *pd;
456             int16_t slen;
457             if ((unsigned) len > (PSTM_MAX_SIZE * sizeof(pstm_digit)))
458             {
459             uint16_t excess = len - (PSTM_MAX_SIZE * sizeof(pstm_digit));
460             len -= excess;
461             buf += excess;
462             }
463             a->used = (len + sizeof(pstm_digit) - 1) / sizeof(pstm_digit);
464             if (a->alloc < a->used)
465             {
466             if (pstm_grow(a, a->used) != PSTM_OKAY)
467             {
468             return PSTM_MEM;
469             }
470             }
471             pd = (unsigned char *) a->dp;
472             /* read the bytes in */
473             /* these loops need len to go negative */
474             slen = (int16_t) len;
475             # ifdef ENDIAN_BIG
476             {
477             /* Use Duff's device to unroll the loop. */
478             uint16_t idx = (slen - 1) & ~3;
479             switch (slen % 4)
480             {
481             case 0: do
482             {
483             pd[idx + 0] = *buf++;
484             case 3: pd[idx + 1] = *buf++;
485             case 2: pd[idx + 2] = *buf++;
486             case 1: pd[idx + 3] = *buf++;
487             idx -= 4;
488             }
489             while ((slen -= 4) > 0);
490             }
491             }
492             # else
493             for (slen -= 1; slen >= 0; slen -= 1)
494             {
495             pd[slen] = *buf++;
496             }
497             # endif
498             }
499             # else
500             /* Big enough based on the len? */
501 39792           a->used = (((len / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
502 19896           / DIGIT_BIT) + 2;
503              
504 19896 50         if (a->alloc < a->used)
505             {
506 0 0         if (pstm_grow(a, a->used) != PSTM_OKAY)
507             {
508 0           return PSTM_MEM;
509             }
510             }
511             /* read the bytes in */
512 2762870 100         for (; len > 0; len--)
513             {
514 2742974 50         if (pstm_mul_2d(a, 8, a) != PSTM_OKAY)
515             {
516 0           return PS_MEM_FAIL;
517             }
518 2742974           a->dp[0] |= *buf++;
519 2742974           a->used += 1;
520             }
521             # endif
522              
523 19896           pstm_clamp(a);
524 19896           return PS_SUCCESS;
525             }
526              
527             /******************************************************************************/
528             /**
529             Read a pstm_int from an ASN.1 formatted buffer.
530             */
531 7686           int32_t pstm_read_asn(psPool_t *pool, const unsigned char **pp, psSize_t len,
532             pstm_int *a)
533             {
534 7686           const unsigned char *p = *pp;
535             psSize_t vlen;
536              
537 15242 50         if (len < 1 || *(p++) != ASN_INTEGER ||
538 15112 50         getAsnLength(&p, len - 1, &vlen) < 0 || (len - 1) < vlen)
539             {
540 130           return PS_PARSE_FAIL;
541             }
542             /* Make a smart size since we know the length */
543 7556 50         if (pstm_init_for_read_unsigned_bin(pool, a, vlen) != PSTM_OKAY)
544             {
545 0           return PS_MEM_FAIL;
546             }
547 7556 50         if (pstm_read_unsigned_bin(a, p, vlen) != 0)
548             {
549 0           pstm_clear(a);
550             psTraceCrypto("pstm_read_asn failed\n");
551 0           return PS_PARSE_FAIL;
552             }
553 7556           *pp = p + vlen;
554 7686           return PS_SUCCESS;
555             }
556              
557             # if defined USE_ECC || defined USE_DH || defined USE_CERT_GEN
558              
559             /******************************************************************************/
560             /**
561             Add a digit to an int.
562             c = a + b;
563             @param[in] pool Memory pool
564             @param[in] a Big integer operand
565             @param[in] b Big digit operand
566             @param[out] c Big integer result
567             @return < 0 on failure
568             */
569 862516           int32_t pstm_add_d(psPool_t *pool, const pstm_int *a, pstm_digit b, pstm_int *c)
570             {
571             pstm_int tmp;
572             int32_t res;
573              
574 862516 50         if (pstm_init_size(pool, &tmp, sizeof(pstm_digit)) != PSTM_OKAY)
575             {
576 0           return PS_MEM_FAIL;
577             }
578 862516           pstm_set(&tmp, b);
579 862516           res = pstm_add(a, &tmp, c);
580 862516           pstm_clear(&tmp);
581 862516           return res;
582             }
583              
584             # endif /* defined USE_ECC || defined USE_DH || defined USE_CERT_GEN */
585             # ifdef USE_ECC
586              
587             /* chars used in radix (base) conversions */
588             const static unsigned char pstm_s_rmap[64] =
589             "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
590              
591             /**
592             Read in data to a big integer from a given radix (base).
593             Typically the radix is 16 to read in hex data.
594             To read in binary data, use pstm_read_unsigned_bin().
595              
596             @param[in] pool Memory pool
597             @param[out] a Big integer to load data into.
598             @param[in] buf ASCII buffer containing 'len' bytes of data to read.
599             @param[in] len Number of bytes to read from 'buf'.
600             @param[in] radix Base of data in each byte of buf. Typically 16.
601             @return < 0 on failure
602             */
603 6665           int32_t pstm_read_radix(psPool_t *pool, pstm_int *a,
604             const char *buf, psSize_t len, uint8_t radix)
605             {
606             int32_t y, neg;
607             unsigned char ch;
608              
609             /* make sure the radix is ok */
610 6665 50         if (radix < 2 || radix > 64)
    50          
611             {
612 0           return PS_ARG_FAIL;
613             }
614              
615             /* if the leading digit is a minus set the sign to negative. */
616 6665 50         if (*buf == '-')
617             {
618 0           ++buf; len--;
619 0           neg = PSTM_NEG;
620             }
621             else
622             {
623 6665           neg = PSTM_ZPOS;
624             }
625              
626             /* set the integer to the default of zero */
627 6665           pstm_zero(a);
628              
629             /* process each digit of the string */
630 869181 100         while (len > 0)
631             {
632             /*
633             if the radix < 36 the conversion is case insensitive this allows
634             numbers like 1AB and 1ab to represent the same value [e.g. in hex].
635             */
636 868964 50         ch = ((radix < 36) ?
637 868964           (unsigned char) toupper((unsigned char) *buf) :
638 0           (unsigned char) *buf);
639 11830276 100         for (y = 0; y < 64; y++)
640             {
641 11823828 100         if (ch == pstm_s_rmap[y])
642             {
643 862516           break;
644             }
645             }
646              
647             /*
648             if the char was found in the map and is less than the given
649             radix, add it to the number, otherwise exit the loop.
650             */
651 868964 100         if (y < radix)
652             {
653 862516           pstm_mul_d(a, (pstm_digit) radix, a);
654 862516           pstm_add_d(pool, a, (pstm_digit) y, a);
655             }
656             else
657             {
658 6448           break;
659             }
660 862516           ++buf; len--;
661             }
662              
663             /* set the sign only if a != 0 */
664 6665 50         if (pstm_iszero(a) != PS_TRUE)
665             {
666 6665           a->sign = neg;
667             }
668 6665           return PS_SUCCESS;
669             }
670             # endif /* USE_ECC */
671              
672             /******************************************************************************/
673              
674 67007           uint16_t pstm_count_bits(const pstm_int *a)
675             {
676             int16 r;
677             pstm_digit q;
678              
679 67007 50         if (a->used == 0)
680             {
681 0           return 0;
682             }
683              
684             /* get number of digits and add that */
685 67007           r = (a->used - 1) * DIGIT_BIT;
686              
687             /* take the last digit and count the bits in it */
688 67007           q = a->dp[a->used - 1];
689 2385123 100         while (q > ((pstm_digit) 0))
690             {
691 2318116           ++r;
692 2318116           q >>= ((pstm_digit) 1);
693             }
694 67007           return r;
695             }
696              
697             /******************************************************************************/
698              
699 20679           uint16_t pstm_unsigned_bin_size(const pstm_int *a)
700             {
701 20679           psSize_t size = pstm_count_bits(a);
702              
703 20679           return size / 8 + ((size & 7) != 0 ? 1 : 0);
704             }
705              
706             /******************************************************************************/
707             /**
708             a = b, where b is a single digit.
709             */
710 885502           void pstm_set(pstm_int *a, pstm_digit b)
711             {
712 885502           pstm_zero(a);
713 885502           a->dp[0] = b;
714 885502           a->used = a->dp[0] ? 1 : 0;
715 885502           }
716              
717             /******************************************************************************/
718             /**
719             Right shift 'a' by 'b' digits.
720             @note This is not a bit shift.
721             */
722 0           void pstm_rshd(pstm_int *a, uint16_t b)
723             {
724             uint16_t y;
725              
726             /* too many digits just zero and return */
727 0 0         if (b >= a->used)
728             {
729 0           pstm_zero(a);
730 0           return;
731             }
732             /* shift */
733 0 0         for (y = 0; y < a->used - b; y++)
734             {
735 0           a->dp[y] = a->dp[y + b];
736             }
737             /* zero the rest */
738 0 0         for (; y < a->used; y++)
739             {
740 0           a->dp[y] = 0;
741             }
742             /* decrement count */
743 0           a->used -= b;
744 0           pstm_clamp(a);
745             }
746              
747             /******************************************************************************/
748             /**
749             Left shift 'a' by 'b' digits.
750             This will grow 'a', and possibly cause a reallocation of memory.
751             @note This is not a bit shift.
752             */
753 24330           int32_t pstm_lshd(pstm_int *a, uint16_t b)
754             {
755             uint16_t x;
756             int32_t res;
757              
758             /* If its less than zero return. */
759 24330 50         if (b <= 0)
760             {
761 0           return PSTM_OKAY;
762             }
763             /* Grow to fit the new digits. */
764 24330 50         if (a->alloc < a->used + b)
765             {
766 0 0         if ((res = pstm_grow(a, a->used + b)) != PSTM_OKAY)
767             {
768 0           return res;
769             }
770             }
771              
772             {
773             register pstm_digit *top, *bottom;
774             /* Increment the used by the shift amount then copy upwards. */
775 24330           a->used += b;
776             /* top */
777 24330           top = a->dp + a->used - 1;
778             /* base */
779 24330           bottom = a->dp + a->used - 1 - b;
780             /*
781             This is implemented using a sliding window except the window goes the
782             other way around. Copying from the bottom to the top.
783             */
784 169347 100         for (x = a->used - 1; x >= b; x--)
785             {
786 145017           *top-- = *bottom--;
787             }
788             /* zero the lower digits */
789 24330           top = a->dp;
790 265800 100         for (x = 0; x < b; x++)
791             {
792 241470           *top++ = 0;
793             }
794             }
795 24330           return PSTM_OKAY;
796             }
797              
798             /******************************************************************************/
799             /**
800             a = 2**b.
801             */
802 2142           int32_t pstm_2expt(pstm_int *a, int16_t b)
803             {
804             uint16_t z;
805              
806             /* zero a as per default */
807 2142           pstm_zero(a);
808              
809 2142 50         if (b < 0)
810             {
811 0           return PSTM_OKAY;
812             }
813              
814 2142           z = b / DIGIT_BIT;
815 2142 50         if (z >= PSTM_MAX_SIZE)
816             {
817 0           return PS_LIMIT_FAIL;
818             }
819              
820             /* set the used count of where the bit will go */
821 2142           a->used = z + 1;
822              
823 2142 50         if (a->used > a->alloc)
824             {
825 0 0         if (pstm_grow(a, a->used) != PSTM_OKAY)
826             {
827 0           return PS_MEM_FAIL;
828             }
829             }
830              
831             /* put the single bit in its place */
832 2142           a->dp[z] = ((pstm_digit) 1) << (b % DIGIT_BIT);
833 2142           return PSTM_OKAY;
834             }
835              
836             /******************************************************************************/
837             /**
838             b = a * 2.
839             Implements multiplication as a left shift of all digits of 1 bit.
840             */
841 119842           int32_t pstm_mul_2(const pstm_int *a, pstm_int *b)
842             {
843             int32 res;
844             int16 x, oldused;
845              
846             /* grow to accomodate result */
847 119842 50         if (b->alloc < a->used + 1)
848             {
849 0 0         if ((res = pstm_grow(b, a->used + 1)) != PSTM_OKAY)
850             {
851 0           return res;
852             }
853             }
854 119842           oldused = b->used;
855 119842           b->used = a->used;
856              
857             {
858             register pstm_digit r, rr, *tmpa, *tmpb;
859              
860             /* alias for source */
861 119842           tmpa = a->dp;
862              
863             /* alias for dest */
864 119842           tmpb = b->dp;
865              
866             /* carry */
867 119842           r = 0;
868 256810 100         for (x = 0; x < a->used; x++)
869             {
870             /*
871             get what will be the *next* carry bit from the
872             MSB of the current digit
873             */
874 136968           rr = *tmpa >> ((pstm_digit) (DIGIT_BIT - 1));
875             /* now shift up this digit, add in the carry [from the previous] */
876 136968           *tmpb++ = ((*tmpa++ << ((pstm_digit) 1)) | r);
877             /*
878             copy the carry that would be from the source
879             digit into the next iteration
880             */
881 136968           r = rr;
882             }
883              
884             /* new leading digit? */
885 119842 100         if (r != 0 && b->used != (PSTM_MAX_SIZE - 1))
    50          
886             {
887             /* add a MSB which is always 1 at this point */
888 2           *tmpb = 1;
889 2           ++(b->used);
890             }
891             /* now zero any excess digits on the destination that we didn't write to */
892 119842           tmpb = b->dp + b->used;
893 119842 50         for (x = b->used; x < oldused; x++)
894             {
895 0           *tmpb++ = 0;
896             }
897             }
898 119842           b->sign = a->sign;
899 119842           return PSTM_OKAY;
900             }
901              
902             /******************************************************************************/
903             /**
904             unsigned subtraction ||a|| must be >= ||b|| ALWAYS.
905             c = a - b.
906             */
907 20061548           int32_t pstm_sub_s(const pstm_int *a, const pstm_int *b, pstm_int *c)
908             {
909             int16 oldbused, oldused;
910             int32 x;
911             pstm_word t;
912              
913 20061548 50         if (b->used > a->used)
914             {
915 0           return PS_LIMIT_FAIL;
916             }
917 20061548 50         if (c->alloc < a->used)
918             {
919 0 0         if ((x = pstm_grow(c, a->used)) != PSTM_OKAY)
920             {
921 0           return x;
922             }
923             }
924 20061548           oldused = c->used;
925 20061548           oldbused = b->used;
926 20061548           c->used = a->used;
927 20061548           t = 0;
928 232097284 100         for (x = 0; x < oldbused; x++)
929             {
930 212035736           t = ((pstm_word) a->dp[x]) - (((pstm_word) b->dp[x]) + t);
931 212035736           c->dp[x] = (pstm_digit) t;
932 212035736           t = (t >> DIGIT_BIT) & 1;
933             }
934 21130634 100         for (; x < a->used; x++)
935             {
936 1069086           t = ((pstm_word) a->dp[x]) - t;
937 1069086           c->dp[x] = (pstm_digit) t;
938 1069086           t = (t >> DIGIT_BIT);
939             }
940 20080482 100         for (; x < oldused; x++)
941             {
942 18934           c->dp[x] = 0;
943             }
944 20061548           pstm_clamp(c);
945 20061548           return PSTM_OKAY;
946             }
947              
948             /******************************************************************************/
949              
950             /**
951             Unsigned addition of two big integers.
952             c = a + b;
953             @param[in] pool Memory pool
954             @param[in] a Big integer operand
955             @param[in] b Big integer operand
956             @param[out] c Big integer result
957             @return < 0 on failure
958             */
959 13005772           static int32_t s_pstm_add(const pstm_int *a, const pstm_int *b, pstm_int *c)
960             {
961             int16 x, y, oldused;
962             register pstm_word t, adp, bdp;
963              
964 13005772           y = a->used;
965 13005772 100         if (b->used > y)
966             {
967 38528           y = b->used;
968             }
969 13005772           oldused = c->used;
970 13005772           c->used = y;
971              
972 13005772 50         if (c->used > c->alloc)
973             {
974 0 0         if (pstm_grow(c, c->used) != PSTM_OKAY)
975             {
976 0           return PS_MEM_FAIL;
977             }
978             }
979              
980 13005772           t = 0;
981 137887622 100         for (x = 0; x < y; x++)
982             {
983 124881850 50         if (a->used <= x)
984             {
985 0           adp = 0;
986             }
987             else
988             {
989 124881850           adp = (pstm_word) a->dp[x];
990             }
991 124881850 100         if (b->used <= x)
992             {
993 26565442           bdp = 0;
994             }
995             else
996             {
997 98316408           bdp = (pstm_word) b->dp[x];
998             }
999 124881850           t += (adp) + (bdp);
1000 124881850           c->dp[x] = (pstm_digit) t;
1001 124881850           t >>= DIGIT_BIT;
1002             }
1003 13005772 100         if (t != 0 && x < PSTM_MAX_SIZE)
    50          
1004             {
1005 3429 50         if (c->used == c->alloc)
1006             {
1007 0 0         if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY)
1008             {
1009 0           return PS_MEM_FAIL;
1010             }
1011             }
1012 3429           c->dp[c->used++] = (pstm_digit) t;
1013 3429           ++x;
1014             }
1015              
1016 13005772           c->used = x;
1017 13005806 100         for (; x < oldused; x++)
1018             {
1019 34           c->dp[x] = 0;
1020             }
1021 13005772           pstm_clamp(c);
1022 13005772           return PSTM_OKAY;
1023             }
1024              
1025              
1026             /******************************************************************************/
1027             /**
1028             Signed subtraction. a and b can be any value.
1029             c = a - b.
1030             */
1031 17327199           int32_t pstm_sub(const pstm_int *a, const pstm_int *b, pstm_int *c)
1032             {
1033             int32 res;
1034             int16 sa, sb;
1035              
1036 17327199           sa = a->sign;
1037 17327199           sb = b->sign;
1038              
1039 17327199 100         if (sa != sb)
1040             {
1041             /*
1042             subtract a negative from a positive, OR a positive from a negative.
1043             For both, ADD their magnitudes, and use the sign of the first number.
1044             */
1045 769020           c->sign = sa;
1046 769020 50         if ((res = s_pstm_add(a, b, c)) != PSTM_OKAY)
1047             {
1048 0           return res;
1049             }
1050             }
1051             else
1052             {
1053             /*
1054             subtract a positive from a positive, OR a negative from a negative.
1055             First, take the difference between their magnitudes, then...
1056             */
1057 16558179 100         if (pstm_cmp_mag(a, b) != PSTM_LT)
1058             {
1059             /* Copy the sign from the first */
1060 12291989           c->sign = sa;
1061             /* The first has a larger or equal magnitude */
1062 12291989 50         if ((res = pstm_sub_s(a, b, c)) != PSTM_OKAY)
1063             {
1064 0           return res;
1065             }
1066             }
1067             else
1068             {
1069             /* The result has the _opposite_ sign from the first number. */
1070 4266190           c->sign = (sa == PSTM_ZPOS) ? PSTM_NEG : PSTM_ZPOS;
1071             /* The second has a larger magnitude */
1072 4266190 50         if ((res = pstm_sub_s(b, a, c)) != PSTM_OKAY)
1073             {
1074 0           return res;
1075             }
1076             }
1077             }
1078 17327199           return PS_SUCCESS;
1079             }
1080              
1081             /******************************************************************************/
1082             /**
1083             Signed subtraction of a digit.
1084             c = a - b, where b is a digit
1085             */
1086 0           int32_t pstm_sub_d(psPool_t *pool, const pstm_int *a, pstm_digit b, pstm_int *c)
1087             {
1088             pstm_int tmp;
1089             int32_t res;
1090              
1091 0 0         if (pstm_init_size(pool, &tmp, sizeof(pstm_digit)) != PSTM_OKAY)
1092             {
1093 0           return PS_MEM_FAIL;
1094             }
1095 0           pstm_set(&tmp, b);
1096 0           res = pstm_sub(a, &tmp, c);
1097 0           pstm_clear(&tmp);
1098 0           return res;
1099             }
1100              
1101             /******************************************************************************/
1102             /**
1103             Sets up the montgomery reduction.
1104             Fast inversion mod 2**k
1105             Based on the fact that
1106             XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
1107             => 2*X*A - X*X*A*A = 1
1108             => 2*(1) - (1) = 1
1109             @param[in] a
1110             @param[out] rho
1111             */
1112 2142           int32_t pstm_montgomery_setup(const pstm_int *a, pstm_digit *rho)
1113             {
1114             pstm_digit x, b;
1115              
1116 2142           b = a->dp[0];
1117              
1118 2142 50         if ((b & 1) == 0)
1119             {
1120             psTraceCrypto("pstm_montogomery_setup failure\n");
1121 0           return PS_ARG_FAIL;
1122             }
1123              
1124 2142           x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
1125 2142           x *= 2 - b * x; /* here x*a==1 mod 2**8 */
1126 2142           x *= 2 - b * x; /* here x*a==1 mod 2**16 */
1127 2142           x *= 2 - b * x; /* here x*a==1 mod 2**32 */
1128             # ifdef PSTM_64BIT
1129 2142           x *= 2 - b * x; /* here x*a==1 mod 2**64 */
1130             # endif
1131             /* rho = -1/m mod b */
1132 2142           *rho = (pstm_digit) (((pstm_word) 1 << ((pstm_word) DIGIT_BIT)) -
1133             ((pstm_word) x));
1134 2142           return PSTM_OKAY;
1135             }
1136              
1137             /******************************************************************************/
1138             /**
1139             a = B**n mod b, without division or multiplication.
1140             Useful for normalizing numbers in a Montgomery system.
1141             */
1142 2142           int32_t pstm_montgomery_calc_normalization(pstm_int *a, const pstm_int *b)
1143             {
1144             uint16_t x, bits;
1145              
1146             /* how many bits of last digit does b use */
1147 2142           bits = pstm_count_bits(b) % DIGIT_BIT;
1148 2142 100         if (!bits)
1149             {
1150 2           bits = DIGIT_BIT;
1151             }
1152              
1153             /* compute A = B^(n-1) * 2^(bits-1) */
1154 2142 50         if (b->used > 1)
1155             {
1156 2142 50         if ((x = pstm_2expt(a, (b->used - 1) * DIGIT_BIT + bits - 1)) !=
1157             PSTM_OKAY)
1158             {
1159 0           return x;
1160             }
1161             }
1162             else
1163             {
1164 0           pstm_set(a, 1);
1165 0           bits = 1;
1166             }
1167             /* now compute C = A * B mod b */
1168 121984 100         for (x = bits - 1; x < (uint16_t) DIGIT_BIT; x++)
1169             {
1170 119842 50         if (pstm_mul_2(a, a) != PSTM_OKAY)
1171             {
1172 0           return PS_MEM_FAIL;
1173             }
1174 119842 100         if (pstm_cmp_mag(a, b) != PSTM_LT)
1175             {
1176 2142 50         if (pstm_sub_s(a, b, a) != PSTM_OKAY)
1177             {
1178 0           return PS_MEM_FAIL;
1179             }
1180             }
1181             }
1182 2142           return PSTM_OKAY;
1183             }
1184              
1185             /******************************************************************************/
1186             /*
1187             c = a * 2**d
1188             */
1189 2775864           static int32_t pstm_mul_2d(const pstm_int *a, int16_t b, pstm_int *c)
1190             {
1191             pstm_digit carry, carrytmp, shift;
1192             uint16_t x;
1193              
1194             /* copy it */
1195 2775864 50         if (pstm_copy(a, c) != PSTM_OKAY)
1196             {
1197 0           return PS_MEM_FAIL;
1198             }
1199              
1200             /* handle whole digits */
1201 2775864 100         if (b >= DIGIT_BIT)
1202             {
1203 24330 50         if (pstm_lshd(c, b / DIGIT_BIT) != PSTM_OKAY)
1204             {
1205 0           return PS_MEM_FAIL;
1206             }
1207             }
1208 2775864           b %= DIGIT_BIT;
1209              
1210             /* shift the digits */
1211 2775864 100         if (b != 0)
1212             {
1213 2775758           carry = 0;
1214 2775758           shift = DIGIT_BIT - b;
1215 45638350 100         for (x = 0; x < c->used; x++)
1216             {
1217 42862592           carrytmp = c->dp[x] >> shift;
1218 42862592           c->dp[x] = (c->dp[x] << b) + carry;
1219 42862592           carry = carrytmp;
1220             }
1221             /* store last carry if room */
1222 2775758 100         if (carry && x < PSTM_MAX_SIZE)
    50          
1223             {
1224 3612 50         if (c->used == c->alloc)
1225             {
1226 0 0         if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY)
1227             {
1228 0           return PS_MEM_FAIL;
1229             }
1230             }
1231 3612           c->dp[c->used++] = carry;
1232             }
1233             }
1234 2775864           pstm_clamp(c);
1235 2775864           return PSTM_OKAY;
1236             }
1237              
1238             /******************************************************************************/
1239             /*
1240             c = a mod 2**b
1241             */
1242 0           static int32_t pstm_mod_2d(const pstm_int *a, int16_t b, pstm_int *c)
1243             {
1244             uint16_t x;
1245              
1246             /* zero if count less than or equal to zero */
1247 0 0         if (b <= 0)
1248             {
1249 0           pstm_zero(c);
1250 0           return PSTM_OKAY;
1251             }
1252              
1253             /* get copy of input */
1254 0 0         if (pstm_copy(a, c) != PSTM_OKAY)
1255             {
1256 0           return PS_MEM_FAIL;
1257             }
1258              
1259             /* if 2**d is larger than we just return */
1260 0 0         if (b >= (DIGIT_BIT * a->used))
1261             {
1262 0           return PSTM_OKAY;
1263             }
1264              
1265             /* zero digits above the last digit of the modulus */
1266 0 0         for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++)
1267             {
1268 0           c->dp[x] = 0;
1269             }
1270             /* clear the digit that is not completely outside/inside the modulus */
1271 0           c->dp[b / DIGIT_BIT] &= ~((pstm_digit) 0) >> (DIGIT_BIT - b);
1272 0           pstm_clamp(c);
1273 0           return PSTM_OKAY;
1274             }
1275              
1276              
1277             /******************************************************************************/
1278             /*
1279             c = a * b
1280             */
1281 862516           int32_t pstm_mul_d(const pstm_int *a, const pstm_digit b, pstm_int *c)
1282             {
1283             pstm_word w;
1284             int32 res;
1285             int16 x, oldused;
1286              
1287 862516 50         if (c->alloc < a->used + 1)
1288             {
1289 0 0         if ((res = pstm_grow(c, a->used + 1)) != PSTM_OKAY)
1290             {
1291 0           return res;
1292             }
1293             }
1294 862516           oldused = c->used;
1295 862516           c->used = a->used;
1296 862516           c->sign = a->sign;
1297 862516           w = 0;
1298 4740280 100         for (x = 0; x < a->used; x++)
1299             {
1300 3877764           w = ((pstm_word) a->dp[x]) * ((pstm_word) b) + w;
1301 3877764           c->dp[x] = (pstm_digit) w;
1302 3877764           w = w >> DIGIT_BIT;
1303             }
1304 862516 100         if (w != 0 && (a->used != PSTM_MAX_SIZE))
    50          
1305             {
1306 52615           c->dp[c->used++] = (pstm_digit) w;
1307 52615           ++x;
1308             }
1309 862516 50         for (; x < oldused; x++)
1310             {
1311 0           c->dp[x] = 0;
1312             }
1313 862516           pstm_clamp(c);
1314 862516           return PSTM_OKAY;
1315             }
1316              
1317             /******************************************************************************/
1318             /**
1319             c = a / 2**b, d = remainder.
1320             @param[in] pool Memory pool
1321             @param[in] a Numerator
1322             @param[in] b Exponent of the denominator
1323             @param[out] c The result
1324             @param[out] d If non-NULL, the remainder of the division is stored here
1325             @return < 0 on failure
1326             */
1327 18025665           int32_t pstm_div_2d(psPool_t *pool, const pstm_int *a, int16_t b, pstm_int *c,
1328             pstm_int *d)
1329             {
1330             pstm_digit D, r, rr;
1331             int32 res;
1332             int16 x;
1333              
1334             /* if the shift count is <= 0 then we do no work */
1335 18025665 50         if (b <= 0)
1336             {
1337 0 0         if (pstm_copy(a, c) != PSTM_OKAY)
1338             {
1339 0           return PS_MEM_FAIL;
1340             }
1341 0 0         if (d != NULL)
1342             {
1343 0           pstm_zero(d);
1344             }
1345 0           return PSTM_OKAY;
1346             }
1347             /* copy */
1348 18025665 50         if (pstm_copy(a, c) != PSTM_OKAY)
1349             {
1350 0           res = PS_MEM_FAIL;
1351 0           goto LBL_DONE;
1352             }
1353              
1354             /* shift by as many digits in the bit count */
1355 18025665 50         if (b >= (int16_t) DIGIT_BIT)
1356             {
1357 0           pstm_rshd(c, b / DIGIT_BIT);
1358             }
1359              
1360             /* shift any bit count < DIGIT_BIT */
1361 18025665           D = (pstm_digit) (b % DIGIT_BIT);
1362 18025665 50         if (D != 0)
1363             {
1364             register pstm_digit *tmpc, mask, shift;
1365              
1366             /* mask */
1367 18025665           mask = (((pstm_digit) 1) << D) - 1;
1368              
1369             /* shift for lsb */
1370 18025665           shift = DIGIT_BIT - D;
1371              
1372             /* alias */
1373 18025665           tmpc = c->dp + (c->used - 1);
1374              
1375             /* carry */
1376 18025665           r = 0;
1377 233084784 100         for (x = c->used - 1; x >= 0; x--)
1378             {
1379             /* get the lower bits of this word in a temp */
1380 215059119           rr = *tmpc & mask;
1381              
1382             /* shift the current word and mix in the carry bits from previous */
1383 215059119           *tmpc = (*tmpc >> D) | (r << shift);
1384 215059119           --tmpc;
1385              
1386             /* set the carry to the carry bits of the current word above */
1387 215059119           r = rr;
1388             }
1389             }
1390 18025665           pstm_clamp(c);
1391              
1392 18025665           res = PSTM_OKAY;
1393             LBL_DONE:
1394             /* set the remainder */
1395 18025665 50         if (d != NULL)
1396             {
1397 0 0         if (pstm_mod_2d(a, b, d) != PSTM_OKAY)
1398             {
1399 0           res = PS_MEM_FAIL;
1400             }
1401             }
1402 18025665           return res;
1403             }
1404              
1405             /******************************************************************************/
1406             /**
1407             b = a / 2.
1408             Implemented as a right shift of one bit.
1409             */
1410 4502187           int32_t pstm_div_2(const pstm_int *a, pstm_int *b)
1411             {
1412             int16 x, oldused;
1413              
1414 4502187 50         if (b->alloc < a->used)
1415             {
1416 0 0         if (pstm_grow(b, a->used) != PSTM_OKAY)
1417             {
1418 0           return PS_MEM_FAIL;
1419             }
1420             }
1421 4502187           oldused = b->used;
1422 4502187           b->used = a->used;
1423             {
1424             register pstm_digit r, rr, *tmpa, *tmpb;
1425              
1426             /* source alias */
1427 4502187           tmpa = a->dp + b->used - 1;
1428              
1429             /* dest alias */
1430 4502187           tmpb = b->dp + b->used - 1;
1431              
1432             /* carry */
1433 4502187           r = 0;
1434 37177229 100         for (x = b->used - 1; x >= 0; x--)
1435             {
1436             /* get the carry for the next iteration */
1437 32675042           rr = *tmpa & 1;
1438              
1439             /* shift the current digit, add in carry and store */
1440 32675042           *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1441              
1442             /* forward carry to next iteration */
1443 32675042           r = rr;
1444             }
1445              
1446             /* zero excess digits */
1447 4502187           tmpb = b->dp + b->used;
1448 4502187 50         for (x = b->used; x < oldused; x++)
1449             {
1450 0           *tmpb++ = 0;
1451             }
1452             }
1453 4502187           b->sign = a->sign;
1454 4502187           pstm_clamp(b);
1455 4502187           return PSTM_OKAY;
1456             }
1457              
1458             /******************************************************************************/
1459             /*
1460             Creates "a" then copies b into it
1461             */
1462 27302           int32_t pstm_init_copy(psPool_t *pool, pstm_int *a, const pstm_int *b,
1463             uint8_t toSqr)
1464             {
1465             int32_t res;
1466             uint16_t x;
1467              
1468 27302 50         if (a == b)
1469             {
1470 0           return PSTM_OKAY;
1471             }
1472 27302           x = b->alloc;
1473              
1474 27302 50         if (toSqr)
1475             {
1476             /*
1477             Smart-size: Increasing size of a if b->used is roughly half
1478             of b->alloc because usage has shown that a lot of these copies
1479             go on to be squared and need these extra digits
1480             */
1481 0 0         if ((b->used * 2) + 2 >= x)
1482             {
1483 0           x = (b->used * 2) + 3;
1484             }
1485             }
1486 27302 50         if ((res = pstm_init_size(pool, a, x)) != PSTM_OKAY)
1487             {
1488 0           return res;
1489             }
1490 27302           return pstm_copy(b, a);
1491             }
1492              
1493             /******************************************************************************/
1494             /*
1495             With some compilers, we have seen issues linking with the builtin
1496             64 bit division routine. The issues with either manifest in a failure
1497             to find 'udivdi3' at link time, or a runtime invalid instruction fault
1498             during an RSA operation.
1499             The routine below divides a 64 bit unsigned int by a 32 bit unsigned int
1500             explicitly, rather than using the division operation
1501             The 64 bit result is placed in the 'numerator' parameter
1502             The 32 bit mod (remainder) of the division is the return parameter
1503             Based on implementations by:
1504             Copyright (C) 2003 Bernardo Innocenti
1505             Copyright (C) 1999 Hewlett-Packard Co
1506             Copyright (C) 1999 David Mosberger-Tang
1507             */
1508             # if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1509             static uint32 psDiv64(uint64 *numerator, uint32 denominator)
1510             {
1511             uint64 rem = *numerator;
1512             uint64 b = denominator;
1513             uint64 res = 0;
1514             uint64 d = 1;
1515             uint32 high = rem >> 32;
1516              
1517             if (high >= denominator)
1518             {
1519             high /= denominator;
1520             res = (uint64) high << 32;
1521             rem -= (uint64) (high * denominator) << 32;
1522             }
1523             while ((int64) b > 0 && b < rem)
1524             {
1525             b = b + b;
1526             d = d + d;
1527             }
1528             do
1529             {
1530             if (rem >= b)
1531             {
1532             rem -= b;
1533             res += d;
1534             }
1535             b >>= 1;
1536             d >>= 1;
1537             }
1538             while (d);
1539             *numerator = res;
1540             return rem;
1541             }
1542             # endif /* USE_MATRIX_DIV64 */
1543              
1544             # if defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1545             typedef unsigned long uint128 __attribute__ ((mode(TI)));
1546             static uint64 psDiv128(uint128 *numerator, uint64 denominator)
1547             {
1548             uint128 rem = *numerator;
1549             uint128 b = denominator;
1550             uint128 res = 0;
1551             uint128 d = 1;
1552             uint64 high = rem >> 64;
1553              
1554             if (high >= denominator)
1555             {
1556             high /= denominator;
1557             res = (uint128) high << 64;
1558             rem -= (uint128) (high * denominator) << 64;
1559             }
1560             while ((uint128) b > 0 && b < rem)
1561             {
1562             b = b + b;
1563             d = d + d;
1564             }
1565             do
1566             {
1567             if (rem >= b)
1568             {
1569             rem -= b;
1570             res += d;
1571             }
1572             b >>= 1;
1573             d >>= 1;
1574             }
1575             while (d);
1576             *numerator = res;
1577             return rem;
1578             }
1579             # endif /* USE_MATRIX_DIV128 */
1580              
1581             # ifndef PSTM_LARGE_DIV
1582              
1583             /* This version of division uses short & small function, but offers
1584             bit worse performance than some others. */
1585 18587           int32_t pstm_div(psPool_t *pool, const pstm_int *a, const pstm_int *b,
1586             pstm_int *c, pstm_int *d)
1587             {
1588             pstm_int ta, tb, tq, q;
1589             int res, n, n2;
1590              
1591             /* is divisor zero ? */
1592 18587 50         if (pstm_iszero(b) == PSTM_YES)
1593             {
1594 0           return PS_LIMIT_FAIL;
1595             }
1596              
1597             /* if a < b then q=0, r = a */
1598 18587 100         if (pstm_cmp_mag(a, b) == PSTM_LT)
1599             {
1600 2142 50         if (d != NULL)
1601             {
1602 2142           res = pstm_copy(a, d);
1603             }
1604             else
1605             {
1606 0           res = PSTM_OKAY;
1607             }
1608 2142 50         if (c != NULL)
1609             {
1610 0           pstm_zero(c);
1611             }
1612 2142           return res;
1613             }
1614              
1615             /* init our temps */
1616 16445           res = pstm_init(pool, &ta);
1617 16445 50         if (res != PSTM_OKAY)
1618             {
1619 0           return res;
1620             }
1621 16445           res = pstm_init(pool, &tb);
1622 16445 50         if (res != PSTM_OKAY)
1623             {
1624 0           pstm_clear(&ta);
1625 0           return res;
1626             }
1627 16445           res = pstm_init(pool, &tq);
1628 16445 50         if (res != PSTM_OKAY)
1629             {
1630 0           pstm_clear(&ta);
1631 0           pstm_clear(&tb);
1632 0           return res;
1633             }
1634 16445           res = pstm_init(pool, &q);
1635 16445 50         if (res != PSTM_OKAY)
1636             {
1637 0           pstm_clear(&ta);
1638 0           pstm_clear(&tb);
1639 0           pstm_clear(&tq);
1640 0           return res;
1641             }
1642              
1643 16445           pstm_set(&tq, 1);
1644 16445           n = pstm_count_bits(a) - pstm_count_bits(b);
1645 16445 50         if (((res = pstm_abs(a, &ta)) != PSTM_OKAY) ||
    50          
1646 16445 50         ((res = pstm_abs(b, &tb)) != PSTM_OKAY) ||
1647 16445 50         ((res = pstm_mul_2d(&tb, n, &tb)) != PSTM_OKAY) ||
1648 16445           ((res = pstm_mul_2d(&tq, n, &tq)) != PSTM_OKAY))
1649             {
1650             goto LBL_ERR;
1651             }
1652              
1653 8241102 100         while (n-- >= 0)
1654             {
1655 8224657 100         if (pstm_cmp(&tb, &ta) != PSTM_GT)
1656             {
1657 4152279 50         if (((res = pstm_sub(&ta, &tb, &ta)) != PSTM_OKAY) ||
    50          
1658             ((res = pstm_add(&q, &tq, &q)) != PSTM_OKAY))
1659             {
1660             goto LBL_ERR;
1661             }
1662             }
1663 8224657 50         if (((res = pstm_div_2d(pool, &tb, 1, &tb, NULL)) !=
1664 8224657 50         PSTM_OKAY) ||
1665             ((res = pstm_div_2d(pool, &tq, 1, &tq, NULL)) !=
1666             PSTM_OKAY))
1667             {
1668             goto LBL_ERR;
1669             }
1670             }
1671              
1672             /* now q == quotient and ta == remainder */
1673 16445           n = a->sign;
1674 16445           n2 = (a->sign == b->sign) ? PSTM_ZPOS : PSTM_NEG;
1675 16445 50         if (c != NULL)
1676             {
1677 0           pstm_exch(c, &q);
1678 0 0         c->sign = (pstm_iszero(c) == PSTM_YES) ? PSTM_ZPOS : n2;
1679             }
1680 16445 50         if (d != NULL)
1681             {
1682 16445           pstm_exch(d, &ta);
1683 16445 50         d->sign = (pstm_iszero(d) == PSTM_YES) ? PSTM_ZPOS : n;
1684             }
1685             LBL_ERR:
1686 16445           pstm_clear(&ta);
1687 16445           pstm_clear(&tb);
1688 16445           pstm_clear(&tq);
1689 16445           pstm_clear(&q);
1690 18587           return res;
1691             }
1692              
1693             # else /* defined PSTM_LARGE_DIV */
1694              
1695             /* This noticed is accompanied with this function to discourage its use. */
1696             # warning "This function has been noticed to give wrong results for some inputs."
1697              
1698             /******************************************************************************/
1699             /**
1700             c = a / b, d = remainder.
1701              
1702             @param[in] pool Memory pool
1703             @param[in] a Numerator
1704             @param[in] b Denominator
1705             @param[out] c The result
1706             @param[out] d If non-NULL, the remainder of the division is stored here
1707             @return < 0 on failure
1708              
1709             a/b => cb + d == a
1710             */
1711             int32_t pstm_div(psPool_t *pool, const pstm_int *a, const pstm_int *b,
1712             pstm_int *c, pstm_int *d)
1713             {
1714             pstm_int q, x, y, t1, t2;
1715             int32_t res;
1716             int16 n, t, i, norm, neg;
1717              
1718             /* is divisor zero ? */
1719             if (pstm_iszero(b) == 1)
1720             {
1721             return PS_LIMIT_FAIL;
1722             }
1723              
1724             /* if a < b then q=0, r = a */
1725             if (pstm_cmp_mag(a, b) == PSTM_LT)
1726             {
1727             if (d != NULL)
1728             {
1729             if (pstm_copy(a, d) != PSTM_OKAY)
1730             {
1731             return PS_MEM_FAIL;
1732             }
1733             }
1734             if (c != NULL)
1735             {
1736             pstm_zero(c);
1737             }
1738             return PSTM_OKAY;
1739             }
1740             /* Smart-size inits */
1741             if ((res = pstm_init_size(pool, &t1, a->alloc)) != PSTM_OKAY)
1742             {
1743             return res;
1744             }
1745             if ((res = pstm_init_size(pool, &t2, 3)) != PSTM_OKAY)
1746             {
1747             goto LBL_T1;
1748             }
1749             if ((res = pstm_init_copy(pool, &x, a, 0)) != PSTM_OKAY)
1750             {
1751             goto LBL_T2;
1752             }
1753             /* Used to be an init_copy on b but pstm_grow was always hit with triple
1754             size */
1755             if ((res = pstm_init_size(pool, &y, b->used * 3)) != PSTM_OKAY)
1756             {
1757             goto LBL_X;
1758             }
1759             if ((res = pstm_copy(b, &y)) != PSTM_OKAY)
1760             {
1761             goto LBL_Y;
1762             }
1763              
1764             /* fix the sign */
1765             neg = (a->sign == b->sign) ? PSTM_ZPOS : PSTM_NEG;
1766             x.sign = y.sign = PSTM_ZPOS;
1767              
1768             /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1769             norm = pstm_count_bits(&y) % DIGIT_BIT;
1770             if (norm < (int16_t) (DIGIT_BIT - 1))
1771             {
1772             norm = (DIGIT_BIT - 1) - norm;
1773             if ((res = pstm_mul_2d(&x, norm, &x)) != PSTM_OKAY)
1774             {
1775             goto LBL_Y;
1776             }
1777             if ((res = pstm_mul_2d(&y, norm, &y)) != PSTM_OKAY)
1778             {
1779             goto LBL_Y;
1780             }
1781             }
1782             else
1783             {
1784             norm = 0;
1785             }
1786              
1787             /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1788             n = x.used - 1;
1789             t = y.used - 1;
1790              
1791             if ((res = pstm_init_size(pool, &q, (n - t) + 1)) != PSTM_OKAY)
1792             {
1793             goto LBL_Y;
1794             }
1795             q.used = (n - t) + 1;
1796              
1797             /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1798             if ((res = pstm_lshd(&y, n - t)) != PSTM_OKAY) /* y = y*b**{n-t} */
1799             {
1800             goto LBL_Q;
1801             }
1802              
1803             while (pstm_cmp(&x, &y) != PSTM_LT)
1804             {
1805             ++(q.dp[n - t]);
1806             if ((res = pstm_sub(&x, &y, &x)) != PSTM_OKAY)
1807             {
1808             goto LBL_Q;
1809             }
1810             }
1811              
1812             /* reset y by shifting it back down */
1813             pstm_rshd(&y, n - t);
1814              
1815             /* step 3. for i from n down to (t + 1) */
1816             for (i = n; i >= (t + 1); i--)
1817             {
1818             if (i > x.used)
1819             {
1820             continue;
1821             }
1822              
1823             /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1824             * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1825             if (x.dp[i] == y.dp[t])
1826             {
1827             q.dp[i - t - 1] = (pstm_digit) ((((pstm_word) 1) << DIGIT_BIT) - 1);
1828             }
1829             else
1830             {
1831             pstm_word tmp;
1832             tmp = ((pstm_word) x.dp[i]) << ((pstm_word) DIGIT_BIT);
1833             tmp |= ((pstm_word) x.dp[i - 1]);
1834             # if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1835             psDiv64(&tmp, y.dp[t]);
1836             # elif defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1837             psDiv128(&tmp, y.dp[t]);
1838             # else
1839             tmp /= ((pstm_word) y.dp[t]);
1840             # endif /* USE_MATRIX_DIV64 */
1841             q.dp[i - t - 1] = (pstm_digit) (tmp);
1842             }
1843              
1844             /* while (q{i-t-1} * (yt * b + y{t-1})) >
1845             xi * b**2 + xi-1 * b + xi-2
1846              
1847             do q{i-t-1} -= 1;
1848             */
1849             q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
1850             do
1851             {
1852             q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
1853              
1854             /* find left hand */
1855             pstm_zero(&t1);
1856             t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1857             t1.dp[1] = y.dp[t];
1858             t1.used = 2;
1859             if ((res = pstm_mul_d(&t1, q.dp[i - t - 1], &t1)) != PSTM_OKAY)
1860             {
1861             goto LBL_Q;
1862             }
1863              
1864             /* find right hand */
1865             t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1866             t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1867             t2.dp[2] = x.dp[i];
1868             t2.used = 3;
1869             }
1870             while (pstm_cmp_mag(&t1, &t2) == PSTM_GT);
1871              
1872             /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1873             if ((res = pstm_mul_d(&y, q.dp[i - t - 1], &t1)) != PSTM_OKAY)
1874             {
1875             goto LBL_Q;
1876             }
1877              
1878             if ((res = pstm_lshd(&t1, i - t - 1)) != PSTM_OKAY)
1879             {
1880             goto LBL_Q;
1881             }
1882              
1883             if ((res = pstm_sub(&x, &t1, &x)) != PSTM_OKAY)
1884             {
1885             goto LBL_Q;
1886             }
1887              
1888             /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1889             if (x.sign == PSTM_NEG)
1890             {
1891             if ((res = pstm_copy(&y, &t1)) != PSTM_OKAY)
1892             {
1893             goto LBL_Q;
1894             }
1895             if ((res = pstm_lshd(&t1, i - t - 1)) != PSTM_OKAY)
1896             {
1897             goto LBL_Q;
1898             }
1899             if ((res = pstm_add(&x, &t1, &x)) != PSTM_OKAY)
1900             {
1901             goto LBL_Q;
1902             }
1903             q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
1904             }
1905             }
1906             /*
1907             now q is the quotient and x is the remainder (which we have to normalize)
1908             */
1909             /* get sign before writing to c */
1910             x.sign = x.used == 0 ? PSTM_ZPOS : a->sign;
1911              
1912             if (c != NULL)
1913             {
1914             pstm_clamp(&q);
1915             if (pstm_copy(&q, c) != PSTM_OKAY)
1916             {
1917             res = PS_MEM_FAIL;
1918             goto LBL_Q;
1919             }
1920             c->sign = neg;
1921             }
1922              
1923             if (d != NULL)
1924             {
1925             if ((res = pstm_div_2d(pool, &x, norm, &x, NULL)) != PSTM_OKAY)
1926             {
1927             goto LBL_Q;
1928             }
1929             /*
1930             the following is a kludge, essentially we were seeing the right
1931             remainder but with excess digits that should have been zero
1932             */
1933             for (i = b->used; i < x.used; i++)
1934             {
1935             x.dp[i] = 0;
1936             }
1937             pstm_clamp(&x);
1938             if (pstm_copy(&x, d) != PSTM_OKAY)
1939             {
1940             res = PS_MEM_FAIL;
1941             goto LBL_Q;
1942             }
1943             }
1944              
1945             res = PSTM_OKAY;
1946              
1947             LBL_Q: pstm_clear(&q);
1948             LBL_Y: pstm_clear(&y);
1949             LBL_X: pstm_clear(&x);
1950             LBL_T2: pstm_clear(&t2);
1951             LBL_T1: pstm_clear(&t1);
1952              
1953             return res;
1954             }
1955             # endif /* PSTM_LARGE_DIV */
1956              
1957             /******************************************************************************/
1958             /*
1959             Swap the elements of two integers, for cases where you can't simply swap
1960             the pstm_int pointers around
1961             */
1962 34406           void pstm_exch(pstm_int *a, pstm_int *b)
1963             {
1964             pstm_int t;
1965              
1966 34406           t = *a;
1967 34406           *a = *b;
1968 34406           *b = t;
1969 34406           }
1970              
1971             /******************************************************************************/
1972             /*
1973             c = a mod b, 0 <= c < b
1974             */
1975 18587           int32_t pstm_mod(psPool_t *pool, const pstm_int *a, const pstm_int *b, pstm_int *c)
1976             {
1977             pstm_int t;
1978             int32_t err;
1979              
1980             /* Smart-size */
1981 18587 50         if ((err = pstm_init_size(pool, &t, b->alloc)) != PSTM_OKAY)
1982             {
1983 0           return err;
1984             }
1985 18587 50         if ((err = pstm_div(pool, a, b, NULL, &t)) != PSTM_OKAY)
1986             {
1987 0           pstm_clear(&t);
1988 0           return err;
1989             }
1990 18587 100         if (t.sign != b->sign)
1991             {
1992 626           err = pstm_add(&t, b, c);
1993             }
1994             else
1995             {
1996 17961           pstm_exch(&t, c);
1997             }
1998 18587           pstm_clear(&t);
1999 18587           return err;
2000             }
2001              
2002             # if defined USE_MATRIX_RSA || defined USE_MATRIX_ECC || defined USE_MATRIX_DH
2003             /******************************************************************************/
2004             /*
2005             d = a * b (mod c)
2006             */
2007 7573           int32_t pstm_mulmod(psPool_t *pool, const pstm_int *a, const pstm_int *b,
2008             const pstm_int *c, pstm_int *d)
2009             {
2010             int32_t res;
2011             psSize_t size;
2012             pstm_int tmp;
2013              
2014             /*
2015             Smart-size pstm_inits. d is an output that is influenced by this local 't'
2016             so don't shrink 'd' if it wants to becuase this will lead to an pstm_grow
2017             in RSA operations
2018             */
2019 7573           size = a->used + b->used + 1;
2020 7573 100         if ((a == d) && (size < a->alloc))
    50          
2021             {
2022 1147           size = a->alloc;
2023             }
2024 7573 50         if ((res = pstm_init_size(pool, &tmp, size)) != PSTM_OKAY)
2025             {
2026 0           return res;
2027             }
2028 7573 50         if ((res = pstm_mul_comba(pool, a, b, &tmp, NULL, 0)) != PSTM_OKAY)
2029             {
2030 0           pstm_clear(&tmp);
2031 0           return res;
2032             }
2033 7573           res = pstm_mod(pool, &tmp, c, d);
2034 7573           pstm_clear(&tmp);
2035 7573           return res;
2036             }
2037              
2038             /******************************************************************************/
2039             /*
2040             * y = g**x (mod p)
2041             * Some restrictions...
2042             * x must be positive and < p
2043             * p must be positive, odd, and [512,1024,1536,2048,3072,4096] bits
2044             */
2045 5648           int32_t pstm_exptmod(psPool_t *pool, const pstm_int *G, const pstm_int *X,
2046             const pstm_int *P, pstm_int *Y)
2047             {
2048             pstm_int M[32], res; /* Keep this winsize based: (1 << max_winsize) */
2049             pstm_digit buf, mp;
2050             pstm_digit *paD;
2051             int32 err, bitbuf;
2052             int16 bitcpy, bitcnt, mode, digidx, x, y, winsize;
2053             uint32 paDlen;
2054              
2055 5648           x = pstm_count_bits(P);
2056 5648 50         switch (x)
2057             {
2058             case 512:
2059             case 1024:
2060             case 1536:
2061             case 2048:
2062             case 3072:
2063             case 4096:
2064 5648           break;
2065             default:
2066             psTraceIntCrypto("pstm_exptmod prime size failed: %hu\n", x);
2067 0           return -1;
2068             }
2069             # ifdef USE_CONSTANT_TIME_MODEXP
2070 5648 50         if (P->dp[0] & 1)
2071             {
2072             pstmnt_word Base[512 / sizeof(pstmnt_word)];
2073             pstmnt_word Mod[512 / sizeof(pstmnt_word)];
2074             pstmnt_word Temp[512 / sizeof(pstmnt_word) * 7 + 1];
2075 5648           pstmnt_word mp = pstmnt_neg_small_inv(pstmnt_const_ptr(P));
2076              
2077              
2078 5648           memset(Base, 0, sizeof(Base));
2079 5648 100         if (G->used <= P->used)
2080             {
2081 3354           memcpy(Base, pstmnt_const_ptr(G), pstmnt_size_bytes(G));
2082             }
2083             else
2084             {
2085             /* Base > P -> have to compute Base % P to get the actual base. */
2086             int32_t err;
2087             pstm_int tmp_int;
2088 2294 50         if ((err = pstm_init_size(pool, &tmp_int, P->used)) != PSTM_OKAY)
2089             {
2090 0           return err;
2091             }
2092 2294           err = pstm_mod(pool, G, P, &tmp_int);
2093 2294           memcpy(Base, pstmnt_const_ptr(&tmp_int),
2094 2294           pstmnt_size_bytes(&tmp_int));
2095 2294           pstm_clear(&tmp_int);
2096 2294 50         if (err != PSTM_OKAY)
2097             {
2098 2294           return err;
2099             }
2100             }
2101 5648           memcpy(Mod, pstmnt_const_ptr(P), pstmnt_size_bytes(P));
2102              
2103 5648           Y->used = P->used;
2104 5648 50         if (Y->used > Y->alloc)
2105             {
2106 0 0         if (pstm_grow(Y, Y->used) != PSTM_OKAY)
2107             {
2108 0           return PS_MEM_FAIL;
2109             }
2110             }
2111              
2112             /* Use constant time variant. */
2113 5648           pstmnt_montgomery_input(Base, Mod, Temp,
2114             pstmnt_ptr(Y), pstmnt_size(P), mp);
2115 11296           pstmnt_mod_exp_montgomery_skip(pstmnt_const_ptr(Y), pstmnt_const_ptr(X),
2116             pstmnt_ptr(Y), 0,
2117 5648           pstm_count_bits(X), Mod, Temp,
2118             mp, pstmnt_size(P));
2119 5648           pstmnt_montgomery_output(pstmnt_const_ptr(Y), pstmnt_ptr(Y), Mod, Temp,
2120             pstmnt_size(P), mp);
2121 5648           pstm_clamp(Y);
2122 5648           memset_s(Base, sizeof(Base), 0, sizeof(Base));
2123 5648           memset_s(Mod, sizeof(Mod), 0, sizeof(Mod));
2124 5648           memset_s(Temp, sizeof(Temp), 0, sizeof(Temp));
2125 5648           return PSTM_OKAY;
2126             }
2127             # endif /* USE_CONSTANT_TIME_MODEXP */
2128             /* set window size from what user set as optimization */
2129 0           x = pstm_count_bits(X);
2130 0 0         if (x < 50)
2131             {
2132 0           winsize = 2;
2133             }
2134             else
2135             {
2136 0           winsize = PS_EXPTMOD_WINSIZE;
2137             }
2138              
2139             /* now setup montgomery */
2140 0 0         if ((err = pstm_montgomery_setup(P, &mp)) != PSTM_OKAY)
2141             {
2142 0           return err;
2143             }
2144              
2145             /* setup result */
2146 0 0         if ((err = pstm_init_size(pool, &res, (P->used * 2) + 1)) != PSTM_OKAY)
2147             {
2148 0           return err;
2149             }
2150             /*
2151             create M table
2152             The M table contains powers of the input base, e.g. M[x] = G^x mod P
2153             The first half of the table is not computed though except for M[0] and M[1]
2154             */
2155             /* now we need R mod m */
2156 0 0         if ((err = pstm_montgomery_calc_normalization(&res, P)) != PSTM_OKAY)
2157             {
2158 0           goto LBL_RES;
2159             }
2160             /*
2161             init M array
2162             init first cell
2163             */
2164 0 0         if ((err = pstm_init_size(pool, &M[1], res.used)) != PSTM_OKAY)
2165             {
2166 0           goto LBL_RES;
2167             }
2168              
2169             /* now set M[1] to G * R mod m */
2170 0 0         if (pstm_cmp_mag(P, G) != PSTM_GT)
2171             {
2172             /* G > P so we reduce it first */
2173 0 0         if ((err = pstm_mod(pool, G, P, &M[1])) != PSTM_OKAY)
2174             {
2175 0           goto LBL_M;
2176             }
2177             }
2178             else
2179             {
2180 0 0         if ((err = pstm_copy(G, &M[1])) != PSTM_OKAY)
2181             {
2182 0           goto LBL_M;
2183             }
2184             }
2185 0 0         if ((err = pstm_mulmod(pool, &M[1], &res, P, &M[1])) != PSTM_OKAY)
2186             {
2187 0           goto LBL_M;
2188             }
2189             /* Pre-allocated digit. Used for mul, sqr, AND reduce */
2190 0           paDlen = ((M[1].used + 3) * 2) * sizeof(pstm_digit);
2191 0 0         if ((paD = psMalloc(pool, paDlen)) == NULL)
2192             {
2193 0           err = PS_MEM_FAIL;
2194 0           goto LBL_M;
2195             }
2196             /* compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times */
2197 0 0         if (pstm_init_copy(pool, &M[1 << (winsize - 1)], &M[1], 1) != PSTM_OKAY)
2198             {
2199 0           err = PS_MEM_FAIL;
2200 0           goto LBL_PAD;
2201             }
2202 0 0         for (x = 0; x < (winsize - 1); x++)
2203             {
2204 0 0         if ((err = pstm_sqr_comba(pool, &M[1 << (winsize - 1)],
2205 0           &M[1 << (winsize - 1)], paD, paDlen)) != PSTM_OKAY)
2206             {
2207 0           goto LBL_PAD;
2208             }
2209 0 0         if ((err = pstm_montgomery_reduce(pool, &M[1 << (winsize - 1)], P, mp,
2210             paD, paDlen)) != PSTM_OKAY)
2211             {
2212 0           goto LBL_PAD;
2213             }
2214             }
2215             /* now init the second half of the array */
2216 0 0         for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++)
2217             {
2218 0 0         if ((err = pstm_init_size(pool, &M[x], M[1 << (winsize - 1)].alloc + 1))
2219             != PSTM_OKAY)
2220             {
2221 0 0         for (y = 1 << (winsize - 1); y < x; y++)
2222             {
2223 0           pstm_clear(&M[y]);
2224             }
2225 0           goto LBL_PAD;
2226             }
2227             }
2228              
2229             /* create upper table */
2230 0 0         for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++)
2231             {
2232 0 0         if ((err = pstm_mul_comba(pool, &M[x - 1], &M[1], &M[x], paD, paDlen))
2233             != PSTM_OKAY)
2234             {
2235 0           goto LBL_MARRAY;
2236             }
2237 0 0         if ((err = pstm_montgomery_reduce(pool, &M[x], P, mp, paD, paDlen)) !=
2238             PSTM_OKAY)
2239             {
2240 0           goto LBL_MARRAY;
2241             }
2242             }
2243              
2244             /* set initial mode and bit cnt */
2245 0           mode = 0;
2246 0           bitcnt = 1;
2247 0           buf = 0;
2248 0           digidx = X->used - 1;
2249 0           bitcpy = 0;
2250 0           bitbuf = 0;
2251              
2252             for (;; )
2253             {
2254             /* grab next digit as required */
2255 0 0         if (--bitcnt == 0)
2256             {
2257             /* if digidx == -1 we are out of digits so break */
2258 0 0         if (digidx == -1)
2259             {
2260 0           break;
2261             }
2262             /* read next digit and reset bitcnt */
2263 0           buf = X->dp[digidx--];
2264 0           bitcnt = (int32) DIGIT_BIT;
2265             }
2266              
2267             /* grab the next msb from the exponent */
2268 0           y = (pstm_digit) (buf >> (DIGIT_BIT - 1)) & 1;
2269 0           buf <<= (pstm_digit) 1;
2270             /*
2271             If the bit is zero and mode == 0 then we ignore it.
2272             These represent the leading zero bits before the first 1 bit
2273             in the exponent. Technically this opt is not required but it
2274             does lower the # of trivial squaring/reductions used
2275             */
2276 0 0         if (mode == 0 && y == 0)
    0          
2277             {
2278 0           continue;
2279             }
2280              
2281             /* if the bit is zero and mode == 1 then we square */
2282 0 0         if (mode == 1 && y == 0)
    0          
2283             {
2284 0 0         if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
2285             PSTM_OKAY)
2286             {
2287 0           goto LBL_MARRAY;
2288             }
2289 0 0         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
2290             != PSTM_OKAY)
2291             {
2292 0           goto LBL_MARRAY;
2293             }
2294 0           continue;
2295             }
2296              
2297             /* else we add it to the window */
2298 0           bitbuf |= (y << (winsize - ++bitcpy));
2299 0           mode = 2;
2300              
2301 0 0         if (bitcpy == winsize)
2302             {
2303             /* ok window is filled so square as required and mul square first */
2304 0 0         for (x = 0; x < winsize; x++)
2305             {
2306 0 0         if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
2307             PSTM_OKAY)
2308             {
2309 0           goto LBL_MARRAY;
2310             }
2311 0 0         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
2312             paDlen)) != PSTM_OKAY)
2313             {
2314 0           goto LBL_MARRAY;
2315             }
2316             }
2317              
2318             /* then multiply */
2319 0 0         if ((err = pstm_mul_comba(pool, &res, &M[bitbuf], &res, paD,
2320             paDlen)) != PSTM_OKAY)
2321             {
2322 0           goto LBL_MARRAY;
2323             }
2324 0 0         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
2325             != PSTM_OKAY)
2326             {
2327 0           goto LBL_MARRAY;
2328             }
2329              
2330             /* empty window and reset */
2331 0           bitcpy = 0;
2332 0           bitbuf = 0;
2333 0           mode = 1;
2334             }
2335 0           }
2336              
2337             /* if bits remain then square/multiply */
2338 0 0         if (mode == 2 && bitcpy > 0)
    0          
2339             {
2340             /* square then multiply if the bit is set */
2341 0 0         for (x = 0; x < bitcpy; x++)
2342             {
2343 0 0         if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
2344             PSTM_OKAY)
2345             {
2346 0           goto LBL_MARRAY;
2347             }
2348 0 0         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
2349             != PSTM_OKAY)
2350             {
2351 0           goto LBL_MARRAY;
2352             }
2353              
2354             /* get next bit of the window */
2355 0           bitbuf <<= 1;
2356 0 0         if ((bitbuf & (1 << winsize)) != 0)
2357             {
2358             /* then multiply */
2359 0 0         if ((err = pstm_mul_comba(pool, &res, &M[1], &res, paD, paDlen))
2360             != PSTM_OKAY)
2361             {
2362 0           goto LBL_MARRAY;
2363             }
2364 0 0         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
2365             paDlen)) != PSTM_OKAY)
2366             {
2367 0           goto LBL_MARRAY;
2368             }
2369             }
2370             }
2371             }
2372             /*
2373             Fix up result if Montgomery reduction is used recall that any value in a
2374             Montgomery system is actually multiplied by R mod n. So we have to reduce
2375             one more time to cancel out the factor of R.
2376             */
2377 0 0         if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen)) !=
2378             PSTM_OKAY)
2379             {
2380 0           goto LBL_MARRAY;
2381             }
2382             /* swap res with Y */
2383 0 0         if ((err = pstm_copy(&res, Y)) != PSTM_OKAY)
2384             {
2385 0           goto LBL_MARRAY;
2386             }
2387 0           err = PSTM_OKAY;
2388             LBL_MARRAY:
2389 0 0         for (x = 1 << (winsize - 1); x < (1 << winsize); x++)
2390             {
2391 0           pstm_clear(&M[x]);
2392             }
2393 0           LBL_PAD: psFree(paD, pool);
2394 0           LBL_M: pstm_clear(&M[1]);
2395 0           LBL_RES: pstm_clear(&res);
2396              
2397 5648           return err;
2398             }
2399             # endif /* USE_MATRIX_RSA || USE_MATRIX_ECC || USE_MATRIX_DH */
2400              
2401             /******************************************************************************/
2402             /**
2403             c = a + b.
2404             */
2405 15736342           int32_t pstm_add(const pstm_int *a, const pstm_int *b, pstm_int *c)
2406             {
2407             int32_t res;
2408             uint8_t sa, sb;
2409              
2410             /* get sign of both inputs */
2411 15736342           sa = a->sign;
2412 15736342           sb = b->sign;
2413              
2414             /* handle two cases, not four */
2415 15736342 100         if (sa == sb)
2416             {
2417             /* both positive or both negative, add their mags, copy the sign */
2418 12236752           c->sign = sa;
2419 12236752 50         if ((res = s_pstm_add(a, b, c)) != PSTM_OKAY)
2420             {
2421 0           return res;
2422             }
2423             }
2424             else
2425             {
2426             /*
2427             one positive, the other negative
2428             subtract the one with the greater magnitude from the one of the lesser
2429             magnitude. The result gets the sign of the one with the greater mag.
2430             */
2431 3499590 100         if (pstm_cmp_mag(a, b) == PSTM_LT)
2432             {
2433 3499353           c->sign = sb;
2434 3499353 50         if ((res = pstm_sub_s(b, a, c)) != PSTM_OKAY)
2435             {
2436 0           return res;
2437             }
2438             }
2439             else
2440             {
2441 237           c->sign = sa;
2442 237 50         if ((res = pstm_sub_s(a, b, c)) != PSTM_OKAY)
2443             {
2444 0           return res;
2445             }
2446             }
2447             }
2448 15736342           return PS_SUCCESS;
2449             }
2450              
2451             /******************************************************************************/
2452             /*
2453             No reverse. Useful in some of the EIP-154 PKA stuff where special byte
2454             order seems to come into play more often
2455             */
2456 0           int32_t pstm_to_unsigned_bin_nr(psPool_t *pool, const pstm_int *a, unsigned char *b)
2457             {
2458             int32_t res;
2459             uint16_t x;
2460 0           pstm_int t = { 0 };
2461              
2462 0 0         if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY)
2463             {
2464 0           return res;
2465             }
2466              
2467 0           x = 0;
2468 0 0         while (pstm_iszero(&t) == 0)
2469             {
2470 0           b[x++] = (unsigned char) (t.dp[0] & 255);
2471 0 0         if ((res = pstm_div_2d(pool, &t, 8, &t, NULL)) != PSTM_OKAY)
2472             {
2473 0           pstm_clear(&t);
2474 0           return res;
2475             }
2476             }
2477 0           pstm_clear(&t);
2478 0           return PS_SUCCESS;
2479             }
2480              
2481             /******************************************************************************/
2482             /*
2483             reverse an array, used for unsigned bin code
2484             */
2485 11040           static void pstm_reverse(unsigned char *s, psSize_t len)
2486             {
2487             uint16_t ix, iy;
2488             unsigned char t;
2489              
2490 11040 50         if (len == 0)
2491             {
2492 0           return;
2493             }
2494 11040           ix = 0;
2495 11040           iy = len - 1;
2496 795547 100         while (ix < iy)
2497             {
2498 784507           t = s[ix];
2499 784507           s[ix] = s[iy];
2500 784507           s[iy] = t;
2501 784507           ++ix;
2502 784507           --iy;
2503             }
2504             }
2505              
2506             /******************************************************************************/
2507             /**
2508             Write a pstm format integer to a raw binary format.
2509             @return < 0 on failure. PS_SUCCESS on success.
2510             */
2511 11040           int32_t pstm_to_unsigned_bin(psPool_t *pool, const pstm_int *a, unsigned char *b)
2512              
2513             {
2514             int32_t res;
2515             uint16_t x;
2516 11040           pstm_int t = { 0 };
2517              
2518 11040 50         if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY)
2519             {
2520 0           return res;
2521             }
2522 11040           x = 0;
2523 1587391 100         while (pstm_iszero(&t) == 0)
2524             {
2525 1576351           b[x++] = (unsigned char) (t.dp[0] & 255);
2526 1576351 50         if ((res = pstm_div_2d(pool, &t, 8, &t, NULL)) != PSTM_OKAY)
2527             {
2528 0           pstm_clear(&t);
2529 0           return res;
2530             }
2531             }
2532 11040           pstm_reverse(b, x);
2533 11040           pstm_clear(&t);
2534 11040           return PS_SUCCESS;
2535             }
2536              
2537             /******************************************************************************/
2538             /**
2539             c = 1/a (mod b).
2540              
2541             @note Slow version supporting an even 'b'. Should call pstm_invmod() and let it
2542             decide which to use.
2543              
2544             Need invmod for ECC and also private key loading for hardware crypto
2545             in cases where dQ > dP. The values must be switched and a new qP must be
2546             calculated using this function
2547             */
2548 0           static int32_t pstm_invmod_slow(psPool_t *pool, const pstm_int *a,
2549             const pstm_int *b, pstm_int *c)
2550             {
2551             pstm_int x, y, u, v, A, B, C, D;
2552             int32 res;
2553              
2554             /* b cannot be negative */
2555 0 0         if (b->sign == PSTM_NEG || pstm_iszero(b) == 1)
    0          
2556             {
2557 0           return PS_LIMIT_FAIL;
2558             }
2559              
2560             /* init temps */
2561 0 0         if (pstm_init_size(pool, &x, b->used) != PSTM_OKAY)
2562             {
2563 0           return PS_MEM_FAIL;
2564             }
2565              
2566             /* x = a, y = b */
2567 0 0         if ((res = pstm_mod(pool, a, b, &x)) != PSTM_OKAY)
2568             {
2569 0           goto LBL_X;
2570             }
2571              
2572 0 0         if (pstm_init_copy(pool, &y, b, 0) != PSTM_OKAY)
2573             {
2574 0           goto LBL_X;
2575             }
2576              
2577             /* 2. [modified] if x,y are both even then return an error! */
2578 0 0         if (pstm_iseven(&x) == 1 && pstm_iseven(&y) == 1)
    0          
    0          
    0          
    0          
    0          
2579             {
2580 0           res = PS_FAILURE;
2581 0           goto LBL_Y;
2582             }
2583              
2584             /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2585 0 0         if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY)
2586             {
2587 0           goto LBL_Y;
2588             }
2589 0 0         if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY)
2590             {
2591 0           goto LBL_U;
2592             }
2593              
2594 0 0         if ((res = pstm_init_size(pool, &A, sizeof(pstm_digit))) != PSTM_OKAY)
2595             {
2596 0           goto LBL_V;
2597             }
2598              
2599 0 0         if ((res = pstm_init_size(pool, &D, sizeof(pstm_digit))) != PSTM_OKAY)
2600             {
2601 0           goto LBL_A;
2602             }
2603 0           pstm_set(&A, 1);
2604 0           pstm_set(&D, 1);
2605              
2606 0 0         if ((res = pstm_init(pool, &B)) != PSTM_OKAY)
2607             {
2608 0           goto LBL_D;
2609             }
2610 0 0         if ((res = pstm_init(pool, &C)) != PSTM_OKAY)
2611             {
2612 0           goto LBL_B;
2613             }
2614              
2615             top:
2616             /* 4. while u is even do */
2617 0 0         while (pstm_iseven(&u) == 1)
    0          
    0          
2618             {
2619             /* 4.1 u = u/2 */
2620 0 0         if ((res = pstm_div_2(&u, &u)) != PSTM_OKAY)
2621             {
2622 0           goto LBL_C;
2623             }
2624              
2625             /* 4.2 if A or B is odd then */
2626 0 0         if (pstm_isodd(&A) == 1 || pstm_isodd(&B) == 1)
    0          
    0          
    0          
    0          
    0          
2627             {
2628             /* A = (A+y)/2, B = (B-x)/2 */
2629 0 0         if ((res = pstm_add(&A, &y, &A)) != PSTM_OKAY)
2630             {
2631 0           goto LBL_C;
2632             }
2633 0 0         if ((res = pstm_sub(&B, &x, &B)) != PSTM_OKAY)
2634             {
2635 0           goto LBL_C;
2636             }
2637             }
2638             /* A = A/2, B = B/2 */
2639 0 0         if ((res = pstm_div_2(&A, &A)) != PSTM_OKAY)
2640             {
2641 0           goto LBL_C;
2642             }
2643 0 0         if ((res = pstm_div_2(&B, &B)) != PSTM_OKAY)
2644             {
2645 0           goto LBL_C;
2646             }
2647             }
2648              
2649             /* 5. while v is even do */
2650 0 0         while (pstm_iseven(&v) == 1)
    0          
    0          
2651             {
2652             /* 5.1 v = v/2 */
2653 0 0         if ((res = pstm_div_2(&v, &v)) != PSTM_OKAY)
2654             {
2655 0           goto LBL_C;
2656             }
2657              
2658             /* 5.2 if C or D is odd then */
2659 0 0         if (pstm_isodd(&C) == 1 || pstm_isodd(&D) == 1)
    0          
    0          
    0          
    0          
    0          
2660             {
2661             /* C = (C+y)/2, D = (D-x)/2 */
2662 0 0         if ((res = pstm_add(&C, &y, &C)) != PSTM_OKAY)
2663             {
2664 0           goto LBL_C;
2665             }
2666 0 0         if ((res = pstm_sub(&D, &x, &D)) != PSTM_OKAY)
2667             {
2668 0           goto LBL_C;
2669             }
2670             }
2671             /* C = C/2, D = D/2 */
2672 0 0         if ((res = pstm_div_2(&C, &C)) != PSTM_OKAY)
2673             {
2674 0           goto LBL_C;
2675             }
2676 0 0         if ((res = pstm_div_2(&D, &D)) != PSTM_OKAY)
2677             {
2678 0           goto LBL_C;
2679             }
2680             }
2681              
2682             /* 6. if u >= v then */
2683 0 0         if (pstm_cmp(&u, &v) != PSTM_LT)
2684             {
2685             /* u = u - v, A = A - C, B = B - D */
2686 0 0         if ((res = pstm_sub(&u, &v, &u)) != PSTM_OKAY)
2687             {
2688 0           goto LBL_C;
2689             }
2690 0 0         if ((res = pstm_sub(&A, &C, &A)) != PSTM_OKAY)
2691             {
2692 0           goto LBL_C;
2693             }
2694 0 0         if ((res = pstm_sub(&B, &D, &B)) != PSTM_OKAY)
2695             {
2696 0           goto LBL_C;
2697             }
2698             }
2699             else
2700             {
2701             /* v - v - u, C = C - A, D = D - B */
2702 0 0         if ((res = pstm_sub(&v, &u, &v)) != PSTM_OKAY)
2703             {
2704 0           goto LBL_C;
2705             }
2706 0 0         if ((res = pstm_sub(&C, &A, &C)) != PSTM_OKAY)
2707             {
2708 0           goto LBL_C;
2709             }
2710 0 0         if ((res = pstm_sub(&D, &B, &D)) != PSTM_OKAY)
2711             {
2712 0           goto LBL_C;
2713             }
2714             }
2715              
2716             /* if not zero goto step 4 */
2717 0 0         if (pstm_iszero(&u) == 0)
2718             {
2719 0           goto top;
2720             }
2721              
2722             /* now a = C, b = D, gcd == g*v */
2723              
2724             /* if v != 1 then there is no inverse */
2725 0 0         if (pstm_cmp_d(&v, 1) != PSTM_EQ)
2726             {
2727 0           res = PS_FAILURE;
2728 0           goto LBL_C;
2729             }
2730              
2731             /* if its too low */
2732 0 0         while (pstm_cmp_d(&C, 0) == PSTM_LT)
2733             {
2734 0 0         if ((res = pstm_add(&C, b, &C)) != PSTM_OKAY)
2735             {
2736 0           goto LBL_C;
2737             }
2738             }
2739              
2740             /* too big */
2741 0 0         while (pstm_cmp_mag(&C, b) != PSTM_LT)
2742             {
2743 0 0         if ((res = pstm_sub(&C, b, &C)) != PSTM_OKAY)
2744             {
2745 0           goto LBL_C;
2746             }
2747             }
2748              
2749             /* C is now the inverse */
2750 0 0         if ((res = pstm_copy(&C, c)) != PSTM_OKAY)
2751             {
2752 0           goto LBL_C;
2753             }
2754 0           res = PSTM_OKAY;
2755              
2756 0           LBL_C: pstm_clear(&C);
2757 0           LBL_D: pstm_clear(&D);
2758 0           LBL_B: pstm_clear(&B);
2759 0           LBL_A: pstm_clear(&A);
2760 0           LBL_V: pstm_clear(&v);
2761 0           LBL_U: pstm_clear(&u);
2762 0           LBL_Y: pstm_clear(&y);
2763 0           LBL_X: pstm_clear(&x);
2764              
2765 0           return res;
2766             }
2767              
2768             /**
2769             c = 1/a (mod b).
2770             This code is for for odd 'b'. pstm_invmod_slow() will be called if 'b'
2771             is even.
2772             */
2773 2142           int32_t pstm_invmod(psPool_t *pool, const pstm_int *a, const pstm_int *b, pstm_int *c)
2774             {
2775             pstm_int x, y, u, v, B, D;
2776             int32 res;
2777             uint16 neg, sanity;
2778              
2779             /* 2. [modified] b must be odd */
2780 2142 50         if (pstm_iseven(b) == 1)
    50          
    50          
2781             {
2782 0           return pstm_invmod_slow(pool, a, b, c);
2783             }
2784              
2785             /* x == modulus, y == value to invert */
2786 2142 50         if ((res = pstm_init_copy(pool, &x, b, 0)) != PSTM_OKAY)
2787             {
2788 0           return res;
2789             }
2790              
2791 2142 50         if ((res = pstm_init_size(pool, &y, a->alloc)) != PSTM_OKAY)
2792             {
2793 0           goto LBL_X;
2794             }
2795              
2796             /* we need y = |a| */
2797 2142 50         if ((res = pstm_abs(a, &y)) != PSTM_OKAY)
2798             {
2799 0           goto LBL_X;
2800             }
2801              
2802             /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2803 2142 50         if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY)
2804             {
2805 0           goto LBL_Y;
2806             }
2807 2142 50         if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY)
2808             {
2809 0           goto LBL_U;
2810             }
2811 2142 50         if ((res = pstm_init(pool, &B)) != PSTM_OKAY)
2812             {
2813 0           goto LBL_V;
2814             }
2815 2142 50         if ((res = pstm_init(pool, &D)) != PSTM_OKAY)
2816             {
2817 0           goto LBL_B;
2818             }
2819              
2820 2142           pstm_set(&D, 1);
2821              
2822 2142           sanity = 0;
2823             top:
2824             /* 4. while u is even do */
2825 1580818 50         while (pstm_iseven(&u) == 1)
    100          
    100          
2826             {
2827             /* 4.1 u = u/2 */
2828 786098 50         if ((res = pstm_div_2(&u, &u)) != PSTM_OKAY)
2829             {
2830 0           goto LBL_D;
2831             }
2832              
2833             /* 4.2 if B is odd then */
2834 786098 50         if (pstm_isodd(&B) == 1)
    100          
    100          
2835             {
2836 422198 50         if ((res = pstm_sub(&B, &x, &B)) != PSTM_OKAY)
2837             {
2838 0           goto LBL_D;
2839             }
2840             }
2841             /* B = B/2 */
2842 786098 50         if ((res = pstm_div_2(&B, &B)) != PSTM_OKAY)
2843             {
2844 0           goto LBL_D;
2845             }
2846             }
2847              
2848             /* 5. while v is even do */
2849 1581745 50         while (pstm_iseven(&v) == 1)
    100          
    100          
2850             {
2851             /* 5.1 v = v/2 */
2852 787025 50         if ((res = pstm_div_2(&v, &v)) != PSTM_OKAY)
2853             {
2854 0           goto LBL_D;
2855             }
2856             /* 5.2 if D is odd then */
2857 787025 50         if (pstm_isodd(&D) == 1)
    100          
    100          
2858             {
2859             /* D = (D-x)/2 */
2860 681074 50         if ((res = pstm_sub(&D, &x, &D)) != PSTM_OKAY)
2861             {
2862 0           goto LBL_D;
2863             }
2864             }
2865             /* D = D/2 */
2866 787025 50         if ((res = pstm_div_2(&D, &D)) != PSTM_OKAY)
2867             {
2868 0           goto LBL_D;
2869             }
2870             }
2871              
2872             /* 6. if u >= v then */
2873 794720 100         if (pstm_cmp(&u, &v) != PSTM_LT)
2874             {
2875             /* u = u - v, B = B - D */
2876 402891 50         if ((res = pstm_sub(&u, &v, &u)) != PSTM_OKAY)
2877             {
2878 0           goto LBL_D;
2879             }
2880 402891 50         if ((res = pstm_sub(&B, &D, &B)) != PSTM_OKAY)
2881             {
2882 0           goto LBL_D;
2883             }
2884             }
2885             else
2886             {
2887             /* v - v - u, D = D - B */
2888 391829 50         if ((res = pstm_sub(&v, &u, &v)) != PSTM_OKAY)
2889             {
2890 0           goto LBL_D;
2891             }
2892 391829 50         if ((res = pstm_sub(&D, &B, &D)) != PSTM_OKAY)
2893             {
2894 0           goto LBL_D;
2895             }
2896             }
2897              
2898             /* if not zero goto step 4 */
2899 794720 50         if (sanity++ > 4096)
2900             {
2901 0           res = PS_LIMIT_FAIL;
2902 0           goto LBL_D;
2903             }
2904 794720 100         if (pstm_iszero(&u) == 0)
2905             {
2906 792578           goto top;
2907             }
2908              
2909             /* now a = C, b = D, gcd == g*v */
2910              
2911             /* if v != 1 then there is no inverse */
2912 2142 50         if (pstm_cmp_d(&v, 1) != PSTM_EQ)
2913             {
2914 0           res = PS_FAILURE;
2915 0           goto LBL_D;
2916             }
2917              
2918             /* b is now the inverse */
2919 2142           neg = a->sign;
2920 3363 100         while (D.sign == PSTM_NEG)
2921             {
2922 1221 50         if ((res = pstm_add(&D, b, &D)) != PSTM_OKAY)
2923             {
2924 0           goto LBL_D;
2925             }
2926             }
2927 2142 50         if ((res = pstm_copy(&D, c)) != PSTM_OKAY)
2928             {
2929 0           goto LBL_D;
2930             }
2931 2142           c->sign = neg;
2932 2142           res = PSTM_OKAY;
2933              
2934 2142           LBL_D: pstm_clear(&D);
2935 2142           LBL_B: pstm_clear(&B);
2936 2142           LBL_V: pstm_clear(&v);
2937 2142           LBL_U: pstm_clear(&u);
2938 2142           LBL_Y: pstm_clear(&y);
2939 2142           LBL_X: pstm_clear(&x);
2940 2142           return res;
2941             }
2942              
2943             /******************************************************************************/
2944              
2945             #endif /* USE_MATRIX_RSA || USE_MATRIX_ECC || USE_MATRIX_DH || USE_CL_RSA || USE_CL_DH || USE_QUICK_ASSIST_RSA || USE_QUICK_ASSIST_ECC */
2946              
2947             /******************************************************************************/
2948