Mercurial > dropbear
comparison libtommath/mtest/mpi.c @ 1436:60fc6476e044
Update to libtommath v1.0
author | Matt Johnston <matt@ucc.asn.au> |
---|---|
date | Sat, 24 Jun 2017 22:37:14 +0800 |
parents | 5ff8218bcee9 |
children |
comparison
equal
deleted
inserted
replaced
1435:f849a5ca2efc | 1436:60fc6476e044 |
---|---|
4 by Michael J. Fromberger <[email protected]> | 4 by Michael J. Fromberger <[email protected]> |
5 Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved | 5 Copyright (C) 1998 Michael J. Fromberger, All Rights Reserved |
6 | 6 |
7 Arbitrary precision integer arithmetic library | 7 Arbitrary precision integer arithmetic library |
8 | 8 |
9 $Id: mpi.c,v 1.2 2005/05/05 14:38:47 tom Exp $ | 9 $Id$ |
10 */ | 10 */ |
11 | 11 |
12 #include "mpi.h" | 12 #include "mpi.h" |
13 #include <stdlib.h> | 13 #include <stdlib.h> |
14 #include <string.h> | 14 #include <string.h> |
20 #define DIAG(T,V) {fprintf(stderr,T);mp_print(V,stderr);fputc('\n',stderr);} | 20 #define DIAG(T,V) {fprintf(stderr,T);mp_print(V,stderr);fputc('\n',stderr);} |
21 #else | 21 #else |
22 #define DIAG(T,V) | 22 #define DIAG(T,V) |
23 #endif | 23 #endif |
24 | 24 |
25 /* | 25 /* |
26 If MP_LOGTAB is not defined, use the math library to compute the | 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. | 27 logarithms on the fly. Otherwise, use the static table below. |
28 Pick which works best for your system. | 28 Pick which works best for your system. |
29 */ | 29 */ |
30 #if MP_LOGTAB | 30 #if MP_LOGTAB |
31 | 31 |
32 /* {{{ s_logv_2[] - log table for 2 in various bases */ | 32 /* {{{ s_logv_2[] - log table for 2 in various bases */ |
33 | 33 |
34 /* | 34 /* |
35 A table of the logs of 2 for various bases (the 0 and 1 entries of | 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). | 36 this table are meaningless and should not be referenced). |
37 | 37 |
38 This table is used to compute output lengths for the mp_toradix() | 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) | 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 | 40 digits, we estimate the output size by taking the least integer |
41 greater than log_r(n), where: | 41 greater than log_r(n), where: |
42 | 42 |
43 log_r(n) = log_2(n) * log_r(2) | 43 log_r(n) = log_2(n) * log_r(2) |
44 | 44 |
45 This table, therefore, is a table of log_r(2) for 2 <= r <= 36, | 45 This table, therefore, is a table of log_r(2) for 2 <= r <= 36, |
46 which are the output bases supported. | 46 which are the output bases supported. |
47 */ | 47 */ |
48 | 48 |
49 #include "logtab.h" | 49 #include "logtab.h" |
50 | 50 |
51 /* }}} */ | 51 /* }}} */ |
102 }; | 102 }; |
103 | 103 |
104 /* Value to digit maps for radix conversion */ | 104 /* Value to digit maps for radix conversion */ |
105 | 105 |
106 /* s_dmap_1 - standard digits and letters */ | 106 /* s_dmap_1 - standard digits and letters */ |
107 static const char *s_dmap_1 = | 107 static const char *s_dmap_1 = |
108 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/"; | 108 "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz+/"; |
109 | 109 |
110 #if 0 | 110 #if 0 |
111 /* s_dmap_2 - base64 ordering for digits */ | 111 /* s_dmap_2 - base64 ordering for digits */ |
112 static const char *s_dmap_2 = | 112 static const char *s_dmap_2 = |
115 | 115 |
116 /* }}} */ | 116 /* }}} */ |
117 | 117 |
118 /* {{{ Static function declarations */ | 118 /* {{{ Static function declarations */ |
119 | 119 |
120 /* | 120 /* |
121 If MP_MACRO is false, these will be defined as actual functions; | 121 If MP_MACRO is false, these will be defined as actual functions; |
122 otherwise, suitable macro definitions will be used. This works | 122 otherwise, suitable macro definitions will be used. This works |
123 around the fact that ANSI C89 doesn't support an 'inline' keyword | 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 | 124 (although I hear C9x will ... about bloody time). At present, the |
125 macro definitions are identical to the function bodies, but they'll | 125 macro definitions are identical to the function bodies, but they'll |
256 } | 256 } |
257 | 257 |
258 return MP_OKAY; | 258 return MP_OKAY; |
259 | 259 |
260 CLEANUP: | 260 CLEANUP: |
261 while(--pos >= 0) | 261 while(--pos >= 0) |
262 mp_clear(&mp[pos]); | 262 mp_clear(&mp[pos]); |
263 | 263 |
264 return res; | 264 return res; |
265 | 265 |
266 } /* end mp_init_array() */ | 266 } /* end mp_init_array() */ |
353 usual | 353 usual |
354 */ | 354 */ |
355 if(ALLOC(to) >= USED(from)) { | 355 if(ALLOC(to) >= USED(from)) { |
356 s_mp_setz(DIGITS(to) + USED(from), 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)); | 357 s_mp_copy(DIGITS(from), DIGITS(to), USED(from)); |
358 | 358 |
359 } else { | 359 } else { |
360 if((tmp = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL) | 360 if((tmp = s_mp_alloc(USED(from), sizeof(mp_digit))) == NULL) |
361 return MP_MEM; | 361 return MP_MEM; |
362 | 362 |
363 s_mp_copy(DIGITS(from), tmp, USED(from)); | 363 s_mp_copy(DIGITS(from), tmp, USED(from)); |
443 | 443 |
444 void mp_clear_array(mp_int mp[], int count) | 444 void mp_clear_array(mp_int mp[], int count) |
445 { | 445 { |
446 ARGCHK(mp != NULL && count > 0, MP_BADARG); | 446 ARGCHK(mp != NULL && count > 0, MP_BADARG); |
447 | 447 |
448 while(--count >= 0) | 448 while(--count >= 0) |
449 mp_clear(&mp[count]); | 449 mp_clear(&mp[count]); |
450 | 450 |
451 } /* end mp_clear_array() */ | 451 } /* end mp_clear_array() */ |
452 | 452 |
453 /* }}} */ | 453 /* }}} */ |
454 | 454 |
455 /* {{{ mp_zero(mp) */ | 455 /* {{{ mp_zero(mp) */ |
456 | 456 |
457 /* | 457 /* |
458 mp_zero(mp) | 458 mp_zero(mp) |
459 | 459 |
460 Set mp to zero. Does not change the allocated size of the structure, | 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) | 461 and therefore cannot fail (except on a bad argument, which we ignore) |
462 */ | 462 */ |
463 void mp_zero(mp_int *mp) | 463 void mp_zero(mp_int *mp) |
504 for(ix = sizeof(long) - 1; ix >= 0; ix--) { | 504 for(ix = sizeof(long) - 1; ix >= 0; ix--) { |
505 | 505 |
506 if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY) | 506 if((res = s_mp_mul_2d(mp, CHAR_BIT)) != MP_OKAY) |
507 return res; | 507 return res; |
508 | 508 |
509 res = s_mp_add_d(mp, | 509 res = s_mp_add_d(mp, |
510 (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX)); | 510 (mp_digit)((v >> (ix * CHAR_BIT)) & UCHAR_MAX)); |
511 if(res != MP_OKAY) | 511 if(res != MP_OKAY) |
512 return res; | 512 return res; |
513 | 513 |
514 } | 514 } |
839 ARGCHK(a != NULL && b != NULL, MP_BADARG); | 839 ARGCHK(a != NULL && b != NULL, MP_BADARG); |
840 | 840 |
841 if((res = mp_copy(a, b)) != MP_OKAY) | 841 if((res = mp_copy(a, b)) != MP_OKAY) |
842 return res; | 842 return res; |
843 | 843 |
844 if(s_mp_cmp_d(b, 0) == MP_EQ) | 844 if(s_mp_cmp_d(b, 0) == MP_EQ) |
845 SIGN(b) = MP_ZPOS; | 845 SIGN(b) = MP_ZPOS; |
846 else | 846 else |
847 SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG; | 847 SIGN(b) = (SIGN(b) == MP_NEG) ? MP_ZPOS : MP_NEG; |
848 | 848 |
849 return MP_OKAY; | 849 return MP_OKAY; |
850 | 850 |
851 } /* end mp_neg() */ | 851 } /* end mp_neg() */ |
868 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); | 868 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); |
869 | 869 |
870 if(SIGN(a) == SIGN(b)) { /* same sign: add values, keep sign */ | 870 if(SIGN(a) == SIGN(b)) { /* same sign: add values, keep sign */ |
871 | 871 |
872 /* Commutativity of addition lets us do this in either order, | 872 /* Commutativity of addition lets us do this in either order, |
873 so we avoid having to use a temporary even if the result | 873 so we avoid having to use a temporary even if the result |
874 is supposed to replace the output | 874 is supposed to replace the output |
875 */ | 875 */ |
876 if(c == b) { | 876 if(c == b) { |
877 if((res = s_mp_add(c, a)) != MP_OKAY) | 877 if((res = s_mp_add(c, a)) != MP_OKAY) |
878 return res; | 878 return res; |
879 } else { | 879 } else { |
880 if(c != a && (res = mp_copy(a, c)) != MP_OKAY) | 880 if(c != a && (res = mp_copy(a, c)) != MP_OKAY) |
881 return res; | 881 return res; |
882 | 882 |
883 if((res = s_mp_add(c, b)) != MP_OKAY) | 883 if((res = s_mp_add(c, b)) != MP_OKAY) |
884 return res; | 884 return res; |
885 } | 885 } |
886 | 886 |
887 } else if((cmp = s_mp_cmp(a, b)) > 0) { /* different sign: a > b */ | 887 } else if((cmp = s_mp_cmp(a, b)) > 0) { /* different sign: a > b */ |
888 | 888 |
889 /* If the output is going to be clobbered, we will use a temporary | 889 /* If the output is going to be clobbered, we will use a temporary |
890 variable; otherwise, we'll do it without touching the memory | 890 variable; otherwise, we'll do it without touching the memory |
891 allocator at all, if possible | 891 allocator at all, if possible |
892 */ | 892 */ |
893 if(c == b) { | 893 if(c == b) { |
894 mp_int tmp; | 894 mp_int tmp; |
895 | 895 |
1017 } | 1017 } |
1018 s_mp_exch(&tmp, c); | 1018 s_mp_exch(&tmp, c); |
1019 mp_clear(&tmp); | 1019 mp_clear(&tmp); |
1020 | 1020 |
1021 } else { | 1021 } else { |
1022 if(c != b && ((res = mp_copy(b, c)) != MP_OKAY)) | 1022 if(c != b && ((res = mp_copy(b, c)) != MP_OKAY)) |
1023 return res; | 1023 return res; |
1024 | 1024 |
1025 if((res = s_mp_sub(c, a)) != MP_OKAY) | 1025 if((res = s_mp_sub(c, a)) != MP_OKAY) |
1026 return res; | 1026 return res; |
1027 } | 1027 } |
1064 return res; | 1064 return res; |
1065 | 1065 |
1066 if((res = s_mp_mul(c, b)) != MP_OKAY) | 1066 if((res = s_mp_mul(c, b)) != MP_OKAY) |
1067 return res; | 1067 return res; |
1068 } | 1068 } |
1069 | 1069 |
1070 if(sgn == MP_ZPOS || s_mp_cmp_d(c, 0) == MP_EQ) | 1070 if(sgn == MP_ZPOS || s_mp_cmp_d(c, 0) == MP_EQ) |
1071 SIGN(c) = MP_ZPOS; | 1071 SIGN(c) = MP_ZPOS; |
1072 else | 1072 else |
1073 SIGN(c) = sgn; | 1073 SIGN(c) = sgn; |
1074 | 1074 |
1075 return MP_OKAY; | 1075 return MP_OKAY; |
1076 | 1076 |
1077 } /* end mp_mul() */ | 1077 } /* end mp_mul() */ |
1078 | 1078 |
1079 /* }}} */ | 1079 /* }}} */ |
1158 if(r) { | 1158 if(r) { |
1159 if((res = mp_copy(a, r)) != MP_OKAY) | 1159 if((res = mp_copy(a, r)) != MP_OKAY) |
1160 return res; | 1160 return res; |
1161 } | 1161 } |
1162 | 1162 |
1163 if(q) | 1163 if(q) |
1164 mp_zero(q); | 1164 mp_zero(q); |
1165 | 1165 |
1166 return MP_OKAY; | 1166 return MP_OKAY; |
1167 | 1167 |
1168 } else if(cmp == 0) { | 1168 } else if(cmp == 0) { |
1204 SIGN(&qtmp) = MP_ZPOS; | 1204 SIGN(&qtmp) = MP_ZPOS; |
1205 if(s_mp_cmp_d(&rtmp, 0) == MP_EQ) | 1205 if(s_mp_cmp_d(&rtmp, 0) == MP_EQ) |
1206 SIGN(&rtmp) = MP_ZPOS; | 1206 SIGN(&rtmp) = MP_ZPOS; |
1207 | 1207 |
1208 /* Copy output, if it is needed */ | 1208 /* Copy output, if it is needed */ |
1209 if(q) | 1209 if(q) |
1210 s_mp_exch(&qtmp, q); | 1210 s_mp_exch(&qtmp, q); |
1211 | 1211 |
1212 if(r) | 1212 if(r) |
1213 s_mp_exch(&rtmp, r); | 1213 s_mp_exch(&rtmp, r); |
1214 | 1214 |
1215 CLEANUP: | 1215 CLEANUP: |
1216 mp_clear(&rtmp); | 1216 mp_clear(&rtmp); |
1217 mp_clear(&qtmp); | 1217 mp_clear(&qtmp); |
1262 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c) | 1262 mp_err mp_expt(mp_int *a, mp_int *b, mp_int *c) |
1263 { | 1263 { |
1264 mp_int s, x; | 1264 mp_int s, x; |
1265 mp_err res; | 1265 mp_err res; |
1266 mp_digit d; | 1266 mp_digit d; |
1267 int dig, bit; | 1267 unsigned int bit, dig; |
1268 | 1268 |
1269 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); | 1269 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); |
1270 | 1270 |
1271 if(mp_cmp_z(b) < 0) | 1271 if(mp_cmp_z(b) < 0) |
1272 return MP_RANGE; | 1272 return MP_RANGE; |
1284 d = DIGIT(b, dig); | 1284 d = DIGIT(b, dig); |
1285 | 1285 |
1286 /* Loop over bits of each non-maximal digit */ | 1286 /* Loop over bits of each non-maximal digit */ |
1287 for(bit = 0; bit < DIGIT_BIT; bit++) { | 1287 for(bit = 0; bit < DIGIT_BIT; bit++) { |
1288 if(d & 1) { | 1288 if(d & 1) { |
1289 if((res = s_mp_mul(&s, &x)) != MP_OKAY) | 1289 if((res = s_mp_mul(&s, &x)) != MP_OKAY) |
1290 goto CLEANUP; | 1290 goto CLEANUP; |
1291 } | 1291 } |
1292 | 1292 |
1293 d >>= 1; | 1293 d >>= 1; |
1294 | 1294 |
1295 if((res = s_mp_sqr(&x)) != MP_OKAY) | 1295 if((res = s_mp_sqr(&x)) != MP_OKAY) |
1296 goto CLEANUP; | 1296 goto CLEANUP; |
1297 } | 1297 } |
1298 } | 1298 } |
1299 | 1299 |
1309 d >>= 1; | 1309 d >>= 1; |
1310 | 1310 |
1311 if((res = s_mp_sqr(&x)) != MP_OKAY) | 1311 if((res = s_mp_sqr(&x)) != MP_OKAY) |
1312 goto CLEANUP; | 1312 goto CLEANUP; |
1313 } | 1313 } |
1314 | 1314 |
1315 if(mp_iseven(b)) | 1315 if(mp_iseven(b)) |
1316 SIGN(&s) = SIGN(a); | 1316 SIGN(&s) = SIGN(a); |
1317 | 1317 |
1318 res = mp_copy(&s, c); | 1318 res = mp_copy(&s, c); |
1319 | 1319 |
1360 if(SIGN(m) == MP_NEG) | 1360 if(SIGN(m) == MP_NEG) |
1361 return MP_RANGE; | 1361 return MP_RANGE; |
1362 | 1362 |
1363 /* | 1363 /* |
1364 If |a| > m, we need to divide to get the remainder and take the | 1364 If |a| > m, we need to divide to get the remainder and take the |
1365 absolute value. | 1365 absolute value. |
1366 | 1366 |
1367 If |a| < m, we don't need to do any division, just copy and adjust | 1367 If |a| < m, we don't need to do any division, just copy and adjust |
1368 the sign (if a is negative). | 1368 the sign (if a is negative). |
1369 | 1369 |
1370 If |a| == m, we can simply set the result to zero. | 1370 If |a| == m, we can simply set the result to zero. |
1374 that |a| != m, so we do those first. | 1374 that |a| != m, so we do those first. |
1375 */ | 1375 */ |
1376 if((mag = s_mp_cmp(a, m)) > 0) { | 1376 if((mag = s_mp_cmp(a, m)) > 0) { |
1377 if((res = mp_div(a, m, NULL, c)) != MP_OKAY) | 1377 if((res = mp_div(a, m, NULL, c)) != MP_OKAY) |
1378 return res; | 1378 return res; |
1379 | 1379 |
1380 if(SIGN(c) == MP_NEG) { | 1380 if(SIGN(c) == MP_NEG) { |
1381 if((res = mp_add(c, m, c)) != MP_OKAY) | 1381 if((res = mp_add(c, m, c)) != MP_OKAY) |
1382 return res; | 1382 return res; |
1383 } | 1383 } |
1384 | 1384 |
1389 if(mp_cmp_z(a) < 0) { | 1389 if(mp_cmp_z(a) < 0) { |
1390 if((res = mp_add(c, m, c)) != MP_OKAY) | 1390 if((res = mp_add(c, m, c)) != MP_OKAY) |
1391 return res; | 1391 return res; |
1392 | 1392 |
1393 } | 1393 } |
1394 | 1394 |
1395 } else { | 1395 } else { |
1396 mp_zero(c); | 1396 mp_zero(c); |
1397 | 1397 |
1398 } | 1398 } |
1399 | 1399 |
1462 /* Cannot take square root of a negative value */ | 1462 /* Cannot take square root of a negative value */ |
1463 if(SIGN(a) == MP_NEG) | 1463 if(SIGN(a) == MP_NEG) |
1464 return MP_RANGE; | 1464 return MP_RANGE; |
1465 | 1465 |
1466 /* Special cases for zero and one, trivial */ | 1466 /* Special cases for zero and one, trivial */ |
1467 if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ) | 1467 if(mp_cmp_d(a, 0) == MP_EQ || mp_cmp_d(a, 1) == MP_EQ) |
1468 return mp_copy(a, b); | 1468 return mp_copy(a, b); |
1469 | 1469 |
1470 /* Initialize the temporaries we'll use below */ | 1470 /* Initialize the temporaries we'll use below */ |
1471 if((res = mp_init_size(&t, USED(a))) != MP_OKAY) | 1471 if((res = mp_init_size(&t, USED(a))) != MP_OKAY) |
1472 return res; | 1472 return res; |
1473 | 1473 |
1474 /* Compute an initial guess for the iteration as a itself */ | 1474 /* Compute an initial guess for the iteration as a itself */ |
1506 s_mp_exch(&x, b); | 1506 s_mp_exch(&x, b); |
1507 | 1507 |
1508 CLEANUP: | 1508 CLEANUP: |
1509 mp_clear(&x); | 1509 mp_clear(&x); |
1510 X: | 1510 X: |
1511 mp_clear(&t); | 1511 mp_clear(&t); |
1512 | 1512 |
1513 return res; | 1513 return res; |
1514 | 1514 |
1515 } /* end mp_sqrt() */ | 1515 } /* end mp_sqrt() */ |
1516 | 1516 |
1624 mp_exptmod(a, b, m, c) | 1624 mp_exptmod(a, b, m, c) |
1625 | 1625 |
1626 Compute c = (a ** b) mod m. Uses a standard square-and-multiply | 1626 Compute c = (a ** b) mod m. Uses a standard square-and-multiply |
1627 method with modular reductions at each step. (This is basically the | 1627 method with modular reductions at each step. (This is basically the |
1628 same code as mp_expt(), except for the addition of the reductions) | 1628 same code as mp_expt(), except for the addition of the reductions) |
1629 | 1629 |
1630 The modular reductions are done using Barrett's algorithm (see | 1630 The modular reductions are done using Barrett's algorithm (see |
1631 s_mp_reduce() below for details) | 1631 s_mp_reduce() below for details) |
1632 */ | 1632 */ |
1633 | 1633 |
1634 mp_err mp_exptmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c) | 1634 mp_err mp_exptmod(mp_int *a, mp_int *b, mp_int *m, mp_int *c) |
1635 { | 1635 { |
1636 mp_int s, x, mu; | 1636 mp_int s, x, mu; |
1637 mp_err res; | 1637 mp_err res; |
1638 mp_digit d, *db = DIGITS(b); | 1638 mp_digit d, *db = DIGITS(b); |
1639 mp_size ub = USED(b); | 1639 mp_size ub = USED(b); |
1640 int dig, bit; | 1640 unsigned int bit, dig; |
1641 | 1641 |
1642 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); | 1642 ARGCHK(a != NULL && b != NULL && c != NULL, MP_BADARG); |
1643 | 1643 |
1644 if(mp_cmp_z(b) < 0 || mp_cmp_z(m) <= 0) | 1644 if(mp_cmp_z(b) < 0 || mp_cmp_z(m) <= 0) |
1645 return MP_RANGE; | 1645 return MP_RANGE; |
1653 goto MU; | 1653 goto MU; |
1654 | 1654 |
1655 mp_set(&s, 1); | 1655 mp_set(&s, 1); |
1656 | 1656 |
1657 /* mu = b^2k / m */ | 1657 /* mu = b^2k / m */ |
1658 s_mp_add_d(&mu, 1); | 1658 s_mp_add_d(&mu, 1); |
1659 s_mp_lshd(&mu, 2 * USED(m)); | 1659 s_mp_lshd(&mu, 2 * USED(m)); |
1660 if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY) | 1660 if((res = mp_div(&mu, m, &mu, NULL)) != MP_OKAY) |
1661 goto CLEANUP; | 1661 goto CLEANUP; |
1662 | 1662 |
1663 /* Loop over digits of b in ascending order, except highest order */ | 1663 /* Loop over digits of b in ascending order, except highest order */ |
1864 { | 1864 { |
1865 mp_int tmp; | 1865 mp_int tmp; |
1866 int out; | 1866 int out; |
1867 | 1867 |
1868 ARGCHK(a != NULL, MP_EQ); | 1868 ARGCHK(a != NULL, MP_EQ); |
1869 | 1869 |
1870 mp_init(&tmp); mp_set_int(&tmp, z); | 1870 mp_init(&tmp); mp_set_int(&tmp, z); |
1871 out = mp_cmp(a, &tmp); | 1871 out = mp_cmp(a, &tmp); |
1872 mp_clear(&tmp); | 1872 mp_clear(&tmp); |
1873 | 1873 |
1874 return out; | 1874 return out; |
1951 | 1951 |
1952 /* Initialize t */ | 1952 /* Initialize t */ |
1953 if(mp_isodd(&u)) { | 1953 if(mp_isodd(&u)) { |
1954 if((res = mp_copy(&v, &t)) != MP_OKAY) | 1954 if((res = mp_copy(&v, &t)) != MP_OKAY) |
1955 goto CLEANUP; | 1955 goto CLEANUP; |
1956 | 1956 |
1957 /* t = -v */ | 1957 /* t = -v */ |
1958 if(SIGN(&v) == MP_ZPOS) | 1958 if(SIGN(&v) == MP_ZPOS) |
1959 SIGN(&t) = MP_NEG; | 1959 SIGN(&t) = MP_NEG; |
1960 else | 1960 else |
1961 SIGN(&t) = MP_ZPOS; | 1961 SIGN(&t) = MP_ZPOS; |
1962 | 1962 |
1963 } else { | 1963 } else { |
1964 if((res = mp_copy(&u, &t)) != MP_OKAY) | 1964 if((res = mp_copy(&u, &t)) != MP_OKAY) |
1965 goto CLEANUP; | 1965 goto CLEANUP; |
1966 | 1966 |
1967 } | 1967 } |
2150 if(x) | 2150 if(x) |
2151 if((res = mp_copy(&C, x)) != MP_OKAY) goto CLEANUP; | 2151 if((res = mp_copy(&C, x)) != MP_OKAY) goto CLEANUP; |
2152 | 2152 |
2153 if(y) | 2153 if(y) |
2154 if((res = mp_copy(&D, y)) != MP_OKAY) goto CLEANUP; | 2154 if((res = mp_copy(&D, y)) != MP_OKAY) goto CLEANUP; |
2155 | 2155 |
2156 if(g) | 2156 if(g) |
2157 if((res = mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP; | 2157 if((res = mp_mul(&gx, &v, g)) != MP_OKAY) goto CLEANUP; |
2158 | 2158 |
2159 break; | 2159 break; |
2160 } | 2160 } |
2253 /*------------------------------------------------------------------------*/ | 2253 /*------------------------------------------------------------------------*/ |
2254 /* {{{ More I/O Functions */ | 2254 /* {{{ More I/O Functions */ |
2255 | 2255 |
2256 /* {{{ mp_read_signed_bin(mp, str, len) */ | 2256 /* {{{ mp_read_signed_bin(mp, str, len) */ |
2257 | 2257 |
2258 /* | 2258 /* |
2259 mp_read_signed_bin(mp, str, len) | 2259 mp_read_signed_bin(mp, str, len) |
2260 | 2260 |
2261 Read in a raw value (base 256) into the given mp_int | 2261 Read in a raw value (base 256) into the given mp_int |
2262 */ | 2262 */ |
2263 | 2263 |
2330 return res; | 2330 return res; |
2331 | 2331 |
2332 if((res = mp_add_d(mp, str[ix], mp)) != MP_OKAY) | 2332 if((res = mp_add_d(mp, str[ix], mp)) != MP_OKAY) |
2333 return res; | 2333 return res; |
2334 } | 2334 } |
2335 | 2335 |
2336 return MP_OKAY; | 2336 return MP_OKAY; |
2337 | 2337 |
2338 } /* end mp_read_unsigned_bin() */ | 2338 } /* end mp_read_unsigned_bin() */ |
2339 | 2339 |
2340 /* }}} */ | 2340 /* }}} */ |
2341 | 2341 |
2342 /* {{{ mp_unsigned_bin_size(mp) */ | 2342 /* {{{ mp_unsigned_bin_size(mp) */ |
2343 | 2343 |
2344 int mp_unsigned_bin_size(mp_int *mp) | 2344 int mp_unsigned_bin_size(mp_int *mp) |
2345 { | 2345 { |
2346 mp_digit topdig; | 2346 mp_digit topdig; |
2347 int count; | 2347 int count; |
2348 | 2348 |
2349 ARGCHK(mp != NULL, 0); | 2349 ARGCHK(mp != NULL, 0); |
2385 return MP_OKAY; | 2385 return MP_OKAY; |
2386 } | 2386 } |
2387 | 2387 |
2388 /* Generate digits in reverse order */ | 2388 /* Generate digits in reverse order */ |
2389 while(dp < end) { | 2389 while(dp < end) { |
2390 int ix; | 2390 unsigned int ix; |
2391 | 2391 |
2392 d = *dp; | 2392 d = *dp; |
2393 for(ix = 0; ix < sizeof(mp_digit); ++ix) { | 2393 for(ix = 0; ix < sizeof(mp_digit); ++ix) { |
2394 *spos = d & UCHAR_MAX; | 2394 *spos = d & UCHAR_MAX; |
2395 d >>= CHAR_BIT; | 2395 d >>= CHAR_BIT; |
2438 ++len; | 2438 ++len; |
2439 d >>= 1; | 2439 d >>= 1; |
2440 } | 2440 } |
2441 | 2441 |
2442 return len; | 2442 return len; |
2443 | 2443 |
2444 } /* end mp_count_bits() */ | 2444 } /* end mp_count_bits() */ |
2445 | 2445 |
2446 /* }}} */ | 2446 /* }}} */ |
2447 | 2447 |
2448 /* {{{ mp_read_radix(mp, str, radix) */ | 2448 /* {{{ mp_read_radix(mp, str, radix) */ |
2460 { | 2460 { |
2461 int ix = 0, val = 0; | 2461 int ix = 0, val = 0; |
2462 mp_err res; | 2462 mp_err res; |
2463 mp_sign sig = MP_ZPOS; | 2463 mp_sign sig = MP_ZPOS; |
2464 | 2464 |
2465 ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX, | 2465 ARGCHK(mp != NULL && str != NULL && radix >= 2 && radix <= MAX_RADIX, |
2466 MP_BADARG); | 2466 MP_BADARG); |
2467 | 2467 |
2468 mp_zero(mp); | 2468 mp_zero(mp); |
2469 | 2469 |
2470 /* Skip leading non-digit characters until a digit or '-' or '+' */ | 2470 /* Skip leading non-digit characters until a digit or '-' or '+' */ |
2471 while(str[ix] && | 2471 while(str[ix] && |
2472 (s_mp_tovalue(str[ix], radix) < 0) && | 2472 (s_mp_tovalue(str[ix], radix) < 0) && |
2473 str[ix] != '-' && | 2473 str[ix] != '-' && |
2474 str[ix] != '+') { | 2474 str[ix] != '+') { |
2475 ++ix; | 2475 ++ix; |
2476 } | 2476 } |
2477 | 2477 |
2523 /* {{{ mp_value_radix_size(num, qty, radix) */ | 2523 /* {{{ mp_value_radix_size(num, qty, radix) */ |
2524 | 2524 |
2525 /* num = number of digits | 2525 /* num = number of digits |
2526 qty = number of bits per digit | 2526 qty = number of bits per digit |
2527 radix = target base | 2527 radix = target base |
2528 | 2528 |
2529 Return the number of digits in the specified radix that would be | 2529 Return the number of digits in the specified radix that would be |
2530 needed to express 'num' digits of 'qty' bits each. | 2530 needed to express 'num' digits of 'qty' bits each. |
2531 */ | 2531 */ |
2532 int mp_value_radix_size(int num, int qty, int radix) | 2532 int mp_value_radix_size(int num, int qty, int radix) |
2533 { | 2533 { |
2539 | 2539 |
2540 /* }}} */ | 2540 /* }}} */ |
2541 | 2541 |
2542 /* {{{ mp_toradix(mp, str, radix) */ | 2542 /* {{{ mp_toradix(mp, str, radix) */ |
2543 | 2543 |
2544 mp_err mp_toradix(mp_int *mp, unsigned char *str, int radix) | 2544 mp_err mp_toradix(mp_int *mp, char *str, int radix) |
2545 { | 2545 { |
2546 int ix, pos = 0; | 2546 int ix, pos = 0; |
2547 | 2547 |
2548 ARGCHK(mp != NULL && str != NULL, MP_BADARG); | 2548 ARGCHK(mp != NULL && str != NULL, MP_BADARG); |
2549 ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE); | 2549 ARGCHK(radix > 1 && radix <= MAX_RADIX, MP_RANGE); |
2585 str[pos--] = '\0'; | 2585 str[pos--] = '\0'; |
2586 | 2586 |
2587 /* Reverse the digits and sign indicator */ | 2587 /* Reverse the digits and sign indicator */ |
2588 ix = 0; | 2588 ix = 0; |
2589 while(ix < pos) { | 2589 while(ix < pos) { |
2590 char tmp = str[ix]; | 2590 char _tmp = str[ix]; |
2591 | 2591 |
2592 str[ix] = str[pos]; | 2592 str[ix] = str[pos]; |
2593 str[pos] = tmp; | 2593 str[pos] = _tmp; |
2594 ++ix; | 2594 ++ix; |
2595 --pos; | 2595 --pos; |
2596 } | 2596 } |
2597 | 2597 |
2598 mp_clear(&tmp); | 2598 mp_clear(&tmp); |
2599 } | 2599 } |
2600 | 2600 |
2601 return MP_OKAY; | 2601 return MP_OKAY; |
2602 | 2602 |
2804 | 2804 |
2805 /* {{{ Arithmetic helpers */ | 2805 /* {{{ Arithmetic helpers */ |
2806 | 2806 |
2807 /* {{{ s_mp_lshd(mp, p) */ | 2807 /* {{{ s_mp_lshd(mp, p) */ |
2808 | 2808 |
2809 /* | 2809 /* |
2810 Shift mp leftward by p digits, growing if needed, and zero-filling | 2810 Shift mp leftward by p digits, growing if needed, and zero-filling |
2811 the in-shifted digits at the right end. This is a convenient | 2811 the in-shifted digits at the right end. This is a convenient |
2812 alternative to multiplication by powers of the radix | 2812 alternative to multiplication by powers of the radix |
2813 */ | 2813 */ |
2814 | 2814 |
2815 mp_err s_mp_lshd(mp_int *mp, mp_size p) | 2815 mp_err s_mp_lshd(mp_int *mp, mp_size p) |
2816 { | 2816 { |
2817 mp_err res; | 2817 mp_err res; |
2818 mp_size pos; | 2818 mp_size pos; |
2819 mp_digit *dp; | 2819 mp_digit *dp; |
2820 int ix; | 2820 int ix; |
2821 | 2821 |
2822 if(p == 0) | 2822 if(p == 0) |
2823 return MP_OKAY; | 2823 return MP_OKAY; |
2824 | 2824 |
2825 if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY) | 2825 if((res = s_mp_pad(mp, USED(mp) + p)) != MP_OKAY) |
2827 | 2827 |
2828 pos = USED(mp) - 1; | 2828 pos = USED(mp) - 1; |
2829 dp = DIGITS(mp); | 2829 dp = DIGITS(mp); |
2830 | 2830 |
2831 /* Shift all the significant figures over as needed */ | 2831 /* Shift all the significant figures over as needed */ |
2832 for(ix = pos - p; ix >= 0; ix--) | 2832 for(ix = pos - p; ix >= 0; ix--) |
2833 dp[ix + p] = dp[ix]; | 2833 dp[ix + p] = dp[ix]; |
2834 | 2834 |
2835 /* Fill the bottom digits with zeroes */ | 2835 /* Fill the bottom digits with zeroes */ |
2836 for(ix = 0; ix < p; ix++) | 2836 for(ix = 0; (unsigned)ix < p; ix++) |
2837 dp[ix] = 0; | 2837 dp[ix] = 0; |
2838 | 2838 |
2839 return MP_OKAY; | 2839 return MP_OKAY; |
2840 | 2840 |
2841 } /* end s_mp_lshd() */ | 2841 } /* end s_mp_lshd() */ |
2842 | 2842 |
2843 /* }}} */ | 2843 /* }}} */ |
2844 | 2844 |
2845 /* {{{ s_mp_rshd(mp, p) */ | 2845 /* {{{ s_mp_rshd(mp, p) */ |
2846 | 2846 |
2847 /* | 2847 /* |
2848 Shift mp rightward by p digits. Maintains the invariant that | 2848 Shift mp rightward by p digits. Maintains the invariant that |
2849 digits above the precision are all zero. Digits shifted off the | 2849 digits above the precision are all zero. Digits shifted off the |
2850 end are lost. Cannot fail. | 2850 end are lost. Cannot fail. |
2851 */ | 2851 */ |
2852 | 2852 |
2896 | 2896 |
2897 /* {{{ s_mp_mul_2(mp) */ | 2897 /* {{{ s_mp_mul_2(mp) */ |
2898 | 2898 |
2899 mp_err s_mp_mul_2(mp_int *mp) | 2899 mp_err s_mp_mul_2(mp_int *mp) |
2900 { | 2900 { |
2901 int ix; | 2901 unsigned int ix; |
2902 mp_digit kin = 0, kout, *dp = DIGITS(mp); | 2902 mp_digit kin = 0, kout, *dp = DIGITS(mp); |
2903 mp_err res; | 2903 mp_err res; |
2904 | 2904 |
2905 /* Shift digits leftward by 1 bit */ | 2905 /* Shift digits leftward by 1 bit */ |
2906 for(ix = 0; ix < USED(mp); ix++) { | 2906 for(ix = 0; ix < USED(mp); ix++) { |
2968 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d) | 2968 mp_err s_mp_mul_2d(mp_int *mp, mp_digit d) |
2969 { | 2969 { |
2970 mp_err res; | 2970 mp_err res; |
2971 mp_digit save, next, mask, *dp; | 2971 mp_digit save, next, mask, *dp; |
2972 mp_size used; | 2972 mp_size used; |
2973 int ix; | 2973 unsigned int ix; |
2974 | 2974 |
2975 if((res = s_mp_lshd(mp, d / DIGIT_BIT)) != MP_OKAY) | 2975 if((res = s_mp_lshd(mp, d / DIGIT_BIT)) != MP_OKAY) |
2976 return res; | 2976 return res; |
2977 | 2977 |
2978 dp = DIGITS(mp); used = USED(mp); | 2978 dp = DIGITS(mp); used = USED(mp); |
3052 accomplish by multiplying a and b by a constant. This constant is | 3052 accomplish by multiplying a and b by a constant. This constant is |
3053 returned (so that it can be divided back out of the remainder at the | 3053 returned (so that it can be divided back out of the remainder at the |
3054 end of the division process). | 3054 end of the division process). |
3055 | 3055 |
3056 We multiply by the smallest power of 2 that gives us a leading digit | 3056 We multiply by the smallest power of 2 that gives us a leading digit |
3057 at least half the radix. By choosing a power of 2, we simplify the | 3057 at least half the radix. By choosing a power of 2, we simplify the |
3058 multiplication and division steps to simple shifts. | 3058 multiplication and division steps to simple shifts. |
3059 */ | 3059 */ |
3060 mp_digit s_mp_norm(mp_int *a, mp_int *b) | 3060 mp_digit s_mp_norm(mp_int *a, mp_int *b) |
3061 { | 3061 { |
3062 mp_digit t, d = 0; | 3062 mp_digit t, d = 0; |
3064 t = DIGIT(b, USED(b) - 1); | 3064 t = DIGIT(b, USED(b) - 1); |
3065 while(t < (RADIX / 2)) { | 3065 while(t < (RADIX / 2)) { |
3066 t <<= 1; | 3066 t <<= 1; |
3067 ++d; | 3067 ++d; |
3068 } | 3068 } |
3069 | 3069 |
3070 if(d != 0) { | 3070 if(d != 0) { |
3071 s_mp_mul_2d(a, d); | 3071 s_mp_mul_2d(a, d); |
3072 s_mp_mul_2d(b, d); | 3072 s_mp_mul_2d(b, d); |
3073 } | 3073 } |
3074 | 3074 |
3186 | 3186 |
3187 /* If there is a precision increase, take care of it here; the above | 3187 /* If there is a precision increase, take care of it here; the above |
3188 test guarantees we have enough storage to do this safely. | 3188 test guarantees we have enough storage to do this safely. |
3189 */ | 3189 */ |
3190 if(k) { | 3190 if(k) { |
3191 dp[max] = k; | 3191 dp[max] = k; |
3192 USED(a) = max + 1; | 3192 USED(a) = max + 1; |
3193 } | 3193 } |
3194 | 3194 |
3195 s_mp_clamp(a); | 3195 s_mp_clamp(a); |
3196 | 3196 |
3197 return MP_OKAY; | 3197 return MP_OKAY; |
3198 | 3198 |
3199 } /* end s_mp_mul_d() */ | 3199 } /* end s_mp_mul_d() */ |
3200 | 3200 |
3201 /* }}} */ | 3201 /* }}} */ |
3202 | 3202 |
3203 /* {{{ s_mp_div_d(mp, d, r) */ | 3203 /* {{{ s_mp_div_d(mp, d, r) */ |
3287 *pa++ = ACCUM(w); | 3287 *pa++ = ACCUM(w); |
3288 w = CARRYOUT(w); | 3288 w = CARRYOUT(w); |
3289 } | 3289 } |
3290 | 3290 |
3291 /* If we run out of 'b' digits before we're actually done, make | 3291 /* If we run out of 'b' digits before we're actually done, make |
3292 sure the carries get propagated upward... | 3292 sure the carries get propagated upward... |
3293 */ | 3293 */ |
3294 used = USED(a); | 3294 used = USED(a); |
3295 while(w && ix < used) { | 3295 while(w && ix < used) { |
3296 w += *pa; | 3296 w += *pa; |
3297 *pa++ = ACCUM(w); | 3297 *pa++ = ACCUM(w); |
3349 } | 3349 } |
3350 | 3350 |
3351 /* Clobber any leading zeroes we created */ | 3351 /* Clobber any leading zeroes we created */ |
3352 s_mp_clamp(a); | 3352 s_mp_clamp(a); |
3353 | 3353 |
3354 /* | 3354 /* |
3355 If there was a borrow out, then |b| > |a| in violation | 3355 If there was a borrow out, then |b| > |a| in violation |
3356 of our input invariant. We've already done the work, | 3356 of our input invariant. We've already done the work, |
3357 but we'll at least complain about it... | 3357 but we'll at least complain about it... |
3358 */ | 3358 */ |
3359 if(w) | 3359 if(w) |
3385 #ifndef SHRT_MUL | 3385 #ifndef SHRT_MUL |
3386 s_mp_mul(&q, m); | 3386 s_mp_mul(&q, m); |
3387 s_mp_mod_2d(&q, (mp_digit)(DIGIT_BIT * (um + 1))); | 3387 s_mp_mod_2d(&q, (mp_digit)(DIGIT_BIT * (um + 1))); |
3388 #else | 3388 #else |
3389 s_mp_mul_dig(&q, m, um + 1); | 3389 s_mp_mul_dig(&q, m, um + 1); |
3390 #endif | 3390 #endif |
3391 | 3391 |
3392 /* x = x - q */ | 3392 /* x = x - q */ |
3393 if((res = mp_sub(x, &q, x)) != MP_OKAY) | 3393 if((res = mp_sub(x, &q, x)) != MP_OKAY) |
3394 goto CLEANUP; | 3394 goto CLEANUP; |
3395 | 3395 |
3439 | 3439 |
3440 /* Outer loop: Digits of b */ | 3440 /* Outer loop: Digits of b */ |
3441 | 3441 |
3442 pb = DIGITS(b); | 3442 pb = DIGITS(b); |
3443 for(ix = 0; ix < ub; ++ix, ++pb) { | 3443 for(ix = 0; ix < ub; ++ix, ++pb) { |
3444 if(*pb == 0) | 3444 if(*pb == 0) |
3445 continue; | 3445 continue; |
3446 | 3446 |
3447 /* Inner product: Digits of a */ | 3447 /* Inner product: Digits of a */ |
3448 pa = DIGITS(a); | 3448 pa = DIGITS(a); |
3449 for(jx = 0; jx < ua; ++jx, ++pa) { | 3449 for(jx = 0; jx < ua; ++jx, ++pa) { |
3478 mp_digit *pa, *pt; | 3478 mp_digit *pa, *pt; |
3479 | 3479 |
3480 for(ix = 0; ix < len; ++ix, ++b) { | 3480 for(ix = 0; ix < len; ++ix, ++b) { |
3481 if(*b == 0) | 3481 if(*b == 0) |
3482 continue; | 3482 continue; |
3483 | 3483 |
3484 pa = a; | 3484 pa = a; |
3485 for(jx = 0; jx < len; ++jx, ++pa) { | 3485 for(jx = 0; jx < len; ++jx, ++pa) { |
3486 pt = out + ix + jx; | 3486 pt = out + ix + jx; |
3487 w = *b * *pa + k + *pt; | 3487 w = *b * *pa + k + *pt; |
3488 *pt = ACCUM(w); | 3488 *pt = ACCUM(w); |
3545 overflow, we have to check explicitly for overflow conditions | 3545 overflow, we have to check explicitly for overflow conditions |
3546 before they happen. | 3546 before they happen. |
3547 */ | 3547 */ |
3548 for(jx = ix + 1, pa2 = DIGITS(a) + jx; jx < used; ++jx, ++pa2) { | 3548 for(jx = ix + 1, pa2 = DIGITS(a) + jx; jx < used; ++jx, ++pa2) { |
3549 mp_word u = 0, v; | 3549 mp_word u = 0, v; |
3550 | 3550 |
3551 /* Store this in a temporary to avoid indirections later */ | 3551 /* Store this in a temporary to avoid indirections later */ |
3552 pt = pbt + ix + jx; | 3552 pt = pbt + ix + jx; |
3553 | 3553 |
3554 /* Compute the multiplicative step */ | 3554 /* Compute the multiplicative step */ |
3555 w = *pa1 * *pa2; | 3555 w = *pa1 * *pa2; |
3566 | 3566 |
3567 /* Compute the additive step */ | 3567 /* Compute the additive step */ |
3568 v = *pt + k; | 3568 v = *pt + k; |
3569 | 3569 |
3570 /* If we do not already have an overflow carry, check to see | 3570 /* If we do not already have an overflow carry, check to see |
3571 if the addition will cause one, and set the carry out if so | 3571 if the addition will cause one, and set the carry out if so |
3572 */ | 3572 */ |
3573 u |= ((MP_WORD_MAX - v) < w); | 3573 u |= ((MP_WORD_MAX - v) < w); |
3574 | 3574 |
3575 /* Add in the rest, again ignoring overflow */ | 3575 /* Add in the rest, again ignoring overflow */ |
3576 w += v; | 3576 w += v; |
3590 k = CARRYOUT(k); | 3590 k = CARRYOUT(k); |
3591 | 3591 |
3592 /* If we are carrying out, propagate the carry to the next digit | 3592 /* If we are carrying out, propagate the carry to the next digit |
3593 in the output. This may cascade, so we have to be somewhat | 3593 in the output. This may cascade, so we have to be somewhat |
3594 circumspect -- but we will have enough precision in the output | 3594 circumspect -- but we will have enough precision in the output |
3595 that we won't overflow | 3595 that we won't overflow |
3596 */ | 3596 */ |
3597 kx = 1; | 3597 kx = 1; |
3598 while(k) { | 3598 while(k) { |
3599 k = pbt[ix + jx + kx] + 1; | 3599 k = pbt[ix + jx + kx] + 1; |
3600 pbt[ix + jx + kx] = ACCUM(k); | 3600 pbt[ix + jx + kx] = ACCUM(k); |
3662 ix = USED(a) - 1; | 3662 ix = USED(a) - 1; |
3663 | 3663 |
3664 while(ix >= 0) { | 3664 while(ix >= 0) { |
3665 /* Find a partial substring of a which is at least b */ | 3665 /* Find a partial substring of a which is at least b */ |
3666 while(s_mp_cmp(&rem, b) < 0 && ix >= 0) { | 3666 while(s_mp_cmp(&rem, b) < 0 && ix >= 0) { |
3667 if((res = s_mp_lshd(&rem, 1)) != MP_OKAY) | 3667 if((res = s_mp_lshd(&rem, 1)) != MP_OKAY) |
3668 goto CLEANUP; | 3668 goto CLEANUP; |
3669 | 3669 |
3670 if((res = s_mp_lshd(", 1)) != MP_OKAY) | 3670 if((res = s_mp_lshd(", 1)) != MP_OKAY) |
3671 goto CLEANUP; | 3671 goto CLEANUP; |
3672 | 3672 |
3674 s_mp_clamp(&rem); | 3674 s_mp_clamp(&rem); |
3675 --ix; | 3675 --ix; |
3676 } | 3676 } |
3677 | 3677 |
3678 /* If we didn't find one, we're finished dividing */ | 3678 /* If we didn't find one, we're finished dividing */ |
3679 if(s_mp_cmp(&rem, b) < 0) | 3679 if(s_mp_cmp(&rem, b) < 0) |
3680 break; | 3680 break; |
3681 | 3681 |
3682 /* Compute a guess for the next quotient digit */ | 3682 /* Compute a guess for the next quotient digit */ |
3683 q = DIGIT(&rem, USED(&rem) - 1); | 3683 q = DIGIT(&rem, USED(&rem) - 1); |
3684 if(q <= DIGIT(b, USED(b) - 1) && USED(&rem) > 1) | 3684 if(q <= DIGIT(b, USED(b) - 1) && USED(&rem) > 1) |
3685 q = (q << DIGIT_BIT) | DIGIT(&rem, USED(&rem) - 2); | 3685 q = (q << DIGIT_BIT) | DIGIT(&rem, USED(&rem) - 2); |
3693 /* See what that multiplies out to */ | 3693 /* See what that multiplies out to */ |
3694 mp_copy(b, &t); | 3694 mp_copy(b, &t); |
3695 if((res = s_mp_mul_d(&t, q)) != MP_OKAY) | 3695 if((res = s_mp_mul_d(&t, q)) != MP_OKAY) |
3696 goto CLEANUP; | 3696 goto CLEANUP; |
3697 | 3697 |
3698 /* | 3698 /* |
3699 If it's too big, back it off. We should not have to do this | 3699 If it's too big, back it off. We should not have to do this |
3700 more than once, or, in rare cases, twice. Knuth describes a | 3700 more than once, or, in rare cases, twice. Knuth describes a |
3701 method by which this could be reduced to a maximum of once, but | 3701 method by which this could be reduced to a maximum of once, but |
3702 I didn't implement that here. | 3702 I didn't implement that here. |
3703 */ | 3703 */ |
3717 */ | 3717 */ |
3718 DIGIT(", 0) = q; | 3718 DIGIT(", 0) = q; |
3719 } | 3719 } |
3720 | 3720 |
3721 /* Denormalize remainder */ | 3721 /* Denormalize remainder */ |
3722 if(d != 0) | 3722 if(d != 0) |
3723 s_mp_div_2d(&rem, d); | 3723 s_mp_div_2d(&rem, d); |
3724 | 3724 |
3725 s_mp_clamp("); | 3725 s_mp_clamp("); |
3726 s_mp_clamp(&rem); | 3726 s_mp_clamp(&rem); |
3727 | 3727 |
3728 /* Copy quotient back to output */ | 3728 /* Copy quotient back to output */ |
3729 s_mp_exch(", a); | 3729 s_mp_exch(", a); |
3730 | 3730 |
3731 /* Copy remainder back to output */ | 3731 /* Copy remainder back to output */ |
3732 s_mp_exch(&rem, b); | 3732 s_mp_exch(&rem, b); |
3733 | 3733 |
3734 CLEANUP: | 3734 CLEANUP: |
3735 mp_clear(&rem); | 3735 mp_clear(&rem); |
3755 bit = k % DIGIT_BIT; | 3755 bit = k % DIGIT_BIT; |
3756 | 3756 |
3757 mp_zero(a); | 3757 mp_zero(a); |
3758 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY) | 3758 if((res = s_mp_pad(a, dig + 1)) != MP_OKAY) |
3759 return res; | 3759 return res; |
3760 | 3760 |
3761 DIGIT(a, dig) |= (1 << bit); | 3761 DIGIT(a, dig) |= (1 << bit); |
3762 | 3762 |
3763 return MP_OKAY; | 3763 return MP_OKAY; |
3764 | 3764 |
3765 } /* end s_mp_2expt() */ | 3765 } /* end s_mp_2expt() */ |
3813 mp_digit *ap = DIGITS(a); | 3813 mp_digit *ap = DIGITS(a); |
3814 | 3814 |
3815 if(ua > 1) | 3815 if(ua > 1) |
3816 return MP_GT; | 3816 return MP_GT; |
3817 | 3817 |
3818 if(*ap < d) | 3818 if(*ap < d) |
3819 return MP_LT; | 3819 return MP_LT; |
3820 else if(*ap > d) | 3820 else if(*ap > d) |
3821 return MP_GT; | 3821 return MP_GT; |
3822 else | 3822 else |
3823 return MP_EQ; | 3823 return MP_EQ; |
3855 | 3855 |
3856 --dp; --ix; | 3856 --dp; --ix; |
3857 } | 3857 } |
3858 | 3858 |
3859 return ((uv - 1) * DIGIT_BIT) + extra; | 3859 return ((uv - 1) * DIGIT_BIT) + extra; |
3860 } | 3860 } |
3861 | 3861 |
3862 return -1; | 3862 return -1; |
3863 | 3863 |
3864 } /* end s_mp_ispow2() */ | 3864 } /* end s_mp_ispow2() */ |
3865 | 3865 |
3899 expected to know what you're up to. | 3899 expected to know what you're up to. |
3900 */ | 3900 */ |
3901 int s_mp_tovalue(char ch, int r) | 3901 int s_mp_tovalue(char ch, int r) |
3902 { | 3902 { |
3903 int val, xch; | 3903 int val, xch; |
3904 | 3904 |
3905 if(r > 36) | 3905 if(r > 36) |
3906 xch = ch; | 3906 xch = ch; |
3907 else | 3907 else |
3908 xch = toupper(ch); | 3908 xch = toupper(ch); |
3909 | 3909 |
3915 val = xch - 'a' + 36; | 3915 val = xch - 'a' + 36; |
3916 else if(xch == '+') | 3916 else if(xch == '+') |
3917 val = 62; | 3917 val = 62; |
3918 else if(xch == '/') | 3918 else if(xch == '/') |
3919 val = 63; | 3919 val = 63; |
3920 else | 3920 else |
3921 return -1; | 3921 return -1; |
3922 | 3922 |
3923 if(val < 0 || val >= r) | 3923 if(val < 0 || val >= r) |
3924 return -1; | 3924 return -1; |
3925 | 3925 |
3937 the value in the given radix. | 3937 the value in the given radix. |
3938 | 3938 |
3939 The results may be odd if you use a radix < 2 or > 64, you are | 3939 The results may be odd if you use a radix < 2 or > 64, you are |
3940 expected to know what you're doing. | 3940 expected to know what you're doing. |
3941 */ | 3941 */ |
3942 | 3942 |
3943 char s_mp_todigit(int val, int r, int low) | 3943 char s_mp_todigit(int val, int r, int low) |
3944 { | 3944 { |
3945 char ch; | 3945 char ch; |
3946 | 3946 |
3947 if(val < 0 || val >= r) | 3947 if(val < 0 || val >= r) |
3958 | 3958 |
3959 /* }}} */ | 3959 /* }}} */ |
3960 | 3960 |
3961 /* {{{ s_mp_outlen(bits, radix) */ | 3961 /* {{{ s_mp_outlen(bits, radix) */ |
3962 | 3962 |
3963 /* | 3963 /* |
3964 Return an estimate for how long a string is needed to hold a radix | 3964 Return an estimate for how long a string is needed to hold a radix |
3965 r representation of a number with 'bits' significant bits. | 3965 r representation of a number with 'bits' significant bits. |
3966 | 3966 |
3967 Does not include space for a sign or a NUL terminator. | 3967 Does not include space for a sign or a NUL terminator. |
3968 */ | 3968 */ |
3978 | 3978 |
3979 /*------------------------------------------------------------------------*/ | 3979 /*------------------------------------------------------------------------*/ |
3980 /* HERE THERE BE DRAGONS */ | 3980 /* HERE THERE BE DRAGONS */ |
3981 /* crc==4242132123, version==2, Sat Feb 02 06:43:52 2002 */ | 3981 /* crc==4242132123, version==2, Sat Feb 02 06:43:52 2002 */ |
3982 | 3982 |
3983 /* $Source: /cvs/libtom/libtommath/mtest/mpi.c,v $ */ | 3983 /* $Source$ */ |
3984 /* $Revision: 1.2 $ */ | 3984 /* $Revision$ */ |
3985 /* $Date: 2005/05/05 14:38:47 $ */ | 3985 /* $Date$ */ |