1 /* imath version 1.3 */
4 Purpose: Arbitrary precision integer arithmetic routines.
5 Author: M. J. Fromberger <http://www.dartmouth.edu/~sting/>
6 Info: Id: imath.c 21 2006-04-02 18:58:36Z sting
8 Copyright (C) 2002 Michael J. Fromberger, All Rights Reserved.
10 Permission is hereby granted, free of charge, to any person
11 obtaining a copy of this software and associated documentation files
12 (the "Software"), to deal in the Software without restriction,
13 including without limitation the rights to use, copy, modify, merge,
14 publish, distribute, sublicense, and/or sell copies of the Software,
15 and to permit persons to whom the Software is furnished to do so,
16 subject to the following conditions:
18 The above copyright notice and this permission notice shall be
19 included in all copies or substantial portions of the Software.
21 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
22 EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
23 MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
24 NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS
25 BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
26 ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
27 CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
37 #define assert(TEST) Assert(TEST)
38 #define TRACEABLE_CLAMP 0
39 #define TRACEABLE_FREE 0
43 const mp_result MP_OK
= 0; /* no error, all is well */
44 const mp_result MP_FALSE
= 0; /* boolean false */
45 const mp_result MP_TRUE
= -1; /* boolean true */
46 const mp_result MP_MEMORY
= -2; /* out of memory */
47 const mp_result MP_RANGE
= -3; /* argument out of range */
48 const mp_result MP_UNDEF
= -4; /* result undefined */
49 const mp_result MP_TRUNC
= -5; /* output truncated */
50 const mp_result MP_BADARG
= -6; /* invalid null argument */
52 const mp_sign MP_NEG
= 1; /* value is strictly negative */
53 const mp_sign MP_ZPOS
= 0; /* value is non-negative */
55 static const char *s_unknown_err
= "unknown result code";
56 static const char *s_error_msg
[] = {
60 "argument out of range",
63 "invalid null argument",
69 /* Optional library flags */
70 #define MP_CAP_DIGITS 1 /* flag bit to capitalize letter digits */
72 /* Argument checking macros
73 Use CHECK() where a return value is required; NRCHECK() elsewhere */
74 #define CHECK(TEST) assert(TEST)
75 #define NRCHECK(TEST) assert(TEST)
77 /* {{{ Logarithm table for computing output sizes */
79 /* The ith entry of this table gives the value of log_i(2).
81 An integer value n requires ceil(log_i(n)) digits to be represented
82 in base i. Since it is easy to compute lg(n), by counting bits, we
83 can compute log_i(n) = lg(n) * log_i(2).
85 static const double s_log2
[] = {
86 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* 0 1 2 3 */
87 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */
88 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */
89 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */
90 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */
91 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */
92 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */
93 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */
94 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */
95 0.193426404, 0.191958720, 0.190551412, 0.189200360, /* 36 37 38 39 */
96 0.187901825, 0.186652411, 0.185449023, 0.184288833, /* 40 41 42 43 */
97 0.183169251, 0.182087900, 0.181042597, 0.180031327, /* 44 45 46 47 */
98 0.179052232, 0.178103594, 0.177183820, 0.176291434, /* 48 49 50 51 */
99 0.175425064, 0.174583430, 0.173765343, 0.172969690, /* 52 53 54 55 */
100 0.172195434, 0.171441601, 0.170707280, 0.169991616, /* 56 57 58 59 */
101 0.169293808, 0.168613099, 0.167948779, 0.167300179, /* 60 61 62 63 */
106 /* {{{ Various macros */
108 /* Return the number of digits needed to represent a static value */
109 #define MP_VALUE_DIGITS(V) \
110 ((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit))
112 /* Round precision P to nearest word boundary */
113 #define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2)))
115 /* Set array P of S digits to zero */
117 do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P);memset(p__,0,i__);}while(0)
119 /* Copy S digits from array P to array Q */
120 #define COPY(P, Q, S) \
121 do{mp_size i__=(S)*sizeof(mp_digit);mp_digit *p__=(P),*q__=(Q);\
122 memcpy(q__,p__,i__);}while(0)
124 /* Reverse N elements of type T in array A */
125 #define REV(T, A, N) \
126 do{T *u_=(A),*v_=u_+(N)-1;while(u_<v_){T xch=*u_;*u_++=*v_;*v_--=xch;}}while(0)
129 #define CLAMP(Z) s_clamp(Z)
132 do{mp_int z_=(Z);mp_size uz_=MP_USED(z_);mp_digit *dz_=MP_DIGITS(z_)+uz_-1;\
133 while(uz_ > 1 && (*dz_-- == 0)) --uz_;MP_USED(z_)=uz_;}while(0)
138 #define MIN(A, B) ((B)<(A)?(B):(A))
139 #define MAX(A, B) ((B)>(A)?(B):(A))
140 #define SWAP(T, A, B) do{T t_=(A);A=(B);B=t_;}while(0)
142 #define TEMP(K) (temp + (K))
143 #define SETUP(E, C) \
144 do{if((res = (E)) != MP_OK) goto CLEANUP; ++(C);}while(0)
147 (((Z)->used==1&&(Z)->digits[0]==0)?0:((Z)->sign==MP_NEG)?-1:1)
149 #define UMUL(X, Y, Z) \
150 do{mp_size ua_=MP_USED(X),ub_=MP_USED(Y);mp_size o_=ua_+ub_;\
151 ZERO(MP_DIGITS(Z),o_);\
152 (void) s_kmul(MP_DIGITS(X),MP_DIGITS(Y),MP_DIGITS(Z),ua_,ub_);\
153 MP_USED(Z)=o_;CLAMP(Z);}while(0)
156 do{mp_size ua_=MP_USED(X),o_=ua_+ua_;ZERO(MP_DIGITS(Z),o_);\
157 (void) s_ksqr(MP_DIGITS(X),MP_DIGITS(Z),ua_);MP_USED(Z)=o_;CLAMP(Z);}while(0)
159 #define UPPER_HALF(W) ((mp_word)((W) >> MP_DIGIT_BIT))
160 #define LOWER_HALF(W) ((mp_digit)(W))
161 #define HIGH_BIT_SET(W) ((W) >> (MP_WORD_BIT - 1))
162 #define ADD_WILL_OVERFLOW(W, V) ((MP_WORD_MAX - (V)) < (W))
166 /* Default number of digits allocated to a new mp_int */
167 static mp_size default_precision
= 64;
169 /* Minimum number of digits to invoke recursive multiply */
170 static mp_size multiply_threshold
= 32;
172 /* Default library configuration flags */
173 static mp_word mp_flags
= MP_CAP_DIGITS
;
175 /* Allocate a buffer of (at least) num digits, or return
176 NULL if that couldn't be done. */
177 static mp_digit
*s_alloc(mp_size num
);
180 static void s_free(void *ptr
);
182 #define s_free(P) px_free(P)
185 /* Insure that z has at least min digits allocated, resizing if
186 necessary. Returns true if successful, false if out of memory. */
187 static int s_pad(mp_int z
, mp_size min
);
189 /* Normalize by removing leading zeroes (except when z = 0) */
191 static void s_clamp(mp_int z
);
194 /* Fill in a "fake" mp_int on the stack with a given value */
195 static void s_fake(mp_int z
, int value
, mp_digit vbuf
[]);
197 /* Compare two runs of digits of given length, returns <0, 0, >0 */
198 static int s_cdig(mp_digit
* da
, mp_digit
* db
, mp_size len
);
200 /* Pack the unsigned digits of v into array t */
201 static int s_vpack(int v
, mp_digit t
[]);
203 /* Compare magnitudes of a and b, returns <0, 0, >0 */
204 static int s_ucmp(mp_int a
, mp_int b
);
206 /* Compare magnitudes of a and v, returns <0, 0, >0 */
207 static int s_vcmp(mp_int a
, int v
);
209 /* Unsigned magnitude addition; assumes dc is big enough.
210 Carry out is returned (no memory allocated). */
211 static mp_digit
s_uadd(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
212 mp_size size_a
, mp_size size_b
);
214 /* Unsigned magnitude subtraction. Assumes dc is big enough. */
215 static void s_usub(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
216 mp_size size_a
, mp_size size_b
);
218 /* Unsigned recursive multiplication. Assumes dc is big enough. */
219 static int s_kmul(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
220 mp_size size_a
, mp_size size_b
);
222 /* Unsigned magnitude multiplication. Assumes dc is big enough. */
223 static void s_umul(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
224 mp_size size_a
, mp_size size_b
);
226 /* Unsigned recursive squaring. Assumes dc is big enough. */
227 static int s_ksqr(mp_digit
* da
, mp_digit
* dc
, mp_size size_a
);
229 /* Unsigned magnitude squaring. Assumes dc is big enough. */
230 static void s_usqr(mp_digit
* da
, mp_digit
* dc
, mp_size size_a
);
232 /* Single digit addition. Assumes a is big enough. */
233 static void s_dadd(mp_int a
, mp_digit b
);
235 /* Single digit multiplication. Assumes a is big enough. */
236 static void s_dmul(mp_int a
, mp_digit b
);
238 /* Single digit multiplication on buffers; assumes dc is big enough. */
239 static void s_dbmul(mp_digit
* da
, mp_digit b
, mp_digit
* dc
,
242 /* Single digit division. Replaces a with the quotient,
243 returns the remainder. */
244 static mp_digit
s_ddiv(mp_int a
, mp_digit b
);
246 /* Quick division by a power of 2, replaces z (no allocation) */
247 static void s_qdiv(mp_int z
, mp_size p2
);
249 /* Quick remainder by a power of 2, replaces z (no allocation) */
250 static void s_qmod(mp_int z
, mp_size p2
);
252 /* Quick multiplication by a power of 2, replaces z.
253 Allocates if necessary; returns false in case this fails. */
254 static int s_qmul(mp_int z
, mp_size p2
);
256 /* Quick subtraction from a power of 2, replaces z.
257 Allocates if necessary; returns false in case this fails. */
258 static int s_qsub(mp_int z
, mp_size p2
);
260 /* Return maximum k such that 2^k divides z. */
261 static int s_dp2k(mp_int z
);
263 /* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */
264 static int s_isp2(mp_int z
);
266 /* Set z to 2^k. May allocate; returns false in case this fails. */
267 static int s_2expt(mp_int z
, int k
);
269 /* Normalize a and b for division, returns normalization constant */
270 static int s_norm(mp_int a
, mp_int b
);
272 /* Compute constant mu for Barrett reduction, given modulus m, result
273 replaces z, m is untouched. */
274 static mp_result
s_brmu(mp_int z
, mp_int m
);
276 /* Reduce a modulo m, using Barrett's algorithm. */
277 static int s_reduce(mp_int x
, mp_int m
, mp_int mu
, mp_int q1
, mp_int q2
);
279 /* Modular exponentiation, using Barrett reduction */
280 static mp_result
s_embar(mp_int a
, mp_int b
, mp_int m
, mp_int mu
, mp_int c
);
282 /* Unsigned magnitude division. Assumes |a| > |b|. Allocates
283 temporaries; overwrites a with quotient, b with remainder. */
284 static mp_result
s_udiv(mp_int a
, mp_int b
);
286 /* Compute the number of digits in radix r required to represent the
287 given value. Does not account for sign flags, terminators, etc. */
288 static int s_outlen(mp_int z
, mp_size r
);
290 /* Guess how many digits of precision will be needed to represent a
291 radix r value of the specified number of digits. Returns a value
292 guaranteed to be no smaller than the actual number required. */
293 static mp_size
s_inlen(int len
, mp_size r
);
295 /* Convert a character to a digit value in radix r, or
296 -1 if out of range */
297 static int s_ch2val(char c
, int r
);
299 /* Convert a digit value to a character */
300 static char s_val2ch(int v
, int caps
);
302 /* Take 2's complement of a buffer in place */
303 static void s_2comp(unsigned char *buf
, int len
);
305 /* Convert a value to binary, ignoring sign. On input, *limpos is the
306 bound on how many bytes should be written to buf; on output, *limpos
307 is set to the number of bytes actually written. */
308 static mp_result
s_tobin(mp_int z
, unsigned char *buf
, int *limpos
, int pad
);
311 /* Dump a representation of the mp_int to standard output */
312 void s_print(char *tag
, mp_int z
);
313 void s_print_buf(char *tag
, mp_digit
* buf
, mp_size num
);
316 /* {{{ get_default_precision() */
319 mp_get_default_precision(void)
321 return default_precision
;
326 /* {{{ mp_set_default_precision(s) */
329 mp_set_default_precision(mp_size s
)
333 default_precision
= (mp_size
) ROUND_PREC(s
);
338 /* {{{ mp_get_multiply_threshold() */
341 mp_get_multiply_threshold(void)
343 return multiply_threshold
;
348 /* {{{ mp_set_multiply_threshold(s) */
351 mp_set_multiply_threshold(mp_size s
)
353 multiply_threshold
= s
;
358 /* {{{ mp_int_init(z) */
361 mp_int_init(mp_int z
)
363 return mp_int_init_size(z
, default_precision
);
368 /* {{{ mp_int_alloc() */
373 mp_int out
= px_alloc(sizeof(mpz_t
));
386 /* {{{ mp_int_init_size(z, prec) */
389 mp_int_init_size(mp_int z
, mp_size prec
)
393 prec
= (mp_size
) ROUND_PREC(prec
);
394 prec
= MAX(prec
, default_precision
);
396 if ((MP_DIGITS(z
) = s_alloc(prec
)) == NULL
)
402 MP_SIGN(z
) = MP_ZPOS
;
409 /* {{{ mp_int_init_copy(z, old) */
412 mp_int_init_copy(mp_int z
, mp_int old
)
418 CHECK(z
!= NULL
&& old
!= NULL
);
421 target
= MAX(uold
, default_precision
);
423 if ((res
= mp_int_init_size(z
, target
)) != MP_OK
)
427 MP_SIGN(z
) = MP_SIGN(old
);
428 COPY(MP_DIGITS(old
), MP_DIGITS(z
), uold
);
435 /* {{{ mp_int_init_value(z, value) */
438 mp_int_init_value(mp_int z
, int value
)
444 if ((res
= mp_int_init(z
)) != MP_OK
)
447 return mp_int_set_value(z
, value
);
452 /* {{{ mp_int_set_value(z, value) */
455 mp_int_set_value(mp_int z
, int value
)
461 /* How many digits to copy */
462 ndig
= (mp_size
) MP_VALUE_DIGITS(value
);
467 MP_USED(z
) = (mp_size
) s_vpack(value
, MP_DIGITS(z
));
468 MP_SIGN(z
) = (value
< 0) ? MP_NEG
: MP_ZPOS
;
475 /* {{{ mp_int_clear(z) */
478 mp_int_clear(mp_int z
)
483 if (MP_DIGITS(z
) != NULL
)
485 s_free(MP_DIGITS(z
));
492 /* {{{ mp_int_free(z) */
495 mp_int_free(mp_int z
)
499 if (z
->digits
!= NULL
)
507 /* {{{ mp_int_copy(a, c) */
510 mp_int_copy(mp_int a
, mp_int c
)
512 CHECK(a
!= NULL
&& c
!= NULL
);
516 mp_size ua
= MP_USED(a
);
528 MP_SIGN(c
) = MP_SIGN(a
);
536 /* {{{ mp_int_swap(a, c) */
539 mp_int_swap(mp_int a
, mp_int c
)
552 /* {{{ mp_int_zero(z) */
555 mp_int_zero(mp_int z
)
561 MP_SIGN(z
) = MP_ZPOS
;
566 /* {{{ mp_int_abs(a, c) */
569 mp_int_abs(mp_int a
, mp_int c
)
573 CHECK(a
!= NULL
&& c
!= NULL
);
575 if ((res
= mp_int_copy(a
, c
)) != MP_OK
)
578 MP_SIGN(c
) = MP_ZPOS
;
584 /* {{{ mp_int_neg(a, c) */
587 mp_int_neg(mp_int a
, mp_int c
)
591 CHECK(a
!= NULL
&& c
!= NULL
);
593 if ((res
= mp_int_copy(a
, c
)) != MP_OK
)
597 MP_SIGN(c
) = 1 - MP_SIGN(a
);
604 /* {{{ mp_int_add(a, b, c) */
607 mp_int_add(mp_int a
, mp_int b
, mp_int c
)
614 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
621 if (MP_SIGN(a
) == MP_SIGN(b
))
623 /* Same sign -- add magnitudes, preserve sign of addends */
629 carry
= s_uadd(MP_DIGITS(a
), MP_DIGITS(b
), MP_DIGITS(c
), ua
, ub
);
634 if (!s_pad(c
, max
+ 1))
637 c
->digits
[max
] = carry
;
642 MP_SIGN(c
) = MP_SIGN(a
);
647 /* Different signs -- subtract magnitudes, preserve sign of greater */
650 int cmp
= s_ucmp(a
, b
); /* magnitude comparision, sign ignored */
652 /* Set x to max(a, b), y to min(a, b) to simplify later code */
664 if (!s_pad(c
, MP_USED(x
)))
667 /* Subtract smaller from larger */
668 s_usub(MP_DIGITS(x
), MP_DIGITS(y
), MP_DIGITS(c
), MP_USED(x
), MP_USED(y
));
669 MP_USED(c
) = MP_USED(x
);
672 /* Give result the sign of the larger */
673 MP_SIGN(c
) = MP_SIGN(x
);
681 /* {{{ mp_int_add_value(a, value, c) */
684 mp_int_add_value(mp_int a
, int value
, mp_int c
)
687 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
689 s_fake(&vtmp
, value
, vbuf
);
691 return mp_int_add(a
, &vtmp
, c
);
696 /* {{{ mp_int_sub(a, b, c) */
699 mp_int_sub(mp_int a
, mp_int b
, mp_int c
)
706 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
713 if (MP_SIGN(a
) != MP_SIGN(b
))
715 /* Different signs -- add magnitudes and keep sign of a */
721 carry
= s_uadd(MP_DIGITS(a
), MP_DIGITS(b
), MP_DIGITS(c
), ua
, ub
);
726 if (!s_pad(c
, max
+ 1))
729 c
->digits
[max
] = carry
;
734 MP_SIGN(c
) = MP_SIGN(a
);
739 /* Same signs -- subtract magnitudes */
743 int cmp
= s_ucmp(a
, b
);
761 if (MP_SIGN(a
) == MP_NEG
&& cmp
!= 0)
764 s_usub(MP_DIGITS(x
), MP_DIGITS(y
), MP_DIGITS(c
), MP_USED(x
), MP_USED(y
));
765 MP_USED(c
) = MP_USED(x
);
776 /* {{{ mp_int_sub_value(a, value, c) */
779 mp_int_sub_value(mp_int a
, int value
, mp_int c
)
782 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
784 s_fake(&vtmp
, value
, vbuf
);
786 return mp_int_sub(a
, &vtmp
, c
);
791 /* {{{ mp_int_mul(a, b, c) */
794 mp_int_mul(mp_int a
, mp_int b
, mp_int c
)
803 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
805 /* If either input is zero, we can shortcut multiplication */
806 if (mp_int_compare_zero(a
) == 0 || mp_int_compare_zero(b
) == 0)
812 /* Output is positive if inputs have same sign, otherwise negative */
813 osign
= (MP_SIGN(a
) == MP_SIGN(b
)) ? MP_ZPOS
: MP_NEG
;
816 * If the output is not equal to any of the inputs, we'll write the
817 * results there directly; otherwise, allocate a temporary space.
823 if (c
== a
|| c
== b
)
825 p
= ROUND_PREC(osize
);
826 p
= MAX(p
, default_precision
);
828 if ((out
= s_alloc(p
)) == NULL
)
833 if (!s_pad(c
, osize
))
840 if (!s_kmul(MP_DIGITS(a
), MP_DIGITS(b
), out
, ua
, ub
))
844 * If we allocated a new buffer, get rid of whatever memory c was already
845 * using, and fix up its fields to reflect that.
847 if (out
!= MP_DIGITS(c
))
849 s_free(MP_DIGITS(c
));
854 MP_USED(c
) = osize
; /* might not be true, but we'll fix it ... */
855 CLAMP(c
); /* ... right here */
863 /* {{{ mp_int_mul_value(a, value, c) */
866 mp_int_mul_value(mp_int a
, int value
, mp_int c
)
869 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
871 s_fake(&vtmp
, value
, vbuf
);
873 return mp_int_mul(a
, &vtmp
, c
);
878 /* {{{ mp_int_mul_pow2(a, p2, c) */
881 mp_int_mul_pow2(mp_int a
, int p2
, mp_int c
)
885 CHECK(a
!= NULL
&& c
!= NULL
&& p2
>= 0);
887 if ((res
= mp_int_copy(a
, c
)) != MP_OK
)
890 if (s_qmul(c
, (mp_size
) p2
))
898 /* {{{ mp_int_sqr(a, c) */
901 mp_int_sqr(mp_int a
, mp_int c
)
907 CHECK(a
!= NULL
&& c
!= NULL
);
909 /* Get a temporary buffer big enough to hold the result */
910 osize
= (mp_size
) 2 *MP_USED(a
);
914 p
= ROUND_PREC(osize
);
915 p
= MAX(p
, default_precision
);
917 if ((out
= s_alloc(p
)) == NULL
)
922 if (!s_pad(c
, osize
))
929 s_ksqr(MP_DIGITS(a
), out
, MP_USED(a
));
932 * Get rid of whatever memory c was already using, and fix up its fields
933 * to reflect the new digit array it's using
935 if (out
!= MP_DIGITS(c
))
937 s_free(MP_DIGITS(c
));
942 MP_USED(c
) = osize
; /* might not be true, but we'll fix it ... */
943 CLAMP(c
); /* ... right here */
944 MP_SIGN(c
) = MP_ZPOS
;
951 /* {{{ mp_int_div(a, b, q, r) */
954 mp_int_div(mp_int a
, mp_int b
, mp_int q
, mp_int r
)
959 mp_result res
= MP_OK
;
963 mp_sign sa
= MP_SIGN(a
),
966 CHECK(a
!= NULL
&& b
!= NULL
&& q
!= r
);
970 else if ((cmp
= s_ucmp(a
, b
)) < 0)
973 * If |a| < |b|, no division is required: q = 0, r = a
975 if (r
&& (res
= mp_int_copy(a
, r
)) != MP_OK
)
986 * If |a| = |b|, no division is required: q = 1 or -1, r = 0
1004 * When |a| > |b|, real division is required. We need someplace to store
1005 * quotient and remainder, but q and r are allowed to be NULL or to
1006 * overlap with the inputs.
1008 if ((lg
= s_isp2(b
)) < 0)
1010 if (q
&& b
!= q
&& (res
= mp_int_copy(a
, q
)) == MP_OK
)
1017 SETUP(mp_int_init_copy(TEMP(last
), a
), last
);
1020 if (r
&& a
!= r
&& (res
= mp_int_copy(b
, r
)) == MP_OK
)
1027 SETUP(mp_int_init_copy(TEMP(last
), b
), last
);
1030 if ((res
= s_udiv(qout
, rout
)) != MP_OK
)
1035 if (q
&& (res
= mp_int_copy(a
, q
)) != MP_OK
)
1037 if (r
&& (res
= mp_int_copy(a
, r
)) != MP_OK
)
1041 s_qdiv(q
, (mp_size
) lg
);
1044 s_qmod(r
, (mp_size
) lg
);
1048 /* Recompute signs for output */
1052 if (CMPZ(rout
) == 0)
1053 MP_SIGN(rout
) = MP_ZPOS
;
1057 MP_SIGN(qout
) = (sa
== sb
) ? MP_ZPOS
: MP_NEG
;
1058 if (CMPZ(qout
) == 0)
1059 MP_SIGN(qout
) = MP_ZPOS
;
1062 if (q
&& (res
= mp_int_copy(qout
, q
)) != MP_OK
)
1064 if (r
&& (res
= mp_int_copy(rout
, r
)) != MP_OK
)
1069 mp_int_clear(TEMP(last
));
1076 /* {{{ mp_int_mod(a, m, c) */
1079 mp_int_mod(mp_int a
, mp_int m
, mp_int c
)
1087 if ((res
= mp_int_init(&tmp
)) != MP_OK
)
1097 if ((res
= mp_int_div(a
, m
, NULL
, out
)) != MP_OK
)
1101 res
= mp_int_add(out
, m
, c
);
1103 res
= mp_int_copy(out
, c
);
1115 /* {{{ mp_int_div_value(a, value, q, r) */
1118 mp_int_div_value(mp_int a
, int value
, mp_int q
, int *r
)
1122 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
1125 if ((res
= mp_int_init(&rtmp
)) != MP_OK
)
1127 s_fake(&vtmp
, value
, vbuf
);
1129 if ((res
= mp_int_div(a
, &vtmp
, q
, &rtmp
)) != MP_OK
)
1133 (void) mp_int_to_int(&rtmp
, r
); /* can't fail */
1136 mp_int_clear(&rtmp
);
1142 /* {{{ mp_int_div_pow2(a, p2, q, r) */
1145 mp_int_div_pow2(mp_int a
, int p2
, mp_int q
, mp_int r
)
1147 mp_result res
= MP_OK
;
1149 CHECK(a
!= NULL
&& p2
>= 0 && q
!= r
);
1151 if (q
!= NULL
&& (res
= mp_int_copy(a
, q
)) == MP_OK
)
1152 s_qdiv(q
, (mp_size
) p2
);
1154 if (res
== MP_OK
&& r
!= NULL
&& (res
= mp_int_copy(a
, r
)) == MP_OK
)
1155 s_qmod(r
, (mp_size
) p2
);
1162 /* {{{ mp_int_expt(a, b, c) */
1165 mp_int_expt(mp_int a
, int b
, mp_int c
)
1169 unsigned int v
= abs(b
);
1171 CHECK(b
>= 0 && c
!= NULL
);
1173 if ((res
= mp_int_init_copy(&t
, a
)) != MP_OK
)
1176 (void) mp_int_set_value(c
, 1);
1181 if ((res
= mp_int_mul(c
, &t
, c
)) != MP_OK
)
1189 if ((res
= mp_int_sqr(&t
, &t
)) != MP_OK
)
1200 /* {{{ mp_int_expt_value(a, b, c) */
1203 mp_int_expt_value(int a
, int b
, mp_int c
)
1207 unsigned int v
= abs(b
);
1209 CHECK(b
>= 0 && c
!= NULL
);
1211 if ((res
= mp_int_init_value(&t
, a
)) != MP_OK
)
1214 (void) mp_int_set_value(c
, 1);
1219 if ((res
= mp_int_mul(c
, &t
, c
)) != MP_OK
)
1227 if ((res
= mp_int_sqr(&t
, &t
)) != MP_OK
)
1238 /* {{{ mp_int_compare(a, b) */
1241 mp_int_compare(mp_int a
, mp_int b
)
1245 CHECK(a
!= NULL
&& b
!= NULL
);
1248 if (sa
== MP_SIGN(b
))
1250 int cmp
= s_ucmp(a
, b
);
1253 * If they're both zero or positive, the normal comparison applies; if
1254 * both negative, the sense is reversed.
1273 /* {{{ mp_int_compare_unsigned(a, b) */
1276 mp_int_compare_unsigned(mp_int a
, mp_int b
)
1278 NRCHECK(a
!= NULL
&& b
!= NULL
);
1280 return s_ucmp(a
, b
);
1285 /* {{{ mp_int_compare_zero(z) */
1288 mp_int_compare_zero(mp_int z
)
1292 if (MP_USED(z
) == 1 && z
->digits
[0] == 0)
1294 else if (MP_SIGN(z
) == MP_ZPOS
)
1302 /* {{{ mp_int_compare_value(z, value) */
1305 mp_int_compare_value(mp_int z
, int value
)
1307 mp_sign vsign
= (value
< 0) ? MP_NEG
: MP_ZPOS
;
1312 if (vsign
== MP_SIGN(z
))
1314 cmp
= s_vcmp(z
, value
);
1316 if (vsign
== MP_ZPOS
)
1332 /* {{{ mp_int_exptmod(a, b, m, c) */
1335 mp_int_exptmod(mp_int a
, mp_int b
, mp_int m
, mp_int c
)
1343 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
&& m
!= NULL
);
1345 /* Zero moduli and negative exponents are not considered. */
1352 SETUP(mp_int_init_size(TEMP(0), 2 * um
), last
);
1353 SETUP(mp_int_init_size(TEMP(1), 2 * um
), last
);
1355 if (c
== b
|| c
== m
)
1357 SETUP(mp_int_init_size(TEMP(2), 2 * um
), last
);
1365 if ((res
= mp_int_mod(a
, m
, TEMP(0))) != MP_OK
)
1368 if ((res
= s_brmu(TEMP(1), m
)) != MP_OK
)
1371 if ((res
= s_embar(TEMP(0), b
, m
, TEMP(1), s
)) != MP_OK
)
1374 res
= mp_int_copy(s
, c
);
1378 mp_int_clear(TEMP(last
));
1385 /* {{{ mp_int_exptmod_evalue(a, value, m, c) */
1388 mp_int_exptmod_evalue(mp_int a
, int value
, mp_int m
, mp_int c
)
1391 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
1393 s_fake(&vtmp
, value
, vbuf
);
1395 return mp_int_exptmod(a
, &vtmp
, m
, c
);
1400 /* {{{ mp_int_exptmod_bvalue(v, b, m, c) */
1403 mp_int_exptmod_bvalue(int value
, mp_int b
,
1407 mp_digit vbuf
[MP_VALUE_DIGITS(value
)];
1409 s_fake(&vtmp
, value
, vbuf
);
1411 return mp_int_exptmod(&vtmp
, b
, m
, c
);
1416 /* {{{ mp_int_exptmod_known(a, b, m, mu, c) */
1419 mp_int_exptmod_known(mp_int a
, mp_int b
, mp_int m
, mp_int mu
, mp_int c
)
1427 CHECK(a
&& b
&& m
&& c
);
1429 /* Zero moduli and negative exponents are not considered. */
1436 SETUP(mp_int_init_size(TEMP(0), 2 * um
), last
);
1438 if (c
== b
|| c
== m
)
1440 SETUP(mp_int_init_size(TEMP(1), 2 * um
), last
);
1448 if ((res
= mp_int_mod(a
, m
, TEMP(0))) != MP_OK
)
1451 if ((res
= s_embar(TEMP(0), b
, m
, mu
, s
)) != MP_OK
)
1454 res
= mp_int_copy(s
, c
);
1458 mp_int_clear(TEMP(last
));
1465 /* {{{ mp_int_redux_const(m, c) */
1468 mp_int_redux_const(mp_int m
, mp_int c
)
1470 CHECK(m
!= NULL
&& c
!= NULL
&& m
!= c
);
1472 return s_brmu(c
, m
);
1477 /* {{{ mp_int_invmod(a, m, c) */
1480 mp_int_invmod(mp_int a
, mp_int m
, mp_int c
)
1487 CHECK(a
!= NULL
&& m
!= NULL
&& c
!= NULL
);
1489 if (CMPZ(a
) == 0 || CMPZ(m
) <= 0)
1492 sa
= MP_SIGN(a
); /* need this for the result later */
1494 for (last
= 0; last
< 2; ++last
)
1495 if ((res
= mp_int_init(TEMP(last
))) != MP_OK
)
1498 if ((res
= mp_int_egcd(a
, m
, TEMP(0), TEMP(1), NULL
)) != MP_OK
)
1501 if (mp_int_compare_value(TEMP(0), 1) != 0)
1507 /* It is first necessary to constrain the value to the proper range */
1508 if ((res
= mp_int_mod(TEMP(1), m
, TEMP(1))) != MP_OK
)
1512 * Now, if 'a' was originally negative, the value we have is actually the
1513 * magnitude of the negative representative; to get the positive value we
1514 * have to subtract from the modulus. Otherwise, the value is okay as it
1518 res
= mp_int_sub(m
, TEMP(1), c
);
1520 res
= mp_int_copy(TEMP(1), c
);
1524 mp_int_clear(TEMP(last
));
1531 /* {{{ mp_int_gcd(a, b, c) */
1533 /* Binary GCD algorithm due to Josef Stein, 1961 */
1535 mp_int_gcd(mp_int a
, mp_int b
, mp_int c
)
1545 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
);
1549 if (ca
== 0 && cb
== 0)
1552 return mp_int_abs(b
, c
);
1554 return mp_int_abs(a
, c
);
1556 if ((res
= mp_int_init(&t
)) != MP_OK
)
1558 if ((res
= mp_int_init_copy(&u
, a
)) != MP_OK
)
1560 if ((res
= mp_int_init_copy(&v
, b
)) != MP_OK
)
1563 MP_SIGN(&u
) = MP_ZPOS
;
1564 MP_SIGN(&v
) = MP_ZPOS
;
1566 { /* Divide out common factors of 2 from u and v */
1567 int div2_u
= s_dp2k(&u
),
1568 div2_v
= s_dp2k(&v
);
1570 k
= MIN(div2_u
, div2_v
);
1571 s_qdiv(&u
, (mp_size
) k
);
1572 s_qdiv(&v
, (mp_size
) k
);
1575 if (mp_int_is_odd(&u
))
1577 if ((res
= mp_int_neg(&v
, &t
)) != MP_OK
)
1582 if ((res
= mp_int_copy(&u
, &t
)) != MP_OK
)
1588 s_qdiv(&t
, s_dp2k(&t
));
1592 if ((res
= mp_int_copy(&t
, &u
)) != MP_OK
)
1597 if ((res
= mp_int_neg(&t
, &v
)) != MP_OK
)
1601 if ((res
= mp_int_sub(&u
, &v
, &t
)) != MP_OK
)
1608 if ((res
= mp_int_abs(&u
, c
)) != MP_OK
)
1610 if (!s_qmul(c
, (mp_size
) k
))
1615 V
: mp_int_clear(&u
);
1616 U
: mp_int_clear(&t
);
1623 /* {{{ mp_int_egcd(a, b, c, x, y) */
1625 /* This is the binary GCD algorithm again, but this time we keep track
1626 of the elementary matrix operations as we go, so we can get values
1627 x and y satisfying c = ax + by.
1630 mp_int_egcd(mp_int a
, mp_int b
, mp_int c
,
1640 CHECK(a
!= NULL
&& b
!= NULL
&& c
!= NULL
&&
1641 (x
!= NULL
|| y
!= NULL
));
1645 if (ca
== 0 && cb
== 0)
1649 if ((res
= mp_int_abs(b
, c
)) != MP_OK
)
1652 (void) mp_int_set_value(y
, 1);
1657 if ((res
= mp_int_abs(a
, c
)) != MP_OK
)
1659 (void) mp_int_set_value(x
, 1);
1665 * Initialize temporaries: A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7
1667 for (last
= 0; last
< 4; ++last
)
1669 if ((res
= mp_int_init(TEMP(last
))) != MP_OK
)
1672 TEMP(0)->digits
[0] = 1;
1673 TEMP(3)->digits
[0] = 1;
1675 SETUP(mp_int_init_copy(TEMP(4), a
), last
);
1676 SETUP(mp_int_init_copy(TEMP(5), b
), last
);
1678 /* We will work with absolute values here */
1679 MP_SIGN(TEMP(4)) = MP_ZPOS
;
1680 MP_SIGN(TEMP(5)) = MP_ZPOS
;
1682 { /* Divide out common factors of 2 from u and v */
1683 int div2_u
= s_dp2k(TEMP(4)),
1684 div2_v
= s_dp2k(TEMP(5));
1686 k
= MIN(div2_u
, div2_v
);
1691 SETUP(mp_int_init_copy(TEMP(6), TEMP(4)), last
);
1692 SETUP(mp_int_init_copy(TEMP(7), TEMP(5)), last
);
1696 while (mp_int_is_even(TEMP(4)))
1700 if (mp_int_is_odd(TEMP(0)) || mp_int_is_odd(TEMP(1)))
1702 if ((res
= mp_int_add(TEMP(0), TEMP(7), TEMP(0))) != MP_OK
)
1704 if ((res
= mp_int_sub(TEMP(1), TEMP(6), TEMP(1))) != MP_OK
)
1712 while (mp_int_is_even(TEMP(5)))
1716 if (mp_int_is_odd(TEMP(2)) || mp_int_is_odd(TEMP(3)))
1718 if ((res
= mp_int_add(TEMP(2), TEMP(7), TEMP(2))) != MP_OK
)
1720 if ((res
= mp_int_sub(TEMP(3), TEMP(6), TEMP(3))) != MP_OK
)
1728 if (mp_int_compare(TEMP(4), TEMP(5)) >= 0)
1730 if ((res
= mp_int_sub(TEMP(4), TEMP(5), TEMP(4))) != MP_OK
)
1732 if ((res
= mp_int_sub(TEMP(0), TEMP(2), TEMP(0))) != MP_OK
)
1734 if ((res
= mp_int_sub(TEMP(1), TEMP(3), TEMP(1))) != MP_OK
)
1739 if ((res
= mp_int_sub(TEMP(5), TEMP(4), TEMP(5))) != MP_OK
)
1741 if ((res
= mp_int_sub(TEMP(2), TEMP(0), TEMP(2))) != MP_OK
)
1743 if ((res
= mp_int_sub(TEMP(3), TEMP(1), TEMP(3))) != MP_OK
)
1747 if (CMPZ(TEMP(4)) == 0)
1749 if (x
&& (res
= mp_int_copy(TEMP(2), x
)) != MP_OK
)
1751 if (y
&& (res
= mp_int_copy(TEMP(3), y
)) != MP_OK
)
1755 if (!s_qmul(TEMP(5), k
))
1761 res
= mp_int_copy(TEMP(5), c
);
1770 mp_int_clear(TEMP(last
));
1777 /* {{{ mp_int_divisible_value(a, v) */
1780 mp_int_divisible_value(mp_int a
, int v
)
1784 if (mp_int_div_value(a
, v
, NULL
, &rem
) != MP_OK
)
1792 /* {{{ mp_int_is_pow2(z) */
1795 mp_int_is_pow2(mp_int z
)
1804 /* {{{ mp_int_sqrt(a, c) */
1807 mp_int_sqrt(mp_int a
, mp_int c
)
1809 mp_result res
= MP_OK
;
1813 CHECK(a
!= NULL
&& c
!= NULL
);
1815 /* The square root of a negative value does not exist in the integers. */
1816 if (MP_SIGN(a
) == MP_NEG
)
1819 SETUP(mp_int_init_copy(TEMP(last
), a
), last
);
1820 SETUP(mp_int_init(TEMP(last
)), last
);
1824 if ((res
= mp_int_sqr(TEMP(0), TEMP(1))) != MP_OK
)
1827 if (mp_int_compare_unsigned(a
, TEMP(1)) == 0)
1830 if ((res
= mp_int_copy(a
, TEMP(1))) != MP_OK
)
1832 if ((res
= mp_int_div(TEMP(1), TEMP(0), TEMP(1), NULL
)) != MP_OK
)
1834 if ((res
= mp_int_add(TEMP(0), TEMP(1), TEMP(1))) != MP_OK
)
1836 if ((res
= mp_int_div_pow2(TEMP(1), 1, TEMP(1), NULL
)) != MP_OK
)
1839 if (mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0)
1841 if ((res
= mp_int_sub_value(TEMP(0), 1, TEMP(0))) != MP_OK
)
1843 if (mp_int_compare_unsigned(TEMP(0), TEMP(1)) == 0)
1846 if ((res
= mp_int_copy(TEMP(1), TEMP(0))) != MP_OK
)
1850 res
= mp_int_copy(TEMP(0), c
);
1854 mp_int_clear(TEMP(last
));
1861 /* {{{ mp_int_to_int(z, out) */
1864 mp_int_to_int(mp_int z
, int *out
)
1866 unsigned int uv
= 0;
1873 /* Make sure the value is representable as an int */
1875 if ((sz
== MP_ZPOS
&& mp_int_compare_value(z
, INT_MAX
) > 0) ||
1876 mp_int_compare_value(z
, INT_MIN
) < 0)
1880 dz
= MP_DIGITS(z
) + uz
- 1;
1884 uv
<<= MP_DIGIT_BIT
/ 2;
1885 uv
= (uv
<< (MP_DIGIT_BIT
/ 2)) | *dz
--;
1890 *out
= (sz
== MP_NEG
) ? -(int) uv
: (int) uv
;
1897 /* {{{ mp_int_to_string(z, radix, str, limit) */
1900 mp_int_to_string(mp_int z
, mp_size radix
,
1901 char *str
, int limit
)
1906 CHECK(z
!= NULL
&& str
!= NULL
&& limit
>= 2);
1908 if (radix
< MP_MIN_RADIX
|| radix
> MP_MAX_RADIX
)
1913 *str
++ = s_val2ch(0, mp_flags
& MP_CAP_DIGITS
);
1921 if ((res
= mp_int_init_copy(&tmp
, z
)) != MP_OK
)
1924 if (MP_SIGN(z
) == MP_NEG
)
1931 /* Generate digits in reverse order until finished or limit reached */
1932 for ( /* */ ; limit
> 0; --limit
)
1936 if ((cmp
= CMPZ(&tmp
)) == 0)
1939 d
= s_ddiv(&tmp
, (mp_digit
) radix
);
1940 *str
++ = s_val2ch(d
, mp_flags
& MP_CAP_DIGITS
);
1944 /* Put digits back in correct output order */
1965 /* {{{ mp_int_string_len(z, radix) */
1968 mp_int_string_len(mp_int z
, mp_size radix
)
1974 if (radix
< MP_MIN_RADIX
|| radix
> MP_MAX_RADIX
)
1977 len
= s_outlen(z
, radix
) + 1; /* for terminator */
1979 /* Allow for sign marker on negatives */
1980 if (MP_SIGN(z
) == MP_NEG
)
1988 /* {{{ mp_int_read_string(z, radix, *str) */
1990 /* Read zero-terminated string into z */
1992 mp_int_read_string(mp_int z
, mp_size radix
, const char *str
)
1994 return mp_int_read_cstring(z
, radix
, str
, NULL
);
2000 /* {{{ mp_int_read_cstring(z, radix, *str, **end) */
2003 mp_int_read_cstring(mp_int z
, mp_size radix
, const char *str
, char **end
)
2007 CHECK(z
!= NULL
&& str
!= NULL
);
2009 if (radix
< MP_MIN_RADIX
|| radix
> MP_MAX_RADIX
)
2012 /* Skip leading whitespace */
2013 while (isspace((unsigned char) *str
))
2016 /* Handle leading sign tag (+/-, positive default) */
2020 MP_SIGN(z
) = MP_NEG
;
2024 ++str
; /* fallthrough */
2026 MP_SIGN(z
) = MP_ZPOS
;
2030 /* Skip leading zeroes */
2031 while ((ch
= s_ch2val(*str
, radix
)) == 0)
2034 /* Make sure there is enough space for the value */
2035 if (!s_pad(z
, s_inlen(strlen(str
), radix
)))
2041 while (*str
!= '\0' && ((ch
= s_ch2val(*str
, radix
)) >= 0))
2043 s_dmul(z
, (mp_digit
) radix
);
2044 s_dadd(z
, (mp_digit
) ch
);
2050 /* Override sign for zero, even if negative specified. */
2052 MP_SIGN(z
) = MP_ZPOS
;
2055 *end
= (char *) str
;
2058 * Return a truncation error if the string has unprocessed characters
2059 * remaining, so the caller can tell if the whole string was done
2069 /* {{{ mp_int_count_bits(z) */
2072 mp_int_count_bits(mp_int z
)
2081 if (uz
== 1 && z
->digits
[0] == 0)
2085 nbits
= uz
* MP_DIGIT_BIT
;
2099 /* {{{ mp_int_to_binary(z, buf, limit) */
2102 mp_int_to_binary(mp_int z
, unsigned char *buf
, int limit
)
2104 static const int PAD_FOR_2C
= 1;
2109 CHECK(z
!= NULL
&& buf
!= NULL
);
2111 res
= s_tobin(z
, buf
, &limpos
, PAD_FOR_2C
);
2113 if (MP_SIGN(z
) == MP_NEG
)
2114 s_2comp(buf
, limpos
);
2121 /* {{{ mp_int_read_binary(z, buf, len) */
2124 mp_int_read_binary(mp_int z
, unsigned char *buf
, int len
)
2131 CHECK(z
!= NULL
&& buf
!= NULL
&& len
> 0);
2133 /* Figure out how many digits are needed to represent this value */
2134 need
= ((len
* CHAR_BIT
) + (MP_DIGIT_BIT
- 1)) / MP_DIGIT_BIT
;
2135 if (!s_pad(z
, need
))
2141 * If the high-order bit is set, take the 2's complement before reading
2142 * the value (it will be restored afterward)
2144 if (buf
[0] >> (CHAR_BIT
- 1))
2146 MP_SIGN(z
) = MP_NEG
;
2151 for (tmp
= buf
, i
= len
; i
> 0; --i
, ++tmp
)
2153 s_qmul(z
, (mp_size
) CHAR_BIT
);
2157 /* Restore 2's complement if we took it before */
2158 if (MP_SIGN(z
) == MP_NEG
)
2166 /* {{{ mp_int_binary_len(z) */
2169 mp_int_binary_len(mp_int z
)
2171 mp_result res
= mp_int_count_bits(z
);
2172 int bytes
= mp_int_unsigned_len(z
);
2177 bytes
= (res
+ (CHAR_BIT
- 1)) / CHAR_BIT
;
2180 * If the highest-order bit falls exactly on a byte boundary, we need to
2181 * pad with an extra byte so that the sign will be read correctly when
2182 * reading it back in.
2184 if (bytes
* CHAR_BIT
== res
)
2192 /* {{{ mp_int_to_unsigned(z, buf, limit) */
2195 mp_int_to_unsigned(mp_int z
, unsigned char *buf
, int limit
)
2197 static const int NO_PADDING
= 0;
2199 CHECK(z
!= NULL
&& buf
!= NULL
);
2201 return s_tobin(z
, buf
, &limit
, NO_PADDING
);
2206 /* {{{ mp_int_read_unsigned(z, buf, len) */
2209 mp_int_read_unsigned(mp_int z
, unsigned char *buf
, int len
)
2216 CHECK(z
!= NULL
&& buf
!= NULL
&& len
> 0);
2218 /* Figure out how many digits are needed to represent this value */
2219 need
= ((len
* CHAR_BIT
) + (MP_DIGIT_BIT
- 1)) / MP_DIGIT_BIT
;
2220 if (!s_pad(z
, need
))
2226 for (tmp
= buf
, i
= len
; i
> 0; --i
, ++tmp
)
2228 (void) s_qmul(z
, CHAR_BIT
);
2237 /* {{{ mp_int_unsigned_len(z) */
2240 mp_int_unsigned_len(mp_int z
)
2242 mp_result res
= mp_int_count_bits(z
);
2248 bytes
= (res
+ (CHAR_BIT
- 1)) / CHAR_BIT
;
2255 /* {{{ mp_error_string(res) */
2258 mp_error_string(mp_result res
)
2263 return s_unknown_err
;
2266 for (ix
= 0; ix
< res
&& s_error_msg
[ix
] != NULL
; ++ix
)
2269 if (s_error_msg
[ix
] != NULL
)
2270 return s_error_msg
[ix
];
2272 return s_unknown_err
;
2277 /*------------------------------------------------------------------------*/
2278 /* Private functions for internal use. These make assumptions. */
2280 /* {{{ s_alloc(num) */
2283 s_alloc(mp_size num
)
2285 mp_digit
*out
= px_alloc(num
* sizeof(mp_digit
));
2287 assert(out
!= NULL
); /* for debugging */
2294 /* {{{ s_realloc(old, num) */
2297 s_realloc(mp_digit
* old
, mp_size num
)
2299 mp_digit
*new = px_realloc(old
, num
* sizeof(mp_digit
));
2301 assert(new != NULL
); /* for debugging */
2308 /* {{{ s_free(ptr) */
2320 /* {{{ s_pad(z, min) */
2323 s_pad(mp_int z
, mp_size min
)
2325 if (MP_ALLOC(z
) < min
)
2327 mp_size nsize
= ROUND_PREC(min
);
2328 mp_digit
*tmp
= s_realloc(MP_DIGITS(z
), nsize
);
2334 MP_ALLOC(z
) = nsize
;
2342 /* {{{ s_clamp(z) */
2348 mp_size uz
= MP_USED(z
);
2349 mp_digit
*zd
= MP_DIGITS(z
) + uz
- 1;
2351 while (uz
> 1 && (*zd
-- == 0))
2360 /* {{{ s_fake(z, value, vbuf) */
2363 s_fake(mp_int z
, int value
, mp_digit vbuf
[])
2365 mp_size uv
= (mp_size
) s_vpack(value
, vbuf
);
2368 z
->alloc
= MP_VALUE_DIGITS(value
);
2369 z
->sign
= (value
< 0) ? MP_NEG
: MP_ZPOS
;
2375 /* {{{ s_cdig(da, db, len) */
2378 s_cdig(mp_digit
* da
, mp_digit
* db
, mp_size len
)
2380 mp_digit
*dat
= da
+ len
- 1,
2381 *dbt
= db
+ len
- 1;
2383 for ( /* */ ; len
!= 0; --len
, --dat
, --dbt
)
2387 else if (*dat
< *dbt
)
2396 /* {{{ s_vpack(v, t[]) */
2399 s_vpack(int v
, mp_digit t
[])
2401 unsigned int uv
= (unsigned int) ((v
< 0) ? -v
: v
);
2410 t
[ndig
++] = (mp_digit
) uv
;
2411 uv
>>= MP_DIGIT_BIT
/ 2;
2412 uv
>>= MP_DIGIT_BIT
/ 2;
2421 /* {{{ s_ucmp(a, b) */
2424 s_ucmp(mp_int a
, mp_int b
)
2426 mp_size ua
= MP_USED(a
),
2434 return s_cdig(MP_DIGITS(a
), MP_DIGITS(b
), ua
);
2439 /* {{{ s_vcmp(a, v) */
2442 s_vcmp(mp_int a
, int v
)
2444 mp_digit vdig
[MP_VALUE_DIGITS(v
)];
2446 mp_size ua
= MP_USED(a
);
2448 ndig
= s_vpack(v
, vdig
);
2455 return s_cdig(MP_DIGITS(a
), vdig
, ndig
);
2460 /* {{{ s_uadd(da, db, dc, size_a, size_b) */
2463 s_uadd(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
2464 mp_size size_a
, mp_size size_b
)
2469 /* Insure that da is the longer of the two to simplify later code */
2470 if (size_b
> size_a
)
2472 SWAP(mp_digit
*, da
, db
);
2473 SWAP(mp_size
, size_a
, size_b
);
2476 /* Add corresponding digits until the shorter number runs out */
2477 for (pos
= 0; pos
< size_b
; ++pos
, ++da
, ++db
, ++dc
)
2479 w
= w
+ (mp_word
) * da
+ (mp_word
) * db
;
2480 *dc
= LOWER_HALF(w
);
2484 /* Propagate carries as far as necessary */
2485 for ( /* */ ; pos
< size_a
; ++pos
, ++da
, ++dc
)
2489 *dc
= LOWER_HALF(w
);
2493 /* Return carry out */
2494 return (mp_digit
) w
;
2499 /* {{{ s_usub(da, db, dc, size_a, size_b) */
2502 s_usub(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
2503 mp_size size_a
, mp_size size_b
)
2508 /* We assume that |a| >= |b| so this should definitely hold */
2509 assert(size_a
>= size_b
);
2511 /* Subtract corresponding digits and propagate borrow */
2512 for (pos
= 0; pos
< size_b
; ++pos
, ++da
, ++db
, ++dc
)
2514 w
= ((mp_word
) MP_DIGIT_MAX
+ 1 + /* MP_RADIX */
2515 (mp_word
) * da
) - w
- (mp_word
) * db
;
2517 *dc
= LOWER_HALF(w
);
2518 w
= (UPPER_HALF(w
) == 0);
2521 /* Finish the subtraction for remaining upper digits of da */
2522 for ( /* */ ; pos
< size_a
; ++pos
, ++da
, ++dc
)
2524 w
= ((mp_word
) MP_DIGIT_MAX
+ 1 + /* MP_RADIX */
2525 (mp_word
) * da
) - w
;
2527 *dc
= LOWER_HALF(w
);
2528 w
= (UPPER_HALF(w
) == 0);
2531 /* If there is a borrow out at the end, it violates the precondition */
2537 /* {{{ s_kmul(da, db, dc, size_a, size_b) */
2540 s_kmul(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
2541 mp_size size_a
, mp_size size_b
)
2545 /* Make sure b is the smaller of the two input values */
2546 if (size_b
> size_a
)
2548 SWAP(mp_digit
*, da
, db
);
2549 SWAP(mp_size
, size_a
, size_b
);
2553 * Insure that the bottom is the larger half in an odd-length split; the
2554 * code below relies on this being true.
2556 bot_size
= (size_a
+ 1) / 2;
2559 * If the values are big enough to bother with recursion, use the
2560 * Karatsuba algorithm to compute the product; otherwise use the normal
2561 * multiplication algorithm
2563 if (multiply_threshold
&&
2564 size_a
>= multiply_threshold
&&
2573 mp_digit
*a_top
= da
+ bot_size
;
2574 mp_digit
*b_top
= db
+ bot_size
;
2576 mp_size at_size
= size_a
- bot_size
;
2577 mp_size bt_size
= size_b
- bot_size
;
2578 mp_size buf_size
= 2 * bot_size
;
2581 * Do a single allocation for all three temporary buffers needed; each
2582 * buffer must be big enough to hold the product of two bottom halves,
2583 * and one buffer needs space for the completed product; twice the
2586 if ((t1
= s_alloc(4 * buf_size
)) == NULL
)
2590 ZERO(t1
, 4 * buf_size
);
2593 * t1 and t2 are initially used as temporaries to compute the inner
2594 * product (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0
2596 carry
= s_uadd(da
, a_top
, t1
, bot_size
, at_size
); /* t1 = a1 + a0 */
2597 t1
[bot_size
] = carry
;
2599 carry
= s_uadd(db
, b_top
, t2
, bot_size
, bt_size
); /* t2 = b1 + b0 */
2600 t2
[bot_size
] = carry
;
2602 (void) s_kmul(t1
, t2
, t3
, bot_size
+ 1, bot_size
+ 1); /* t3 = t1 * t2 */
2605 * Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so
2606 * that we're left with only the pieces we want: t3 = a1b0 + a0b1
2608 ZERO(t1
, bot_size
+ 1);
2609 ZERO(t2
, bot_size
+ 1);
2610 (void) s_kmul(da
, db
, t1
, bot_size
, bot_size
); /* t1 = a0 * b0 */
2611 (void) s_kmul(a_top
, b_top
, t2
, at_size
, bt_size
); /* t2 = a1 * b1 */
2613 /* Subtract out t1 and t2 to get the inner product */
2614 s_usub(t3
, t1
, t3
, buf_size
+ 2, buf_size
);
2615 s_usub(t3
, t2
, t3
, buf_size
+ 2, buf_size
);
2617 /* Assemble the output value */
2618 COPY(t1
, dc
, buf_size
);
2619 (void) s_uadd(t3
, dc
+ bot_size
, dc
+ bot_size
,
2620 buf_size
+ 1, buf_size
+ 1);
2622 (void) s_uadd(t2
, dc
+ 2 * bot_size
, dc
+ 2 * bot_size
,
2623 buf_size
, buf_size
);
2625 s_free(t1
); /* note t2 and t3 are just internal pointers
2630 s_umul(da
, db
, dc
, size_a
, size_b
);
2638 /* {{{ s_umul(da, db, dc, size_a, size_b) */
2641 s_umul(mp_digit
* da
, mp_digit
* db
, mp_digit
* dc
,
2642 mp_size size_a
, mp_size size_b
)
2648 for (a
= 0; a
< size_a
; ++a
, ++dc
, ++da
)
2657 for (b
= 0; b
< size_b
; ++b
, ++dbt
, ++dct
)
2659 w
= (mp_word
) * da
* (mp_word
) * dbt
+ w
+ (mp_word
) * dct
;
2661 *dct
= LOWER_HALF(w
);
2665 *dct
= (mp_digit
) w
;
2671 /* {{{ s_ksqr(da, dc, size_a) */
2674 s_ksqr(mp_digit
* da
, mp_digit
* dc
, mp_size size_a
)
2676 if (multiply_threshold
&& size_a
> multiply_threshold
)
2678 mp_size bot_size
= (size_a
+ 1) / 2;
2679 mp_digit
*a_top
= da
+ bot_size
;
2683 mp_size at_size
= size_a
- bot_size
;
2684 mp_size buf_size
= 2 * bot_size
;
2686 if ((t1
= s_alloc(4 * buf_size
)) == NULL
)
2690 ZERO(t1
, 4 * buf_size
);
2692 (void) s_ksqr(da
, t1
, bot_size
); /* t1 = a0 ^ 2 */
2693 (void) s_ksqr(a_top
, t2
, at_size
); /* t2 = a1 ^ 2 */
2695 (void) s_kmul(da
, a_top
, t3
, bot_size
, at_size
); /* t3 = a0 * a1 */
2697 /* Quick multiply t3 by 2, shifting left (can't overflow) */
2700 top
= bot_size
+ at_size
;
2704 for (i
= 0; i
< top
; ++i
)
2707 w
= (w
<< 1) | save
;
2708 t3
[i
] = LOWER_HALF(w
);
2709 save
= UPPER_HALF(w
);
2711 t3
[i
] = LOWER_HALF(save
);
2714 /* Assemble the output value */
2715 COPY(t1
, dc
, 2 * bot_size
);
2716 (void) s_uadd(t3
, dc
+ bot_size
, dc
+ bot_size
,
2717 buf_size
+ 1, buf_size
+ 1);
2719 (void) s_uadd(t2
, dc
+ 2 * bot_size
, dc
+ 2 * bot_size
,
2720 buf_size
, buf_size
);
2722 px_free(t1
); /* note that t2 and t2 are internal pointers
2728 s_usqr(da
, dc
, size_a
);
2736 /* {{{ s_usqr(da, dc, size_a) */
2739 s_usqr(mp_digit
* da
, mp_digit
* dc
, mp_size size_a
)
2745 for (i
= 0; i
< size_a
; ++i
, dc
+= 2, ++da
)
2753 /* Take care of the first digit, no rollover */
2754 w
= (mp_word
) * dat
* (mp_word
) * dat
+ (mp_word
) * dct
;
2755 *dct
= LOWER_HALF(w
);
2760 for (j
= i
+ 1; j
< size_a
; ++j
, ++dat
, ++dct
)
2762 mp_word t
= (mp_word
) * da
* (mp_word
) * dat
;
2763 mp_word u
= w
+ (mp_word
) * dct
,
2766 /* Check if doubling t will overflow a word */
2767 if (HIGH_BIT_SET(t
))
2772 /* Check if adding u to w will overflow a word */
2773 if (ADD_WILL_OVERFLOW(w
, u
))
2778 *dct
= LOWER_HALF(w
);
2782 w
+= MP_DIGIT_MAX
; /* MP_RADIX */
2788 *dct
= (mp_digit
) w
;
2789 while ((w
= UPPER_HALF(w
)) != 0)
2793 *dct
= LOWER_HALF(w
);
2802 /* {{{ s_dadd(a, b) */
2805 s_dadd(mp_int a
, mp_digit b
)
2808 mp_digit
*da
= MP_DIGITS(a
);
2809 mp_size ua
= MP_USED(a
);
2811 w
= (mp_word
) * da
+ b
;
2812 *da
++ = LOWER_HALF(w
);
2815 for (ua
-= 1; ua
> 0; --ua
, ++da
)
2817 w
= (mp_word
) * da
+ w
;
2819 *da
= LOWER_HALF(w
);
2832 /* {{{ s_dmul(a, b) */
2835 s_dmul(mp_int a
, mp_digit b
)
2838 mp_digit
*da
= MP_DIGITS(a
);
2839 mp_size ua
= MP_USED(a
);
2843 w
= (mp_word
) * da
* b
+ w
;
2844 *da
++ = LOWER_HALF(w
);
2858 /* {{{ s_dbmul(da, b, dc, size_a) */
2861 s_dbmul(mp_digit
* da
, mp_digit b
, mp_digit
* dc
, mp_size size_a
)
2867 w
= (mp_word
) * da
++ * (mp_word
) b
+ w
;
2869 *dc
++ = LOWER_HALF(w
);
2875 *dc
= LOWER_HALF(w
);
2880 /* {{{ s_ddiv(da, d, dc, size_a) */
2883 s_ddiv(mp_int a
, mp_digit b
)
2887 mp_size ua
= MP_USED(a
);
2888 mp_digit
*da
= MP_DIGITS(a
) + ua
- 1;
2890 for ( /* */ ; ua
> 0; --ua
, --da
)
2892 w
= (w
<< MP_DIGIT_BIT
) | *da
;
2904 *da
= (mp_digit
) qdigit
;
2908 return (mp_digit
) w
;
2913 /* {{{ s_qdiv(z, p2) */
2916 s_qdiv(mp_int z
, mp_size p2
)
2918 mp_size ndig
= p2
/ MP_DIGIT_BIT
,
2919 nbits
= p2
% MP_DIGIT_BIT
;
2920 mp_size uz
= MP_USED(z
);
2937 for (mark
= ndig
; mark
< uz
; ++mark
)
2940 MP_USED(z
) = uz
- ndig
;
2948 mp_size up
= MP_DIGIT_BIT
- nbits
;
2951 dz
= MP_DIGITS(z
) + uz
- 1;
2953 for ( /* */ ; uz
> 0; --uz
, --dz
)
2957 *dz
= (*dz
>> nbits
) | (d
<< up
);
2964 if (MP_USED(z
) == 1 && z
->digits
[0] == 0)
2965 MP_SIGN(z
) = MP_ZPOS
;
2970 /* {{{ s_qmod(z, p2) */
2973 s_qmod(mp_int z
, mp_size p2
)
2975 mp_size start
= p2
/ MP_DIGIT_BIT
+ 1,
2976 rest
= p2
% MP_DIGIT_BIT
;
2977 mp_size uz
= MP_USED(z
);
2978 mp_digit mask
= (1 << rest
) - 1;
2983 z
->digits
[start
- 1] &= mask
;
2990 /* {{{ s_qmul(z, p2) */
2993 s_qmul(mp_int z
, mp_size p2
)
3008 need
= p2
/ MP_DIGIT_BIT
;
3009 rest
= p2
% MP_DIGIT_BIT
;
3012 * Figure out if we need an extra digit at the top end; this occurs if the
3013 * topmost `rest' bits of the high-order digit of z are not zero, meaning
3014 * they will be shifted off the end if not preserved
3019 mp_digit
*dz
= MP_DIGITS(z
) + uz
- 1;
3021 if ((*dz
>> (MP_DIGIT_BIT
- rest
)) != 0)
3025 if (!s_pad(z
, uz
+ need
+ extra
))
3029 * If we need to shift by whole digits, do that in one pass, then to back
3030 * and shift by partial digits.
3034 from
= MP_DIGITS(z
) + uz
- 1;
3037 for (i
= 0; i
< uz
; ++i
)
3040 ZERO(MP_DIGITS(z
), need
);
3047 for (i
= need
, from
= MP_DIGITS(z
) + need
; i
< uz
; ++i
, ++from
)
3049 mp_digit save
= *from
;
3051 *from
= (*from
<< rest
) | (d
>> (MP_DIGIT_BIT
- rest
));
3055 d
>>= (MP_DIGIT_BIT
- rest
);
3071 /* {{{ s_qsub(z, p2) */
3073 /* Subtract |z| from 2^p2, assuming 2^p2 > |z|, and set z to be positive */
3075 s_qsub(mp_int z
, mp_size p2
)
3077 mp_digit hi
= (1 << (p2
% MP_DIGIT_BIT
)),
3079 mp_size tdig
= (p2
/ MP_DIGIT_BIT
),
3083 if (!s_pad(z
, tdig
+ 1))
3086 for (pos
= 0, zp
= MP_DIGITS(z
); pos
< tdig
; ++pos
, ++zp
)
3088 w
= ((mp_word
) MP_DIGIT_MAX
+ 1) - w
- (mp_word
) * zp
;
3090 *zp
= LOWER_HALF(w
);
3091 w
= UPPER_HALF(w
) ? 0 : 1;
3094 w
= ((mp_word
) MP_DIGIT_MAX
+ 1 + hi
) - w
- (mp_word
) * zp
;
3095 *zp
= LOWER_HALF(w
);
3097 assert(UPPER_HALF(w
) != 0); /* no borrow out should be possible */
3099 MP_SIGN(z
) = MP_ZPOS
;
3113 mp_digit
*dp
= MP_DIGITS(z
),
3116 if (MP_USED(z
) == 1 && *dp
== 0)
3126 while ((d
& 1) == 0)
3142 mp_size uz
= MP_USED(z
),
3144 mp_digit
*dz
= MP_DIGITS(z
),
3169 /* {{{ s_2expt(z, k) */
3172 s_2expt(mp_int z
, int k
)
3178 ndig
= (k
+ MP_DIGIT_BIT
) / MP_DIGIT_BIT
;
3179 rest
= k
% MP_DIGIT_BIT
;
3181 if (!s_pad(z
, ndig
))
3186 *(dz
+ ndig
- 1) = (1 << rest
);
3194 /* {{{ s_norm(a, b) */
3197 s_norm(mp_int a
, mp_int b
)
3199 mp_digit d
= b
->digits
[MP_USED(b
) - 1];
3202 while (d
< (mp_digit
) ((mp_digit
) 1 << (MP_DIGIT_BIT
- 1)))
3203 { /* d < (MP_RADIX / 2) */
3208 /* These multiplications can't fail */
3211 (void) s_qmul(a
, (mp_size
) k
);
3212 (void) s_qmul(b
, (mp_size
) k
);
3220 /* {{{ s_brmu(z, m) */
3223 s_brmu(mp_int z
, mp_int m
)
3225 mp_size um
= MP_USED(m
) * 2;
3230 s_2expt(z
, MP_DIGIT_BIT
* um
);
3231 return mp_int_div(z
, m
, z
, NULL
);
3236 /* {{{ s_reduce(x, m, mu, q1, q2) */
3239 s_reduce(mp_int x
, mp_int m
, mp_int mu
, mp_int q1
, mp_int q2
)
3241 mp_size um
= MP_USED(m
),
3245 umb_p1
= (um
+ 1) * MP_DIGIT_BIT
;
3246 umb_m1
= (um
- 1) * MP_DIGIT_BIT
;
3248 if (mp_int_copy(x
, q1
) != MP_OK
)
3251 /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */
3256 /* Set x = x mod b^(k+1) */
3260 * Now, q is a guess for the quotient a / m. Compute x - q * m mod
3261 * b^(k+1), replacing x. This may be off by a factor of 2m, but no more
3266 (void) mp_int_sub(x
, q1
, x
); /* can't fail */
3269 * The result may be < 0; if it is, add b^(k+1) to pin it in the proper
3272 if ((CMPZ(x
) < 0) && !s_qsub(x
, umb_p1
))
3276 * If x > m, we need to back it off until it is in range. This will be
3277 * required at most twice.
3279 if (mp_int_compare(x
, m
) >= 0)
3280 (void) mp_int_sub(x
, m
, x
);
3281 if (mp_int_compare(x
, m
) >= 0)
3282 (void) mp_int_sub(x
, m
, x
);
3284 /* At this point, x has been properly reduced. */
3290 /* {{{ s_embar(a, b, m, mu, c) */
3292 /* Perform modular exponentiation using Barrett's method, where mu is
3293 the reduction constant for m. Assumes a < m, b > 0. */
3295 s_embar(mp_int a
, mp_int b
, mp_int m
, mp_int mu
, mp_int c
)
3307 dbt
= db
+ MP_USED(b
) - 1;
3310 SETUP(mp_int_init_size(TEMP(last
), 2 * umu
), last
);
3312 (void) mp_int_set_value(c
, 1);
3314 /* Take care of low-order digits */
3319 for (d
= *db
, i
= MP_DIGIT_BIT
; i
> 0; --i
, d
>>= 1)
3323 /* The use of a second temporary avoids allocation */
3324 UMUL(c
, a
, TEMP(0));
3325 if (!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2)))
3330 mp_int_copy(TEMP(0), c
);
3335 assert(MP_SIGN(TEMP(0)) == MP_ZPOS
);
3336 if (!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2)))
3341 assert(MP_SIGN(TEMP(0)) == MP_ZPOS
);
3342 mp_int_copy(TEMP(0), a
);
3350 /* Take care of highest-order digit */
3356 UMUL(c
, a
, TEMP(0));
3357 if (!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2)))
3362 mp_int_copy(TEMP(0), c
);
3370 if (!s_reduce(TEMP(0), m
, mu
, TEMP(1), TEMP(2)))
3375 (void) mp_int_copy(TEMP(0), a
);
3380 mp_int_clear(TEMP(last
));
3387 /* {{{ s_udiv(a, b) */
3389 /* Precondition: a >= b and b > 0
3390 Postcondition: a' = a / b, b' = a % b
3393 s_udiv(mp_int a
, mp_int b
)
3403 mp_result res
= MP_OK
;
3407 /* Force signs to positive */
3408 MP_SIGN(a
) = MP_ZPOS
;
3409 MP_SIGN(b
) = MP_ZPOS
;
3411 /* Normalize, per Knuth */
3416 btop
= b
->digits
[ub
- 1];
3417 if ((res
= mp_int_init_size(&q
, ua
)) != MP_OK
)
3419 if ((res
= mp_int_init_size(&t
, ua
+ 1)) != MP_OK
)
3423 r
.digits
= da
+ ua
- 1; /* The contents of r are shared with a */
3426 r
.alloc
= MP_ALLOC(a
);
3427 ZERO(t
.digits
, t
.alloc
);
3429 /* Solve for quotient digits, store in q.digits in reverse order */
3430 while (r
.digits
>= da
)
3432 assert(qpos
<= q
.alloc
);
3434 if (s_ucmp(b
, &r
) > 0)
3440 q
.digits
[qpos
++] = 0;
3446 mp_word pfx
= r
.digits
[r
.used
- 1];
3449 if (r
.used
> 1 && (pfx
< btop
|| r
.digits
[r
.used
- 2] == 0))
3451 pfx
<<= MP_DIGIT_BIT
/ 2;
3452 pfx
<<= MP_DIGIT_BIT
/ 2;
3453 pfx
|= r
.digits
[r
.used
- 2];
3456 qdigit
= pfx
/ btop
;
3457 if (qdigit
> MP_DIGIT_MAX
)
3460 s_dbmul(MP_DIGITS(b
), (mp_digit
) qdigit
, t
.digits
, ub
);
3463 while (s_ucmp(&t
, &r
) > 0)
3466 (void) mp_int_sub(&t
, b
, &t
); /* cannot fail */
3469 s_usub(r
.digits
, t
.digits
, r
.digits
, r
.used
, t
.used
);
3472 q
.digits
[qpos
++] = (mp_digit
) qdigit
;
3473 ZERO(t
.digits
, t
.used
);
3478 /* Put quotient digits in the correct order, and discard extra zeroes */
3480 REV(mp_digit
, q
.digits
, qpos
);
3483 /* Denormalize the remainder */
3488 mp_int_copy(a
, b
); /* ok: 0 <= r < b */
3489 mp_int_copy(&q
, a
); /* ok: q <= a */
3499 /* {{{ s_outlen(z, r) */
3501 /* Precondition: 2 <= r < 64 */
3503 s_outlen(mp_int z
, mp_size r
)
3508 bits
= mp_int_count_bits(z
);
3509 raw
= (double) bits
*s_log2
[r
];
3511 return (int) (raw
+ 0.999999);
3516 /* {{{ s_inlen(len, r) */
3519 s_inlen(int len
, mp_size r
)
3521 double raw
= (double) len
/ s_log2
[r
];
3522 mp_size bits
= (mp_size
) (raw
+ 0.5);
3524 return (mp_size
) ((bits
+ (MP_DIGIT_BIT
- 1)) / MP_DIGIT_BIT
);
3529 /* {{{ s_ch2val(c, r) */
3532 s_ch2val(char c
, int r
)
3536 if (isdigit((unsigned char) c
))
3538 else if (r
> 10 && isalpha((unsigned char) c
))
3539 out
= toupper((unsigned char) c
) - 'A' + 10;
3543 return (out
>= r
) ? -1 : out
;
3548 /* {{{ s_val2ch(v, caps) */
3551 s_val2ch(int v
, int caps
)
3559 char out
= (v
- 10) + 'a';
3562 return toupper((unsigned char) out
);
3570 /* {{{ s_2comp(buf, len) */
3573 s_2comp(unsigned char *buf
, int len
)
3576 unsigned short s
= 1;
3578 for (i
= len
- 1; i
>= 0; --i
)
3580 unsigned char c
= ~buf
[i
];
3589 /* last carry out is ignored */
3594 /* {{{ s_tobin(z, buf, *limpos) */
3597 s_tobin(mp_int z
, unsigned char *buf
, int *limpos
, int pad
)
3606 while (uz
> 0 && pos
< limit
)
3611 for (i
= sizeof(mp_digit
); i
> 0 && pos
< limit
; --i
)
3613 buf
[pos
++] = (unsigned char) d
;
3616 /* Don't write leading zeroes */
3617 if (d
== 0 && uz
== 1)
3618 i
= 0; /* exit loop without signaling truncation */
3621 /* Detect truncation (loop exited with pos >= limit) */
3628 if (pad
!= 0 && (buf
[pos
- 1] >> (CHAR_BIT
- 1)))
3636 /* Digits are in reverse order, fix that */
3637 REV(unsigned char, buf
, pos
);
3639 /* Return the number of bytes actually written */
3642 return (uz
== 0) ? MP_OK
: MP_TRUNC
;
3647 /* {{{ s_print(tag, z) */
3651 s_print(char *tag
, mp_int z
)
3655 fprintf(stderr
, "%s: %c ", tag
,
3656 (MP_SIGN(z
) == MP_NEG
) ? '-' : '+');
3658 for (i
= MP_USED(z
) - 1; i
>= 0; --i
)
3659 fprintf(stderr
, "%0*X", (int) (MP_DIGIT_BIT
/ 4), z
->digits
[i
]);
3661 fputc('\n', stderr
);
3666 s_print_buf(char *tag
, mp_digit
* buf
, mp_size num
)
3670 fprintf(stderr
, "%s: ", tag
);
3672 for (i
= num
- 1; i
>= 0; --i
)
3673 fprintf(stderr
, "%0*X", (int) (MP_DIGIT_BIT
/ 4), buf
[i
]);
3675 fputc('\n', stderr
);
3681 /* HERE THERE BE DRAGONS */