2 * Copyright (C) 2017 Denys Vlasenko
4 * Licensed under GPLv2, see file LICENSE in this source tree.
8 /* The file is taken almost verbatim from matrixssl-3-7-2b-open/crypto/math/.
9 * Changes are flagged with //bbox
14 * @version 33ef80f (HEAD, tag: MATRIXSSL-3-7-2-OPEN, tag: MATRIXSSL-3-7-2-COMM, origin/master, origin/HEAD, master)
16 * Multiprecision number implementation.
19 * Copyright (c) 2013-2015 INSIDE Secure Corporation
20 * Copyright (c) PeerSec Networks, 2002-2011
23 * The latest version of this code is available at http://www.matrixssl.org
25 * This software is open source; you can redistribute it and/or modify
26 * it under the terms of the GNU General Public License as published by
27 * the Free Software Foundation; either version 2 of the License, or
28 * (at your option) any later version.
30 * This General Public License does NOT permit incorporating this software
31 * into proprietary programs. If you are unable to comply with the GPL, a
32 * commercial license for this software may be purchased from INSIDE at
33 * http://www.insidesecure.com/eng/Company/Locations
35 * This program is distributed in WITHOUT ANY WARRANTY; without even the
36 * implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
37 * See the GNU General Public License for more details.
39 * You should have received a copy of the GNU General Public License
40 * along with this program; if not, write to the Free Software
41 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
42 * http://www.gnu.org/copyleft/gpl.html
44 /******************************************************************************/
47 //#include "../cryptoApi.h"
51 static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c); //bbox: was int16 b
52 #define pstm_mul_2d(a, b, c) (pstm_mul_2d(a, b, c), PSTM_OKAY)
54 /******************************************************************************/
56 init an pstm_int for a given size
59 #define pstm_init_size(pool, a, size) \
60 pstm_init_size( a, size)
61 int32 FAST_FUNC pstm_init_size(psPool_t *pool, pstm_int * a, uint32 size)
69 a->dp = xzalloc(sizeof (pstm_digit) * size);//bbox
70 //bbox a->pool = pool;
78 // for (x = 0; x < size; x++) {
84 #define pstm_init_size(pool, a, size) (pstm_init_size(a, size), PSTM_OKAY)
86 /******************************************************************************/
91 #define pstm_init(pool, a) \
93 static int32 pstm_init(psPool_t *pool, pstm_int * a)
98 allocate memory required and clear it
100 a->dp = xzalloc(sizeof (pstm_digit) * PSTM_DEFAULT_INIT);//bbox
102 set the digits to zero
105 // for (i = 0; i < PSTM_DEFAULT_INIT; i++) {
109 set the used to zero, allocated digits to the default precision and sign
112 //bbox a->pool = pool;
114 a->alloc = PSTM_DEFAULT_INIT;
120 #define pstm_init(pool, a) (pstm_init(a), PSTM_OKAY)
122 /******************************************************************************/
127 int32 FAST_FUNC pstm_grow(pstm_int * a, int size)
129 int i; //bbox: was int16
133 If the alloc size is smaller alloc more ram.
135 if (a->alloc < size) {
137 Reallocate the array a->dp
139 We store the return in a temporary variable in case the operation
140 failed we don't want to overwrite the dp member of a.
142 tmp = xrealloc(a->dp, sizeof (pstm_digit) * size);//bbox
144 reallocation succeeded so set a->dp
152 for (; i < a->alloc; i++) {
158 #define pstm_grow(a, size) (pstm_grow(a, size), PSTM_OKAY)
160 /******************************************************************************/
162 copy, b = a (b must be pre-allocated)
165 int32 pstm_copy(pstm_int * a, pstm_int * b)
170 If dst == src do nothing
178 if (b->alloc < a->used) {
179 if ((res = pstm_grow (b, a->used)) != PSTM_OKAY) {
184 Zero b and copy the parameters over
187 register pstm_digit *tmpa, *tmpb;
189 /* pointer aliases */
196 /* copy all the digits */
197 for (n = 0; n < a->used; n++) {
201 /* clear high digits */
202 for (; n < b->used; n++) {
207 copy used count and sign
213 #define pstm_copy(a, b) (pstm_copy(a, b), PSTM_OKAY)
215 /******************************************************************************/
219 This is used to ensure that leading zero digits are trimed and the
220 leading "used" digit will be non-zero. Typically very fast. Also fixes
221 the sign if there are no more leading digits
223 void FAST_FUNC pstm_clamp(pstm_int * a)
225 /* decrease used while the most significant digit is zero. */
226 while (a->used > 0 && a->dp[a->used - 1] == 0) {
229 /* reset the sign flag if used == 0 */
235 /******************************************************************************/
239 void FAST_FUNC pstm_clear(pstm_int * a)
243 only do anything if a hasn't been freed previously
245 if (a != NULL && a->dp != NULL) {
247 first zero the digits
249 for (i = 0; i < a->used; i++) {
253 psFree (a->dp, a->pool);
255 reset members to make debugging easier
258 a->alloc = a->used = 0;
263 /******************************************************************************/
268 void pstm_clear_multi(pstm_int *mp0, pstm_int *mp1, pstm_int *mp2,
269 pstm_int *mp3, pstm_int *mp4, pstm_int *mp5,
270 pstm_int *mp6, pstm_int *mp7)
272 int32 n; /* Number of ok inits */
274 pstm_int *tempArray[9];
286 for (n = 0; tempArray[n] != NULL; n++) {
287 if ((tempArray[n] != NULL) && (tempArray[n]->dp != NULL)) {
288 pstm_clear(tempArray[n]);
294 /******************************************************************************/
298 static void pstm_zero(pstm_int * a)
307 for (n = 0; n < a->alloc; n++) {
313 /******************************************************************************/
315 Compare maginitude of two ints (unsigned).
317 int32 FAST_FUNC pstm_cmp_mag(pstm_int * a, pstm_int * b)
319 int n; //bbox: was int16
320 pstm_digit *tmpa, *tmpb;
323 compare based on # of non-zero digits
325 if (a->used > b->used) {
329 if (a->used < b->used) {
334 tmpa = a->dp + (a->used - 1);
337 tmpb = b->dp + (a->used - 1);
340 compare based on digits
342 for (n = 0; n < a->used; ++n, --tmpa, --tmpb) {
353 /******************************************************************************/
355 Compare two ints (signed)
357 int32 FAST_FUNC pstm_cmp(pstm_int * a, pstm_int * b)
360 compare based on sign
362 if (a->sign != b->sign) {
363 if (a->sign == PSTM_NEG) {
372 if (a->sign == PSTM_NEG) {
373 /* if negative compare opposite direction */
374 return pstm_cmp_mag(b, a);
376 return pstm_cmp_mag(a, b);
380 /******************************************************************************/
382 pstm_ints can be initialized more precisely when they will populated
383 using pstm_read_unsigned_bin since the length of the byte stream is known
385 int32 FAST_FUNC pstm_init_for_read_unsigned_bin(psPool_t *pool, pstm_int *a, uint32 len)
389 Need to set this based on how many words max it will take to store the bin.
391 1 to round up for the remainder of this integer math
392 1 for the initial carry of '1' bits that fall between DIGIT_BIT and 8
394 size = (((len / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
396 return pstm_init_size(pool, a, size);
400 /******************************************************************************/
402 Reads a unsigned char array into pstm_int format. User should have
403 called pstm_init_for_read_unsigned_bin first. There is some grow logic
404 here if the default pstm_init was used but we don't really want to hit it.
406 int32 FAST_FUNC pstm_read_unsigned_bin(pstm_int *a, unsigned char *b, int32 c)
412 If we know the endianness of this architecture, and we're using
413 32-bit pstm_digits, we can optimize this
415 #if (defined(ENDIAN_LITTLE) || defined(ENDIAN_BIG)) && !defined(PSTM_64BIT)
416 /* But not for both simultaneously */
417 #if defined(ENDIAN_LITTLE) && defined(ENDIAN_BIG)
418 #error Both ENDIAN_LITTLE and ENDIAN_BIG defined.
422 if ((unsigned)c > (PSTM_MAX_SIZE * sizeof(pstm_digit))) {
423 uint32 excess = c - (PSTM_MAX_SIZE * sizeof(pstm_digit));
427 a->used = ((c + sizeof(pstm_digit) - 1)/sizeof(pstm_digit));
428 if (a->alloc < a->used) {
429 if (pstm_grow(a, a->used) != PSTM_OKAY) {
433 pd = (unsigned char *)a->dp;
434 /* read the bytes in */
437 /* Use Duff's device to unroll the loop. */
438 int32 idx = (c - 1) & ~3;
440 case 0: do { pd[idx+0] = *b++;
441 case 3: pd[idx+1] = *b++;
442 case 2: pd[idx+2] = *b++;
443 case 1: pd[idx+3] = *b++;
445 } while ((c -= 4) > 0);
449 for (c -= 1; c >= 0; c -= 1) {
455 /* Big enough based on the len? */
456 a->used = (((c / sizeof(pstm_digit)) * (sizeof(pstm_digit) * CHAR_BIT))
459 if (a->alloc < a->used) {
460 if (pstm_grow(a, a->used) != PSTM_OKAY) {
464 /* read the bytes in */
466 if (pstm_mul_2d (a, 8, a) != PSTM_OKAY) {
478 /******************************************************************************/
481 static int pstm_count_bits(pstm_int * a)
483 int r; //bbox: was int16
490 /* get number of digits and add that */
491 r = (a->used - 1) * DIGIT_BIT;
493 /* take the last digit and count the bits in it */
494 q = a->dp[a->used - 1];
495 while (q > ((pstm_digit) 0)) {
497 q >>= ((pstm_digit) 1);
502 /******************************************************************************/
503 int32 FAST_FUNC pstm_unsigned_bin_size(pstm_int *a)
505 int32 size = pstm_count_bits (a);
506 return (size / 8 + ((size & 7) != 0 ? 1 : 0));
509 /******************************************************************************/
510 static void pstm_set(pstm_int *a, pstm_digit b)
514 a->used = a->dp[0] ? 1 : 0;
517 /******************************************************************************/
521 static void pstm_rshd(pstm_int *a, int x)
523 int y; //bbox: was int16
525 /* too many digits just zero and return */
532 for (y = 0; y < a->used - x; y++) {
533 a->dp[y] = a->dp[y+x];
537 for (; y < a->used; y++) {
541 /* decrement count */
546 /******************************************************************************/
548 Shift left a certain amount of digits.
551 static int32 pstm_lshd(pstm_int * a, int b)
553 int x; //bbox: was int16
557 If its less than zero return.
563 Grow to fit the new digits.
565 if (a->alloc < a->used + b) {
566 if ((res = pstm_grow (a, a->used + b)) != PSTM_OKAY) {
572 register pstm_digit *top, *bottom;
574 Increment the used by the shift amount then copy upwards.
579 top = a->dp + a->used - 1;
582 bottom = a->dp + a->used - 1 - b;
584 This is implemented using a sliding window except the window goes the
585 other way around. Copying from the bottom to the top.
587 for (x = a->used - 1; x >= b; x--) {
591 /* zero the lower digits */
593 for (x = 0; x < b; x++) {
599 #define pstm_lshd(a, b) (pstm_lshd(a, b), PSTM_OKAY)
601 /******************************************************************************/
605 static int32 pstm_2expt(pstm_int *a, int b)
607 int z; //bbox: was int16
609 /* zero a as per default */
617 if (z >= PSTM_MAX_SIZE) {
618 return PS_LIMIT_FAIL;
621 /* set the used count of where the bit will go */
624 if (a->used > a->alloc) {
625 if (pstm_grow(a, a->used) != PSTM_OKAY) {
630 /* put the single bit in its place */
631 a->dp[z] = ((pstm_digit)1) << (b % DIGIT_BIT);
635 /******************************************************************************/
639 int32 FAST_FUNC pstm_mul_2(pstm_int * a, pstm_int * b)
642 int x, oldused; //bbox: was int16
645 grow to accomodate result
647 if (b->alloc < a->used + 1) {
648 if ((res = pstm_grow (b, a->used + 1)) != PSTM_OKAY) {
656 register pstm_digit r, rr, *tmpa, *tmpb;
658 /* alias for source */
666 for (x = 0; x < a->used; x++) {
668 get what will be the *next* carry bit from the
669 MSB of the current digit
671 rr = *tmpa >> ((pstm_digit)(DIGIT_BIT - 1));
673 now shift up this digit, add in the carry [from the previous]
675 *tmpb++ = ((*tmpa++ << ((pstm_digit)1)) | r);
677 copy the carry that would be from the source
678 digit into the next iteration
683 /* new leading digit? */
684 if (r != 0 && b->used != (PSTM_MAX_SIZE-1)) {
685 /* add a MSB which is always 1 at this point */
690 now zero any excess digits on the destination that we didn't write to
692 tmpb = b->dp + b->used;
693 for (x = b->used; x < oldused; x++) {
701 /******************************************************************************/
703 unsigned subtraction ||a|| >= ||b|| ALWAYS!
705 int32 FAST_FUNC s_pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
707 int oldbused, oldused; //bbox: was int16
711 if (b->used > a->used) {
712 return PS_LIMIT_FAIL;
714 if (c->alloc < a->used) {
715 if ((x = pstm_grow (c, a->used)) != PSTM_OKAY) {
724 for (x = 0; x < oldbused; x++) {
725 t = ((pstm_word)a->dp[x]) - (((pstm_word)b->dp[x]) + t);
726 c->dp[x] = (pstm_digit)t;
727 t = (t >> DIGIT_BIT)&1;
729 for (; x < a->used; x++) {
730 t = ((pstm_word)a->dp[x]) - t;
731 c->dp[x] = (pstm_digit)t;
732 t = (t >> DIGIT_BIT);
734 for (; x < oldused; x++) {
741 /******************************************************************************/
745 static int32 s_pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
747 int x, y, oldused; //bbox: was int16
748 register pstm_word t, adp, bdp;
757 if (c->used > c->alloc) {
758 if (pstm_grow(c, c->used) != PSTM_OKAY) {
764 for (x = 0; x < y; x++) {
768 adp = (pstm_word)a->dp[x];
773 bdp = (pstm_word)b->dp[x];
776 c->dp[x] = (pstm_digit)t;
779 if (t != 0 && x < PSTM_MAX_SIZE) {
780 if (c->used == c->alloc) {
781 if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
785 c->dp[c->used++] = (pstm_digit)t;
790 for (; x < oldused; x++) {
798 /******************************************************************************/
802 int32 FAST_FUNC pstm_sub(pstm_int *a, pstm_int *b, pstm_int *c)
805 int sa, sb; //bbox: was int16
812 subtract a negative from a positive, OR a positive from a negative.
813 For both, ADD their magnitudes, and use the sign of the first number.
816 if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
821 subtract a positive from a positive, OR a negative from a negative.
822 First, take the difference between their magnitudes, then...
824 if (pstm_cmp_mag (a, b) != PSTM_LT) {
825 /* Copy the sign from the first */
827 /* The first has a larger or equal magnitude */
828 if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
832 /* The result has the _opposite_ sign from the first number. */
833 c->sign = (sa == PSTM_ZPOS) ? PSTM_NEG : PSTM_ZPOS;
834 /* The second has a larger magnitude */
835 if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
843 /******************************************************************************/
848 int32 pstm_sub_d(psPool_t *pool, pstm_int *a, pstm_digit b, pstm_int *c)
853 if (pstm_init_size(pool, &tmp, sizeof(pstm_digit)) != PSTM_OKAY) {
857 res = pstm_sub(a, &tmp, c);
863 /******************************************************************************/
865 setups the montgomery reduction
867 static int32 pstm_montgomery_setup(pstm_int *a, pstm_digit *rho)
872 fast inversion mod 2**k
873 Based on the fact that
874 XA = 1 (mod 2**n) => (X(2-XA)) A = 1 (mod 2**2n)
875 => 2*X*A - X*X*A*A = 1
881 psTraceCrypto("pstm_montogomery_setup failure\n");
885 x = (((b + 2) & 4) << 1) + b; /* here x*a==1 mod 2**4 */
886 x *= 2 - b * x; /* here x*a==1 mod 2**8 */
887 x *= 2 - b * x; /* here x*a==1 mod 2**16 */
888 x *= 2 - b * x; /* here x*a==1 mod 2**32 */
890 x *= 2 - b * x; /* here x*a==1 mod 2**64 */
892 /* rho = -1/m mod b */
893 *rho = (pstm_digit)(((pstm_word) 1 << ((pstm_word) DIGIT_BIT)) -
898 /******************************************************************************/
900 * computes a = B**n mod b without division or multiplication useful for
901 * normalizing numbers in a Montgomery system.
903 static int32 pstm_montgomery_calc_normalization(pstm_int *a, pstm_int *b)
906 int bits; //bbox: was int16
908 /* how many bits of last digit does b use */
909 bits = pstm_count_bits (b) % DIGIT_BIT;
910 if (!bits) bits = DIGIT_BIT;
912 /* compute A = B^(n-1) * 2^(bits-1) */
914 if ((x = pstm_2expt (a, (b->used - 1) * DIGIT_BIT + bits - 1)) !=
923 /* now compute C = A * B mod b */
924 for (x = bits - 1; x < (int32)DIGIT_BIT; x++) {
925 if (pstm_mul_2 (a, a) != PSTM_OKAY) {
928 if (pstm_cmp_mag (a, b) != PSTM_LT) {
929 if (s_pstm_sub (a, b, a) != PSTM_OKAY) {
937 /******************************************************************************/
942 static int32 pstm_mul_2d(pstm_int *a, int b, pstm_int *c)
944 pstm_digit carry, carrytmp, shift;
945 int x; //bbox: was int16
948 if (pstm_copy(a, c) != PSTM_OKAY) {
952 /* handle whole digits */
953 if (b >= DIGIT_BIT) {
954 if (pstm_lshd(c, b/DIGIT_BIT) != PSTM_OKAY) {
960 /* shift the digits */
963 shift = DIGIT_BIT - b;
964 for (x = 0; x < c->used; x++) {
965 carrytmp = c->dp[x] >> shift;
966 c->dp[x] = (c->dp[x] << b) + carry;
969 /* store last carry if room */
970 if (carry && x < PSTM_MAX_SIZE) {
971 if (c->used == c->alloc) {
972 if (pstm_grow(c, c->alloc + 1) != PSTM_OKAY) {
976 c->dp[c->used++] = carry;
982 #define pstm_mul_2d(a, b, c) (pstm_mul_2d(a, b, c), PSTM_OKAY)
984 /******************************************************************************/
989 static int32 pstm_mod_2d(pstm_int *a, int b, pstm_int *c) //bbox: was int16 b
991 int x; //bbox: was int16
993 /* zero if count less than or equal to zero */
999 /* get copy of input */
1000 if (pstm_copy(a, c) != PSTM_OKAY) {
1004 /* if 2**d is larger than we just return */
1005 if (b >= (DIGIT_BIT * a->used)) {
1009 /* zero digits above the last digit of the modulus */
1010 for (x = (b / DIGIT_BIT) + ((b % DIGIT_BIT) == 0 ? 0 : 1); x < c->used; x++)
1014 /* clear the digit that is not completely outside/inside the modulus */
1015 c->dp[b / DIGIT_BIT] &= ~((pstm_digit)0) >> (DIGIT_BIT - b);
1019 #define pstm_mod_2d(a, b, c) (pstm_mod_2d(a, b, c), PSTM_OKAY)
1022 /******************************************************************************/
1027 static int32 pstm_mul_d(pstm_int *a, pstm_digit b, pstm_int *c)
1031 int x, oldused; //bbox: was int16
1033 if (c->alloc < a->used + 1) {
1034 if ((res = pstm_grow (c, a->used + 1)) != PSTM_OKAY) {
1042 for (x = 0; x < a->used; x++) {
1043 w = ((pstm_word)a->dp[x]) * ((pstm_word)b) + w;
1044 c->dp[x] = (pstm_digit)w;
1047 if (w != 0 && (a->used != PSTM_MAX_SIZE)) {
1048 c->dp[c->used++] = (pstm_digit)w;
1051 for (; x < oldused; x++) {
1057 #define pstm_mul_d(a, b, c) (pstm_mul_d(a, b, c), PSTM_OKAY)
1059 /******************************************************************************/
1064 #define pstm_div_2d(pool, a, b, c, d) \
1065 pstm_div_2d( a, b, c, d)
1066 static int32 pstm_div_2d(psPool_t *pool, pstm_int *a, int b, pstm_int *c,
1069 pstm_digit D, r, rr;
1071 int x; //bbox: was int16
1074 /* if the shift count is <= 0 then we do no work */
1076 if (pstm_copy (a, c) != PSTM_OKAY) {
1085 /* get the remainder */
1087 if (pstm_init(pool, &t) != PSTM_OKAY) {
1090 if (pstm_mod_2d (a, b, &t) != PSTM_OKAY) {
1097 if (pstm_copy(a, c) != PSTM_OKAY) {
1102 /* shift by as many digits in the bit count */
1103 if (b >= (int32)DIGIT_BIT) {
1104 pstm_rshd (c, b / DIGIT_BIT);
1107 /* shift any bit count < DIGIT_BIT */
1108 D = (pstm_digit) (b % DIGIT_BIT);
1110 register pstm_digit *tmpc, mask, shift;
1113 mask = (((pstm_digit)1) << D) - 1;
1116 shift = DIGIT_BIT - D;
1119 tmpc = c->dp + (c->used - 1);
1123 for (x = c->used - 1; x >= 0; x--) {
1124 /* get the lower bits of this word in a temp */
1127 /* shift the current word and mix in the carry bits from previous */
1128 *tmpc = (*tmpc >> D) | (r << shift);
1131 /* set the carry to the carry bits of the current word above */
1140 if (pstm_copy(&t, d) != PSTM_OKAY) {
1148 #define pstm_div_2d(pool, a, b, c, d) (pstm_div_2d(a, b, c, d), PSTM_OKAY)
1150 /******************************************************************************/
1155 int32 pstm_div_2(pstm_int * a, pstm_int * b)
1157 int x, oldused; //bbox: was int16
1159 if (b->alloc < a->used) {
1160 if (pstm_grow(b, a->used) != PSTM_OKAY) {
1167 register pstm_digit r, rr, *tmpa, *tmpb;
1170 tmpa = a->dp + b->used - 1;
1173 tmpb = b->dp + b->used - 1;
1177 for (x = b->used - 1; x >= 0; x--) {
1178 /* get the carry for the next iteration */
1181 /* shift the current digit, add in carry and store */
1182 *tmpb-- = (*tmpa-- >> 1) | (r << (DIGIT_BIT - 1));
1184 /* forward carry to next iteration */
1188 /* zero excess digits */
1189 tmpb = b->dp + b->used;
1190 for (x = b->used; x < oldused; x++) {
1200 /******************************************************************************/
1202 Creates "a" then copies b into it
1204 #undef pstm_init_copy
1205 #define pstm_init_copy(pool, a, b, toSqr) \
1206 pstm_init_copy( a, b, toSqr)
1207 static int32 pstm_init_copy(psPool_t *pool, pstm_int * a, pstm_int * b, int toSqr)
1209 int x; //bbox: was int16
1219 Smart-size: Increasing size of a if b->used is roughly half
1220 of b->alloc because usage has shown that a lot of these copies
1221 go on to be squared and need these extra digits
1223 if ((b->used * 2) + 2 >= x) {
1224 x = (b->used * 2) + 3;
1227 if ((res = pstm_init_size(pool, a, x)) != PSTM_OKAY) {
1230 return pstm_copy(b, a);
1232 #undef pstm_init_copy
1233 #define pstm_init_copy(pool, a, b, toSqr) (pstm_init_copy(a, b, toSqr), PSTM_OKAY)
1235 /******************************************************************************/
1237 With some compilers, we have seen issues linking with the builtin
1238 64 bit division routine. The issues with either manifest in a failure
1239 to find 'udivdi3' at link time, or a runtime invalid instruction fault
1240 during an RSA operation.
1241 The routine below divides a 64 bit unsigned int by a 32 bit unsigned int
1242 explicitly, rather than using the division operation
1243 The 64 bit result is placed in the 'numerator' parameter
1244 The 32 bit mod (remainder) of the division is the return parameter
1245 Based on implementations by:
1246 Copyright (C) 2003 Bernardo Innocenti <bernie@develer.com>
1247 Copyright (C) 1999 Hewlett-Packard Co
1248 Copyright (C) 1999 David Mosberger-Tang <davidm@hpl.hp.com>
1250 #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1251 static uint32 psDiv64(uint64 *numerator, uint32 denominator)
1253 uint64 rem = *numerator;
1254 uint64 b = denominator;
1257 uint32 high = rem >> 32;
1259 if (high >= denominator) {
1260 high /= denominator;
1261 res = (uint64) high << 32;
1262 rem -= (uint64) (high * denominator) << 32;
1264 while ((int64)b > 0 && b < rem) {
1279 #endif /* USE_MATRIX_DIV64 */
1281 #if defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1282 typedef unsigned long uint128 __attribute__ ((mode(TI)));
1283 static uint64 psDiv128(uint128 *numerator, uint64 denominator)
1285 uint128 rem = *numerator;
1286 uint128 b = denominator;
1289 uint64 high = rem >> 64;
1291 if (high >= denominator) {
1292 high /= denominator;
1293 res = (uint128) high << 64;
1294 rem -= (uint128) (high * denominator) << 64;
1296 while ((uint128)b > 0 && b < rem) {
1311 #endif /* USE_MATRIX_DIV128 */
1313 /******************************************************************************/
1317 static int32 pstm_div(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
1320 pstm_int q, x, y, t1, t2;
1322 int n, t, i, norm, neg; //bbox: was int16
1324 /* is divisor zero ? */
1325 if (pstm_iszero (b) == 1) {
1326 return PS_LIMIT_FAIL;
1329 /* if a < b then q=0, r = a */
1330 if (pstm_cmp_mag (a, b) == PSTM_LT) {
1332 if (pstm_copy(a, d) != PSTM_OKAY) {
1344 if ((res = pstm_init_size(pool, &t1, a->alloc)) != PSTM_OKAY) {
1347 if ((res = pstm_init_size(pool, &t2, 3)) != PSTM_OKAY) {
1350 if ((res = pstm_init_copy(pool, &x, a, 0)) != PSTM_OKAY) {
1354 Used to be an init_copy on b but pstm_grow was always hit with triple size
1356 if ((res = pstm_init_size(pool, &y, b->used * 3)) != PSTM_OKAY) {
1359 if ((res = pstm_copy(b, &y)) != PSTM_OKAY) {
1364 neg = (a->sign == b->sign) ? PSTM_ZPOS : PSTM_NEG;
1365 x.sign = y.sign = PSTM_ZPOS;
1367 /* normalize both x and y, ensure that y >= b/2, [b == 2**DIGIT_BIT] */
1368 norm = pstm_count_bits(&y) % DIGIT_BIT;
1369 if (norm < (int32)(DIGIT_BIT-1)) {
1370 norm = (DIGIT_BIT-1) - norm;
1371 if ((res = pstm_mul_2d(&x, norm, &x)) != PSTM_OKAY) {
1374 if ((res = pstm_mul_2d(&y, norm, &y)) != PSTM_OKAY) {
1381 /* note hac does 0 based, so if used==5 then its 0,1,2,3,4, e.g. use 4 */
1385 if ((res = pstm_init_size(pool, &q, n - t + 1)) != PSTM_OKAY) {
1390 /* while (x >= y*b**n-t) do { q[n-t] += 1; x -= y*b**{n-t} } */
1391 if ((res = pstm_lshd(&y, n - t)) != PSTM_OKAY) { /* y = y*b**{n-t} */
1395 while (pstm_cmp (&x, &y) != PSTM_LT) {
1397 if ((res = pstm_sub(&x, &y, &x)) != PSTM_OKAY) {
1402 /* reset y by shifting it back down */
1403 pstm_rshd (&y, n - t);
1405 /* step 3. for i from n down to (t + 1) */
1406 for (i = n; i >= (t + 1); i--) {
1411 /* step 3.1 if xi == yt then set q{i-t-1} to b-1,
1412 * otherwise set q{i-t-1} to (xi*b + x{i-1})/yt */
1413 if (x.dp[i] == y.dp[t]) {
1414 q.dp[i - t - 1] = (pstm_digit)((((pstm_word)1) << DIGIT_BIT) - 1);
1417 tmp = ((pstm_word) x.dp[i]) << ((pstm_word) DIGIT_BIT);
1418 tmp |= ((pstm_word) x.dp[i - 1]);
1419 #if defined(USE_MATRIX_DIV64) && defined(PSTM_32BIT)
1420 psDiv64(&tmp, y.dp[t]);
1421 #elif defined(USE_MATRIX_DIV128) && defined(PSTM_64BIT)
1422 psDiv128(&tmp, y.dp[t]);
1424 tmp /= ((pstm_word) y.dp[t]);
1425 #endif /* USE_MATRIX_DIV64 */
1426 q.dp[i - t - 1] = (pstm_digit) (tmp);
1429 /* while (q{i-t-1} * (yt * b + y{t-1})) >
1430 xi * b**2 + xi-1 * b + xi-2
1434 q.dp[i - t - 1] = (q.dp[i - t - 1] + 1);
1436 q.dp[i - t - 1] = (q.dp[i - t - 1] - 1);
1438 /* find left hand */
1440 t1.dp[0] = (t - 1 < 0) ? 0 : y.dp[t - 1];
1443 if ((res = pstm_mul_d (&t1, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
1447 /* find right hand */
1448 t2.dp[0] = (i - 2 < 0) ? 0 : x.dp[i - 2];
1449 t2.dp[1] = (i - 1 < 0) ? 0 : x.dp[i - 1];
1452 } while (pstm_cmp_mag(&t1, &t2) == PSTM_GT);
1454 /* step 3.3 x = x - q{i-t-1} * y * b**{i-t-1} */
1455 if ((res = pstm_mul_d(&y, q.dp[i - t - 1], &t1)) != PSTM_OKAY) {
1459 if ((res = pstm_lshd(&t1, i - t - 1)) != PSTM_OKAY) {
1463 if ((res = pstm_sub(&x, &t1, &x)) != PSTM_OKAY) {
1467 /* if x < 0 then { x = x + y*b**{i-t-1}; q{i-t-1} -= 1; } */
1468 if (x.sign == PSTM_NEG) {
1469 if ((res = pstm_copy(&y, &t1)) != PSTM_OKAY) {
1472 if ((res = pstm_lshd (&t1, i - t - 1)) != PSTM_OKAY) {
1475 if ((res = pstm_add (&x, &t1, &x)) != PSTM_OKAY) {
1478 q.dp[i - t - 1] = q.dp[i - t - 1] - 1;
1482 now q is the quotient and x is the remainder (which we have to normalize)
1484 /* get sign before writing to c */
1485 x.sign = x.used == 0 ? PSTM_ZPOS : a->sign;
1489 if (pstm_copy (&q, c) != PSTM_OKAY) {
1497 if ((res = pstm_div_2d (pool, &x, norm, &x, NULL)) != PSTM_OKAY) {
1501 the following is a kludge, essentially we were seeing the right
1502 remainder but with excess digits that should have been zero
1504 for (i = b->used; i < x.used; i++) {
1508 if (pstm_copy (&x, d) != PSTM_OKAY) {
1516 LBL_Q:pstm_clear (&q);
1517 LBL_Y:pstm_clear (&y);
1518 LBL_X:pstm_clear (&x);
1519 LBL_T2:pstm_clear (&t2);
1520 LBL_T1:pstm_clear (&t1);
1525 /******************************************************************************/
1527 Swap the elements of two integers, for cases where you can't simply swap
1528 the pstm_int pointers around
1530 static void pstm_exch(pstm_int * a, pstm_int * b)
1539 /******************************************************************************/
1541 c = a mod b, 0 <= c < b
1543 static int32 pstm_mod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
1550 if ((err = pstm_init_size(pool, &t, b->alloc)) != PSTM_OKAY) {
1553 if ((err = pstm_div(pool, a, b, NULL, &t)) != PSTM_OKAY) {
1557 if (t.sign != b->sign) {
1558 err = pstm_add(&t, b, c);
1566 /******************************************************************************/
1570 int32 FAST_FUNC pstm_mulmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c,
1574 int size; //bbox: was int16
1578 Smart-size pstm_inits. d is an output that is influenced by this local 't'
1579 so don't shrink 'd' if it wants to becuase this will lead to an pstm_grow
1582 size = a->used + b->used + 1;
1583 if ((a == d) && (size < a->alloc)) {
1586 if ((res = pstm_init_size(pool, &tmp, size)) != PSTM_OKAY) {
1589 if ((res = pstm_mul_comba(pool, a, b, &tmp, NULL, 0)) != PSTM_OKAY) {
1593 res = pstm_mod(pool, &tmp, c, d);
1598 /******************************************************************************/
1601 * Some restrictions... x must be positive and < b
1603 int32 FAST_FUNC pstm_exptmod(psPool_t *pool, pstm_int *G, pstm_int *X, pstm_int *P,
1606 pstm_int M[32], res; /* Keep this winsize based: (1 << max_winsize) */
1610 int bitcpy, bitcnt, mode, digidx, x, y, winsize; //bbox: was int16
1613 /* set window size from what user set as optimization */
1614 x = pstm_count_bits(X);
1618 winsize = PS_EXPTMOD_WINSIZE;
1621 /* now setup montgomery */
1622 if ((err = pstm_montgomery_setup (P, &mp)) != PSTM_OKAY) {
1627 if ((err = pstm_init_size(pool, &res, (P->used * 2) + 1)) != PSTM_OKAY) {
1632 The M table contains powers of the input base, e.g. M[x] = G^x mod P
1633 The first half of the table is not computed though except for M[0] and M[1]
1635 /* now we need R mod m */
1636 if ((err = pstm_montgomery_calc_normalization (&res, P)) != PSTM_OKAY) {
1643 if ((err = pstm_init_size(pool, &M[1], res.used)) != PSTM_OKAY) {
1647 /* now set M[1] to G * R mod m */
1648 if (pstm_cmp_mag(P, G) != PSTM_GT) {
1649 /* G > P so we reduce it first */
1650 if ((err = pstm_mod(pool, G, P, &M[1])) != PSTM_OKAY) {
1654 if ((err = pstm_copy(G, &M[1])) != PSTM_OKAY) {
1658 if ((err = pstm_mulmod (pool, &M[1], &res, P, &M[1])) != PSTM_OKAY) {
1662 Pre-allocated digit. Used for mul, sqr, AND reduce
1664 paDlen = ((M[1].used + 3) * 2) * sizeof(pstm_digit);
1665 paD = xzalloc(paDlen);//bbox
1667 compute the value at M[1<<(winsize-1)] by squaring M[1] (winsize-1) times
1669 if (pstm_init_copy(pool, &M[1 << (winsize - 1)], &M[1], 1) != PSTM_OKAY) {
1673 for (x = 0; x < (winsize - 1); x++) {
1674 if ((err = pstm_sqr_comba (pool, &M[1 << (winsize - 1)],
1675 &M[1 << (winsize - 1)], paD, paDlen)) != PSTM_OKAY) {
1678 if ((err = pstm_montgomery_reduce(pool, &M[1 << (winsize - 1)], P, mp,
1679 paD, paDlen)) != PSTM_OKAY) {
1684 now init the second half of the array
1686 for (x = (1<<(winsize-1)) + 1; x < (1 << winsize); x++) {
1687 if ((err = pstm_init_size(pool, &M[x], M[1<<(winsize-1)].alloc + 1))
1689 for (y = 1<<(winsize-1); y < x; y++) {
1696 /* create upper table */
1697 for (x = (1 << (winsize - 1)) + 1; x < (1 << winsize); x++) {
1698 if ((err = pstm_mul_comba(pool, &M[x - 1], &M[1], &M[x], paD, paDlen))
1702 if ((err = pstm_montgomery_reduce(pool, &M[x], P, mp, paD, paDlen)) !=
1708 /* set initial mode and bit cnt */
1712 digidx = X->used - 1;
1717 /* grab next digit as required */
1718 if (--bitcnt == 0) {
1719 /* if digidx == -1 we are out of digits so break */
1723 /* read next digit and reset bitcnt */
1724 buf = X->dp[digidx--];
1725 bitcnt = (int32)DIGIT_BIT;
1728 /* grab the next msb from the exponent */
1729 y = (pstm_digit)(buf >> (DIGIT_BIT - 1)) & 1;
1730 buf <<= (pstm_digit)1;
1732 If the bit is zero and mode == 0 then we ignore it.
1733 These represent the leading zero bits before the first 1 bit
1734 in the exponent. Technically this opt is not required but it
1735 does lower the # of trivial squaring/reductions used
1737 if (mode == 0 && y == 0) {
1741 /* if the bit is zero and mode == 1 then we square */
1742 if (mode == 1 && y == 0) {
1743 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1747 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1754 /* else we add it to the window */
1755 bitbuf |= (y << (winsize - ++bitcpy));
1758 if (bitcpy == winsize) {
1759 /* ok window is filled so square as required and mul square first */
1760 for (x = 0; x < winsize; x++) {
1761 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1765 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
1766 paDlen)) != PSTM_OKAY) {
1772 if ((err = pstm_mul_comba(pool, &res, &M[bitbuf], &res, paD,
1773 paDlen)) != PSTM_OKAY) {
1776 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1781 /* empty window and reset */
1788 /* if bits remain then square/multiply */
1789 if (mode == 2 && bitcpy > 0) {
1790 /* square then multiply if the bit is set */
1791 for (x = 0; x < bitcpy; x++) {
1792 if ((err = pstm_sqr_comba(pool, &res, &res, paD, paDlen)) !=
1796 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen))
1801 /* get next bit of the window */
1803 if ((bitbuf & (1 << winsize)) != 0) {
1805 if ((err = pstm_mul_comba(pool, &res, &M[1], &res, paD, paDlen))
1809 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD,
1810 paDlen)) != PSTM_OKAY) {
1817 Fix up result if Montgomery reduction is used recall that any value in a
1818 Montgomery system is actually multiplied by R mod n. So we have to reduce
1819 one more time to cancel out the factor of R.
1821 if ((err = pstm_montgomery_reduce(pool, &res, P, mp, paD, paDlen)) !=
1825 /* swap res with Y */
1826 if ((err = pstm_copy (&res, Y)) != PSTM_OKAY) {
1831 for (x = 1<<(winsize-1); x < (1 << winsize); x++) {
1834 LBL_PAD:psFree(paD, pool);
1835 LBL_M: pstm_clear(&M[1]);
1836 LBL_RES:pstm_clear(&res);
1840 /******************************************************************************/
1844 int32 FAST_FUNC pstm_add(pstm_int *a, pstm_int *b, pstm_int *c)
1847 int sa, sb; //bbox: was int16
1849 /* get sign of both inputs */
1853 /* handle two cases, not four */
1855 /* both positive or both negative, add their mags, copy the sign */
1857 if ((res = s_pstm_add (a, b, c)) != PSTM_OKAY) {
1862 one positive, the other negative
1863 subtract the one with the greater magnitude from the one of the lesser
1864 magnitude. The result gets the sign of the one with the greater mag.
1866 if (pstm_cmp_mag (a, b) == PSTM_LT) {
1868 if ((res = s_pstm_sub (b, a, c)) != PSTM_OKAY) {
1873 if ((res = s_pstm_sub (a, b, c)) != PSTM_OKAY) {
1881 /******************************************************************************/
1883 reverse an array, used for radix code
1885 static void pstm_reverse (unsigned char *s, int len) //bbox: was int16 len
1900 /******************************************************************************/
1902 No reverse. Useful in some of the EIP-154 PKA stuff where special byte
1903 order seems to come into play more often
1906 int32 pstm_to_unsigned_bin_nr(psPool_t *pool, pstm_int *a, unsigned char *b)
1909 int x; //bbox: was int16
1912 if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
1917 while (pstm_iszero (&t) == 0) {
1918 b[x++] = (unsigned char) (t.dp[0] & 255);
1919 if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
1928 /******************************************************************************/
1932 int32 FAST_FUNC pstm_to_unsigned_bin(psPool_t *pool, pstm_int *a, unsigned char *b)
1935 int x; //bbox: was int16
1938 if ((res = pstm_init_copy(pool, &t, a, 0)) != PSTM_OKAY) {
1943 while (pstm_iszero (&t) == 0) {
1944 b[x++] = (unsigned char) (t.dp[0] & 255);
1945 if ((res = pstm_div_2d (pool, &t, 8, &t, NULL)) != PSTM_OKAY) {
1950 pstm_reverse (b, x);
1956 /******************************************************************************/
1958 compare against a single digit
1960 static int32 pstm_cmp_d(pstm_int *a, pstm_digit b)
1962 /* compare based on sign */
1963 if ((b && a->used == 0) || a->sign == PSTM_NEG) {
1967 /* compare based on magnitude */
1972 /* compare the only digit of a to b */
1975 } else if (a->dp[0] < b) {
1983 Need invmod for ECC and also private key loading for hardware crypto
1984 in cases where dQ > dP. The values must be switched and a new qP must be
1985 calculated using this function
1988 #define pstm_invmod_slow(pool, a, b, c) \
1989 pstm_invmod_slow( a, b, c)
1990 static int32 pstm_invmod_slow(psPool_t *pool, pstm_int * a, pstm_int * b,
1993 pstm_int x, y, u, v, A, B, C, D;
1996 /* b cannot be negative */
1997 if (b->sign == PSTM_NEG || pstm_iszero(b) == 1) {
1998 return PS_LIMIT_FAIL;
2002 if (pstm_init_size(pool, &x, b->used) != PSTM_OKAY) {
2007 if ((res = pstm_mod(pool, a, b, &x)) != PSTM_OKAY) {
2011 if (pstm_init_copy(pool, &y, b, 0) != PSTM_OKAY) {
2015 /* 2. [modified] if x,y are both even then return an error! */
2016 if (pstm_iseven (&x) == 1 && pstm_iseven (&y) == 1) {
2021 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2022 if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
2025 if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
2029 if ((res = pstm_init_size(pool, &A, sizeof(pstm_digit))) != PSTM_OKAY) {
2033 if ((res = pstm_init_size(pool, &D, sizeof(pstm_digit))) != PSTM_OKAY) {
2039 if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
2042 if ((res = pstm_init(pool, &C)) != PSTM_OKAY) {
2047 /* 4. while u is even do */
2048 while (pstm_iseven (&u) == 1) {
2050 if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
2054 /* 4.2 if A or B is odd then */
2055 if (pstm_isodd (&A) == 1 || pstm_isodd (&B) == 1) {
2056 /* A = (A+y)/2, B = (B-x)/2 */
2057 if ((res = pstm_add (&A, &y, &A)) != PSTM_OKAY) {
2060 if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
2064 /* A = A/2, B = B/2 */
2065 if ((res = pstm_div_2 (&A, &A)) != PSTM_OKAY) {
2068 if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
2073 /* 5. while v is even do */
2074 while (pstm_iseven (&v) == 1) {
2076 if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
2080 /* 5.2 if C or D is odd then */
2081 if (pstm_isodd (&C) == 1 || pstm_isodd (&D) == 1) {
2082 /* C = (C+y)/2, D = (D-x)/2 */
2083 if ((res = pstm_add (&C, &y, &C)) != PSTM_OKAY) {
2086 if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
2090 /* C = C/2, D = D/2 */
2091 if ((res = pstm_div_2 (&C, &C)) != PSTM_OKAY) {
2094 if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
2099 /* 6. if u >= v then */
2100 if (pstm_cmp (&u, &v) != PSTM_LT) {
2101 /* u = u - v, A = A - C, B = B - D */
2102 if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
2105 if ((res = pstm_sub (&A, &C, &A)) != PSTM_OKAY) {
2108 if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
2112 /* v - v - u, C = C - A, D = D - B */
2113 if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
2116 if ((res = pstm_sub (&C, &A, &C)) != PSTM_OKAY) {
2119 if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
2124 /* if not zero goto step 4 */
2125 if (pstm_iszero (&u) == 0)
2128 /* now a = C, b = D, gcd == g*v */
2130 /* if v != 1 then there is no inverse */
2131 if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
2136 /* if its too low */
2137 while (pstm_cmp_d(&C, 0) == PSTM_LT) {
2138 if ((res = pstm_add(&C, b, &C)) != PSTM_OKAY) {
2144 while (pstm_cmp_mag(&C, b) != PSTM_LT) {
2145 if ((res = pstm_sub(&C, b, &C)) != PSTM_OKAY) {
2150 /* C is now the inverse */
2151 if ((res = pstm_copy(&C, c)) != PSTM_OKAY) {
2156 LBL_C: pstm_clear(&C);
2157 LBL_D: pstm_clear(&D);
2158 LBL_B: pstm_clear(&B);
2159 LBL_A: pstm_clear(&A);
2160 LBL_V: pstm_clear(&v);
2161 LBL_U: pstm_clear(&u);
2162 LBL_Y: pstm_clear(&y);
2163 LBL_X: pstm_clear(&x);
2168 /* c = 1/a (mod b) for odd b only */
2169 int32 pstm_invmod(psPool_t *pool, pstm_int *a, pstm_int *b, pstm_int *c)
2171 pstm_int x, y, u, v, B, D;
2173 int neg, sanity; //bbox: was uint16
2175 /* 2. [modified] b must be odd */
2176 if (pstm_iseven (b) == 1) {
2177 return pstm_invmod_slow(pool, a,b,c);
2180 /* x == modulus, y == value to invert */
2181 if ((res = pstm_init_copy(pool, &x, b, 0)) != PSTM_OKAY) {
2185 if ((res = pstm_init_size(pool, &y, a->alloc)) != PSTM_OKAY) {
2189 /* we need y = |a| */
2192 /* 3. u=x, v=y, A=1, B=0, C=0,D=1 */
2193 if ((res = pstm_init_copy(pool, &u, &x, 0)) != PSTM_OKAY) {
2196 if ((res = pstm_init_copy(pool, &v, &y, 0)) != PSTM_OKAY) {
2199 if ((res = pstm_init(pool, &B)) != PSTM_OKAY) {
2202 if ((res = pstm_init(pool, &D)) != PSTM_OKAY) {
2210 /* 4. while u is even do */
2211 while (pstm_iseven (&u) == 1) {
2213 if ((res = pstm_div_2 (&u, &u)) != PSTM_OKAY) {
2217 /* 4.2 if B is odd then */
2218 if (pstm_isodd (&B) == 1) {
2219 if ((res = pstm_sub (&B, &x, &B)) != PSTM_OKAY) {
2224 if ((res = pstm_div_2 (&B, &B)) != PSTM_OKAY) {
2229 /* 5. while v is even do */
2230 while (pstm_iseven (&v) == 1) {
2232 if ((res = pstm_div_2 (&v, &v)) != PSTM_OKAY) {
2235 /* 5.2 if D is odd then */
2236 if (pstm_isodd (&D) == 1) {
2238 if ((res = pstm_sub (&D, &x, &D)) != PSTM_OKAY) {
2243 if ((res = pstm_div_2 (&D, &D)) != PSTM_OKAY) {
2248 /* 6. if u >= v then */
2249 if (pstm_cmp (&u, &v) != PSTM_LT) {
2250 /* u = u - v, B = B - D */
2251 if ((res = pstm_sub (&u, &v, &u)) != PSTM_OKAY) {
2254 if ((res = pstm_sub (&B, &D, &B)) != PSTM_OKAY) {
2258 /* v - v - u, D = D - B */
2259 if ((res = pstm_sub (&v, &u, &v)) != PSTM_OKAY) {
2262 if ((res = pstm_sub (&D, &B, &D)) != PSTM_OKAY) {
2267 /* if not zero goto step 4 */
2268 if (sanity++ > 1000) {
2269 res = PS_LIMIT_FAIL;
2272 if (pstm_iszero (&u) == 0) {
2276 /* now a = C, b = D, gcd == g*v */
2278 /* if v != 1 then there is no inverse */
2279 if (pstm_cmp_d (&v, 1) != PSTM_EQ) {
2284 /* b is now the inverse */
2286 while (D.sign == PSTM_NEG) {
2287 if ((res = pstm_add (&D, b, &D)) != PSTM_OKAY) {
2291 if ((res = pstm_copy (&D, c)) != PSTM_OKAY) {
2297 LBL_D: pstm_clear(&D);
2298 LBL_B: pstm_clear(&B);
2299 LBL_V: pstm_clear(&v);
2300 LBL_U: pstm_clear(&u);
2301 LBL_Y: pstm_clear(&y);
2302 LBL_X: pstm_clear(&x);
2307 #endif /* !DISABLE_PSTM */
2308 /******************************************************************************/