security/nss/lib/freebl/mpi/mpi.c@b8a032363ba2
security/nss/lib/freebl/mpi/mpi.c
Thu, 22 Jan 2015 13:21:57 +0100
- author
- Michael Schloh von Bennewitz <michael@schloh.com>
- date
- Thu, 22 Jan 2015 13:21:57 +0100
- branch
- TOR_BUG_9701
- changeset 15
- b8a032363ba2
- permissions
- -rw-r--r--
Incorporate requested changes from Mozilla in review:
https://bugzilla.mozilla.org/show_bug.cgi?id=1123480#c6
1 /*
2 * mpi.c
3 *
4 * Arbitrary precision integer arithmetic library
5 *
6 * This Source Code Form is subject to the terms of the Mozilla Public
7 * License, v. 2.0. If a copy of the MPL was not distributed with this
8 * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
10 #include "mpi-priv.h"
11 #if defined(OSF1)
12 #include <c_asm.h>
13 #endif
15 #if defined(__arm__) && \
16 ((defined(__thumb__) && !defined(__thumb2__)) || defined(__ARM_ARCH_3__))
17 /* 16-bit thumb or ARM v3 doesn't work inlined assember version */
18 #undef MP_ASSEMBLY_MULTIPLY
19 #undef MP_ASSEMBLY_SQUARE
20 #endif
22 #if MP_LOGTAB
23 /*
24 A table of the logs of 2 for various bases (the 0 and 1 entries of
25 this table are meaningless and should not be referenced).
27 This table is used to compute output lengths for the mp_toradix()
28 function. Since a number n in radix r takes up about log_r(n)
29 digits, we estimate the output size by taking the least integer
30 greater than log_r(n), where:
32 log_r(n) = log_2(n) * log_r(2)
34 This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
35 which are the output bases supported.
36 */
37 #include "logtab.h"
38 #endif
40 /* {{{ Constant strings */
42 /* Constant strings returned by mp_strerror() */
43 static const char *mp_err_string[] = {
44 "unknown result code", /* say what? */
45 "boolean true", /* MP_OKAY, MP_YES */
46 "boolean false", /* MP_NO */
47 "out of memory", /* MP_MEM */
48 "argument out of range", /* MP_RANGE */
49 "invalid input parameter", /* MP_BADARG */
50 "result is undefined" /* MP_UNDEF */
51 };
53 /* Value to digit maps for radix conversion */
55 /* s_dmap_1 - standard digits and letters */
56 static const char *s_dmap_1 =
57 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
59 /* }}} */
61 unsigned long mp_allocs;
62 unsigned long mp_frees;
63 unsigned long mp_copies;
65 /* {{{ Default precision manipulation */
67 /* Default precision for newly created mp_int's */
68 static mp_size s_mp_defprec = MP_DEFPREC;
70 mp_size mp_get_prec(void)
71 {
72 return s_mp_defprec;
74 } /* end mp_get_prec() */
76 void mp_set_prec(mp_size prec)
77 {
78 if(prec == 0)
79 s_mp_defprec = MP_DEFPREC;
80 else
81 s_mp_defprec = prec;
83 } /* end mp_set_prec() */
85 /* }}} */
87 /*------------------------------------------------------------------------*/
88 /* {{{ mp_init(mp) */
90 /*
91 mp_init(mp)
93 Initialize a new zero-valued mp_int. Returns MP_OKAY if successful,
94 MP_MEM if memory could not be allocated for the structure.
95 */
97 mp_err mp_init(mp_int *mp)
98 {
99 return mp_init_size(mp, s_mp_defprec);
101 } /* end mp_init() */
103 /* }}} */
105 /* {{{ mp_init_size(mp, prec) */
107 /*
108 mp_init_size(mp, prec)
110 Initialize a new zero-valued mp_int with at least the given
111 precision; returns MP_OKAY if successful, or MP_MEM if memory could
112 not be allocated for the structure.
113 */
115 mp_err mp_init_size(mp_int *mp, mp_size prec)
116 {
117 ARGCHK(mp != NULL && prec > 0, MP_BADARG);
119 prec = MP_ROUNDUP(prec, s_mp_defprec);
120 if((DIGITS(mp) = s_mp_alloc(prec, sizeof(mp_digit))) == NULL)
121 return MP_MEM;
123 SIGN(mp) = ZPOS;
124 USED(mp) = 1;
125 ALLOC(mp) = prec;
127 return MP_OKAY;
129 } /* end mp_init_size() */
131 /* }}} */
133 /* {{{ mp_init_copy(mp, from) */
135 /*
136 mp_init_copy(mp, from)
138 Initialize mp as an exact copy of from. Returns MP_OKAY if
139 successful, MP_MEM if memory could not be allocated for the new
140 structure.
141 */
143 mp_err mp_init_copy(mp_int *mp, const mp_int *from)
144 {
145 ARGCHK(mp != NULL && from != NULL, MP_BADARG);
147 if(mp == from)
148 return MP_OKAY;
150 if((DIGITS(mp) = s_mp_alloc(ALLOC(from), sizeof(mp_digit))) == NULL)
151 return MP_MEM;
153 s_mp_copy(DIGITS(from), DIGITS(mp), USED(from));
154 USED(mp) = USED(from);
155 ALLOC(mp) = ALLOC(from);
156 SIGN(mp) = SIGN(from);
158 return MP_OKAY;
160 } /* end mp_init_copy() */
162 /* }}} */
164 /* {{{ mp_copy(from, to) */
166 /*
167 mp_copy(from, to)
169 Copies the mp_int 'from' to the mp_int 'to'. It is presumed that
170 'to' has already been initialized (if not, use mp_init_copy()
171 instead). If 'from' and 'to' are identical, nothing happens.
172 */
174 mp_err mp_copy(const mp_int *from, mp_int *to)
175 {
176 ARGCHK(from != NULL && to != NULL, MP_BADARG);
178 if(from == to)
179 return MP_OKAY;
181 { /* copy */
182 mp_digit *tmp;
184 /*
185 If the allocated buffer in 'to' already has enough space to hold
186 all the used digits of 'from', we'll re-use it to avoid hitting
187 the memory allocater more than necessary; otherwise, we'd have
188 to grow anyway, so we just allocate a hunk and make the copy as
189 usual
190 */
191 if(ALLOC(to) >= USED(from)) {
192 s_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));
193 s_mp_copy(DIGITS(from), DIGITS(to), USED(from));
195 } else {
196 if((tmp = s_mp_alloc(ALLOC(from), sizeof(mp_digit))) == NULL)
197 return MP_MEM;
199 s_mp_copy(DIGITS(from), tmp, USED(from));
201 if(DIGITS(to) != NULL) {
202 #if MP_CRYPTO
203 s_mp_setz(DIGITS(to), ALLOC(to));
204 #endif
205 s_mp_free(DIGITS(to));
206 }
208 DIGITS(to) = tmp;
209 ALLOC(to) = ALLOC(from);
210 }
212 /* Copy the precision and sign from the original */
213 USED(to) = USED(from);
214 SIGN(to) = SIGN(from);
215 } /* end copy */
217 return MP_OKAY;
219 } /* end mp_copy() */
221 /* }}} */
223 /* {{{ mp_exch(mp1, mp2) */
225 /*
226 mp_exch(mp1, mp2)
228 Exchange mp1 and mp2 without allocating any intermediate memory
229 (well, unless you count the stack space needed for this call and the
230 locals it creates...). This cannot fail.
231 */
233 void mp_exch(mp_int *mp1, mp_int *mp2)
234 {
235 #if MP_ARGCHK == 2
236 assert(mp1 != NULL && mp2 != NULL);
237 #else
238 if(mp1 == NULL || mp2 == NULL)
239 return;
240 #endif
242 s_mp_exch(mp1, mp2);
244 } /* end mp_exch() */
246 /* }}} */
248 /* {{{ mp_clear(mp) */
250 /*
251 mp_clear(mp)
253 Release the storage used by an mp_int, and void its fields so that
254 if someone calls mp_clear() again for the same int later, we won't
255 get tollchocked.
256 */
258 void mp_clear(mp_int *mp)
259 {
260 if(mp == NULL)
261 return;
263 if(DIGITS(mp) != NULL) {
264 #if MP_CRYPTO
265 s_mp_setz(DIGITS(mp), ALLOC(mp));
266 #endif
267 s_mp_free(DIGITS(mp));
268 DIGITS(mp) = NULL;
269 }
271 USED(mp) = 0;
272 ALLOC(mp) = 0;
274 } /* end mp_clear() */
276 /* }}} */
278 /* {{{ mp_zero(mp) */
280 /*
281 mp_zero(mp)
283 Set mp to zero. Does not change the allocated size of the structure,
284 and therefore cannot fail (except on a bad argument, which we ignore)
285 */
286 void mp_zero(mp_int *mp)
287 {
288 if(mp == NULL)
289 return;
291 s_mp_setz(DIGITS(mp), ALLOC(mp));
292 USED(mp) = 1;
293 SIGN(mp) = ZPOS;
295 } /* end mp_zero() */
297 /* }}} */
299 /* {{{ mp_set(mp, d) */
301 void mp_set(mp_int *mp, mp_digit d)
302 {
303 if(mp == NULL)
304 return;
306 mp_zero(mp);
307 DIGIT(mp, 0) = d;
309 } /* end mp_set() */
311 /* }}} */
313 /* {{{ mp_set_int(mp, z) */
315 mp_err mp_set_int(mp_int *mp, long z)
316 {
317 int ix;
318 unsigned long v = labs(z);
319 mp_err res;
321 ARGCHK(mp != NULL, MP_BADARG);
323 mp_zero(mp);
324 if(z == 0)
325 return MP_OKAY; /* shortcut for zero */
327 if (sizeof v <= sizeof(mp_digit)) {
328 DIGIT(mp,0) = v;
329 } else {
330 for (ix = sizeof(long) - 1; ix >= 0; ix--) {
331 if ((res = s_mp_mul_d(mp, (UCHAR_MAX + 1))) != MP_OKAY)
332 return res;
334 res = s_mp_add_d(mp, (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));
335 if (res != MP_OKAY)
336 return res;
337 }
338 }
339 if(z < 0)
340 SIGN(mp) = NEG;
342 return MP_OKAY;
344 } /* end mp_set_int() */
346 /* }}} */
348 /* {{{ mp_set_ulong(mp, z) */
350 mp_err mp_set_ulong(mp_int *mp, unsigned long z)
351 {
352 int ix;
353 mp_err res;
355 ARGCHK(mp != NULL, MP_BADARG);
357 mp_zero(mp);
358 if(z == 0)
359 return MP_OKAY; /* shortcut for zero */
361 if (sizeof z <= sizeof(mp_digit)) {
362 DIGIT(mp,0) = z;
363 } else {
364 for (ix = sizeof(long) - 1; ix >= 0; ix--) {
365 if ((res = s_mp_mul_d(mp, (UCHAR_MAX + 1))) != MP_OKAY)
366 return res;
368 res = s_mp_add_d(mp, (mp_digit)((z >> (ix * CHAR_BIT)) & UCHAR_MAX));
369 if (res != MP_OKAY)
370 return res;
371 }
372 }
373 return MP_OKAY;
374 } /* end mp_set_ulong() */
376 /* }}} */
378 /*------------------------------------------------------------------------*/
379 /* {{{ Digit arithmetic */
381 /* {{{ mp_add_d(a, d, b) */
383 /*
384 mp_add_d(a, d, b)
386 Compute the sum b = a + d, for a single digit d. Respects the sign of
387 its primary addend (single digits are unsigned anyway).
388 */
390 mp_err mp_add_d(const mp_int *a, mp_digit d, mp_int *b)
391 {
392 mp_int tmp;
393 mp_err res;
395 ARGCHK(a != NULL && b != NULL, MP_BADARG);
397 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
398 return res;
400 if(SIGN(&tmp) == ZPOS) {
401 if((res = s_mp_add_d(&tmp, d)) != MP_OKAY)
402 goto CLEANUP;
403 } else if(s_mp_cmp_d(&tmp, d) >= 0) {
404 if((res = s_mp_sub_d(&tmp, d)) != MP_OKAY)
405 goto CLEANUP;
406 } else {
407 mp_neg(&tmp, &tmp);
409 DIGIT(&tmp, 0) = d - DIGIT(&tmp, 0);
410 }
412 if(s_mp_cmp_d(&tmp, 0) == 0)
413 SIGN(&tmp) = ZPOS;
415 s_mp_exch(&tmp, b);
417 CLEANUP:
418 mp_clear(&tmp);
419 return res;
421 } /* end mp_add_d() */
423 /* }}} */
425 /* {{{ mp_sub_d(a, d, b) */
427 /*
428 mp_sub_d(a, d, b)
430 Compute the difference b = a - d, for a single digit d. Respects the
431 sign of its subtrahend (single digits are unsigned anyway).
432 */
434 mp_err mp_sub_d(const mp_int *a, mp_digit d, mp_int *b)
435 {
436 mp_int tmp;
437 mp_err res;
439 ARGCHK(a != NULL && b != NULL, MP_BADARG);
441 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
442 return res;
444 if(SIGN(&tmp) == NEG) {
445 if((res = s_mp_add_d(&tmp, d)) != MP_OKAY)
446 goto CLEANUP;
447 } else if(s_mp_cmp_d(&tmp, d) >= 0) {
448 if((res = s_mp_sub_d(&tmp, d)) != MP_OKAY)
449 goto CLEANUP;
450 } else {
451 mp_neg(&tmp, &tmp);
453 DIGIT(&tmp, 0) = d - DIGIT(&tmp, 0);
454 SIGN(&tmp) = NEG;
455 }
457 if(s_mp_cmp_d(&tmp, 0) == 0)
458 SIGN(&tmp) = ZPOS;
460 s_mp_exch(&tmp, b);
462 CLEANUP:
463 mp_clear(&tmp);
464 return res;
466 } /* end mp_sub_d() */
468 /* }}} */
470 /* {{{ mp_mul_d(a, d, b) */
472 /*
473 mp_mul_d(a, d, b)
475 Compute the product b = a * d, for a single digit d. Respects the sign
476 of its multiplicand (single digits are unsigned anyway)
477 */
479 mp_err mp_mul_d(const mp_int *a, mp_digit d, mp_int *b)
480 {
481 mp_err res;
483 ARGCHK(a != NULL && b != NULL, MP_BADARG);
485 if(d == 0) {
486 mp_zero(b);
487 return MP_OKAY;
488 }
490 if((res = mp_copy(a, b)) != MP_OKAY)
491 return res;
493 res = s_mp_mul_d(b, d);
495 return res;
497 } /* end mp_mul_d() */
499 /* }}} */
501 /* {{{ mp_mul_2(a, c) */
503 mp_err mp_mul_2(const mp_int *a, mp_int *c)
504 {
505 mp_err res;
507 ARGCHK(a != NULL && c != NULL, MP_BADARG);
509 if((res = mp_copy(a, c)) != MP_OKAY)
510 return res;
512 return s_mp_mul_2(c);
514 } /* end mp_mul_2() */
516 /* }}} */
518 /* {{{ mp_div_d(a, d, q, r) */
520 /*
521 mp_div_d(a, d, q, r)
523 Compute the quotient q = a / d and remainder r = a mod d, for a
524 single digit d. Respects the sign of its divisor (single digits are
525 unsigned anyway).
526 */
528 mp_err mp_div_d(const mp_int *a, mp_digit d, mp_int *q, mp_digit *r)
529 {
530 mp_err res;
531 mp_int qp;
532 mp_digit rem;
533 int pow;
535 ARGCHK(a != NULL, MP_BADARG);
537 if(d == 0)
538 return MP_RANGE;
540 /* Shortcut for powers of two ... */
541 if((pow = s_mp_ispow2d(d)) >= 0) {
542 mp_digit mask;
544 mask = ((mp_digit)1 << pow) - 1;
545 rem = DIGIT(a, 0) & mask;
547 if(q) {
548 mp_copy(a, q);
549 s_mp_div_2d(q, pow);
550 }
552 if(r)
553 *r = rem;
555 return MP_OKAY;
556 }
558 if((res = mp_init_copy(&qp, a)) != MP_OKAY)
559 return res;
561 res = s_mp_div_d(&qp, d, &rem);
563 if(s_mp_cmp_d(&qp, 0) == 0)
564 SIGN(q) = ZPOS;
566 if(r)
567 *r = rem;
569 if(q)
570 s_mp_exch(&qp, q);
572 mp_clear(&qp);
573 return res;
575 } /* end mp_div_d() */
577 /* }}} */
579 /* {{{ mp_div_2(a, c) */
581 /*
582 mp_div_2(a, c)
584 Compute c = a / 2, disregarding the remainder.
585 */
587 mp_err mp_div_2(const mp_int *a, mp_int *c)
588 {
589 mp_err res;
591 ARGCHK(a != NULL && c != NULL, MP_BADARG);
593 if((res = mp_copy(a, c)) != MP_OKAY)
594 return res;
596 s_mp_div_2(c);
598 return MP_OKAY;
600 } /* end mp_div_2() */
602 /* }}} */
604 /* {{{ mp_expt_d(a, d, b) */
606 mp_err mp_expt_d(const mp_int *a, mp_digit d, mp_int *c)
607 {
608 mp_int s, x;
609 mp_err res;
611 ARGCHK(a != NULL && c != NULL, MP_BADARG);
613 if((res = mp_init(&s)) != MP_OKAY)
614 return res;
615 if((res = mp_init_copy(&x, a)) != MP_OKAY)
616 goto X;
618 DIGIT(&s, 0) = 1;
620 while(d != 0) {
621 if(d & 1) {
622 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
623 goto CLEANUP;
624 }
626 d /= 2;
628 if((res = s_mp_sqr(&x)) != MP_OKAY)
629 goto CLEANUP;
630 }
632 s_mp_exch(&s, c);
634 CLEANUP:
635 mp_clear(&x);
636 X:
637 mp_clear(&s);
639 return res;
641 } /* end mp_expt_d() */
643 /* }}} */
645 /* }}} */
647 /*------------------------------------------------------------------------*/
648 /* {{{ Full arithmetic */
650 /* {{{ mp_abs(a, b) */
652 /*
653 mp_abs(a, b)
655 Compute b = |a|. 'a' and 'b' may be identical.
656 */
658 mp_err mp_abs(const mp_int *a, mp_int *b)
659 {
660 mp_err res;
662 ARGCHK(a != NULL && b != NULL, MP_BADARG);
664 if((res = mp_copy(a, b)) != MP_OKAY)
665 return res;
667 SIGN(b) = ZPOS;
669 return MP_OKAY;
671 } /* end mp_abs() */
673 /* }}} */
675 /* {{{ mp_neg(a, b) */
677 /*
678 mp_neg(a, b)
680 Compute b = -a. 'a' and 'b' may be identical.
681 */
683 mp_err mp_neg(const mp_int *a, mp_int *b)
684 {
685 mp_err res;
687 ARGCHK(a != NULL && b != NULL, MP_BADARG);
689 if((res = mp_copy(a, b)) != MP_OKAY)
690 return res;
692 if(s_mp_cmp_d(b, 0) == MP_EQ)
693 SIGN(b) = ZPOS;
694 else
695 SIGN(b) = (SIGN(b) == NEG) ? ZPOS : NEG;
697 return MP_OKAY;
699 } /* end mp_neg() */
701 /* }}} */
703 /* {{{ mp_add(a, b, c) */
705 /*
706 mp_add(a, b, c)
708 Compute c = a + b. All parameters may be identical.
709 */
711 mp_err mp_add(const mp_int *a, const mp_int *b, mp_int *c)
712 {
713 mp_err res;
715 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
717 if(SIGN(a) == SIGN(b)) { /* same sign: add values, keep sign */
718 MP_CHECKOK( s_mp_add_3arg(a, b, c) );
719 } else if(s_mp_cmp(a, b) >= 0) { /* different sign: |a| >= |b| */
720 MP_CHECKOK( s_mp_sub_3arg(a, b, c) );
721 } else { /* different sign: |a| < |b| */
722 MP_CHECKOK( s_mp_sub_3arg(b, a, c) );
723 }
725 if (s_mp_cmp_d(c, 0) == MP_EQ)
726 SIGN(c) = ZPOS;
728 CLEANUP:
729 return res;
731 } /* end mp_add() */
733 /* }}} */
735 /* {{{ mp_sub(a, b, c) */
737 /*
738 mp_sub(a, b, c)
740 Compute c = a - b. All parameters may be identical.
741 */
743 mp_err mp_sub(const mp_int *a, const mp_int *b, mp_int *c)
744 {
745 mp_err res;
746 int magDiff;
748 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
750 if (a == b) {
751 mp_zero(c);
752 return MP_OKAY;
753 }
755 if (MP_SIGN(a) != MP_SIGN(b)) {
756 MP_CHECKOK( s_mp_add_3arg(a, b, c) );
757 } else if (!(magDiff = s_mp_cmp(a, b))) {
758 mp_zero(c);
759 res = MP_OKAY;
760 } else if (magDiff > 0) {
761 MP_CHECKOK( s_mp_sub_3arg(a, b, c) );
762 } else {
763 MP_CHECKOK( s_mp_sub_3arg(b, a, c) );
764 MP_SIGN(c) = !MP_SIGN(a);
765 }
767 if (s_mp_cmp_d(c, 0) == MP_EQ)
768 MP_SIGN(c) = MP_ZPOS;
770 CLEANUP:
771 return res;
773 } /* end mp_sub() */
775 /* }}} */
777 /* {{{ mp_mul(a, b, c) */
779 /*
780 mp_mul(a, b, c)
782 Compute c = a * b. All parameters may be identical.
783 */
784 mp_err mp_mul(const mp_int *a, const mp_int *b, mp_int * c)
785 {
786 mp_digit *pb;
787 mp_int tmp;
788 mp_err res;
789 mp_size ib;
790 mp_size useda, usedb;
792 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
794 if (a == c) {
795 if ((res = mp_init_copy(&tmp, a)) != MP_OKAY)
796 return res;
797 if (a == b)
798 b = &tmp;
799 a = &tmp;
800 } else if (b == c) {
801 if ((res = mp_init_copy(&tmp, b)) != MP_OKAY)
802 return res;
803 b = &tmp;
804 } else {
805 MP_DIGITS(&tmp) = 0;
806 }
808 if (MP_USED(a) < MP_USED(b)) {
809 const mp_int *xch = b; /* switch a and b, to do fewer outer loops */
810 b = a;
811 a = xch;
812 }
814 MP_USED(c) = 1; MP_DIGIT(c, 0) = 0;
815 if((res = s_mp_pad(c, USED(a) + USED(b))) != MP_OKAY)
816 goto CLEANUP;
818 #ifdef NSS_USE_COMBA
819 if ((MP_USED(a) == MP_USED(b)) && IS_POWER_OF_2(MP_USED(b))) {
820 if (MP_USED(a) == 4) {
821 s_mp_mul_comba_4(a, b, c);
822 goto CLEANUP;
823 }
824 if (MP_USED(a) == 8) {
825 s_mp_mul_comba_8(a, b, c);
826 goto CLEANUP;
827 }
828 if (MP_USED(a) == 16) {
829 s_mp_mul_comba_16(a, b, c);
830 goto CLEANUP;
831 }
832 if (MP_USED(a) == 32) {
833 s_mp_mul_comba_32(a, b, c);
834 goto CLEANUP;
835 }
836 }
837 #endif
839 pb = MP_DIGITS(b);
840 s_mpv_mul_d(MP_DIGITS(a), MP_USED(a), *pb++, MP_DIGITS(c));
842 /* Outer loop: Digits of b */
843 useda = MP_USED(a);
844 usedb = MP_USED(b);
845 for (ib = 1; ib < usedb; ib++) {
846 mp_digit b_i = *pb++;
848 /* Inner product: Digits of a */
849 if (b_i)
850 s_mpv_mul_d_add(MP_DIGITS(a), useda, b_i, MP_DIGITS(c) + ib);
851 else
852 MP_DIGIT(c, ib + useda) = b_i;
853 }
855 s_mp_clamp(c);
857 if(SIGN(a) == SIGN(b) || s_mp_cmp_d(c, 0) == MP_EQ)
858 SIGN(c) = ZPOS;
859 else
860 SIGN(c) = NEG;
862 CLEANUP:
863 mp_clear(&tmp);
864 return res;
865 } /* end mp_mul() */
867 /* }}} */
869 /* {{{ mp_sqr(a, sqr) */
871 #if MP_SQUARE
872 /*
873 Computes the square of a. This can be done more
874 efficiently than a general multiplication, because many of the
875 computation steps are redundant when squaring. The inner product
876 step is a bit more complicated, but we save a fair number of
877 iterations of the multiplication loop.
878 */
880 /* sqr = a^2; Caller provides both a and tmp; */
881 mp_err mp_sqr(const mp_int *a, mp_int *sqr)
882 {
883 mp_digit *pa;
884 mp_digit d;
885 mp_err res;
886 mp_size ix;
887 mp_int tmp;
888 int count;
890 ARGCHK(a != NULL && sqr != NULL, MP_BADARG);
892 if (a == sqr) {
893 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
894 return res;
895 a = &tmp;
896 } else {
897 DIGITS(&tmp) = 0;
898 res = MP_OKAY;
899 }
901 ix = 2 * MP_USED(a);
902 if (ix > MP_ALLOC(sqr)) {
903 MP_USED(sqr) = 1;
904 MP_CHECKOK( s_mp_grow(sqr, ix) );
905 }
906 MP_USED(sqr) = ix;
907 MP_DIGIT(sqr, 0) = 0;
909 #ifdef NSS_USE_COMBA
910 if (IS_POWER_OF_2(MP_USED(a))) {
911 if (MP_USED(a) == 4) {
912 s_mp_sqr_comba_4(a, sqr);
913 goto CLEANUP;
914 }
915 if (MP_USED(a) == 8) {
916 s_mp_sqr_comba_8(a, sqr);
917 goto CLEANUP;
918 }
919 if (MP_USED(a) == 16) {
920 s_mp_sqr_comba_16(a, sqr);
921 goto CLEANUP;
922 }
923 if (MP_USED(a) == 32) {
924 s_mp_sqr_comba_32(a, sqr);
925 goto CLEANUP;
926 }
927 }
928 #endif
930 pa = MP_DIGITS(a);
931 count = MP_USED(a) - 1;
932 if (count > 0) {
933 d = *pa++;
934 s_mpv_mul_d(pa, count, d, MP_DIGITS(sqr) + 1);
935 for (ix = 3; --count > 0; ix += 2) {
936 d = *pa++;
937 s_mpv_mul_d_add(pa, count, d, MP_DIGITS(sqr) + ix);
938 } /* for(ix ...) */
939 MP_DIGIT(sqr, MP_USED(sqr)-1) = 0; /* above loop stopped short of this. */
941 /* now sqr *= 2 */
942 s_mp_mul_2(sqr);
943 } else {
944 MP_DIGIT(sqr, 1) = 0;
945 }
947 /* now add the squares of the digits of a to sqr. */
948 s_mpv_sqr_add_prop(MP_DIGITS(a), MP_USED(a), MP_DIGITS(sqr));
950 SIGN(sqr) = ZPOS;
951 s_mp_clamp(sqr);
953 CLEANUP:
954 mp_clear(&tmp);
955 return res;
957 } /* end mp_sqr() */
958 #endif
960 /* }}} */
962 /* {{{ mp_div(a, b, q, r) */
964 /*
965 mp_div(a, b, q, r)
967 Compute q = a / b and r = a mod b. Input parameters may be re-used
968 as output parameters. If q or r is NULL, that portion of the
969 computation will be discarded (although it will still be computed)
970 */
971 mp_err mp_div(const mp_int *a, const mp_int *b, mp_int *q, mp_int *r)
972 {
973 mp_err res;
974 mp_int *pQ, *pR;
975 mp_int qtmp, rtmp, btmp;
976 int cmp;
977 mp_sign signA;
978 mp_sign signB;
980 ARGCHK(a != NULL && b != NULL, MP_BADARG);
982 signA = MP_SIGN(a);
983 signB = MP_SIGN(b);
985 if(mp_cmp_z(b) == MP_EQ)
986 return MP_RANGE;
988 DIGITS(&qtmp) = 0;
989 DIGITS(&rtmp) = 0;
990 DIGITS(&btmp) = 0;
992 /* Set up some temporaries... */
993 if (!r || r == a || r == b) {
994 MP_CHECKOK( mp_init_copy(&rtmp, a) );
995 pR = &rtmp;
996 } else {
997 MP_CHECKOK( mp_copy(a, r) );
998 pR = r;
999 }
1001 if (!q || q == a || q == b) {
1002 MP_CHECKOK( mp_init_size(&qtmp, MP_USED(a)) );
1003 pQ = &qtmp;
1004 } else {
1005 MP_CHECKOK( s_mp_pad(q, MP_USED(a)) );
1006 pQ = q;
1007 mp_zero(pQ);
1008 }
1010 /*
1011 If |a| <= |b|, we can compute the solution without division;
1012 otherwise, we actually do the work required.
1013 */
1014 if ((cmp = s_mp_cmp(a, b)) <= 0) {
1015 if (cmp) {
1016 /* r was set to a above. */
1017 mp_zero(pQ);
1018 } else {
1019 mp_set(pQ, 1);
1020 mp_zero(pR);
1021 }
1022 } else {
1023 MP_CHECKOK( mp_init_copy(&btmp, b) );
1024 MP_CHECKOK( s_mp_div(pR, &btmp, pQ) );
1025 }
1027 /* Compute the signs for the output */
1028 MP_SIGN(pR) = signA; /* Sr = Sa */
1029 /* Sq = ZPOS if Sa == Sb */ /* Sq = NEG if Sa != Sb */
1030 MP_SIGN(pQ) = (signA == signB) ? ZPOS : NEG;
1032 if(s_mp_cmp_d(pQ, 0) == MP_EQ)
1033 SIGN(pQ) = ZPOS;
1034 if(s_mp_cmp_d(pR, 0) == MP_EQ)
1035 SIGN(pR) = ZPOS;
1037 /* Copy output, if it is needed */
1038 if(q && q != pQ)
1039 s_mp_exch(pQ, q);
1041 if(r && r != pR)
1042 s_mp_exch(pR, r);
1044 CLEANUP:
1045 mp_clear(&btmp);
1046 mp_clear(&rtmp);
1047 mp_clear(&qtmp);
1049 return res;
1051 } /* end mp_div() */
1053 /* }}} */
1055 /* {{{ mp_div_2d(a, d, q, r) */
1057 mp_err mp_div_2d(const mp_int *a, mp_digit d, mp_int *q, mp_int *r)
1058 {
1059 mp_err res;
1061 ARGCHK(a != NULL, MP_BADARG);
1063 if(q) {
1064 if((res = mp_copy(a, q)) != MP_OKAY)
1065 return res;
1066 }
1067 if(r) {
1068 if((res = mp_copy(a, r)) != MP_OKAY)
1069 return res;
1070 }
1071 if(q) {
1072 s_mp_div_2d(q, d);
1073 }
1074 if(r) {
1075 s_mp_mod_2d(r, d);
1076 }
1078 return MP_OKAY;
1080 } /* end mp_div_2d() */
1082 /* }}} */
1084 /* {{{ mp_expt(a, b, c) */
1086 /*
1087 mp_expt(a, b, c)
1089 Compute c = a ** b, that is, raise a to the b power. Uses a
1090 standard iterative square-and-multiply technique.
1091 */
1093 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c)
1094 {
1095 mp_int s, x;
1096 mp_err res;
1097 mp_digit d;
1098 int dig, bit;
1100 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1102 if(mp_cmp_z(b) < 0)
1103 return MP_RANGE;
1105 if((res = mp_init(&s)) != MP_OKAY)
1106 return res;
1108 mp_set(&s, 1);
1110 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1111 goto X;
1113 /* Loop over low-order digits in ascending order */
1114 for(dig = 0; dig < (USED(b) - 1); dig++) {
1115 d = DIGIT(b, dig);
1117 /* Loop over bits of each non-maximal digit */
1118 for(bit = 0; bit < DIGIT_BIT; bit++) {
1119 if(d & 1) {
1120 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1121 goto CLEANUP;
1122 }
1124 d >>= 1;
1126 if((res = s_mp_sqr(&x)) != MP_OKAY)
1127 goto CLEANUP;
1128 }
1129 }
1131 /* Consider now the last digit... */
1132 d = DIGIT(b, dig);
1134 while(d) {
1135 if(d & 1) {
1136 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1137 goto CLEANUP;
1138 }
1140 d >>= 1;
1142 if((res = s_mp_sqr(&x)) != MP_OKAY)
1143 goto CLEANUP;
1144 }
1146 if(mp_iseven(b))
1147 SIGN(&s) = SIGN(a);
1149 res = mp_copy(&s, c);
1151 CLEANUP:
1152 mp_clear(&x);
1153 X:
1154 mp_clear(&s);
1156 return res;
1158 } /* end mp_expt() */
1160 /* }}} */
1162 /* {{{ mp_2expt(a, k) */
1164 /* Compute a = 2^k */
1166 mp_err mp_2expt(mp_int *a, mp_digit k)
1167 {
1168 ARGCHK(a != NULL, MP_BADARG);
1170 return s_mp_2expt(a, k);
1172 } /* end mp_2expt() */
1174 /* }}} */
1176 /* {{{ mp_mod(a, m, c) */
1178 /*
1179 mp_mod(a, m, c)
1181 Compute c = a (mod m). Result will always be 0 <= c < m.
1182 */
1184 mp_err mp_mod(const mp_int *a, const mp_int *m, mp_int *c)
1185 {
1186 mp_err res;
1187 int mag;
1189 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
1191 if(SIGN(m) == NEG)
1192 return MP_RANGE;
1194 /*
1195 If |a| > m, we need to divide to get the remainder and take the
1196 absolute value.
1198 If |a| < m, we don't need to do any division, just copy and adjust
1199 the sign (if a is negative).
1201 If |a| == m, we can simply set the result to zero.
1203 This order is intended to minimize the average path length of the
1204 comparison chain on common workloads -- the most frequent cases are
1205 that |a| != m, so we do those first.
1206 */
1207 if((mag = s_mp_cmp(a, m)) > 0) {
1208 if((res = mp_div(a, m, NULL, c)) != MP_OKAY)
1209 return res;
1211 if(SIGN(c) == NEG) {
1212 if((res = mp_add(c, m, c)) != MP_OKAY)
1213 return res;
1214 }
1216 } else if(mag < 0) {
1217 if((res = mp_copy(a, c)) != MP_OKAY)
1218 return res;
1220 if(mp_cmp_z(a) < 0) {
1221 if((res = mp_add(c, m, c)) != MP_OKAY)
1222 return res;
1224 }
1226 } else {
1227 mp_zero(c);
1229 }
1231 return MP_OKAY;
1233 } /* end mp_mod() */
1235 /* }}} */
1237 /* {{{ mp_mod_d(a, d, c) */
1239 /*
1240 mp_mod_d(a, d, c)
1242 Compute c = a (mod d). Result will always be 0 <= c < d
1243 */
1244 mp_err mp_mod_d(const mp_int *a, mp_digit d, mp_digit *c)
1245 {
1246 mp_err res;
1247 mp_digit rem;
1249 ARGCHK(a != NULL && c != NULL, MP_BADARG);
1251 if(s_mp_cmp_d(a, d) > 0) {
1252 if((res = mp_div_d(a, d, NULL, &rem)) != MP_OKAY)
1253 return res;
1255 } else {
1256 if(SIGN(a) == NEG)
1257 rem = d - DIGIT(a, 0);
1258 else
1259 rem = DIGIT(a, 0);
1260 }
1262 if(c)
1263 *c = rem;
1265 return MP_OKAY;
1267 } /* end mp_mod_d() */
1269 /* }}} */
1271 /* {{{ mp_sqrt(a, b) */
1273 /*
1274 mp_sqrt(a, b)
1276 Compute the integer square root of a, and store the result in b.
1277 Uses an integer-arithmetic version of Newton's iterative linear
1278 approximation technique to determine this value; the result has the
1279 following two properties:
1281 b^2 <= a
1282 (b+1)^2 >= a
1284 It is a range error to pass a negative value.
1285 */
1286 mp_err mp_sqrt(const mp_int *a, mp_int *b)
1287 {
1288 mp_int x, t;
1289 mp_err res;
1290 mp_size used;
1292 ARGCHK(a != NULL && b != NULL, MP_BADARG);
1294 /* Cannot take square root of a negative value */
1295 if(SIGN(a) == NEG)
1296 return MP_RANGE;
1298 /* Special cases for zero and one, trivial */
1299 if(mp_cmp_d(a, 1) <= 0)
1300 return mp_copy(a, b);
1302 /* Initialize the temporaries we'll use below */
1303 if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
1304 return res;
1306 /* Compute an initial guess for the iteration as a itself */
1307 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1308 goto X;
1310 used = MP_USED(&x);
1311 if (used > 1) {
1312 s_mp_rshd(&x, used / 2);
1313 }
1315 for(;;) {
1316 /* t = (x * x) - a */
1317 mp_copy(&x, &t); /* can't fail, t is big enough for original x */
1318 if((res = mp_sqr(&t, &t)) != MP_OKAY ||
1319 (res = mp_sub(&t, a, &t)) != MP_OKAY)
1320 goto CLEANUP;
1322 /* t = t / 2x */
1323 s_mp_mul_2(&x);
1324 if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)
1325 goto CLEANUP;
1326 s_mp_div_2(&x);
1328 /* Terminate the loop, if the quotient is zero */
1329 if(mp_cmp_z(&t) == MP_EQ)
1330 break;
1332 /* x = x - t */
1333 if((res = mp_sub(&x, &t, &x)) != MP_OKAY)
1334 goto CLEANUP;
1336 }
1338 /* Copy result to output parameter */
1339 mp_sub_d(&x, 1, &x);
1340 s_mp_exch(&x, b);
1342 CLEANUP:
1343 mp_clear(&x);
1344 X:
1345 mp_clear(&t);
1347 return res;
1349 } /* end mp_sqrt() */
1351 /* }}} */
1353 /* }}} */
1355 /*------------------------------------------------------------------------*/
1356 /* {{{ Modular arithmetic */
1358 #if MP_MODARITH
1359 /* {{{ mp_addmod(a, b, m, c) */
1361 /*
1362 mp_addmod(a, b, m, c)
1364 Compute c = (a + b) mod m
1365 */
1367 mp_err mp_addmod(const mp_int *a, const mp_int *b, const mp_int *m, mp_int *c)
1368 {
1369 mp_err res;
1371 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1373 if((res = mp_add(a, b, c)) != MP_OKAY)
1374 return res;
1375 if((res = mp_mod(c, m, c)) != MP_OKAY)
1376 return res;
1378 return MP_OKAY;
1380 }
1382 /* }}} */
1384 /* {{{ mp_submod(a, b, m, c) */
1386 /*
1387 mp_submod(a, b, m, c)
1389 Compute c = (a - b) mod m
1390 */
1392 mp_err mp_submod(const mp_int *a, const mp_int *b, const mp_int *m, mp_int *c)
1393 {
1394 mp_err res;
1396 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1398 if((res = mp_sub(a, b, c)) != MP_OKAY)
1399 return res;
1400 if((res = mp_mod(c, m, c)) != MP_OKAY)
1401 return res;
1403 return MP_OKAY;
1405 }
1407 /* }}} */
1409 /* {{{ mp_mulmod(a, b, m, c) */
1411 /*
1412 mp_mulmod(a, b, m, c)
1414 Compute c = (a * b) mod m
1415 */
1417 mp_err mp_mulmod(const mp_int *a, const mp_int *b, const mp_int *m, mp_int *c)
1418 {
1419 mp_err res;
1421 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1423 if((res = mp_mul(a, b, c)) != MP_OKAY)
1424 return res;
1425 if((res = mp_mod(c, m, c)) != MP_OKAY)
1426 return res;
1428 return MP_OKAY;
1430 }
1432 /* }}} */
1434 /* {{{ mp_sqrmod(a, m, c) */
1436 #if MP_SQUARE
1437 mp_err mp_sqrmod(const mp_int *a, const mp_int *m, mp_int *c)
1438 {
1439 mp_err res;
1441 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
1443 if((res = mp_sqr(a, c)) != MP_OKAY)
1444 return res;
1445 if((res = mp_mod(c, m, c)) != MP_OKAY)
1446 return res;
1448 return MP_OKAY;
1450 } /* end mp_sqrmod() */
1451 #endif
1453 /* }}} */
1455 /* {{{ s_mp_exptmod(a, b, m, c) */
1457 /*
1458 s_mp_exptmod(a, b, m, c)
1460 Compute c = (a ** b) mod m. Uses a standard square-and-multiply
1461 method with modular reductions at each step. (This is basically the
1462 same code as mp_expt(), except for the addition of the reductions)
1464 The modular reductions are done using Barrett's algorithm (see
1465 s_mp_reduce() below for details)
1466 */
1468 mp_err s_mp_exptmod(const mp_int *a, const mp_int *b, const mp_int *m, mp_int *c)
1469 {
1470 mp_int s, x, mu;
1471 mp_err res;
1472 mp_digit d;
1473 int dig, bit;
1475 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1477 if(mp_cmp_z(b) < 0 || mp_cmp_z(m) <= 0)
1478 return MP_RANGE;
1480 if((res = mp_init(&s)) != MP_OKAY)
1481 return res;
1482 if((res = mp_init_copy(&x, a)) != MP_OKAY ||
1483 (res = mp_mod(&x, m, &x)) != MP_OKAY)
1484 goto X;
1485 if((res = mp_init(&mu)) != MP_OKAY)
1486 goto MU;
1488 mp_set(&s, 1);
1490 /* mu = b^2k / m */
1491 s_mp_add_d(&mu, 1);
1492 s_mp_lshd(&mu, 2 * USED(m));
1493 if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY)
1494 goto CLEANUP;
1496 /* Loop over digits of b in ascending order, except highest order */
1497 for(dig = 0; dig < (USED(b) - 1); dig++) {
1498 d = DIGIT(b, dig);
1500 /* Loop over the bits of the lower-order digits */
1501 for(bit = 0; bit < DIGIT_BIT; bit++) {
1502 if(d & 1) {
1503 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1504 goto CLEANUP;
1505 if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)
1506 goto CLEANUP;
1507 }
1509 d >>= 1;
1511 if((res = s_mp_sqr(&x)) != MP_OKAY)
1512 goto CLEANUP;
1513 if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)
1514 goto CLEANUP;
1515 }
1516 }
1518 /* Now do the last digit... */
1519 d = DIGIT(b, dig);
1521 while(d) {
1522 if(d & 1) {
1523 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1524 goto CLEANUP;
1525 if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)
1526 goto CLEANUP;
1527 }
1529 d >>= 1;
1531 if((res = s_mp_sqr(&x)) != MP_OKAY)
1532 goto CLEANUP;
1533 if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)
1534 goto CLEANUP;
1535 }
1537 s_mp_exch(&s, c);
1539 CLEANUP:
1540 mp_clear(&mu);
1541 MU:
1542 mp_clear(&x);
1543 X:
1544 mp_clear(&s);
1546 return res;
1548 } /* end s_mp_exptmod() */
1550 /* }}} */
1552 /* {{{ mp_exptmod_d(a, d, m, c) */
1554 mp_err mp_exptmod_d(const mp_int *a, mp_digit d, const mp_int *m, mp_int *c)
1555 {
1556 mp_int s, x;
1557 mp_err res;
1559 ARGCHK(a != NULL && c != NULL, MP_BADARG);
1561 if((res = mp_init(&s)) != MP_OKAY)
1562 return res;
1563 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1564 goto X;
1566 mp_set(&s, 1);
1568 while(d != 0) {
1569 if(d & 1) {
1570 if((res = s_mp_mul(&s, &x)) != MP_OKAY ||
1571 (res = mp_mod(&s, m, &s)) != MP_OKAY)
1572 goto CLEANUP;
1573 }
1575 d /= 2;
1577 if((res = s_mp_sqr(&x)) != MP_OKAY ||
1578 (res = mp_mod(&x, m, &x)) != MP_OKAY)
1579 goto CLEANUP;
1580 }
1582 s_mp_exch(&s, c);
1584 CLEANUP:
1585 mp_clear(&x);
1586 X:
1587 mp_clear(&s);
1589 return res;
1591 } /* end mp_exptmod_d() */
1593 /* }}} */
1594 #endif /* if MP_MODARITH */
1596 /* }}} */
1598 /*------------------------------------------------------------------------*/
1599 /* {{{ Comparison functions */
1601 /* {{{ mp_cmp_z(a) */
1603 /*
1604 mp_cmp_z(a)
1606 Compare a <=> 0. Returns <0 if a<0, 0 if a=0, >0 if a>0.
1607 */
1609 int mp_cmp_z(const mp_int *a)
1610 {
1611 if(SIGN(a) == NEG)
1612 return MP_LT;
1613 else if(USED(a) == 1 && DIGIT(a, 0) == 0)
1614 return MP_EQ;
1615 else
1616 return MP_GT;
1618 } /* end mp_cmp_z() */
1620 /* }}} */
1622 /* {{{ mp_cmp_d(a, d) */
1624 /*
1625 mp_cmp_d(a, d)
1627 Compare a <=> d. Returns <0 if a<d, 0 if a=d, >0 if a>d
1628 */
1630 int mp_cmp_d(const mp_int *a, mp_digit d)
1631 {
1632 ARGCHK(a != NULL, MP_EQ);
1634 if(SIGN(a) == NEG)
1635 return MP_LT;
1637 return s_mp_cmp_d(a, d);
1639 } /* end mp_cmp_d() */
1641 /* }}} */
1643 /* {{{ mp_cmp(a, b) */
1645 int mp_cmp(const mp_int *a, const mp_int *b)
1646 {
1647 ARGCHK(a != NULL && b != NULL, MP_EQ);
1649 if(SIGN(a) == SIGN(b)) {
1650 int mag;
1652 if((mag = s_mp_cmp(a, b)) == MP_EQ)
1653 return MP_EQ;
1655 if(SIGN(a) == ZPOS)
1656 return mag;
1657 else
1658 return -mag;
1660 } else if(SIGN(a) == ZPOS) {
1661 return MP_GT;
1662 } else {
1663 return MP_LT;
1664 }
1666 } /* end mp_cmp() */
1668 /* }}} */
1670 /* {{{ mp_cmp_mag(a, b) */
1672 /*
1673 mp_cmp_mag(a, b)
1675 Compares |a| <=> |b|, and returns an appropriate comparison result
1676 */
1678 int mp_cmp_mag(mp_int *a, mp_int *b)
1679 {
1680 ARGCHK(a != NULL && b != NULL, MP_EQ);
1682 return s_mp_cmp(a, b);
1684 } /* end mp_cmp_mag() */
1686 /* }}} */
1688 /* {{{ mp_cmp_int(a, z) */
1690 /*
1691 This just converts z to an mp_int, and uses the existing comparison
1692 routines. This is sort of inefficient, but it's not clear to me how
1693 frequently this wil get used anyway. For small positive constants,
1694 you can always use mp_cmp_d(), and for zero, there is mp_cmp_z().
1695 */
1696 int mp_cmp_int(const mp_int *a, long z)
1697 {
1698 mp_int tmp;
1699 int out;
1701 ARGCHK(a != NULL, MP_EQ);
1703 mp_init(&tmp); mp_set_int(&tmp, z);
1704 out = mp_cmp(a, &tmp);
1705 mp_clear(&tmp);
1707 return out;
1709 } /* end mp_cmp_int() */
1711 /* }}} */
1713 /* {{{ mp_isodd(a) */
1715 /*
1716 mp_isodd(a)
1718 Returns a true (non-zero) value if a is odd, false (zero) otherwise.
1719 */
1720 int mp_isodd(const mp_int *a)
1721 {
1722 ARGCHK(a != NULL, 0);
1724 return (int)(DIGIT(a, 0) & 1);
1726 } /* end mp_isodd() */
1728 /* }}} */
1730 /* {{{ mp_iseven(a) */
1732 int mp_iseven(const mp_int *a)
1733 {
1734 return !mp_isodd(a);
1736 } /* end mp_iseven() */
1738 /* }}} */
1740 /* }}} */
1742 /*------------------------------------------------------------------------*/
1743 /* {{{ Number theoretic functions */
1745 #if MP_NUMTH
1746 /* {{{ mp_gcd(a, b, c) */
1748 /*
1749 Like the old mp_gcd() function, except computes the GCD using the
1750 binary algorithm due to Josef Stein in 1961 (via Knuth).
1751 */
1752 mp_err mp_gcd(mp_int *a, mp_int *b, mp_int *c)
1753 {
1754 mp_err res;
1755 mp_int u, v, t;
1756 mp_size k = 0;
1758 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1760 if(mp_cmp_z(a) == MP_EQ && mp_cmp_z(b) == MP_EQ)
1761 return MP_RANGE;
1762 if(mp_cmp_z(a) == MP_EQ) {
1763 return mp_copy(b, c);
1764 } else if(mp_cmp_z(b) == MP_EQ) {
1765 return mp_copy(a, c);
1766 }
1768 if((res = mp_init(&t)) != MP_OKAY)
1769 return res;
1770 if((res = mp_init_copy(&u, a)) != MP_OKAY)
1771 goto U;
1772 if((res = mp_init_copy(&v, b)) != MP_OKAY)
1773 goto V;
1775 SIGN(&u) = ZPOS;
1776 SIGN(&v) = ZPOS;
1778 /* Divide out common factors of 2 until at least 1 of a, b is even */
1779 while(mp_iseven(&u) && mp_iseven(&v)) {
1780 s_mp_div_2(&u);
1781 s_mp_div_2(&v);
1782 ++k;
1783 }
1785 /* Initialize t */
1786 if(mp_isodd(&u)) {
1787 if((res = mp_copy(&v, &t)) != MP_OKAY)
1788 goto CLEANUP;
1790 /* t = -v */
1791 if(SIGN(&v) == ZPOS)
1792 SIGN(&t) = NEG;
1793 else
1794 SIGN(&t) = ZPOS;
1796 } else {
1797 if((res = mp_copy(&u, &t)) != MP_OKAY)
1798 goto CLEANUP;
1800 }
1802 for(;;) {
1803 while(mp_iseven(&t)) {
1804 s_mp_div_2(&t);
1805 }
1807 if(mp_cmp_z(&t) == MP_GT) {
1808 if((res = mp_copy(&t, &u)) != MP_OKAY)
1809 goto CLEANUP;
1811 } else {
1812 if((res = mp_copy(&t, &v)) != MP_OKAY)
1813 goto CLEANUP;
1815 /* v = -t */
1816 if(SIGN(&t) == ZPOS)
1817 SIGN(&v) = NEG;
1818 else
1819 SIGN(&v) = ZPOS;
1820 }
1822 if((res = mp_sub(&u, &v, &t)) != MP_OKAY)
1823 goto CLEANUP;
1825 if(s_mp_cmp_d(&t, 0) == MP_EQ)
1826 break;
1827 }
1829 s_mp_2expt(&v, k); /* v = 2^k */
1830 res = mp_mul(&u, &v, c); /* c = u * v */
1832 CLEANUP:
1833 mp_clear(&v);
1834 V:
1835 mp_clear(&u);
1836 U:
1837 mp_clear(&t);
1839 return res;
1841 } /* end mp_gcd() */
1843 /* }}} */
1845 /* {{{ mp_lcm(a, b, c) */
1847 /* We compute the least common multiple using the rule:
1849 ab = [a, b](a, b)
1851 ... by computing the product, and dividing out the gcd.
1852 */
1854 mp_err mp_lcm(mp_int *a, mp_int *b, mp_int *c)
1855 {
1856 mp_int gcd, prod;
1857 mp_err res;
1859 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1861 /* Set up temporaries */
1862 if((res = mp_init(&gcd)) != MP_OKAY)
1863 return res;
1864 if((res = mp_init(&prod)) != MP_OKAY)
1865 goto GCD;
1867 if((res = mp_mul(a, b, &prod)) != MP_OKAY)
1868 goto CLEANUP;
1869 if((res = mp_gcd(a, b, &gcd)) != MP_OKAY)
1870 goto CLEANUP;
1872 res = mp_div(&prod, &gcd, c, NULL);
1874 CLEANUP:
1875 mp_clear(&prod);
1876 GCD:
1877 mp_clear(&gcd);
1879 return res;
1881 } /* end mp_lcm() */
1883 /* }}} */
1885 /* {{{ mp_xgcd(a, b, g, x, y) */
1887 /*
1888 mp_xgcd(a, b, g, x, y)
1890 Compute g = (a, b) and values x and y satisfying Bezout's identity
1891 (that is, ax + by = g). This uses the binary extended GCD algorithm
1892 based on the Stein algorithm used for mp_gcd()
1893 See algorithm 14.61 in Handbook of Applied Cryptogrpahy.
1894 */
1896 mp_err mp_xgcd(const mp_int *a, const mp_int *b, mp_int *g, mp_int *x, mp_int *y)
1897 {
1898 mp_int gx, xc, yc, u, v, A, B, C, D;
1899 mp_int *clean[9];
1900 mp_err res;
1901 int last = -1;
1903 if(mp_cmp_z(b) == 0)
1904 return MP_RANGE;
1906 /* Initialize all these variables we need */
1907 MP_CHECKOK( mp_init(&u) );
1908 clean[++last] = &u;
1909 MP_CHECKOK( mp_init(&v) );
1910 clean[++last] = &v;
1911 MP_CHECKOK( mp_init(&gx) );
1912 clean[++last] = &gx;
1913 MP_CHECKOK( mp_init(&A) );
1914 clean[++last] = &A;
1915 MP_CHECKOK( mp_init(&B) );
1916 clean[++last] = &B;
1917 MP_CHECKOK( mp_init(&C) );
1918 clean[++last] = &C;
1919 MP_CHECKOK( mp_init(&D) );
1920 clean[++last] = &D;
1921 MP_CHECKOK( mp_init_copy(&xc, a) );
1922 clean[++last] = &xc;
1923 mp_abs(&xc, &xc);
1924 MP_CHECKOK( mp_init_copy(&yc, b) );
1925 clean[++last] = &yc;
1926 mp_abs(&yc, &yc);
1928 mp_set(&gx, 1);
1930 /* Divide by two until at least one of them is odd */
1931 while(mp_iseven(&xc) && mp_iseven(&yc)) {
1932 mp_size nx = mp_trailing_zeros(&xc);
1933 mp_size ny = mp_trailing_zeros(&yc);
1934 mp_size n = MP_MIN(nx, ny);
1935 s_mp_div_2d(&xc,n);
1936 s_mp_div_2d(&yc,n);
1937 MP_CHECKOK( s_mp_mul_2d(&gx,n) );
1938 }
1940 mp_copy(&xc, &u);
1941 mp_copy(&yc, &v);
1942 mp_set(&A, 1); mp_set(&D, 1);
1944 /* Loop through binary GCD algorithm */
1945 do {
1946 while(mp_iseven(&u)) {
1947 s_mp_div_2(&u);
1949 if(mp_iseven(&A) && mp_iseven(&B)) {
1950 s_mp_div_2(&A); s_mp_div_2(&B);
1951 } else {
1952 MP_CHECKOK( mp_add(&A, &yc, &A) );
1953 s_mp_div_2(&A);
1954 MP_CHECKOK( mp_sub(&B, &xc, &B) );
1955 s_mp_div_2(&B);
1956 }
1957 }
1959 while(mp_iseven(&v)) {
1960 s_mp_div_2(&v);
1962 if(mp_iseven(&C) && mp_iseven(&D)) {
1963 s_mp_div_2(&C); s_mp_div_2(&D);
1964 } else {
1965 MP_CHECKOK( mp_add(&C, &yc, &C) );
1966 s_mp_div_2(&C);
1967 MP_CHECKOK( mp_sub(&D, &xc, &D) );
1968 s_mp_div_2(&D);
1969 }
1970 }
1972 if(mp_cmp(&u, &v) >= 0) {
1973 MP_CHECKOK( mp_sub(&u, &v, &u) );
1974 MP_CHECKOK( mp_sub(&A, &C, &A) );
1975 MP_CHECKOK( mp_sub(&B, &D, &B) );
1976 } else {
1977 MP_CHECKOK( mp_sub(&v, &u, &v) );
1978 MP_CHECKOK( mp_sub(&C, &A, &C) );
1979 MP_CHECKOK( mp_sub(&D, &B, &D) );
1980 }
1981 } while (mp_cmp_z(&u) != 0);
1983 /* copy results to output */
1984 if(x)
1985 MP_CHECKOK( mp_copy(&C, x) );
1987 if(y)
1988 MP_CHECKOK( mp_copy(&D, y) );
1990 if(g)
1991 MP_CHECKOK( mp_mul(&gx, &v, g) );
1993 CLEANUP:
1994 while(last >= 0)
1995 mp_clear(clean[last--]);
1997 return res;
1999 } /* end mp_xgcd() */
2001 /* }}} */
2003 mp_size mp_trailing_zeros(const mp_int *mp)
2004 {
2005 mp_digit d;
2006 mp_size n = 0;
2007 int ix;
2009 if (!mp || !MP_DIGITS(mp) || !mp_cmp_z(mp))
2010 return n;
2012 for (ix = 0; !(d = MP_DIGIT(mp,ix)) && (ix < MP_USED(mp)); ++ix)
2013 n += MP_DIGIT_BIT;
2014 if (!d)
2015 return 0; /* shouldn't happen, but ... */
2016 #if !defined(MP_USE_UINT_DIGIT)
2017 if (!(d & 0xffffffffU)) {
2018 d >>= 32;
2019 n += 32;
2020 }
2021 #endif
2022 if (!(d & 0xffffU)) {
2023 d >>= 16;
2024 n += 16;
2025 }
2026 if (!(d & 0xffU)) {
2027 d >>= 8;
2028 n += 8;
2029 }
2030 if (!(d & 0xfU)) {
2031 d >>= 4;
2032 n += 4;
2033 }
2034 if (!(d & 0x3U)) {
2035 d >>= 2;
2036 n += 2;
2037 }
2038 if (!(d & 0x1U)) {
2039 d >>= 1;
2040 n += 1;
2041 }
2042 #if MP_ARGCHK == 2
2043 assert(0 != (d & 1));
2044 #endif
2045 return n;
2046 }
2048 /* Given a and prime p, computes c and k such that a*c == 2**k (mod p).
2049 ** Returns k (positive) or error (negative).
2050 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)
2051 ** by Richard Schroeppel (a.k.a. Captain Nemo).
2052 */
2053 mp_err s_mp_almost_inverse(const mp_int *a, const mp_int *p, mp_int *c)
2054 {
2055 mp_err res;
2056 mp_err k = 0;
2057 mp_int d, f, g;
2059 ARGCHK(a && p && c, MP_BADARG);
2061 MP_DIGITS(&d) = 0;
2062 MP_DIGITS(&f) = 0;
2063 MP_DIGITS(&g) = 0;
2064 MP_CHECKOK( mp_init(&d) );
2065 MP_CHECKOK( mp_init_copy(&f, a) ); /* f = a */
2066 MP_CHECKOK( mp_init_copy(&g, p) ); /* g = p */
2068 mp_set(c, 1);
2069 mp_zero(&d);
2071 if (mp_cmp_z(&f) == 0) {
2072 res = MP_UNDEF;
2073 } else
2074 for (;;) {
2075 int diff_sign;
2076 while (mp_iseven(&f)) {
2077 mp_size n = mp_trailing_zeros(&f);
2078 if (!n) {
2079 res = MP_UNDEF;
2080 goto CLEANUP;
2081 }
2082 s_mp_div_2d(&f, n);
2083 MP_CHECKOK( s_mp_mul_2d(&d, n) );
2084 k += n;
2085 }
2086 if (mp_cmp_d(&f, 1) == MP_EQ) { /* f == 1 */
2087 res = k;
2088 break;
2089 }
2090 diff_sign = mp_cmp(&f, &g);
2091 if (diff_sign < 0) { /* f < g */
2092 s_mp_exch(&f, &g);
2093 s_mp_exch(c, &d);
2094 } else if (diff_sign == 0) { /* f == g */
2095 res = MP_UNDEF; /* a and p are not relatively prime */
2096 break;
2097 }
2098 if ((MP_DIGIT(&f,0) % 4) == (MP_DIGIT(&g,0) % 4)) {
2099 MP_CHECKOK( mp_sub(&f, &g, &f) ); /* f = f - g */
2100 MP_CHECKOK( mp_sub(c, &d, c) ); /* c = c - d */
2101 } else {
2102 MP_CHECKOK( mp_add(&f, &g, &f) ); /* f = f + g */
2103 MP_CHECKOK( mp_add(c, &d, c) ); /* c = c + d */
2104 }
2105 }
2106 if (res >= 0) {
2107 while (MP_SIGN(c) != MP_ZPOS) {
2108 MP_CHECKOK( mp_add(c, p, c) );
2109 }
2110 res = k;
2111 }
2113 CLEANUP:
2114 mp_clear(&d);
2115 mp_clear(&f);
2116 mp_clear(&g);
2117 return res;
2118 }
2120 /* Compute T = (P ** -1) mod MP_RADIX. Also works for 16-bit mp_digits.
2121 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)
2122 ** by Richard Schroeppel (a.k.a. Captain Nemo).
2123 */
2124 mp_digit s_mp_invmod_radix(mp_digit P)
2125 {
2126 mp_digit T = P;
2127 T *= 2 - (P * T);
2128 T *= 2 - (P * T);
2129 T *= 2 - (P * T);
2130 T *= 2 - (P * T);
2131 #if !defined(MP_USE_UINT_DIGIT)
2132 T *= 2 - (P * T);
2133 T *= 2 - (P * T);
2134 #endif
2135 return T;
2136 }
2138 /* Given c, k, and prime p, where a*c == 2**k (mod p),
2139 ** Compute x = (a ** -1) mod p. This is similar to Montgomery reduction.
2140 ** This technique from the paper "Fast Modular Reciprocals" (unpublished)
2141 ** by Richard Schroeppel (a.k.a. Captain Nemo).
2142 */
2143 mp_err s_mp_fixup_reciprocal(const mp_int *c, const mp_int *p, int k, mp_int *x)
2144 {
2145 int k_orig = k;
2146 mp_digit r;
2147 mp_size ix;
2148 mp_err res;
2150 if (mp_cmp_z(c) < 0) { /* c < 0 */
2151 MP_CHECKOK( mp_add(c, p, x) ); /* x = c + p */
2152 } else {
2153 MP_CHECKOK( mp_copy(c, x) ); /* x = c */
2154 }
2156 /* make sure x is large enough */
2157 ix = MP_HOWMANY(k, MP_DIGIT_BIT) + MP_USED(p) + 1;
2158 ix = MP_MAX(ix, MP_USED(x));
2159 MP_CHECKOK( s_mp_pad(x, ix) );
2161 r = 0 - s_mp_invmod_radix(MP_DIGIT(p,0));
2163 for (ix = 0; k > 0; ix++) {
2164 int j = MP_MIN(k, MP_DIGIT_BIT);
2165 mp_digit v = r * MP_DIGIT(x, ix);
2166 if (j < MP_DIGIT_BIT) {
2167 v &= ((mp_digit)1 << j) - 1; /* v = v mod (2 ** j) */
2168 }
2169 s_mp_mul_d_add_offset(p, v, x, ix); /* x += p * v * (RADIX ** ix) */
2170 k -= j;
2171 }
2172 s_mp_clamp(x);
2173 s_mp_div_2d(x, k_orig);
2174 res = MP_OKAY;
2176 CLEANUP:
2177 return res;
2178 }
2180 /* compute mod inverse using Schroeppel's method, only if m is odd */
2181 mp_err s_mp_invmod_odd_m(const mp_int *a, const mp_int *m, mp_int *c)
2182 {
2183 int k;
2184 mp_err res;
2185 mp_int x;
2187 ARGCHK(a && m && c, MP_BADARG);
2189 if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
2190 return MP_RANGE;
2191 if (mp_iseven(m))
2192 return MP_UNDEF;
2194 MP_DIGITS(&x) = 0;
2196 if (a == c) {
2197 if ((res = mp_init_copy(&x, a)) != MP_OKAY)
2198 return res;
2199 if (a == m)
2200 m = &x;
2201 a = &x;
2202 } else if (m == c) {
2203 if ((res = mp_init_copy(&x, m)) != MP_OKAY)
2204 return res;
2205 m = &x;
2206 } else {
2207 MP_DIGITS(&x) = 0;
2208 }
2210 MP_CHECKOK( s_mp_almost_inverse(a, m, c) );
2211 k = res;
2212 MP_CHECKOK( s_mp_fixup_reciprocal(c, m, k, c) );
2213 CLEANUP:
2214 mp_clear(&x);
2215 return res;
2216 }
2218 /* Known good algorithm for computing modular inverse. But slow. */
2219 mp_err mp_invmod_xgcd(const mp_int *a, const mp_int *m, mp_int *c)
2220 {
2221 mp_int g, x;
2222 mp_err res;
2224 ARGCHK(a && m && c, MP_BADARG);
2226 if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
2227 return MP_RANGE;
2229 MP_DIGITS(&g) = 0;
2230 MP_DIGITS(&x) = 0;
2231 MP_CHECKOK( mp_init(&x) );
2232 MP_CHECKOK( mp_init(&g) );
2234 MP_CHECKOK( mp_xgcd(a, m, &g, &x, NULL) );
2236 if (mp_cmp_d(&g, 1) != MP_EQ) {
2237 res = MP_UNDEF;
2238 goto CLEANUP;
2239 }
2241 res = mp_mod(&x, m, c);
2242 SIGN(c) = SIGN(a);
2244 CLEANUP:
2245 mp_clear(&x);
2246 mp_clear(&g);
2248 return res;
2249 }
2251 /* modular inverse where modulus is 2**k. */
2252 /* c = a**-1 mod 2**k */
2253 mp_err s_mp_invmod_2d(const mp_int *a, mp_size k, mp_int *c)
2254 {
2255 mp_err res;
2256 mp_size ix = k + 4;
2257 mp_int t0, t1, val, tmp, two2k;
2259 static const mp_digit d2 = 2;
2260 static const mp_int two = { MP_ZPOS, 1, 1, (mp_digit *)&d2 };
2262 if (mp_iseven(a))
2263 return MP_UNDEF;
2264 if (k <= MP_DIGIT_BIT) {
2265 mp_digit i = s_mp_invmod_radix(MP_DIGIT(a,0));
2266 if (k < MP_DIGIT_BIT)
2267 i &= ((mp_digit)1 << k) - (mp_digit)1;
2268 mp_set(c, i);
2269 return MP_OKAY;
2270 }
2271 MP_DIGITS(&t0) = 0;
2272 MP_DIGITS(&t1) = 0;
2273 MP_DIGITS(&val) = 0;
2274 MP_DIGITS(&tmp) = 0;
2275 MP_DIGITS(&two2k) = 0;
2276 MP_CHECKOK( mp_init_copy(&val, a) );
2277 s_mp_mod_2d(&val, k);
2278 MP_CHECKOK( mp_init_copy(&t0, &val) );
2279 MP_CHECKOK( mp_init_copy(&t1, &t0) );
2280 MP_CHECKOK( mp_init(&tmp) );
2281 MP_CHECKOK( mp_init(&two2k) );
2282 MP_CHECKOK( s_mp_2expt(&two2k, k) );
2283 do {
2284 MP_CHECKOK( mp_mul(&val, &t1, &tmp) );
2285 MP_CHECKOK( mp_sub(&two, &tmp, &tmp) );
2286 MP_CHECKOK( mp_mul(&t1, &tmp, &t1) );
2287 s_mp_mod_2d(&t1, k);
2288 while (MP_SIGN(&t1) != MP_ZPOS) {
2289 MP_CHECKOK( mp_add(&t1, &two2k, &t1) );
2290 }
2291 if (mp_cmp(&t1, &t0) == MP_EQ)
2292 break;
2293 MP_CHECKOK( mp_copy(&t1, &t0) );
2294 } while (--ix > 0);
2295 if (!ix) {
2296 res = MP_UNDEF;
2297 } else {
2298 mp_exch(c, &t1);
2299 }
2301 CLEANUP:
2302 mp_clear(&t0);
2303 mp_clear(&t1);
2304 mp_clear(&val);
2305 mp_clear(&tmp);
2306 mp_clear(&two2k);
2307 return res;
2308 }
2310 mp_err s_mp_invmod_even_m(const mp_int *a, const mp_int *m, mp_int *c)
2311 {
2312 mp_err res;
2313 mp_size k;
2314 mp_int oddFactor, evenFactor; /* factors of the modulus */
2315 mp_int oddPart, evenPart; /* parts to combine via CRT. */
2316 mp_int C2, tmp1, tmp2;
2318 /*static const mp_digit d1 = 1; */
2319 /*static const mp_int one = { MP_ZPOS, 1, 1, (mp_digit *)&d1 }; */
2321 if ((res = s_mp_ispow2(m)) >= 0) {
2322 k = res;
2323 return s_mp_invmod_2d(a, k, c);
2324 }
2325 MP_DIGITS(&oddFactor) = 0;
2326 MP_DIGITS(&evenFactor) = 0;
2327 MP_DIGITS(&oddPart) = 0;
2328 MP_DIGITS(&evenPart) = 0;
2329 MP_DIGITS(&C2) = 0;
2330 MP_DIGITS(&tmp1) = 0;
2331 MP_DIGITS(&tmp2) = 0;
2333 MP_CHECKOK( mp_init_copy(&oddFactor, m) ); /* oddFactor = m */
2334 MP_CHECKOK( mp_init(&evenFactor) );
2335 MP_CHECKOK( mp_init(&oddPart) );
2336 MP_CHECKOK( mp_init(&evenPart) );
2337 MP_CHECKOK( mp_init(&C2) );
2338 MP_CHECKOK( mp_init(&tmp1) );
2339 MP_CHECKOK( mp_init(&tmp2) );
2341 k = mp_trailing_zeros(m);
2342 s_mp_div_2d(&oddFactor, k);
2343 MP_CHECKOK( s_mp_2expt(&evenFactor, k) );
2345 /* compute a**-1 mod oddFactor. */
2346 MP_CHECKOK( s_mp_invmod_odd_m(a, &oddFactor, &oddPart) );
2347 /* compute a**-1 mod evenFactor, where evenFactor == 2**k. */
2348 MP_CHECKOK( s_mp_invmod_2d( a, k, &evenPart) );
2350 /* Use Chinese Remainer theorem to compute a**-1 mod m. */
2351 /* let m1 = oddFactor, v1 = oddPart,
2352 * let m2 = evenFactor, v2 = evenPart.
2353 */
2355 /* Compute C2 = m1**-1 mod m2. */
2356 MP_CHECKOK( s_mp_invmod_2d(&oddFactor, k, &C2) );
2358 /* compute u = (v2 - v1)*C2 mod m2 */
2359 MP_CHECKOK( mp_sub(&evenPart, &oddPart, &tmp1) );
2360 MP_CHECKOK( mp_mul(&tmp1, &C2, &tmp2) );
2361 s_mp_mod_2d(&tmp2, k);
2362 while (MP_SIGN(&tmp2) != MP_ZPOS) {
2363 MP_CHECKOK( mp_add(&tmp2, &evenFactor, &tmp2) );
2364 }
2366 /* compute answer = v1 + u*m1 */
2367 MP_CHECKOK( mp_mul(&tmp2, &oddFactor, c) );
2368 MP_CHECKOK( mp_add(&oddPart, c, c) );
2369 /* not sure this is necessary, but it's low cost if not. */
2370 MP_CHECKOK( mp_mod(c, m, c) );
2372 CLEANUP:
2373 mp_clear(&oddFactor);
2374 mp_clear(&evenFactor);
2375 mp_clear(&oddPart);
2376 mp_clear(&evenPart);
2377 mp_clear(&C2);
2378 mp_clear(&tmp1);
2379 mp_clear(&tmp2);
2380 return res;
2381 }
2384 /* {{{ mp_invmod(a, m, c) */
2386 /*
2387 mp_invmod(a, m, c)
2389 Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
2390 This is equivalent to the question of whether (a, m) = 1. If not,
2391 MP_UNDEF is returned, and there is no inverse.
2392 */
2394 mp_err mp_invmod(const mp_int *a, const mp_int *m, mp_int *c)
2395 {
2397 ARGCHK(a && m && c, MP_BADARG);
2399 if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
2400 return MP_RANGE;
2402 if (mp_isodd(m)) {
2403 return s_mp_invmod_odd_m(a, m, c);
2404 }
2405 if (mp_iseven(a))
2406 return MP_UNDEF; /* not invertable */
2408 return s_mp_invmod_even_m(a, m, c);
2410 } /* end mp_invmod() */
2412 /* }}} */
2413 #endif /* if MP_NUMTH */
2415 /* }}} */
2417 /*------------------------------------------------------------------------*/
2418 /* {{{ mp_print(mp, ofp) */
2420 #if MP_IOFUNC
2421 /*
2422 mp_print(mp, ofp)
2424 Print a textual representation of the given mp_int on the output
2425 stream 'ofp'. Output is generated using the internal radix.
2426 */
2428 void mp_print(mp_int *mp, FILE *ofp)
2429 {
2430 int ix;
2432 if(mp == NULL || ofp == NULL)
2433 return;
2435 fputc((SIGN(mp) == NEG) ? '-' : '+', ofp);
2437 for(ix = USED(mp) - 1; ix >= 0; ix--) {
2438 fprintf(ofp, DIGIT_FMT, DIGIT(mp, ix));
2439 }
2441 } /* end mp_print() */
2443 #endif /* if MP_IOFUNC */
2445 /* }}} */
2447 /*------------------------------------------------------------------------*/
2448 /* {{{ More I/O Functions */
2450 /* {{{ mp_read_raw(mp, str, len) */
2452 /*
2453 mp_read_raw(mp, str, len)
2455 Read in a raw value (base 256) into the given mp_int
2456 */
2458 mp_err mp_read_raw(mp_int *mp, char *str, int len)
2459 {
2460 int ix;
2461 mp_err res;
2462 unsigned char *ustr = (unsigned char *)str;
2464 ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
2466 mp_zero(mp);
2468 /* Get sign from first byte */
2469 if(ustr[0])
2470 SIGN(mp) = NEG;
2471 else
2472 SIGN(mp) = ZPOS;
2474 /* Read the rest of the digits */
2475 for(ix = 1; ix < len; ix++) {
2476 if((res = mp_mul_d(mp, 256, mp)) != MP_OKAY)
2477 return res;
2478 if((res = mp_add_d(mp, ustr[ix], mp)) != MP_OKAY)
2479 return res;
2480 }
2482 return MP_OKAY;
2484 } /* end mp_read_raw() */
2486 /* }}} */
2488 /* {{{ mp_raw_size(mp) */
2490 int mp_raw_size(mp_int *mp)
2491 {
2492 ARGCHK(mp != NULL, 0);
2494 return (USED(mp) * sizeof(mp_digit)) + 1;
2496 } /* end mp_raw_size() */
2498 /* }}} */
2500 /* {{{ mp_toraw(mp, str) */
2502 mp_err mp_toraw(mp_int *mp, char *str)
2503 {
2504 int ix, jx, pos = 1;
2506 ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2508 str[0] = (char)SIGN(mp);
2510 /* Iterate over each digit... */
2511 for(ix = USED(mp) - 1; ix >= 0; ix--) {
2512 mp_digit d = DIGIT(mp, ix);
2514 /* Unpack digit bytes, high order first */
2515 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {
2516 str[pos++] = (char)(d >> (jx * CHAR_BIT));
2517 }
2518 }
2520 return MP_OKAY;
2522 } /* end mp_toraw() */
2524 /* }}} */
2526 /* {{{ mp_read_radix(mp, str, radix) */
2528 /*
2529 mp_read_radix(mp, str, radix)
2531 Read an integer from the given string, and set mp to the resulting
2532 value. The input is presumed to be in base 10. Leading non-digit
2533 characters are ignored, and the function reads until a non-digit
2534 character or the end of the string.
2535 */
2537 mp_err mp_read_radix(mp_int *mp, const char *str, int radix)
2538 {
2539 int ix = 0, val = 0;
2540 mp_err res;
2541 mp_sign sig = ZPOS;
2543 ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX,
2544 MP_BADARG);
2546 mp_zero(mp);
2548 /* Skip leading non-digit characters until a digit or '-' or '+' */
2549 while(str[ix] &&
2550 (s_mp_tovalue(str[ix], radix) < 0) &&
2551 str[ix] != '-' &&
2552 str[ix] != '+') {
2553 ++ix;
2554 }
2556 if(str[ix] == '-') {
2557 sig = NEG;
2558 ++ix;
2559 } else if(str[ix] == '+') {
2560 sig = ZPOS; /* this is the default anyway... */
2561 ++ix;
2562 }
2564 while((val = s_mp_tovalue(str[ix], radix)) >= 0) {
2565 if((res = s_mp_mul_d(mp, radix)) != MP_OKAY)
2566 return res;
2567 if((res = s_mp_add_d(mp, val)) != MP_OKAY)
2568 return res;
2569 ++ix;
2570 }
2572 if(s_mp_cmp_d(mp, 0) == MP_EQ)
2573 SIGN(mp) = ZPOS;
2574 else
2575 SIGN(mp) = sig;
2577 return MP_OKAY;
2579 } /* end mp_read_radix() */
2581 mp_err mp_read_variable_radix(mp_int *a, const char * str, int default_radix)
2582 {
2583 int radix = default_radix;
2584 int cx;
2585 mp_sign sig = ZPOS;
2586 mp_err res;
2588 /* Skip leading non-digit characters until a digit or '-' or '+' */
2589 while ((cx = *str) != 0 &&
2590 (s_mp_tovalue(cx, radix) < 0) &&
2591 cx != '-' &&
2592 cx != '+') {
2593 ++str;
2594 }
2596 if (cx == '-') {
2597 sig = NEG;
2598 ++str;
2599 } else if (cx == '+') {
2600 sig = ZPOS; /* this is the default anyway... */
2601 ++str;
2602 }
2604 if (str[0] == '0') {
2605 if ((str[1] | 0x20) == 'x') {
2606 radix = 16;
2607 str += 2;
2608 } else {
2609 radix = 8;
2610 str++;
2611 }
2612 }
2613 res = mp_read_radix(a, str, radix);
2614 if (res == MP_OKAY) {
2615 MP_SIGN(a) = (s_mp_cmp_d(a, 0) == MP_EQ) ? ZPOS : sig;
2616 }
2617 return res;
2618 }
2620 /* }}} */
2622 /* {{{ mp_radix_size(mp, radix) */
2624 int mp_radix_size(mp_int *mp, int radix)
2625 {
2626 int bits;
2628 if(!mp || radix < 2 || radix > MAX_RADIX)
2629 return 0;
2631 bits = USED(mp) * DIGIT_BIT - 1;
2633 return s_mp_outlen(bits, radix);
2635 } /* end mp_radix_size() */
2637 /* }}} */
2639 /* {{{ mp_toradix(mp, str, radix) */
2641 mp_err mp_toradix(mp_int *mp, char *str, int radix)
2642 {
2643 int ix, pos = 0;
2645 ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2646 ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);
2648 if(mp_cmp_z(mp) == MP_EQ) {
2649 str[0] = '0';
2650 str[1] = '\0';
2651 } else {
2652 mp_err res;
2653 mp_int tmp;
2654 mp_sign sgn;
2655 mp_digit rem, rdx = (mp_digit)radix;
2656 char ch;
2658 if((res = mp_init_copy(&tmp, mp)) != MP_OKAY)
2659 return res;
2661 /* Save sign for later, and take absolute value */
2662 sgn = SIGN(&tmp); SIGN(&tmp) = ZPOS;
2664 /* Generate output digits in reverse order */
2665 while(mp_cmp_z(&tmp) != 0) {
2666 if((res = mp_div_d(&tmp, rdx, &tmp, &rem)) != MP_OKAY) {
2667 mp_clear(&tmp);
2668 return res;
2669 }
2671 /* Generate digits, use capital letters */
2672 ch = s_mp_todigit(rem, radix, 0);
2674 str[pos++] = ch;
2675 }
2677 /* Add - sign if original value was negative */
2678 if(sgn == NEG)
2679 str[pos++] = '-';
2681 /* Add trailing NUL to end the string */
2682 str[pos--] = '\0';
2684 /* Reverse the digits and sign indicator */
2685 ix = 0;
2686 while(ix < pos) {
2687 char tmp = str[ix];
2689 str[ix] = str[pos];
2690 str[pos] = tmp;
2691 ++ix;
2692 --pos;
2693 }
2695 mp_clear(&tmp);
2696 }
2698 return MP_OKAY;
2700 } /* end mp_toradix() */
2702 /* }}} */
2704 /* {{{ mp_tovalue(ch, r) */
2706 int mp_tovalue(char ch, int r)
2707 {
2708 return s_mp_tovalue(ch, r);
2710 } /* end mp_tovalue() */
2712 /* }}} */
2714 /* }}} */
2716 /* {{{ mp_strerror(ec) */
2718 /*
2719 mp_strerror(ec)
2721 Return a string describing the meaning of error code 'ec'. The
2722 string returned is allocated in static memory, so the caller should
2723 not attempt to modify or free the memory associated with this
2724 string.
2725 */
2726 const char *mp_strerror(mp_err ec)
2727 {
2728 int aec = (ec < 0) ? -ec : ec;
2730 /* Code values are negative, so the senses of these comparisons
2731 are accurate */
2732 if(ec < MP_LAST_CODE || ec > MP_OKAY) {
2733 return mp_err_string[0]; /* unknown error code */
2734 } else {
2735 return mp_err_string[aec + 1];
2736 }
2738 } /* end mp_strerror() */
2740 /* }}} */
2742 /*========================================================================*/
2743 /*------------------------------------------------------------------------*/
2744 /* Static function definitions (internal use only) */
2746 /* {{{ Memory management */
2748 /* {{{ s_mp_grow(mp, min) */
2750 /* Make sure there are at least 'min' digits allocated to mp */
2751 mp_err s_mp_grow(mp_int *mp, mp_size min)
2752 {
2753 if(min > ALLOC(mp)) {
2754 mp_digit *tmp;
2756 /* Set min to next nearest default precision block size */
2757 min = MP_ROUNDUP(min, s_mp_defprec);
2759 if((tmp = s_mp_alloc(min, sizeof(mp_digit))) == NULL)
2760 return MP_MEM;
2762 s_mp_copy(DIGITS(mp), tmp, USED(mp));
2764 #if MP_CRYPTO
2765 s_mp_setz(DIGITS(mp), ALLOC(mp));
2766 #endif
2767 s_mp_free(DIGITS(mp));
2768 DIGITS(mp) = tmp;
2769 ALLOC(mp) = min;
2770 }
2772 return MP_OKAY;
2774 } /* end s_mp_grow() */
2776 /* }}} */
2778 /* {{{ s_mp_pad(mp, min) */
2780 /* Make sure the used size of mp is at least 'min', growing if needed */
2781 mp_err s_mp_pad(mp_int *mp, mp_size min)
2782 {
2783 if(min > USED(mp)) {
2784 mp_err res;
2786 /* Make sure there is room to increase precision */
2787 if (min > ALLOC(mp)) {
2788 if ((res = s_mp_grow(mp, min)) != MP_OKAY)
2789 return res;
2790 } else {
2791 s_mp_setz(DIGITS(mp) + USED(mp), min - USED(mp));
2792 }
2794 /* Increase precision; should already be 0-filled */
2795 USED(mp) = min;
2796 }
2798 return MP_OKAY;
2800 } /* end s_mp_pad() */
2802 /* }}} */
2804 /* {{{ s_mp_setz(dp, count) */
2806 #if MP_MACRO == 0
2807 /* Set 'count' digits pointed to by dp to be zeroes */
2808 void s_mp_setz(mp_digit *dp, mp_size count)
2809 {
2810 #if MP_MEMSET == 0
2811 int ix;
2813 for(ix = 0; ix < count; ix++)
2814 dp[ix] = 0;
2815 #else
2816 memset(dp, 0, count * sizeof(mp_digit));
2817 #endif
2819 } /* end s_mp_setz() */
2820 #endif
2822 /* }}} */
2824 /* {{{ s_mp_copy(sp, dp, count) */
2826 #if MP_MACRO == 0
2827 /* Copy 'count' digits from sp to dp */
2828 void s_mp_copy(const mp_digit *sp, mp_digit *dp, mp_size count)
2829 {
2830 #if MP_MEMCPY == 0
2831 int ix;
2833 for(ix = 0; ix < count; ix++)
2834 dp[ix] = sp[ix];
2835 #else
2836 memcpy(dp, sp, count * sizeof(mp_digit));
2837 #endif
2838 ++mp_copies;
2840 } /* end s_mp_copy() */
2841 #endif
2843 /* }}} */
2845 /* {{{ s_mp_alloc(nb, ni) */
2847 #if MP_MACRO == 0
2848 /* Allocate ni records of nb bytes each, and return a pointer to that */
2849 void *s_mp_alloc(size_t nb, size_t ni)
2850 {
2851 ++mp_allocs;
2852 return calloc(nb, ni);
2854 } /* end s_mp_alloc() */
2855 #endif
2857 /* }}} */
2859 /* {{{ s_mp_free(ptr) */
2861 #if MP_MACRO == 0
2862 /* Free the memory pointed to by ptr */
2863 void s_mp_free(void *ptr)
2864 {
2865 if(ptr) {
2866 ++mp_frees;
2867 free(ptr);
2868 }
2869 } /* end s_mp_free() */
2870 #endif
2872 /* }}} */
2874 /* {{{ s_mp_clamp(mp) */
2876 #if MP_MACRO == 0
2877 /* Remove leading zeroes from the given value */
2878 void s_mp_clamp(mp_int *mp)
2879 {
2880 mp_size used = MP_USED(mp);
2881 while (used > 1 && DIGIT(mp, used - 1) == 0)
2882 --used;
2883 MP_USED(mp) = used;
2884 } /* end s_mp_clamp() */
2885 #endif
2887 /* }}} */
2889 /* {{{ s_mp_exch(a, b) */
2891 /* Exchange the data for a and b; (b, a) = (a, b) */
2892 void s_mp_exch(mp_int *a, mp_int *b)
2893 {
2894 mp_int tmp;
2896 tmp = *a;
2897 *a = *b;
2898 *b = tmp;
2900 } /* end s_mp_exch() */
2902 /* }}} */
2904 /* }}} */
2906 /* {{{ Arithmetic helpers */
2908 /* {{{ s_mp_lshd(mp, p) */
2910 /*
2911 Shift mp leftward by p digits, growing if needed, and zero-filling
2912 the in-shifted digits at the right end. This is a convenient
2913 alternative to multiplication by powers of the radix
2914 */
2916 mp_err s_mp_lshd(mp_int *mp, mp_size p)
2917 {
2918 mp_err res;
2919 mp_size pos;
2920 int ix;
2922 if(p == 0)
2923 return MP_OKAY;
2925 if (MP_USED(mp) == 1 && MP_DIGIT(mp, 0) == 0)
2926 return MP_OKAY;
2928 if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY)
2929 return res;
2931 pos = USED(mp) - 1;
2933 /* Shift all the significant figures over as needed */
2934 for(ix = pos - p; ix >= 0; ix--)
2935 DIGIT(mp, ix + p) = DIGIT(mp, ix);
2937 /* Fill the bottom digits with zeroes */
2938 for(ix = 0; ix < p; ix++)
2939 DIGIT(mp, ix) = 0;
2941 return MP_OKAY;
2943 } /* end s_mp_lshd() */
2945 /* }}} */
2947 /* {{{ s_mp_mul_2d(mp, d) */
2949 /*
2950 Multiply the integer by 2^d, where d is a number of bits. This
2951 amounts to a bitwise shift of the value.
2952 */
2953 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d)
2954 {
2955 mp_err res;
2956 mp_digit dshift, bshift;
2957 mp_digit mask;
2959 ARGCHK(mp != NULL, MP_BADARG);
2961 dshift = d / MP_DIGIT_BIT;
2962 bshift = d % MP_DIGIT_BIT;
2963 /* bits to be shifted out of the top word */
2964 mask = ((mp_digit)~0 << (MP_DIGIT_BIT - bshift));
2965 mask &= MP_DIGIT(mp, MP_USED(mp) - 1);
2967 if (MP_OKAY != (res = s_mp_pad(mp, MP_USED(mp) + dshift + (mask != 0) )))
2968 return res;
2970 if (dshift && MP_OKAY != (res = s_mp_lshd(mp, dshift)))
2971 return res;
2973 if (bshift) {
2974 mp_digit *pa = MP_DIGITS(mp);
2975 mp_digit *alim = pa + MP_USED(mp);
2976 mp_digit prev = 0;
2978 for (pa += dshift; pa < alim; ) {
2979 mp_digit x = *pa;
2980 *pa++ = (x << bshift) | prev;
2981 prev = x >> (DIGIT_BIT - bshift);
2982 }
2983 }
2985 s_mp_clamp(mp);
2986 return MP_OKAY;
2987 } /* end s_mp_mul_2d() */
2989 /* {{{ s_mp_rshd(mp, p) */
2991 /*
2992 Shift mp rightward by p digits. Maintains the invariant that
2993 digits above the precision are all zero. Digits shifted off the
2994 end are lost. Cannot fail.
2995 */
2997 void s_mp_rshd(mp_int *mp, mp_size p)
2998 {
2999 mp_size ix;
3000 mp_digit *src, *dst;
3002 if(p == 0)
3003 return;
3005 /* Shortcut when all digits are to be shifted off */
3006 if(p >= USED(mp)) {
3007 s_mp_setz(DIGITS(mp), ALLOC(mp));
3008 USED(mp) = 1;
3009 SIGN(mp) = ZPOS;
3010 return;
3011 }
3013 /* Shift all the significant figures over as needed */
3014 dst = MP_DIGITS(mp);
3015 src = dst + p;
3016 for (ix = USED(mp) - p; ix > 0; ix--)
3017 *dst++ = *src++;
3019 MP_USED(mp) -= p;
3020 /* Fill the top digits with zeroes */
3021 while (p-- > 0)
3022 *dst++ = 0;
3024 #if 0
3025 /* Strip off any leading zeroes */
3026 s_mp_clamp(mp);
3027 #endif
3029 } /* end s_mp_rshd() */
3031 /* }}} */
3033 /* {{{ s_mp_div_2(mp) */
3035 /* Divide by two -- take advantage of radix properties to do it fast */
3036 void s_mp_div_2(mp_int *mp)
3037 {
3038 s_mp_div_2d(mp, 1);
3040 } /* end s_mp_div_2() */
3042 /* }}} */
3044 /* {{{ s_mp_mul_2(mp) */
3046 mp_err s_mp_mul_2(mp_int *mp)
3047 {
3048 mp_digit *pd;
3049 int ix, used;
3050 mp_digit kin = 0;
3052 /* Shift digits leftward by 1 bit */
3053 used = MP_USED(mp);
3054 pd = MP_DIGITS(mp);
3055 for (ix = 0; ix < used; ix++) {
3056 mp_digit d = *pd;
3057 *pd++ = (d << 1) | kin;
3058 kin = (d >> (DIGIT_BIT - 1));
3059 }
3061 /* Deal with rollover from last digit */
3062 if (kin) {
3063 if (ix >= ALLOC(mp)) {
3064 mp_err res;
3065 if((res = s_mp_grow(mp, ALLOC(mp) + 1)) != MP_OKAY)
3066 return res;
3067 }
3069 DIGIT(mp, ix) = kin;
3070 USED(mp) += 1;
3071 }
3073 return MP_OKAY;
3075 } /* end s_mp_mul_2() */
3077 /* }}} */
3079 /* {{{ s_mp_mod_2d(mp, d) */
3081 /*
3082 Remainder the integer by 2^d, where d is a number of bits. This
3083 amounts to a bitwise AND of the value, and does not require the full
3084 division code
3085 */
3086 void s_mp_mod_2d(mp_int *mp, mp_digit d)
3087 {
3088 mp_size ndig = (d / DIGIT_BIT), nbit = (d % DIGIT_BIT);
3089 mp_size ix;
3090 mp_digit dmask;
3092 if(ndig >= USED(mp))
3093 return;
3095 /* Flush all the bits above 2^d in its digit */
3096 dmask = ((mp_digit)1 << nbit) - 1;
3097 DIGIT(mp, ndig) &= dmask;
3099 /* Flush all digits above the one with 2^d in it */
3100 for(ix = ndig + 1; ix < USED(mp); ix++)
3101 DIGIT(mp, ix) = 0;
3103 s_mp_clamp(mp);
3105 } /* end s_mp_mod_2d() */
3107 /* }}} */
3109 /* {{{ s_mp_div_2d(mp, d) */
3111 /*
3112 Divide the integer by 2^d, where d is a number of bits. This
3113 amounts to a bitwise shift of the value, and does not require the
3114 full division code (used in Barrett reduction, see below)
3115 */
3116 void s_mp_div_2d(mp_int *mp, mp_digit d)
3117 {
3118 int ix;
3119 mp_digit save, next, mask;
3121 s_mp_rshd(mp, d / DIGIT_BIT);
3122 d %= DIGIT_BIT;
3123 if (d) {
3124 mask = ((mp_digit)1 << d) - 1;
3125 save = 0;
3126 for(ix = USED(mp) - 1; ix >= 0; ix--) {
3127 next = DIGIT(mp, ix) & mask;
3128 DIGIT(mp, ix) = (DIGIT(mp, ix) >> d) | (save << (DIGIT_BIT - d));
3129 save = next;
3130 }
3131 }
3132 s_mp_clamp(mp);
3134 } /* end s_mp_div_2d() */
3136 /* }}} */
3138 /* {{{ s_mp_norm(a, b, *d) */
3140 /*
3141 s_mp_norm(a, b, *d)
3143 Normalize a and b for division, where b is the divisor. In order
3144 that we might make good guesses for quotient digits, we want the
3145 leading digit of b to be at least half the radix, which we
3146 accomplish by multiplying a and b by a power of 2. The exponent
3147 (shift count) is placed in *pd, so that the remainder can be shifted
3148 back at the end of the division process.
3149 */
3151 mp_err s_mp_norm(mp_int *a, mp_int *b, mp_digit *pd)
3152 {
3153 mp_digit d;
3154 mp_digit mask;
3155 mp_digit b_msd;
3156 mp_err res = MP_OKAY;
3158 d = 0;
3159 mask = DIGIT_MAX & ~(DIGIT_MAX >> 1); /* mask is msb of digit */
3160 b_msd = DIGIT(b, USED(b) - 1);
3161 while (!(b_msd & mask)) {
3162 b_msd <<= 1;
3163 ++d;
3164 }
3166 if (d) {
3167 MP_CHECKOK( s_mp_mul_2d(a, d) );
3168 MP_CHECKOK( s_mp_mul_2d(b, d) );
3169 }
3171 *pd = d;
3172 CLEANUP:
3173 return res;
3175 } /* end s_mp_norm() */
3177 /* }}} */
3179 /* }}} */
3181 /* {{{ Primitive digit arithmetic */
3183 /* {{{ s_mp_add_d(mp, d) */
3185 /* Add d to |mp| in place */
3186 mp_err s_mp_add_d(mp_int *mp, mp_digit d) /* unsigned digit addition */
3187 {
3188 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3189 mp_word w, k = 0;
3190 mp_size ix = 1;
3192 w = (mp_word)DIGIT(mp, 0) + d;
3193 DIGIT(mp, 0) = ACCUM(w);
3194 k = CARRYOUT(w);
3196 while(ix < USED(mp) && k) {
3197 w = (mp_word)DIGIT(mp, ix) + k;
3198 DIGIT(mp, ix) = ACCUM(w);
3199 k = CARRYOUT(w);
3200 ++ix;
3201 }
3203 if(k != 0) {
3204 mp_err res;
3206 if((res = s_mp_pad(mp, USED(mp) + 1)) != MP_OKAY)
3207 return res;
3209 DIGIT(mp, ix) = (mp_digit)k;
3210 }
3212 return MP_OKAY;
3213 #else
3214 mp_digit * pmp = MP_DIGITS(mp);
3215 mp_digit sum, mp_i, carry = 0;
3216 mp_err res = MP_OKAY;
3217 int used = (int)MP_USED(mp);
3219 mp_i = *pmp;
3220 *pmp++ = sum = d + mp_i;
3221 carry = (sum < d);
3222 while (carry && --used > 0) {
3223 mp_i = *pmp;
3224 *pmp++ = sum = carry + mp_i;
3225 carry = !sum;
3226 }
3227 if (carry && !used) {
3228 /* mp is growing */
3229 used = MP_USED(mp);
3230 MP_CHECKOK( s_mp_pad(mp, used + 1) );
3231 MP_DIGIT(mp, used) = carry;
3232 }
3233 CLEANUP:
3234 return res;
3235 #endif
3236 } /* end s_mp_add_d() */
3238 /* }}} */
3240 /* {{{ s_mp_sub_d(mp, d) */
3242 /* Subtract d from |mp| in place, assumes |mp| > d */
3243 mp_err s_mp_sub_d(mp_int *mp, mp_digit d) /* unsigned digit subtract */
3244 {
3245 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3246 mp_word w, b = 0;
3247 mp_size ix = 1;
3249 /* Compute initial subtraction */
3250 w = (RADIX + (mp_word)DIGIT(mp, 0)) - d;
3251 b = CARRYOUT(w) ? 0 : 1;
3252 DIGIT(mp, 0) = ACCUM(w);
3254 /* Propagate borrows leftward */
3255 while(b && ix < USED(mp)) {
3256 w = (RADIX + (mp_word)DIGIT(mp, ix)) - b;
3257 b = CARRYOUT(w) ? 0 : 1;
3258 DIGIT(mp, ix) = ACCUM(w);
3259 ++ix;
3260 }
3262 /* Remove leading zeroes */
3263 s_mp_clamp(mp);
3265 /* If we have a borrow out, it's a violation of the input invariant */
3266 if(b)
3267 return MP_RANGE;
3268 else
3269 return MP_OKAY;
3270 #else
3271 mp_digit *pmp = MP_DIGITS(mp);
3272 mp_digit mp_i, diff, borrow;
3273 mp_size used = MP_USED(mp);
3275 mp_i = *pmp;
3276 *pmp++ = diff = mp_i - d;
3277 borrow = (diff > mp_i);
3278 while (borrow && --used) {
3279 mp_i = *pmp;
3280 *pmp++ = diff = mp_i - borrow;
3281 borrow = (diff > mp_i);
3282 }
3283 s_mp_clamp(mp);
3284 return (borrow && !used) ? MP_RANGE : MP_OKAY;
3285 #endif
3286 } /* end s_mp_sub_d() */
3288 /* }}} */
3290 /* {{{ s_mp_mul_d(a, d) */
3292 /* Compute a = a * d, single digit multiplication */
3293 mp_err s_mp_mul_d(mp_int *a, mp_digit d)
3294 {
3295 mp_err res;
3296 mp_size used;
3297 int pow;
3299 if (!d) {
3300 mp_zero(a);
3301 return MP_OKAY;
3302 }
3303 if (d == 1)
3304 return MP_OKAY;
3305 if (0 <= (pow = s_mp_ispow2d(d))) {
3306 return s_mp_mul_2d(a, (mp_digit)pow);
3307 }
3309 used = MP_USED(a);
3310 MP_CHECKOK( s_mp_pad(a, used + 1) );
3312 s_mpv_mul_d(MP_DIGITS(a), used, d, MP_DIGITS(a));
3314 s_mp_clamp(a);
3316 CLEANUP:
3317 return res;
3319 } /* end s_mp_mul_d() */
3321 /* }}} */
3323 /* {{{ s_mp_div_d(mp, d, r) */
3325 /*
3326 s_mp_div_d(mp, d, r)
3328 Compute the quotient mp = mp / d and remainder r = mp mod d, for a
3329 single digit d. If r is null, the remainder will be discarded.
3330 */
3332 mp_err s_mp_div_d(mp_int *mp, mp_digit d, mp_digit *r)
3333 {
3334 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
3335 mp_word w = 0, q;
3336 #else
3337 mp_digit w, q;
3338 #endif
3339 int ix;
3340 mp_err res;
3341 mp_int quot;
3342 mp_int rem;
3344 if(d == 0)
3345 return MP_RANGE;
3346 if (d == 1) {
3347 if (r)
3348 *r = 0;
3349 return MP_OKAY;
3350 }
3351 /* could check for power of 2 here, but mp_div_d does that. */
3352 if (MP_USED(mp) == 1) {
3353 mp_digit n = MP_DIGIT(mp,0);
3354 mp_digit rem;
3356 q = n / d;
3357 rem = n % d;
3358 MP_DIGIT(mp,0) = q;
3359 if (r)
3360 *r = rem;
3361 return MP_OKAY;
3362 }
3364 MP_DIGITS(&rem) = 0;
3365 MP_DIGITS(") = 0;
3366 /* Make room for the quotient */
3367 MP_CHECKOK( mp_init_size(", USED(mp)) );
3369 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
3370 for(ix = USED(mp) - 1; ix >= 0; ix--) {
3371 w = (w << DIGIT_BIT) | DIGIT(mp, ix);
3373 if(w >= d) {
3374 q = w / d;
3375 w = w % d;
3376 } else {
3377 q = 0;
3378 }
3380 s_mp_lshd(", 1);
3381 DIGIT(", 0) = (mp_digit)q;
3382 }
3383 #else
3384 {
3385 mp_digit p;
3386 #if !defined(MP_ASSEMBLY_DIV_2DX1D)
3387 mp_digit norm;
3388 #endif
3390 MP_CHECKOK( mp_init_copy(&rem, mp) );
3392 #if !defined(MP_ASSEMBLY_DIV_2DX1D)
3393 MP_DIGIT(", 0) = d;
3394 MP_CHECKOK( s_mp_norm(&rem, ", &norm) );
3395 if (norm)
3396 d <<= norm;
3397 MP_DIGIT(", 0) = 0;
3398 #endif
3400 p = 0;
3401 for (ix = USED(&rem) - 1; ix >= 0; ix--) {
3402 w = DIGIT(&rem, ix);
3404 if (p) {
3405 MP_CHECKOK( s_mpv_div_2dx1d(p, w, d, &q, &w) );
3406 } else if (w >= d) {
3407 q = w / d;
3408 w = w % d;
3409 } else {
3410 q = 0;
3411 }
3413 MP_CHECKOK( s_mp_lshd(", 1) );
3414 DIGIT(", 0) = q;
3415 p = w;
3416 }
3417 #if !defined(MP_ASSEMBLY_DIV_2DX1D)
3418 if (norm)
3419 w >>= norm;
3420 #endif
3421 }
3422 #endif
3424 /* Deliver the remainder, if desired */
3425 if(r)
3426 *r = (mp_digit)w;
3428 s_mp_clamp(");
3429 mp_exch(", mp);
3430 CLEANUP:
3431 mp_clear(");
3432 mp_clear(&rem);
3434 return res;
3435 } /* end s_mp_div_d() */
3437 /* }}} */
3440 /* }}} */
3442 /* {{{ Primitive full arithmetic */
3444 /* {{{ s_mp_add(a, b) */
3446 /* Compute a = |a| + |b| */
3447 mp_err s_mp_add(mp_int *a, const mp_int *b) /* magnitude addition */
3448 {
3449 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3450 mp_word w = 0;
3451 #else
3452 mp_digit d, sum, carry = 0;
3453 #endif
3454 mp_digit *pa, *pb;
3455 mp_size ix;
3456 mp_size used;
3457 mp_err res;
3459 /* Make sure a has enough precision for the output value */
3460 if((USED(b) > USED(a)) && (res = s_mp_pad(a, USED(b))) != MP_OKAY)
3461 return res;
3463 /*
3464 Add up all digits up to the precision of b. If b had initially
3465 the same precision as a, or greater, we took care of it by the
3466 padding step above, so there is no problem. If b had initially
3467 less precision, we'll have to make sure the carry out is duly
3468 propagated upward among the higher-order digits of the sum.
3469 */
3470 pa = MP_DIGITS(a);
3471 pb = MP_DIGITS(b);
3472 used = MP_USED(b);
3473 for(ix = 0; ix < used; ix++) {
3474 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3475 w = w + *pa + *pb++;
3476 *pa++ = ACCUM(w);
3477 w = CARRYOUT(w);
3478 #else
3479 d = *pa;
3480 sum = d + *pb++;
3481 d = (sum < d); /* detect overflow */
3482 *pa++ = sum += carry;
3483 carry = d + (sum < carry); /* detect overflow */
3484 #endif
3485 }
3487 /* If we run out of 'b' digits before we're actually done, make
3488 sure the carries get propagated upward...
3489 */
3490 used = MP_USED(a);
3491 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3492 while (w && ix < used) {
3493 w = w + *pa;
3494 *pa++ = ACCUM(w);
3495 w = CARRYOUT(w);
3496 ++ix;
3497 }
3498 #else
3499 while (carry && ix < used) {
3500 sum = carry + *pa;
3501 *pa++ = sum;
3502 carry = !sum;
3503 ++ix;
3504 }
3505 #endif
3507 /* If there's an overall carry out, increase precision and include
3508 it. We could have done this initially, but why touch the memory
3509 allocator unless we're sure we have to?
3510 */
3511 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3512 if (w) {
3513 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)
3514 return res;
3516 DIGIT(a, ix) = (mp_digit)w;
3517 }
3518 #else
3519 if (carry) {
3520 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)
3521 return res;
3523 DIGIT(a, used) = carry;
3524 }
3525 #endif
3527 return MP_OKAY;
3528 } /* end s_mp_add() */
3530 /* }}} */
3532 /* Compute c = |a| + |b| */ /* magnitude addition */
3533 mp_err s_mp_add_3arg(const mp_int *a, const mp_int *b, mp_int *c)
3534 {
3535 mp_digit *pa, *pb, *pc;
3536 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3537 mp_word w = 0;
3538 #else
3539 mp_digit sum, carry = 0, d;
3540 #endif
3541 mp_size ix;
3542 mp_size used;
3543 mp_err res;
3545 MP_SIGN(c) = MP_SIGN(a);
3546 if (MP_USED(a) < MP_USED(b)) {
3547 const mp_int *xch = a;
3548 a = b;
3549 b = xch;
3550 }
3552 /* Make sure a has enough precision for the output value */
3553 if (MP_OKAY != (res = s_mp_pad(c, MP_USED(a))))
3554 return res;
3556 /*
3557 Add up all digits up to the precision of b. If b had initially
3558 the same precision as a, or greater, we took care of it by the
3559 exchange step above, so there is no problem. If b had initially
3560 less precision, we'll have to make sure the carry out is duly
3561 propagated upward among the higher-order digits of the sum.
3562 */
3563 pa = MP_DIGITS(a);
3564 pb = MP_DIGITS(b);
3565 pc = MP_DIGITS(c);
3566 used = MP_USED(b);
3567 for (ix = 0; ix < used; ix++) {
3568 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3569 w = w + *pa++ + *pb++;
3570 *pc++ = ACCUM(w);
3571 w = CARRYOUT(w);
3572 #else
3573 d = *pa++;
3574 sum = d + *pb++;
3575 d = (sum < d); /* detect overflow */
3576 *pc++ = sum += carry;
3577 carry = d + (sum < carry); /* detect overflow */
3578 #endif
3579 }
3581 /* If we run out of 'b' digits before we're actually done, make
3582 sure the carries get propagated upward...
3583 */
3584 for (used = MP_USED(a); ix < used; ++ix) {
3585 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3586 w = w + *pa++;
3587 *pc++ = ACCUM(w);
3588 w = CARRYOUT(w);
3589 #else
3590 *pc++ = sum = carry + *pa++;
3591 carry = (sum < carry);
3592 #endif
3593 }
3595 /* If there's an overall carry out, increase precision and include
3596 it. We could have done this initially, but why touch the memory
3597 allocator unless we're sure we have to?
3598 */
3599 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3600 if (w) {
3601 if((res = s_mp_pad(c, used + 1)) != MP_OKAY)
3602 return res;
3604 DIGIT(c, used) = (mp_digit)w;
3605 ++used;
3606 }
3607 #else
3608 if (carry) {
3609 if((res = s_mp_pad(c, used + 1)) != MP_OKAY)
3610 return res;
3612 DIGIT(c, used) = carry;
3613 ++used;
3614 }
3615 #endif
3616 MP_USED(c) = used;
3617 return MP_OKAY;
3618 }
3619 /* {{{ s_mp_add_offset(a, b, offset) */
3621 /* Compute a = |a| + ( |b| * (RADIX ** offset) ) */
3622 mp_err s_mp_add_offset(mp_int *a, mp_int *b, mp_size offset)
3623 {
3624 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3625 mp_word w, k = 0;
3626 #else
3627 mp_digit d, sum, carry = 0;
3628 #endif
3629 mp_size ib;
3630 mp_size ia;
3631 mp_size lim;
3632 mp_err res;
3634 /* Make sure a has enough precision for the output value */
3635 lim = MP_USED(b) + offset;
3636 if((lim > USED(a)) && (res = s_mp_pad(a, lim)) != MP_OKAY)
3637 return res;
3639 /*
3640 Add up all digits up to the precision of b. If b had initially
3641 the same precision as a, or greater, we took care of it by the
3642 padding step above, so there is no problem. If b had initially
3643 less precision, we'll have to make sure the carry out is duly
3644 propagated upward among the higher-order digits of the sum.
3645 */
3646 lim = USED(b);
3647 for(ib = 0, ia = offset; ib < lim; ib++, ia++) {
3648 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3649 w = (mp_word)DIGIT(a, ia) + DIGIT(b, ib) + k;
3650 DIGIT(a, ia) = ACCUM(w);
3651 k = CARRYOUT(w);
3652 #else
3653 d = MP_DIGIT(a, ia);
3654 sum = d + MP_DIGIT(b, ib);
3655 d = (sum < d);
3656 MP_DIGIT(a,ia) = sum += carry;
3657 carry = d + (sum < carry);
3658 #endif
3659 }
3661 /* If we run out of 'b' digits before we're actually done, make
3662 sure the carries get propagated upward...
3663 */
3664 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3665 for (lim = MP_USED(a); k && (ia < lim); ++ia) {
3666 w = (mp_word)DIGIT(a, ia) + k;
3667 DIGIT(a, ia) = ACCUM(w);
3668 k = CARRYOUT(w);
3669 }
3670 #else
3671 for (lim = MP_USED(a); carry && (ia < lim); ++ia) {
3672 d = MP_DIGIT(a, ia);
3673 MP_DIGIT(a,ia) = sum = d + carry;
3674 carry = (sum < d);
3675 }
3676 #endif
3678 /* If there's an overall carry out, increase precision and include
3679 it. We could have done this initially, but why touch the memory
3680 allocator unless we're sure we have to?
3681 */
3682 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_ADD_WORD)
3683 if(k) {
3684 if((res = s_mp_pad(a, USED(a) + 1)) != MP_OKAY)
3685 return res;
3687 DIGIT(a, ia) = (mp_digit)k;
3688 }
3689 #else
3690 if (carry) {
3691 if((res = s_mp_pad(a, lim + 1)) != MP_OKAY)
3692 return res;
3694 DIGIT(a, lim) = carry;
3695 }
3696 #endif
3697 s_mp_clamp(a);
3699 return MP_OKAY;
3701 } /* end s_mp_add_offset() */
3703 /* }}} */
3705 /* {{{ s_mp_sub(a, b) */
3707 /* Compute a = |a| - |b|, assumes |a| >= |b| */
3708 mp_err s_mp_sub(mp_int *a, const mp_int *b) /* magnitude subtract */
3709 {
3710 mp_digit *pa, *pb, *limit;
3711 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3712 mp_sword w = 0;
3713 #else
3714 mp_digit d, diff, borrow = 0;
3715 #endif
3717 /*
3718 Subtract and propagate borrow. Up to the precision of b, this
3719 accounts for the digits of b; after that, we just make sure the
3720 carries get to the right place. This saves having to pad b out to
3721 the precision of a just to make the loops work right...
3722 */
3723 pa = MP_DIGITS(a);
3724 pb = MP_DIGITS(b);
3725 limit = pb + MP_USED(b);
3726 while (pb < limit) {
3727 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3728 w = w + *pa - *pb++;
3729 *pa++ = ACCUM(w);
3730 w >>= MP_DIGIT_BIT;
3731 #else
3732 d = *pa;
3733 diff = d - *pb++;
3734 d = (diff > d); /* detect borrow */
3735 if (borrow && --diff == MP_DIGIT_MAX)
3736 ++d;
3737 *pa++ = diff;
3738 borrow = d;
3739 #endif
3740 }
3741 limit = MP_DIGITS(a) + MP_USED(a);
3742 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3743 while (w && pa < limit) {
3744 w = w + *pa;
3745 *pa++ = ACCUM(w);
3746 w >>= MP_DIGIT_BIT;
3747 }
3748 #else
3749 while (borrow && pa < limit) {
3750 d = *pa;
3751 *pa++ = diff = d - borrow;
3752 borrow = (diff > d);
3753 }
3754 #endif
3756 /* Clobber any leading zeroes we created */
3757 s_mp_clamp(a);
3759 /*
3760 If there was a borrow out, then |b| > |a| in violation
3761 of our input invariant. We've already done the work,
3762 but we'll at least complain about it...
3763 */
3764 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3765 return w ? MP_RANGE : MP_OKAY;
3766 #else
3767 return borrow ? MP_RANGE : MP_OKAY;
3768 #endif
3769 } /* end s_mp_sub() */
3771 /* }}} */
3773 /* Compute c = |a| - |b|, assumes |a| >= |b| */ /* magnitude subtract */
3774 mp_err s_mp_sub_3arg(const mp_int *a, const mp_int *b, mp_int *c)
3775 {
3776 mp_digit *pa, *pb, *pc;
3777 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3778 mp_sword w = 0;
3779 #else
3780 mp_digit d, diff, borrow = 0;
3781 #endif
3782 int ix, limit;
3783 mp_err res;
3785 MP_SIGN(c) = MP_SIGN(a);
3787 /* Make sure a has enough precision for the output value */
3788 if (MP_OKAY != (res = s_mp_pad(c, MP_USED(a))))
3789 return res;
3791 /*
3792 Subtract and propagate borrow. Up to the precision of b, this
3793 accounts for the digits of b; after that, we just make sure the
3794 carries get to the right place. This saves having to pad b out to
3795 the precision of a just to make the loops work right...
3796 */
3797 pa = MP_DIGITS(a);
3798 pb = MP_DIGITS(b);
3799 pc = MP_DIGITS(c);
3800 limit = MP_USED(b);
3801 for (ix = 0; ix < limit; ++ix) {
3802 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3803 w = w + *pa++ - *pb++;
3804 *pc++ = ACCUM(w);
3805 w >>= MP_DIGIT_BIT;
3806 #else
3807 d = *pa++;
3808 diff = d - *pb++;
3809 d = (diff > d);
3810 if (borrow && --diff == MP_DIGIT_MAX)
3811 ++d;
3812 *pc++ = diff;
3813 borrow = d;
3814 #endif
3815 }
3816 for (limit = MP_USED(a); ix < limit; ++ix) {
3817 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3818 w = w + *pa++;
3819 *pc++ = ACCUM(w);
3820 w >>= MP_DIGIT_BIT;
3821 #else
3822 d = *pa++;
3823 *pc++ = diff = d - borrow;
3824 borrow = (diff > d);
3825 #endif
3826 }
3828 /* Clobber any leading zeroes we created */
3829 MP_USED(c) = ix;
3830 s_mp_clamp(c);
3832 /*
3833 If there was a borrow out, then |b| > |a| in violation
3834 of our input invariant. We've already done the work,
3835 but we'll at least complain about it...
3836 */
3837 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_SUB_WORD)
3838 return w ? MP_RANGE : MP_OKAY;
3839 #else
3840 return borrow ? MP_RANGE : MP_OKAY;
3841 #endif
3842 }
3843 /* {{{ s_mp_mul(a, b) */
3845 /* Compute a = |a| * |b| */
3846 mp_err s_mp_mul(mp_int *a, const mp_int *b)
3847 {
3848 return mp_mul(a, b, a);
3849 } /* end s_mp_mul() */
3851 /* }}} */
3853 #if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)
3854 /* This trick works on Sparc V8 CPUs with the Workshop compilers. */
3855 #define MP_MUL_DxD(a, b, Phi, Plo) \
3856 { unsigned long long product = (unsigned long long)a * b; \
3857 Plo = (mp_digit)product; \
3858 Phi = (mp_digit)(product >> MP_DIGIT_BIT); }
3859 #elif defined(OSF1)
3860 #define MP_MUL_DxD(a, b, Phi, Plo) \
3861 { Plo = asm ("mulq %a0, %a1, %v0", a, b);\
3862 Phi = asm ("umulh %a0, %a1, %v0", a, b); }
3863 #else
3864 #define MP_MUL_DxD(a, b, Phi, Plo) \
3865 { mp_digit a0b1, a1b0; \
3866 Plo = (a & MP_HALF_DIGIT_MAX) * (b & MP_HALF_DIGIT_MAX); \
3867 Phi = (a >> MP_HALF_DIGIT_BIT) * (b >> MP_HALF_DIGIT_BIT); \
3868 a0b1 = (a & MP_HALF_DIGIT_MAX) * (b >> MP_HALF_DIGIT_BIT); \
3869 a1b0 = (a >> MP_HALF_DIGIT_BIT) * (b & MP_HALF_DIGIT_MAX); \
3870 a1b0 += a0b1; \
3871 Phi += a1b0 >> MP_HALF_DIGIT_BIT; \
3872 if (a1b0 < a0b1) \
3873 Phi += MP_HALF_RADIX; \
3874 a1b0 <<= MP_HALF_DIGIT_BIT; \
3875 Plo += a1b0; \
3876 if (Plo < a1b0) \
3877 ++Phi; \
3878 }
3879 #endif
3881 #if !defined(MP_ASSEMBLY_MULTIPLY)
3882 /* c = a * b */
3883 void s_mpv_mul_d(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *c)
3884 {
3885 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
3886 mp_digit d = 0;
3888 /* Inner product: Digits of a */
3889 while (a_len--) {
3890 mp_word w = ((mp_word)b * *a++) + d;
3891 *c++ = ACCUM(w);
3892 d = CARRYOUT(w);
3893 }
3894 *c = d;
3895 #else
3896 mp_digit carry = 0;
3897 while (a_len--) {
3898 mp_digit a_i = *a++;
3899 mp_digit a0b0, a1b1;
3901 MP_MUL_DxD(a_i, b, a1b1, a0b0);
3903 a0b0 += carry;
3904 if (a0b0 < carry)
3905 ++a1b1;
3906 *c++ = a0b0;
3907 carry = a1b1;
3908 }
3909 *c = carry;
3910 #endif
3911 }
3913 /* c += a * b */
3914 void s_mpv_mul_d_add(const mp_digit *a, mp_size a_len, mp_digit b,
3915 mp_digit *c)
3916 {
3917 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
3918 mp_digit d = 0;
3920 /* Inner product: Digits of a */
3921 while (a_len--) {
3922 mp_word w = ((mp_word)b * *a++) + *c + d;
3923 *c++ = ACCUM(w);
3924 d = CARRYOUT(w);
3925 }
3926 *c = d;
3927 #else
3928 mp_digit carry = 0;
3929 while (a_len--) {
3930 mp_digit a_i = *a++;
3931 mp_digit a0b0, a1b1;
3933 MP_MUL_DxD(a_i, b, a1b1, a0b0);
3935 a0b0 += carry;
3936 if (a0b0 < carry)
3937 ++a1b1;
3938 a0b0 += a_i = *c;
3939 if (a0b0 < a_i)
3940 ++a1b1;
3941 *c++ = a0b0;
3942 carry = a1b1;
3943 }
3944 *c = carry;
3945 #endif
3946 }
3948 /* Presently, this is only used by the Montgomery arithmetic code. */
3949 /* c += a * b */
3950 void s_mpv_mul_d_add_prop(const mp_digit *a, mp_size a_len, mp_digit b, mp_digit *c)
3951 {
3952 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
3953 mp_digit d = 0;
3955 /* Inner product: Digits of a */
3956 while (a_len--) {
3957 mp_word w = ((mp_word)b * *a++) + *c + d;
3958 *c++ = ACCUM(w);
3959 d = CARRYOUT(w);
3960 }
3962 while (d) {
3963 mp_word w = (mp_word)*c + d;
3964 *c++ = ACCUM(w);
3965 d = CARRYOUT(w);
3966 }
3967 #else
3968 mp_digit carry = 0;
3969 while (a_len--) {
3970 mp_digit a_i = *a++;
3971 mp_digit a0b0, a1b1;
3973 MP_MUL_DxD(a_i, b, a1b1, a0b0);
3975 a0b0 += carry;
3976 if (a0b0 < carry)
3977 ++a1b1;
3979 a0b0 += a_i = *c;
3980 if (a0b0 < a_i)
3981 ++a1b1;
3983 *c++ = a0b0;
3984 carry = a1b1;
3985 }
3986 while (carry) {
3987 mp_digit c_i = *c;
3988 carry += c_i;
3989 *c++ = carry;
3990 carry = carry < c_i;
3991 }
3992 #endif
3993 }
3994 #endif
3996 #if defined(MP_USE_UINT_DIGIT) && defined(MP_USE_LONG_LONG_MULTIPLY)
3997 /* This trick works on Sparc V8 CPUs with the Workshop compilers. */
3998 #define MP_SQR_D(a, Phi, Plo) \
3999 { unsigned long long square = (unsigned long long)a * a; \
4000 Plo = (mp_digit)square; \
4001 Phi = (mp_digit)(square >> MP_DIGIT_BIT); }
4002 #elif defined(OSF1)
4003 #define MP_SQR_D(a, Phi, Plo) \
4004 { Plo = asm ("mulq %a0, %a0, %v0", a);\
4005 Phi = asm ("umulh %a0, %a0, %v0", a); }
4006 #else
4007 #define MP_SQR_D(a, Phi, Plo) \
4008 { mp_digit Pmid; \
4009 Plo = (a & MP_HALF_DIGIT_MAX) * (a & MP_HALF_DIGIT_MAX); \
4010 Phi = (a >> MP_HALF_DIGIT_BIT) * (a >> MP_HALF_DIGIT_BIT); \
4011 Pmid = (a & MP_HALF_DIGIT_MAX) * (a >> MP_HALF_DIGIT_BIT); \
4012 Phi += Pmid >> (MP_HALF_DIGIT_BIT - 1); \
4013 Pmid <<= (MP_HALF_DIGIT_BIT + 1); \
4014 Plo += Pmid; \
4015 if (Plo < Pmid) \
4016 ++Phi; \
4017 }
4018 #endif
4020 #if !defined(MP_ASSEMBLY_SQUARE)
4021 /* Add the squares of the digits of a to the digits of b. */
4022 void s_mpv_sqr_add_prop(const mp_digit *pa, mp_size a_len, mp_digit *ps)
4023 {
4024 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_MUL_WORD)
4025 mp_word w;
4026 mp_digit d;
4027 mp_size ix;
4029 w = 0;
4030 #define ADD_SQUARE(n) \
4031 d = pa[n]; \
4032 w += (d * (mp_word)d) + ps[2*n]; \
4033 ps[2*n] = ACCUM(w); \
4034 w = (w >> DIGIT_BIT) + ps[2*n+1]; \
4035 ps[2*n+1] = ACCUM(w); \
4036 w = (w >> DIGIT_BIT)
4038 for (ix = a_len; ix >= 4; ix -= 4) {
4039 ADD_SQUARE(0);
4040 ADD_SQUARE(1);
4041 ADD_SQUARE(2);
4042 ADD_SQUARE(3);
4043 pa += 4;
4044 ps += 8;
4045 }
4046 if (ix) {
4047 ps += 2*ix;
4048 pa += ix;
4049 switch (ix) {
4050 case 3: ADD_SQUARE(-3); /* FALLTHRU */
4051 case 2: ADD_SQUARE(-2); /* FALLTHRU */
4052 case 1: ADD_SQUARE(-1); /* FALLTHRU */
4053 case 0: break;
4054 }
4055 }
4056 while (w) {
4057 w += *ps;
4058 *ps++ = ACCUM(w);
4059 w = (w >> DIGIT_BIT);
4060 }
4061 #else
4062 mp_digit carry = 0;
4063 while (a_len--) {
4064 mp_digit a_i = *pa++;
4065 mp_digit a0a0, a1a1;
4067 MP_SQR_D(a_i, a1a1, a0a0);
4069 /* here a1a1 and a0a0 constitute a_i ** 2 */
4070 a0a0 += carry;
4071 if (a0a0 < carry)
4072 ++a1a1;
4074 /* now add to ps */
4075 a0a0 += a_i = *ps;
4076 if (a0a0 < a_i)
4077 ++a1a1;
4078 *ps++ = a0a0;
4079 a1a1 += a_i = *ps;
4080 carry = (a1a1 < a_i);
4081 *ps++ = a1a1;
4082 }
4083 while (carry) {
4084 mp_digit s_i = *ps;
4085 carry += s_i;
4086 *ps++ = carry;
4087 carry = carry < s_i;
4088 }
4089 #endif
4090 }
4091 #endif
4093 #if (defined(MP_NO_MP_WORD) || defined(MP_NO_DIV_WORD)) \
4094 && !defined(MP_ASSEMBLY_DIV_2DX1D)
4095 /*
4096 ** Divide 64-bit (Nhi,Nlo) by 32-bit divisor, which must be normalized
4097 ** so its high bit is 1. This code is from NSPR.
4098 */
4099 mp_err s_mpv_div_2dx1d(mp_digit Nhi, mp_digit Nlo, mp_digit divisor,
4100 mp_digit *qp, mp_digit *rp)
4101 {
4102 mp_digit d1, d0, q1, q0;
4103 mp_digit r1, r0, m;
4105 d1 = divisor >> MP_HALF_DIGIT_BIT;
4106 d0 = divisor & MP_HALF_DIGIT_MAX;
4107 r1 = Nhi % d1;
4108 q1 = Nhi / d1;
4109 m = q1 * d0;
4110 r1 = (r1 << MP_HALF_DIGIT_BIT) | (Nlo >> MP_HALF_DIGIT_BIT);
4111 if (r1 < m) {
4112 q1--, r1 += divisor;
4113 if (r1 >= divisor && r1 < m) {
4114 q1--, r1 += divisor;
4115 }
4116 }
4117 r1 -= m;
4118 r0 = r1 % d1;
4119 q0 = r1 / d1;
4120 m = q0 * d0;
4121 r0 = (r0 << MP_HALF_DIGIT_BIT) | (Nlo & MP_HALF_DIGIT_MAX);
4122 if (r0 < m) {
4123 q0--, r0 += divisor;
4124 if (r0 >= divisor && r0 < m) {
4125 q0--, r0 += divisor;
4126 }
4127 }
4128 if (qp)
4129 *qp = (q1 << MP_HALF_DIGIT_BIT) | q0;
4130 if (rp)
4131 *rp = r0 - m;
4132 return MP_OKAY;
4133 }
4134 #endif
4136 #if MP_SQUARE
4137 /* {{{ s_mp_sqr(a) */
4139 mp_err s_mp_sqr(mp_int *a)
4140 {
4141 mp_err res;
4142 mp_int tmp;
4144 if((res = mp_init_size(&tmp, 2 * USED(a))) != MP_OKAY)
4145 return res;
4146 res = mp_sqr(a, &tmp);
4147 if (res == MP_OKAY) {
4148 s_mp_exch(&tmp, a);
4149 }
4150 mp_clear(&tmp);
4151 return res;
4152 }
4154 /* }}} */
4155 #endif
4157 /* {{{ s_mp_div(a, b) */
4159 /*
4160 s_mp_div(a, b)
4162 Compute a = a / b and b = a mod b. Assumes b > a.
4163 */
4165 mp_err s_mp_div(mp_int *rem, /* i: dividend, o: remainder */
4166 mp_int *div, /* i: divisor */
4167 mp_int *quot) /* i: 0; o: quotient */
4168 {
4169 mp_int part, t;
4170 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
4171 mp_word q_msd;
4172 #else
4173 mp_digit q_msd;
4174 #endif
4175 mp_err res;
4176 mp_digit d;
4177 mp_digit div_msd;
4178 int ix;
4180 if(mp_cmp_z(div) == 0)
4181 return MP_RANGE;
4183 DIGITS(&t) = 0;
4184 /* Shortcut if divisor is power of two */
4185 if((ix = s_mp_ispow2(div)) >= 0) {
4186 MP_CHECKOK( mp_copy(rem, quot) );
4187 s_mp_div_2d(quot, (mp_digit)ix);
4188 s_mp_mod_2d(rem, (mp_digit)ix);
4190 return MP_OKAY;
4191 }
4193 MP_SIGN(rem) = ZPOS;
4194 MP_SIGN(div) = ZPOS;
4196 /* A working temporary for division */
4197 MP_CHECKOK( mp_init_size(&t, MP_ALLOC(rem)));
4199 /* Normalize to optimize guessing */
4200 MP_CHECKOK( s_mp_norm(rem, div, &d) );
4202 part = *rem;
4204 /* Perform the division itself...woo! */
4205 MP_USED(quot) = MP_ALLOC(quot);
4207 /* Find a partial substring of rem which is at least div */
4208 /* If we didn't find one, we're finished dividing */
4209 while (MP_USED(rem) > MP_USED(div) || s_mp_cmp(rem, div) >= 0) {
4210 int i;
4211 int unusedRem;
4213 unusedRem = MP_USED(rem) - MP_USED(div);
4214 MP_DIGITS(&part) = MP_DIGITS(rem) + unusedRem;
4215 MP_ALLOC(&part) = MP_ALLOC(rem) - unusedRem;
4216 MP_USED(&part) = MP_USED(div);
4217 if (s_mp_cmp(&part, div) < 0) {
4218 -- unusedRem;
4219 #if MP_ARGCHK == 2
4220 assert(unusedRem >= 0);
4221 #endif
4222 -- MP_DIGITS(&part);
4223 ++ MP_USED(&part);
4224 ++ MP_ALLOC(&part);
4225 }
4227 /* Compute a guess for the next quotient digit */
4228 q_msd = MP_DIGIT(&part, MP_USED(&part) - 1);
4229 div_msd = MP_DIGIT(div, MP_USED(div) - 1);
4230 if (q_msd >= div_msd) {
4231 q_msd = 1;
4232 } else if (MP_USED(&part) > 1) {
4233 #if !defined(MP_NO_MP_WORD) && !defined(MP_NO_DIV_WORD)
4234 q_msd = (q_msd << MP_DIGIT_BIT) | MP_DIGIT(&part, MP_USED(&part) - 2);
4235 q_msd /= div_msd;
4236 if (q_msd == RADIX)
4237 --q_msd;
4238 #else
4239 mp_digit r;
4240 MP_CHECKOK( s_mpv_div_2dx1d(q_msd, MP_DIGIT(&part, MP_USED(&part) - 2),
4241 div_msd, &q_msd, &r) );
4242 #endif
4243 } else {
4244 q_msd = 0;
4245 }
4246 #if MP_ARGCHK == 2
4247 assert(q_msd > 0); /* This case should never occur any more. */
4248 #endif
4249 if (q_msd <= 0)
4250 break;
4252 /* See what that multiplies out to */
4253 mp_copy(div, &t);
4254 MP_CHECKOK( s_mp_mul_d(&t, (mp_digit)q_msd) );
4256 /*
4257 If it's too big, back it off. We should not have to do this
4258 more than once, or, in rare cases, twice. Knuth describes a
4259 method by which this could be reduced to a maximum of once, but
4260 I didn't implement that here.
4261 * When using s_mpv_div_2dx1d, we may have to do this 3 times.
4262 */
4263 for (i = 4; s_mp_cmp(&t, &part) > 0 && i > 0; --i) {
4264 --q_msd;
4265 s_mp_sub(&t, div); /* t -= div */
4266 }
4267 if (i < 0) {
4268 res = MP_RANGE;
4269 goto CLEANUP;
4270 }
4272 /* At this point, q_msd should be the right next digit */
4273 MP_CHECKOK( s_mp_sub(&part, &t) ); /* part -= t */
4274 s_mp_clamp(rem);
4276 /*
4277 Include the digit in the quotient. We allocated enough memory
4278 for any quotient we could ever possibly get, so we should not
4279 have to check for failures here
4280 */
4281 MP_DIGIT(quot, unusedRem) = (mp_digit)q_msd;
4282 }
4284 /* Denormalize remainder */
4285 if (d) {
4286 s_mp_div_2d(rem, d);
4287 }
4289 s_mp_clamp(quot);
4291 CLEANUP:
4292 mp_clear(&t);
4294 return res;
4296 } /* end s_mp_div() */
4299 /* }}} */
4301 /* {{{ s_mp_2expt(a, k) */
4303 mp_err s_mp_2expt(mp_int *a, mp_digit k)
4304 {
4305 mp_err res;
4306 mp_size dig, bit;
4308 dig = k / DIGIT_BIT;
4309 bit = k % DIGIT_BIT;
4311 mp_zero(a);
4312 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY)
4313 return res;
4315 DIGIT(a, dig) |= ((mp_digit)1 << bit);
4317 return MP_OKAY;
4319 } /* end s_mp_2expt() */
4321 /* }}} */
4323 /* {{{ s_mp_reduce(x, m, mu) */
4325 /*
4326 Compute Barrett reduction, x (mod m), given a precomputed value for
4327 mu = b^2k / m, where b = RADIX and k = #digits(m). This should be
4328 faster than straight division, when many reductions by the same
4329 value of m are required (such as in modular exponentiation). This
4330 can nearly halve the time required to do modular exponentiation,
4331 as compared to using the full integer divide to reduce.
4333 This algorithm was derived from the _Handbook of Applied
4334 Cryptography_ by Menezes, Oorschot and VanStone, Ch. 14,
4335 pp. 603-604.
4336 */
4338 mp_err s_mp_reduce(mp_int *x, const mp_int *m, const mp_int *mu)
4339 {
4340 mp_int q;
4341 mp_err res;
4343 if((res = mp_init_copy(&q, x)) != MP_OKAY)
4344 return res;
4346 s_mp_rshd(&q, USED(m) - 1); /* q1 = x / b^(k-1) */
4347 s_mp_mul(&q, mu); /* q2 = q1 * mu */
4348 s_mp_rshd(&q, USED(m) + 1); /* q3 = q2 / b^(k+1) */
4350 /* x = x mod b^(k+1), quick (no division) */
4351 s_mp_mod_2d(x, DIGIT_BIT * (USED(m) + 1));
4353 /* q = q * m mod b^(k+1), quick (no division) */
4354 s_mp_mul(&q, m);
4355 s_mp_mod_2d(&q, DIGIT_BIT * (USED(m) + 1));
4357 /* x = x - q */
4358 if((res = mp_sub(x, &q, x)) != MP_OKAY)
4359 goto CLEANUP;
4361 /* If x < 0, add b^(k+1) to it */
4362 if(mp_cmp_z(x) < 0) {
4363 mp_set(&q, 1);
4364 if((res = s_mp_lshd(&q, USED(m) + 1)) != MP_OKAY)
4365 goto CLEANUP;
4366 if((res = mp_add(x, &q, x)) != MP_OKAY)
4367 goto CLEANUP;
4368 }
4370 /* Back off if it's too big */
4371 while(mp_cmp(x, m) >= 0) {
4372 if((res = s_mp_sub(x, m)) != MP_OKAY)
4373 break;
4374 }
4376 CLEANUP:
4377 mp_clear(&q);
4379 return res;
4381 } /* end s_mp_reduce() */
4383 /* }}} */
4385 /* }}} */
4387 /* {{{ Primitive comparisons */
4389 /* {{{ s_mp_cmp(a, b) */
4391 /* Compare |a| <=> |b|, return 0 if equal, <0 if a<b, >0 if a>b */
4392 int s_mp_cmp(const mp_int *a, const mp_int *b)
4393 {
4394 mp_size used_a = MP_USED(a);
4395 {
4396 mp_size used_b = MP_USED(b);
4398 if (used_a > used_b)
4399 goto IS_GT;
4400 if (used_a < used_b)
4401 goto IS_LT;
4402 }
4403 {
4404 mp_digit *pa, *pb;
4405 mp_digit da = 0, db = 0;
4407 #define CMP_AB(n) if ((da = pa[n]) != (db = pb[n])) goto done
4409 pa = MP_DIGITS(a) + used_a;
4410 pb = MP_DIGITS(b) + used_a;
4411 while (used_a >= 4) {
4412 pa -= 4;
4413 pb -= 4;
4414 used_a -= 4;
4415 CMP_AB(3);
4416 CMP_AB(2);
4417 CMP_AB(1);
4418 CMP_AB(0);
4419 }
4420 while (used_a-- > 0 && ((da = *--pa) == (db = *--pb)))
4421 /* do nothing */;
4422 done:
4423 if (da > db)
4424 goto IS_GT;
4425 if (da < db)
4426 goto IS_LT;
4427 }
4428 return MP_EQ;
4429 IS_LT:
4430 return MP_LT;
4431 IS_GT:
4432 return MP_GT;
4433 } /* end s_mp_cmp() */
4435 /* }}} */
4437 /* {{{ s_mp_cmp_d(a, d) */
4439 /* Compare |a| <=> d, return 0 if equal, <0 if a<d, >0 if a>d */
4440 int s_mp_cmp_d(const mp_int *a, mp_digit d)
4441 {
4442 if(USED(a) > 1)
4443 return MP_GT;
4445 if(DIGIT(a, 0) < d)
4446 return MP_LT;
4447 else if(DIGIT(a, 0) > d)
4448 return MP_GT;
4449 else
4450 return MP_EQ;
4452 } /* end s_mp_cmp_d() */
4454 /* }}} */
4456 /* {{{ s_mp_ispow2(v) */
4458 /*
4459 Returns -1 if the value is not a power of two; otherwise, it returns
4460 k such that v = 2^k, i.e. lg(v).
4461 */
4462 int s_mp_ispow2(const mp_int *v)
4463 {
4464 mp_digit d;
4465 int extra = 0, ix;
4467 ix = MP_USED(v) - 1;
4468 d = MP_DIGIT(v, ix); /* most significant digit of v */
4470 extra = s_mp_ispow2d(d);
4471 if (extra < 0 || ix == 0)
4472 return extra;
4474 while (--ix >= 0) {
4475 if (DIGIT(v, ix) != 0)
4476 return -1; /* not a power of two */
4477 extra += MP_DIGIT_BIT;
4478 }
4480 return extra;
4482 } /* end s_mp_ispow2() */
4484 /* }}} */
4486 /* {{{ s_mp_ispow2d(d) */
4488 int s_mp_ispow2d(mp_digit d)
4489 {
4490 if ((d != 0) && ((d & (d-1)) == 0)) { /* d is a power of 2 */
4491 int pow = 0;
4492 #if defined (MP_USE_UINT_DIGIT)
4493 if (d & 0xffff0000U)
4494 pow += 16;
4495 if (d & 0xff00ff00U)
4496 pow += 8;
4497 if (d & 0xf0f0f0f0U)
4498 pow += 4;
4499 if (d & 0xccccccccU)
4500 pow += 2;
4501 if (d & 0xaaaaaaaaU)
4502 pow += 1;
4503 #elif defined(MP_USE_LONG_LONG_DIGIT)
4504 if (d & 0xffffffff00000000ULL)
4505 pow += 32;
4506 if (d & 0xffff0000ffff0000ULL)
4507 pow += 16;
4508 if (d & 0xff00ff00ff00ff00ULL)
4509 pow += 8;
4510 if (d & 0xf0f0f0f0f0f0f0f0ULL)
4511 pow += 4;
4512 if (d & 0xccccccccccccccccULL)
4513 pow += 2;
4514 if (d & 0xaaaaaaaaaaaaaaaaULL)
4515 pow += 1;
4516 #elif defined(MP_USE_LONG_DIGIT)
4517 if (d & 0xffffffff00000000UL)
4518 pow += 32;
4519 if (d & 0xffff0000ffff0000UL)
4520 pow += 16;
4521 if (d & 0xff00ff00ff00ff00UL)
4522 pow += 8;
4523 if (d & 0xf0f0f0f0f0f0f0f0UL)
4524 pow += 4;
4525 if (d & 0xccccccccccccccccUL)
4526 pow += 2;
4527 if (d & 0xaaaaaaaaaaaaaaaaUL)
4528 pow += 1;
4529 #else
4530 #error "unknown type for mp_digit"
4531 #endif
4532 return pow;
4533 }
4534 return -1;
4536 } /* end s_mp_ispow2d() */
4538 /* }}} */
4540 /* }}} */
4542 /* {{{ Primitive I/O helpers */
4544 /* {{{ s_mp_tovalue(ch, r) */
4546 /*
4547 Convert the given character to its digit value, in the given radix.
4548 If the given character is not understood in the given radix, -1 is
4549 returned. Otherwise the digit's numeric value is returned.
4551 The results will be odd if you use a radix < 2 or > 62, you are
4552 expected to know what you're up to.
4553 */
4554 int s_mp_tovalue(char ch, int r)
4555 {
4556 int val, xch;
4558 if(r > 36)
4559 xch = ch;
4560 else
4561 xch = toupper(ch);
4563 if(isdigit(xch))
4564 val = xch - '0';
4565 else if(isupper(xch))
4566 val = xch - 'A' + 10;
4567 else if(islower(xch))
4568 val = xch - 'a' + 36;
4569 else if(xch == '+')
4570 val = 62;
4571 else if(xch == '/')
4572 val = 63;
4573 else
4574 return -1;
4576 if(val < 0 || val >= r)
4577 return -1;
4579 return val;
4581 } /* end s_mp_tovalue() */
4583 /* }}} */
4585 /* {{{ s_mp_todigit(val, r, low) */
4587 /*
4588 Convert val to a radix-r digit, if possible. If val is out of range
4589 for r, returns zero. Otherwise, returns an ASCII character denoting
4590 the value in the given radix.
4592 The results may be odd if you use a radix < 2 or > 64, you are
4593 expected to know what you're doing.
4594 */
4596 char s_mp_todigit(mp_digit val, int r, int low)
4597 {
4598 char ch;
4600 if(val >= r)
4601 return 0;
4603 ch = s_dmap_1[val];
4605 if(r <= 36 && low)
4606 ch = tolower(ch);
4608 return ch;
4610 } /* end s_mp_todigit() */
4612 /* }}} */
4614 /* {{{ s_mp_outlen(bits, radix) */
4616 /*
4617 Return an estimate for how long a string is needed to hold a radix
4618 r representation of a number with 'bits' significant bits, plus an
4619 extra for a zero terminator (assuming C style strings here)
4620 */
4621 int s_mp_outlen(int bits, int r)
4622 {
4623 return (int)((double)bits * LOG_V_2(r) + 1.5) + 1;
4625 } /* end s_mp_outlen() */
4627 /* }}} */
4629 /* }}} */
4631 /* {{{ mp_read_unsigned_octets(mp, str, len) */
4632 /* mp_read_unsigned_octets(mp, str, len)
4633 Read in a raw value (base 256) into the given mp_int
4634 No sign bit, number is positive. Leading zeros ignored.
4635 */
4637 mp_err
4638 mp_read_unsigned_octets(mp_int *mp, const unsigned char *str, mp_size len)
4639 {
4640 int count;
4641 mp_err res;
4642 mp_digit d;
4644 ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
4646 mp_zero(mp);
4648 count = len % sizeof(mp_digit);
4649 if (count) {
4650 for (d = 0; count-- > 0; --len) {
4651 d = (d << 8) | *str++;
4652 }
4653 MP_DIGIT(mp, 0) = d;
4654 }
4656 /* Read the rest of the digits */
4657 for(; len > 0; len -= sizeof(mp_digit)) {
4658 for (d = 0, count = sizeof(mp_digit); count > 0; --count) {
4659 d = (d << 8) | *str++;
4660 }
4661 if (MP_EQ == mp_cmp_z(mp)) {
4662 if (!d)
4663 continue;
4664 } else {
4665 if((res = s_mp_lshd(mp, 1)) != MP_OKAY)
4666 return res;
4667 }
4668 MP_DIGIT(mp, 0) = d;
4669 }
4670 return MP_OKAY;
4671 } /* end mp_read_unsigned_octets() */
4672 /* }}} */
4674 /* {{{ mp_unsigned_octet_size(mp) */
4675 int
4676 mp_unsigned_octet_size(const mp_int *mp)
4677 {
4678 int bytes;
4679 int ix;
4680 mp_digit d = 0;
4682 ARGCHK(mp != NULL, MP_BADARG);
4683 ARGCHK(MP_ZPOS == SIGN(mp), MP_BADARG);
4685 bytes = (USED(mp) * sizeof(mp_digit));
4687 /* subtract leading zeros. */
4688 /* Iterate over each digit... */
4689 for(ix = USED(mp) - 1; ix >= 0; ix--) {
4690 d = DIGIT(mp, ix);
4691 if (d)
4692 break;
4693 bytes -= sizeof(d);
4694 }
4695 if (!bytes)
4696 return 1;
4698 /* Have MSD, check digit bytes, high order first */
4699 for(ix = sizeof(mp_digit) - 1; ix >= 0; ix--) {
4700 unsigned char x = (unsigned char)(d >> (ix * CHAR_BIT));
4701 if (x)
4702 break;
4703 --bytes;
4704 }
4705 return bytes;
4706 } /* end mp_unsigned_octet_size() */
4707 /* }}} */
4709 /* {{{ mp_to_unsigned_octets(mp, str) */
4710 /* output a buffer of big endian octets no longer than specified. */
4711 mp_err
4712 mp_to_unsigned_octets(const mp_int *mp, unsigned char *str, mp_size maxlen)
4713 {
4714 int ix, pos = 0;
4715 int bytes;
4717 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);
4719 bytes = mp_unsigned_octet_size(mp);
4720 ARGCHK(bytes >= 0 && bytes <= maxlen, MP_BADARG);
4722 /* Iterate over each digit... */
4723 for(ix = USED(mp) - 1; ix >= 0; ix--) {
4724 mp_digit d = DIGIT(mp, ix);
4725 int jx;
4727 /* Unpack digit bytes, high order first */
4728 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {
4729 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));
4730 if (!pos && !x) /* suppress leading zeros */
4731 continue;
4732 str[pos++] = x;
4733 }
4734 }
4735 if (!pos)
4736 str[pos++] = 0;
4737 return pos;
4738 } /* end mp_to_unsigned_octets() */
4739 /* }}} */
4741 /* {{{ mp_to_signed_octets(mp, str) */
4742 /* output a buffer of big endian octets no longer than specified. */
4743 mp_err
4744 mp_to_signed_octets(const mp_int *mp, unsigned char *str, mp_size maxlen)
4745 {
4746 int ix, pos = 0;
4747 int bytes;
4749 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);
4751 bytes = mp_unsigned_octet_size(mp);
4752 ARGCHK(bytes >= 0 && bytes <= maxlen, MP_BADARG);
4754 /* Iterate over each digit... */
4755 for(ix = USED(mp) - 1; ix >= 0; ix--) {
4756 mp_digit d = DIGIT(mp, ix);
4757 int jx;
4759 /* Unpack digit bytes, high order first */
4760 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {
4761 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));
4762 if (!pos) {
4763 if (!x) /* suppress leading zeros */
4764 continue;
4765 if (x & 0x80) { /* add one leading zero to make output positive. */
4766 ARGCHK(bytes + 1 <= maxlen, MP_BADARG);
4767 if (bytes + 1 > maxlen)
4768 return MP_BADARG;
4769 str[pos++] = 0;
4770 }
4771 }
4772 str[pos++] = x;
4773 }
4774 }
4775 if (!pos)
4776 str[pos++] = 0;
4777 return pos;
4778 } /* end mp_to_signed_octets() */
4779 /* }}} */
4781 /* {{{ mp_to_fixlen_octets(mp, str) */
4782 /* output a buffer of big endian octets exactly as long as requested. */
4783 mp_err
4784 mp_to_fixlen_octets(const mp_int *mp, unsigned char *str, mp_size length)
4785 {
4786 int ix, pos = 0;
4787 int bytes;
4789 ARGCHK(mp != NULL && str != NULL && !SIGN(mp), MP_BADARG);
4791 bytes = mp_unsigned_octet_size(mp);
4792 ARGCHK(bytes >= 0 && bytes <= length, MP_BADARG);
4794 /* place any needed leading zeros */
4795 for (;length > bytes; --length) {
4796 *str++ = 0;
4797 }
4799 /* Iterate over each digit... */
4800 for(ix = USED(mp) - 1; ix >= 0; ix--) {
4801 mp_digit d = DIGIT(mp, ix);
4802 int jx;
4804 /* Unpack digit bytes, high order first */
4805 for(jx = sizeof(mp_digit) - 1; jx >= 0; jx--) {
4806 unsigned char x = (unsigned char)(d >> (jx * CHAR_BIT));
4807 if (!pos && !x) /* suppress leading zeros */
4808 continue;
4809 str[pos++] = x;
4810 }
4811 }
4812 if (!pos)
4813 str[pos++] = 0;
4814 return MP_OKAY;
4815 } /* end mp_to_fixlen_octets() */
4816 /* }}} */
4819 /*------------------------------------------------------------------------*/
4820 /* HERE THERE BE DRAGONS */
