Mercurial > pidgin3 / file comparison

comparison: src/protocols/sametime/meanwhile/mpi/mpi.c

src/protocols/sametime/meanwhile/mpi/mpi.c

changeset 12957
9af807a5c9e7
parent 12956
39a4efae983c
child 12958
706645a0b944
equal deleted inserted replaced
12956:39a4efae983c 12957:9af807a5c9e7
1 /*
2 mpi.c
3
4 by Michael J. Fromberger <http://www.dartmouth.edu/~sting/>
5 Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved
6
7 Arbitrary precision integer arithmetic library
8
9 $Id: mpi.c 14563 2005-11-29 23:31:40Z taliesein $
10 */
11
12 #include "mpi.h"
13 #include <stdlib.h>
14 #include <string.h>
15 #include <ctype.h>
16
17 #if MP_DEBUG
18 #include <stdio.h>
19
20 #define DIAG(T,V) {fprintf(stderr,T);mp_print(V,stderr);fputc('\n',stderr);}
21 #else
22 #define DIAG(T,V)
23 #endif
24
25 /*
26 If MP_LOGTAB is not defined, use the math library to compute the
27 logarithms on the fly. Otherwise, use the static table below.
28 Pick which works best for your system.
29 */
30 #if MP_LOGTAB
31
32 /* {{{ s_logv_2[] - log table for 2 in various bases */
33
34 /*
35 A table of the logs of 2 for various bases (the 0 and 1 entries of
36 this table are meaningless and should not be referenced).
37
38 This table is used to compute output lengths for the mp_toradix()
39 function. Since a number n in radix r takes up about log_r(n)
40 digits, we estimate the output size by taking the least integer
41 greater than log_r(n), where:
42
43 log_r(n) = log_2(n) * log_r(2)
44
45 This table, therefore, is a table of log_r(2) for 2 <= r <= 36,
46 which are the output bases supported.
47 */
48
49 #include "logtab.h"
50
51 /* }}} */
52 #define LOG_V_2(R) s_logv_2[(R)]
53
54 #else
55
56 #include <math.h>
57 #define LOG_V_2(R) (log(2.0)/log(R))
58
59 #endif
60
61 /* Default precision for newly created mp_int's */
62 static unsigned int s_mp_defprec = MP_DEFPREC;
63
64 /* {{{ Digit arithmetic macros */
65
66 /*
67 When adding and multiplying digits, the results can be larger than
68 can be contained in an mp_digit. Thus, an mp_word is used. These
69 macros mask off the upper and lower digits of the mp_word (the
70 mp_word may be more than 2 mp_digits wide, but we only concern
71 ourselves with the low-order 2 mp_digits)
72
73 If your mp_word DOES have more than 2 mp_digits, you need to
74 uncomment the first line, and comment out the second.
75 */
76
77 /* #define CARRYOUT(W) (((W)>>DIGIT_BIT)&MP_DIGIT_MAX) */
78 #define CARRYOUT(W) ((W)>>DIGIT_BIT)
79 #define ACCUM(W) ((W)&MP_DIGIT_MAX)
80
81 /* }}} */
82
83 /* {{{ Comparison constants */
84
85 #define MP_LT -1
86 #define MP_EQ 0
87 #define MP_GT 1
88
89 /* }}} */
90
91 /* {{{ Constant strings */
92
93 /* Constant strings returned by mp_strerror() */
94 static const char *mp_err_string[] = {
95 "unknown result code", /* say what? */
96 "boolean true", /* MP_OKAY, MP_YES */
97 "boolean false", /* MP_NO */
98 "out of memory", /* MP_MEM */
99 "argument out of range", /* MP_RANGE */
100 "invalid input parameter", /* MP_BADARG */
101 "result is undefined" /* MP_UNDEF */
102 };
103
104 /* Value to digit maps for radix conversion */
105
106 /* s_dmap_1 - standard digits and letters */
107 static const char *s_dmap_1 =
108 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/";
109
110 #if 0
111 /* s_dmap_2 - base64 ordering for digits */
112 static const char *s_dmap_2 =
113 "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
114 #endif
115
116 /* }}} */
117
118 /* {{{ Static function declarations */
119
120 /*
121 If MP_MACRO is false, these will be defined as actual functions;
122 otherwise, suitable macro definitions will be used. This works
123 around the fact that ANSI C89 doesn't support an 'inline' keyword
124 (although I hear C9x will ... about bloody time). At present, the
125 macro definitions are identical to the function bodies, but they'll
126 expand in place, instead of generating a function call.
127
128 I chose these particular functions to be made into macros because
129 some profiling showed they are called a lot on a typical workload,
130 and yet they are primarily housekeeping.
131 */
132 #if MP_MACRO == 0
133 void s_mp_setz(mp_digit *dp, mp_size count); /* zero digits */
134 void s_mp_copy(mp_digit *sp, mp_digit *dp, mp_size count); /* copy */
135 void *s_mp_alloc(size_t nb, size_t ni); /* general allocator */
136 void s_mp_free(void *ptr); /* general free function */
137 #else
138
139 /* Even if these are defined as macros, we need to respect the settings
140 of the MP_MEMSET and MP_MEMCPY configuration options...
141 */
142 #if MP_MEMSET == 0
143 #define s_mp_setz(dp, count) \
144 {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=0;}
145 #else
146 #define s_mp_setz(dp, count) memset(dp, 0, (count) * sizeof(mp_digit))
147 #endif /* MP_MEMSET */
148
149 #if MP_MEMCPY == 0
150 #define s_mp_copy(sp, dp, count) \
151 {int ix;for(ix=0;ix<(count);ix++)(dp)[ix]=(sp)[ix];}
152 #else
153 #define s_mp_copy(sp, dp, count) memcpy(dp, sp, (count) * sizeof(mp_digit))
154 #endif /* MP_MEMCPY */
155
156 #define s_mp_alloc(nb, ni) calloc(nb, ni)
157 #define s_mp_free(ptr) {if(ptr) free(ptr);}
158 #endif /* MP_MACRO */
159
160 mp_err s_mp_grow(mp_int *mp, mp_size min); /* increase allocated size */
161 mp_err s_mp_pad(mp_int *mp, mp_size min); /* left pad with zeroes */
162
163 void s_mp_clamp(mp_int *mp); /* clip leading zeroes */
164
165 void s_mp_exch(mp_int *a, mp_int *b); /* swap a and b in place */
166
167 mp_err s_mp_lshd(mp_int *mp, mp_size p); /* left-shift by p digits */
168 void s_mp_rshd(mp_int *mp, mp_size p); /* right-shift by p digits */
169 void s_mp_div_2d(mp_int *mp, mp_digit d); /* divide by 2^d in place */
170 void s_mp_mod_2d(mp_int *mp, mp_digit d); /* modulo 2^d in place */
171 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d); /* multiply by 2^d in place*/
172 void s_mp_div_2(mp_int *mp); /* divide by 2 in place */
173 mp_err s_mp_mul_2(mp_int *mp); /* multiply by 2 in place */
174 mp_digit s_mp_norm(mp_int *a, mp_int *b); /* normalize for division */
175 mp_err s_mp_add_d(mp_int *mp, mp_digit d); /* unsigned digit addition */
176 mp_err s_mp_sub_d(mp_int *mp, mp_digit d); /* unsigned digit subtract */
177 mp_err s_mp_mul_d(mp_int *mp, mp_digit d); /* unsigned digit multiply */
178 mp_err s_mp_div_d(mp_int *mp, mp_digit d, mp_digit *r);
179 /* unsigned digit divide */
180 mp_err s_mp_reduce(mp_int *x, mp_int *m, mp_int *mu);
181 /* Barrett reduction */
182 mp_err s_mp_add(mp_int *a, mp_int *b); /* magnitude addition */
183 mp_err s_mp_sub(mp_int *a, mp_int *b); /* magnitude subtract */
184 mp_err s_mp_mul(mp_int *a, mp_int *b); /* magnitude multiply */
185 #if 0
186 void s_mp_kmul(mp_digit *a, mp_digit *b, mp_digit *out, mp_size len);
187 /* multiply buffers in place */
188 #endif
189 #if MP_SQUARE
190 mp_err s_mp_sqr(mp_int *a); /* magnitude square */
191 #else
192 #define s_mp_sqr(a) s_mp_mul(a, a)
193 #endif
194 mp_err s_mp_div(mp_int *a, mp_int *b); /* magnitude divide */
195 mp_err s_mp_2expt(mp_int *a, mp_digit k); /* a = 2^k */
196 int s_mp_cmp(mp_int *a, mp_int *b); /* magnitude comparison */
197 int s_mp_cmp_d(mp_int *a, mp_digit d); /* magnitude digit compare */
198 int s_mp_ispow2(mp_int *v); /* is v a power of 2? */
199 int s_mp_ispow2d(mp_digit d); /* is d a power of 2? */
200
201 int s_mp_tovalue(char ch, int r); /* convert ch to value */
202 char s_mp_todigit(int val, int r, int low); /* convert val to digit */
203 int s_mp_outlen(int bits, int r); /* output length in bytes */
204
205 /* }}} */
206
207 /* {{{ Default precision manipulation */
208
209 unsigned int mp_get_prec(void)
210 {
211 return s_mp_defprec;
212
213 } /* end mp_get_prec() */
214
215 void mp_set_prec(unsigned int prec)
216 {
217 if(prec == 0)
218 s_mp_defprec = MP_DEFPREC;
219 else
220 s_mp_defprec = prec;
221
222 } /* end mp_set_prec() */
223
224 /* }}} */
225
226 /*------------------------------------------------------------------------*/
227 /* {{{ mp_init(mp) */
228
229 /*
230 mp_init(mp)
231
232 Initialize a new zero-valued mp_int. Returns MP_OKAY if successful,
233 MP_MEM if memory could not be allocated for the structure.
234 */
235
236 mp_err mp_init(mp_int *mp)
237 {
238 return mp_init_size(mp, s_mp_defprec);
239
240 } /* end mp_init() */
241
242 /* }}} */
243
244 /* {{{ mp_init_array(mp[], count) */
245
246 mp_err mp_init_array(mp_int mp[], int count)
247 {
248 mp_err res;
249 int pos;
250
251 ARGCHK(mp !=NULL && count > 0, MP_BADARG);
252
253 for(pos = 0; pos < count; ++pos) {
254 if((res = mp_init(&mp[pos])) != MP_OKAY)
255 goto CLEANUP;
256 }
257
258 return MP_OKAY;
259
260 CLEANUP:
261 while(--pos >= 0)
262 mp_clear(&mp[pos]);
263
264 return res;
265
266 } /* end mp_init_array() */
267
268 /* }}} */
269
270 /* {{{ mp_init_size(mp, prec) */
271
272 /*
273 mp_init_size(mp, prec)
274
275 Initialize a new zero-valued mp_int with at least the given
276 precision; returns MP_OKAY if successful, or MP_MEM if memory could
277 not be allocated for the structure.
278 */
279
280 mp_err mp_init_size(mp_int *mp, mp_size prec)
281 {
282 ARGCHK(mp != NULL && prec > 0, MP_BADARG);
283
284 if((DIGITS(mp) = s_mp_alloc(prec, sizeof(mp_digit))) == NULL)
285 return MP_MEM;
286
287 SIGN(mp) = MP_ZPOS;
288 USED(mp) = 1;
289 ALLOC(mp) = prec;
290
291 return MP_OKAY;
292
293 } /* end mp_init_size() */
294
295 /* }}} */
296
297 /* {{{ mp_init_copy(mp, from) */
298
299 /*
300 mp_init_copy(mp, from)
301
302 Initialize mp as an exact copy of from. Returns MP_OKAY if
303 successful, MP_MEM if memory could not be allocated for the new
304 structure.
305 */
306
307 mp_err mp_init_copy(mp_int *mp, mp_int *from)
308 {
309 ARGCHK(mp != NULL && from != NULL, MP_BADARG);
310
311 if(mp == from)
312 return MP_OKAY;
313
314 if((DIGITS(mp) = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL)
315 return MP_MEM;
316
317 s_mp_copy(DIGITS(from), DIGITS(mp), USED(from));
318 USED(mp) = USED(from);
319 ALLOC(mp) = USED(from);
320 SIGN(mp) = SIGN(from);
321
322 return MP_OKAY;
323
324 } /* end mp_init_copy() */
325
326 /* }}} */
327
328 /* {{{ mp_copy(from, to) */
329
330 /*
331 mp_copy(from, to)
332
333 Copies the mp_int 'from' to the mp_int 'to'. It is presumed that
334 'to' has already been initialized (if not, use mp_init_copy()
335 instead). If 'from' and 'to' are identical, nothing happens.
336 */
337
338 mp_err mp_copy(mp_int *from, mp_int *to)
339 {
340 ARGCHK(from != NULL && to != NULL, MP_BADARG);
341
342 if(from == to)
343 return MP_OKAY;
344
345 { /* copy */
346 mp_digit *tmp;
347
348 /*
349 If the allocated buffer in 'to' already has enough space to hold
350 all the used digits of 'from', we'll re-use it to avoid hitting
351 the memory allocater more than necessary; otherwise, we'd have
352 to grow anyway, so we just allocate a hunk and make the copy as
353 usual
354 */
355 if(ALLOC(to) >= USED(from)) {
356 s_mp_setz(DIGITS(to) + USED(from), ALLOC(to) - USED(from));
357 s_mp_copy(DIGITS(from), DIGITS(to), USED(from));
358
359 } else {
360 if((tmp = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL)
361 return MP_MEM;
362
363 s_mp_copy(DIGITS(from), tmp, USED(from));
364
365 if(DIGITS(to) != NULL) {
366 #if MP_CRYPTO
367 s_mp_setz(DIGITS(to), ALLOC(to));
368 #endif
369 s_mp_free(DIGITS(to));
370 }
371
372 DIGITS(to) = tmp;
373 ALLOC(to) = USED(from);
374 }
375
376 /* Copy the precision and sign from the original */
377 USED(to) = USED(from);
378 SIGN(to) = SIGN(from);
379 } /* end copy */
380
381 return MP_OKAY;
382
383 } /* end mp_copy() */
384
385 /* }}} */
386
387 /* {{{ mp_exch(mp1, mp2) */
388
389 /*
390 mp_exch(mp1, mp2)
391
392 Exchange mp1 and mp2 without allocating any intermediate memory
393 (well, unless you count the stack space needed for this call and the
394 locals it creates...). This cannot fail.
395 */
396
397 void mp_exch(mp_int *mp1, mp_int *mp2)
398 {
399 #if MP_ARGCHK == 2
400 assert(mp1 != NULL && mp2 != NULL);
401 #else
402 if(mp1 == NULL || mp2 == NULL)
403 return;
404 #endif
405
406 s_mp_exch(mp1, mp2);
407
408 } /* end mp_exch() */
409
410 /* }}} */
411
412 /* {{{ mp_clear(mp) */
413
414 /*
415 mp_clear(mp)
416
417 Release the storage used by an mp_int, and void its fields so that
418 if someone calls mp_clear() again for the same int later, we won't
419 get tollchocked.
420 */
421
422 void mp_clear(mp_int *mp)
423 {
424 if(mp == NULL)
425 return;
426
427 if(DIGITS(mp) != NULL) {
428 #if MP_CRYPTO
429 s_mp_setz(DIGITS(mp), ALLOC(mp));
430 #endif
431 s_mp_free(DIGITS(mp));
432 DIGITS(mp) = NULL;
433 }
434
435 USED(mp) = 0;
436 ALLOC(mp) = 0;
437
438 } /* end mp_clear() */
439
440 /* }}} */
441
442 /* {{{ mp_clear_array(mp[], count) */
443
444 void mp_clear_array(mp_int mp[], int count)
445 {
446 ARGCHK(mp != NULL && count > 0, MP_BADARG);
447
448 while(--count >= 0)
449 mp_clear(&mp[count]);
450
451 } /* end mp_clear_array() */
452
453 /* }}} */
454
455 /* {{{ mp_zero(mp) */
456
457 /*
458 mp_zero(mp)
459
460 Set mp to zero. Does not change the allocated size of the structure,
461 and therefore cannot fail (except on a bad argument, which we ignore)
462 */
463 void mp_zero(mp_int *mp)
464 {
465 if(mp == NULL)
466 return;
467
468 s_mp_setz(DIGITS(mp), ALLOC(mp));
469 USED(mp) = 1;
470 SIGN(mp) = MP_ZPOS;
471
472 } /* end mp_zero() */
473
474 /* }}} */
475
476 /* {{{ mp_set(mp, d) */
477
478 void mp_set(mp_int *mp, mp_digit d)
479 {
480 if(mp == NULL)
481 return;
482
483 mp_zero(mp);
484 DIGIT(mp, 0) = d;
485
486 } /* end mp_set() */
487
488 /* }}} */
489
490 /* {{{ mp_set_int(mp, z) */
491
492 mp_err mp_set_int(mp_int *mp, long z)
493 {
494 int ix;
495 unsigned long v = abs(z);
496 mp_err res;
497
498 ARGCHK(mp != NULL, MP_BADARG);
499
500 mp_zero(mp);
501 if(z == 0)
502 return MP_OKAY; /* shortcut for zero */
503
504 for(ix = sizeof(long) - 1; ix >= 0; ix--) {
505
506 if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
507 return res;
508
509 res = s_mp_add_d(mp,
510 (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX));
511 if(res != MP_OKAY)
512 return res;
513
514 }
515
516 if(z < 0)
517 SIGN(mp) = MP_NEG;
518
519 return MP_OKAY;
520
521 } /* end mp_set_int() */
522
523 /* }}} */
524
525 /*------------------------------------------------------------------------*/
526 /* {{{ Digit arithmetic */
527
528 /* {{{ mp_add_d(a, d, b) */
529
530 /*
531 mp_add_d(a, d, b)
532
533 Compute the sum b = a + d, for a single digit d. Respects the sign of
534 its primary addend (single digits are unsigned anyway).
535 */
536
537 mp_err mp_add_d(mp_int *a, mp_digit d, mp_int *b)
538 {
539 mp_err res = MP_OKAY;
540
541 ARGCHK(a != NULL && b != NULL, MP_BADARG);
542
543 if((res = mp_copy(a, b)) != MP_OKAY)
544 return res;
545
546 if(SIGN(b) == MP_ZPOS) {
547 res = s_mp_add_d(b, d);
548 } else if(s_mp_cmp_d(b, d) >= 0) {
549 res = s_mp_sub_d(b, d);
550 } else {
551 SIGN(b) = MP_ZPOS;
552
553 DIGIT(b, 0) = d - DIGIT(b, 0);
554 }
555
556 return res;
557
558 } /* end mp_add_d() */
559
560 /* }}} */
561
562 /* {{{ mp_sub_d(a, d, b) */
563
564 /*
565 mp_sub_d(a, d, b)
566
567 Compute the difference b = a - d, for a single digit d. Respects the
568 sign of its subtrahend (single digits are unsigned anyway).
569 */
570
571 mp_err mp_sub_d(mp_int *a, mp_digit d, mp_int *b)
572 {
573 mp_err res;
574
575 ARGCHK(a != NULL && b != NULL, MP_BADARG);
576
577 if((res = mp_copy(a, b)) != MP_OKAY)
578 return res;
579
580 if(SIGN(b) == MP_NEG) {
581 if((res = s_mp_add_d(b, d)) != MP_OKAY)
582 return res;
583
584 } else if(s_mp_cmp_d(b, d) >= 0) {
585 if((res = s_mp_sub_d(b, d)) != MP_OKAY)
586 return res;
587
588 } else {
589 mp_neg(b, b);
590
591 DIGIT(b, 0) = d - DIGIT(b, 0);
592 SIGN(b) = MP_NEG;
593 }
594
595 if(s_mp_cmp_d(b, 0) == 0)
596 SIGN(b) = MP_ZPOS;
597
598 return MP_OKAY;
599
600 } /* end mp_sub_d() */
601
602 /* }}} */
603
604 /* {{{ mp_mul_d(a, d, b) */
605
606 /*
607 mp_mul_d(a, d, b)
608
609 Compute the product b = a * d, for a single digit d. Respects the sign
610 of its multiplicand (single digits are unsigned anyway)
611 */
612
613 mp_err mp_mul_d(mp_int *a, mp_digit d, mp_int *b)
614 {
615 mp_err res;
616
617 ARGCHK(a != NULL && b != NULL, MP_BADARG);
618
619 if(d == 0) {
620 mp_zero(b);
621 return MP_OKAY;
622 }
623
624 if((res = mp_copy(a, b)) != MP_OKAY)
625 return res;
626
627 res = s_mp_mul_d(b, d);
628
629 return res;
630
631 } /* end mp_mul_d() */
632
633 /* }}} */
634
635 /* {{{ mp_mul_2(a, c) */
636
637 mp_err mp_mul_2(mp_int *a, mp_int *c)
638 {
639 mp_err res;
640
641 ARGCHK(a != NULL && c != NULL, MP_BADARG);
642
643 if((res = mp_copy(a, c)) != MP_OKAY)
644 return res;
645
646 return s_mp_mul_2(c);
647
648 } /* end mp_mul_2() */
649
650 /* }}} */
651
652 /* {{{ mp_div_d(a, d, q, r) */
653
654 /*
655 mp_div_d(a, d, q, r)
656
657 Compute the quotient q = a / d and remainder r = a mod d, for a
658 single digit d. Respects the sign of its divisor (single digits are
659 unsigned anyway).
660 */
661
662 mp_err mp_div_d(mp_int *a, mp_digit d, mp_int *q, mp_digit *r)
663 {
664 mp_err res;
665 mp_digit rem;
666 int pow;
667
668 ARGCHK(a != NULL, MP_BADARG);
669
670 if(d == 0)
671 return MP_RANGE;
672
673 /* Shortcut for powers of two ... */
674 if((pow = s_mp_ispow2d(d)) >= 0) {
675 mp_digit mask;
676
677 mask = (1 << pow) - 1;
678 rem = DIGIT(a, 0) & mask;
679
680 if(q) {
681 mp_copy(a, q);
682 s_mp_div_2d(q, pow);
683 }
684
685 if(r)
686 *r = rem;
687
688 return MP_OKAY;
689 }
690
691 /*
692 If the quotient is actually going to be returned, we'll try to
693 avoid hitting the memory allocator by copying the dividend into it
694 and doing the division there. This can't be any _worse_ than
695 always copying, and will sometimes be better (since it won't make
696 another copy)
697
698 If it's not going to be returned, we need to allocate a temporary
699 to hold the quotient, which will just be discarded.
700 */
701 if(q) {
702 if((res = mp_copy(a, q)) != MP_OKAY)
703 return res;
704
705 res = s_mp_div_d(q, d, &rem);
706 if(s_mp_cmp_d(q, 0) == MP_EQ)
707 SIGN(q) = MP_ZPOS;
708
709 } else {
710 mp_int qp;
711
712 if((res = mp_init_copy(&qp, a)) != MP_OKAY)
713 return res;
714
715 res = s_mp_div_d(&qp, d, &rem);
716 if(s_mp_cmp_d(&qp, 0) == 0)
717 SIGN(&qp) = MP_ZPOS;
718
719 mp_clear(&qp);
720 }
721
722 if(r)
723 *r = rem;
724
725 return res;
726
727 } /* end mp_div_d() */
728
729 /* }}} */
730
731 /* {{{ mp_div_2(a, c) */
732
733 /*
734 mp_div_2(a, c)
735
736 Compute c = a / 2, disregarding the remainder.
737 */
738
739 mp_err mp_div_2(mp_int *a, mp_int *c)
740 {
741 mp_err res;
742
743 ARGCHK(a != NULL && c != NULL, MP_BADARG);
744
745 if((res = mp_copy(a, c)) != MP_OKAY)
746 return res;
747
748 s_mp_div_2(c);
749
750 return MP_OKAY;
751
752 } /* end mp_div_2() */
753
754 /* }}} */
755
756 /* {{{ mp_expt_d(a, d, b) */
757
758 mp_err mp_expt_d(mp_int *a, mp_digit d, mp_int *c)
759 {
760 mp_int s, x;
761 mp_err res;
762 mp_sign cs = MP_ZPOS;
763
764 ARGCHK(a != NULL && c != NULL, MP_BADARG);
765
766 if((res = mp_init(&s)) != MP_OKAY)
767 return res;
768 if((res = mp_init_copy(&x, a)) != MP_OKAY)
769 goto X;
770
771 DIGIT(&s, 0) = 1;
772
773 if((d % 2) == 1)
774 cs = SIGN(a);
775
776 while(d != 0) {
777 if(d & 1) {
778 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
779 goto CLEANUP;
780 }
781
782 d >>= 1;
783
784 if((res = s_mp_sqr(&x)) != MP_OKAY)
785 goto CLEANUP;
786 }
787
788 SIGN(&s) = cs;
789
790 s_mp_exch(&s, c);
791
792 CLEANUP:
793 mp_clear(&x);
794 X:
795 mp_clear(&s);
796
797 return res;
798
799 } /* end mp_expt_d() */
800
801 /* }}} */
802
803 /* }}} */
804
805 /*------------------------------------------------------------------------*/
806 /* {{{ Full arithmetic */
807
808 /* {{{ mp_abs(a, b) */
809
810 /*
811 mp_abs(a, b)
812
813 Compute b = |a|. 'a' and 'b' may be identical.
814 */
815
816 mp_err mp_abs(mp_int *a, mp_int *b)
817 {
818 mp_err res;
819
820 ARGCHK(a != NULL && b != NULL, MP_BADARG);
821
822 if((res = mp_copy(a, b)) != MP_OKAY)
823 return res;
824
825 SIGN(b) = MP_ZPOS;
826
827 return MP_OKAY;
828
829 } /* end mp_abs() */
830
831 /* }}} */
832
833 /* {{{ mp_neg(a, b) */
834
835 /*
836 mp_neg(a, b)
837
838 Compute b = -a. 'a' and 'b' may be identical.
839 */
840
841 mp_err mp_neg(mp_int *a, mp_int *b)
842 {
843 mp_err res;
844
845 ARGCHK(a != NULL && b != NULL, MP_BADARG);
846
847 if((res = mp_copy(a, b)) != MP_OKAY)
848 return res;
849
850 if(s_mp_cmp_d(b, 0) == MP_EQ)
851 SIGN(b) = MP_ZPOS;
852 else
853 SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG;
854
855 return MP_OKAY;
856
857 } /* end mp_neg() */
858
859 /* }}} */
860
861 /* {{{ mp_add(a, b, c) */
862
863 /*
864 mp_add(a, b, c)
865
866 Compute c = a + b. All parameters may be identical.
867 */
868
869 mp_err mp_add(mp_int *a, mp_int *b, mp_int *c)
870 {
871 mp_err res;
872 int cmp;
873
874 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
875
876 if(SIGN(a) == SIGN(b)) { /* same sign: add values, keep sign */
877
878 /* Commutativity of addition lets us do this in either order,
879 so we avoid having to use a temporary even if the result
880 is supposed to replace the output
881 */
882 if(c == b) {
883 if((res = s_mp_add(c, a)) != MP_OKAY)
884 return res;
885 } else {
886 if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
887 return res;
888
889 if((res = s_mp_add(c, b)) != MP_OKAY)
890 return res;
891 }
892
893 } else if((cmp = s_mp_cmp(a, b)) > 0) { /* different sign: a > b */
894
895 /* If the output is going to be clobbered, we will use a temporary
896 variable; otherwise, we'll do it without touching the memory
897 allocator at all, if possible
898 */
899 if(c == b) {
900 mp_int tmp;
901
902 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
903 return res;
904 if((res = s_mp_sub(&tmp, b)) != MP_OKAY) {
905 mp_clear(&tmp);
906 return res;
907 }
908
909 s_mp_exch(&tmp, c);
910 mp_clear(&tmp);
911
912 } else {
913
914 if(c != a && (res = mp_copy(a, c)) != MP_OKAY)
915 return res;
916 if((res = s_mp_sub(c, b)) != MP_OKAY)
917 return res;
918
919 }
920
921 } else if(cmp == 0) { /* different sign, a == b */
922
923 mp_zero(c);
924 return MP_OKAY;
925
926 } else { /* different sign: a < b */
927
928 /* See above... */
929 if(c == a) {
930 mp_int tmp;
931
932 if((res = mp_init_copy(&tmp, b)) != MP_OKAY)
933 return res;
934 if((res = s_mp_sub(&tmp, a)) != MP_OKAY) {
935 mp_clear(&tmp);
936 return res;
937 }
938
939 s_mp_exch(&tmp, c);
940 mp_clear(&tmp);
941
942 } else {
943
944 if(c != b && (res = mp_copy(b, c)) != MP_OKAY)
945 return res;
946 if((res = s_mp_sub(c, a)) != MP_OKAY)
947 return res;
948
949 }
950 }
951
952 if(USED(c) == 1 && DIGIT(c, 0) == 0)
953 SIGN(c) = MP_ZPOS;
954
955 return MP_OKAY;
956
957 } /* end mp_add() */
958
959 /* }}} */
960
961 /* {{{ mp_sub(a, b, c) */
962
963 /*
964 mp_sub(a, b, c)
965
966 Compute c = a - b. All parameters may be identical.
967 */
968
969 mp_err mp_sub(mp_int *a, mp_int *b, mp_int *c)
970 {
971 mp_err res;
972 int cmp;
973
974 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
975
976 if(SIGN(a) != SIGN(b)) {
977 if(c == a) {
978 if((res = s_mp_add(c, b)) != MP_OKAY)
979 return res;
980 } else {
981 if(c != b && ((res = mp_copy(b, c)) != MP_OKAY))
982 return res;
983 if((res = s_mp_add(c, a)) != MP_OKAY)
984 return res;
985 SIGN(c) = SIGN(a);
986 }
987
988 } else if((cmp = s_mp_cmp(a, b)) > 0) { /* Same sign, a > b */
989 if(c == b) {
990 mp_int tmp;
991
992 if((res = mp_init_copy(&tmp, a)) != MP_OKAY)
993 return res;
994 if((res = s_mp_sub(&tmp, b)) != MP_OKAY) {
995 mp_clear(&tmp);
996 return res;
997 }
998 s_mp_exch(&tmp, c);
999 mp_clear(&tmp);
1000
1001 } else {
1002 if(c != a && ((res = mp_copy(a, c)) != MP_OKAY))
1003 return res;
1004
1005 if((res = s_mp_sub(c, b)) != MP_OKAY)
1006 return res;
1007 }
1008
1009 } else if(cmp == 0) { /* Same sign, equal magnitude */
1010 mp_zero(c);
1011 return MP_OKAY;
1012
1013 } else { /* Same sign, b > a */
1014 if(c == a) {
1015 mp_int tmp;
1016
1017 if((res = mp_init_copy(&tmp, b)) != MP_OKAY)
1018 return res;
1019
1020 if((res = s_mp_sub(&tmp, a)) != MP_OKAY) {
1021 mp_clear(&tmp);
1022 return res;
1023 }
1024 s_mp_exch(&tmp, c);
1025 mp_clear(&tmp);
1026
1027 } else {
1028 if(c != b && ((res = mp_copy(b, c)) != MP_OKAY))
1029 return res;
1030
1031 if((res = s_mp_sub(c, a)) != MP_OKAY)
1032 return res;
1033 }
1034
1035 SIGN(c) = !SIGN(b);
1036 }
1037
1038 if(USED(c) == 1 && DIGIT(c, 0) == 0)
1039 SIGN(c) = MP_ZPOS;
1040
1041 return MP_OKAY;
1042
1043 } /* end mp_sub() */
1044
1045 /* }}} */
1046
1047 /* {{{ mp_mul(a, b, c) */
1048
1049 /*
1050 mp_mul(a, b, c)
1051
1052 Compute c = a * b. All parameters may be identical.
1053 */
1054
1055 mp_err mp_mul(mp_int *a, mp_int *b, mp_int *c)
1056 {
1057 mp_err res;
1058 mp_sign sgn;
1059
1060 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1061
1062 sgn = (SIGN(a) == SIGN(b)) ? MP_ZPOS : MP_NEG;
1063
1064 if(c == b) {
1065 if((res = s_mp_mul(c, a)) != MP_OKAY)
1066 return res;
1067
1068 } else {
1069 if((res = mp_copy(a, c)) != MP_OKAY)
1070 return res;
1071
1072 if((res = s_mp_mul(c, b)) != MP_OKAY)
1073 return res;
1074 }
1075
1076 if(sgn == MP_ZPOS || s_mp_cmp_d(c, 0) == MP_EQ)
1077 SIGN(c) = MP_ZPOS;
1078 else
1079 SIGN(c) = sgn;
1080
1081 return MP_OKAY;
1082
1083 } /* end mp_mul() */
1084
1085 /* }}} */
1086
1087 /* {{{ mp_mul_2d(a, d, c) */
1088
1089 /*
1090 mp_mul_2d(a, d, c)
1091
1092 Compute c = a * 2^d. a may be the same as c.
1093 */
1094
1095 mp_err mp_mul_2d(mp_int *a, mp_digit d, mp_int *c)
1096 {
1097 mp_err res;
1098
1099 ARGCHK(a != NULL && c != NULL, MP_BADARG);
1100
1101 if((res = mp_copy(a, c)) != MP_OKAY)
1102 return res;
1103
1104 if(d == 0)
1105 return MP_OKAY;
1106
1107 return s_mp_mul_2d(c, d);
1108
1109 } /* end mp_mul() */
1110
1111 /* }}} */
1112
1113 /* {{{ mp_sqr(a, b) */
1114
1115 #if MP_SQUARE
1116 mp_err mp_sqr(mp_int *a, mp_int *b)
1117 {
1118 mp_err res;
1119
1120 ARGCHK(a != NULL && b != NULL, MP_BADARG);
1121
1122 if((res = mp_copy(a, b)) != MP_OKAY)
1123 return res;
1124
1125 if((res = s_mp_sqr(b)) != MP_OKAY)
1126 return res;
1127
1128 SIGN(b) = MP_ZPOS;
1129
1130 return MP_OKAY;
1131
1132 } /* end mp_sqr() */
1133 #endif
1134
1135 /* }}} */
1136
1137 /* {{{ mp_div(a, b, q, r) */
1138
1139 /*
1140 mp_div(a, b, q, r)
1141
1142 Compute q = a / b and r = a mod b. Input parameters may be re-used
1143 as output parameters. If q or r is NULL, that portion of the
1144 computation will be discarded (although it will still be computed)
1145
1146 Pay no attention to the hacker behind the curtain.
1147 */
1148
1149 mp_err mp_div(mp_int *a, mp_int *b, mp_int *q, mp_int *r)
1150 {
1151 mp_err res;
1152 mp_int qtmp, rtmp;
1153 int cmp;
1154
1155 ARGCHK(a != NULL && b != NULL, MP_BADARG);
1156
1157 if(mp_cmp_z(b) == MP_EQ)
1158 return MP_RANGE;
1159
1160 /* If a <= b, we can compute the solution without division, and
1161 avoid any memory allocation
1162 */
1163 if((cmp = s_mp_cmp(a, b)) < 0) {
1164 if(r) {
1165 if((res = mp_copy(a, r)) != MP_OKAY)
1166 return res;
1167 }
1168
1169 if(q)
1170 mp_zero(q);
1171
1172 return MP_OKAY;
1173
1174 } else if(cmp == 0) {
1175
1176 /* Set quotient to 1, with appropriate sign */
1177 if(q) {
1178 int qneg = (SIGN(a) != SIGN(b));
1179
1180 mp_set(q, 1);
1181 if(qneg)
1182 SIGN(q) = MP_NEG;
1183 }
1184
1185 if(r)
1186 mp_zero(r);
1187
1188 return MP_OKAY;
1189 }
1190
1191 /* If we get here, it means we actually have to do some division */
1192
1193 /* Set up some temporaries... */
1194 if((res = mp_init_copy(&qtmp, a)) != MP_OKAY)
1195 return res;
1196 if((res = mp_init_copy(&rtmp, b)) != MP_OKAY)
1197 goto CLEANUP;
1198
1199 if((res = s_mp_div(&qtmp, &rtmp)) != MP_OKAY)
1200 goto CLEANUP;
1201
1202 /* Compute the signs for the output */
1203 SIGN(&rtmp) = SIGN(a); /* Sr = Sa */
1204 if(SIGN(a) == SIGN(b))
1205 SIGN(&qtmp) = MP_ZPOS; /* Sq = MP_ZPOS if Sa = Sb */
1206 else
1207 SIGN(&qtmp) = MP_NEG; /* Sq = MP_NEG if Sa != Sb */
1208
1209 if(s_mp_cmp_d(&qtmp, 0) == MP_EQ)
1210 SIGN(&qtmp) = MP_ZPOS;
1211 if(s_mp_cmp_d(&rtmp, 0) == MP_EQ)
1212 SIGN(&rtmp) = MP_ZPOS;
1213
1214 /* Copy output, if it is needed */
1215 if(q)
1216 s_mp_exch(&qtmp, q);
1217
1218 if(r)
1219 s_mp_exch(&rtmp, r);
1220
1221 CLEANUP:
1222 mp_clear(&rtmp);
1223 mp_clear(&qtmp);
1224
1225 return res;
1226
1227 } /* end mp_div() */
1228
1229 /* }}} */
1230
1231 /* {{{ mp_div_2d(a, d, q, r) */
1232
1233 mp_err mp_div_2d(mp_int *a, mp_digit d, mp_int *q, mp_int *r)
1234 {
1235 mp_err res;
1236
1237 ARGCHK(a != NULL, MP_BADARG);
1238
1239 if(q) {
1240 if((res = mp_copy(a, q)) != MP_OKAY)
1241 return res;
1242
1243 s_mp_div_2d(q, d);
1244 }
1245
1246 if(r) {
1247 if((res = mp_copy(a, r)) != MP_OKAY)
1248 return res;
1249
1250 s_mp_mod_2d(r, d);
1251 }
1252
1253 return MP_OKAY;
1254
1255 } /* end mp_div_2d() */
1256
1257 /* }}} */
1258
1259 /* {{{ mp_expt(a, b, c) */
1260
1261 /*
1262 mp_expt(a, b, c)
1263
1264 Compute c = a ** b, that is, raise a to the b power. Uses a
1265 standard iterative square-and-multiply technique.
1266 */
1267
1268 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c)
1269 {
1270 mp_int s, x;
1271 mp_err res;
1272 mp_digit d;
1273 int dig, bit;
1274
1275 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1276
1277 if(mp_cmp_z(b) < 0)
1278 return MP_RANGE;
1279
1280 if((res = mp_init(&s)) != MP_OKAY)
1281 return res;
1282
1283 mp_set(&s, 1);
1284
1285 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1286 goto X;
1287
1288 /* Loop over low-order digits in ascending order */
1289 for(dig = 0; dig < (USED(b) - 1); dig++) {
1290 d = DIGIT(b, dig);
1291
1292 /* Loop over bits of each non-maximal digit */
1293 for(bit = 0; bit < DIGIT_BIT; bit++) {
1294 if(d & 1) {
1295 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1296 goto CLEANUP;
1297 }
1298
1299 d >>= 1;
1300
1301 if((res = s_mp_sqr(&x)) != MP_OKAY)
1302 goto CLEANUP;
1303 }
1304 }
1305
1306 /* Consider now the last digit... */
1307 d = DIGIT(b, dig);
1308
1309 while(d) {
1310 if(d & 1) {
1311 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1312 goto CLEANUP;
1313 }
1314
1315 d >>= 1;
1316
1317 if((res = s_mp_sqr(&x)) != MP_OKAY)
1318 goto CLEANUP;
1319 }
1320
1321 if(mp_iseven(b))
1322 SIGN(&s) = SIGN(a);
1323
1324 res = mp_copy(&s, c);
1325
1326 CLEANUP:
1327 mp_clear(&x);
1328 X:
1329 mp_clear(&s);
1330
1331 return res;
1332
1333 } /* end mp_expt() */
1334
1335 /* }}} */
1336
1337 /* {{{ mp_2expt(a, k) */
1338
1339 /* Compute a = 2^k */
1340
1341 mp_err mp_2expt(mp_int *a, mp_digit k)
1342 {
1343 ARGCHK(a != NULL, MP_BADARG);
1344
1345 return s_mp_2expt(a, k);
1346
1347 } /* end mp_2expt() */
1348
1349 /* }}} */
1350
1351 /* {{{ mp_mod(a, m, c) */
1352
1353 /*
1354 mp_mod(a, m, c)
1355
1356 Compute c = a (mod m). Result will always be 0 <= c < m.
1357 */
1358
1359 mp_err mp_mod(mp_int *a, mp_int *m, mp_int *c)
1360 {
1361 mp_err res;
1362 int mag;
1363
1364 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
1365
1366 if(SIGN(m) == MP_NEG)
1367 return MP_RANGE;
1368
1369 /*
1370 If |a| > m, we need to divide to get the remainder and take the
1371 absolute value.
1372
1373 If |a| < m, we don't need to do any division, just copy and adjust
1374 the sign (if a is negative).
1375
1376 If |a| == m, we can simply set the result to zero.
1377
1378 This order is intended to minimize the average path length of the
1379 comparison chain on common workloads -- the most frequent cases are
1380 that |a| != m, so we do those first.
1381 */
1382 if((mag = s_mp_cmp(a, m)) > 0) {
1383 if((res = mp_div(a, m, NULL, c)) != MP_OKAY)
1384 return res;
1385
1386 if(SIGN(c) == MP_NEG) {
1387 if((res = mp_add(c, m, c)) != MP_OKAY)
1388 return res;
1389 }
1390
1391 } else if(mag < 0) {
1392 if((res = mp_copy(a, c)) != MP_OKAY)
1393 return res;
1394
1395 if(mp_cmp_z(a) < 0) {
1396 if((res = mp_add(c, m, c)) != MP_OKAY)
1397 return res;
1398
1399 }
1400
1401 } else {
1402 mp_zero(c);
1403
1404 }
1405
1406 return MP_OKAY;
1407
1408 } /* end mp_mod() */
1409
1410 /* }}} */
1411
1412 /* {{{ mp_mod_d(a, d, c) */
1413
1414 /*
1415 mp_mod_d(a, d, c)
1416
1417 Compute c = a (mod d). Result will always be 0 <= c < d
1418 */
1419 mp_err mp_mod_d(mp_int *a, mp_digit d, mp_digit *c)
1420 {
1421 mp_err res;
1422 mp_digit rem;
1423
1424 ARGCHK(a != NULL && c != NULL, MP_BADARG);
1425
1426 if(s_mp_cmp_d(a, d) > 0) {
1427 if((res = mp_div_d(a, d, NULL, &rem)) != MP_OKAY)
1428 return res;
1429
1430 } else {
1431 if(SIGN(a) == MP_NEG)
1432 rem = d - DIGIT(a, 0);
1433 else
1434 rem = DIGIT(a, 0);
1435 }
1436
1437 if(c)
1438 *c = rem;
1439
1440 return MP_OKAY;
1441
1442 } /* end mp_mod_d() */
1443
1444 /* }}} */
1445
1446 /* {{{ mp_sqrt(a, b) */
1447
1448 /*
1449 mp_sqrt(a, b)
1450
1451 Compute the integer square root of a, and store the result in b.
1452 Uses an integer-arithmetic version of Newton's iterative linear
1453 approximation technique to determine this value; the result has the
1454 following two properties:
1455
1456 b^2 <= a
1457 (b+1)^2 >= a
1458
1459 It is a range error to pass a negative value.
1460 */
1461 mp_err mp_sqrt(mp_int *a, mp_int *b)
1462 {
1463 mp_int x, t;
1464 mp_err res;
1465
1466 ARGCHK(a != NULL && b != NULL, MP_BADARG);
1467
1468 /* Cannot take square root of a negative value */
1469 if(SIGN(a) == MP_NEG)
1470 return MP_RANGE;
1471
1472 /* Special cases for zero and one, trivial */
1473 if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ)
1474 return mp_copy(a, b);
1475
1476 /* Initialize the temporaries we'll use below */
1477 if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
1478 return res;
1479
1480 /* Compute an initial guess for the iteration as a itself */
1481 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1482 goto X;
1483
1484 for(;;) {
1485 /* t = (x * x) - a */
1486 mp_copy(&x, &t); /* can't fail, t is big enough for original x */
1487 if((res = mp_sqr(&t, &t)) != MP_OKAY ||
1488 (res = mp_sub(&t, a, &t)) != MP_OKAY)
1489 goto CLEANUP;
1490
1491 /* t = t / 2x */
1492 s_mp_mul_2(&x);
1493 if((res = mp_div(&t, &x, &t, NULL)) != MP_OKAY)
1494 goto CLEANUP;
1495 s_mp_div_2(&x);
1496
1497 /* Terminate the loop, if the quotient is zero */
1498 if(mp_cmp_z(&t) == MP_EQ)
1499 break;
1500
1501 /* x = x - t */
1502 if((res = mp_sub(&x, &t, &x)) != MP_OKAY)
1503 goto CLEANUP;
1504
1505 }
1506
1507 /* Copy result to output parameter */
1508 mp_sub_d(&x, 1, &x);
1509 s_mp_exch(&x, b);
1510
1511 CLEANUP:
1512 mp_clear(&x);
1513 X:
1514 mp_clear(&t);
1515
1516 return res;
1517
1518 } /* end mp_sqrt() */
1519
1520 /* }}} */
1521
1522 /* }}} */
1523
1524 /*------------------------------------------------------------------------*/
1525 /* {{{ Modular arithmetic */
1526
1527 #if MP_MODARITH
1528 /* {{{ mp_addmod(a, b, m, c) */
1529
1530 /*
1531 mp_addmod(a, b, m, c)
1532
1533 Compute c = (a + b) mod m
1534 */
1535
1536 mp_err mp_addmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1537 {
1538 mp_err res;
1539
1540 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1541
1542 if((res = mp_add(a, b, c)) != MP_OKAY)
1543 return res;
1544 if((res = mp_mod(c, m, c)) != MP_OKAY)
1545 return res;
1546
1547 return MP_OKAY;
1548
1549 }
1550
1551 /* }}} */
1552
1553 /* {{{ mp_submod(a, b, m, c) */
1554
1555 /*
1556 mp_submod(a, b, m, c)
1557
1558 Compute c = (a - b) mod m
1559 */
1560
1561 mp_err mp_submod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1562 {
1563 mp_err res;
1564
1565 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1566
1567 if((res = mp_sub(a, b, c)) != MP_OKAY)
1568 return res;
1569 if((res = mp_mod(c, m, c)) != MP_OKAY)
1570 return res;
1571
1572 return MP_OKAY;
1573
1574 }
1575
1576 /* }}} */
1577
1578 /* {{{ mp_mulmod(a, b, m, c) */
1579
1580 /*
1581 mp_mulmod(a, b, m, c)
1582
1583 Compute c = (a * b) mod m
1584 */
1585
1586 mp_err mp_mulmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1587 {
1588 mp_err res;
1589
1590 ARGCHK(a != NULL && b != NULL && m != NULL && c != NULL, MP_BADARG);
1591
1592 if((res = mp_mul(a, b, c)) != MP_OKAY)
1593 return res;
1594 if((res = mp_mod(c, m, c)) != MP_OKAY)
1595 return res;
1596
1597 return MP_OKAY;
1598
1599 }
1600
1601 /* }}} */
1602
1603 /* {{{ mp_sqrmod(a, m, c) */
1604
1605 #if MP_SQUARE
1606 mp_err mp_sqrmod(mp_int *a, mp_int *m, mp_int *c)
1607 {
1608 mp_err res;
1609
1610 ARGCHK(a != NULL && m != NULL && c != NULL, MP_BADARG);
1611
1612 if((res = mp_sqr(a, c)) != MP_OKAY)
1613 return res;
1614 if((res = mp_mod(c, m, c)) != MP_OKAY)
1615 return res;
1616
1617 return MP_OKAY;
1618
1619 } /* end mp_sqrmod() */
1620 #endif
1621
1622 /* }}} */
1623
1624 /* {{{ mp_exptmod(a, b, m, c) */
1625
1626 /*
1627 mp_exptmod(a, b, m, c)
1628
1629 Compute c = (a ** b) mod m. Uses a standard square-and-multiply
1630 method with modular reductions at each step. (This is basically the
1631 same code as mp_expt(), except for the addition of the reductions)
1632
1633 The modular reductions are done using Barrett's algorithm (see
1634 s_mp_reduce() below for details)
1635 */
1636
1637 mp_err mp_exptmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c)
1638 {
1639 mp_int s, x, mu;
1640 mp_err res;
1641 mp_digit d, *db = DIGITS(b);
1642 mp_size ub = USED(b);
1643 int dig, bit;
1644
1645 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1646
1647 if(mp_cmp_z(b) < 0 || mp_cmp_z(m) <= 0)
1648 return MP_RANGE;
1649
1650 if((res = mp_init(&s)) != MP_OKAY)
1651 return res;
1652 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1653 goto X;
1654 if((res = mp_mod(&x, m, &x)) != MP_OKAY ||
1655 (res = mp_init(&mu)) != MP_OKAY)
1656 goto MU;
1657
1658 mp_set(&s, 1);
1659
1660 /* mu = b^2k / m */
1661 s_mp_add_d(&mu, 1);
1662 s_mp_lshd(&mu, 2 * USED(m));
1663 if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY)
1664 goto CLEANUP;
1665
1666 /* Loop over digits of b in ascending order, except highest order */
1667 for(dig = 0; dig < (ub - 1); dig++) {
1668 d = *db++;
1669
1670 /* Loop over the bits of the lower-order digits */
1671 for(bit = 0; bit < DIGIT_BIT; bit++) {
1672 if(d & 1) {
1673 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1674 goto CLEANUP;
1675 if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)
1676 goto CLEANUP;
1677 }
1678
1679 d >>= 1;
1680
1681 if((res = s_mp_sqr(&x)) != MP_OKAY)
1682 goto CLEANUP;
1683 if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)
1684 goto CLEANUP;
1685 }
1686 }
1687
1688 /* Now do the last digit... */
1689 d = *db;
1690
1691 while(d) {
1692 if(d & 1) {
1693 if((res = s_mp_mul(&s, &x)) != MP_OKAY)
1694 goto CLEANUP;
1695 if((res = s_mp_reduce(&s, m, &mu)) != MP_OKAY)
1696 goto CLEANUP;
1697 }
1698
1699 d >>= 1;
1700
1701 if((res = s_mp_sqr(&x)) != MP_OKAY)
1702 goto CLEANUP;
1703 if((res = s_mp_reduce(&x, m, &mu)) != MP_OKAY)
1704 goto CLEANUP;
1705 }
1706
1707 s_mp_exch(&s, c);
1708
1709 CLEANUP:
1710 mp_clear(&mu);
1711 MU:
1712 mp_clear(&x);
1713 X:
1714 mp_clear(&s);
1715
1716 return res;
1717
1718 } /* end mp_exptmod() */
1719
1720 /* }}} */
1721
1722 /* {{{ mp_exptmod_d(a, d, m, c) */
1723
1724 mp_err mp_exptmod_d(mp_int *a, mp_digit d, mp_int *m, mp_int *c)
1725 {
1726 mp_int s, x;
1727 mp_err res;
1728
1729 ARGCHK(a != NULL && c != NULL, MP_BADARG);
1730
1731 if((res = mp_init(&s)) != MP_OKAY)
1732 return res;
1733 if((res = mp_init_copy(&x, a)) != MP_OKAY)
1734 goto X;
1735
1736 mp_set(&s, 1);
1737
1738 while(d != 0) {
1739 if(d & 1) {
1740 if((res = s_mp_mul(&s, &x)) != MP_OKAY ||
1741 (res = mp_mod(&s, m, &s)) != MP_OKAY)
1742 goto CLEANUP;
1743 }
1744
1745 d /= 2;
1746
1747 if((res = s_mp_sqr(&x)) != MP_OKAY ||
1748 (res = mp_mod(&x, m, &x)) != MP_OKAY)
1749 goto CLEANUP;
1750 }
1751
1752 s_mp_exch(&s, c);
1753
1754 CLEANUP:
1755 mp_clear(&x);
1756 X:
1757 mp_clear(&s);
1758
1759 return res;
1760
1761 } /* end mp_exptmod_d() */
1762
1763 /* }}} */
1764 #endif /* if MP_MODARITH */
1765
1766 /* }}} */
1767
1768 /*------------------------------------------------------------------------*/
1769 /* {{{ Comparison functions */
1770
1771 /* {{{ mp_cmp_z(a) */
1772
1773 /*
1774 mp_cmp_z(a)
1775
1776 Compare a <=> 0. Returns <0 if a<0, 0 if a=0, >0 if a>0.
1777 */
1778
1779 int mp_cmp_z(mp_int *a)
1780 {
1781 if(SIGN(a) == MP_NEG)
1782 return MP_LT;
1783 else if(USED(a) == 1 && DIGIT(a, 0) == 0)
1784 return MP_EQ;
1785 else
1786 return MP_GT;
1787
1788 } /* end mp_cmp_z() */
1789
1790 /* }}} */
1791
1792 /* {{{ mp_cmp_d(a, d) */
1793
1794 /*
1795 mp_cmp_d(a, d)
1796
1797 Compare a <=> d. Returns <0 if a<d, 0 if a=d, >0 if a>d
1798 */
1799
1800 int mp_cmp_d(mp_int *a, mp_digit d)
1801 {
1802 ARGCHK(a != NULL, MP_EQ);
1803
1804 if(SIGN(a) == MP_NEG)
1805 return MP_LT;
1806
1807 return s_mp_cmp_d(a, d);
1808
1809 } /* end mp_cmp_d() */
1810
1811 /* }}} */
1812
1813 /* {{{ mp_cmp(a, b) */
1814
1815 int mp_cmp(mp_int *a, mp_int *b)
1816 {
1817 ARGCHK(a != NULL && b != NULL, MP_EQ);
1818
1819 if(SIGN(a) == SIGN(b)) {
1820 int mag;
1821
1822 if((mag = s_mp_cmp(a, b)) == MP_EQ)
1823 return MP_EQ;
1824
1825 if(SIGN(a) == MP_ZPOS)
1826 return mag;
1827 else
1828 return -mag;
1829
1830 } else if(SIGN(a) == MP_ZPOS) {
1831 return MP_GT;
1832 } else {
1833 return MP_LT;
1834 }
1835
1836 } /* end mp_cmp() */
1837
1838 /* }}} */
1839
1840 /* {{{ mp_cmp_mag(a, b) */
1841
1842 /*
1843 mp_cmp_mag(a, b)
1844
1845 Compares |a| <=> |b|, and returns an appropriate comparison result
1846 */
1847
1848 int mp_cmp_mag(mp_int *a, mp_int *b)
1849 {
1850 ARGCHK(a != NULL && b != NULL, MP_EQ);
1851
1852 return s_mp_cmp(a, b);
1853
1854 } /* end mp_cmp_mag() */
1855
1856 /* }}} */
1857
1858 /* {{{ mp_cmp_int(a, z) */
1859
1860 /*
1861 This just converts z to an mp_int, and uses the existing comparison
1862 routines. This is sort of inefficient, but it's not clear to me how
1863 frequently this wil get used anyway. For small positive constants,
1864 you can always use mp_cmp_d(), and for zero, there is mp_cmp_z().
1865 */
1866 int mp_cmp_int(mp_int *a, long z)
1867 {
1868 mp_int tmp;
1869 int out;
1870
1871 ARGCHK(a != NULL, MP_EQ);
1872
1873 mp_init(&tmp); mp_set_int(&tmp, z);
1874 out = mp_cmp(a, &tmp);
1875 mp_clear(&tmp);
1876
1877 return out;
1878
1879 } /* end mp_cmp_int() */
1880
1881 /* }}} */
1882
1883 /* {{{ mp_isodd(a) */
1884
1885 /*
1886 mp_isodd(a)
1887
1888 Returns a true (non-zero) value if a is odd, false (zero) otherwise.
1889 */
1890 int mp_isodd(mp_int *a)
1891 {
1892 ARGCHK(a != NULL, 0);
1893
1894 return (DIGIT(a, 0) & 1);
1895
1896 } /* end mp_isodd() */
1897
1898 /* }}} */
1899
1900 /* {{{ mp_iseven(a) */
1901
1902 int mp_iseven(mp_int *a)
1903 {
1904 return !mp_isodd(a);
1905
1906 } /* end mp_iseven() */
1907
1908 /* }}} */
1909
1910 /* }}} */
1911
1912 /*------------------------------------------------------------------------*/
1913 /* {{{ Number theoretic functions */
1914
1915 #if MP_NUMTH
1916 /* {{{ mp_gcd(a, b, c) */
1917
1918 /*
1919 Like the old mp_gcd() function, except computes the GCD using the
1920 binary algorithm due to Josef Stein in 1961 (via Knuth).
1921 */
1922 mp_err mp_gcd(mp_int *a, mp_int *b, mp_int *c)
1923 {
1924 mp_err res;
1925 mp_int u, v, t;
1926 mp_size k = 0;
1927
1928 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
1929
1930 if(mp_cmp_z(a) == MP_EQ && mp_cmp_z(b) == MP_EQ)
1931 return MP_RANGE;
1932 if(mp_cmp_z(a) == MP_EQ) {
1933 if((res = mp_copy(b, c)) != MP_OKAY)
1934 return res;
1935 SIGN(c) = MP_ZPOS; return MP_OKAY;
1936 } else if(mp_cmp_z(b) == MP_EQ) {
1937 if((res = mp_copy(a, c)) != MP_OKAY)
1938 return res;
1939 SIGN(c) = MP_ZPOS; return MP_OKAY;
1940 }
1941
1942 if((res = mp_init(&t)) != MP_OKAY)
1943 return res;
1944 if((res = mp_init_copy(&u, a)) != MP_OKAY)
1945 goto U;
1946 if((res = mp_init_copy(&v, b)) != MP_OKAY)
1947 goto V;
1948
1949 SIGN(&u) = MP_ZPOS;
1950 SIGN(&v) = MP_ZPOS;
1951
1952 /* Divide out common factors of 2 until at least 1 of a, b is even */
1953 while(mp_iseven(&u) && mp_iseven(&v)) {
1954 s_mp_div_2(&u);
1955 s_mp_div_2(&v);
1956 ++k;
1957 }
1958
1959 /* Initialize t */
1960 if(mp_isodd(&u)) {
1961 if((res = mp_copy(&v, &t)) != MP_OKAY)
1962 goto CLEANUP;
1963
1964 /* t = -v */
1965 if(SIGN(&v) == MP_ZPOS)
1966 SIGN(&t) = MP_NEG;
1967 else
1968 SIGN(&t) = MP_ZPOS;
1969
1970 } else {
1971 if((res = mp_copy(&u, &t)) != MP_OKAY)
1972 goto CLEANUP;
1973
1974 }
1975
1976 for(;;) {
1977 while(mp_iseven(&t)) {
1978 s_mp_div_2(&t);
1979 }
1980
1981 if(mp_cmp_z(&t) == MP_GT) {
1982 if((res = mp_copy(&t, &u)) != MP_OKAY)
1983 goto CLEANUP;
1984
1985 } else {
1986 if((res = mp_copy(&t, &v)) != MP_OKAY)
1987 goto CLEANUP;
1988
1989 /* v = -t */
1990 if(SIGN(&t) == MP_ZPOS)
1991 SIGN(&v) = MP_NEG;
1992 else
1993 SIGN(&v) = MP_ZPOS;
1994 }
1995
1996 if((res = mp_sub(&u, &v, &t)) != MP_OKAY)
1997 goto CLEANUP;
1998
1999 if(s_mp_cmp_d(&t, 0) == MP_EQ)
2000 break;
2001 }
2002
2003 s_mp_2expt(&v, k); /* v = 2^k */
2004 res = mp_mul(&u, &v, c); /* c = u * v */
2005
2006 CLEANUP:
2007 mp_clear(&v);
2008 V:
2009 mp_clear(&u);
2010 U:
2011 mp_clear(&t);
2012
2013 return res;
2014
2015 } /* end mp_bgcd() */
2016
2017 /* }}} */
2018
2019 /* {{{ mp_lcm(a, b, c) */
2020
2021 /* We compute the least common multiple using the rule:
2022
2023 ab = [a, b](a, b)
2024
2025 ... by computing the product, and dividing out the gcd.
2026 */
2027
2028 mp_err mp_lcm(mp_int *a, mp_int *b, mp_int *c)
2029 {
2030 mp_int gcd, prod;
2031 mp_err res;
2032
2033 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG);
2034
2035 /* Set up temporaries */
2036 if((res = mp_init(&gcd)) != MP_OKAY)
2037 return res;
2038 if((res = mp_init(&prod)) != MP_OKAY)
2039 goto GCD;
2040
2041 if((res = mp_mul(a, b, &prod)) != MP_OKAY)
2042 goto CLEANUP;
2043 if((res = mp_gcd(a, b, &gcd)) != MP_OKAY)
2044 goto CLEANUP;
2045
2046 res = mp_div(&prod, &gcd, c, NULL);
2047
2048 CLEANUP:
2049 mp_clear(&prod);
2050 GCD:
2051 mp_clear(&gcd);
2052
2053 return res;
2054
2055 } /* end mp_lcm() */
2056
2057 /* }}} */
2058
2059 /* {{{ mp_xgcd(a, b, g, x, y) */
2060
2061 /*
2062 mp_xgcd(a, b, g, x, y)
2063
2064 Compute g = (a, b) and values x and y satisfying Bezout's identity
2065 (that is, ax + by = g). This uses the extended binary GCD algorithm
2066 based on the Stein algorithm used for mp_gcd()
2067 */
2068
2069 mp_err mp_xgcd(mp_int *a, mp_int *b, mp_int *g, mp_int *x, mp_int *y)
2070 {
2071 mp_int gx, xc, yc, u, v, A, B, C, D;
2072 mp_int *clean[9];
2073 mp_err res;
2074 int last = -1;
2075
2076 if(mp_cmp_z(b) == 0)
2077 return MP_RANGE;
2078
2079 /* Initialize all these variables we need */
2080 if((res = mp_init(&u)) != MP_OKAY) goto CLEANUP;
2081 clean[++last] = &u;
2082 if((res = mp_init(&v)) != MP_OKAY) goto CLEANUP;
2083 clean[++last] = &v;
2084 if((res = mp_init(&gx)) != MP_OKAY) goto CLEANUP;
2085 clean[++last] = &gx;
2086 if((res = mp_init(&A)) != MP_OKAY) goto CLEANUP;
2087 clean[++last] = &A;
2088 if((res = mp_init(&B)) != MP_OKAY) goto CLEANUP;
2089 clean[++last] = &B;
2090 if((res = mp_init(&C)) != MP_OKAY) goto CLEANUP;
2091 clean[++last] = &C;
2092 if((res = mp_init(&D)) != MP_OKAY) goto CLEANUP;
2093 clean[++last] = &D;
2094 if((res = mp_init_copy(&xc, a)) != MP_OKAY) goto CLEANUP;
2095 clean[++last] = &xc;
2096 mp_abs(&xc, &xc);
2097 if((res = mp_init_copy(&yc, b)) != MP_OKAY) goto CLEANUP;
2098 clean[++last] = &yc;
2099 mp_abs(&yc, &yc);
2100
2101 mp_set(&gx, 1);
2102
2103 /* Divide by two until at least one of them is even */
2104 while(mp_iseven(&xc) && mp_iseven(&yc)) {
2105 s_mp_div_2(&xc);
2106 s_mp_div_2(&yc);
2107 if((res = s_mp_mul_2(&gx)) != MP_OKAY)
2108 goto CLEANUP;
2109 }
2110
2111 mp_copy(&xc, &u);
2112 mp_copy(&yc, &v);
2113 mp_set(&A, 1); mp_set(&D, 1);
2114
2115 /* Loop through binary GCD algorithm */
2116 for(;;) {
2117 while(mp_iseven(&u)) {
2118 s_mp_div_2(&u);
2119
2120 if(mp_iseven(&A) && mp_iseven(&B)) {
2121 s_mp_div_2(&A); s_mp_div_2(&B);
2122 } else {
2123 if((res = mp_add(&A, &yc, &A)) != MP_OKAY) goto CLEANUP;
2124 s_mp_div_2(&A);
2125 if((res = mp_sub(&B, &xc, &B)) != MP_OKAY) goto CLEANUP;
2126 s_mp_div_2(&B);
2127 }
2128 }
2129
2130 while(mp_iseven(&v)) {
2131 s_mp_div_2(&v);
2132
2133 if(mp_iseven(&C) && mp_iseven(&D)) {
2134 s_mp_div_2(&C); s_mp_div_2(&D);
2135 } else {
2136 if((res = mp_add(&C, &yc, &C)) != MP_OKAY) goto CLEANUP;
2137 s_mp_div_2(&C);
2138 if((res = mp_sub(&D, &xc, &D)) != MP_OKAY) goto CLEANUP;
2139 s_mp_div_2(&D);
2140 }
2141 }
2142
2143 if(mp_cmp(&u, &v) >= 0) {
2144 if((res = mp_sub(&u, &v, &u)) != MP_OKAY) goto CLEANUP;
2145 if((res = mp_sub(&A, &C, &A)) != MP_OKAY) goto CLEANUP;
2146 if((res = mp_sub(&B, &D, &B)) != MP_OKAY) goto CLEANUP;
2147
2148 } else {
2149 if((res = mp_sub(&v, &u, &v)) != MP_OKAY) goto CLEANUP;
2150 if((res = mp_sub(&C, &A, &C)) != MP_OKAY) goto CLEANUP;
2151 if((res = mp_sub(&D, &B, &D)) != MP_OKAY) goto CLEANUP;
2152
2153 }
2154
2155 /* If we're done, copy results to output */
2156 if(mp_cmp_z(&u) == 0) {
2157 if(x)
2158 if((res = mp_copy(&C, x)) != MP_OKAY) goto CLEANUP;
2159
2160 if(y)
2161 if((res = mp_copy(&D, y)) != MP_OKAY) goto CLEANUP;
2162
2163 if(g)
2164 if((res = mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP;
2165
2166 break;
2167 }
2168 }
2169
2170 CLEANUP:
2171 while(last >= 0)
2172 mp_clear(clean[last--]);
2173
2174 return res;
2175
2176 } /* end mp_xgcd() */
2177
2178 /* }}} */
2179
2180 /* {{{ mp_invmod(a, m, c) */
2181
2182 /*
2183 mp_invmod(a, m, c)
2184
2185 Compute c = a^-1 (mod m), if there is an inverse for a (mod m).
2186 This is equivalent to the question of whether (a, m) = 1. If not,
2187 MP_UNDEF is returned, and there is no inverse.
2188 */
2189
2190 mp_err mp_invmod(mp_int *a, mp_int *m, mp_int *c)
2191 {
2192 mp_int g, x;
2193 mp_sign sa;
2194 mp_err res;
2195
2196 ARGCHK(a && m && c, MP_BADARG);
2197
2198 if(mp_cmp_z(a) == 0 || mp_cmp_z(m) == 0)
2199 return MP_RANGE;
2200
2201 sa = SIGN(a);
2202
2203 if((res = mp_init(&g)) != MP_OKAY)
2204 return res;
2205 if((res = mp_init(&x)) != MP_OKAY)
2206 goto X;
2207
2208 if((res = mp_xgcd(a, m, &g, &x, NULL)) != MP_OKAY)
2209 goto CLEANUP;
2210
2211 if(mp_cmp_d(&g, 1) != MP_EQ) {
2212 res = MP_UNDEF;
2213 goto CLEANUP;
2214 }
2215
2216 res = mp_mod(&x, m, c);
2217 SIGN(c) = sa;
2218
2219 CLEANUP:
2220 mp_clear(&x);
2221 X:
2222 mp_clear(&g);
2223
2224 return res;
2225
2226 } /* end mp_invmod() */
2227
2228 /* }}} */
2229 #endif /* if MP_NUMTH */
2230
2231 /* }}} */
2232
2233 /*------------------------------------------------------------------------*/
2234 /* {{{ mp_print(mp, ofp) */
2235
2236 #if MP_IOFUNC
2237 /*
2238 mp_print(mp, ofp)
2239
2240 Print a textual representation of the given mp_int on the output
2241 stream 'ofp'. Output is generated using the internal radix.
2242 */
2243
2244 void mp_print(mp_int *mp, FILE *ofp)
2245 {
2246 int ix;
2247
2248 if(mp == NULL || ofp == NULL)
2249 return;
2250
2251 fputc((SIGN(mp) == MP_NEG) ? '-' : '+', ofp);
2252
2253 for(ix = USED(mp) - 1; ix >= 0; ix--) {
2254 fprintf(ofp, DIGIT_FMT, DIGIT(mp, ix));
2255 }
2256
2257 } /* end mp_print() */
2258
2259 #endif /* if MP_IOFUNC */
2260
2261 /* }}} */
2262
2263 /*------------------------------------------------------------------------*/
2264 /* {{{ More I/O Functions */
2265
2266 /* {{{ mp_read_signed_bin(mp, str, len) */
2267
2268 /*
2269 mp_read_signed_bin(mp, str, len)
2270
2271 Read in a raw value (base 256) into the given mp_int
2272 */
2273
2274 mp_err mp_read_signed_bin(mp_int *mp, unsigned char *str, int len)
2275 {
2276 mp_err res;
2277
2278 ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
2279
2280 if((res = mp_read_unsigned_bin(mp, str + 1, len - 1)) == MP_OKAY) {
2281 /* Get sign from first byte */
2282 if(str[0])
2283 SIGN(mp) = MP_NEG;
2284 else
2285 SIGN(mp) = MP_ZPOS;
2286 }
2287
2288 return res;
2289
2290 } /* end mp_read_signed_bin() */
2291
2292 /* }}} */
2293
2294 /* {{{ mp_signed_bin_size(mp) */
2295
2296 int mp_signed_bin_size(mp_int *mp)
2297 {
2298 ARGCHK(mp != NULL, 0);
2299
2300 return mp_unsigned_bin_size(mp) + 1;
2301
2302 } /* end mp_signed_bin_size() */
2303
2304 /* }}} */
2305
2306 /* {{{ mp_to_signed_bin(mp, str) */
2307
2308 mp_err mp_to_signed_bin(mp_int *mp, unsigned char *str)
2309 {
2310 ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2311
2312 /* Caller responsible for allocating enough memory (use mp_raw_size(mp)) */
2313 str[0] = (char)SIGN(mp);
2314
2315 return mp_to_unsigned_bin(mp, str + 1);
2316
2317 } /* end mp_to_signed_bin() */
2318
2319 /* }}} */
2320
2321 /* {{{ mp_read_unsigned_bin(mp, str, len) */
2322
2323 /*
2324 mp_read_unsigned_bin(mp, str, len)
2325
2326 Read in an unsigned value (base 256) into the given mp_int
2327 */
2328
2329 mp_err mp_read_unsigned_bin(mp_int *mp, unsigned char *str, int len)
2330 {
2331 int ix;
2332 mp_err res;
2333
2334 ARGCHK(mp != NULL && str != NULL && len > 0, MP_BADARG);
2335
2336 mp_zero(mp);
2337
2338 for(ix = 0; ix < len; ix++) {
2339 if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY)
2340 return res;
2341
2342 if((res = mp_add_d(mp, str[ix], mp)) != MP_OKAY)
2343 return res;
2344 }
2345
2346 return MP_OKAY;
2347
2348 } /* end mp_read_unsigned_bin() */
2349
2350 /* }}} */
2351
2352 /* {{{ mp_unsigned_bin_size(mp) */
2353
2354 int mp_unsigned_bin_size(mp_int *mp)
2355 {
2356 mp_digit topdig;
2357 int count;
2358
2359 ARGCHK(mp != NULL, 0);
2360
2361 /* Special case for the value zero */
2362 if(USED(mp) == 1 && DIGIT(mp, 0) == 0)
2363 return 1;
2364
2365 count = (USED(mp) - 1) * sizeof(mp_digit);
2366 topdig = DIGIT(mp, USED(mp) - 1);
2367
2368 while(topdig != 0) {
2369 ++count;
2370 topdig >>= CHAR_BIT;
2371 }
2372
2373 return count;
2374
2375 } /* end mp_unsigned_bin_size() */
2376
2377 /* }}} */
2378
2379 /* {{{ mp_to_unsigned_bin(mp, str) */
2380
2381 mp_err mp_to_unsigned_bin(mp_int *mp, unsigned char *str)
2382 {
2383 mp_digit *dp, *end, d;
2384 unsigned char *spos;
2385
2386 ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2387
2388 dp = DIGITS(mp);
2389 end = dp + USED(mp) - 1;
2390 spos = str;
2391
2392 /* Special case for zero, quick test */
2393 if(dp == end && *dp == 0) {
2394 *str = '\0';
2395 return MP_OKAY;
2396 }
2397
2398 /* Generate digits in reverse order */
2399 while(dp < end) {
2400 int ix;
2401
2402 d = *dp;
2403 for(ix = 0; ix < sizeof(mp_digit); ++ix) {
2404 *spos = d & UCHAR_MAX;
2405 d >>= CHAR_BIT;
2406 ++spos;
2407 }
2408
2409 ++dp;
2410 }
2411
2412 /* Now handle last digit specially, high order zeroes are not written */
2413 d = *end;
2414 while(d != 0) {
2415 *spos = d & UCHAR_MAX;
2416 d >>= CHAR_BIT;
2417 ++spos;
2418 }
2419
2420 /* Reverse everything to get digits in the correct order */
2421 while(--spos > str) {
2422 unsigned char t = *str;
2423 *str = *spos;
2424 *spos = t;
2425
2426 ++str;
2427 }
2428
2429 return MP_OKAY;
2430
2431 } /* end mp_to_unsigned_bin() */
2432
2433 /* }}} */
2434
2435 /* {{{ mp_count_bits(mp) */
2436
2437 int mp_count_bits(mp_int *mp)
2438 {
2439 int len;
2440 mp_digit d;
2441
2442 ARGCHK(mp != NULL, MP_BADARG);
2443
2444 len = DIGIT_BIT * (USED(mp) - 1);
2445 d = DIGIT(mp, USED(mp) - 1);
2446
2447 while(d != 0) {
2448 ++len;
2449 d >>= 1;
2450 }
2451
2452 return len;
2453
2454 } /* end mp_count_bits() */
2455
2456 /* }}} */
2457
2458 /* {{{ mp_read_radix(mp, str, radix) */
2459
2460 /*
2461 mp_read_radix(mp, str, radix)
2462
2463 Read an integer from the given string, and set mp to the resulting
2464 value. The input is presumed to be in base 10. Leading non-digit
2465 characters are ignored, and the function reads until a non-digit
2466 character or the end of the string.
2467 */
2468
2469 mp_err mp_read_radix(mp_int *mp, unsigned char *str, int radix)
2470 {
2471 int ix = 0, val = 0;
2472 mp_err res;
2473 mp_sign sig = MP_ZPOS;
2474
2475 ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX,
2476 MP_BADARG);
2477
2478 mp_zero(mp);
2479
2480 /* Skip leading non-digit characters until a digit or '-' or '+' */
2481 while(str[ix] &&
2482 (s_mp_tovalue(str[ix], radix) < 0) &&
2483 str[ix] != '-' &&
2484 str[ix] != '+') {
2485 ++ix;
2486 }
2487
2488 if(str[ix] == '-') {
2489 sig = MP_NEG;
2490 ++ix;
2491 } else if(str[ix] == '+') {
2492 sig = MP_ZPOS; /* this is the default anyway... */
2493 ++ix;
2494 }
2495
2496 while((val = s_mp_tovalue(str[ix], radix)) >= 0) {
2497 if((res = s_mp_mul_d(mp, radix)) != MP_OKAY)
2498 return res;
2499 if((res = s_mp_add_d(mp, val)) != MP_OKAY)
2500 return res;
2501 ++ix;
2502 }
2503
2504 if(s_mp_cmp_d(mp, 0) == MP_EQ)
2505 SIGN(mp) = MP_ZPOS;
2506 else
2507 SIGN(mp) = sig;
2508
2509 return MP_OKAY;
2510
2511 } /* end mp_read_radix() */
2512
2513 /* }}} */
2514
2515 /* {{{ mp_radix_size(mp, radix) */
2516
2517 int mp_radix_size(mp_int *mp, int radix)
2518 {
2519 int len;
2520 ARGCHK(mp != NULL, 0);
2521
2522 len = s_mp_outlen(mp_count_bits(mp), radix) + 1; /* for NUL terminator */
2523
2524 if(mp_cmp_z(mp) < 0)
2525 ++len; /* for sign */
2526
2527 return len;
2528
2529 } /* end mp_radix_size() */
2530
2531 /* }}} */
2532
2533 /* {{{ mp_value_radix_size(num, qty, radix) */
2534
2535 /* num = number of digits
2536 qty = number of bits per digit
2537 radix = target base
2538
2539 Return the number of digits in the specified radix that would be
2540 needed to express 'num' digits of 'qty' bits each.
2541 */
2542 int mp_value_radix_size(int num, int qty, int radix)
2543 {
2544 ARGCHK(num >= 0 && qty > 0 && radix >= 2 && radix <= MAX_RADIX, 0);
2545
2546 return s_mp_outlen(num * qty, radix);
2547
2548 } /* end mp_value_radix_size() */
2549
2550 /* }}} */
2551
2552 /* {{{ mp_toradix(mp, str, radix) */
2553
2554 mp_err mp_toradix(mp_int *mp, unsigned char *str, int radix)
2555 {
2556 int ix, pos = 0;
2557
2558 ARGCHK(mp != NULL && str != NULL, MP_BADARG);
2559 ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE);
2560
2561 if(mp_cmp_z(mp) == MP_EQ) {
2562 str[0] = '0';
2563 str[1] = '\0';
2564 } else {
2565 mp_err res;
2566 mp_int tmp;
2567 mp_sign sgn;
2568 mp_digit rem, rdx = (mp_digit)radix;
2569 char ch;
2570
2571 if((res = mp_init_copy(&tmp, mp)) != MP_OKAY)
2572 return res;
2573
2574 /* Save sign for later, and take absolute value */
2575 sgn = SIGN(&tmp); SIGN(&tmp) = MP_ZPOS;
2576
2577 /* Generate output digits in reverse order */
2578 while(mp_cmp_z(&tmp) != 0) {
2579 if((res = s_mp_div_d(&tmp, rdx, &rem)) != MP_OKAY) {
2580 mp_clear(&tmp);
2581 return res;
2582 }
2583
2584 /* Generate digits, use capital letters */
2585 ch = s_mp_todigit(rem, radix, 0);
2586
2587 str[pos++] = ch;
2588 }
2589
2590 /* Add - sign if original value was negative */
2591 if(sgn == MP_NEG)
2592 str[pos++] = '-';
2593
2594 /* Add trailing NUL to end the string */
2595 str[pos--] = '\0';
2596
2597 /* Reverse the digits and sign indicator */
2598 ix = 0;
2599 while(ix < pos) {
2600 char tmp = str[ix];
2601
2602 str[ix] = str[pos];
2603 str[pos] = tmp;
2604 ++ix;
2605 --pos;
2606 }
2607
2608 mp_clear(&tmp);
2609 }
2610
2611 return MP_OKAY;
2612
2613 } /* end mp_toradix() */
2614
2615 /* }}} */
2616
2617 /* {{{ mp_char2value(ch, r) */
2618
2619 int mp_char2value(char ch, int r)
2620 {
2621 return s_mp_tovalue(ch, r);
2622
2623 } /* end mp_tovalue() */
2624
2625 /* }}} */
2626
2627 /* }}} */
2628
2629 /* {{{ mp_strerror(ec) */
2630
2631 /*
2632 mp_strerror(ec)
2633
2634 Return a string describing the meaning of error code 'ec'. The
2635 string returned is allocated in static memory, so the caller should
2636 not attempt to modify or free the memory associated with this
2637 string.
2638 */
2639 const char *mp_strerror(mp_err ec)
2640 {
2641 int aec = (ec < 0) ? -ec : ec;
2642
2643 /* Code values are negative, so the senses of these comparisons
2644 are accurate */
2645 if(ec < MP_LAST_CODE || ec > MP_OKAY) {
2646 return mp_err_string[0]; /* unknown error code */
2647 } else {
2648 return mp_err_string[aec + 1];
2649 }
2650
2651 } /* end mp_strerror() */
2652
2653 /* }}} */
2654
2655 /*========================================================================*/
2656 /*------------------------------------------------------------------------*/
2657 /* Static function definitions (internal use only) */
2658
2659 /* {{{ Memory management */
2660
2661 /* {{{ s_mp_grow(mp, min) */
2662
2663 /* Make sure there are at least 'min' digits allocated to mp */
2664 mp_err s_mp_grow(mp_int *mp, mp_size min)
2665 {
2666 if(min > ALLOC(mp)) {
2667 mp_digit *tmp;
2668
2669 /* Set min to next nearest default precision block size */
2670 min = ((min + (s_mp_defprec - 1)) / s_mp_defprec) * s_mp_defprec;
2671
2672 if((tmp = s_mp_alloc(min, sizeof(mp_digit))) == NULL)
2673 return MP_MEM;
2674
2675 s_mp_copy(DIGITS(mp), tmp, USED(mp));
2676
2677 #if MP_CRYPTO
2678 s_mp_setz(DIGITS(mp), ALLOC(mp));
2679 #endif
2680 s_mp_free(DIGITS(mp));
2681 DIGITS(mp) = tmp;
2682 ALLOC(mp) = min;
2683 }
2684
2685 return MP_OKAY;
2686
2687 } /* end s_mp_grow() */
2688
2689 /* }}} */
2690
2691 /* {{{ s_mp_pad(mp, min) */
2692
2693 /* Make sure the used size of mp is at least 'min', growing if needed */
2694 mp_err s_mp_pad(mp_int *mp, mp_size min)
2695 {
2696 if(min > USED(mp)) {
2697 mp_err res;
2698
2699 /* Make sure there is room to increase precision */
2700 if(min > ALLOC(mp) && (res = s_mp_grow(mp, min)) != MP_OKAY)
2701 return res;
2702
2703 /* Increase precision; should already be 0-filled */
2704 USED(mp) = min;
2705 }
2706
2707 return MP_OKAY;
2708
2709 } /* end s_mp_pad() */
2710
2711 /* }}} */
2712
2713 /* {{{ s_mp_setz(dp, count) */
2714
2715 #if MP_MACRO == 0
2716 /* Set 'count' digits pointed to by dp to be zeroes */
2717 void s_mp_setz(mp_digit *dp, mp_size count)
2718 {
2719 #if MP_MEMSET == 0
2720 int ix;
2721
2722 for(ix = 0; ix < count; ix++)
2723 dp[ix] = 0;
2724 #else
2725 memset(dp, 0, count * sizeof(mp_digit));
2726 #endif
2727
2728 } /* end s_mp_setz() */
2729 #endif
2730
2731 /* }}} */
2732
2733 /* {{{ s_mp_copy(sp, dp, count) */
2734
2735 #if MP_MACRO == 0
2736 /* Copy 'count' digits from sp to dp */
2737 void s_mp_copy(mp_digit *sp, mp_digit *dp, mp_size count)
2738 {
2739 #if MP_MEMCPY == 0
2740 int ix;
2741
2742 for(ix = 0; ix < count; ix++)
2743 dp[ix] = sp[ix];
2744 #else
2745 memcpy(dp, sp, count * sizeof(mp_digit));
2746 #endif
2747
2748 } /* end s_mp_copy() */
2749 #endif
2750
2751 /* }}} */
2752
2753 /* {{{ s_mp_alloc(nb, ni) */
2754
2755 #if MP_MACRO == 0
2756 /* Allocate ni records of nb bytes each, and return a pointer to that */
2757 void *s_mp_alloc(size_t nb, size_t ni)
2758 {
2759 return calloc(nb, ni);
2760
2761 } /* end s_mp_alloc() */
2762 #endif
2763
2764 /* }}} */
2765
2766 /* {{{ s_mp_free(ptr) */
2767
2768 #if MP_MACRO == 0
2769 /* Free the memory pointed to by ptr */
2770 void s_mp_free(void *ptr)
2771 {
2772 if(ptr)
2773 free(ptr);
2774
2775 } /* end s_mp_free() */
2776 #endif
2777
2778 /* }}} */
2779
2780 /* {{{ s_mp_clamp(mp) */
2781
2782 /* Remove leading zeroes from the given value */
2783 void s_mp_clamp(mp_int *mp)
2784 {
2785 mp_size du = USED(mp);
2786 mp_digit *zp = DIGITS(mp) + du - 1;
2787
2788 while(du > 1 && !*zp--)
2789 --du;
2790
2791 if(du == 1 && *zp == 0)
2792 SIGN(mp) = MP_ZPOS;
2793
2794 USED(mp) = du;
2795
2796 } /* end s_mp_clamp() */
2797
2798
2799 /* }}} */
2800
2801 /* {{{ s_mp_exch(a, b) */
2802
2803 /* Exchange the data for a and b; (b, a) = (a, b) */
2804 void s_mp_exch(mp_int *a, mp_int *b)
2805 {
2806 mp_int tmp;
2807
2808 tmp = *a;
2809 *a = *b;
2810 *b = tmp;
2811
2812 } /* end s_mp_exch() */
2813
2814 /* }}} */
2815
2816 /* }}} */
2817
2818 /* {{{ Arithmetic helpers */
2819
2820 /* {{{ s_mp_lshd(mp, p) */
2821
2822 /*
2823 Shift mp leftward by p digits, growing if needed, and zero-filling
2824 the in-shifted digits at the right end. This is a convenient
2825 alternative to multiplication by powers of the radix
2826 */
2827
2828 mp_err s_mp_lshd(mp_int *mp, mp_size p)
2829 {
2830 mp_err res;
2831 mp_size pos;
2832 mp_digit *dp;
2833 int ix;
2834
2835 if(p == 0)
2836 return MP_OKAY;
2837
2838 if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY)
2839 return res;
2840
2841 pos = USED(mp) - 1;
2842 dp = DIGITS(mp);
2843
2844 /* Shift all the significant figures over as needed */
2845 for(ix = pos - p; ix >= 0; ix--)
2846 dp[ix + p] = dp[ix];
2847
2848 /* Fill the bottom digits with zeroes */
2849 for(ix = 0; ix < p; ix++)
2850 dp[ix] = 0;
2851
2852 return MP_OKAY;
2853
2854 } /* end s_mp_lshd() */
2855
2856 /* }}} */
2857
2858 /* {{{ s_mp_rshd(mp, p) */
2859
2860 /*
2861 Shift mp rightward by p digits. Maintains the invariant that
2862 digits above the precision are all zero. Digits shifted off the
2863 end are lost. Cannot fail.
2864 */
2865
2866 void s_mp_rshd(mp_int *mp, mp_size p)
2867 {
2868 mp_size ix;
2869 mp_digit *dp;
2870
2871 if(p == 0)
2872 return;
2873
2874 /* Shortcut when all digits are to be shifted off */
2875 if(p >= USED(mp)) {
2876 s_mp_setz(DIGITS(mp), ALLOC(mp));
2877 USED(mp) = 1;
2878 SIGN(mp) = MP_ZPOS;
2879 return;
2880 }
2881
2882 /* Shift all the significant figures over as needed */
2883 dp = DIGITS(mp);
2884 for(ix = p; ix < USED(mp); ix++)
2885 dp[ix - p] = dp[ix];
2886
2887 /* Fill the top digits with zeroes */
2888 ix -= p;
2889 while(ix < USED(mp))
2890 dp[ix++] = 0;
2891
2892 /* Strip off any leading zeroes */
2893 s_mp_clamp(mp);
2894
2895 } /* end s_mp_rshd() */
2896
2897 /* }}} */
2898
2899 /* {{{ s_mp_div_2(mp) */
2900
2901 /* Divide by two -- take advantage of radix properties to do it fast */
2902 void s_mp_div_2(mp_int *mp)
2903 {
2904 s_mp_div_2d(mp, 1);
2905
2906 } /* end s_mp_div_2() */
2907
2908 /* }}} */
2909
2910 /* {{{ s_mp_mul_2(mp) */
2911
2912 mp_err s_mp_mul_2(mp_int *mp)
2913 {
2914 int ix;
2915 mp_digit kin = 0, kout, *dp = DIGITS(mp);
2916 mp_err res;
2917
2918 /* Shift digits leftward by 1 bit */
2919 for(ix = 0; ix < USED(mp); ix++) {
2920 kout = (dp[ix] >> (DIGIT_BIT - 1)) & 1;
2921 dp[ix] = (dp[ix] << 1) | kin;
2922
2923 kin = kout;
2924 }
2925
2926 /* Deal with rollover from last digit */
2927 if(kin) {
2928 if(ix >= ALLOC(mp)) {
2929 if((res = s_mp_grow(mp, ALLOC(mp) + 1)) != MP_OKAY)
2930 return res;
2931 dp = DIGITS(mp);
2932 }
2933
2934 dp[ix] = kin;
2935 USED(mp) += 1;
2936 }
2937
2938 return MP_OKAY;
2939
2940 } /* end s_mp_mul_2() */
2941
2942 /* }}} */
2943
2944 /* {{{ s_mp_mod_2d(mp, d) */
2945
2946 /*
2947 Remainder the integer by 2^d, where d is a number of bits. This
2948 amounts to a bitwise AND of the value, and does not require the full
2949 division code
2950 */
2951 void s_mp_mod_2d(mp_int *mp, mp_digit d)
2952 {
2953 unsigned int ndig = (d / DIGIT_BIT), nbit = (d % DIGIT_BIT);
2954 unsigned int ix;
2955 mp_digit dmask, *dp = DIGITS(mp);
2956
2957 if(ndig >= USED(mp))
2958 return;
2959
2960 /* Flush all the bits above 2^d in its digit */
2961 dmask = (1 << nbit) - 1;
2962 dp[ndig] &= dmask;
2963
2964 /* Flush all digits above the one with 2^d in it */
2965 for(ix = ndig + 1; ix < USED(mp); ix++)
2966 dp[ix] = 0;
2967
2968 s_mp_clamp(mp);
2969
2970 } /* end s_mp_mod_2d() */
2971
2972 /* }}} */
2973
2974 /* {{{ s_mp_mul_2d(mp, d) */
2975
2976 /*
2977 Multiply by the integer 2^d, where d is a number of bits. This
2978 amounts to a bitwise shift of the value, and does not require the
2979 full multiplication code.
2980 */
2981 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d)
2982 {
2983 mp_err res;
2984 mp_digit save, next, mask, *dp;
2985 mp_size used;
2986 int ix;
2987
2988 if((res = s_mp_lshd(mp, d / DIGIT_BIT)) != MP_OKAY)
2989 return res;
2990
2991 dp = DIGITS(mp); used = USED(mp);
2992 d %= DIGIT_BIT;
2993
2994 mask = (1 << d) - 1;
2995
2996 /* If the shift requires another digit, make sure we've got one to
2997 work with */
2998 if((dp[used - 1] >> (DIGIT_BIT - d)) & mask) {
2999 if((res = s_mp_grow(mp, used + 1)) != MP_OKAY)
3000 return res;
3001 dp = DIGITS(mp);
3002 }
3003
3004 /* Do the shifting... */
3005 save = 0;
3006 for(ix = 0; ix < used; ix++) {
3007 next = (dp[ix] >> (DIGIT_BIT - d)) & mask;
3008 dp[ix] = (dp[ix] << d) | save;
3009 save = next;
3010 }
3011
3012 /* If, at this point, we have a nonzero carryout into the next
3013 digit, we'll increase the size by one digit, and store it...
3014 */
3015 if(save) {
3016 dp[used] = save;
3017 USED(mp) += 1;
3018 }
3019
3020 s_mp_clamp(mp);
3021 return MP_OKAY;
3022
3023 } /* end s_mp_mul_2d() */
3024
3025 /* }}} */
3026
3027 /* {{{ s_mp_div_2d(mp, d) */
3028
3029 /*
3030 Divide the integer by 2^d, where d is a number of bits. This
3031 amounts to a bitwise shift of the value, and does not require the
3032 full division code (used in Barrett reduction, see below)
3033 */
3034 void s_mp_div_2d(mp_int *mp, mp_digit d)
3035 {
3036 int ix;
3037 mp_digit save, next, mask, *dp = DIGITS(mp);
3038
3039 s_mp_rshd(mp, d / DIGIT_BIT);
3040 d %= DIGIT_BIT;
3041
3042 mask = (1 << d) - 1;
3043
3044 save = 0;
3045 for(ix = USED(mp) - 1; ix >= 0; ix--) {
3046 next = dp[ix] & mask;
3047 dp[ix] = (dp[ix] >> d) | (save << (DIGIT_BIT - d));
3048 save = next;
3049 }
3050
3051 s_mp_clamp(mp);
3052
3053 } /* end s_mp_div_2d() */
3054
3055 /* }}} */
3056
3057 /* {{{ s_mp_norm(a, b) */
3058
3059 /*
3060 s_mp_norm(a, b)
3061
3062 Normalize a and b for division, where b is the divisor. In order
3063 that we might make good guesses for quotient digits, we want the
3064 leading digit of b to be at least half the radix, which we
3065 accomplish by multiplying a and b by a constant. This constant is
3066 returned (so that it can be divided back out of the remainder at the
3067 end of the division process).
3068
3069 We multiply by the smallest power of 2 that gives us a leading digit
3070 at least half the radix. By choosing a power of 2, we simplify the
3071 multiplication and division steps to simple shifts.
3072 */
3073 mp_digit s_mp_norm(mp_int *a, mp_int *b)
3074 {
3075 mp_digit t, d = 0;
3076
3077 t = DIGIT(b, USED(b) - 1);
3078 while(t < (RADIX / 2)) {
3079 t <<= 1;
3080 ++d;
3081 }
3082
3083 if(d != 0) {
3084 s_mp_mul_2d(a, d);
3085 s_mp_mul_2d(b, d);
3086 }
3087
3088 return d;
3089
3090 } /* end s_mp_norm() */
3091
3092 /* }}} */
3093
3094 /* }}} */
3095
3096 /* {{{ Primitive digit arithmetic */
3097
3098 /* {{{ s_mp_add_d(mp, d) */
3099
3100 /* Add d to |mp| in place */
3101 mp_err s_mp_add_d(mp_int *mp, mp_digit d) /* unsigned digit addition */
3102 {
3103 mp_word w, k = 0;
3104 mp_size ix = 1, used = USED(mp);
3105 mp_digit *dp = DIGITS(mp);
3106
3107 w = dp[0] + d;
3108 dp[0] = ACCUM(w);
3109 k = CARRYOUT(w);
3110
3111 while(ix < used && k) {
3112 w = dp[ix] + k;
3113 dp[ix] = ACCUM(w);
3114 k = CARRYOUT(w);
3115 ++ix;
3116 }
3117
3118 if(k != 0) {
3119 mp_err res;
3120
3121 if((res = s_mp_pad(mp, USED(mp) + 1)) != MP_OKAY)
3122 return res;
3123
3124 DIGIT(mp, ix) = k;
3125 }
3126
3127 return MP_OKAY;
3128
3129 } /* end s_mp_add_d() */
3130
3131 /* }}} */
3132
3133 /* {{{ s_mp_sub_d(mp, d) */
3134
3135 /* Subtract d from |mp| in place, assumes |mp| > d */
3136 mp_err s_mp_sub_d(mp_int *mp, mp_digit d) /* unsigned digit subtract */
3137 {
3138 mp_word w, b = 0;
3139 mp_size ix = 1, used = USED(mp);
3140 mp_digit *dp = DIGITS(mp);
3141
3142 /* Compute initial subtraction */
3143 w = (RADIX + dp[0]) - d;
3144 b = CARRYOUT(w) ? 0 : 1;
3145 dp[0] = ACCUM(w);
3146
3147 /* Propagate borrows leftward */
3148 while(b && ix < used) {
3149 w = (RADIX + dp[ix]) - b;
3150 b = CARRYOUT(w) ? 0 : 1;
3151 dp[ix] = ACCUM(w);
3152 ++ix;
3153 }
3154
3155 /* Remove leading zeroes */
3156 s_mp_clamp(mp);
3157
3158 /* If we have a borrow out, it's a violation of the input invariant */
3159 if(b)
3160 return MP_RANGE;
3161 else
3162 return MP_OKAY;
3163
3164 } /* end s_mp_sub_d() */
3165
3166 /* }}} */
3167
3168 /* {{{ s_mp_mul_d(a, d) */
3169
3170 /* Compute a = a * d, single digit multiplication */
3171 mp_err s_mp_mul_d(mp_int *a, mp_digit d)
3172 {
3173 mp_word w, k = 0;
3174 mp_size ix, max;
3175 mp_err res;
3176 mp_digit *dp = DIGITS(a);
3177
3178 /*
3179 Single-digit multiplication will increase the precision of the
3180 output by at most one digit. However, we can detect when this
3181 will happen -- if the high-order digit of a, times d, gives a
3182 two-digit result, then the precision of the result will increase;
3183 otherwise it won't. We use this fact to avoid calling s_mp_pad()
3184 unless absolutely necessary.
3185 */
3186 max = USED(a);
3187 w = dp[max - 1] * d;
3188 if(CARRYOUT(w) != 0) {
3189 if((res = s_mp_pad(a, max + 1)) != MP_OKAY)
3190 return res;
3191 dp = DIGITS(a);
3192 }
3193
3194 for(ix = 0; ix < max; ix++) {
3195 w = (dp[ix] * d) + k;
3196 dp[ix] = ACCUM(w);
3197 k = CARRYOUT(w);
3198 }
3199
3200 /* If there is a precision increase, take care of it here; the above
3201 test guarantees we have enough storage to do this safely.
3202 */
3203 if(k) {
3204 dp[max] = k;
3205 USED(a) = max + 1;
3206 }
3207
3208 s_mp_clamp(a);
3209
3210 return MP_OKAY;
3211
3212 } /* end s_mp_mul_d() */
3213
3214 /* }}} */
3215
3216 /* {{{ s_mp_div_d(mp, d, r) */
3217
3218 /*
3219 s_mp_div_d(mp, d, r)
3220
3221 Compute the quotient mp = mp / d and remainder r = mp mod d, for a
3222 single digit d. If r is null, the remainder will be discarded.
3223 */
3224
3225 mp_err s_mp_div_d(mp_int *mp, mp_digit d, mp_digit *r)
3226 {
3227 mp_word w = 0, t;
3228 mp_int quot;
3229 mp_err res;
3230 mp_digit *dp = DIGITS(mp), *qp;
3231 int ix;
3232
3233 if(d == 0)
3234 return MP_RANGE;
3235
3236 /* Make room for the quotient */
3237 if((res = mp_init_size(&quot, USED(mp))) != MP_OKAY)
3238 return res;
3239
3240 USED(&quot) = USED(mp); /* so clamping will work below */
3241 qp = DIGITS(&quot);
3242
3243 /* Divide without subtraction */
3244 for(ix = USED(mp) - 1; ix >= 0; ix--) {
3245 w = (w << DIGIT_BIT) | dp[ix];
3246
3247 if(w >= d) {
3248 t = w / d;
3249 w = w % d;
3250 } else {
3251 t = 0;
3252 }
3253
3254 qp[ix] = t;
3255 }
3256
3257 /* Deliver the remainder, if desired */
3258 if(r)
3259 *r = w;
3260
3261 s_mp_clamp(&quot);
3262 mp_exch(&quot, mp);
3263 mp_clear(&quot);
3264
3265 return MP_OKAY;
3266
3267 } /* end s_mp_div_d() */
3268
3269 /* }}} */
3270
3271 /* }}} */
3272
3273 /* {{{ Primitive full arithmetic */
3274
3275 /* {{{ s_mp_add(a, b) */
3276
3277 /* Compute a = |a| + |b| */
3278 mp_err s_mp_add(mp_int *a, mp_int *b) /* magnitude addition */
3279 {
3280 mp_word w = 0;
3281 mp_digit *pa, *pb;
3282 mp_size ix, used = USED(b);
3283 mp_err res;
3284
3285 /* Make sure a has enough precision for the output value */
3286 if((used > USED(a)) && (res = s_mp_pad(a, used)) != MP_OKAY)
3287 return res;
3288
3289 /*
3290 Add up all digits up to the precision of b. If b had initially
3291 the same precision as a, or greater, we took care of it by the
3292 padding step above, so there is no problem. If b had initially
3293 less precision, we'll have to make sure the carry out is duly
3294 propagated upward among the higher-order digits of the sum.
3295 */
3296 pa = DIGITS(a);
3297 pb = DIGITS(b);
3298 for(ix = 0; ix < used; ++ix) {
3299 w += *pa + *pb++;
3300 *pa++ = ACCUM(w);
3301 w = CARRYOUT(w);
3302 }
3303
3304 /* If we run out of 'b' digits before we're actually done, make
3305 sure the carries get propagated upward...
3306 */
3307 used = USED(a);
3308 while(w && ix < used) {
3309 w += *pa;
3310 *pa++ = ACCUM(w);
3311 w = CARRYOUT(w);
3312 ++ix;
3313 }
3314
3315 /* If there's an overall carry out, increase precision and include
3316 it. We could have done this initially, but why touch the memory
3317 allocator unless we're sure we have to?
3318 */
3319 if(w) {
3320 if((res = s_mp_pad(a, used + 1)) != MP_OKAY)
3321 return res;
3322
3323 DIGIT(a, ix) = w; /* pa may not be valid after s_mp_pad() call */
3324 }
3325
3326 return MP_OKAY;
3327
3328 } /* end s_mp_add() */
3329
3330 /* }}} */
3331
3332 /* {{{ s_mp_sub(a, b) */
3333
3334 /* Compute a = |a| - |b|, assumes |a| >= |b| */
3335 mp_err s_mp_sub(mp_int *a, mp_int *b) /* magnitude subtract */
3336 {
3337 mp_word w = 0;
3338 mp_digit *pa, *pb;
3339 mp_size ix, used = USED(b);
3340
3341 /*
3342 Subtract and propagate borrow. Up to the precision of b, this
3343 accounts for the digits of b; after that, we just make sure the
3344 carries get to the right place. This saves having to pad b out to
3345 the precision of a just to make the loops work right...
3346 */
3347 pa = DIGITS(a);
3348 pb = DIGITS(b);
3349
3350 for(ix = 0; ix < used; ++ix) {
3351 w = (RADIX + *pa) - w - *pb++;
3352 *pa++ = ACCUM(w);
3353 w = CARRYOUT(w) ? 0 : 1;
3354 }
3355
3356 used = USED(a);
3357 while(ix < used) {
3358 w = RADIX + *pa - w;
3359 *pa++ = ACCUM(w);
3360 w = CARRYOUT(w) ? 0 : 1;
3361 ++ix;
3362 }
3363
3364 /* Clobber any leading zeroes we created */
3365 s_mp_clamp(a);
3366
3367 /*
3368 If there was a borrow out, then |b| > |a| in violation
3369 of our input invariant. We've already done the work,
3370 but we'll at least complain about it...
3371 */
3372 if(w)
3373 return MP_RANGE;
3374 else
3375 return MP_OKAY;
3376
3377 } /* end s_mp_sub() */
3378
3379 /* }}} */
3380
3381 /* {{{ s_mp_mul(a, b) */
3382
3383 /* Compute a = |a| * |b| */
3384 mp_err s_mp_mul(mp_int *a, mp_int *b)
3385 {
3386 mp_word w, k = 0;
3387 mp_int tmp;
3388 mp_err res;
3389 mp_size ix, jx, ua = USED(a), ub = USED(b);
3390 mp_digit *pa, *pb, *pt, *pbt;
3391
3392 if((res = mp_init_size(&tmp, ua + ub)) != MP_OKAY)
3393 return res;
3394
3395 /* This has the effect of left-padding with zeroes... */
3396 USED(&tmp) = ua + ub;
3397
3398 /* We're going to need the base value each iteration */
3399 pbt = DIGITS(&tmp);
3400
3401 /* Outer loop: Digits of b */
3402
3403 pb = DIGITS(b);
3404 for(ix = 0; ix < ub; ++ix, ++pb) {
3405 if(*pb == 0)
3406 continue;
3407
3408 /* Inner product: Digits of a */
3409 pa = DIGITS(a);
3410 for(jx = 0; jx < ua; ++jx, ++pa) {
3411 pt = pbt + ix + jx;
3412 w = *pb * *pa + k + *pt;
3413 *pt = ACCUM(w);
3414 k = CARRYOUT(w);
3415 }
3416
3417 pbt[ix + jx] = k;
3418 k = 0;
3419 }
3420
3421 s_mp_clamp(&tmp);
3422 s_mp_exch(&tmp, a);
3423
3424 mp_clear(&tmp);
3425
3426 return MP_OKAY;
3427
3428 } /* end s_mp_mul() */
3429
3430 /* }}} */
3431
3432 /* {{{ s_mp_kmul(a, b, out, len) */
3433
3434 #if 0
3435 void s_mp_kmul(mp_digit *a, mp_digit *b, mp_digit *out, mp_size len)
3436 {
3437 mp_word w, k = 0;
3438 mp_size ix, jx;
3439 mp_digit *pa, *pt;
3440
3441 for(ix = 0; ix < len; ++ix, ++b) {
3442 if(*b == 0)
3443 continue;
3444
3445 pa = a;
3446 for(jx = 0; jx < len; ++jx, ++pa) {
3447 pt = out + ix + jx;
3448 w = *b * *pa + k + *pt;
3449 *pt = ACCUM(w);
3450 k = CARRYOUT(w);
3451 }
3452
3453 out[ix + jx] = k;
3454 k = 0;
3455 }
3456
3457 } /* end s_mp_kmul() */
3458 #endif
3459
3460 /* }}} */
3461
3462 /* {{{ s_mp_sqr(a) */
3463
3464 /*
3465 Computes the square of a, in place. This can be done more
3466 efficiently than a general multiplication, because many of the
3467 computation steps are redundant when squaring. The inner product
3468 step is a bit more complicated, but we save a fair number of
3469 iterations of the multiplication loop.
3470 */
3471 #if MP_SQUARE
3472 mp_err s_mp_sqr(mp_int *a)
3473 {
3474 mp_word w, k = 0;
3475 mp_int tmp;
3476 mp_err res;
3477 mp_size ix, jx, kx, used = USED(a);
3478 mp_digit *pa1, *pa2, *pt, *pbt;
3479
3480 if((res = mp_init_size(&tmp, 2 * used)) != MP_OKAY)
3481 return res;
3482
3483 /* Left-pad with zeroes */
3484 USED(&tmp) = 2 * used;
3485
3486 /* We need the base value each time through the loop */
3487 pbt = DIGITS(&tmp);
3488
3489 pa1 = DIGITS(a);
3490 for(ix = 0; ix < used; ++ix, ++pa1) {
3491 if(*pa1 == 0)
3492 continue;
3493
3494 w = DIGIT(&tmp, ix + ix) + (*pa1 * *pa1);
3495
3496 pbt[ix + ix] = ACCUM(w);
3497 k = CARRYOUT(w);
3498
3499 /*
3500 The inner product is computed as:
3501
3502 (C, S) = t[i,j] + 2 a[i] a[j] + C
3503
3504 This can overflow what can be represented in an mp_word, and
3505 since C arithmetic does not provide any way to check for
3506 overflow, we have to check explicitly for overflow conditions
3507 before they happen.
3508 */
3509 for(jx = ix + 1, pa2 = DIGITS(a) + jx; jx < used; ++jx, ++pa2) {
3510 mp_word u = 0, v;
3511
3512 /* Store this in a temporary to avoid indirections later */
3513 pt = pbt + ix + jx;
3514
3515 /* Compute the multiplicative step */
3516 w = *pa1 * *pa2;
3517
3518 /* If w is more than half MP_WORD_MAX, the doubling will
3519 overflow, and we need to record a carry out into the next
3520 word */
3521 u = (w >> (MP_WORD_BIT - 1)) & 1;
3522
3523 /* Double what we've got, overflow will be ignored as defined
3524 for C arithmetic (we've already noted if it is to occur)
3525 */
3526 w *= 2;
3527
3528 /* Compute the additive step */
3529 v = *pt + k;
3530
3531 /* If we do not already have an overflow carry, check to see
3532 if the addition will cause one, and set the carry out if so
3533 */
3534 u |= ((MP_WORD_MAX - v) < w);
3535
3536 /* Add in the rest, again ignoring overflow */
3537 w += v;
3538
3539 /* Set the i,j digit of the output */
3540 *pt = ACCUM(w);
3541
3542 /* Save carry information for the next iteration of the loop.
3543 This is why k must be an mp_word, instead of an mp_digit */
3544 k = CARRYOUT(w) | (u << DIGIT_BIT);
3545
3546 } /* for(jx ...) */
3547
3548 /* Set the last digit in the cycle and reset the carry */
3549 k = DIGIT(&tmp, ix + jx) + k;
3550 pbt[ix + jx] = ACCUM(k);
3551 k = CARRYOUT(k);
3552
3553 /* If we are carrying out, propagate the carry to the next digit
3554 in the output. This may cascade, so we have to be somewhat
3555 circumspect -- but we will have enough precision in the output
3556 that we won't overflow
3557 */
3558 kx = 1;
3559 while(k) {
3560 k = pbt[ix + jx + kx] + 1;
3561 pbt[ix + jx + kx] = ACCUM(k);
3562 k = CARRYOUT(k);
3563 ++kx;
3564 }
3565 } /* for(ix ...) */
3566
3567 s_mp_clamp(&tmp);
3568 s_mp_exch(&tmp, a);
3569
3570 mp_clear(&tmp);
3571
3572 return MP_OKAY;
3573
3574 } /* end s_mp_sqr() */
3575 #endif
3576
3577 /* }}} */
3578
3579 /* {{{ s_mp_div(a, b) */
3580
3581 /*
3582 s_mp_div(a, b)
3583
3584 Compute a = a / b and b = a mod b. Assumes b > a.
3585 */
3586
3587 mp_err s_mp_div(mp_int *a, mp_int *b)
3588 {
3589 mp_int quot, rem, t;
3590 mp_word q;
3591 mp_err res;
3592 mp_digit d;
3593 int ix;
3594
3595 if(mp_cmp_z(b) == 0)
3596 return MP_RANGE;
3597
3598 /* Shortcut if b is power of two */
3599 if((ix = s_mp_ispow2(b)) >= 0) {
3600 mp_copy(a, b); /* need this for remainder */
3601 s_mp_div_2d(a, (mp_digit)ix);
3602 s_mp_mod_2d(b, (mp_digit)ix);
3603
3604 return MP_OKAY;
3605 }
3606
3607 /* Allocate space to store the quotient */
3608 if((res = mp_init_size(&quot, USED(a))) != MP_OKAY)
3609 return res;
3610
3611 /* A working temporary for division */
3612 if((res = mp_init_size(&t, USED(a))) != MP_OKAY)
3613 goto T;
3614
3615 /* Allocate space for the remainder */
3616 if((res = mp_init_size(&rem, USED(a))) != MP_OKAY)
3617 goto REM;
3618
3619 /* Normalize to optimize guessing */
3620 d = s_mp_norm(a, b);
3621
3622 /* Perform the division itself...woo! */
3623 ix = USED(a) - 1;
3624
3625 while(ix >= 0) {
3626 /* Find a partial substring of a which is at least b */
3627 while(s_mp_cmp(&rem, b) < 0 && ix >= 0) {
3628 if((res = s_mp_lshd(&rem, 1)) != MP_OKAY)
3629 goto CLEANUP;
3630
3631 if((res = s_mp_lshd(&quot, 1)) != MP_OKAY)
3632 goto CLEANUP;
3633
3634 DIGIT(&rem, 0) = DIGIT(a, ix);
3635 s_mp_clamp(&rem);
3636 --ix;
3637 }
3638
3639 /* If we didn't find one, we're finished dividing */
3640 if(s_mp_cmp(&rem, b) < 0)
3641 break;
3642
3643 /* Compute a guess for the next quotient digit */
3644 q = DIGIT(&rem, USED(&rem) - 1);
3645 if(q <= DIGIT(b, USED(b) - 1) && USED(&rem) > 1)
3646 q = (q << DIGIT_BIT) | DIGIT(&rem, USED(&rem) - 2);
3647
3648 q /= DIGIT(b, USED(b) - 1);
3649
3650 /* The guess can be as much as RADIX + 1 */
3651 if(q >= RADIX)
3652 q = RADIX - 1;
3653
3654 /* See what that multiplies out to */
3655 mp_copy(b, &t);
3656 if((res = s_mp_mul_d(&t, q)) != MP_OKAY)
3657 goto CLEANUP;
3658
3659 /*
3660 If it's too big, back it off. We should not have to do this
3661 more than once, or, in rare cases, twice. Knuth describes a
3662 method by which this could be reduced to a maximum of once, but
3663 I didn't implement that here.
3664 */
3665 while(s_mp_cmp(&t, &rem) > 0) {
3666 --q;
3667 s_mp_sub(&t, b);
3668 }
3669
3670 /* At this point, q should be the right next digit */
3671 if((res = s_mp_sub(&rem, &t)) != MP_OKAY)
3672 goto CLEANUP;
3673
3674 /*
3675 Include the digit in the quotient. We allocated enough memory
3676 for any quotient we could ever possibly get, so we should not
3677 have to check for failures here
3678 */
3679 DIGIT(&quot, 0) = q;
3680 }
3681
3682 /* Denormalize remainder */
3683 if(d != 0)
3684 s_mp_div_2d(&rem, d);
3685
3686 s_mp_clamp(&quot);
3687 s_mp_clamp(&rem);
3688
3689 /* Copy quotient back to output */
3690 s_mp_exch(&quot, a);
3691
3692 /* Copy remainder back to output */
3693 s_mp_exch(&rem, b);
3694
3695 CLEANUP:
3696 mp_clear(&rem);
3697 REM:
3698 mp_clear(&t);
3699 T:
3700 mp_clear(&quot);
3701
3702 return res;
3703
3704 } /* end s_mp_div() */
3705
3706 /* }}} */
3707
3708 /* {{{ s_mp_2expt(a, k) */
3709
3710 mp_err s_mp_2expt(mp_int *a, mp_digit k)
3711 {
3712 mp_err res;
3713 mp_size dig, bit;
3714
3715 dig = k / DIGIT_BIT;
3716 bit = k % DIGIT_BIT;
3717
3718 mp_zero(a);
3719 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY)
3720 return res;
3721
3722 DIGIT(a, dig) |= (1 << bit);
3723
3724 return MP_OKAY;
3725
3726 } /* end s_mp_2expt() */
3727
3728 /* }}} */
3729
3730 /* {{{ s_mp_reduce(x, m, mu) */
3731
3732 /*
3733 Compute Barrett reduction, x (mod m), given a precomputed value for
3734 mu = b^2k / m, where b = RADIX and k = #digits(m). This should be
3735 faster than straight division, when many reductions by the same
3736 value of m are required (such as in modular exponentiation). This
3737 can nearly halve the time required to do modular exponentiation,
3738 as compared to using the full integer divide to reduce.
3739
3740 This algorithm was derived from the _Handbook of Applied
3741 Cryptography_ by Menezes, Oorschot and VanStone, Ch. 14,
3742 pp. 603-604.
3743 */
3744
3745 mp_err s_mp_reduce(mp_int *x, mp_int *m, mp_int *mu)
3746 {
3747 mp_int q;
3748 mp_err res;
3749 mp_size um = USED(m);
3750
3751 if((res = mp_init_copy(&q, x)) != MP_OKAY)
3752 return res;
3753
3754 s_mp_rshd(&q, um - 1); /* q1 = x / b^(k-1) */
3755 s_mp_mul(&q, mu); /* q2 = q1 * mu */
3756 s_mp_rshd(&q, um + 1); /* q3 = q2 / b^(k+1) */
3757
3758 /* x = x mod b^(k+1), quick (no division) */
3759 s_mp_mod_2d(x, DIGIT_BIT * (um + 1));
3760
3761 /* q = q * m mod b^(k+1), quick (no division) */
3762 s_mp_mul(&q, m);
3763 s_mp_mod_2d(&q, DIGIT_BIT * (um + 1));
3764
3765 /* x = x - q */
3766 if((res = mp_sub(x, &q, x)) != MP_OKAY)
3767 goto CLEANUP;
3768
3769 /* If x < 0, add b^(k+1) to it */
3770 if(mp_cmp_z(x) < 0) {
3771 mp_set(&q, 1);
3772 if((res = s_mp_lshd(&q, um + 1)) != MP_OKAY)
3773 goto CLEANUP;
3774 if((res = mp_add(x, &q, x)) != MP_OKAY)
3775 goto CLEANUP;
3776 }
3777
3778 /* Back off if it's too big */
3779 while(mp_cmp(x, m) >= 0) {
3780 if((res = s_mp_sub(x, m)) != MP_OKAY)
3781 break;
3782 }
3783
3784 CLEANUP:
3785 mp_clear(&q);
3786
3787 return res;
3788
3789 } /* end s_mp_reduce() */
3790
3791 /* }}} */
3792
3793 /* }}} */
3794
3795 /* {{{ Primitive comparisons */
3796
3797 /* {{{ s_mp_cmp(a, b) */
3798
3799 /* Compare |a| <=> |b|, return 0 if equal, <0 if a<b, >0 if a>b */
3800 int s_mp_cmp(mp_int *a, mp_int *b)
3801 {
3802 mp_size ua = USED(a), ub = USED(b);
3803
3804 if(ua > ub)
3805 return MP_GT;
3806 else if(ua < ub)
3807 return MP_LT;
3808 else {
3809 int ix = ua - 1;
3810 mp_digit *ap = DIGITS(a) + ix, *bp = DIGITS(b) + ix;
3811
3812 while(ix >= 0) {
3813 if(*ap > *bp)
3814 return MP_GT;
3815 else if(*ap < *bp)
3816 return MP_LT;
3817
3818 --ap; --bp; --ix;
3819 }
3820
3821 return MP_EQ;
3822 }
3823
3824 } /* end s_mp_cmp() */
3825
3826 /* }}} */
3827
3828 /* {{{ s_mp_cmp_d(a, d) */
3829
3830 /* Compare |a| <=> d, return 0 if equal, <0 if a<d, >0 if a>d */
3831 int s_mp_cmp_d(mp_int *a, mp_digit d)
3832 {
3833 mp_size ua = USED(a);
3834 mp_digit *ap = DIGITS(a);
3835
3836 if(ua > 1)
3837 return MP_GT;
3838
3839 if(*ap < d)
3840 return MP_LT;
3841 else if(*ap > d)
3842 return MP_GT;
3843 else
3844 return MP_EQ;
3845
3846 } /* end s_mp_cmp_d() */
3847
3848 /* }}} */
3849
3850 /* {{{ s_mp_ispow2(v) */
3851
3852 /*
3853 Returns -1 if the value is not a power of two; otherwise, it returns
3854 k such that v = 2^k, i.e. lg(v).
3855 */
3856 int s_mp_ispow2(mp_int *v)
3857 {
3858 mp_digit d, *dp;
3859 mp_size uv = USED(v);
3860 int extra = 0, ix;
3861
3862 d = DIGIT(v, uv - 1); /* most significant digit of v */
3863
3864 while(d && ((d & 1) == 0)) {
3865 d >>= 1;
3866 ++extra;
3867 }
3868
3869 if(d == 1) {
3870 ix = uv - 2;
3871 dp = DIGITS(v) + ix;
3872
3873 while(ix >= 0) {
3874 if(*dp)
3875 return -1; /* not a power of two */
3876
3877 --dp; --ix;
3878 }
3879
3880 return ((uv - 1) * DIGIT_BIT) + extra;
3881 }
3882
3883 return -1;
3884
3885 } /* end s_mp_ispow2() */
3886
3887 /* }}} */
3888
3889 /* {{{ s_mp_ispow2d(d) */
3890
3891 int s_mp_ispow2d(mp_digit d)
3892 {
3893 int pow = 0;
3894
3895 while((d & 1) == 0) {
3896 ++pow; d >>= 1;
3897 }
3898
3899 if(d == 1)
3900 return pow;
3901
3902 return -1;
3903
3904 } /* end s_mp_ispow2d() */
3905
3906 /* }}} */
3907
3908 /* }}} */
3909
3910 /* {{{ Primitive I/O helpers */
3911
3912 /* {{{ s_mp_tovalue(ch, r) */
3913
3914 /*
3915 Convert the given character to its digit value, in the given radix.
3916 If the given character is not understood in the given radix, -1 is
3917 returned. Otherwise the digit's numeric value is returned.
3918
3919 The results will be odd if you use a radix < 2 or > 62, you are
3920 expected to know what you're up to.
3921 */
3922 int s_mp_tovalue(char ch, int r)
3923 {
3924 int val, xch;
3925
3926 if(r > 36)
3927 xch = ch;
3928 else
3929 xch = toupper(ch);
3930
3931 if(isdigit(xch))
3932 val = xch - '0';
3933 else if(isupper(xch))
3934 val = xch - 'A' + 10;
3935 else if(islower(xch))
3936 val = xch - 'a' + 36;
3937 else if(xch == '+')
3938 val = 62;
3939 else if(xch == '/')
3940 val = 63;
3941 else
3942 return -1;
3943
3944 if(val < 0 || val >= r)
3945 return -1;
3946
3947 return val;
3948
3949 } /* end s_mp_tovalue() */
3950
3951 /* }}} */
3952
3953 /* {{{ s_mp_todigit(val, r, low) */
3954
3955 /*
3956 Convert val to a radix-r digit, if possible. If val is out of range
3957 for r, returns zero. Otherwise, returns an ASCII character denoting
3958 the value in the given radix.
3959
3960 The results may be odd if you use a radix < 2 or > 64, you are
3961 expected to know what you're doing.
3962 */
3963
3964 char s_mp_todigit(int val, int r, int low)
3965 {
3966 char ch;
3967
3968 if(val < 0 || val >= r)
3969 return 0;
3970
3971 ch = s_dmap_1[val];
3972
3973 if(r <= 36 && low)
3974 ch = tolower(ch);
3975
3976 return ch;
3977
3978 } /* end s_mp_todigit() */
3979
3980 /* }}} */
3981
3982 /* {{{ s_mp_outlen(bits, radix) */
3983
3984 /*
3985 Return an estimate for how long a string is needed to hold a radix
3986 r representation of a number with 'bits' significant bits.
3987
3988 Does not include space for a sign or a NUL terminator.
3989 */
3990 int s_mp_outlen(int bits, int r)
3991 {
3992 return (int)((double)bits * LOG_V_2(r) + 0.5);
3993
3994 } /* end s_mp_outlen() */
3995
3996 /* }}} */
3997
3998 /* }}} */
3999
4000 /*------------------------------------------------------------------------*/
4001 /* HERE THERE BE DRAGONS */

Mercurial Repository: pidgin3

mercurial