3 Purpose: Arbitrary precision integer arithmetic routines.
4 Author: M. J. Fromberger <http://spinning-yarns.org/michael/>
5 Info: $Id: imath.c 826 2009-02-11 16:21:04Z 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
43 #define STATIC /* public */
50 const mp_result MP_OK
= 0; /* no error, all is well */
51 const mp_result MP_FALSE
= 0; /* boolean false */
52 const mp_result MP_TRUE
= -1; /* boolean true */
53 const mp_result MP_MEMORY
= -2; /* out of memory */
54 const mp_result MP_RANGE
= -3; /* argument out of range */
55 const mp_result MP_UNDEF
= -4; /* result undefined */
56 const mp_result MP_TRUNC
= -5; /* output truncated */
57 const mp_result MP_BADARG
= -6; /* invalid null argument */
58 const mp_result MP_MINERR
= -6;
60 const mp_sign MP_NEG
= 1; /* value is strictly negative */
61 const mp_sign MP_ZPOS
= 0; /* value is non-negative */
63 STATIC
const char *s_unknown_err
= "unknown result code";
64 STATIC
const char *s_error_msg
[] = {
68 "argument out of range",
77 /* Argument checking macros
78 Use CHECK() where a return value is required; NRCHECK() elsewhere */
79 #define CHECK(TEST) assert(TEST)
80 #define NRCHECK(TEST) assert(TEST)
82 /* {{{ Logarithm table for computing output sizes */
84 /* The ith entry of this table gives the value of log_i(2).
86 An integer value n requires ceil(log_i(n)) digits to be represented
87 in base i. Since it is easy to compute lg(n), by counting bits, we
88 can compute log_i(n) = lg(n) * log_i(2).
90 The use of this table eliminates a dependency upon linkage against
91 the standard math libraries.
93 STATIC
const double s_log2
[] = {
94 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
95 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
96 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
97 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */
98 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */
99 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */
100 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */
101 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */
102 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */
103 0.193426404, /* 36 */
107 /* {{{ Various macros */
109 /* Return the number of digits needed to represent a static value */
110 #define MP_VALUE_DIGITS(V) \
111 ((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
113 /* Round precision P to nearest word boundary */
114 #define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
116 /* Set array P of S digits to zero */
118 do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P);memset(p__,0,i__);}while(0)
120 /* Copy S digits from array P to array Q */
121 #define COPY(P, Q, S) \
122 do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P),*q__=(Q);\
123 memcpy(q__,p__,i__);}while(0)
125 /* Reverse N elements of type T in array A */
126 #define REV(T, A, N) \
127 do{T *u_=(A),*v_=u_+(N)-1;while(u_<v_){T xch=*u_;*u_++=*v_;*v_--=xch;}}while(0)
130 do{mp_int z_=(Z);mp_size uz_=MP_USED(z_);mp_digit *dz_=MP_DIGITS(z_)+uz_-1;\
131 while(uz_ > 1 && (*dz_-- == 0)) --uz_;MP_USED(z_)=uz_;}while(0)
133 /* Select min/max. Do not provide expressions for which multiple
134 evaluation would be problematic, e.g. x++ */
135 #define MIN(A, B) ((B)<(A)?(B):(A))
136 #define MAX(A, B) ((B)>(A)?(B):(A))
138 /* Exchange lvalues A and B of type T, e.g.
139 SWAP(int, x, y) where x and y are variables of type int. */
140 #define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
142 /* Used to set up and access simple temp stacks within functions. */
143 #define TEMP(K) (temp + (K))
144 #define SETUP(E, C) \
145 do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
147 /* Compare value to zero. */
149 (((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
151 /* Multiply X by Y into Z, ignoring signs. Requires that Z have
152 enough storage preallocated to hold the result. */
153 #define UMUL(X, Y, Z) \
154 do{mp_size ua_=MP_USED(X),ub_=MP_USED(Y);mp_size o_=ua_+ub_;\
155 ZERO(MP_DIGITS(Z),o_);\
156 (void) s_kmul(MP_DIGITS(X),MP_DIGITS(Y),MP_DIGITS(Z),ua_,ub_);\
157 MP_USED(Z)=o_;CLAMP(Z);}while(0)
159 /* Square X into Z. Requires that Z have enough storage to hold the
162 do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
163 (void) s_ksqr(MP_DIGITS(X),MP_DIGITS(Z),ua_);MP_USED(Z)=o_;CLAMP(Z);}while(0)
165 #define UPPER_HALF(W) ((mp_word)((W) >> MP_DIGIT_BIT))
166 #define LOWER_HALF(W) ((mp_digit)(W))
167 #define HIGH_BIT_SET(W) ((W) >> (MP_WORD_BIT - 1))
168 #define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
171 /* {{{ Default configuration settings */
173 /* Default number of digits allocated to a new mp_int */
175 mp_size default_precision
= MP_DEFAULT_PREC
;
177 STATIC
const mp_size default_precision
= MP_DEFAULT_PREC
;
180 /* Minimum number of digits to invoke recursive multiply */
182 mp_size multiply_threshold
= MP_MULT_THRESH
;
184 STATIC
const mp_size multiply_threshold
= MP_MULT_THRESH
;
189 /* Allocate a buffer of (at least) num digits, or return
190 NULL if that couldn't be done. */
191 STATIC mp_digit
*s_alloc(mp_size num
);
193 /* Release a buffer of digits allocated by s_alloc(). */
194 STATIC
void s_free(void *ptr
);
196 /* Insure that z has at least min digits allocated, resizing if
197 necessary. Returns true if successful, false if out of memory. */
198 STATIC
int s_pad(mp_int z
, mp_size min
);
200 /* Fill in a "fake" mp_int on the stack with a given value */
201 STATIC
void s_fake(mp_int z
, mp_small value
, mp_digit vbuf
[]);
203 /* Compare two runs of digits of given length, returns <0, 0, >0 */
204 STATIC
int s_cdig(mp_digit
*da
, mp_digit
*db
, mp_size len
);
206 /* Pack the unsigned digits of v into array t */
207 STATIC
int s_vpack(mp_small v
, mp_digit t
[]);
209 /* Compare magnitudes of a and b, returns <0, 0, >0 */
210 STATIC
int s_ucmp(mp_int a
, mp_int b
);
212 /* Compare magnitudes of a and v, returns <0, 0, >0 */
213 STATIC
int s_vcmp(mp_int a
, mp_small v
);
215 /* Unsigned magnitude addition; assumes dc is big enough.
216 Carry out is returned (no memory allocated). */
217 STATIC mp_digit
s_uadd(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
218 mp_size size_a
, mp_size size_b
);
220 /* Unsigned magnitude subtraction. Assumes dc is big enough. */
221 STATIC
void s_usub(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
222 mp_size size_a
, mp_size size_b
);
224 /* Unsigned recursive multiplication. Assumes dc is big enough. */
225 STATIC
int s_kmul(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
226 mp_size size_a
, mp_size size_b
);
228 /* Unsigned magnitude multiplication. Assumes dc is big enough. */
229 STATIC
void s_umul(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
230 mp_size size_a
, mp_size size_b
);
232 /* Unsigned recursive squaring. Assumes dc is big enough. */
233 STATIC
int s_ksqr(mp_digit
*da
, mp_digit
*dc
, mp_size size_a
);
235 /* Unsigned magnitude squaring. Assumes dc is big enough. */
236 STATIC
void s_usqr(mp_digit
*da
, mp_digit
*dc
, mp_size size_a
);
238 /* Single digit addition. Assumes a is big enough. */
239 STATIC
void s_dadd(mp_int a
, mp_digit b
);
241 /* Single digit multiplication. Assumes a is big enough. */
242 STATIC
void s_dmul(mp_int a
, mp_digit b
);
244 /* Single digit multiplication on buffers; assumes dc is big enough. */
245 STATIC
void s_dbmul(mp_digit
*da
, mp_digit b
, mp_digit
*dc
,
248 /* Single digit division. Replaces a with the quotient,
249 returns the remainder. */
250 STATIC mp_digit
s_ddiv(mp_int a
, mp_digit b
);
252 /* Quick division by a power of 2, replaces z (no allocation) */
253 STATIC
void s_qdiv(mp_int z
, mp_size p2
);
255 /* Quick remainder by a power of 2, replaces z (no allocation) */
256 STATIC
void s_qmod(mp_int z
, mp_size p2
);
258 /* Quick multiplication by a power of 2, replaces z.
259 Allocates if necessary; returns false in case this fails. */
260 STATIC
int s_qmul(mp_int z
, mp_size p2
);
262 /* Quick subtraction from a power of 2, replaces z.
263 Allocates if necessary; returns false in case this fails. */
264 STATIC
int s_qsub(mp_int z
, mp_size p2
);
266 /* Return maximum k such that 2^k divides z. */
267 STATIC
int s_dp2k(mp_int z
);
269 /* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
270 STATIC
int s_isp2(mp_int z
);
272 /* Set z to 2^k. May allocate; returns false in case this fails. */
273 STATIC
int s_2expt(mp_int z
, mp_small k
);
275 /* Normalize a and b for division, returns normalization constant */
276 STATIC
int s_norm(mp_int a
, mp_int b
);
278 /* Compute constant mu for Barrett reduction, given modulus m, result
279 replaces z, m is untouched. */
280 STATIC mp_result
s_brmu(mp_int z
, mp_int m
);
282 /* Reduce a modulo m, using Barrett's algorithm. */
283 STATIC
int s_reduce(mp_int x
, mp_int m
, mp_int mu
, mp_int q1
, mp_int q2
);
285 /* Modular exponentiation, using Barrett reduction */
286 STATIC mp_result
s_embar(mp_int a
, mp_int b
, mp_int m
, mp_int mu
, mp_int c
);
288 /* Unsigned magnitude division. Assumes |a| > |b|. Allocates
289 temporaries; overwrites a with quotient, b with remainder. */
290 STATIC mp_result
s_udiv(mp_int a
, mp_int b
);
292 /* Compute the number of digits in radix r required to represent the
293 given value. Does not account for sign flags, terminators, etc. */
294 STATIC
int s_outlen(mp_int z
, mp_size r
);
296 /* Guess how many digits of precision will be needed to represent a
297 radix r value of the specified number of digits. Returns a value
298 guaranteed to be no smaller than the actual number required. */
299 STATIC mp_size
s_inlen(int len
, mp_size r
);
301 /* Convert a character to a digit value in radix r, or
302 -1 if out of range */
303 STATIC
int s_ch2val(char c
, int r
);
305 /* Convert a digit value to a character */
306 STATIC
char s_val2ch(int v
, int caps
);
308 /* Take 2's complement of a buffer in place */
309 STATIC
void s_2comp(unsigned char *buf
, int len
);
311 /* Convert a value to binary, ignoring sign. On input, *limpos is the
312 bound on how many bytes should be written to buf; on output, *limpos
313 is set to the number of bytes actually written. */
314 STATIC mp_result
s_tobin(mp_int z
, unsigned char *buf
, int *limpos
, int pad
);
317 /* Dump a representation of the mp_int to standard output */
318 void s_print(char *tag
, mp_int z
);
319 void s_print_buf(char *tag
, mp_digit
*buf
, mp_size num
);
322 /* {{{ mp_int_init(z) */
324 mp_result
mp_int_init(mp_int z
)
330 z
->digits
= &(z
->single
);
340 /* {{{ mp_int_alloc() */
342 mp_int
mp_int_alloc(void)
344 mp_int out
= malloc(sizeof(mpz_t
));
354 /* {{{ mp_int_init_size(z, prec) */
356 mp_result
mp_int_init_size(mp_int z
, mp_size prec
)
361 prec
= default_precision
;
363 return mp_int_init(z
);
365 prec
= (mp_size
) ROUND_PREC(prec
);
367 if((MP_DIGITS(z
) = s_alloc(prec
)) == NULL
)
373 MP_SIGN(z
) = MP_ZPOS
;
380 /* {{{ mp_int_init_copy(z, old) */
382 mp_result
mp_int_init_copy(mp_int z
, mp_int old
)
387 CHECK(z
!= NULL
&& old
!= NULL
);
394 mp_size target
= MAX(uold
, default_precision
);
396 if((res
= mp_int_init_size(z
, target
)) != MP_OK
)
401 MP_SIGN(z
) = MP_SIGN(old
);
402 COPY(MP_DIGITS(old
), MP_DIGITS(z
), uold
);
409 /* {{{ mp_int_init_value(z, value) */
411 mp_result
mp_int_init_value(mp_int z
, mp_small value
)
414 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
416 s_fake(&vtmp
, value
, vbuf
);
417 return mp_int_init_copy(z
, &vtmp
);
422 /* {{{ mp_int_set_value(z, value) */
424 mp_result
mp_int_set_value(mp_int z
, mp_small value
)
427 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
429 s_fake(&vtmp
, value
, vbuf
);
430 return mp_int_copy(&vtmp
, z
);
435 /* {{{ mp_int_clear(z) */
437 void mp_int_clear(mp_int z
)
442 if(MP_DIGITS(z
) != NULL
) {
443 if((void *) MP_DIGITS(z
) != (void *) z
)
444 s_free(MP_DIGITS(z
));
452 /* {{{ mp_int_free(z) */
454 void mp_int_free(mp_int z
)
459 free(z
); /* note: NOT s_free() */
464 /* {{{ mp_int_copy(a, c) */
466 mp_result
mp_int_copy(mp_int a
, mp_int c
)
468 CHECK(a
!= NULL
&& c
!= NULL
);
471 mp_size ua
= MP_USED(a
);
477 da
= MP_DIGITS(a
); dc
= MP_DIGITS(c
);
481 MP_SIGN(c
) = MP_SIGN(a
);
489 /* {{{ mp_int_swap(a, c) */
491 void mp_int_swap(mp_int a
, mp_int c
)
503 /* {{{ mp_int_zero(z) */
505 void mp_int_zero(mp_int z
)
511 MP_SIGN(z
) = MP_ZPOS
;
516 /* {{{ mp_int_abs(a, c) */
518 mp_result
mp_int_abs(mp_int a
, mp_int c
)
522 CHECK(a
!= NULL
&& c
!= NULL
);
524 if((res
= mp_int_copy(a
, c
)) != MP_OK
)
527 MP_SIGN(c
) = MP_ZPOS
;
533 /* {{{ mp_int_neg(a, c) */
535 mp_result
mp_int_neg(mp_int a
, mp_int c
)
539 CHECK(a
!= NULL
&& c
!= NULL
);
541 if((res
= mp_int_copy(a
, c
)) != MP_OK
)
545 MP_SIGN(c
) = 1 - MP_SIGN(a
);
552 /* {{{ mp_int_add(a, b, c) */
554 mp_result
mp_int_add(mp_int a
, mp_int b
, mp_int c
)
558 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
560 ua
= MP_USED(a
); ub
= MP_USED(b
);
563 if(MP_SIGN(a
) == MP_SIGN(b
)) {
564 /* Same sign -- add magnitudes, preserve sign of addends */
571 carry
= s_uadd(MP_DIGITS(a
), MP_DIGITS(b
), MP_DIGITS(c
), ua
, ub
);
575 if(!s_pad(c
, max
+ 1))
578 c
->digits
[max
] = carry
;
583 MP_SIGN(c
) = MP_SIGN(a
);
587 /* Different signs -- subtract magnitudes, preserve sign of greater */
589 int cmp
= s_ucmp(a
, b
); /* magnitude comparision, sign ignored */
591 /* Set x to max(a, b), y to min(a, b) to simplify later code.
592 A special case yields zero for equal magnitudes.
605 if(!s_pad(c
, MP_USED(x
)))
608 /* Subtract smaller from larger */
609 s_usub(MP_DIGITS(x
), MP_DIGITS(y
), MP_DIGITS(c
), MP_USED(x
), MP_USED(y
));
610 MP_USED(c
) = MP_USED(x
);
613 /* Give result the sign of the larger */
614 MP_SIGN(c
) = MP_SIGN(x
);
622 /* {{{ mp_int_add_value(a, value, c) */
624 mp_result
mp_int_add_value(mp_int a
, mp_small value
, mp_int c
)
627 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
629 s_fake(&vtmp
, value
, vbuf
);
631 return mp_int_add(a
, &vtmp
, c
);
636 /* {{{ mp_int_sub(a, b, c) */
638 mp_result
mp_int_sub(mp_int a
, mp_int b
, mp_int c
)
642 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
644 ua
= MP_USED(a
); ub
= MP_USED(b
);
647 if(MP_SIGN(a
) != MP_SIGN(b
)) {
648 /* Different signs -- add magnitudes and keep sign of a */
655 carry
= s_uadd(MP_DIGITS(a
), MP_DIGITS(b
), MP_DIGITS(c
), ua
, ub
);
659 if(!s_pad(c
, max
+ 1))
662 c
->digits
[max
] = carry
;
667 MP_SIGN(c
) = MP_SIGN(a
);
671 /* Same signs -- subtract magnitudes */
674 int cmp
= s_ucmp(a
, b
);
680 x
= a
; y
= b
; osign
= MP_ZPOS
;
683 x
= b
; y
= a
; osign
= MP_NEG
;
686 if(MP_SIGN(a
) == MP_NEG
&& cmp
!= 0)
689 s_usub(MP_DIGITS(x
), MP_DIGITS(y
), MP_DIGITS(c
), MP_USED(x
), MP_USED(y
));
690 MP_USED(c
) = MP_USED(x
);
701 /* {{{ mp_int_sub_value(a, value, c) */
703 mp_result
mp_int_sub_value(mp_int a
, mp_small value
, mp_int c
)
706 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
708 s_fake(&vtmp
, value
, vbuf
);
710 return mp_int_sub(a
, &vtmp
, c
);
715 /* {{{ mp_int_mul(a, b, c) */
717 mp_result
mp_int_mul(mp_int a
, mp_int b
, mp_int c
)
720 mp_size osize
, ua
, ub
, p
= 0;
723 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
725 /* If either input is zero, we can shortcut multiplication */
726 if(mp_int_compare_zero(a
) == 0 || mp_int_compare_zero(b
) == 0) {
731 /* Output is positive if inputs have same sign, otherwise negative */
732 osign
= (MP_SIGN(a
) == MP_SIGN(b
)) ? MP_ZPOS
: MP_NEG
;
734 /* If the output is not identical to any of the inputs, we'll write
735 the results directly; otherwise, allocate a temporary space. */
736 ua
= MP_USED(a
); ub
= MP_USED(b
);
738 osize
= 4 * ((osize
+ 1) / 2);
740 if(c
== a
|| c
== b
) {
741 p
= ROUND_PREC(osize
);
742 p
= MAX(p
, default_precision
);
744 if((out
= s_alloc(p
)) == NULL
)
755 if(!s_kmul(MP_DIGITS(a
), MP_DIGITS(b
), out
, ua
, ub
))
758 /* If we allocated a new buffer, get rid of whatever memory c was
759 already using, and fix up its fields to reflect that.
761 if(out
!= MP_DIGITS(c
)) {
762 if((void *) MP_DIGITS(c
) != (void *) c
)
763 s_free(MP_DIGITS(c
));
768 MP_USED(c
) = osize
; /* might not be true, but we'll fix it ... */
769 CLAMP(c
); /* ... right here */
777 /* {{{ mp_int_mul_value(a, value, c) */
779 mp_result
mp_int_mul_value(mp_int a
, mp_small value
, mp_int c
)
782 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
784 s_fake(&vtmp
, value
, vbuf
);
786 return mp_int_mul(a
, &vtmp
, c
);
791 /* {{{ mp_int_mul_pow2(a, p2, c) */
793 mp_result
mp_int_mul_pow2(mp_int a
, mp_small p2
, mp_int c
)
796 CHECK(a
!= NULL
&& c
!= NULL
&& p2
>= 0);
798 if((res
= mp_int_copy(a
, c
)) != MP_OK
)
801 if(s_qmul(c
, (mp_size
) p2
))
809 /* {{{ mp_int_sqr(a, c) */
811 mp_result
mp_int_sqr(mp_int a
, mp_int c
)
814 mp_size osize
, p
= 0;
816 CHECK(a
!= NULL
&& c
!= NULL
);
818 /* Get a temporary buffer big enough to hold the result */
819 osize
= (mp_size
) 4 * ((MP_USED(a
) + 1) / 2);
821 p
= ROUND_PREC(osize
);
822 p
= MAX(p
, default_precision
);
824 if((out
= s_alloc(p
)) == NULL
)
835 s_ksqr(MP_DIGITS(a
), out
, MP_USED(a
));
837 /* Get rid of whatever memory c was already using, and fix up its
838 fields to reflect the new digit array it's using
840 if(out
!= MP_DIGITS(c
)) {
841 if((void *) MP_DIGITS(c
) != (void *) c
)
842 s_free(MP_DIGITS(c
));
847 MP_USED(c
) = osize
; /* might not be true, but we'll fix it ... */
848 CLAMP(c
); /* ... right here */
849 MP_SIGN(c
) = MP_ZPOS
;
856 /* {{{ mp_int_div(a, b, q, r) */
858 mp_result
mp_int_div(mp_int a
, mp_int b
, mp_int q
, mp_int r
)
860 int cmp
, last
= 0, lg
;
861 mp_result res
= MP_OK
;
864 mp_sign sa
= MP_SIGN(a
), sb
= MP_SIGN(b
);
866 CHECK(a
!= NULL
&& b
!= NULL
&& q
!= r
);
870 else if((cmp
= s_ucmp(a
, b
)) < 0) {
871 /* If |a| < |b|, no division is required:
874 if(r
&& (res
= mp_int_copy(a
, r
)) != MP_OK
)
883 /* If |a| = |b|, no division is required:
900 /* When |a| > |b|, real division is required. We need someplace to
901 store quotient and remainder, but q and r are allowed to be NULL
902 or to overlap with the inputs.
904 if((lg
= s_isp2(b
)) < 0) {
906 if((res
= mp_int_copy(a
, q
)) != MP_OK
)
913 SETUP(mp_int_init_copy(TEMP(last
), a
), last
);
917 if((res
= mp_int_copy(b
, r
)) != MP_OK
)
924 SETUP(mp_int_init_copy(TEMP(last
), b
), last
);
927 if((res
= s_udiv(qout
, rout
)) != MP_OK
) goto CLEANUP
;
930 if(q
&& (res
= mp_int_copy(a
, q
)) != MP_OK
) goto CLEANUP
;
931 if(r
&& (res
= mp_int_copy(a
, r
)) != MP_OK
) goto CLEANUP
;
933 if(q
) s_qdiv(q
, (mp_size
) lg
); qout
= q
;
934 if(r
) s_qmod(r
, (mp_size
) lg
); rout
= r
;
937 /* Recompute signs for output */
941 MP_SIGN(rout
) = MP_ZPOS
;
944 MP_SIGN(qout
) = (sa
== sb
) ? MP_ZPOS
: MP_NEG
;
946 MP_SIGN(qout
) = MP_ZPOS
;
949 if(q
&& (res
= mp_int_copy(qout
, q
)) != MP_OK
) goto CLEANUP
;
950 if(r
&& (res
= mp_int_copy(rout
, r
)) != MP_OK
) goto CLEANUP
;
954 mp_int_clear(TEMP(last
));
961 /* {{{ mp_int_mod(a, m, c) */
963 mp_result
mp_int_mod(mp_int a
, mp_int m
, mp_int c
)
977 if((res
= mp_int_div(a
, m
, NULL
, out
)) != MP_OK
)
981 res
= mp_int_add(out
, m
, c
);
983 res
= mp_int_copy(out
, c
);
994 /* {{{ mp_int_div_value(a, value, q, r) */
996 mp_result
mp_int_div_value(mp_int a
, mp_small value
, mp_int q
, mp_small
*r
)
999 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
1003 s_fake(&vtmp
, value
, vbuf
);
1005 if((res
= mp_int_div(a
, &vtmp
, q
, &rtmp
)) != MP_OK
)
1009 (void) mp_int_to_int(&rtmp
, r
); /* can't fail */
1012 mp_int_clear(&rtmp
);
1018 /* {{{ mp_int_div_pow2(a, p2, q, r) */
1020 mp_result
mp_int_div_pow2(mp_int a
, mp_small p2
, mp_int q
, mp_int r
)
1022 mp_result res
= MP_OK
;
1024 CHECK(a
!= NULL
&& p2
>= 0 && q
!= r
);
1026 if(q
!= NULL
&& (res
= mp_int_copy(a
, q
)) == MP_OK
)
1027 s_qdiv(q
, (mp_size
) p2
);
1029 if(res
== MP_OK
&& r
!= NULL
&& (res
= mp_int_copy(a
, r
)) == MP_OK
)
1030 s_qmod(r
, (mp_size
) p2
);
1037 /* {{{ mp_int_expt(a, b, c) */
1039 mp_result
mp_int_expt(mp_int a
, mp_small b
, mp_int c
)
1043 unsigned int v
= abs(b
);
1045 CHECK(b
>= 0 && c
!= NULL
);
1047 if((res
= mp_int_init_copy(&t
, a
)) != MP_OK
)
1050 (void) mp_int_set_value(c
, 1);
1053 if((res
= mp_int_mul(c
, &t
, c
)) != MP_OK
)
1060 if((res
= mp_int_sqr(&t
, &t
)) != MP_OK
)
1071 /* {{{ mp_int_expt_value(a, b, c) */
1073 mp_result
mp_int_expt_value(mp_small a
, mp_small b
, mp_int c
)
1077 unsigned int v
= abs(b
);
1079 CHECK(b
>= 0 && c
!= NULL
);
1081 if((res
= mp_int_init_value(&t
, a
)) != MP_OK
)
1084 (void) mp_int_set_value(c
, 1);
1087 if((res
= mp_int_mul(c
, &t
, c
)) != MP_OK
)
1094 if((res
= mp_int_sqr(&t
, &t
)) != MP_OK
)
1105 /* {{{ mp_int_compare(a, b) */
1107 int mp_int_compare(mp_int a
, mp_int b
)
1111 CHECK(a
!= NULL
&& b
!= NULL
);
1114 if(sa
== MP_SIGN(b
)) {
1115 int cmp
= s_ucmp(a
, b
);
1117 /* If they're both zero or positive, the normal comparison
1118 applies; if both negative, the sense is reversed. */
1135 /* {{{ mp_int_compare_unsigned(a, b) */
1137 int mp_int_compare_unsigned(mp_int a
, mp_int b
)
1139 NRCHECK(a
!= NULL
&& b
!= NULL
);
1141 return s_ucmp(a
, b
);
1146 /* {{{ mp_int_compare_zero(z) */
1148 int mp_int_compare_zero(mp_int z
)
1152 if(MP_USED(z
) == 1 && z
->digits
[0] == 0)
1154 else if(MP_SIGN(z
) == MP_ZPOS
)
1162 /* {{{ mp_int_compare_value(z, value) */
1164 int mp_int_compare_value(mp_int z
, mp_small value
)
1166 mp_sign vsign
= (value
< 0) ? MP_NEG
: MP_ZPOS
;
1171 if(vsign
== MP_SIGN(z
)) {
1172 cmp
= s_vcmp(z
, value
);
1174 if(vsign
== MP_ZPOS
)
1189 /* {{{ mp_int_exptmod(a, b, m, c) */
1191 mp_result
mp_int_exptmod(mp_int a
, mp_int b
, mp_int m
, mp_int c
)
1199 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
&& m
!= NULL
);
1201 /* Zero moduli and negative exponents are not considered. */
1208 SETUP(mp_int_init_size(TEMP(0), 2 * um
), last
);
1209 SETUP(mp_int_init_size(TEMP(1), 2 * um
), last
);
1211 if(c
== b
|| c
== m
) {
1212 SETUP(mp_int_init_size(TEMP(2), 2 * um
), last
);
1219 if((res
= mp_int_mod(a
, m
, TEMP(0))) != MP_OK
) goto CLEANUP
;
1221 if((res
= s_brmu(TEMP(1), m
)) != MP_OK
) goto CLEANUP
;
1223 if((res
= s_embar(TEMP(0), b
, m
, TEMP(1), s
)) != MP_OK
)
1226 res
= mp_int_copy(s
, c
);
1230 mp_int_clear(TEMP(last
));
1237 /* {{{ mp_int_exptmod_evalue(a, value, m, c) */
1239 mp_result
mp_int_exptmod_evalue(mp_int a
, mp_small value
, mp_int m
, mp_int c
)
1242 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
1244 s_fake(&vtmp
, value
, vbuf
);
1246 return mp_int_exptmod(a
, &vtmp
, m
, c
);
1251 /* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
1253 mp_result
mp_int_exptmod_bvalue(mp_small value
, mp_int b
,
1257 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
1259 s_fake(&vtmp
, value
, vbuf
);
1261 return mp_int_exptmod(&vtmp
, b
, m
, c
);
1266 /* {{{ mp_int_exptmod_known(a, b, m, mu, c) */
1268 mp_result
mp_int_exptmod_known(mp_int a
, mp_int b
, mp_int m
, mp_int mu
, mp_int c
)
1276 CHECK(a
&& b
&& m
&& c
);
1278 /* Zero moduli and negative exponents are not considered. */
1285 SETUP(mp_int_init_size(TEMP(0), 2 * um
), last
);
1287 if(c
== b
|| c
== m
) {
1288 SETUP(mp_int_init_size(TEMP(1), 2 * um
), last
);
1295 if((res
= mp_int_mod(a
, m
, TEMP(0))) != MP_OK
) goto CLEANUP
;
1297 if((res
= s_embar(TEMP(0), b
, m
, mu
, s
)) != MP_OK
)
1300 res
= mp_int_copy(s
, c
);
1304 mp_int_clear(TEMP(last
));
1311 /* {{{ mp_int_redux_const(m, c) */
1313 mp_result
mp_int_redux_const(mp_int m
, mp_int c
)
1315 CHECK(m
!= NULL
&& c
!= NULL
&& m
!= c
);
1317 return s_brmu(c
, m
);
1322 /* {{{ mp_int_invmod(a, m, c) */
1324 mp_result
mp_int_invmod(mp_int a
, mp_int m
, mp_int c
)
1331 CHECK(a
!= NULL
&& m
!= NULL
&& c
!= NULL
);
1333 if(CMPZ(a
) == 0 || CMPZ(m
) <= 0)
1336 sa
= MP_SIGN(a
); /* need this for the result later */
1338 for(last
= 0; last
< 2; ++last
)
1339 mp_int_init(TEMP(last
));
1341 if((res
= mp_int_egcd(a
, m
, TEMP(0), TEMP(1), NULL
)) != MP_OK
)
1344 if(mp_int_compare_value(TEMP(0), 1) != 0) {
1349 /* It is first necessary to constrain the value to the proper range */
1350 if((res
= mp_int_mod(TEMP(1), m
, TEMP(1))) != MP_OK
)
1353 /* Now, if 'a' was originally negative, the value we have is
1354 actually the magnitude of the negative representative; to get the
1355 positive value we have to subtract from the modulus. Otherwise,
1356 the value is okay as it stands.
1359 res
= mp_int_sub(m
, TEMP(1), c
);
1361 res
= mp_int_copy(TEMP(1), c
);
1365 mp_int_clear(TEMP(last
));
1372 /* {{{ mp_int_gcd(a, b, c) */
1374 /* Binary GCD algorithm due to Josef Stein, 1961 */
1375 mp_result
mp_int_gcd(mp_int a
, mp_int b
, mp_int c
)
1381 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
1385 if(ca
== 0 && cb
== 0)
1388 return mp_int_abs(b
, c
);
1390 return mp_int_abs(a
, c
);
1393 if((res
= mp_int_init_copy(&u
, a
)) != MP_OK
)
1395 if((res
= mp_int_init_copy(&v
, b
)) != MP_OK
)
1398 MP_SIGN(&u
) = MP_ZPOS
; MP_SIGN(&v
) = MP_ZPOS
;
1400 { /* Divide out common factors of 2 from u and v */
1401 int div2_u
= s_dp2k(&u
), div2_v
= s_dp2k(&v
);
1403 k
= MIN(div2_u
, div2_v
);
1404 s_qdiv(&u
, (mp_size
) k
);
1405 s_qdiv(&v
, (mp_size
) k
);
1408 if(mp_int_is_odd(&u
)) {
1409 if((res
= mp_int_neg(&v
, &t
)) != MP_OK
)
1413 if((res
= mp_int_copy(&u
, &t
)) != MP_OK
)
1418 s_qdiv(&t
, s_dp2k(&t
));
1421 if((res
= mp_int_copy(&t
, &u
)) != MP_OK
)
1425 if((res
= mp_int_neg(&t
, &v
)) != MP_OK
)
1429 if((res
= mp_int_sub(&u
, &v
, &t
)) != MP_OK
)
1436 if((res
= mp_int_abs(&u
, c
)) != MP_OK
)
1438 if(!s_qmul(c
, (mp_size
) k
))
1443 V
: mp_int_clear(&u
);
1444 U
: mp_int_clear(&t
);
1451 /* {{{ mp_int_egcd(a, b, c, x, y) */
1453 /* This is the binary GCD algorithm again, but this time we keep track
1454 of the elementary matrix operations as we go, so we can get values
1455 x and y satisfying c = ax + by.
1457 mp_result
mp_int_egcd(mp_int a
, mp_int b
, mp_int c
,
1460 int k
, last
= 0, ca
, cb
;
1464 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
&&
1465 (x
!= NULL
|| y
!= NULL
));
1469 if(ca
== 0 && cb
== 0)
1472 if((res
= mp_int_abs(b
, c
)) != MP_OK
) return res
;
1473 mp_int_zero(x
); (void) mp_int_set_value(y
, 1); return MP_OK
;
1476 if((res
= mp_int_abs(a
, c
)) != MP_OK
) return res
;
1477 (void) mp_int_set_value(x
, 1); mp_int_zero(y
); return MP_OK
;
1480 /* Initialize temporaries:
1481 A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */
1482 for(last
= 0; last
< 4; ++last
)
1483 mp_int_init(TEMP(last
));
1484 TEMP(0)->digits
[0] = 1;
1485 TEMP(3)->digits
[0] = 1;
1487 SETUP(mp_int_init_copy(TEMP(4), a
), last
);
1488 SETUP(mp_int_init_copy(TEMP(5), b
), last
);
1490 /* We will work with absolute values here */
1491 MP_SIGN(TEMP(4)) = MP_ZPOS
;
1492 MP_SIGN(TEMP(5)) = MP_ZPOS
;
1494 { /* Divide out common factors of 2 from u and v */
1495 int div2_u
= s_dp2k(TEMP(4)), div2_v
= s_dp2k(TEMP(5));
1497 k
= MIN(div2_u
, div2_v
);
1502 SETUP(mp_int_init_copy(TEMP(6), TEMP(4)), last
);
1503 SETUP(mp_int_init_copy(TEMP(7), TEMP(5)), last
);
1506 while(mp_int_is_even(TEMP(4))) {
1509 if(mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1))) {
1510 if((res
= mp_int_add(TEMP(0), TEMP(7), TEMP(0))) != MP_OK
)
1512 if((res
= mp_int_sub(TEMP(1), TEMP(6), TEMP(1))) != MP_OK
)
1520 while(mp_int_is_even(TEMP(5))) {
1523 if(mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3))) {
1524 if((res
= mp_int_add(TEMP(2), TEMP(7), TEMP(2))) != MP_OK
)
1526 if((res
= mp_int_sub(TEMP(3), TEMP(6), TEMP(3))) != MP_OK
)
1534 if(mp_int_compare(TEMP(4), TEMP(5)) >= 0) {
1535 if((res
= mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != MP_OK
) goto CLEANUP
;
1536 if((res
= mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != MP_OK
) goto CLEANUP
;
1537 if((res
= mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != MP_OK
) goto CLEANUP
;
1540 if((res
= mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != MP_OK
) goto CLEANUP
;
1541 if((res
= mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK
) goto CLEANUP
;
1542 if((res
= mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != MP_OK
) goto CLEANUP
;
1545 if(CMPZ(TEMP(4)) == 0) {
1546 if(x
&& (res
= mp_int_copy(TEMP(2), x
)) != MP_OK
) goto CLEANUP
;
1547 if(y
&& (res
= mp_int_copy(TEMP(3), y
)) != MP_OK
) goto CLEANUP
;
1549 if(!s_qmul(TEMP(5), k
)) {
1554 res
= mp_int_copy(TEMP(5), c
);
1563 mp_int_clear(TEMP(last
));
1570 /* {{{ mp_int_lcm(a, b, c) */
1572 mp_result
mp_int_lcm(mp_int a
, mp_int b
, mp_int c
)
1577 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
1579 /* Since a * b = gcd(a, b) * lcm(a, b), we can compute
1580 lcm(a, b) = (a / gcd(a, b)) * b.
1582 This formulation insures everything works even if the input
1583 variables share space.
1585 if((res
= mp_int_init(&lcm
)) != MP_OK
)
1587 if((res
= mp_int_gcd(a
, b
, &lcm
)) != MP_OK
)
1589 if((res
= mp_int_div(a
, &lcm
, &lcm
, NULL
)) != MP_OK
)
1591 if((res
= mp_int_mul(&lcm
, b
, &lcm
)) != MP_OK
)
1594 res
= mp_int_copy(&lcm
, c
);
1604 /* {{{ mp_int_divisible_value(a, v) */
1606 int mp_int_divisible_value(mp_int a
, mp_small v
)
1610 if(mp_int_div_value(a
, v
, NULL
, &rem
) != MP_OK
)
1618 /* {{{ mp_int_is_pow2(z) */
1620 int mp_int_is_pow2(mp_int z
)
1629 /* {{{ mp_int_root(a, b, c) */
1631 /* Implementation of Newton's root finding method, based loosely on a
1632 patch contributed by Hal Finkel <half@halssoftware.com>
1633 modified by M. J. Fromberger.
1635 mp_result
mp_int_root(mp_int a
, mp_small b
, mp_int c
)
1637 mp_result res
= MP_OK
;
1642 CHECK(a
!= NULL
&& c
!= NULL
&& b
> 0);
1645 return mp_int_copy(a
, c
);
1647 if(MP_SIGN(a
) == MP_NEG
) {
1649 return MP_UNDEF
; /* root does not exist for negative a with even b */
1654 SETUP(mp_int_init_copy(TEMP(last
), a
), last
);
1655 SETUP(mp_int_init_copy(TEMP(last
), a
), last
);
1656 SETUP(mp_int_init(TEMP(last
)), last
);
1657 SETUP(mp_int_init(TEMP(last
)), last
);
1658 SETUP(mp_int_init(TEMP(last
)), last
);
1660 (void) mp_int_abs(TEMP(0), TEMP(0));
1661 (void) mp_int_abs(TEMP(1), TEMP(1));
1664 if((res
= mp_int_expt(TEMP(1), b
, TEMP(2))) != MP_OK
)
1667 if(mp_int_compare_unsigned(TEMP(2), TEMP(0)) <= 0)
1670 if((res
= mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK
)
1672 if((res
= mp_int_expt(TEMP(1), b
- 1, TEMP(3))) != MP_OK
)
1674 if((res
= mp_int_mul_value(TEMP(3), b
, TEMP(3))) != MP_OK
)
1676 if((res
= mp_int_div(TEMP(2), TEMP(3), TEMP(4), NULL
)) != MP_OK
)
1678 if((res
= mp_int_sub(TEMP(1), TEMP(4), TEMP(4))) != MP_OK
)
1681 if(mp_int_compare_unsigned(TEMP(1), TEMP(4)) == 0) {
1682 if((res
= mp_int_sub_value(TEMP(4), 1, TEMP(4))) != MP_OK
)
1685 if((res
= mp_int_copy(TEMP(4), TEMP(1))) != MP_OK
)
1689 if((res
= mp_int_copy(TEMP(1), c
)) != MP_OK
)
1692 /* If the original value of a was negative, flip the output sign. */
1694 (void) mp_int_neg(c
, c
); /* cannot fail */
1698 mp_int_clear(TEMP(last
));
1705 /* {{{ mp_int_to_int(z, out) */
1707 mp_result
mp_int_to_int(mp_int z
, mp_small
*out
)
1716 /* Make sure the value is representable as an int */
1718 if((sz
== MP_ZPOS
&& mp_int_compare_value(z
, MP_SMALL_MAX
) > 0) ||
1719 mp_int_compare_value(z
, MP_SMALL_MIN
) < 0)
1723 dz
= MP_DIGITS(z
) + uz
- 1;
1726 uv
<<= MP_DIGIT_BIT
/2;
1727 uv
= (uv
<< (MP_DIGIT_BIT
/2)) | *dz
--;
1732 *out
= (sz
== MP_NEG
) ? -(mp_small
)uv
: (mp_small
)uv
;
1739 /* {{{ mp_int_to_uint(z, *out) */
1741 mp_result
mp_int_to_uint(mp_int z
, mp_usmall
*out
)
1750 /* Make sure the value is representable as an int */
1752 if(!(sz
== MP_ZPOS
&& mp_int_compare_value(z
, UINT_MAX
) <= 0))
1756 dz
= MP_DIGITS(z
) + uz
- 1;
1759 uv
<<= MP_DIGIT_BIT
/2;
1760 uv
= (uv
<< (MP_DIGIT_BIT
/2)) | *dz
--;
1772 /* {{{ mp_int_to_string(z, radix, str, limit) */
1774 mp_result
mp_int_to_string(mp_int z
, mp_size radix
,
1775 char *str
, int limit
)
1780 CHECK(z
!= NULL
&& str
!= NULL
&& limit
>= 2);
1782 if(radix
< MP_MIN_RADIX
|| radix
> MP_MAX_RADIX
)
1786 *str
++ = s_val2ch(0, 1);
1792 if((res
= mp_int_init_copy(&tmp
, z
)) != MP_OK
)
1795 if(MP_SIGN(z
) == MP_NEG
) {
1801 /* Generate digits in reverse order until finished or limit reached */
1802 for(/* */; limit
> 0; --limit
) {
1805 if((cmp
= CMPZ(&tmp
)) == 0)
1808 d
= s_ddiv(&tmp
, (mp_digit
)radix
);
1809 *str
++ = s_val2ch(d
, 1);
1813 /* Put digits back in correct output order */
1832 /* {{{ mp_int_string_len(z, radix) */
1834 mp_result
mp_int_string_len(mp_int z
, mp_size radix
)
1840 if(radix
< MP_MIN_RADIX
|| radix
> MP_MAX_RADIX
)
1843 len
= s_outlen(z
, radix
) + 1; /* for terminator */
1845 /* Allow for sign marker on negatives */
1846 if(MP_SIGN(z
) == MP_NEG
)
1854 /* {{{ mp_int_read_string(z, radix, *str) */
1856 /* Read zero-terminated string into z */
1857 mp_result
mp_int_read_string(mp_int z
, mp_size radix
, const char *str
)
1859 return mp_int_read_cstring(z
, radix
, str
, NULL
);
1865 /* {{{ mp_int_read_cstring(z, radix, *str, **end) */
1867 mp_result
mp_int_read_cstring(mp_int z
, mp_size radix
, const char *str
, char **end
)
1871 CHECK(z
!= NULL
&& str
!= NULL
);
1873 if(radix
< MP_MIN_RADIX
|| radix
> MP_MAX_RADIX
)
1876 /* Skip leading whitespace */
1877 while(isspace((int)*str
))
1880 /* Handle leading sign tag (+/-, positive default) */
1883 MP_SIGN(z
) = MP_NEG
;
1887 ++str
; /* fallthrough */
1889 MP_SIGN(z
) = MP_ZPOS
;
1893 /* Skip leading zeroes */
1894 while((ch
= s_ch2val(*str
, radix
)) == 0)
1897 /* Make sure there is enough space for the value */
1898 if(!s_pad(z
, s_inlen(strlen(str
), radix
)))
1901 MP_USED(z
) = 1; z
->digits
[0] = 0;
1903 while(*str
!= '\0' && ((ch
= s_ch2val(*str
, radix
)) >= 0)) {
1904 s_dmul(z
, (mp_digit
)radix
);
1905 s_dadd(z
, (mp_digit
)ch
);
1911 /* Override sign for zero, even if negative specified. */
1913 MP_SIGN(z
) = MP_ZPOS
;
1918 /* Return a truncation error if the string has unprocessed
1919 characters remaining, so the caller can tell if the whole string
1929 /* {{{ mp_int_count_bits(z) */
1931 mp_result
mp_int_count_bits(mp_int z
)
1933 mp_size nbits
= 0, uz
;
1939 if(uz
== 1 && z
->digits
[0] == 0)
1943 nbits
= uz
* MP_DIGIT_BIT
;
1956 /* {{{ mp_int_to_binary(z, buf, limit) */
1958 mp_result
mp_int_to_binary(mp_int z
, unsigned char *buf
, int limit
)
1960 static const int PAD_FOR_2C
= 1;
1965 CHECK(z
!= NULL
&& buf
!= NULL
);
1967 res
= s_tobin(z
, buf
, &limpos
, PAD_FOR_2C
);
1969 if(MP_SIGN(z
) == MP_NEG
)
1970 s_2comp(buf
, limpos
);
1977 /* {{{ mp_int_read_binary(z, buf, len) */
1979 mp_result
mp_int_read_binary(mp_int z
, unsigned char *buf
, int len
)
1985 CHECK(z
!= NULL
&& buf
!= NULL
&& len
> 0);
1987 /* Figure out how many digits are needed to represent this value */
1988 need
= ((len
* CHAR_BIT
) + (MP_DIGIT_BIT
- 1)) / MP_DIGIT_BIT
;
1994 /* If the high-order bit is set, take the 2's complement before
1995 reading the value (it will be restored afterward) */
1996 if(buf
[0] >> (CHAR_BIT
- 1)) {
1997 MP_SIGN(z
) = MP_NEG
;
2002 for(tmp
= buf
, i
= len
; i
> 0; --i
, ++tmp
) {
2003 s_qmul(z
, (mp_size
) CHAR_BIT
);
2007 /* Restore 2's complement if we took it before */
2008 if(MP_SIGN(z
) == MP_NEG
)
2016 /* {{{ mp_int_binary_len(z) */
2018 mp_result
mp_int_binary_len(mp_int z
)
2020 mp_result res
= mp_int_count_bits(z
);
2026 bytes
= (res
+ (CHAR_BIT
- 1)) / CHAR_BIT
;
2028 /* If the highest-order bit falls exactly on a byte boundary, we
2029 need to pad with an extra byte so that the sign will be read
2030 correctly when reading it back in. */
2031 if(bytes
* CHAR_BIT
== res
)
2039 /* {{{ mp_int_to_unsigned(z, buf, limit) */
2041 mp_result
mp_int_to_unsigned(mp_int z
, unsigned char *buf
, int limit
)
2043 static const int NO_PADDING
= 0;
2045 CHECK(z
!= NULL
&& buf
!= NULL
);
2047 return s_tobin(z
, buf
, &limit
, NO_PADDING
);
2052 /* {{{ mp_int_read_unsigned(z, buf, len) */
2054 mp_result
mp_int_read_unsigned(mp_int z
, unsigned char *buf
, int len
)
2060 CHECK(z
!= NULL
&& buf
!= NULL
&& len
> 0);
2062 /* Figure out how many digits are needed to represent this value */
2063 need
= ((len
* CHAR_BIT
) + (MP_DIGIT_BIT
- 1)) / MP_DIGIT_BIT
;
2070 for(tmp
= buf
, i
= len
; i
> 0; --i
, ++tmp
) {
2071 (void) s_qmul(z
, CHAR_BIT
);
2080 /* {{{ mp_int_unsigned_len(z) */
2082 mp_result
mp_int_unsigned_len(mp_int z
)
2084 mp_result res
= mp_int_count_bits(z
);
2090 bytes
= (res
+ (CHAR_BIT
- 1)) / CHAR_BIT
;
2097 /* {{{ mp_error_string(res) */
2099 const char *mp_error_string(mp_result res
)
2103 return s_unknown_err
;
2106 for(ix
= 0; ix
< res
&& s_error_msg
[ix
] != NULL
; ++ix
)
2109 if(s_error_msg
[ix
] != NULL
)
2110 return s_error_msg
[ix
];
2112 return s_unknown_err
;
2117 /*------------------------------------------------------------------------*/
2118 /* Private functions for internal use. These make assumptions. */
2120 /* {{{ s_alloc(num) */
2122 STATIC mp_digit
*s_alloc(mp_size num
)
2124 mp_digit
*out
= malloc(num
* sizeof(mp_digit
));
2126 assert(out
!= NULL
); /* for debugging */
2129 mp_digit v
= (mp_digit
) 0xdeadbeef;
2132 for(ix
= 0; ix
< num
; ++ix
)
2142 /* {{{ s_realloc(old, osize, nsize) */
2144 STATIC mp_digit
*s_realloc(mp_digit
*old
, mp_size osize
, mp_size nsize
)
2147 mp_digit
*new = s_alloc(nsize
);
2150 for(ix
= 0; ix
< nsize
; ++ix
)
2151 new[ix
] = (mp_digit
) 0xdeadbeef;
2153 memcpy(new, old
, osize
* sizeof(mp_digit
));
2155 mp_digit
*new = realloc(old
, nsize
* sizeof(mp_digit
));
2157 assert(new != NULL
); /* for debugging */
2164 /* {{{ s_free(ptr) */
2166 STATIC
void s_free(void *ptr
)
2173 /* {{{ s_pad(z, min) */
2175 STATIC
int s_pad(mp_int z
, mp_size min
)
2177 if(MP_ALLOC(z
) < min
) {
2178 mp_size nsize
= ROUND_PREC(min
);
2181 if((void *)z
->digits
== (void *)z
) {
2182 if((tmp
= s_alloc(nsize
)) == NULL
)
2185 COPY(MP_DIGITS(z
), tmp
, MP_USED(z
));
2187 else if((tmp
= s_realloc(MP_DIGITS(z
), MP_ALLOC(z
), nsize
)) == NULL
)
2191 MP_ALLOC(z
) = nsize
;
2199 /* {{{ s_fake(z, value, vbuf) */
2201 STATIC
void s_fake(mp_int z
, mp_small value
, mp_digit vbuf
[])
2203 mp_size uv
= (mp_size
) s_vpack(value
, vbuf
);
2206 z
->alloc
= MP_VALUE_DIGITS(value
);
2207 z
->sign
= (value
< 0) ? MP_NEG
: MP_ZPOS
;
2213 /* {{{ s_cdig(da, db, len) */
2215 STATIC
int s_cdig(mp_digit
*da
, mp_digit
*db
, mp_size len
)
2217 mp_digit
*dat
= da
+ len
- 1, *dbt
= db
+ len
- 1;
2219 for(/* */; len
!= 0; --len
, --dat
, --dbt
) {
2222 else if(*dat
< *dbt
)
2231 /* {{{ s_vpack(v, t[]) */
2233 STATIC
int s_vpack(mp_small v
, mp_digit t
[])
2235 mp_usmall uv
= (mp_usmall
) ((v
< 0) ? -v
: v
);
2242 t
[ndig
++] = (mp_digit
) uv
;
2243 uv
>>= MP_DIGIT_BIT
/2;
2244 uv
>>= MP_DIGIT_BIT
/2;
2253 /* {{{ s_ucmp(a, b) */
2255 STATIC
int s_ucmp(mp_int a
, mp_int b
)
2257 mp_size ua
= MP_USED(a
), ub
= MP_USED(b
);
2264 return s_cdig(MP_DIGITS(a
), MP_DIGITS(b
), ua
);
2269 /* {{{ s_vcmp(a, v) */
2271 STATIC
int s_vcmp(mp_int a
, mp_small v
)
2273 mp_digit vdig
[MP_VALUE_DIGITS(v
)];
2275 mp_size ua
= MP_USED(a
);
2277 ndig
= s_vpack(v
, vdig
);
2284 return s_cdig(MP_DIGITS(a
), vdig
, ndig
);
2289 /* {{{ s_uadd(da, db, dc, size_a, size_b) */
2291 STATIC mp_digit
s_uadd(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
2292 mp_size size_a
, mp_size size_b
)
2297 /* Insure that da is the longer of the two to simplify later code */
2298 if(size_b
> size_a
) {
2299 SWAP(mp_digit
*, da
, db
);
2300 SWAP(mp_size
, size_a
, size_b
);
2303 /* Add corresponding digits until the shorter number runs out */
2304 for(pos
= 0; pos
< size_b
; ++pos
, ++da
, ++db
, ++dc
) {
2305 w
= w
+ (mp_word
) *da
+ (mp_word
) *db
;
2306 *dc
= LOWER_HALF(w
);
2310 /* Propagate carries as far as necessary */
2311 for(/* */; pos
< size_a
; ++pos
, ++da
, ++dc
) {
2314 *dc
= LOWER_HALF(w
);
2318 /* Return carry out */
2324 /* {{{ s_usub(da, db, dc, size_a, size_b) */
2326 STATIC
void s_usub(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
2327 mp_size size_a
, mp_size size_b
)
2332 /* We assume that |a| >= |b| so this should definitely hold */
2333 assert(size_a
>= size_b
);
2335 /* Subtract corresponding digits and propagate borrow */
2336 for(pos
= 0; pos
< size_b
; ++pos
, ++da
, ++db
, ++dc
) {
2337 w
= ((mp_word
)MP_DIGIT_MAX
+ 1 + /* MP_RADIX */
2338 (mp_word
)*da
) - w
- (mp_word
)*db
;
2340 *dc
= LOWER_HALF(w
);
2341 w
= (UPPER_HALF(w
) == 0);
2344 /* Finish the subtraction for remaining upper digits of da */
2345 for(/* */; pos
< size_a
; ++pos
, ++da
, ++dc
) {
2346 w
= ((mp_word
)MP_DIGIT_MAX
+ 1 + /* MP_RADIX */
2349 *dc
= LOWER_HALF(w
);
2350 w
= (UPPER_HALF(w
) == 0);
2353 /* If there is a borrow out at the end, it violates the precondition */
2359 /* {{{ s_kmul(da, db, dc, size_a, size_b) */
2361 STATIC
int s_kmul(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
2362 mp_size size_a
, mp_size size_b
)
2366 /* Make sure b is the smaller of the two input values */
2367 if(size_b
> size_a
) {
2368 SWAP(mp_digit
*, da
, db
);
2369 SWAP(mp_size
, size_a
, size_b
);
2372 /* Insure that the bottom is the larger half in an odd-length split;
2373 the code below relies on this being true.
2375 bot_size
= (size_a
+ 1) / 2;
2377 /* If the values are big enough to bother with recursion, use the
2378 Karatsuba algorithm to compute the product; otherwise use the
2379 normal multiplication algorithm
2381 if(multiply_threshold
&&
2382 size_a
>= multiply_threshold
&&
2383 size_b
> bot_size
) {
2385 mp_digit
*t1
, *t2
, *t3
, carry
;
2387 mp_digit
*a_top
= da
+ bot_size
;
2388 mp_digit
*b_top
= db
+ bot_size
;
2390 mp_size at_size
= size_a
- bot_size
;
2391 mp_size bt_size
= size_b
- bot_size
;
2392 mp_size buf_size
= 2 * bot_size
;
2394 /* Do a single allocation for all three temporary buffers needed;
2395 each buffer must be big enough to hold the product of two
2396 bottom halves, and one buffer needs space for the completed
2397 product; twice the space is plenty.
2399 if((t1
= s_alloc(4 * buf_size
)) == NULL
) return 0;
2402 ZERO(t1
, 4 * buf_size
);
2404 /* t1 and t2 are initially used as temporaries to compute the inner product
2405 (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
2407 carry
= s_uadd(da
, a_top
, t1
, bot_size
, at_size
); /* t1 = a1 + a0 */
2408 t1
[bot_size
] = carry
;
2410 carry
= s_uadd(db
, b_top
, t2
, bot_size
, bt_size
); /* t2 = b1 + b0 */
2411 t2
[bot_size
] = carry
;
2413 (void) s_kmul(t1
, t2
, t3
, bot_size
+ 1, bot_size
+ 1); /* t3 = t1 * t2 */
2415 /* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that
2416 we're left with only the pieces we want: t3 = a1b0 + a0b1
2420 (void) s_kmul(da
, db
, t1
, bot_size
, bot_size
); /* t1 = a0 * b0 */
2421 (void) s_kmul(a_top
, b_top
, t2
, at_size
, bt_size
); /* t2 = a1 * b1 */
2423 /* Subtract out t1 and t2 to get the inner product */
2424 s_usub(t3
, t1
, t3
, buf_size
+ 2, buf_size
);
2425 s_usub(t3
, t2
, t3
, buf_size
+ 2, buf_size
);
2427 /* Assemble the output value */
2428 COPY(t1
, dc
, buf_size
);
2429 carry
= s_uadd(t3
, dc
+ bot_size
, dc
+ bot_size
,
2430 buf_size
+ 1, buf_size
);
2433 carry
= s_uadd(t2
, dc
+ 2*bot_size
, dc
+ 2*bot_size
,
2434 buf_size
, buf_size
);
2437 s_free(t1
); /* note t2 and t3 are just internal pointers to t1 */
2440 s_umul(da
, db
, dc
, size_a
, size_b
);
2448 /* {{{ s_umul(da, db, dc, size_a, size_b) */
2450 STATIC
void s_umul(mp_digit
*da
, mp_digit
*db
, mp_digit
*dc
,
2451 mp_size size_a
, mp_size size_b
)
2456 for(a
= 0; a
< size_a
; ++a
, ++dc
, ++da
) {
2464 for(b
= 0; b
< size_b
; ++b
, ++dbt
, ++dct
) {
2465 w
= (mp_word
)*da
* (mp_word
)*dbt
+ w
+ (mp_word
)*dct
;
2467 *dct
= LOWER_HALF(w
);
2477 /* {{{ s_ksqr(da, dc, size_a) */
2479 STATIC
int s_ksqr(mp_digit
*da
, mp_digit
*dc
, mp_size size_a
)
2481 if(multiply_threshold
&& size_a
> multiply_threshold
) {
2482 mp_size bot_size
= (size_a
+ 1) / 2;
2483 mp_digit
*a_top
= da
+ bot_size
;
2484 mp_digit
*t1
, *t2
, *t3
, carry
;
2485 mp_size at_size
= size_a
- bot_size
;
2486 mp_size buf_size
= 2 * bot_size
;
2488 if((t1
= s_alloc(4 * buf_size
)) == NULL
) return 0;
2491 ZERO(t1
, 4 * buf_size
);
2493 (void) s_ksqr(da
, t1
, bot_size
); /* t1 = a0 ^ 2 */
2494 (void) s_ksqr(a_top
, t2
, at_size
); /* t2 = a1 ^ 2 */
2496 (void) s_kmul(da
, a_top
, t3
, bot_size
, at_size
); /* t3 = a0 * a1 */
2498 /* Quick multiply t3 by 2, shifting left (can't overflow) */
2500 int i
, top
= bot_size
+ at_size
;
2501 mp_word w
, save
= 0;
2503 for(i
= 0; i
< top
; ++i
) {
2505 w
= (w
<< 1) | save
;
2506 t3
[i
] = LOWER_HALF(w
);
2507 save
= UPPER_HALF(w
);
2509 t3
[i
] = LOWER_HALF(save
);
2512 /* Assemble the output value */
2513 COPY(t1
, dc
, 2 * bot_size
);
2514 carry
= s_uadd(t3
, dc
+ bot_size
, dc
+ bot_size
,
2515 buf_size
+ 1, buf_size
);
2518 carry
= s_uadd(t2
, dc
+ 2*bot_size
, dc
+ 2*bot_size
,
2519 buf_size
, buf_size
);
2522 s_free(t1
); /* note that t2 and t2 are internal pointers only */
2526 s_usqr(da
, dc
, size_a
);
2534 /* {{{ s_usqr(da, dc, size_a) */
2536 STATIC
void s_usqr(mp_digit
*da
, mp_digit
*dc
, mp_size size_a
)
2541 for(i
= 0; i
< size_a
; ++i
, dc
+= 2, ++da
) {
2542 mp_digit
*dct
= dc
, *dat
= da
;
2547 /* Take care of the first digit, no rollover */
2548 w
= (mp_word
)*dat
* (mp_word
)*dat
+ (mp_word
)*dct
;
2549 *dct
= LOWER_HALF(w
);
2553 for(j
= i
+ 1; j
< size_a
; ++j
, ++dat
, ++dct
) {
2554 mp_word t
= (mp_word
)*da
* (mp_word
)*dat
;
2555 mp_word u
= w
+ (mp_word
)*dct
, ov
= 0;
2557 /* Check if doubling t will overflow a word */
2563 /* Check if adding u to w will overflow a word */
2564 if(ADD_WILL_OVERFLOW(w
, u
))
2569 *dct
= LOWER_HALF(w
);
2572 w
+= MP_DIGIT_MAX
; /* MP_RADIX */
2579 while((w
= UPPER_HALF(w
)) != 0) {
2580 ++dct
; w
= w
+ *dct
;
2581 *dct
= LOWER_HALF(w
);
2590 /* {{{ s_dadd(a, b) */
2592 STATIC
void s_dadd(mp_int a
, mp_digit b
)
2595 mp_digit
*da
= MP_DIGITS(a
);
2596 mp_size ua
= MP_USED(a
);
2598 w
= (mp_word
)*da
+ b
;
2599 *da
++ = LOWER_HALF(w
);
2602 for(ua
-= 1; ua
> 0; --ua
, ++da
) {
2603 w
= (mp_word
)*da
+ w
;
2605 *da
= LOWER_HALF(w
);
2617 /* {{{ s_dmul(a, b) */
2619 STATIC
void s_dmul(mp_int a
, mp_digit b
)
2622 mp_digit
*da
= MP_DIGITS(a
);
2623 mp_size ua
= MP_USED(a
);
2626 w
= (mp_word
)*da
* b
+ w
;
2627 *da
++ = LOWER_HALF(w
);
2640 /* {{{ s_dbmul(da, b, dc, size_a) */
2642 STATIC
void s_dbmul(mp_digit
*da
, mp_digit b
, mp_digit
*dc
, mp_size size_a
)
2647 w
= (mp_word
)*da
++ * (mp_word
)b
+ w
;
2649 *dc
++ = LOWER_HALF(w
);
2655 *dc
= LOWER_HALF(w
);
2660 /* {{{ s_ddiv(da, d, dc, size_a) */
2662 STATIC mp_digit
s_ddiv(mp_int a
, mp_digit b
)
2664 mp_word w
= 0, qdigit
;
2665 mp_size ua
= MP_USED(a
);
2666 mp_digit
*da
= MP_DIGITS(a
) + ua
- 1;
2668 for(/* */; ua
> 0; --ua
, --da
) {
2669 w
= (w
<< MP_DIGIT_BIT
) | *da
;
2679 *da
= (mp_digit
)qdigit
;
2688 /* {{{ s_qdiv(z, p2) */
2690 STATIC
void s_qdiv(mp_int z
, mp_size p2
)
2692 mp_size ndig
= p2
/ MP_DIGIT_BIT
, nbits
= p2
% MP_DIGIT_BIT
;
2693 mp_size uz
= MP_USED(z
);
2697 mp_digit
*to
, *from
;
2704 to
= MP_DIGITS(z
); from
= to
+ ndig
;
2706 for(mark
= ndig
; mark
< uz
; ++mark
)
2709 MP_USED(z
) = uz
- ndig
;
2713 mp_digit d
= 0, *dz
, save
;
2714 mp_size up
= MP_DIGIT_BIT
- nbits
;
2717 dz
= MP_DIGITS(z
) + uz
- 1;
2719 for(/* */; uz
> 0; --uz
, --dz
) {
2722 *dz
= (*dz
>> nbits
) | (d
<< up
);
2729 if(MP_USED(z
) == 1 && z
->digits
[0] == 0)
2730 MP_SIGN(z
) = MP_ZPOS
;
2735 /* {{{ s_qmod(z, p2) */
2737 STATIC
void s_qmod(mp_int z
, mp_size p2
)
2739 mp_size start
= p2
/ MP_DIGIT_BIT
+ 1, rest
= p2
% MP_DIGIT_BIT
;
2740 mp_size uz
= MP_USED(z
);
2741 mp_digit mask
= (1 << rest
) - 1;
2745 z
->digits
[start
- 1] &= mask
;
2752 /* {{{ s_qmul(z, p2) */
2754 STATIC
int s_qmul(mp_int z
, mp_size p2
)
2756 mp_size uz
, need
, rest
, extra
, i
;
2757 mp_digit
*from
, *to
, d
;
2763 need
= p2
/ MP_DIGIT_BIT
; rest
= p2
% MP_DIGIT_BIT
;
2765 /* Figure out if we need an extra digit at the top end; this occurs
2766 if the topmost `rest' bits of the high-order digit of z are not
2767 zero, meaning they will be shifted off the end if not preserved */
2770 mp_digit
*dz
= MP_DIGITS(z
) + uz
- 1;
2772 if((*dz
>> (MP_DIGIT_BIT
- rest
)) != 0)
2776 if(!s_pad(z
, uz
+ need
+ extra
))
2779 /* If we need to shift by whole digits, do that in one pass, then
2780 to back and shift by partial digits.
2783 from
= MP_DIGITS(z
) + uz
- 1;
2786 for(i
= 0; i
< uz
; ++i
)
2789 ZERO(MP_DIGITS(z
), need
);
2795 for(i
= need
, from
= MP_DIGITS(z
) + need
; i
< uz
; ++i
, ++from
) {
2796 mp_digit save
= *from
;
2798 *from
= (*from
<< rest
) | (d
>> (MP_DIGIT_BIT
- rest
));
2802 d
>>= (MP_DIGIT_BIT
- rest
);
2817 /* {{{ s_qsub(z, p2) */
2819 /* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z|
2820 The sign of the result is always zero/positive.
2822 STATIC
int s_qsub(mp_int z
, mp_size p2
)
2824 mp_digit hi
= (1 << (p2
% MP_DIGIT_BIT
)), *zp
;
2825 mp_size tdig
= (p2
/ MP_DIGIT_BIT
), pos
;
2828 if(!s_pad(z
, tdig
+ 1))
2831 for(pos
= 0, zp
= MP_DIGITS(z
); pos
< tdig
; ++pos
, ++zp
) {
2832 w
= ((mp_word
) MP_DIGIT_MAX
+ 1) - w
- (mp_word
)*zp
;
2834 *zp
= LOWER_HALF(w
);
2835 w
= UPPER_HALF(w
) ? 0 : 1;
2838 w
= ((mp_word
) MP_DIGIT_MAX
+ 1 + hi
) - w
- (mp_word
)*zp
;
2839 *zp
= LOWER_HALF(w
);
2841 assert(UPPER_HALF(w
) != 0); /* no borrow out should be possible */
2843 MP_SIGN(z
) = MP_ZPOS
;
2853 STATIC
int s_dp2k(mp_int z
)
2856 mp_digit
*dp
= MP_DIGITS(z
), d
;
2858 if(MP_USED(z
) == 1 && *dp
== 0)
2867 while((d
& 1) == 0) {
2879 STATIC
int s_isp2(mp_int z
)
2881 mp_size uz
= MP_USED(z
), k
= 0;
2882 mp_digit
*dz
= MP_DIGITS(z
), d
;
2903 /* {{{ s_2expt(z, k) */
2905 STATIC
int s_2expt(mp_int z
, mp_small k
)
2910 ndig
= (k
+ MP_DIGIT_BIT
) / MP_DIGIT_BIT
;
2911 rest
= k
% MP_DIGIT_BIT
;
2918 *(dz
+ ndig
- 1) = (1 << rest
);
2926 /* {{{ s_norm(a, b) */
2928 STATIC
int s_norm(mp_int a
, mp_int b
)
2930 mp_digit d
= b
->digits
[MP_USED(b
) - 1];
2933 while(d
< (mp_digit
) (1 << (MP_DIGIT_BIT
- 1))) { /* d < (MP_RADIX / 2) */
2938 /* These multiplications can't fail */
2940 (void) s_qmul(a
, (mp_size
) k
);
2941 (void) s_qmul(b
, (mp_size
) k
);
2949 /* {{{ s_brmu(z, m) */
2951 STATIC mp_result
s_brmu(mp_int z
, mp_int m
)
2953 mp_size um
= MP_USED(m
) * 2;
2958 s_2expt(z
, MP_DIGIT_BIT
* um
);
2959 return mp_int_div(z
, m
, z
, NULL
);
2964 /* {{{ s_reduce(x, m, mu, q1, q2) */
2966 STATIC
int s_reduce(mp_int x
, mp_int m
, mp_int mu
, mp_int q1
, mp_int q2
)
2968 mp_size um
= MP_USED(m
), umb_p1
, umb_m1
;
2970 umb_p1
= (um
+ 1) * MP_DIGIT_BIT
;
2971 umb_m1
= (um
- 1) * MP_DIGIT_BIT
;
2973 if(mp_int_copy(x
, q1
) != MP_OK
)
2976 /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
2981 /* Set x = x mod b^(k+1) */
2984 /* Now, q is a guess for the quotient a / m.
2985 Compute x - q * m mod b^(k+1), replacing x. This may be off
2986 by a factor of 2m, but no more than that.
2990 (void) mp_int_sub(x
, q1
, x
); /* can't fail */
2992 /* The result may be < 0; if it is, add b^(k+1) to pin it in the
2994 if((CMPZ(x
) < 0) && !s_qsub(x
, umb_p1
))
2997 /* If x > m, we need to back it off until it is in range.
2998 This will be required at most twice. */
2999 if(mp_int_compare(x
, m
) >= 0) {
3000 (void) mp_int_sub(x
, m
, x
);
3001 if(mp_int_compare(x
, m
) >= 0)
3002 (void) mp_int_sub(x
, m
, x
);
3005 /* At this point, x has been properly reduced. */
3011 /* {{{ s_embar(a, b, m, mu, c) */
3013 /* Perform modular exponentiation using Barrett's method, where mu is
3014 the reduction constant for m. Assumes a < m, b > 0. */
3015 STATIC mp_result
s_embar(mp_int a
, mp_int b
, mp_int m
, mp_int mu
, mp_int c
)
3017 mp_digit
*db
, *dbt
, umu
, d
;
3022 umu
= MP_USED(mu
); db
= MP_DIGITS(b
); dbt
= db
+ MP_USED(b
) - 1;
3025 SETUP(mp_int_init_size(TEMP(last
), 4 * umu
), last
);
3026 ZERO(MP_DIGITS(TEMP(last
- 1)), MP_ALLOC(TEMP(last
- 1)));
3029 (void) mp_int_set_value(c
, 1);
3031 /* Take care of low-order digits */
3035 for(d
= *db
, i
= MP_DIGIT_BIT
; i
> 0; --i
, d
>>= 1) {
3037 /* The use of a second temporary avoids allocation */
3038 UMUL(c
, a
, TEMP(0));
3039 if(!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2))) {
3040 res
= MP_MEMORY
; goto CLEANUP
;
3042 mp_int_copy(TEMP(0), c
);
3047 assert(MP_SIGN(TEMP(0)) == MP_ZPOS
);
3048 if(!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2))) {
3049 res
= MP_MEMORY
; goto CLEANUP
;
3051 assert(MP_SIGN(TEMP(0)) == MP_ZPOS
);
3052 mp_int_copy(TEMP(0), a
);
3060 /* Take care of highest-order digit */
3064 UMUL(c
, a
, TEMP(0));
3065 if(!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2))) {
3066 res
= MP_MEMORY
; goto CLEANUP
;
3068 mp_int_copy(TEMP(0), c
);
3075 if(!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2))) {
3076 res
= MP_MEMORY
; goto CLEANUP
;
3078 (void) mp_int_copy(TEMP(0), a
);
3083 mp_int_clear(TEMP(last
));
3090 /* {{{ s_udiv(a, b) */
3092 /* Precondition: a >= b and b > 0
3093 Postcondition: a' = a / b, b' = a % b
3095 STATIC mp_result
s_udiv(mp_int a
, mp_int b
)
3098 mp_size ua
, ub
, qpos
= 0;
3100 mp_result res
= MP_OK
;
3103 /* Force signs to positive */
3104 MP_SIGN(a
) = MP_ZPOS
;
3105 MP_SIGN(b
) = MP_ZPOS
;
3107 /* Normalize, per Knuth */
3110 ua
= MP_USED(a
); ub
= MP_USED(b
); btop
= b
->digits
[ub
- 1];
3111 if((res
= mp_int_init_size(&q
, ua
)) != MP_OK
) return res
;
3112 if((res
= mp_int_init_size(&t
, ua
+ 1)) != MP_OK
) goto CLEANUP
;
3115 r
.digits
= da
+ ua
- 1; /* The contents of r are shared with a */
3118 r
.alloc
= MP_ALLOC(a
);
3119 ZERO(t
.digits
, t
.alloc
);
3121 /* Solve for quotient digits, store in q.digits in reverse order */
3122 while(r
.digits
>= da
) {
3123 assert(qpos
<= q
.alloc
);
3125 if(s_ucmp(b
, &r
) > 0) {
3129 if(++skip
> 1 && qpos
> 0)
3130 q
.digits
[qpos
++] = 0;
3135 mp_word pfx
= r
.digits
[r
.used
- 1];
3138 if(r
.used
> 1 && pfx
<= btop
) {
3139 pfx
<<= MP_DIGIT_BIT
/ 2;
3140 pfx
<<= MP_DIGIT_BIT
/ 2;
3141 pfx
|= r
.digits
[r
.used
- 2];
3144 qdigit
= pfx
/ btop
;
3145 if(qdigit
> MP_DIGIT_MAX
) {
3146 qdigit
= MP_DIGIT_MAX
;
3149 s_dbmul(MP_DIGITS(b
), (mp_digit
) qdigit
, t
.digits
, ub
);
3150 t
.used
= ub
+ 1; CLAMP(&t
);
3151 while(s_ucmp(&t
, &r
) > 0) {
3153 (void) mp_int_sub(&t
, b
, &t
); /* cannot fail */
3156 s_usub(r
.digits
, t
.digits
, r
.digits
, r
.used
, t
.used
);
3159 q
.digits
[qpos
++] = (mp_digit
) qdigit
;
3160 ZERO(t
.digits
, t
.used
);
3165 /* Put quotient digits in the correct order, and discard extra zeroes */
3167 REV(mp_digit
, q
.digits
, qpos
);
3170 /* Denormalize the remainder */
3175 mp_int_copy(a
, b
); /* ok: 0 <= r < b */
3176 mp_int_copy(&q
, a
); /* ok: q <= a */
3186 /* {{{ s_outlen(z, r) */
3188 STATIC
int s_outlen(mp_int z
, mp_size r
)
3193 assert(r
>= MP_MIN_RADIX
&& r
<= MP_MAX_RADIX
);
3195 bits
= mp_int_count_bits(z
);
3196 raw
= (double)bits
* s_log2
[r
];
3198 return (int)(raw
+ 0.999999);
3203 /* {{{ s_inlen(len, r) */
3205 STATIC mp_size
s_inlen(int len
, mp_size r
)
3207 double raw
= (double)len
/ s_log2
[r
];
3208 mp_size bits
= (mp_size
)(raw
+ 0.5);
3210 return (mp_size
)((bits
+ (MP_DIGIT_BIT
- 1)) / MP_DIGIT_BIT
);
3215 /* {{{ s_ch2val(c, r) */
3217 STATIC
int s_ch2val(char c
, int r
)
3221 if(isdigit((unsigned char) c
))
3223 else if(r
> 10 && isalpha((unsigned char) c
))
3224 out
= toupper(c
) - 'A' + 10;
3228 return (out
>= r
) ? -1 : out
;
3233 /* {{{ s_val2ch(v, caps) */
3235 STATIC
char s_val2ch(int v
, int caps
)
3242 char out
= (v
- 10) + 'a';
3245 return toupper(out
);
3253 /* {{{ s_2comp(buf, len) */
3255 STATIC
void s_2comp(unsigned char *buf
, int len
)
3258 unsigned short s
= 1;
3260 for(i
= len
- 1; i
>= 0; --i
) {
3261 unsigned char c
= ~buf
[i
];
3270 /* last carry out is ignored */
3275 /* {{{ s_tobin(z, buf, *limpos) */
3277 STATIC mp_result
s_tobin(mp_int z
, unsigned char *buf
, int *limpos
, int pad
)
3281 int pos
= 0, limit
= *limpos
;
3283 uz
= MP_USED(z
); dz
= MP_DIGITS(z
);
3284 while(uz
> 0 && pos
< limit
) {
3288 for(i
= sizeof(mp_digit
); i
> 0 && pos
< limit
; --i
) {
3289 buf
[pos
++] = (unsigned char)d
;
3292 /* Don't write leading zeroes */
3293 if(d
== 0 && uz
== 1)
3294 i
= 0; /* exit loop without signaling truncation */
3297 /* Detect truncation (loop exited with pos >= limit) */
3303 if(pad
!= 0 && (buf
[pos
- 1] >> (CHAR_BIT
- 1))) {
3310 /* Digits are in reverse order, fix that */
3311 REV(unsigned char, buf
, pos
);
3313 /* Return the number of bytes actually written */
3316 return (uz
== 0) ? MP_OK
: MP_TRUNC
;
3321 /* {{{ s_print(tag, z) */
3324 void s_print(char *tag
, mp_int z
)
3328 fprintf(stderr
, "%s: %c ", tag
,
3329 (MP_SIGN(z
) == MP_NEG
) ? '-' : '+');
3331 for(i
= MP_USED(z
) - 1; i
>= 0; --i
)
3332 fprintf(stderr
, "%0*X", (int)(MP_DIGIT_BIT
/ 4), z
->digits
[i
]);
3334 fputc('\n', stderr
);
3338 void s_print_buf(char *tag
, mp_digit
*buf
, mp_size num
)
3342 fprintf(stderr
, "%s: ", tag
);
3344 for(i
= num
- 1; i
>= 0; --i
)
3345 fprintf(stderr
, "%0*X", (int)(MP_DIGIT_BIT
/ 4), buf
[i
]);
3347 fputc('\n', stderr
);
3353 /* HERE THERE BE DRAGONS */