3
Purpose: Arbitrary precision integer arithmetic routines.
4
Author: M. J. Fromberger <http://spinning-yarns.org/michael/>
5
Info: $Id: imath.c 645 2008-08-03 04:00:30Z sting $
7
Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved.
9
Permission is hereby granted, free of charge, to any person
10
obtaining a copy of this software and associated documentation files
11
(the "Software"), to deal in the Software without restriction,
12
including without limitation the rights to use, copy, modify, merge,
13
publish, distribute, sublicense, and/or sell copies of the Software,
14
and to permit persons to whom the Software is furnished to do so,
15
subject to the following conditions:
17
The above copyright notice and this permission notice shall be
18
included in all copies or substantial portions of the Software.
20
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
21
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
22
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
23
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
24
BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
25
ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
26
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
48
const mp_result MP_OK = 0; /* no error, all is well */
49
const mp_result MP_FALSE = 0; /* boolean false */
50
const mp_result MP_TRUE = -1; /* boolean true */
51
const mp_result MP_MEMORY = -2; /* out of memory */
52
const mp_result MP_RANGE = -3; /* argument out of range */
53
const mp_result MP_UNDEF = -4; /* result undefined */
54
const mp_result MP_TRUNC = -5; /* output truncated */
55
const mp_result MP_BADARG = -6; /* invalid null argument */
56
const mp_result MP_MINERR = -6;
58
const mp_sign MP_NEG = 1; /* value is strictly negative */
59
const mp_sign MP_ZPOS = 0; /* value is non-negative */
61
static const char *s_unknown_err = "unknown result code";
62
static const char *s_error_msg[] = {
66
"argument out of range",
75
/* Argument checking macros
76
Use CHECK() where a return value is required; NRCHECK() elsewhere */
77
#define CHECK(TEST) assert(TEST)
78
#define NRCHECK(TEST) assert(TEST)
80
/* {{{ Logarithm table for computing output sizes */
82
/* The ith entry of this table gives the value of log_i(2).
84
An integer value n requires ceil(log_i(n)) digits to be represented
85
in base i. Since it is easy to compute lg(n), by counting bits, we
86
can compute log_i(n) = lg(n) * log_i(2).
88
The use of this table eliminates a dependency upon linkage against
89
the standard math libraries.
91
static const double s_log2[] = {
92
0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
93
0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
94
0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
95
0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */
96
0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */
97
0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */
98
0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */
99
0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */
100
0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */
101
0.193426404, /* 36 */
105
/* {{{ Various macros */
107
/* Return the number of digits needed to represent a static value */
108
#define MP_VALUE_DIGITS(V) \
109
((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
111
/* Round precision P to nearest word boundary */
112
#define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
114
/* Set array P of S digits to zero */
116
do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P);memset(p__,0,i__);}while(0)
118
/* Copy S digits from array P to array Q */
119
#define COPY(P, Q, S) \
120
do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P),*q__=(Q);\
121
memcpy(q__,p__,i__);}while(0)
123
/* Reverse N elements of type T in array A */
124
#define REV(T, A, N) \
125
do{T *u_=(A),*v_=u_+(N)-1;while(u_<v_){T xch=*u_;*u_++=*v_;*v_--=xch;}}while(0)
128
do{mp_int z_=(Z);mp_size uz_=MP_USED(z_);mp_digit *dz_=MP_DIGITS(z_)+uz_-1;\
129
while(uz_ > 1 && (*dz_-- == 0)) --uz_;MP_USED(z_)=uz_;}while(0)
131
/* Select min/max. Do not provide expressions for which multiple
132
evaluation would be problematic, e.g. x++ */
133
#define MIN(A, B) ((B)<(A)?(B):(A))
134
#define MAX(A, B) ((B)>(A)?(B):(A))
136
/* Exchange lvalues A and B of type T, e.g.
137
SWAP(int, x, y) where x and y are variables of type int. */
138
#define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
140
/* Used to set up and access simple temp stacks within functions. */
141
#define TEMP(K) (temp + (K))
142
#define SETUP(E, C) \
143
do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
145
/* Compare value to zero. */
147
(((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
149
/* Multiply X by Y into Z, ignoring signs. Requires that Z have
150
enough storage preallocated to hold the result. */
151
#define UMUL(X, Y, Z) \
152
do{mp_size ua_=MP_USED(X),ub_=MP_USED(Y);mp_size o_=ua_+ub_;\
153
ZERO(MP_DIGITS(Z),o_);\
154
(void) s_kmul(MP_DIGITS(X),MP_DIGITS(Y),MP_DIGITS(Z),ua_,ub_);\
155
MP_USED(Z)=o_;CLAMP(Z);}while(0)
157
/* Square X into Z. Requires that Z have enough storage to hold the
160
do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
161
(void) s_ksqr(MP_DIGITS(X),MP_DIGITS(Z),ua_);MP_USED(Z)=o_;CLAMP(Z);}while(0)
163
#define UPPER_HALF(W) ((mp_word)((W) >> MP_DIGIT_BIT))
164
#define LOWER_HALF(W) ((mp_digit)(W))
165
#define HIGH_BIT_SET(W) ((W) >> (MP_WORD_BIT - 1))
166
#define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
169
/* {{{ Default configuration settings */
171
/* Default number of digits allocated to a new mp_int */
173
mp_size default_precision = MP_DEFAULT_PREC;
175
static const mp_size default_precision = MP_DEFAULT_PREC;
178
/* Minimum number of digits to invoke recursive multiply */
180
mp_size multiply_threshold = MP_MULT_THRESH;
182
static const mp_size multiply_threshold = MP_MULT_THRESH;
187
/* Allocate a buffer of (at least) num digits, or return
188
NULL if that couldn't be done. */
189
static mp_digit *s_alloc(mp_size num);
191
/* Release a buffer of digits allocated by s_alloc(). */
192
static void s_free(void *ptr);
194
/* Insure that z has at least min digits allocated, resizing if
195
necessary. Returns true if successful, false if out of memory. */
196
static int s_pad(mp_int z, mp_size min);
198
/* Fill in a "fake" mp_int on the stack with a given value */
199
static void s_fake(mp_int z, mp_small value, mp_digit vbuf[]);
201
/* Compare two runs of digits of given length, returns <0, 0, >0 */
202
static int s_cdig(mp_digit *da, mp_digit *db, mp_size len);
204
/* Pack the unsigned digits of v into array t */
205
static int s_vpack(mp_small v, mp_digit t[]);
207
/* Compare magnitudes of a and b, returns <0, 0, >0 */
208
static int s_ucmp(mp_int a, mp_int b);
210
/* Compare magnitudes of a and v, returns <0, 0, >0 */
211
static int s_vcmp(mp_int a, mp_small v);
213
/* Unsigned magnitude addition; assumes dc is big enough.
214
Carry out is returned (no memory allocated). */
215
static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
216
mp_size size_a, mp_size size_b);
218
/* Unsigned magnitude subtraction. Assumes dc is big enough. */
219
static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
220
mp_size size_a, mp_size size_b);
222
/* Unsigned recursive multiplication. Assumes dc is big enough. */
223
static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
224
mp_size size_a, mp_size size_b);
226
/* Unsigned magnitude multiplication. Assumes dc is big enough. */
227
static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
228
mp_size size_a, mp_size size_b);
230
/* Unsigned recursive squaring. Assumes dc is big enough. */
231
static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a);
233
/* Unsigned magnitude squaring. Assumes dc is big enough. */
234
static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a);
236
/* Single digit addition. Assumes a is big enough. */
237
static void s_dadd(mp_int a, mp_digit b);
239
/* Single digit multiplication. Assumes a is big enough. */
240
static void s_dmul(mp_int a, mp_digit b);
242
/* Single digit multiplication on buffers; assumes dc is big enough. */
243
static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc,
246
/* Single digit division. Replaces a with the quotient,
247
returns the remainder. */
248
static mp_digit s_ddiv(mp_int a, mp_digit b);
250
/* Quick division by a power of 2, replaces z (no allocation) */
251
static void s_qdiv(mp_int z, mp_size p2);
253
/* Quick remainder by a power of 2, replaces z (no allocation) */
254
static void s_qmod(mp_int z, mp_size p2);
256
/* Quick multiplication by a power of 2, replaces z.
257
Allocates if necessary; returns false in case this fails. */
258
static int s_qmul(mp_int z, mp_size p2);
260
/* Quick subtraction from a power of 2, replaces z.
261
Allocates if necessary; returns false in case this fails. */
262
static int s_qsub(mp_int z, mp_size p2);
264
/* Return maximum k such that 2^k divides z. */
265
static int s_dp2k(mp_int z);
267
/* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
268
static int s_isp2(mp_int z);
270
/* Set z to 2^k. May allocate; returns false in case this fails. */
271
static int s_2expt(mp_int z, mp_small k);
273
/* Normalize a and b for division, returns normalization constant */
274
static int s_norm(mp_int a, mp_int b);
276
/* Compute constant mu for Barrett reduction, given modulus m, result
277
replaces z, m is untouched. */
278
static mp_result s_brmu(mp_int z, mp_int m);
280
/* Reduce a modulo m, using Barrett's algorithm. */
281
static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2);
283
/* Modular exponentiation, using Barrett reduction */
284
static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c);
286
/* Unsigned magnitude division. Assumes |a| > |b|. Allocates
287
temporaries; overwrites a with quotient, b with remainder. */
288
static mp_result s_udiv(mp_int a, mp_int b);
290
/* Compute the number of digits in radix r required to represent the
291
given value. Does not account for sign flags, terminators, etc. */
292
static int s_outlen(mp_int z, mp_size r);
294
/* Guess how many digits of precision will be needed to represent a
295
radix r value of the specified number of digits. Returns a value
296
guaranteed to be no smaller than the actual number required. */
297
static mp_size s_inlen(int len, mp_size r);
299
/* Convert a character to a digit value in radix r, or
300
-1 if out of range */
301
static int s_ch2val(char c, int r);
303
/* Convert a digit value to a character */
304
static char s_val2ch(int v, int caps);
306
/* Take 2's complement of a buffer in place */
307
static void s_2comp(unsigned char *buf, int len);
309
/* Convert a value to binary, ignoring sign. On input, *limpos is the
310
bound on how many bytes should be written to buf; on output, *limpos
311
is set to the number of bytes actually written. */
312
static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad);
315
/* Dump a representation of the mp_int to standard output */
316
void s_print(char *tag, mp_int z);
317
void s_print_buf(char *tag, mp_digit *buf, mp_size num);
320
/* {{{ mp_int_init(z) */
322
mp_result mp_int_init(mp_int z)
328
z->digits = &(z->single);
338
/* {{{ mp_int_alloc() */
340
mp_int mp_int_alloc(void)
342
mp_int out = malloc(sizeof(mpz_t));
352
/* {{{ mp_int_init_size(z, prec) */
354
mp_result mp_int_init_size(mp_int z, mp_size prec)
359
prec = default_precision;
361
return mp_int_init(z);
363
prec = (mp_size) ROUND_PREC(prec);
365
if((MP_DIGITS(z) = s_alloc(prec)) == NULL)
371
MP_SIGN(z) = MP_ZPOS;
378
/* {{{ mp_int_init_copy(z, old) */
380
mp_result mp_int_init_copy(mp_int z, mp_int old)
385
CHECK(z != NULL && old != NULL);
392
mp_size target = MAX(uold, default_precision);
394
if((res = mp_int_init_size(z, target)) != MP_OK)
399
MP_SIGN(z) = MP_SIGN(old);
400
COPY(MP_DIGITS(old), MP_DIGITS(z), uold);
407
/* {{{ mp_int_init_value(z, value) */
409
mp_result mp_int_init_value(mp_int z, mp_small value)
412
mp_digit vbuf[MP_VALUE_DIGITS(value)];
414
s_fake(&vtmp, value, vbuf);
415
return mp_int_init_copy(z, &vtmp);
420
/* {{{ mp_int_set_value(z, value) */
422
mp_result mp_int_set_value(mp_int z, mp_small value)
425
mp_digit vbuf[MP_VALUE_DIGITS(value)];
427
s_fake(&vtmp, value, vbuf);
428
return mp_int_copy(&vtmp, z);
433
/* {{{ mp_int_clear(z) */
435
void mp_int_clear(mp_int z)
440
if(MP_DIGITS(z) != NULL) {
441
if((void *) MP_DIGITS(z) != (void *) z)
442
s_free(MP_DIGITS(z));
450
/* {{{ mp_int_free(z) */
452
void mp_int_free(mp_int z)
457
free(z); /* note: NOT s_free() */
462
/* {{{ mp_int_copy(a, c) */
464
mp_result mp_int_copy(mp_int a, mp_int c)
466
CHECK(a != NULL && c != NULL);
469
mp_size ua = MP_USED(a);
475
da = MP_DIGITS(a); dc = MP_DIGITS(c);
479
MP_SIGN(c) = MP_SIGN(a);
487
/* {{{ mp_int_swap(a, c) */
489
void mp_int_swap(mp_int a, mp_int c)
501
/* {{{ mp_int_zero(z) */
503
void mp_int_zero(mp_int z)
509
MP_SIGN(z) = MP_ZPOS;
514
/* {{{ mp_int_abs(a, c) */
516
mp_result mp_int_abs(mp_int a, mp_int c)
520
CHECK(a != NULL && c != NULL);
522
if((res = mp_int_copy(a, c)) != MP_OK)
525
MP_SIGN(c) = MP_ZPOS;
531
/* {{{ mp_int_neg(a, c) */
533
mp_result mp_int_neg(mp_int a, mp_int c)
537
CHECK(a != NULL && c != NULL);
539
if((res = mp_int_copy(a, c)) != MP_OK)
543
MP_SIGN(c) = 1 - MP_SIGN(a);
550
/* {{{ mp_int_add(a, b, c) */
552
mp_result mp_int_add(mp_int a, mp_int b, mp_int c)
554
mp_size ua, ub, uc, max;
556
CHECK(a != NULL && b != NULL && c != NULL);
558
ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
561
if(MP_SIGN(a) == MP_SIGN(b)) {
562
/* Same sign -- add magnitudes, preserve sign of addends */
568
carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
572
if(!s_pad(c, max + 1))
575
c->digits[max] = carry;
580
MP_SIGN(c) = MP_SIGN(a);
584
/* Different signs -- subtract magnitudes, preserve sign of greater */
586
int cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */
588
/* Set x to max(a, b), y to min(a, b) to simplify later code.
589
A special case yields zero for equal magnitudes.
602
if(!s_pad(c, MP_USED(x)))
605
/* Subtract smaller from larger */
606
s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
607
MP_USED(c) = MP_USED(x);
610
/* Give result the sign of the larger */
611
MP_SIGN(c) = MP_SIGN(x);
619
/* {{{ mp_int_add_value(a, value, c) */
621
mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c)
624
mp_digit vbuf[MP_VALUE_DIGITS(value)];
626
s_fake(&vtmp, value, vbuf);
628
return mp_int_add(a, &vtmp, c);
633
/* {{{ mp_int_sub(a, b, c) */
635
mp_result mp_int_sub(mp_int a, mp_int b, mp_int c)
637
mp_size ua, ub, uc, max;
639
CHECK(a != NULL && b != NULL && c != NULL);
641
ua = MP_USED(a); ub = MP_USED(b); uc = MP_USED(c);
644
if(MP_SIGN(a) != MP_SIGN(b)) {
645
/* Different signs -- add magnitudes and keep sign of a */
651
carry = s_uadd(MP_DIGITS(a), MP_DIGITS(b), MP_DIGITS(c), ua, ub);
655
if(!s_pad(c, max + 1))
658
c->digits[max] = carry;
663
MP_SIGN(c) = MP_SIGN(a);
667
/* Same signs -- subtract magnitudes */
670
int cmp = s_ucmp(a, b);
676
x = a; y = b; osign = MP_ZPOS;
679
x = b; y = a; osign = MP_NEG;
682
if(MP_SIGN(a) == MP_NEG && cmp != 0)
685
s_usub(MP_DIGITS(x), MP_DIGITS(y), MP_DIGITS(c), MP_USED(x), MP_USED(y));
686
MP_USED(c) = MP_USED(x);
697
/* {{{ mp_int_sub_value(a, value, c) */
699
mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c)
702
mp_digit vbuf[MP_VALUE_DIGITS(value)];
704
s_fake(&vtmp, value, vbuf);
706
return mp_int_sub(a, &vtmp, c);
711
/* {{{ mp_int_mul(a, b, c) */
713
mp_result mp_int_mul(mp_int a, mp_int b, mp_int c)
716
mp_size osize, ua, ub, p = 0;
719
CHECK(a != NULL && b != NULL && c != NULL);
721
/* If either input is zero, we can shortcut multiplication */
722
if(mp_int_compare_zero(a) == 0 || mp_int_compare_zero(b) == 0) {
727
/* Output is positive if inputs have same sign, otherwise negative */
728
osign = (MP_SIGN(a) == MP_SIGN(b)) ? MP_ZPOS : MP_NEG;
730
/* If the output is not identical to any of the inputs, we'll write
731
the results directly; otherwise, allocate a temporary space. */
732
ua = MP_USED(a); ub = MP_USED(b);
734
osize = 4 * ((osize + 1) / 2);
736
if(c == a || c == b) {
737
p = ROUND_PREC(osize);
738
p = MAX(p, default_precision);
740
if((out = s_alloc(p)) == NULL)
751
if(!s_kmul(MP_DIGITS(a), MP_DIGITS(b), out, ua, ub))
754
/* If we allocated a new buffer, get rid of whatever memory c was
755
already using, and fix up its fields to reflect that.
757
if(out != MP_DIGITS(c)) {
758
if((void *) MP_DIGITS(c) != (void *) c)
759
s_free(MP_DIGITS(c));
764
MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
765
CLAMP(c); /* ... right here */
773
/* {{{ mp_int_mul_value(a, value, c) */
775
mp_result mp_int_mul_value(mp_int a, mp_small value, mp_int c)
778
mp_digit vbuf[MP_VALUE_DIGITS(value)];
780
s_fake(&vtmp, value, vbuf);
782
return mp_int_mul(a, &vtmp, c);
787
/* {{{ mp_int_mul_pow2(a, p2, c) */
789
mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c)
792
CHECK(a != NULL && c != NULL && p2 >= 0);
794
if((res = mp_int_copy(a, c)) != MP_OK)
797
if(s_qmul(c, (mp_size) p2))
805
/* {{{ mp_int_sqr(a, c) */
807
mp_result mp_int_sqr(mp_int a, mp_int c)
810
mp_size osize, p = 0;
812
CHECK(a != NULL && c != NULL);
814
/* Get a temporary buffer big enough to hold the result */
815
osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2);
817
p = ROUND_PREC(osize);
818
p = MAX(p, default_precision);
820
if((out = s_alloc(p)) == NULL)
831
s_ksqr(MP_DIGITS(a), out, MP_USED(a));
833
/* Get rid of whatever memory c was already using, and fix up its
834
fields to reflect the new digit array it's using
836
if(out != MP_DIGITS(c)) {
837
if((void *) MP_DIGITS(c) != (void *) c)
838
s_free(MP_DIGITS(c));
843
MP_USED(c) = osize; /* might not be true, but we'll fix it ... */
844
CLAMP(c); /* ... right here */
845
MP_SIGN(c) = MP_ZPOS;
852
/* {{{ mp_int_div(a, b, q, r) */
854
mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r)
856
int cmp, last = 0, lg;
857
mp_result res = MP_OK;
860
mp_sign sa = MP_SIGN(a), sb = MP_SIGN(b);
862
CHECK(a != NULL && b != NULL && q != r);
866
else if((cmp = s_ucmp(a, b)) < 0) {
867
/* If |a| < |b|, no division is required:
870
if(r && (res = mp_int_copy(a, r)) != MP_OK)
879
/* If |a| = |b|, no division is required:
896
/* When |a| > |b|, real division is required. We need someplace to
897
store quotient and remainder, but q and r are allowed to be NULL
898
or to overlap with the inputs.
900
if((lg = s_isp2(b)) < 0) {
902
if((res = mp_int_copy(a, q)) != MP_OK)
909
SETUP(mp_int_init_copy(TEMP(last), a), last);
913
if((res = mp_int_copy(b, r)) != MP_OK)
920
SETUP(mp_int_init_copy(TEMP(last), b), last);
923
if((res = s_udiv(qout, rout)) != MP_OK) goto CLEANUP;
926
if(q && (res = mp_int_copy(a, q)) != MP_OK) goto CLEANUP;
927
if(r && (res = mp_int_copy(a, r)) != MP_OK) goto CLEANUP;
929
if(q) s_qdiv(q, (mp_size) lg); qout = q;
930
if(r) s_qmod(r, (mp_size) lg); rout = r;
933
/* Recompute signs for output */
937
MP_SIGN(rout) = MP_ZPOS;
940
MP_SIGN(qout) = (sa == sb) ? MP_ZPOS : MP_NEG;
942
MP_SIGN(qout) = MP_ZPOS;
945
if(q && (res = mp_int_copy(qout, q)) != MP_OK) goto CLEANUP;
946
if(r && (res = mp_int_copy(rout, r)) != MP_OK) goto CLEANUP;
950
mp_int_clear(TEMP(last));
957
/* {{{ mp_int_mod(a, m, c) */
959
mp_result mp_int_mod(mp_int a, mp_int m, mp_int c)
973
if((res = mp_int_div(a, m, NULL, out)) != MP_OK)
977
res = mp_int_add(out, m, c);
979
res = mp_int_copy(out, c);
990
/* {{{ mp_int_div_value(a, value, q, r) */
992
mp_result mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r)
995
mp_digit vbuf[MP_VALUE_DIGITS(value)];
999
s_fake(&vtmp, value, vbuf);
1001
if((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK)
1005
(void) mp_int_to_int(&rtmp, r); /* can't fail */
1008
mp_int_clear(&rtmp);
1014
/* {{{ mp_int_div_pow2(a, p2, q, r) */
1016
mp_result mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r)
1018
mp_result res = MP_OK;
1020
CHECK(a != NULL && p2 >= 0 && q != r);
1022
if(q != NULL && (res = mp_int_copy(a, q)) == MP_OK)
1023
s_qdiv(q, (mp_size) p2);
1025
if(res == MP_OK && r != NULL && (res = mp_int_copy(a, r)) == MP_OK)
1026
s_qmod(r, (mp_size) p2);
1033
/* {{{ mp_int_expt(a, b, c) */
1035
mp_result mp_int_expt(mp_int a, mp_small b, mp_int c)
1039
unsigned int v = abs(b);
1041
CHECK(b >= 0 && c != NULL);
1043
if((res = mp_int_init_copy(&t, a)) != MP_OK)
1046
(void) mp_int_set_value(c, 1);
1049
if((res = mp_int_mul(c, &t, c)) != MP_OK)
1056
if((res = mp_int_sqr(&t, &t)) != MP_OK)
1067
/* {{{ mp_int_expt_value(a, b, c) */
1069
mp_result mp_int_expt_value(mp_small a, mp_small b, mp_int c)
1073
unsigned int v = abs(b);
1075
CHECK(b >= 0 && c != NULL);
1077
if((res = mp_int_init_value(&t, a)) != MP_OK)
1080
(void) mp_int_set_value(c, 1);
1083
if((res = mp_int_mul(c, &t, c)) != MP_OK)
1090
if((res = mp_int_sqr(&t, &t)) != MP_OK)
1101
/* {{{ mp_int_compare(a, b) */
1103
int mp_int_compare(mp_int a, mp_int b)
1107
CHECK(a != NULL && b != NULL);
1110
if(sa == MP_SIGN(b)) {
1111
int cmp = s_ucmp(a, b);
1113
/* If they're both zero or positive, the normal comparison
1114
applies; if both negative, the sense is reversed. */
1131
/* {{{ mp_int_compare_unsigned(a, b) */
1133
int mp_int_compare_unsigned(mp_int a, mp_int b)
1135
NRCHECK(a != NULL && b != NULL);
1137
return s_ucmp(a, b);
1142
/* {{{ mp_int_compare_zero(z) */
1144
int mp_int_compare_zero(mp_int z)
1148
if(MP_USED(z) == 1 && z->digits[0] == 0)
1150
else if(MP_SIGN(z) == MP_ZPOS)
1158
/* {{{ mp_int_compare_value(z, value) */
1160
int mp_int_compare_value(mp_int z, mp_small value)
1162
mp_sign vsign = (value < 0) ? MP_NEG : MP_ZPOS;
1167
if(vsign == MP_SIGN(z)) {
1168
cmp = s_vcmp(z, value);
1170
if(vsign == MP_ZPOS)
1185
/* {{{ mp_int_exptmod(a, b, m, c) */
1187
mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c)
1195
CHECK(a != NULL && b != NULL && c != NULL && m != NULL);
1197
/* Zero moduli and negative exponents are not considered. */
1204
SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
1205
SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
1207
if(c == b || c == m) {
1208
SETUP(mp_int_init_size(TEMP(2), 2 * um), last);
1215
if((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;
1217
if((res = s_brmu(TEMP(1), m)) != MP_OK) goto CLEANUP;
1219
if((res = s_embar(TEMP(0), b, m, TEMP(1), s)) != MP_OK)
1222
res = mp_int_copy(s, c);
1226
mp_int_clear(TEMP(last));
1233
/* {{{ mp_int_exptmod_evalue(a, value, m, c) */
1235
mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c)
1238
mp_digit vbuf[MP_VALUE_DIGITS(value)];
1240
s_fake(&vtmp, value, vbuf);
1242
return mp_int_exptmod(a, &vtmp, m, c);
1247
/* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
1249
mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b,
1253
mp_digit vbuf[MP_VALUE_DIGITS(value)];
1255
s_fake(&vtmp, value, vbuf);
1257
return mp_int_exptmod(&vtmp, b, m, c);
1262
/* {{{ mp_int_exptmod_known(a, b, m, mu, c) */
1264
mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
1272
CHECK(a && b && m && c);
1274
/* Zero moduli and negative exponents are not considered. */
1281
SETUP(mp_int_init_size(TEMP(0), 2 * um), last);
1283
if(c == b || c == m) {
1284
SETUP(mp_int_init_size(TEMP(1), 2 * um), last);
1291
if((res = mp_int_mod(a, m, TEMP(0))) != MP_OK) goto CLEANUP;
1293
if((res = s_embar(TEMP(0), b, m, mu, s)) != MP_OK)
1296
res = mp_int_copy(s, c);
1300
mp_int_clear(TEMP(last));
1307
/* {{{ mp_int_redux_const(m, c) */
1309
mp_result mp_int_redux_const(mp_int m, mp_int c)
1311
CHECK(m != NULL && c != NULL && m != c);
1313
return s_brmu(c, m);
1318
/* {{{ mp_int_invmod(a, m, c) */
1320
mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c)
1327
CHECK(a != NULL && m != NULL && c != NULL);
1329
if(CMPZ(a) == 0 || CMPZ(m) <= 0)
1332
sa = MP_SIGN(a); /* need this for the result later */
1334
for(last = 0; last < 2; ++last)
1335
mp_int_init(TEMP(last));
1337
if((res = mp_int_egcd(a, m, TEMP(0), TEMP(1), NULL)) != MP_OK)
1340
if(mp_int_compare_value(TEMP(0), 1) != 0) {
1345
/* It is first necessary to constrain the value to the proper range */
1346
if((res = mp_int_mod(TEMP(1), m, TEMP(1))) != MP_OK)
1349
/* Now, if 'a' was originally negative, the value we have is
1350
actually the magnitude of the negative representative; to get the
1351
positive value we have to subtract from the modulus. Otherwise,
1352
the value is okay as it stands.
1355
res = mp_int_sub(m, TEMP(1), c);
1357
res = mp_int_copy(TEMP(1), c);
1361
mp_int_clear(TEMP(last));
1368
/* {{{ mp_int_gcd(a, b, c) */
1370
/* Binary GCD algorithm due to Josef Stein, 1961 */
1371
mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c)
1377
CHECK(a != NULL && b != NULL && c != NULL);
1381
if(ca == 0 && cb == 0)
1384
return mp_int_abs(b, c);
1386
return mp_int_abs(a, c);
1389
if((res = mp_int_init_copy(&u, a)) != MP_OK)
1391
if((res = mp_int_init_copy(&v, b)) != MP_OK)
1394
MP_SIGN(&u) = MP_ZPOS; MP_SIGN(&v) = MP_ZPOS;
1396
{ /* Divide out common factors of 2 from u and v */
1397
int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v);
1399
k = MIN(div2_u, div2_v);
1400
s_qdiv(&u, (mp_size) k);
1401
s_qdiv(&v, (mp_size) k);
1404
if(mp_int_is_odd(&u)) {
1405
if((res = mp_int_neg(&v, &t)) != MP_OK)
1409
if((res = mp_int_copy(&u, &t)) != MP_OK)
1414
s_qdiv(&t, s_dp2k(&t));
1417
if((res = mp_int_copy(&t, &u)) != MP_OK)
1421
if((res = mp_int_neg(&t, &v)) != MP_OK)
1425
if((res = mp_int_sub(&u, &v, &t)) != MP_OK)
1432
if((res = mp_int_abs(&u, c)) != MP_OK)
1434
if(!s_qmul(c, (mp_size) k))
1439
V: mp_int_clear(&u);
1440
U: mp_int_clear(&t);
1447
/* {{{ mp_int_egcd(a, b, c, x, y) */
1449
/* This is the binary GCD algorithm again, but this time we keep track
1450
of the elementary matrix operations as we go, so we can get values
1451
x and y satisfying c = ax + by.
1453
mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c,
1456
int k, last = 0, ca, cb;
1460
CHECK(a != NULL && b != NULL && c != NULL &&
1461
(x != NULL || y != NULL));
1465
if(ca == 0 && cb == 0)
1468
if((res = mp_int_abs(b, c)) != MP_OK) return res;
1469
mp_int_zero(x); (void) mp_int_set_value(y, 1); return MP_OK;
1472
if((res = mp_int_abs(a, c)) != MP_OK) return res;
1473
(void) mp_int_set_value(x, 1); mp_int_zero(y); return MP_OK;
1476
/* Initialize temporaries:
1477
A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */
1478
for(last = 0; last < 4; ++last)
1479
mp_int_init(TEMP(last));
1480
TEMP(0)->digits[0] = 1;
1481
TEMP(3)->digits[0] = 1;
1483
SETUP(mp_int_init_copy(TEMP(4), a), last);
1484
SETUP(mp_int_init_copy(TEMP(5), b), last);
1486
/* We will work with absolute values here */
1487
MP_SIGN(TEMP(4)) = MP_ZPOS;
1488
MP_SIGN(TEMP(5)) = MP_ZPOS;
1490
{ /* Divide out common factors of 2 from u and v */
1491
int div2_u = s_dp2k(TEMP(4)), div2_v = s_dp2k(TEMP(5));
1493
k = MIN(div2_u, div2_v);
1498
SETUP(mp_int_init_copy(TEMP(6), TEMP(4)), last);
1499
SETUP(mp_int_init_copy(TEMP(7), TEMP(5)), last);
1502
while(mp_int_is_even(TEMP(4))) {
1505
if(mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) {
1506
if((res = mp_int_add(TEMP(0), TEMP(7), TEMP(0))) != MP_OK)
1508
if((res = mp_int_sub(TEMP(1), TEMP(6), TEMP(1))) != MP_OK)
1516
while(mp_int_is_even(TEMP(5))) {
1519
if(mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) {
1520
if((res = mp_int_add(TEMP(2), TEMP(7), TEMP(2))) != MP_OK)
1522
if((res = mp_int_sub(TEMP(3), TEMP(6), TEMP(3))) != MP_OK)
1530
if(mp_int_compare(TEMP(4), TEMP(5)) >= 0) {
1531
if((res = mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != MP_OK) goto CLEANUP;
1532
if((res = mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != MP_OK) goto CLEANUP;
1533
if((res = mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != MP_OK) goto CLEANUP;
1536
if((res = mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != MP_OK) goto CLEANUP;
1537
if((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK) goto CLEANUP;
1538
if((res = mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != MP_OK) goto CLEANUP;
1541
if(CMPZ(TEMP(4)) == 0) {
1542
if(x && (res = mp_int_copy(TEMP(2), x)) != MP_OK) goto CLEANUP;
1543
if(y && (res = mp_int_copy(TEMP(3), y)) != MP_OK) goto CLEANUP;
1545
if(!s_qmul(TEMP(5), k)) {
1550
res = mp_int_copy(TEMP(5), c);
1559
mp_int_clear(TEMP(last));
1566
/* {{{ mp_int_lcm(a, b, c) */
1568
mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c)
1573
CHECK(a != NULL && b != NULL && c != NULL);
1575
/* Since a * b = gcd(a, b) * lcm(a, b), we can compute
1576
lcm(a, b) = (a / gcd(a, b)) * b.
1578
This formulation insures everything works even if the input
1579
variables share space.
1581
if((res = mp_int_init(&lcm)) != MP_OK)
1583
if((res = mp_int_gcd(a, b, &lcm)) != MP_OK)
1585
if((res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK)
1587
if((res = mp_int_mul(&lcm, b, &lcm)) != MP_OK)
1590
res = mp_int_copy(&lcm, c);
1600
/* {{{ mp_int_divisible_value(a, v) */
1602
int mp_int_divisible_value(mp_int a, mp_small v)
1606
if(mp_int_div_value(a, v, NULL, &rem) != MP_OK)
1614
/* {{{ mp_int_is_pow2(z) */
1616
int mp_int_is_pow2(mp_int z)
1625
/* {{{ mp_int_root(a, b, c) */
1627
/* Implementation of Newton's root finding method, based loosely on a
1628
patch contributed by Hal Finkel <half@halssoftware.com>
1629
modified by M. J. Fromberger.
1631
mp_result mp_int_root(mp_int a, mp_small b, mp_int c)
1633
mp_result res = MP_OK;
1638
CHECK(a != NULL && c != NULL && b > 0);
1641
return mp_int_copy(a, c);
1643
if(MP_SIGN(a) == MP_NEG) {
1645
return MP_UNDEF; /* root does not exist for negative a with even b */
1650
SETUP(mp_int_init_copy(TEMP(last), a), last);
1651
SETUP(mp_int_init_copy(TEMP(last), a), last);
1652
SETUP(mp_int_init(TEMP(last)), last);
1653
SETUP(mp_int_init(TEMP(last)), last);
1654
SETUP(mp_int_init(TEMP(last)), last);
1656
(void) mp_int_abs(TEMP(0), TEMP(0));
1657
(void) mp_int_abs(TEMP(1), TEMP(1));
1660
if((res = mp_int_expt(TEMP(1), b, TEMP(2))) != MP_OK)
1663
if(mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0)
1666
if((res = mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK)
1668
if((res = mp_int_expt(TEMP(1), b - 1, TEMP(3))) != MP_OK)
1670
if((res = mp_int_mul_value(TEMP(3), b, TEMP(3))) != MP_OK)
1672
if((res = mp_int_div(TEMP(2), TEMP(3), TEMP(4), NULL)) != MP_OK)
1674
if((res = mp_int_sub(TEMP(1), TEMP(4), TEMP(4))) != MP_OK)
1677
if(mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) {
1678
if((res = mp_int_sub_value(TEMP(4), 1, TEMP(4))) != MP_OK)
1681
if((res = mp_int_copy(TEMP(4), TEMP(1))) != MP_OK)
1685
if((res = mp_int_copy(TEMP(1), c)) != MP_OK)
1688
/* If the original value of a was negative, flip the output sign. */
1690
(void) mp_int_neg(c, c); /* cannot fail */
1694
mp_int_clear(TEMP(last));
1701
/* {{{ mp_int_to_int(z, out) */
1703
mp_result mp_int_to_int(mp_int z, mp_small *out)
1712
/* Make sure the value is representable as an int */
1714
if((sz == MP_ZPOS && mp_int_compare_value(z, MP_SMALL_MAX) > 0) ||
1715
mp_int_compare_value(z, MP_SMALL_MIN) < 0)
1719
dz = MP_DIGITS(z) + uz - 1;
1722
uv <<= MP_DIGIT_BIT/2;
1723
uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1728
*out = (sz == MP_NEG) ? -(mp_small)uv : (mp_small)uv;
1735
/* {{{ mp_int_to_uint(z, *out) */
1737
mp_result mp_int_to_uint(mp_int z, mp_usmall *out)
1746
/* Make sure the value is representable as an int */
1748
if(!(sz == MP_ZPOS && mp_int_compare_value(z, UINT_MAX) <= 0))
1752
dz = MP_DIGITS(z) + uz - 1;
1755
uv <<= MP_DIGIT_BIT/2;
1756
uv = (uv << (MP_DIGIT_BIT/2)) | *dz--;
1768
/* {{{ mp_int_to_string(z, radix, str, limit) */
1770
mp_result mp_int_to_string(mp_int z, mp_size radix,
1771
char *str, int limit)
1776
CHECK(z != NULL && str != NULL && limit >= 2);
1778
if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1782
*str++ = s_val2ch(0, 1);
1788
if((res = mp_int_init_copy(&tmp, z)) != MP_OK)
1791
if(MP_SIGN(z) == MP_NEG) {
1797
/* Generate digits in reverse order until finished or limit reached */
1798
for(/* */; limit > 0; --limit) {
1801
if((cmp = CMPZ(&tmp)) == 0)
1804
d = s_ddiv(&tmp, (mp_digit)radix);
1805
*str++ = s_val2ch(d, 1);
1809
/* Put digits back in correct output order */
1828
/* {{{ mp_int_string_len(z, radix) */
1830
mp_result mp_int_string_len(mp_int z, mp_size radix)
1836
if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1839
len = s_outlen(z, radix) + 1; /* for terminator */
1841
/* Allow for sign marker on negatives */
1842
if(MP_SIGN(z) == MP_NEG)
1850
/* {{{ mp_int_read_string(z, radix, *str) */
1852
/* Read zero-terminated string into z */
1853
mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str)
1855
return mp_int_read_cstring(z, radix, str, NULL);
1861
/* {{{ mp_int_read_cstring(z, radix, *str, **end) */
1863
mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end)
1867
CHECK(z != NULL && str != NULL);
1869
if(radix < MP_MIN_RADIX || radix > MP_MAX_RADIX)
1872
/* Skip leading whitespace */
1873
while(isspace((int)*str))
1876
/* Handle leading sign tag (+/-, positive default) */
1879
MP_SIGN(z) = MP_NEG;
1883
++str; /* fallthrough */
1885
MP_SIGN(z) = MP_ZPOS;
1889
/* Skip leading zeroes */
1890
while((ch = s_ch2val(*str, radix)) == 0)
1893
/* Make sure there is enough space for the value */
1894
if(!s_pad(z, s_inlen(strlen(str), radix)))
1897
MP_USED(z) = 1; z->digits[0] = 0;
1899
while(*str != '\0' && ((ch = s_ch2val(*str, radix)) >= 0)) {
1900
s_dmul(z, (mp_digit)radix);
1901
s_dadd(z, (mp_digit)ch);
1907
/* Override sign for zero, even if negative specified. */
1909
MP_SIGN(z) = MP_ZPOS;
1914
/* Return a truncation error if the string has unprocessed
1915
characters remaining, so the caller can tell if the whole string
1925
/* {{{ mp_int_count_bits(z) */
1927
mp_result mp_int_count_bits(mp_int z)
1929
mp_size nbits = 0, uz;
1935
if(uz == 1 && z->digits[0] == 0)
1939
nbits = uz * MP_DIGIT_BIT;
1952
/* {{{ mp_int_to_binary(z, buf, limit) */
1954
mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit)
1956
static const int PAD_FOR_2C = 1;
1961
CHECK(z != NULL && buf != NULL);
1963
res = s_tobin(z, buf, &limpos, PAD_FOR_2C);
1965
if(MP_SIGN(z) == MP_NEG)
1966
s_2comp(buf, limpos);
1973
/* {{{ mp_int_read_binary(z, buf, len) */
1975
mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len)
1981
CHECK(z != NULL && buf != NULL && len > 0);
1983
/* Figure out how many digits are needed to represent this value */
1984
need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
1990
/* If the high-order bit is set, take the 2's complement before
1991
reading the value (it will be restored afterward) */
1992
if(buf[0] >> (CHAR_BIT - 1)) {
1993
MP_SIGN(z) = MP_NEG;
1998
for(tmp = buf, i = len; i > 0; --i, ++tmp) {
1999
s_qmul(z, (mp_size) CHAR_BIT);
2003
/* Restore 2's complement if we took it before */
2004
if(MP_SIGN(z) == MP_NEG)
2012
/* {{{ mp_int_binary_len(z) */
2014
mp_result mp_int_binary_len(mp_int z)
2016
mp_result res = mp_int_count_bits(z);
2017
int bytes = mp_int_unsigned_len(z);
2022
bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
2024
/* If the highest-order bit falls exactly on a byte boundary, we
2025
need to pad with an extra byte so that the sign will be read
2026
correctly when reading it back in. */
2027
if(bytes * CHAR_BIT == res)
2035
/* {{{ mp_int_to_unsigned(z, buf, limit) */
2037
mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit)
2039
static const int NO_PADDING = 0;
2041
CHECK(z != NULL && buf != NULL);
2043
return s_tobin(z, buf, &limit, NO_PADDING);
2048
/* {{{ mp_int_read_unsigned(z, buf, len) */
2050
mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len)
2056
CHECK(z != NULL && buf != NULL && len > 0);
2058
/* Figure out how many digits are needed to represent this value */
2059
need = ((len * CHAR_BIT) + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT;
2066
for(tmp = buf, i = len; i > 0; --i, ++tmp) {
2067
(void) s_qmul(z, CHAR_BIT);
2076
/* {{{ mp_int_unsigned_len(z) */
2078
mp_result mp_int_unsigned_len(mp_int z)
2080
mp_result res = mp_int_count_bits(z);
2086
bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT;
2093
/* {{{ mp_error_string(res) */
2095
const char *mp_error_string(mp_result res)
2099
return s_unknown_err;
2102
for(ix = 0; ix < res && s_error_msg[ix] != NULL; ++ix)
2105
if(s_error_msg[ix] != NULL)
2106
return s_error_msg[ix];
2108
return s_unknown_err;
2113
/*------------------------------------------------------------------------*/
2114
/* Private functions for internal use. These make assumptions. */
2116
/* {{{ s_alloc(num) */
2118
static mp_digit *s_alloc(mp_size num)
2120
mp_digit *out = malloc(num * sizeof(mp_digit));
2122
assert(out != NULL); /* for debugging */
2125
mp_digit v = (mp_digit) 0xdeadbeef;
2128
for(ix = 0; ix < num; ++ix)
2138
/* {{{ s_realloc(old, osize, nsize) */
2140
static mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize)
2143
mp_digit *new = s_alloc(nsize);
2146
for(ix = 0; ix < nsize; ++ix)
2147
new[ix] = (mp_digit) 0xdeadbeef;
2149
memcpy(new, old, osize * sizeof(mp_digit));
2151
mp_digit *new = realloc(old, nsize * sizeof(mp_digit));
2153
assert(new != NULL); /* for debugging */
2160
/* {{{ s_free(ptr) */
2162
static void s_free(void *ptr)
2169
/* {{{ s_pad(z, min) */
2171
static int s_pad(mp_int z, mp_size min)
2173
if(MP_ALLOC(z) < min) {
2174
mp_size nsize = ROUND_PREC(min);
2177
if((void *)z->digits == (void *)z) {
2178
if((tmp = s_alloc(nsize)) == NULL)
2181
COPY(MP_DIGITS(z), tmp, MP_USED(z));
2183
else if((tmp = s_realloc(MP_DIGITS(z), MP_ALLOC(z), nsize)) == NULL)
2187
MP_ALLOC(z) = nsize;
2195
/* {{{ s_fake(z, value, vbuf) */
2197
static void s_fake(mp_int z, mp_small value, mp_digit vbuf[])
2199
mp_size uv = (mp_size) s_vpack(value, vbuf);
2202
z->alloc = MP_VALUE_DIGITS(value);
2203
z->sign = (value < 0) ? MP_NEG : MP_ZPOS;
2209
/* {{{ s_cdig(da, db, len) */
2211
static int s_cdig(mp_digit *da, mp_digit *db, mp_size len)
2213
mp_digit *dat = da + len - 1, *dbt = db + len - 1;
2215
for(/* */; len != 0; --len, --dat, --dbt) {
2218
else if(*dat < *dbt)
2227
/* {{{ s_vpack(v, t[]) */
2229
static int s_vpack(mp_small v, mp_digit t[])
2231
mp_usmall uv = (mp_usmall) ((v < 0) ? -v : v);
2238
t[ndig++] = (mp_digit) uv;
2239
uv >>= MP_DIGIT_BIT/2;
2240
uv >>= MP_DIGIT_BIT/2;
2249
/* {{{ s_ucmp(a, b) */
2251
static int s_ucmp(mp_int a, mp_int b)
2253
mp_size ua = MP_USED(a), ub = MP_USED(b);
2260
return s_cdig(MP_DIGITS(a), MP_DIGITS(b), ua);
2265
/* {{{ s_vcmp(a, v) */
2267
static int s_vcmp(mp_int a, mp_small v)
2269
mp_digit vdig[MP_VALUE_DIGITS(v)];
2271
mp_size ua = MP_USED(a);
2273
ndig = s_vpack(v, vdig);
2280
return s_cdig(MP_DIGITS(a), vdig, ndig);
2285
/* {{{ s_uadd(da, db, dc, size_a, size_b) */
2287
static mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc,
2288
mp_size size_a, mp_size size_b)
2293
/* Insure that da is the longer of the two to simplify later code */
2294
if(size_b > size_a) {
2295
SWAP(mp_digit *, da, db);
2296
SWAP(mp_size, size_a, size_b);
2299
/* Add corresponding digits until the shorter number runs out */
2300
for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) {
2301
w = w + (mp_word) *da + (mp_word) *db;
2302
*dc = LOWER_HALF(w);
2306
/* Propagate carries as far as necessary */
2307
for(/* */; pos < size_a; ++pos, ++da, ++dc) {
2310
*dc = LOWER_HALF(w);
2314
/* Return carry out */
2320
/* {{{ s_usub(da, db, dc, size_a, size_b) */
2322
static void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc,
2323
mp_size size_a, mp_size size_b)
2328
/* We assume that |a| >= |b| so this should definitely hold */
2329
assert(size_a >= size_b);
2331
/* Subtract corresponding digits and propagate borrow */
2332
for(pos = 0; pos < size_b; ++pos, ++da, ++db, ++dc) {
2333
w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */
2334
(mp_word)*da) - w - (mp_word)*db;
2336
*dc = LOWER_HALF(w);
2337
w = (UPPER_HALF(w) == 0);
2340
/* Finish the subtraction for remaining upper digits of da */
2341
for(/* */; pos < size_a; ++pos, ++da, ++dc) {
2342
w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */
2345
*dc = LOWER_HALF(w);
2346
w = (UPPER_HALF(w) == 0);
2349
/* If there is a borrow out at the end, it violates the precondition */
2355
/* {{{ s_kmul(da, db, dc, size_a, size_b) */
2357
static int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc,
2358
mp_size size_a, mp_size size_b)
2362
/* Make sure b is the smaller of the two input values */
2363
if(size_b > size_a) {
2364
SWAP(mp_digit *, da, db);
2365
SWAP(mp_size, size_a, size_b);
2368
/* Insure that the bottom is the larger half in an odd-length split;
2369
the code below relies on this being true.
2371
bot_size = (size_a + 1) / 2;
2373
/* If the values are big enough to bother with recursion, use the
2374
Karatsuba algorithm to compute the product; otherwise use the
2375
normal multiplication algorithm
2377
if(multiply_threshold &&
2378
size_a >= multiply_threshold &&
2379
size_b > bot_size) {
2381
mp_digit *t1, *t2, *t3, carry;
2383
mp_digit *a_top = da + bot_size;
2384
mp_digit *b_top = db + bot_size;
2386
mp_size at_size = size_a - bot_size;
2387
mp_size bt_size = size_b - bot_size;
2388
mp_size buf_size = 2 * bot_size;
2390
/* Do a single allocation for all three temporary buffers needed;
2391
each buffer must be big enough to hold the product of two
2392
bottom halves, and one buffer needs space for the completed
2393
product; twice the space is plenty.
2395
if((t1 = s_alloc(4 * buf_size)) == NULL) return 0;
2398
ZERO(t1, 4 * buf_size);
2400
/* t1 and t2 are initially used as temporaries to compute the inner product
2401
(a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
2403
carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */
2404
t1[bot_size] = carry;
2406
carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */
2407
t2[bot_size] = carry;
2409
(void) s_kmul(t1, t2, t3, bot_size + 1, bot_size + 1); /* t3 = t1 * t2 */
2411
/* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that
2412
we're left with only the pieces we want: t3 = a1b0 + a0b1
2416
(void) s_kmul(da, db, t1, bot_size, bot_size); /* t1 = a0 * b0 */
2417
(void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */
2419
/* Subtract out t1 and t2 to get the inner product */
2420
s_usub(t3, t1, t3, buf_size + 2, buf_size);
2421
s_usub(t3, t2, t3, buf_size + 2, buf_size);
2423
/* Assemble the output value */
2424
COPY(t1, dc, buf_size);
2425
carry = s_uadd(t3, dc + bot_size, dc + bot_size,
2426
buf_size + 1, buf_size);
2429
carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
2430
buf_size, buf_size);
2433
s_free(t1); /* note t2 and t3 are just internal pointers to t1 */
2436
s_umul(da, db, dc, size_a, size_b);
2444
/* {{{ s_umul(da, db, dc, size_a, size_b) */
2446
static void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc,
2447
mp_size size_a, mp_size size_b)
2452
for(a = 0; a < size_a; ++a, ++dc, ++da) {
2460
for(b = 0; b < size_b; ++b, ++dbt, ++dct) {
2461
w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct;
2463
*dct = LOWER_HALF(w);
2473
/* {{{ s_ksqr(da, dc, size_a) */
2475
static int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2477
if(multiply_threshold && size_a > multiply_threshold) {
2478
mp_size bot_size = (size_a + 1) / 2;
2479
mp_digit *a_top = da + bot_size;
2480
mp_digit *t1, *t2, *t3, carry;
2481
mp_size at_size = size_a - bot_size;
2482
mp_size buf_size = 2 * bot_size;
2484
if((t1 = s_alloc(4 * buf_size)) == NULL) return 0;
2487
ZERO(t1, 4 * buf_size);
2489
(void) s_ksqr(da, t1, bot_size); /* t1 = a0 ^ 2 */
2490
(void) s_ksqr(a_top, t2, at_size); /* t2 = a1 ^ 2 */
2492
(void) s_kmul(da, a_top, t3, bot_size, at_size); /* t3 = a0 * a1 */
2494
/* Quick multiply t3 by 2, shifting left (can't overflow) */
2496
int i, top = bot_size + at_size;
2497
mp_word w, save = 0;
2499
for(i = 0; i < top; ++i) {
2501
w = (w << 1) | save;
2502
t3[i] = LOWER_HALF(w);
2503
save = UPPER_HALF(w);
2505
t3[i] = LOWER_HALF(save);
2508
/* Assemble the output value */
2509
COPY(t1, dc, 2 * bot_size);
2510
carry = s_uadd(t3, dc + bot_size, dc + bot_size,
2511
buf_size + 1, buf_size);
2514
carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size,
2515
buf_size, buf_size);
2518
s_free(t1); /* note that t2 and t2 are internal pointers only */
2522
s_usqr(da, dc, size_a);
2530
/* {{{ s_usqr(da, dc, size_a) */
2532
static void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a)
2537
for(i = 0; i < size_a; ++i, dc += 2, ++da) {
2538
mp_digit *dct = dc, *dat = da;
2543
/* Take care of the first digit, no rollover */
2544
w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct;
2545
*dct = LOWER_HALF(w);
2549
for(j = i + 1; j < size_a; ++j, ++dat, ++dct) {
2550
mp_word t = (mp_word)*da * (mp_word)*dat;
2551
mp_word u = w + (mp_word)*dct, ov = 0;
2553
/* Check if doubling t will overflow a word */
2559
/* Check if adding u to w will overflow a word */
2560
if(ADD_WILL_OVERFLOW(w, u))
2565
*dct = LOWER_HALF(w);
2568
w += MP_DIGIT_MAX; /* MP_RADIX */
2575
while((w = UPPER_HALF(w)) != 0) {
2576
++dct; w = w + *dct;
2577
*dct = LOWER_HALF(w);
2586
/* {{{ s_dadd(a, b) */
2588
static void s_dadd(mp_int a, mp_digit b)
2591
mp_digit *da = MP_DIGITS(a);
2592
mp_size ua = MP_USED(a);
2594
w = (mp_word)*da + b;
2595
*da++ = LOWER_HALF(w);
2598
for(ua -= 1; ua > 0; --ua, ++da) {
2599
w = (mp_word)*da + w;
2601
*da = LOWER_HALF(w);
2613
/* {{{ s_dmul(a, b) */
2615
static void s_dmul(mp_int a, mp_digit b)
2618
mp_digit *da = MP_DIGITS(a);
2619
mp_size ua = MP_USED(a);
2622
w = (mp_word)*da * b + w;
2623
*da++ = LOWER_HALF(w);
2636
/* {{{ s_dbmul(da, b, dc, size_a) */
2638
static void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a)
2643
w = (mp_word)*da++ * (mp_word)b + w;
2645
*dc++ = LOWER_HALF(w);
2651
*dc = LOWER_HALF(w);
2656
/* {{{ s_ddiv(da, d, dc, size_a) */
2658
static mp_digit s_ddiv(mp_int a, mp_digit b)
2660
mp_word w = 0, qdigit;
2661
mp_size ua = MP_USED(a);
2662
mp_digit *da = MP_DIGITS(a) + ua - 1;
2664
for(/* */; ua > 0; --ua, --da) {
2665
w = (w << MP_DIGIT_BIT) | *da;
2675
*da = (mp_digit)qdigit;
2684
/* {{{ s_qdiv(z, p2) */
2686
static void s_qdiv(mp_int z, mp_size p2)
2688
mp_size ndig = p2 / MP_DIGIT_BIT, nbits = p2 % MP_DIGIT_BIT;
2689
mp_size uz = MP_USED(z);
2693
mp_digit *to, *from;
2700
to = MP_DIGITS(z); from = to + ndig;
2702
for(mark = ndig; mark < uz; ++mark)
2705
MP_USED(z) = uz - ndig;
2709
mp_digit d = 0, *dz, save;
2710
mp_size up = MP_DIGIT_BIT - nbits;
2713
dz = MP_DIGITS(z) + uz - 1;
2715
for(/* */; uz > 0; --uz, --dz) {
2718
*dz = (*dz >> nbits) | (d << up);
2725
if(MP_USED(z) == 1 && z->digits[0] == 0)
2726
MP_SIGN(z) = MP_ZPOS;
2731
/* {{{ s_qmod(z, p2) */
2733
static void s_qmod(mp_int z, mp_size p2)
2735
mp_size start = p2 / MP_DIGIT_BIT + 1, rest = p2 % MP_DIGIT_BIT;
2736
mp_size uz = MP_USED(z);
2737
mp_digit mask = (1 << rest) - 1;
2741
z->digits[start - 1] &= mask;
2748
/* {{{ s_qmul(z, p2) */
2750
static int s_qmul(mp_int z, mp_size p2)
2752
mp_size uz, need, rest, extra, i;
2753
mp_digit *from, *to, d;
2759
need = p2 / MP_DIGIT_BIT; rest = p2 % MP_DIGIT_BIT;
2761
/* Figure out if we need an extra digit at the top end; this occurs
2762
if the topmost `rest' bits of the high-order digit of z are not
2763
zero, meaning they will be shifted off the end if not preserved */
2766
mp_digit *dz = MP_DIGITS(z) + uz - 1;
2768
if((*dz >> (MP_DIGIT_BIT - rest)) != 0)
2772
if(!s_pad(z, uz + need + extra))
2775
/* If we need to shift by whole digits, do that in one pass, then
2776
to back and shift by partial digits.
2779
from = MP_DIGITS(z) + uz - 1;
2782
for(i = 0; i < uz; ++i)
2785
ZERO(MP_DIGITS(z), need);
2791
for(i = need, from = MP_DIGITS(z) + need; i < uz; ++i, ++from) {
2792
mp_digit save = *from;
2794
*from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest));
2798
d >>= (MP_DIGIT_BIT - rest);
2813
/* {{{ s_qsub(z, p2) */
2815
/* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
2816
The sign of the result is always zero/positive.
2818
static int s_qsub(mp_int z, mp_size p2)
2820
mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp;
2821
mp_size tdig = (p2 / MP_DIGIT_BIT), pos;
2824
if(!s_pad(z, tdig + 1))
2827
for(pos = 0, zp = MP_DIGITS(z); pos < tdig; ++pos, ++zp) {
2828
w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp;
2830
*zp = LOWER_HALF(w);
2831
w = UPPER_HALF(w) ? 0 : 1;
2834
w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp;
2835
*zp = LOWER_HALF(w);
2837
assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */
2839
MP_SIGN(z) = MP_ZPOS;
2849
static int s_dp2k(mp_int z)
2852
mp_digit *dp = MP_DIGITS(z), d;
2854
if(MP_USED(z) == 1 && *dp == 0)
2863
while((d & 1) == 0) {
2875
static int s_isp2(mp_int z)
2877
mp_size uz = MP_USED(z), k = 0;
2878
mp_digit *dz = MP_DIGITS(z), d;
2899
/* {{{ s_2expt(z, k) */
2901
static int s_2expt(mp_int z, mp_small k)
2906
ndig = (k + MP_DIGIT_BIT) / MP_DIGIT_BIT;
2907
rest = k % MP_DIGIT_BIT;
2914
*(dz + ndig - 1) = (1 << rest);
2922
/* {{{ s_norm(a, b) */
2924
static int s_norm(mp_int a, mp_int b)
2926
mp_digit d = b->digits[MP_USED(b) - 1];
2929
while(d < (mp_digit) (1 << (MP_DIGIT_BIT - 1))) { /* d < (MP_RADIX / 2) */
2934
/* These multiplications can't fail */
2936
(void) s_qmul(a, (mp_size) k);
2937
(void) s_qmul(b, (mp_size) k);
2945
/* {{{ s_brmu(z, m) */
2947
static mp_result s_brmu(mp_int z, mp_int m)
2949
mp_size um = MP_USED(m) * 2;
2954
s_2expt(z, MP_DIGIT_BIT * um);
2955
return mp_int_div(z, m, z, NULL);
2960
/* {{{ s_reduce(x, m, mu, q1, q2) */
2962
static int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2)
2964
mp_size um = MP_USED(m), umb_p1, umb_m1;
2966
umb_p1 = (um + 1) * MP_DIGIT_BIT;
2967
umb_m1 = (um - 1) * MP_DIGIT_BIT;
2969
if(mp_int_copy(x, q1) != MP_OK)
2972
/* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
2977
/* Set x = x mod b^(k+1) */
2980
/* Now, q is a guess for the quotient a / m.
2981
Compute x - q * m mod b^(k+1), replacing x. This may be off
2982
by a factor of 2m, but no more than that.
2986
(void) mp_int_sub(x, q1, x); /* can't fail */
2988
/* The result may be < 0; if it is, add b^(k+1) to pin it in the
2990
if((CMPZ(x) < 0) && !s_qsub(x, umb_p1))
2993
/* If x > m, we need to back it off until it is in range.
2994
This will be required at most twice. */
2995
if(mp_int_compare(x, m) >= 0) {
2996
(void) mp_int_sub(x, m, x);
2997
if(mp_int_compare(x, m) >= 0)
2998
(void) mp_int_sub(x, m, x);
3001
/* At this point, x has been properly reduced. */
3007
/* {{{ s_embar(a, b, m, mu, c) */
3009
/* Perform modular exponentiation using Barrett's method, where mu is
3010
the reduction constant for m. Assumes a < m, b > 0. */
3011
static mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c)
3013
mp_digit *db, *dbt, umu, d;
3018
umu = MP_USED(mu); db = MP_DIGITS(b); dbt = db + MP_USED(b) - 1;
3021
SETUP(mp_int_init_size(TEMP(last), 4 * umu), last);
3022
ZERO(MP_DIGITS(TEMP(last - 1)), MP_ALLOC(TEMP(last - 1)));
3025
(void) mp_int_set_value(c, 1);
3027
/* Take care of low-order digits */
3031
for(d = *db, i = MP_DIGIT_BIT; i > 0; --i, d >>= 1) {
3033
/* The use of a second temporary avoids allocation */
3034
UMUL(c, a, TEMP(0));
3035
if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
3036
res = MP_MEMORY; goto CLEANUP;
3038
mp_int_copy(TEMP(0), c);
3043
assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
3044
if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
3045
res = MP_MEMORY; goto CLEANUP;
3047
assert(MP_SIGN(TEMP(0)) == MP_ZPOS);
3048
mp_int_copy(TEMP(0), a);
3056
/* Take care of highest-order digit */
3060
UMUL(c, a, TEMP(0));
3061
if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
3062
res = MP_MEMORY; goto CLEANUP;
3064
mp_int_copy(TEMP(0), c);
3071
if(!s_reduce(TEMP(0), m, mu, TEMP(1), TEMP(2))) {
3072
res = MP_MEMORY; goto CLEANUP;
3074
(void) mp_int_copy(TEMP(0), a);
3079
mp_int_clear(TEMP(last));
3086
/* {{{ s_udiv(a, b) */
3088
/* Precondition: a >= b and b > 0
3089
Postcondition: a' = a / b, b' = a % b
3091
static mp_result s_udiv(mp_int a, mp_int b)
3094
mp_size ua, ub, qpos = 0;
3096
mp_result res = MP_OK;
3099
/* Force signs to positive */
3100
MP_SIGN(a) = MP_ZPOS;
3101
MP_SIGN(b) = MP_ZPOS;
3103
/* Normalize, per Knuth */
3106
ua = MP_USED(a); ub = MP_USED(b); btop = b->digits[ub - 1];
3107
if((res = mp_int_init_size(&q, ua)) != MP_OK) return res;
3108
if((res = mp_int_init_size(&t, ua + 1)) != MP_OK) goto CLEANUP;
3111
r.digits = da + ua - 1; /* The contents of r are shared with a */
3114
r.alloc = MP_ALLOC(a);
3115
ZERO(t.digits, t.alloc);
3117
/* Solve for quotient digits, store in q.digits in reverse order */
3118
while(r.digits >= da) {
3119
assert(qpos <= q.alloc);
3121
if(s_ucmp(b, &r) > 0) {
3125
if(++skip > 1 && qpos > 0)
3126
q.digits[qpos++] = 0;
3131
mp_word pfx = r.digits[r.used - 1];
3134
if(r.used > 1 && pfx <= btop) {
3135
pfx <<= MP_DIGIT_BIT / 2;
3136
pfx <<= MP_DIGIT_BIT / 2;
3137
pfx |= r.digits[r.used - 2];
3140
qdigit = pfx / btop;
3141
if(qdigit > MP_DIGIT_MAX) {
3142
if(qdigit & MP_DIGIT_MAX)
3143
qdigit = MP_DIGIT_MAX;
3148
s_dbmul(MP_DIGITS(b), (mp_digit) qdigit, t.digits, ub);
3149
t.used = ub + 1; CLAMP(&t);
3150
while(s_ucmp(&t, &r) > 0) {
3152
(void) mp_int_sub(&t, b, &t); /* cannot fail */
3155
s_usub(r.digits, t.digits, r.digits, r.used, t.used);
3158
q.digits[qpos++] = (mp_digit) qdigit;
3159
ZERO(t.digits, t.used);
3164
/* Put quotient digits in the correct order, and discard extra zeroes */
3166
REV(mp_digit, q.digits, qpos);
3169
/* Denormalize the remainder */
3174
mp_int_copy(a, b); /* ok: 0 <= r < b */
3175
mp_int_copy(&q, a); /* ok: q <= a */
3185
/* {{{ s_outlen(z, r) */
3187
static int s_outlen(mp_int z, mp_size r)
3192
assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX);
3194
bits = mp_int_count_bits(z);
3195
raw = (double)bits * s_log2[r];
3197
return (int)(raw + 0.999999);
3202
/* {{{ s_inlen(len, r) */
3204
static mp_size s_inlen(int len, mp_size r)
3206
double raw = (double)len / s_log2[r];
3207
mp_size bits = (mp_size)(raw + 0.5);
3209
return (mp_size)((bits + (MP_DIGIT_BIT - 1)) / MP_DIGIT_BIT);
3214
/* {{{ s_ch2val(c, r) */
3216
static int s_ch2val(char c, int r)
3220
if(isdigit((unsigned char) c))
3222
else if(r > 10 && isalpha((unsigned char) c))
3223
out = toupper(c) - 'A' + 10;
3227
return (out >= r) ? -1 : out;
3232
/* {{{ s_val2ch(v, caps) */
3234
static char s_val2ch(int v, int caps)
3241
char out = (v - 10) + 'a';
3244
return toupper(out);
3252
/* {{{ s_2comp(buf, len) */
3254
static void s_2comp(unsigned char *buf, int len)
3257
unsigned short s = 1;
3259
for(i = len - 1; i >= 0; --i) {
3260
unsigned char c = ~buf[i];
3269
/* last carry out is ignored */
3274
/* {{{ s_tobin(z, buf, *limpos) */
3276
static mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad)
3280
int pos = 0, limit = *limpos;
3282
uz = MP_USED(z); dz = MP_DIGITS(z);
3283
while(uz > 0 && pos < limit) {
3287
for(i = sizeof(mp_digit); i > 0 && pos < limit; --i) {
3288
buf[pos++] = (unsigned char)d;
3291
/* Don't write leading zeroes */
3292
if(d == 0 && uz == 1)
3293
i = 0; /* exit loop without signaling truncation */
3296
/* Detect truncation (loop exited with pos >= limit) */
3302
if(pad != 0 && (buf[pos - 1] >> (CHAR_BIT - 1))) {
3309
/* Digits are in reverse order, fix that */
3310
REV(unsigned char, buf, pos);
3312
/* Return the number of bytes actually written */
3315
return (uz == 0) ? MP_OK : MP_TRUNC;
3320
/* {{{ s_print(tag, z) */
3323
void s_print(char *tag, mp_int z)
3327
fprintf(stderr, "%s: %c ", tag,
3328
(MP_SIGN(z) == MP_NEG) ? '-' : '+');
3330
for(i = MP_USED(z) - 1; i >= 0; --i)
3331
fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), z->digits[i]);
3333
fputc('\n', stderr);
3337
void s_print_buf(char *tag, mp_digit *buf, mp_size num)
3341
fprintf(stderr, "%s: ", tag);
3343
for(i = num - 1; i >= 0; --i)
3344
fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), buf[i]);
3346
fputc('\n', stderr);
3352
/* HERE THERE BE DRAGONS */
added
removed
source4/heimdal/lib/hcrypto/imath/imath.c
