/Users/buildslave/jenkins/sharedspace/clang-stage2-coverage-R@2/llvm/tools/polly/lib/External/isl/imath/imath.c
Line | Count | Source (jump to first uncovered line) |
1 | | /* |
2 | | Name: imath.c |
3 | | Purpose: Arbitrary precision integer arithmetic routines. |
4 | | Author: M. J. Fromberger <http://spinning-yarns.org/michael/> |
5 | | |
6 | | Copyright (C) 2002-2007 Michael J. Fromberger, All Rights Reserved. |
7 | | |
8 | | Permission is hereby granted, free of charge, to any person obtaining a copy |
9 | | of this software and associated documentation files (the "Software"), to deal |
10 | | in the Software without restriction, including without limitation the rights |
11 | | to use, copy, modify, merge, publish, distribute, sublicense, and/or sell |
12 | | copies of the Software, and to permit persons to whom the Software is |
13 | | furnished to do so, subject to the following conditions: |
14 | | |
15 | | The above copyright notice and this permission notice shall be included in |
16 | | all copies or substantial portions of the Software. |
17 | | |
18 | | THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR |
19 | | IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, |
20 | | FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE |
21 | | AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER |
22 | | LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, |
23 | | OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE |
24 | | SOFTWARE. |
25 | | */ |
26 | | |
27 | | #include "imath.h" |
28 | | |
29 | | #if DEBUG |
30 | | #include <stdio.h> |
31 | | #endif |
32 | | |
33 | | #include <stdlib.h> |
34 | | #include <string.h> |
35 | | #include <ctype.h> |
36 | | |
37 | | #include <assert.h> |
38 | | |
39 | | #if DEBUG |
40 | | #define STATIC /* public */ |
41 | | #else |
42 | | #define STATIC static |
43 | | #endif |
44 | | |
45 | | const mp_result MP_OK = 0; /* no error, all is well */ |
46 | | const mp_result MP_FALSE = 0; /* boolean false */ |
47 | | const mp_result MP_TRUE = -1; /* boolean true */ |
48 | | const mp_result MP_MEMORY = -2; /* out of memory */ |
49 | | const mp_result MP_RANGE = -3; /* argument out of range */ |
50 | | const mp_result MP_UNDEF = -4; /* result undefined */ |
51 | | const mp_result MP_TRUNC = -5; /* output truncated */ |
52 | | const mp_result MP_BADARG = -6; /* invalid null argument */ |
53 | | const mp_result MP_MINERR = -6; |
54 | | |
55 | | const mp_sign MP_NEG = 1; /* value is strictly negative */ |
56 | | const mp_sign MP_ZPOS = 0; /* value is non-negative */ |
57 | | |
58 | | STATIC const char *s_unknown_err = "unknown result code"; |
59 | | STATIC const char *s_error_msg[] = { |
60 | | "error code 0", |
61 | | "boolean true", |
62 | | "out of memory", |
63 | | "argument out of range", |
64 | | "result undefined", |
65 | | "output truncated", |
66 | | "invalid argument", |
67 | | NULL |
68 | | }; |
69 | | |
70 | | /* Argument checking macros |
71 | | Use CHECK() where a return value is required; NRCHECK() elsewhere */ |
72 | 439M | #define CHECK(TEST) assert(TEST) |
73 | 117M | #define NRCHECK(TEST) assert(TEST) |
74 | | |
75 | | /* The ith entry of this table gives the value of log_i(2). |
76 | | |
77 | | An integer value n requires ceil(log_i(n)) digits to be represented |
78 | | in base i. Since it is easy to compute lg(n), by counting bits, we |
79 | | can compute log_i(n) = lg(n) * log_i(2). |
80 | | |
81 | | The use of this table eliminates a dependency upon linkage against |
82 | | the standard math libraries. |
83 | | |
84 | | If MP_MAX_RADIX is increased, this table should be expanded too. |
85 | | */ |
86 | | STATIC const double s_log2[] = { |
87 | | 0.000000000, 0.000000000, 1.000000000, 0.630929754, /* (D)(D) 2 3 */ |
88 | | 0.500000000, 0.430676558, 0.386852807, 0.356207187, /* 4 5 6 7 */ |
89 | | 0.333333333, 0.315464877, 0.301029996, 0.289064826, /* 8 9 10 11 */ |
90 | | 0.278942946, 0.270238154, 0.262649535, 0.255958025, /* 12 13 14 15 */ |
91 | | 0.250000000, 0.244650542, 0.239812467, 0.235408913, /* 16 17 18 19 */ |
92 | | 0.231378213, 0.227670249, 0.224243824, 0.221064729, /* 20 21 22 23 */ |
93 | | 0.218104292, 0.215338279, 0.212746054, 0.210309918, /* 24 25 26 27 */ |
94 | | 0.208014598, 0.205846832, 0.203795047, 0.201849087, /* 28 29 30 31 */ |
95 | | 0.200000000, 0.198239863, 0.196561632, 0.194959022, /* 32 33 34 35 */ |
96 | | 0.193426404, /* 36 */ |
97 | | }; |
98 | | |
99 | | |
100 | | |
101 | | /* Return the number of digits needed to represent a static value */ |
102 | 49.2M | #define MP_VALUE_DIGITS(V) \ |
103 | 49.2M | ((sizeof(V)+(sizeof(mp_digit)-1))/sizeof(mp_digit)) |
104 | | |
105 | | /* Round precision P to nearest word boundary */ |
106 | 21.7M | #define ROUND_PREC(P) ((mp_size)(2*(((P)+1)/2))) |
107 | | |
108 | | /* Set array P of S digits to zero */ |
109 | 10.3M | #define ZERO(P, S) \ |
110 | 10.3M | do{ \ |
111 | 10.3M | mp_size i__ = (S) * sizeof(mp_digit); \ |
112 | 10.3M | mp_digit *p__ = (P); \ |
113 | 10.3M | memset(p__, 0, i__); \ |
114 | 10.3M | } while(0) |
115 | | |
116 | | /* Copy S digits from array P to array Q */ |
117 | 132M | #define COPY(P, Q, S) \ |
118 | 119M | do{ \ |
119 | 119M | mp_size i__ = (S) * sizeof(mp_digit); \ |
120 | 119M | mp_digit *p__ = (P), *q__ = (Q); \ |
121 | 119M | memcpy(q__, p__, i__); \ |
122 | 132M | } while(0) |
123 | | |
124 | | /* Reverse N elements of type T in array A */ |
125 | 0 | #define REV(T, A, N) \ |
126 | 0 | do{ \ |
127 | 0 | T *u_ = (A), *v_ = u_ + (N) - 1; \ |
128 | 0 | while (u_ < v_0 ) { \ |
129 | 0 | T xch = *u_; \ |
130 | 0 | *u_++ = *v_; \ |
131 | 0 | *v_-- = xch; \ |
132 | 0 | } \ |
133 | 0 | } while(0) |
134 | | |
135 | 182M | #define CLAMP(Z) \ |
136 | 182M | do{ \ |
137 | 182M | mp_int z_ = (Z); \ |
138 | 182M | mp_size uz_ = MP_USED(z_); \ |
139 | 182M | mp_digit *dz_ = MP_DIGITS(z_) + uz_ -1; \ |
140 | 215M | while (uz_ > 1 && 215M (*dz_-- == 0)147M ) \ |
141 | 33.5M | --uz_; \ |
142 | 182M | MP_USED(z_) = uz_; \ |
143 | 182M | } while(0) |
144 | | |
145 | | /* Select min/max. Do not provide expressions for which multiple |
146 | | evaluation would be problematic, e.g. x++ */ |
147 | 2.76M | #define MIN(A, B) ((B)<(A)?2.76M (B)544k :(A)2.21M ) |
148 | 106M | #define MAX(A, B) ((B)>(A)?106M (B)47.1M :(A)59.1M ) |
149 | | |
150 | | /* Exchange lvalues A and B of type T, e.g. |
151 | | SWAP(int, x, y) where x and y are variables of type int. */ |
152 | 7.71M | #define SWAP(T, A, B) \ |
153 | 7.71M | do{ \ |
154 | 7.71M | T t_ = (A); \ |
155 | 7.71M | A = (B); \ |
156 | 7.71M | B = t_; \ |
157 | 7.71M | } while(0) |
158 | | |
159 | | /* Used to set up and access simple temp stacks within functions. */ |
160 | | #define DECLARE_TEMP(N) \ |
161 | 4.20M | mpz_t temp[(N)]; \ |
162 | 4.20M | int last__ = 0 |
163 | | #define CLEANUP_TEMP() \ |
164 | 5.36M | CLEANUP: \ |
165 | 4.39M | while (--last__ >= 0) \ |
166 | 5.36M | mp_int_clear(1.70M TEMP1.70M (last__)) |
167 | 3.41M | #define TEMP(K) (temp + (K)) |
168 | 1.70M | #define LAST_TEMP() TEMP(last__) |
169 | 1.70M | #define SETUP(E) \ |
170 | 1.70M | do{ \ |
171 | 1.70M | if ((res = (E)) != MP_OK) \ |
172 | 0 | goto CLEANUP; \ |
173 | 1.70M | ++(last__); \ |
174 | 1.70M | } while(0) |
175 | | |
176 | | /* Compare value to zero. */ |
177 | 242M | #define CMPZ(Z) \ |
178 | 242M | (((Z)->used==1&&(Z)->digits[0]==092.2M )?04.60M :((Z)->sign==MP_NEG)?238M -1203M :134.7M ) |
179 | | |
180 | | /* Multiply X by Y into Z, ignoring signs. Requires that Z have |
181 | | enough storage preallocated to hold the result. */ |
182 | 0 | #define UMUL(X, Y, Z) \ |
183 | 0 | do{ \ |
184 | 0 | mp_size ua_ = MP_USED0 (X), ub_ = MP_USED0 (Y); \ |
185 | 0 | mp_size o_ = ua_ + ub_; \ |
186 | 0 | ZERO(MP_DIGITS(Z), o_); \ |
187 | 0 | (void) s_kmul(MP_DIGITS0 (X), MP_DIGITS0 (Y), MP_DIGITS0 (Z), ua_, ub_); \ |
188 | 0 | MP_USED(Z) = o_; \ |
189 | 0 | CLAMP(Z); \ |
190 | 0 | } while(0) |
191 | | |
192 | | /* Square X into Z. Requires that Z have enough storage to hold the |
193 | | result. */ |
194 | 0 | #define USQR(X, Z) \ |
195 | 0 | do{ \ |
196 | 0 | mp_size ua_ = MP_USED(X), o_ = ua_ + ua_; \ |
197 | 0 | ZERO(MP_DIGITS(Z), o_); \ |
198 | 0 | (void) s_ksqr(MP_DIGITS0 (X), MP_DIGITS0 (Z), ua_); \ |
199 | 0 | MP_USED(Z) = o_; \ |
200 | 0 | CLAMP(Z); \ |
201 | 0 | } while(0) |
202 | | |
203 | 272M | #define UPPER_HALF(W) ((mp_word)((W) >> 272M MP_DIGIT_BIT272M )) |
204 | 272M | #define LOWER_HALF(W) ((mp_digit)(W)) |
205 | 0 | #define HIGH_BIT_SET(W) ((W) >> (0 MP_WORD_BIT0 - 1)) |
206 | 0 | #define ADD_WILL_OVERFLOW(W, V) ((0 MP_WORD_MAX0 - (V)) < (W)) |
207 | | |
208 | | |
209 | | |
210 | | /* Default number of digits allocated to a new mp_int */ |
211 | | #if IMATH_TEST |
212 | | mp_size default_precision = MP_DEFAULT_PREC; |
213 | | #else |
214 | | STATIC const mp_size default_precision = MP_DEFAULT_PREC; |
215 | | #endif |
216 | | |
217 | | /* Minimum number of digits to invoke recursive multiply */ |
218 | | #if IMATH_TEST |
219 | | mp_size multiply_threshold = MP_MULT_THRESH; |
220 | | #else |
221 | | STATIC const mp_size multiply_threshold = MP_MULT_THRESH; |
222 | | #endif |
223 | | |
224 | | /* Allocate a buffer of (at least) num digits, or return |
225 | | NULL if that couldn't be done. */ |
226 | | STATIC mp_digit *s_alloc(mp_size num); |
227 | | |
228 | | /* Release a buffer of digits allocated by s_alloc(). */ |
229 | | STATIC void s_free(void *ptr); |
230 | | |
231 | | /* Insure that z has at least min digits allocated, resizing if |
232 | | necessary. Returns true if successful, false if out of memory. */ |
233 | | STATIC int s_pad(mp_int z, mp_size min); |
234 | | |
235 | | /* Fill in a "fake" mp_int on the stack with a given value */ |
236 | | STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[]); |
237 | | STATIC void s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[]); |
238 | | |
239 | | /* Compare two runs of digits of given length, returns <0, 0, >0 */ |
240 | | STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len); |
241 | | |
242 | | /* Pack the unsigned digits of v into array t */ |
243 | | STATIC int s_uvpack(mp_usmall v, mp_digit t[]); |
244 | | |
245 | | /* Compare magnitudes of a and b, returns <0, 0, >0 */ |
246 | | STATIC int s_ucmp(mp_int a, mp_int b); |
247 | | |
248 | | /* Compare magnitudes of a and v, returns <0, 0, >0 */ |
249 | | STATIC int s_vcmp(mp_int a, mp_small v); |
250 | | STATIC int s_uvcmp(mp_int a, mp_usmall uv); |
251 | | |
252 | | /* Unsigned magnitude addition; assumes dc is big enough. |
253 | | Carry out is returned (no memory allocated). */ |
254 | | STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, |
255 | | mp_size size_a, mp_size size_b); |
256 | | |
257 | | /* Unsigned magnitude subtraction. Assumes dc is big enough. */ |
258 | | STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, |
259 | | mp_size size_a, mp_size size_b); |
260 | | |
261 | | /* Unsigned recursive multiplication. Assumes dc is big enough. */ |
262 | | STATIC int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc, |
263 | | mp_size size_a, mp_size size_b); |
264 | | |
265 | | /* Unsigned magnitude multiplication. Assumes dc is big enough. */ |
266 | | STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, |
267 | | mp_size size_a, mp_size size_b); |
268 | | |
269 | | /* Unsigned recursive squaring. Assumes dc is big enough. */ |
270 | | STATIC int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a); |
271 | | |
272 | | /* Unsigned magnitude squaring. Assumes dc is big enough. */ |
273 | | STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a); |
274 | | |
275 | | /* Single digit addition. Assumes a is big enough. */ |
276 | | STATIC void s_dadd(mp_int a, mp_digit b); |
277 | | |
278 | | /* Single digit multiplication. Assumes a is big enough. */ |
279 | | STATIC void s_dmul(mp_int a, mp_digit b); |
280 | | |
281 | | /* Single digit multiplication on buffers; assumes dc is big enough. */ |
282 | | STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, |
283 | | mp_size size_a); |
284 | | |
285 | | /* Single digit division. Replaces a with the quotient, |
286 | | returns the remainder. */ |
287 | | STATIC mp_digit s_ddiv(mp_int a, mp_digit b); |
288 | | |
289 | | /* Quick division by a power of 2, replaces z (no allocation) */ |
290 | | STATIC void s_qdiv(mp_int z, mp_size p2); |
291 | | |
292 | | /* Quick remainder by a power of 2, replaces z (no allocation) */ |
293 | | STATIC void s_qmod(mp_int z, mp_size p2); |
294 | | |
295 | | /* Quick multiplication by a power of 2, replaces z. |
296 | | Allocates if necessary; returns false in case this fails. */ |
297 | | STATIC int s_qmul(mp_int z, mp_size p2); |
298 | | |
299 | | /* Quick subtraction from a power of 2, replaces z. |
300 | | Allocates if necessary; returns false in case this fails. */ |
301 | | STATIC int s_qsub(mp_int z, mp_size p2); |
302 | | |
303 | | /* Return maximum k such that 2^k divides z. */ |
304 | | STATIC int s_dp2k(mp_int z); |
305 | | |
306 | | /* Return k >= 0 such that z = 2^k, or -1 if there is no such k. */ |
307 | | STATIC int s_isp2(mp_int z); |
308 | | |
309 | | /* Set z to 2^k. May allocate; returns false in case this fails. */ |
310 | | STATIC int s_2expt(mp_int z, mp_small k); |
311 | | |
312 | | /* Normalize a and b for division, returns normalization constant */ |
313 | | STATIC int s_norm(mp_int a, mp_int b); |
314 | | |
315 | | /* Compute constant mu for Barrett reduction, given modulus m, result |
316 | | replaces z, m is untouched. */ |
317 | | STATIC mp_result s_brmu(mp_int z, mp_int m); |
318 | | |
319 | | /* Reduce a modulo m, using Barrett's algorithm. */ |
320 | | STATIC int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2); |
321 | | |
322 | | /* Modular exponentiation, using Barrett reduction */ |
323 | | STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c); |
324 | | |
325 | | /* Unsigned magnitude division. Assumes |a| > |b|. Allocates temporaries; |
326 | | overwrites a with quotient, b with remainder. */ |
327 | | STATIC mp_result s_udiv_knuth(mp_int a, mp_int b); |
328 | | |
329 | | /* Compute the number of digits in radix r required to represent the given |
330 | | value. Does not account for sign flags, terminators, etc. */ |
331 | | STATIC int s_outlen(mp_int z, mp_size r); |
332 | | |
333 | | /* Guess how many digits of precision will be needed to represent a radix r |
334 | | value of the specified number of digits. Returns a value guaranteed to be |
335 | | no smaller than the actual number required. */ |
336 | | STATIC mp_size s_inlen(int len, mp_size r); |
337 | | |
338 | | /* Convert a character to a digit value in radix r, or |
339 | | -1 if out of range */ |
340 | | STATIC int s_ch2val(char c, int r); |
341 | | |
342 | | /* Convert a digit value to a character */ |
343 | | STATIC char s_val2ch(int v, int caps); |
344 | | |
345 | | /* Take 2's complement of a buffer in place */ |
346 | | STATIC void s_2comp(unsigned char *buf, int len); |
347 | | |
348 | | /* Convert a value to binary, ignoring sign. On input, *limpos is the bound on |
349 | | how many bytes should be written to buf; on output, *limpos is set to the |
350 | | number of bytes actually written. */ |
351 | | STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad); |
352 | | |
353 | | #if DEBUG |
354 | | /* Dump a representation of the mp_int to standard output */ |
355 | | void s_print(char *tag, mp_int z); |
356 | | void s_print_buf(char *tag, mp_digit *buf, mp_size num); |
357 | | #endif |
358 | | |
359 | | mp_result mp_int_init(mp_int z) |
360 | 38.6M | { |
361 | 38.6M | if (z == NULL) |
362 | 0 | return MP_BADARG; |
363 | 38.6M | |
364 | 38.6M | z->single = 0; |
365 | 38.6M | z->digits = &(z->single); |
366 | 38.6M | z->alloc = 1; |
367 | 38.6M | z->used = 1; |
368 | 38.6M | z->sign = MP_ZPOS; |
369 | 38.6M | |
370 | 38.6M | return MP_OK; |
371 | 38.6M | } |
372 | | |
373 | | mp_int mp_int_alloc(void) |
374 | 31.6M | { |
375 | 31.6M | mp_int out = malloc(sizeof(mpz_t)); |
376 | 31.6M | |
377 | 31.6M | if (out != NULL) |
378 | 31.6M | mp_int_init(out); |
379 | 31.6M | |
380 | 31.6M | return out; |
381 | 31.6M | } |
382 | | |
383 | | mp_result mp_int_init_size(mp_int z, mp_size prec) |
384 | 6.16M | { |
385 | 6.16M | CHECK(z != NULL); |
386 | 6.16M | |
387 | 6.16M | if (prec == 0) |
388 | 0 | prec = default_precision; |
389 | 6.16M | else if (6.16M prec == 16.16M ) |
390 | 334k | return mp_int_init(z); |
391 | 6.16M | else |
392 | 5.82M | prec = (mp_size) 5.82M ROUND_PREC5.82M (prec); |
393 | 6.16M | |
394 | 5.82M | if (5.82M (5.82M MP_DIGITS5.82M (z) = s_alloc(prec)) == NULL) |
395 | 0 | return MP_MEMORY; |
396 | 5.82M | |
397 | 5.82M | z->digits[0] = 0; |
398 | 5.82M | MP_USED(z) = 1; |
399 | 5.82M | MP_ALLOC(z) = prec; |
400 | 5.82M | MP_SIGN(z) = MP_ZPOS; |
401 | 6.16M | |
402 | 6.16M | return MP_OK; |
403 | 6.16M | } |
404 | | |
405 | | mp_result mp_int_init_copy(mp_int z, mp_int old) |
406 | 7.47M | { |
407 | 7.47M | mp_result res; |
408 | 7.47M | mp_size uold; |
409 | 7.47M | |
410 | 7.47M | CHECK(z != NULL && old != NULL); |
411 | 7.47M | |
412 | 7.47M | uold = MP_USED(old); |
413 | 7.47M | if (uold == 17.47M ) { |
414 | 3.02M | mp_int_init(z); |
415 | 3.02M | } |
416 | 4.45M | else { |
417 | 4.45M | mp_size target = MAX(uold, default_precision); |
418 | 4.45M | |
419 | 4.45M | if ((res = mp_int_init_size(z, target)) != MP_OK) |
420 | 0 | return res; |
421 | 7.47M | } |
422 | 7.47M | |
423 | 7.47M | MP_USED7.47M (z) = uold; |
424 | 7.47M | MP_SIGN7.47M (z) = MP_SIGN7.47M (old); |
425 | 7.47M | COPY(MP_DIGITS(old), MP_DIGITS(z), uold); |
426 | 7.47M | |
427 | 7.47M | return MP_OK; |
428 | 7.47M | } |
429 | | |
430 | | mp_result mp_int_init_value(mp_int z, mp_small value) |
431 | 0 | { |
432 | 0 | mpz_t vtmp; |
433 | 0 | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
434 | 0 |
|
435 | 0 | s_fake(&vtmp, value, vbuf); |
436 | 0 | return mp_int_init_copy(z, &vtmp); |
437 | 0 | } |
438 | | |
439 | | mp_result mp_int_init_uvalue(mp_int z, mp_usmall uvalue) |
440 | 170k | { |
441 | 170k | mpz_t vtmp; |
442 | 170k | mp_digit vbuf[MP_VALUE_DIGITS(uvalue)]; |
443 | 170k | |
444 | 170k | s_ufake(&vtmp, uvalue, vbuf); |
445 | 170k | return mp_int_init_copy(z, &vtmp); |
446 | 170k | } |
447 | | |
448 | | mp_result mp_int_set_value(mp_int z, mp_small value) |
449 | 951k | { |
450 | 951k | mpz_t vtmp; |
451 | 951k | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
452 | 951k | |
453 | 951k | s_fake(&vtmp, value, vbuf); |
454 | 951k | return mp_int_copy(&vtmp, z); |
455 | 951k | } |
456 | | |
457 | | mp_result mp_int_set_uvalue(mp_int z, mp_usmall uvalue) |
458 | 131k | { |
459 | 131k | mpz_t vtmp; |
460 | 131k | mp_digit vbuf[MP_VALUE_DIGITS(uvalue)]; |
461 | 131k | |
462 | 131k | s_ufake(&vtmp, uvalue, vbuf); |
463 | 131k | return mp_int_copy(&vtmp, z); |
464 | 131k | } |
465 | | |
466 | | void mp_int_clear(mp_int z) |
467 | 44.4M | { |
468 | 44.4M | if (z == NULL) |
469 | 0 | return; |
470 | 44.4M | |
471 | 44.4M | if (44.4M MP_DIGITS44.4M (z) != NULL44.4M ) { |
472 | 44.4M | if (MP_DIGITS44.4M (z) != &(z->single)44.4M ) |
473 | 19.7M | s_free(19.7M MP_DIGITS19.7M (z)); |
474 | 44.4M | |
475 | 44.4M | MP_DIGITS(z) = NULL; |
476 | 44.4M | } |
477 | 44.4M | } |
478 | | |
479 | | void mp_int_free(mp_int z) |
480 | 31.6M | { |
481 | 31.6M | NRCHECK(z != NULL); |
482 | 31.6M | |
483 | 31.6M | mp_int_clear(z); |
484 | 31.6M | free(z); /* note: NOT s_free() */ |
485 | 31.6M | } |
486 | | |
487 | | mp_result mp_int_copy(mp_int a, mp_int c) |
488 | 103M | { |
489 | 103M | CHECK(a != NULL && c != NULL); |
490 | 103M | |
491 | 103M | if (a != c103M ) { |
492 | 98.2M | mp_size ua = MP_USED(a); |
493 | 98.2M | mp_digit *da, *dc; |
494 | 98.2M | |
495 | 98.2M | if (!s_pad(c, ua)) |
496 | 0 | return MP_MEMORY; |
497 | 98.2M | |
498 | 98.2M | da = 98.2M MP_DIGITS98.2M (a); dc = MP_DIGITS98.2M (c); |
499 | 98.2M | COPY(da, dc, ua); |
500 | 98.2M | |
501 | 98.2M | MP_USED(c) = ua; |
502 | 98.2M | MP_SIGN98.2M (c) = MP_SIGN98.2M (a); |
503 | 98.2M | } |
504 | 103M | |
505 | 103M | return MP_OK; |
506 | 103M | } |
507 | | |
508 | | void mp_int_swap(mp_int a, mp_int c) |
509 | 0 | { |
510 | 0 | if (a != c0 ) { |
511 | 0 | mpz_t tmp = *a; |
512 | 0 |
|
513 | 0 | *a = *c; |
514 | 0 | *c = tmp; |
515 | 0 |
|
516 | 0 | if (MP_DIGITS0 (a) == &(c->single)0 ) |
517 | 0 | MP_DIGITS0 (a) = &(a->single)0 ; |
518 | 0 | if (MP_DIGITS0 (c) == &(a->single)0 ) |
519 | 0 | MP_DIGITS0 (c) = &(c->single)0 ; |
520 | 0 | } |
521 | 0 | } |
522 | | |
523 | | void mp_int_zero(mp_int z) |
524 | 19.5M | { |
525 | 19.5M | NRCHECK(z != NULL); |
526 | 19.5M | |
527 | 19.5M | z->digits[0] = 0; |
528 | 19.5M | MP_USED(z) = 1; |
529 | 19.5M | MP_SIGN(z) = MP_ZPOS; |
530 | 19.5M | } |
531 | | |
532 | | mp_result mp_int_abs(mp_int a, mp_int c) |
533 | 3.50M | { |
534 | 3.50M | mp_result res; |
535 | 3.50M | |
536 | 3.50M | CHECK(a != NULL && c != NULL); |
537 | 3.50M | |
538 | 3.50M | if ((res = mp_int_copy(a, c)) != MP_OK) |
539 | 0 | return res; |
540 | 3.50M | |
541 | 3.50M | MP_SIGN3.50M (c) = MP_ZPOS; |
542 | 3.50M | return MP_OK; |
543 | 3.50M | } |
544 | | |
545 | | mp_result mp_int_neg(mp_int a, mp_int c) |
546 | 70.9M | { |
547 | 70.9M | mp_result res; |
548 | 70.9M | |
549 | 70.9M | CHECK(a != NULL && c != NULL); |
550 | 70.9M | |
551 | 70.9M | if ((res = mp_int_copy(a, c)) != MP_OK) |
552 | 0 | return res; |
553 | 70.9M | |
554 | 70.9M | if (70.9M CMPZ70.9M (c) != 070.9M ) |
555 | 70.8M | MP_SIGN70.8M (c) = 1 - 70.8M MP_SIGN70.8M (a); |
556 | 70.9M | |
557 | 70.9M | return MP_OK; |
558 | 70.9M | } |
559 | | |
560 | | mp_result mp_int_add(mp_int a, mp_int b, mp_int c) |
561 | 13.2M | { |
562 | 13.2M | mp_size ua, ub, uc, max; |
563 | 13.2M | |
564 | 13.2M | CHECK(a != NULL && b != NULL && c != NULL); |
565 | 13.2M | |
566 | 13.2M | ua = MP_USED13.2M (a); ub = MP_USED13.2M (b); uc = MP_USED13.2M (c); |
567 | 13.2M | max = MAX(ua, ub); |
568 | 13.2M | |
569 | 13.2M | if (MP_SIGN13.2M (a) == 13.2M MP_SIGN13.2M (b)) { |
570 | 7.84M | /* Same sign -- add magnitudes, preserve sign of addends */ |
571 | 7.84M | mp_digit carry; |
572 | 7.84M | |
573 | 7.84M | if (!s_pad(c, max)) |
574 | 0 | return MP_MEMORY; |
575 | 7.84M | |
576 | 7.84M | carry = s_uadd(7.84M MP_DIGITS7.84M (a), MP_DIGITS7.84M (b), MP_DIGITS7.84M (c), ua, ub); |
577 | 7.84M | uc = max; |
578 | 7.84M | |
579 | 7.84M | if (carry7.84M ) { |
580 | 47.9k | if (!s_pad(c, max + 1)) |
581 | 0 | return MP_MEMORY; |
582 | 47.9k | |
583 | 47.9k | c->digits[max] = carry; |
584 | 47.9k | ++uc; |
585 | 47.9k | } |
586 | 7.84M | |
587 | 7.84M | MP_USED7.84M (c) = uc; |
588 | 7.84M | MP_SIGN7.84M (c) = MP_SIGN7.84M (a); |
589 | 7.84M | |
590 | 7.84M | } |
591 | 5.43M | else { |
592 | 5.43M | /* Different signs -- subtract magnitudes, preserve sign of greater */ |
593 | 5.43M | mp_int x, y; |
594 | 5.43M | int cmp = s_ucmp(a, b); /* magnitude comparision, sign ignored */ |
595 | 5.43M | |
596 | 5.43M | /* Set x to max(a, b), y to min(a, b) to simplify later code. |
597 | 5.43M | A special case yields zero for equal magnitudes. |
598 | 5.43M | */ |
599 | 5.43M | if (cmp == 05.43M ) { |
600 | 809k | mp_int_zero(c); |
601 | 809k | return MP_OK; |
602 | 809k | } |
603 | 4.62M | else if (4.62M cmp < 04.62M ) { |
604 | 1.34M | x = b; y = a; |
605 | 1.34M | } |
606 | 3.27M | else { |
607 | 3.27M | x = a; y = b; |
608 | 3.27M | } |
609 | 5.43M | |
610 | 4.62M | if (4.62M !s_pad(c, 4.62M MP_USED4.62M (x))) |
611 | 0 | return MP_MEMORY; |
612 | 4.62M | |
613 | 4.62M | /* Subtract smaller from larger */ |
614 | 4.62M | s_usub(4.62M MP_DIGITS4.62M (x), MP_DIGITS4.62M (y), MP_DIGITS4.62M (c), MP_USED4.62M (x), MP_USED4.62M (y)); |
615 | 4.62M | MP_USED4.62M (c) = MP_USED4.62M (x); |
616 | 4.62M | CLAMP(c); |
617 | 4.62M | |
618 | 4.62M | /* Give result the sign of the larger */ |
619 | 4.62M | MP_SIGN4.62M (c) = MP_SIGN4.62M (x); |
620 | 5.43M | } |
621 | 13.2M | |
622 | 12.4M | return MP_OK; |
623 | 13.2M | } |
624 | | |
625 | | mp_result mp_int_add_value(mp_int a, mp_small value, mp_int c) |
626 | 2.89k | { |
627 | 2.89k | mpz_t vtmp; |
628 | 2.89k | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
629 | 2.89k | |
630 | 2.89k | s_fake(&vtmp, value, vbuf); |
631 | 2.89k | |
632 | 2.89k | return mp_int_add(a, &vtmp, c); |
633 | 2.89k | } |
634 | | |
635 | | mp_result mp_int_sub(mp_int a, mp_int b, mp_int c) |
636 | 79.3M | { |
637 | 79.3M | mp_size ua, ub, uc, max; |
638 | 79.3M | |
639 | 79.3M | CHECK(a != NULL && b != NULL && c != NULL); |
640 | 79.3M | |
641 | 79.3M | ua = MP_USED79.3M (a); ub = MP_USED79.3M (b); uc = MP_USED79.3M (c); |
642 | 79.3M | max = MAX(ua, ub); |
643 | 79.3M | |
644 | 79.3M | if (MP_SIGN79.3M (a) != 79.3M MP_SIGN79.3M (b)) { |
645 | 267k | /* Different signs -- add magnitudes and keep sign of a */ |
646 | 267k | mp_digit carry; |
647 | 267k | |
648 | 267k | if (!s_pad(c, max)) |
649 | 0 | return MP_MEMORY; |
650 | 267k | |
651 | 267k | carry = s_uadd(267k MP_DIGITS267k (a), MP_DIGITS267k (b), MP_DIGITS267k (c), ua, ub); |
652 | 267k | uc = max; |
653 | 267k | |
654 | 267k | if (carry267k ) { |
655 | 2.13k | if (!s_pad(c, max + 1)) |
656 | 0 | return MP_MEMORY; |
657 | 2.13k | |
658 | 2.13k | c->digits[max] = carry; |
659 | 2.13k | ++uc; |
660 | 2.13k | } |
661 | 267k | |
662 | 267k | MP_USED267k (c) = uc; |
663 | 267k | MP_SIGN267k (c) = MP_SIGN267k (a); |
664 | 267k | |
665 | 267k | } |
666 | 79.1M | else { |
667 | 79.1M | /* Same signs -- subtract magnitudes */ |
668 | 79.1M | mp_int x, y; |
669 | 79.1M | mp_sign osign; |
670 | 79.1M | int cmp = s_ucmp(a, b); |
671 | 79.1M | |
672 | 79.1M | if (!s_pad(c, max)) |
673 | 0 | return MP_MEMORY; |
674 | 79.1M | |
675 | 79.1M | if (79.1M cmp >= 079.1M ) { |
676 | 13.4M | x = a; y = b; osign = MP_ZPOS; |
677 | 13.4M | } |
678 | 65.6M | else { |
679 | 65.6M | x = b; y = a; osign = MP_NEG; |
680 | 65.6M | } |
681 | 79.1M | |
682 | 79.1M | if (MP_SIGN79.1M (a) == MP_NEG && 79.1M cmp != 0109k ) |
683 | 91.2k | osign = 1 - osign; |
684 | 79.1M | |
685 | 79.1M | s_usub(MP_DIGITS79.1M (x), MP_DIGITS79.1M (y), MP_DIGITS79.1M (c), MP_USED79.1M (x), MP_USED79.1M (y)); |
686 | 79.1M | MP_USED79.1M (c) = MP_USED79.1M (x); |
687 | 79.1M | CLAMP(c); |
688 | 79.1M | |
689 | 79.1M | MP_SIGN(c) = osign; |
690 | 79.1M | } |
691 | 79.3M | |
692 | 79.3M | return MP_OK; |
693 | 79.3M | } |
694 | | |
695 | | mp_result mp_int_sub_value(mp_int a, mp_small value, mp_int c) |
696 | 15.0k | { |
697 | 15.0k | mpz_t vtmp; |
698 | 15.0k | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
699 | 15.0k | |
700 | 15.0k | s_fake(&vtmp, value, vbuf); |
701 | 15.0k | |
702 | 15.0k | return mp_int_sub(a, &vtmp, c); |
703 | 15.0k | } |
704 | | |
705 | | mp_result mp_int_mul(mp_int a, mp_int b, mp_int c) |
706 | 24.6M | { |
707 | 24.6M | mp_digit *out; |
708 | 24.6M | mp_size osize, ua, ub, p = 0; |
709 | 24.6M | mp_sign osign; |
710 | 24.6M | |
711 | 24.6M | CHECK(a != NULL && b != NULL && c != NULL); |
712 | 24.6M | |
713 | 24.6M | /* If either input is zero, we can shortcut multiplication */ |
714 | 24.6M | if (mp_int_compare_zero(a) == 0 || 24.6M mp_int_compare_zero(b) == 019.8M ) { |
715 | 17.2M | mp_int_zero(c); |
716 | 17.2M | return MP_OK; |
717 | 17.2M | } |
718 | 7.48M | |
719 | 7.48M | /* Output is positive if inputs have same sign, otherwise negative */ |
720 | 7.48M | osign = (7.48M MP_SIGN7.48M (a) == MP_SIGN7.48M (b)) ? MP_ZPOS4.51M : MP_NEG2.97M ; |
721 | 7.48M | |
722 | 7.48M | /* If the output is not identical to any of the inputs, we'll write the |
723 | 7.48M | results directly; otherwise, allocate a temporary space. */ |
724 | 7.48M | ua = MP_USED7.48M (a); ub = MP_USED7.48M (b); |
725 | 7.48M | osize = MAX(ua, ub); |
726 | 7.48M | osize = 4 * ((osize + 1) / 2); |
727 | 7.48M | |
728 | 7.48M | if (c == a || 7.48M c == b5.80M ) { |
729 | 1.68M | p = ROUND_PREC(osize); |
730 | 1.68M | p = MAX(p, default_precision); |
731 | 1.68M | |
732 | 1.68M | if ((out = s_alloc(p)) == NULL) |
733 | 0 | return MP_MEMORY; |
734 | 7.48M | } |
735 | 5.80M | else { |
736 | 5.80M | if (!s_pad(c, osize)) |
737 | 0 | return MP_MEMORY; |
738 | 5.80M | |
739 | 5.80M | out = 5.80M MP_DIGITS5.80M (c); |
740 | 5.80M | } |
741 | 7.48M | ZERO7.48M (out, osize); |
742 | 7.48M | |
743 | 7.48M | if (!s_kmul(7.48M MP_DIGITS7.48M (a), MP_DIGITS7.48M (b), out, ua, ub)) |
744 | 0 | return MP_MEMORY; |
745 | 7.48M | |
746 | 7.48M | /* If we allocated a new buffer, get rid of whatever memory c was already |
747 | 7.48M | using, and fix up its fields to reflect that. |
748 | 7.48M | */ |
749 | 7.48M | if (7.48M out != 7.48M MP_DIGITS7.48M (c)) { |
750 | 1.68M | if ((void *) 1.68M MP_DIGITS1.68M (c) != (void *) c) |
751 | 1.43M | s_free(1.43M MP_DIGITS1.43M (c)); |
752 | 1.68M | MP_DIGITS(c) = out; |
753 | 1.68M | MP_ALLOC(c) = p; |
754 | 1.68M | } |
755 | 7.48M | |
756 | 7.48M | MP_USED(c) = osize; /* might not be true, but we'll fix it ... */ |
757 | 7.48M | CLAMP(c); /* ... right here */ |
758 | 7.48M | MP_SIGN(c) = osign; |
759 | 24.6M | |
760 | 24.6M | return MP_OK; |
761 | 24.6M | } |
762 | | |
763 | | mp_result mp_int_mul_value(mp_int a, mp_small value, mp_int c) |
764 | 0 | { |
765 | 0 | mpz_t vtmp; |
766 | 0 | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
767 | 0 |
|
768 | 0 | s_fake(&vtmp, value, vbuf); |
769 | 0 |
|
770 | 0 | return mp_int_mul(a, &vtmp, c); |
771 | 0 | } |
772 | | |
773 | | mp_result mp_int_mul_pow2(mp_int a, mp_small p2, mp_int c) |
774 | 3.94k | { |
775 | 3.94k | mp_result res; |
776 | 3.94k | CHECK(a != NULL && c != NULL && p2 >= 0); |
777 | 3.94k | |
778 | 3.94k | if ((res = mp_int_copy(a, c)) != MP_OK) |
779 | 0 | return res; |
780 | 3.94k | |
781 | 3.94k | if (3.94k s_qmul(c, (mp_size) p2)3.94k ) |
782 | 3.94k | return MP_OK; |
783 | 3.94k | else |
784 | 0 | return MP_MEMORY; |
785 | 0 | } |
786 | | |
787 | | mp_result mp_int_sqr(mp_int a, mp_int c) |
788 | 0 | { |
789 | 0 | mp_digit *out; |
790 | 0 | mp_size osize, p = 0; |
791 | 0 |
|
792 | 0 | CHECK(a != NULL && c != NULL); |
793 | 0 |
|
794 | 0 | /* Get a temporary buffer big enough to hold the result */ |
795 | 0 | osize = (mp_size) 4 * ((MP_USED(a) + 1) / 2); |
796 | 0 | if (a == c0 ) { |
797 | 0 | p = ROUND_PREC(osize); |
798 | 0 | p = MAX(p, default_precision); |
799 | 0 |
|
800 | 0 | if ((out = s_alloc(p)) == NULL) |
801 | 0 | return MP_MEMORY; |
802 | 0 | } |
803 | 0 | else { |
804 | 0 | if (!s_pad(c, osize)) |
805 | 0 | return MP_MEMORY; |
806 | 0 |
|
807 | 0 | out = 0 MP_DIGITS0 (c); |
808 | 0 | } |
809 | 0 | ZERO0 (out, osize); |
810 | 0 |
|
811 | 0 | s_ksqr(MP_DIGITS0 (a), out, MP_USED0 (a)); |
812 | 0 |
|
813 | 0 | /* Get rid of whatever memory c was already using, and fix up its fields to |
814 | 0 | reflect the new digit array it's using |
815 | 0 | */ |
816 | 0 | if (out != 0 MP_DIGITS0 (c)) { |
817 | 0 | if ((void *) 0 MP_DIGITS0 (c) != (void *) c) |
818 | 0 | s_free(0 MP_DIGITS0 (c)); |
819 | 0 | MP_DIGITS(c) = out; |
820 | 0 | MP_ALLOC(c) = p; |
821 | 0 | } |
822 | 0 |
|
823 | 0 | MP_USED(c) = osize; /* might not be true, but we'll fix it ... */ |
824 | 0 | CLAMP(c); /* ... right here */ |
825 | 0 | MP_SIGN(c) = MP_ZPOS; |
826 | 0 |
|
827 | 0 | return MP_OK; |
828 | 0 | } |
829 | | |
830 | | mp_result mp_int_div(mp_int a, mp_int b, mp_int q, mp_int r) |
831 | 4.20M | { |
832 | 4.20M | int cmp, lg; |
833 | 4.20M | mp_result res = MP_OK; |
834 | 4.20M | mp_int qout, rout; |
835 | 4.20M | mp_sign sa = MP_SIGN4.20M (a), sb = MP_SIGN4.20M (b); |
836 | 4.20M | DECLARE_TEMP(2); |
837 | 4.20M | |
838 | 4.20M | CHECK(a != NULL && b != NULL && q != r); |
839 | 4.20M | |
840 | 4.20M | if (CMPZ4.20M (b) == 04.20M ) |
841 | 0 | return MP_UNDEF; |
842 | 4.20M | else if (4.20M (cmp = s_ucmp(a, b)) < 04.20M ) { |
843 | 790k | /* If |a| < |b|, no division is required: |
844 | 790k | q = 0, r = a |
845 | 790k | */ |
846 | 790k | if (r && 790k (res = mp_int_copy(a, r)) != MP_OK83.2k ) |
847 | 0 | return res; |
848 | 790k | |
849 | 790k | if (790k q790k ) |
850 | 782k | mp_int_zero(q); |
851 | 790k | |
852 | 790k | return MP_OK; |
853 | 790k | } |
854 | 3.41M | else if (3.41M cmp == 03.41M ) { |
855 | 727k | /* If |a| = |b|, no division is required: |
856 | 727k | q = 1 or -1, r = 0 |
857 | 727k | */ |
858 | 727k | if (r) |
859 | 5.23k | mp_int_zero(r); |
860 | 727k | |
861 | 727k | if (q727k ) { |
862 | 723k | mp_int_zero(q); |
863 | 723k | q->digits[0] = 1; |
864 | 723k | |
865 | 723k | if (sa != sb) |
866 | 86.0k | MP_SIGN86.0k (q) = MP_NEG86.0k ; |
867 | 723k | } |
868 | 4.20M | |
869 | 4.20M | return MP_OK; |
870 | 4.20M | } |
871 | 2.68M | |
872 | 2.68M | /* When |a| > |b|, real division is required. We need someplace to store |
873 | 2.68M | quotient and remainder, but q and r are allowed to be NULL or to overlap |
874 | 2.68M | with the inputs. |
875 | 2.68M | */ |
876 | 2.68M | if (2.68M (lg = s_isp2(b)) < 02.68M ) { |
877 | 1.71M | if (q && 1.71M b != q1.70M ) { |
878 | 1.59M | if ((res = mp_int_copy(a, q)) != MP_OK) |
879 | 0 | goto CLEANUP; |
880 | 1.59M | else |
881 | 1.59M | qout = q; |
882 | 1.59M | } |
883 | 117k | else { |
884 | 117k | qout = LAST_TEMP(); |
885 | 117k | SETUP(mp_int_init_copy(LAST_TEMP(), a)); |
886 | 117k | } |
887 | 1.71M | |
888 | 1.71M | if (1.71M r && 1.71M a != r124k ) { |
889 | 124k | if ((res = mp_int_copy(b, r)) != MP_OK) |
890 | 0 | goto CLEANUP; |
891 | 124k | else |
892 | 124k | rout = r; |
893 | 124k | } |
894 | 1.59M | else { |
895 | 1.59M | rout = LAST_TEMP(); |
896 | 1.59M | SETUP(mp_int_init_copy(LAST_TEMP(), b)); |
897 | 1.59M | } |
898 | 1.71M | |
899 | 1.71M | if (1.71M (res = s_udiv_knuth(qout, rout)) != MP_OK1.71M ) goto CLEANUP0 ; |
900 | 2.68M | } |
901 | 966k | else { |
902 | 966k | if (q && 966k (res = mp_int_copy(a, q)) != MP_OK919k ) goto CLEANUP0 ; |
903 | 966k | if (966k r && 966k (res = mp_int_copy(a, r)) != MP_OK138k ) goto CLEANUP0 ; |
904 | 966k | |
905 | 966k | if (966k q966k ) s_qdiv(q, (mp_size) lg)919k ; qout = q; |
906 | 966k | if (r966k ) s_qmod(r, (mp_size) lg)138k ; rout = r; |
907 | 966k | } |
908 | 2.68M | |
909 | 2.68M | /* Recompute signs for output */ |
910 | 2.68M | if (2.68M rout2.68M ) { |
911 | 1.85M | MP_SIGN(rout) = sa; |
912 | 1.85M | if (CMPZ1.85M (rout) == 01.85M ) |
913 | 1.76M | MP_SIGN1.76M (rout) = MP_ZPOS1.76M ; |
914 | 1.85M | } |
915 | 2.68M | if (qout2.68M ) { |
916 | 2.63M | MP_SIGN2.63M (qout) = (sa == sb) ? MP_ZPOS1.86M : MP_NEG769k ; |
917 | 2.63M | if (CMPZ2.63M (qout) == 02.63M ) |
918 | 0 | MP_SIGN0 (qout) = MP_ZPOS0 ; |
919 | 2.63M | } |
920 | 2.68M | |
921 | 2.68M | if (q && 2.68M (res = mp_int_copy(qout, q)) != MP_OK2.62M ) goto CLEANUP0 ; |
922 | 2.68M | if (2.68M r && 2.68M (res = mp_int_copy(rout, r)) != MP_OK263k ) goto CLEANUP0 ; |
923 | 2.68M | |
924 | 2.68M | CLEANUP_TEMP2.68M (); |
925 | 2.68M | return res; |
926 | 4.20M | } |
927 | | |
928 | | mp_result mp_int_mod(mp_int a, mp_int m, mp_int c) |
929 | 0 | { |
930 | 0 | mp_result res; |
931 | 0 | mpz_t tmp; |
932 | 0 | mp_int out; |
933 | 0 |
|
934 | 0 | if (m == c0 ) { |
935 | 0 | mp_int_init(&tmp); |
936 | 0 | out = &tmp; |
937 | 0 | } |
938 | 0 | else { |
939 | 0 | out = c; |
940 | 0 | } |
941 | 0 |
|
942 | 0 | if ((res = mp_int_div(a, m, NULL, out)) != MP_OK) |
943 | 0 | goto CLEANUP; |
944 | 0 |
|
945 | 0 | if (0 CMPZ0 (out) < 00 ) |
946 | 0 | res = mp_int_add(out, m, c); |
947 | 0 | else |
948 | 0 | res = mp_int_copy(out, c); |
949 | 0 |
|
950 | 0 | CLEANUP: |
951 | 0 | if (out != c) |
952 | 0 | mp_int_clear(&tmp); |
953 | 0 |
|
954 | 0 | return res; |
955 | 0 | } |
956 | | |
957 | | mp_result mp_int_div_value(mp_int a, mp_small value, mp_int q, mp_small *r) |
958 | 65.1k | { |
959 | 65.1k | mpz_t vtmp, rtmp; |
960 | 65.1k | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
961 | 65.1k | mp_result res; |
962 | 65.1k | |
963 | 65.1k | mp_int_init(&rtmp); |
964 | 65.1k | s_fake(&vtmp, value, vbuf); |
965 | 65.1k | |
966 | 65.1k | if ((res = mp_int_div(a, &vtmp, q, &rtmp)) != MP_OK) |
967 | 0 | goto CLEANUP; |
968 | 65.1k | |
969 | 65.1k | if (65.1k r65.1k ) |
970 | 54.1k | (void) mp_int_to_int(&rtmp, r); /* can't fail */ |
971 | 65.1k | |
972 | 65.1k | CLEANUP: |
973 | 65.1k | mp_int_clear(&rtmp); |
974 | 65.1k | return res; |
975 | 65.1k | } |
976 | | |
977 | | mp_result mp_int_div_pow2(mp_int a, mp_small p2, mp_int q, mp_int r) |
978 | 0 | { |
979 | 0 | mp_result res = MP_OK; |
980 | 0 |
|
981 | 0 | CHECK(a != NULL && p2 >= 0 && q != r); |
982 | 0 |
|
983 | 0 | if (q != NULL && 0 (res = mp_int_copy(a, q)) == MP_OK0 ) |
984 | 0 | s_qdiv(q, (mp_size) p2); |
985 | 0 |
|
986 | 0 | if (res == MP_OK && 0 r != NULL0 && (res = mp_int_copy(a, r)) == MP_OK0 ) |
987 | 0 | s_qmod(r, (mp_size) p2); |
988 | 0 |
|
989 | 0 | return res; |
990 | 0 | } |
991 | | |
992 | | mp_result mp_int_expt(mp_int a, mp_small b, mp_int c) |
993 | 0 | { |
994 | 0 | mpz_t t; |
995 | 0 | mp_result res; |
996 | 0 | unsigned int v = labs(b); |
997 | 0 |
|
998 | 0 | CHECK(c != NULL); |
999 | 0 | if (b < 0) |
1000 | 0 | return MP_RANGE; |
1001 | 0 |
|
1002 | 0 | if (0 (res = mp_int_init_copy(&t, a)) != MP_OK0 ) |
1003 | 0 | return res; |
1004 | 0 |
|
1005 | 0 | (void) mp_int_set_value(c, 1); |
1006 | 0 | while (v != 00 ) { |
1007 | 0 | if (v & 10 ) { |
1008 | 0 | if ((res = mp_int_mul(c, &t, c)) != MP_OK) |
1009 | 0 | goto CLEANUP; |
1010 | 0 | } |
1011 | 0 |
|
1012 | 0 | v >>= 1; |
1013 | 0 | if (v == 00 ) break0 ; |
1014 | 0 |
|
1015 | 0 | if (0 (res = mp_int_sqr(&t, &t)) != MP_OK0 ) |
1016 | 0 | goto CLEANUP; |
1017 | 0 | } |
1018 | 0 |
|
1019 | 0 | CLEANUP: |
1020 | 0 | mp_int_clear(&t); |
1021 | 0 | return res; |
1022 | 0 | } |
1023 | | |
1024 | | mp_result mp_int_expt_value(mp_small a, mp_small b, mp_int c) |
1025 | 0 | { |
1026 | 0 | mpz_t t; |
1027 | 0 | mp_result res; |
1028 | 0 | unsigned int v = labs(b); |
1029 | 0 |
|
1030 | 0 | CHECK(c != NULL); |
1031 | 0 | if (b < 0) |
1032 | 0 | return MP_RANGE; |
1033 | 0 |
|
1034 | 0 | if (0 (res = mp_int_init_value(&t, a)) != MP_OK0 ) |
1035 | 0 | return res; |
1036 | 0 |
|
1037 | 0 | (void) mp_int_set_value(c, 1); |
1038 | 0 | while (v != 00 ) { |
1039 | 0 | if (v & 10 ) { |
1040 | 0 | if ((res = mp_int_mul(c, &t, c)) != MP_OK) |
1041 | 0 | goto CLEANUP; |
1042 | 0 | } |
1043 | 0 |
|
1044 | 0 | v >>= 1; |
1045 | 0 | if (v == 00 ) break0 ; |
1046 | 0 |
|
1047 | 0 | if (0 (res = mp_int_sqr(&t, &t)) != MP_OK0 ) |
1048 | 0 | goto CLEANUP; |
1049 | 0 | } |
1050 | 0 |
|
1051 | 0 | CLEANUP: |
1052 | 0 | mp_int_clear(&t); |
1053 | 0 | return res; |
1054 | 0 | } |
1055 | | |
1056 | | mp_result mp_int_expt_full(mp_int a, mp_int b, mp_int c) |
1057 | 0 | { |
1058 | 0 | mpz_t t; |
1059 | 0 | mp_result res; |
1060 | 0 | unsigned ix, jx; |
1061 | 0 |
|
1062 | 0 | CHECK(a != NULL && b != NULL && c != NULL); |
1063 | 0 | if (MP_SIGN0 (b) == MP_NEG0 ) |
1064 | 0 | return MP_RANGE; |
1065 | 0 |
|
1066 | 0 | if (0 (res = mp_int_init_copy(&t, a)) != MP_OK0 ) |
1067 | 0 | return res; |
1068 | 0 |
|
1069 | 0 | (void) mp_int_set_value(c, 1); |
1070 | 0 | for (ix = 0; ix < 0 MP_USED0 (b); ++ix0 ) { |
1071 | 0 | mp_digit d = b->digits[ix]; |
1072 | 0 |
|
1073 | 0 | for (jx = 0; jx < 0 MP_DIGIT_BIT0 ; ++jx0 ) { |
1074 | 0 | if (d & 10 ) { |
1075 | 0 | if ((res = mp_int_mul(c, &t, c)) != MP_OK) |
1076 | 0 | goto CLEANUP; |
1077 | 0 | } |
1078 | 0 |
|
1079 | 0 | d >>= 1; |
1080 | 0 | if (d == 0 && 0 ix + 1 == 0 MP_USED0 (b)) |
1081 | 0 | break; |
1082 | 0 | if (0 (res = mp_int_sqr(&t, &t)) != MP_OK0 ) |
1083 | 0 | goto CLEANUP; |
1084 | 0 | } |
1085 | 0 | } |
1086 | 0 |
|
1087 | 0 | CLEANUP: |
1088 | 0 | mp_int_clear(&t); |
1089 | 0 | return res; |
1090 | 0 | } |
1091 | | |
1092 | | int mp_int_compare(mp_int a, mp_int b) |
1093 | 262k | { |
1094 | 262k | mp_sign sa; |
1095 | 262k | |
1096 | 262k | CHECK(a != NULL && b != NULL); |
1097 | 262k | |
1098 | 262k | sa = MP_SIGN(a); |
1099 | 262k | if (sa == 262k MP_SIGN262k (b)) { |
1100 | 229k | int cmp = s_ucmp(a, b); |
1101 | 229k | |
1102 | 229k | /* If they're both zero or positive, the normal comparison applies; if both |
1103 | 229k | negative, the sense is reversed. */ |
1104 | 229k | if (sa == MP_ZPOS) |
1105 | 163k | return cmp; |
1106 | 229k | else |
1107 | 66.4k | return -cmp; |
1108 | 262k | |
1109 | 262k | } |
1110 | 33.3k | else { |
1111 | 33.3k | if (sa == MP_ZPOS) |
1112 | 15.1k | return 1; |
1113 | 33.3k | else |
1114 | 18.1k | return -1; |
1115 | 0 | } |
1116 | 262k | } |
1117 | | |
1118 | | int mp_int_compare_unsigned(mp_int a, mp_int b) |
1119 | 3.10M | { |
1120 | 3.10M | NRCHECK(a != NULL && b != NULL); |
1121 | 3.10M | |
1122 | 3.10M | return s_ucmp(a, b); |
1123 | 3.10M | } |
1124 | | |
1125 | | int mp_int_compare_zero(mp_int z) |
1126 | 62.8M | { |
1127 | 62.8M | NRCHECK(z != NULL); |
1128 | 62.8M | |
1129 | 62.8M | if (MP_USED62.8M (z) == 1 && 62.8M z->digits[0] == 035.4M ) |
1130 | 17.5M | return 0; |
1131 | 45.2M | else if (45.2M MP_SIGN45.2M (z) == MP_ZPOS45.2M ) |
1132 | 32.1M | return 1; |
1133 | 45.2M | else |
1134 | 13.1M | return -1; |
1135 | 0 | } |
1136 | | |
1137 | | int mp_int_compare_value(mp_int z, mp_small value) |
1138 | 78.6M | { |
1139 | 78.6M | mp_sign vsign = (value < 0) ? MP_NEG38.7M : MP_ZPOS39.8M ; |
1140 | 78.6M | int cmp; |
1141 | 78.6M | |
1142 | 78.6M | CHECK(z != NULL); |
1143 | 78.6M | |
1144 | 78.6M | if (vsign == 78.6M MP_SIGN78.6M (z)) { |
1145 | 47.9M | cmp = s_vcmp(z, value); |
1146 | 47.9M | |
1147 | 47.9M | return (vsign == MP_ZPOS) ? cmp39.8M : -cmp8.06M ; |
1148 | 47.9M | } |
1149 | 30.7M | else { |
1150 | 30.7M | return (value < 0) ? 130.7M : -132.7k ; |
1151 | 30.7M | } |
1152 | 0 | } |
1153 | | |
1154 | | int mp_int_compare_uvalue(mp_int z, mp_usmall uv) |
1155 | 0 | { |
1156 | 0 | CHECK(z != NULL); |
1157 | 0 |
|
1158 | 0 | if (MP_SIGN0 (z) == MP_NEG0 ) |
1159 | 0 | return -1; |
1160 | 0 | else |
1161 | 0 | return s_uvcmp(z, uv); |
1162 | 0 | } |
1163 | | |
1164 | | mp_result mp_int_exptmod(mp_int a, mp_int b, mp_int m, mp_int c) |
1165 | 0 | { |
1166 | 0 | mp_result res; |
1167 | 0 | mp_size um; |
1168 | 0 | mp_int s; |
1169 | 0 | DECLARE_TEMP(3); |
1170 | 0 |
|
1171 | 0 | CHECK(a != NULL && b != NULL && c != NULL && m != NULL); |
1172 | 0 |
|
1173 | 0 | /* Zero moduli and negative exponents are not considered. */ |
1174 | 0 | if (CMPZ0 (m) == 00 ) |
1175 | 0 | return MP_UNDEF; |
1176 | 0 | if (0 CMPZ0 (b) < 00 ) |
1177 | 0 | return MP_RANGE; |
1178 | 0 |
|
1179 | 0 | um = 0 MP_USED0 (m); |
1180 | 0 | SETUP(mp_int_init_size(TEMP(0), 2 * um)); |
1181 | 0 | SETUP0 (mp_int_init_size(TEMP(1), 2 * um)); |
1182 | 0 |
|
1183 | 0 | if (0 c == b || 0 c == m0 ) { |
1184 | 0 | SETUP(mp_int_init_size(TEMP(2), 2 * um)); |
1185 | 0 | s = 0 TEMP0 (2); |
1186 | 0 | } |
1187 | 0 | else { |
1188 | 0 | s = c; |
1189 | 0 | } |
1190 | 0 |
|
1191 | 0 | if (0 (res = mp_int_mod(a, m, 0 TEMP0 (0))) != MP_OK) goto CLEANUP0 ; |
1192 | 0 |
|
1193 | 0 | if (0 (res = s_brmu(0 TEMP0 (1), m)) != MP_OK) goto CLEANUP0 ; |
1194 | 0 |
|
1195 | 0 | if (0 (res = s_embar(0 TEMP0 (0), b, m, TEMP0 (1), s)) != MP_OK) |
1196 | 0 | goto CLEANUP; |
1197 | 0 |
|
1198 | 0 | res = mp_int_copy(s, c); |
1199 | 0 |
|
1200 | 0 | CLEANUP_TEMP0 (); |
1201 | 0 | return res; |
1202 | 0 | } |
1203 | | |
1204 | | mp_result mp_int_exptmod_evalue(mp_int a, mp_small value, mp_int m, mp_int c) |
1205 | 0 | { |
1206 | 0 | mpz_t vtmp; |
1207 | 0 | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
1208 | 0 |
|
1209 | 0 | s_fake(&vtmp, value, vbuf); |
1210 | 0 |
|
1211 | 0 | return mp_int_exptmod(a, &vtmp, m, c); |
1212 | 0 | } |
1213 | | |
1214 | | mp_result mp_int_exptmod_bvalue(mp_small value, mp_int b, |
1215 | | mp_int m, mp_int c) |
1216 | 0 | { |
1217 | 0 | mpz_t vtmp; |
1218 | 0 | mp_digit vbuf[MP_VALUE_DIGITS(value)]; |
1219 | 0 |
|
1220 | 0 | s_fake(&vtmp, value, vbuf); |
1221 | 0 |
|
1222 | 0 | return mp_int_exptmod(&vtmp, b, m, c); |
1223 | 0 | } |
1224 | | |
1225 | | mp_result mp_int_exptmod_known(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c) |
1226 | 0 | { |
1227 | 0 | mp_result res; |
1228 | 0 | mp_size um; |
1229 | 0 | mp_int s; |
1230 | 0 | DECLARE_TEMP(2); |
1231 | 0 |
|
1232 | 0 | CHECK(a && b && m && c); |
1233 | 0 |
|
1234 | 0 | /* Zero moduli and negative exponents are not considered. */ |
1235 | 0 | if (CMPZ0 (m) == 00 ) |
1236 | 0 | return MP_UNDEF; |
1237 | 0 | if (0 CMPZ0 (b) < 00 ) |
1238 | 0 | return MP_RANGE; |
1239 | 0 |
|
1240 | 0 | um = 0 MP_USED0 (m); |
1241 | 0 | SETUP(mp_int_init_size(TEMP(0), 2 * um)); |
1242 | 0 |
|
1243 | 0 | if (0 c == b || 0 c == m0 ) { |
1244 | 0 | SETUP(mp_int_init_size(TEMP(1), 2 * um)); |
1245 | 0 | s = 0 TEMP0 (1); |
1246 | 0 | } |
1247 | 0 | else { |
1248 | 0 | s = c; |
1249 | 0 | } |
1250 | 0 |
|
1251 | 0 | if (0 (res = mp_int_mod(a, m, 0 TEMP0 (0))) != MP_OK) goto CLEANUP0 ; |
1252 | 0 |
|
1253 | 0 | if (0 (res = s_embar(0 TEMP0 (0), b, m, mu, s)) != MP_OK) |
1254 | 0 | goto CLEANUP; |
1255 | 0 |
|
1256 | 0 | res = mp_int_copy(s, c); |
1257 | 0 |
|
1258 | 0 | CLEANUP_TEMP0 (); |
1259 | 0 | return res; |
1260 | 0 | } |
1261 | | |
1262 | | mp_result mp_int_redux_const(mp_int m, mp_int c) |
1263 | 0 | { |
1264 | 0 | CHECK(m != NULL && c != NULL && m != c); |
1265 | 0 |
|
1266 | 0 | return s_brmu(c, m); |
1267 | 0 | } |
1268 | | |
1269 | | mp_result mp_int_invmod(mp_int a, mp_int m, mp_int c) |
1270 | 0 | { |
1271 | 0 | mp_result res; |
1272 | 0 | mp_sign sa; |
1273 | 0 | DECLARE_TEMP(2); |
1274 | 0 |
|
1275 | 0 | CHECK(a != NULL && m != NULL && c != NULL); |
1276 | 0 |
|
1277 | 0 | if (CMPZ0 (a) == 0 || 0 CMPZ0 (m) <= 00 ) |
1278 | 0 | return MP_RANGE; |
1279 | 0 |
|
1280 | 0 | sa = 0 MP_SIGN0 (a); /* need this for the result later */ |
1281 | 0 |
|
1282 | 0 | for (last__ = 0; last__ < 20 ; ++last__0 ) |
1283 | 0 | mp_int_init(0 LAST_TEMP0 ()); |
1284 | 0 |
|
1285 | 0 | if ((res = mp_int_egcd(a, m, 0 TEMP0 (0), TEMP0 (1), NULL)) != MP_OK) |
1286 | 0 | goto CLEANUP; |
1287 | 0 |
|
1288 | 0 | if (0 mp_int_compare_value(0 TEMP0 (0), 1) != 0) { |
1289 | 0 | res = MP_UNDEF; |
1290 | 0 | goto CLEANUP; |
1291 | 0 | } |
1292 | 0 |
|
1293 | 0 | /* It is first necessary to constrain the value to the proper range */ |
1294 | 0 | if (0 (res = mp_int_mod(0 TEMP0 (1), m, TEMP0 (1))) != MP_OK) |
1295 | 0 | goto CLEANUP; |
1296 | 0 |
|
1297 | 0 | /* Now, if 'a' was originally negative, the value we have is actually the |
1298 | 0 | magnitude of the negative representative; to get the positive value we |
1299 | 0 | have to subtract from the modulus. Otherwise, the value is okay as it |
1300 | 0 | stands. |
1301 | 0 | */ |
1302 | 0 | if (0 sa == MP_NEG0 ) |
1303 | 0 | res = mp_int_sub(m, 0 TEMP0 (1), c); |
1304 | 0 | else |
1305 | 0 | res = mp_int_copy(0 TEMP0 (1), c); |
1306 | 0 |
|
1307 | 0 | CLEANUP_TEMP0 (); |
1308 | 0 | return res; |
1309 | 0 | } |
1310 | | |
1311 | | /* Binary GCD algorithm due to Josef Stein, 1961 */ |
1312 | | mp_result mp_int_gcd(mp_int a, mp_int b, mp_int c) |
1313 | 2.77M | { |
1314 | 2.77M | int ca, cb, k = 0; |
1315 | 2.77M | mpz_t u, v, t; |
1316 | 2.77M | mp_result res; |
1317 | 2.77M | |
1318 | 2.77M | CHECK(a != NULL && b != NULL && c != NULL); |
1319 | 2.77M | |
1320 | 2.77M | ca = CMPZ(a); |
1321 | 2.77M | cb = CMPZ(b); |
1322 | 2.77M | if (ca == 0 && 2.77M cb == 018 ) |
1323 | 0 | return MP_UNDEF; |
1324 | 2.77M | else if (2.77M ca == 02.77M ) |
1325 | 18 | return mp_int_abs(b, c); |
1326 | 2.77M | else if (2.77M cb == 02.77M ) |
1327 | 9.43k | return mp_int_abs(a, c); |
1328 | 2.76M | |
1329 | 2.76M | mp_int_init(&t); |
1330 | 2.76M | if ((res = mp_int_init_copy(&u, a)) != MP_OK) |
1331 | 0 | goto U; |
1332 | 2.76M | if (2.76M (res = mp_int_init_copy(&v, b)) != MP_OK2.76M ) |
1333 | 0 | goto V; |
1334 | 2.76M | |
1335 | 2.76M | MP_SIGN2.76M (&u) = MP_ZPOS; 2.76M MP_SIGN2.76M (&v) = MP_ZPOS; |
1336 | 2.76M | |
1337 | 2.76M | { /* Divide out common factors of 2 from u and v */ |
1338 | 2.76M | int div2_u = s_dp2k(&u), div2_v = s_dp2k(&v); |
1339 | 2.76M | |
1340 | 2.76M | k = MIN(div2_u, div2_v); |
1341 | 2.76M | s_qdiv(&u, (mp_size) k); |
1342 | 2.76M | s_qdiv(&v, (mp_size) k); |
1343 | 2.76M | } |
1344 | 2.76M | |
1345 | 2.76M | if (mp_int_is_odd2.76M (&u)) { |
1346 | 2.21M | if ((res = mp_int_neg(&v, &t)) != MP_OK) |
1347 | 0 | goto CLEANUP; |
1348 | 2.76M | } |
1349 | 544k | else { |
1350 | 544k | if ((res = mp_int_copy(&u, &t)) != MP_OK) |
1351 | 0 | goto CLEANUP; |
1352 | 2.76M | } |
1353 | 2.76M | |
1354 | 2.76M | for (;;) 2.76M { |
1355 | 78.7M | s_qdiv(&t, s_dp2k(&t)); |
1356 | 78.7M | |
1357 | 78.7M | if (CMPZ78.7M (&t) > 078.7M ) { |
1358 | 10.9M | if ((res = mp_int_copy(&t, &u)) != MP_OK) |
1359 | 0 | goto CLEANUP; |
1360 | 78.7M | } |
1361 | 67.8M | else { |
1362 | 67.8M | if ((res = mp_int_neg(&t, &v)) != MP_OK) |
1363 | 0 | goto CLEANUP; |
1364 | 78.7M | } |
1365 | 78.7M | |
1366 | 78.7M | if (78.7M (res = mp_int_sub(&u, &v, &t)) != MP_OK78.7M ) |
1367 | 0 | goto CLEANUP; |
1368 | 78.7M | |
1369 | 78.7M | if (78.7M CMPZ78.7M (&t) == 078.7M ) |
1370 | 2.76M | break; |
1371 | 2.76M | } |
1372 | 2.76M | |
1373 | 2.76M | if (2.76M (res = mp_int_abs(&u, c)) != MP_OK2.76M ) |
1374 | 0 | goto CLEANUP; |
1375 | 2.76M | if (2.76M !s_qmul(c, (mp_size) k)2.76M ) |
1376 | 0 | res = MP_MEMORY; |
1377 | 2.76M | |
1378 | 2.76M | CLEANUP: |
1379 | 2.76M | mp_int_clear(&v); |
1380 | 2.76M | V: mp_int_clear(&u); |
1381 | 2.76M | U: mp_int_clear(&t); |
1382 | 2.76M | |
1383 | 2.76M | return res; |
1384 | 2.77M | } |
1385 | | |
1386 | | /* This is the binary GCD algorithm again, but this time we keep track of the |
1387 | | elementary matrix operations as we go, so we can get values x and y |
1388 | | satisfying c = ax + by. |
1389 | | */ |
1390 | | mp_result mp_int_egcd(mp_int a, mp_int b, mp_int c, |
1391 | | mp_int x, mp_int y) |
1392 | 0 | { |
1393 | 0 | int k, ca, cb; |
1394 | 0 | mp_result res; |
1395 | 0 | DECLARE_TEMP(8); |
1396 | 0 |
|
1397 | 0 | CHECK(a != NULL && b != NULL && c != NULL && |
1398 | 0 | (x != NULL || y != NULL)); |
1399 | 0 |
|
1400 | 0 | ca = CMPZ(a); |
1401 | 0 | cb = CMPZ(b); |
1402 | 0 | if (ca == 0 && 0 cb == 00 ) |
1403 | 0 | return MP_UNDEF; |
1404 | 0 | else if (0 ca == 00 ) { |
1405 | 0 | if ((res = mp_int_abs(b, c)) != MP_OK0 ) return res0 ; |
1406 | 0 | mp_int_zero(x); (void) mp_int_set_value(y, 1); return MP_OK; |
1407 | 0 | } |
1408 | 0 | else if (0 cb == 00 ) { |
1409 | 0 | if ((res = mp_int_abs(a, c)) != MP_OK0 ) return res0 ; |
1410 | 0 | (void) mp_int_set_value(x, 1); mp_int_zero(y); return MP_OK; |
1411 | 0 | } |
1412 | 0 |
|
1413 | 0 | /* Initialize temporaries: |
1414 | 0 | A:0, B:1, C:2, D:3, u:4, v:5, ou:6, ov:7 */ |
1415 | 0 | for (last__ = 0; 0 last__ < 40 ; ++last__0 ) |
1416 | 0 | mp_int_init(0 LAST_TEMP0 ()); |
1417 | 0 | TEMP(0)->digits[0] = 1; |
1418 | 0 | TEMP(3)->digits[0] = 1; |
1419 | 0 |
|
1420 | 0 | SETUP(mp_int_init_copy(TEMP(4), a)); |
1421 | 0 | SETUP0 (mp_int_init_copy(TEMP(5), b)); |
1422 | 0 |
|
1423 | 0 | /* We will work with absolute values here */ |
1424 | 0 | MP_SIGN0 (TEMP(4)) = MP_ZPOS; |
1425 | 0 | MP_SIGN(TEMP(5)) = MP_ZPOS; |
1426 | 0 |
|
1427 | 0 | { /* Divide out common factors of 2 from u and v */ |
1428 | 0 | int div2_u = s_dp2k(TEMP0 (4)), div2_v = s_dp2k(TEMP0 (5)); |
1429 | 0 |
|
1430 | 0 | k = MIN(div2_u, div2_v); |
1431 | 0 | s_qdiv(TEMP(4), k); |
1432 | 0 | s_qdiv(TEMP(5), k); |
1433 | 0 | } |
1434 | 0 |
|
1435 | 0 | SETUP(mp_int_init_copy(TEMP(6), TEMP(4))); |
1436 | 0 | SETUP0 (mp_int_init_copy(TEMP(7), TEMP(5))); |
1437 | 0 |
|
1438 | 0 | for (;;) 0 { |
1439 | 0 | while (mp_int_is_even0 (TEMP(4))) { |
1440 | 0 | s_qdiv(TEMP(4), 1); |
1441 | 0 |
|
1442 | 0 | if (mp_int_is_odd0 (TEMP(0)) || 0 mp_int_is_odd0 (TEMP(1))) { |
1443 | 0 | if ((res = mp_int_add(0 TEMP0 (0), TEMP0 (7), TEMP0 (0))) != MP_OK) |
1444 | 0 | goto CLEANUP; |
1445 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (1), TEMP0 (6), TEMP0 (1))) != MP_OK) |
1446 | 0 | goto CLEANUP; |
1447 | 0 | } |
1448 | 0 |
|
1449 | 0 | s_qdiv(0 TEMP0 (0), 1); |
1450 | 0 | s_qdiv(TEMP(1), 1); |
1451 | 0 | } |
1452 | 0 |
|
1453 | 0 | while (0 mp_int_is_even0 (TEMP(5))) { |
1454 | 0 | s_qdiv(TEMP(5), 1); |
1455 | 0 |
|
1456 | 0 | if (mp_int_is_odd0 (TEMP(2)) || 0 mp_int_is_odd0 (TEMP(3))) { |
1457 | 0 | if ((res = mp_int_add(0 TEMP0 (2), TEMP0 (7), TEMP0 (2))) != MP_OK) |
1458 | 0 | goto CLEANUP; |
1459 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (3), TEMP0 (6), TEMP0 (3))) != MP_OK) |
1460 | 0 | goto CLEANUP; |
1461 | 0 | } |
1462 | 0 |
|
1463 | 0 | s_qdiv(0 TEMP0 (2), 1); |
1464 | 0 | s_qdiv(TEMP(3), 1); |
1465 | 0 | } |
1466 | 0 |
|
1467 | 0 | if (0 mp_int_compare(0 TEMP0 (4), TEMP0 (5)) >= 0) { |
1468 | 0 | if ((res = mp_int_sub(0 TEMP0 (4), TEMP0 (5), TEMP0 (4))) != MP_OK) goto CLEANUP0 ; |
1469 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (0), TEMP0 (2), TEMP0 (0))) != MP_OK) goto CLEANUP0 ; |
1470 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (1), TEMP0 (3), TEMP0 (1))) != MP_OK) goto CLEANUP0 ; |
1471 | 0 | } |
1472 | 0 | else { |
1473 | 0 | if ((res = mp_int_sub(0 TEMP0 (5), TEMP0 (4), TEMP0 (5))) != MP_OK) goto CLEANUP0 ; |
1474 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (2), TEMP0 (0), TEMP0 (2))) != MP_OK) goto CLEANUP0 ; |
1475 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (3), TEMP0 (1), TEMP0 (3))) != MP_OK) goto CLEANUP0 ; |
1476 | 0 | } |
1477 | 0 |
|
1478 | 0 | if (0 CMPZ0 (TEMP(4)) == 00 ) { |
1479 | 0 | if (x && 0 (res = mp_int_copy(0 TEMP0 (2), x)) != MP_OK) goto CLEANUP0 ; |
1480 | 0 | if (0 y && 0 (res = mp_int_copy(0 TEMP0 (3), y)) != MP_OK) goto CLEANUP0 ; |
1481 | 0 | if (0 c0 ) { |
1482 | 0 | if (!s_qmul(0 TEMP0 (5), k)) { |
1483 | 0 | res = MP_MEMORY; |
1484 | 0 | goto CLEANUP; |
1485 | 0 | } |
1486 | 0 |
|
1487 | 0 | res = mp_int_copy(0 TEMP0 (5), c); |
1488 | 0 | } |
1489 | 0 |
|
1490 | 0 | break; |
1491 | 0 | } |
1492 | 0 | } |
1493 | 0 |
|
1494 | 0 | CLEANUP_TEMP0 (); |
1495 | 0 | return res; |
1496 | 0 | } |
1497 | | |
1498 | | mp_result mp_int_lcm(mp_int a, mp_int b, mp_int c) |
1499 | 271k | { |
1500 | 271k | mpz_t lcm; |
1501 | 271k | mp_result res; |
1502 | 271k | |
1503 | 271k | CHECK(a != NULL && b != NULL && c != NULL); |
1504 | 271k | |
1505 | 271k | /* Since a * b = gcd(a, b) * lcm(a, b), we can compute |
1506 | 271k | lcm(a, b) = (a / gcd(a, b)) * b. |
1507 | 271k | |
1508 | 271k | This formulation insures everything works even if the input |
1509 | 271k | variables share space. |
1510 | 271k | */ |
1511 | 271k | if ((res = mp_int_init(&lcm)) != MP_OK) |
1512 | 0 | return res; |
1513 | 271k | if (271k (res = mp_int_gcd(a, b, &lcm)) != MP_OK271k ) |
1514 | 0 | goto CLEANUP; |
1515 | 271k | if (271k (res = mp_int_div(a, &lcm, &lcm, NULL)) != MP_OK271k ) |
1516 | 0 | goto CLEANUP; |
1517 | 271k | if (271k (res = mp_int_mul(&lcm, b, &lcm)) != MP_OK271k ) |
1518 | 0 | goto CLEANUP; |
1519 | 271k | |
1520 | 271k | res = mp_int_copy(&lcm, c); |
1521 | 271k | |
1522 | 271k | CLEANUP: |
1523 | 271k | mp_int_clear(&lcm); |
1524 | 271k | |
1525 | 271k | return res; |
1526 | 271k | } |
1527 | | |
1528 | | int mp_int_divisible_value(mp_int a, mp_small v) |
1529 | 54.1k | { |
1530 | 54.1k | mp_small rem = 0; |
1531 | 54.1k | |
1532 | 54.1k | if (mp_int_div_value(a, v, NULL, &rem) != MP_OK) |
1533 | 0 | return 0; |
1534 | 54.1k | |
1535 | 54.1k | return rem == 0; |
1536 | 54.1k | } |
1537 | | |
1538 | | int mp_int_is_pow2(mp_int z) |
1539 | 0 | { |
1540 | 0 | CHECK(z != NULL); |
1541 | 0 |
|
1542 | 0 | return s_isp2(z); |
1543 | 0 | } |
1544 | | |
1545 | | /* Implementation of Newton's root finding method, based loosely on a patch |
1546 | | contributed by Hal Finkel <half@halssoftware.com> |
1547 | | modified by M. J. Fromberger. |
1548 | | */ |
1549 | | mp_result mp_int_root(mp_int a, mp_small b, mp_int c) |
1550 | 0 | { |
1551 | 0 | mp_result res = MP_OK; |
1552 | 0 | int flips = 0; |
1553 | 0 | DECLARE_TEMP(5); |
1554 | 0 |
|
1555 | 0 | CHECK(a != NULL && c != NULL && b > 0); |
1556 | 0 |
|
1557 | 0 | if (b == 10 ) { |
1558 | 0 | return mp_int_copy(a, c); |
1559 | 0 | } |
1560 | 0 | if (0 MP_SIGN0 (a) == MP_NEG0 ) { |
1561 | 0 | if (b % 2 == 0) |
1562 | 0 | return MP_UNDEF; /* root does not exist for negative a with even b */ |
1563 | 0 | else |
1564 | 0 | flips = 1; |
1565 | 0 | } |
1566 | 0 |
|
1567 | 0 | SETUP0 (mp_int_init_copy(LAST_TEMP(), a)); |
1568 | 0 | SETUP0 (mp_int_init_copy(LAST_TEMP(), a)); |
1569 | 0 | SETUP0 (mp_int_init(LAST_TEMP())); |
1570 | 0 | SETUP0 (mp_int_init(LAST_TEMP())); |
1571 | 0 | SETUP0 (mp_int_init(LAST_TEMP())); |
1572 | 0 |
|
1573 | 0 | (void) mp_int_abs(0 TEMP0 (0), TEMP0 (0)); |
1574 | 0 | (void) mp_int_abs(TEMP0 (1), TEMP0 (1)); |
1575 | 0 |
|
1576 | 0 | for (;;) { |
1577 | 0 | if ((res = mp_int_expt(0 TEMP0 (1), b, TEMP0 (2))) != MP_OK) |
1578 | 0 | goto CLEANUP; |
1579 | 0 |
|
1580 | 0 | if (0 mp_int_compare_unsigned(0 TEMP0 (2), TEMP0 (0)) <= 0) |
1581 | 0 | break; |
1582 | 0 |
|
1583 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (2), TEMP0 (0), TEMP0 (2))) != MP_OK) |
1584 | 0 | goto CLEANUP; |
1585 | 0 | if (0 (res = mp_int_expt(0 TEMP0 (1), b - 1, TEMP0 (3))) != MP_OK) |
1586 | 0 | goto CLEANUP; |
1587 | 0 | if (0 (res = mp_int_mul_value(0 TEMP0 (3), b, TEMP0 (3))) != MP_OK) |
1588 | 0 | goto CLEANUP; |
1589 | 0 | if (0 (res = mp_int_div(0 TEMP0 (2), TEMP0 (3), TEMP0 (4), NULL)) != MP_OK) |
1590 | 0 | goto CLEANUP; |
1591 | 0 | if (0 (res = mp_int_sub(0 TEMP0 (1), TEMP0 (4), TEMP0 (4))) != MP_OK) |
1592 | 0 | goto CLEANUP; |
1593 | 0 |
|
1594 | 0 | if (0 mp_int_compare_unsigned(0 TEMP0 (1), TEMP0 (4)) == 0) { |
1595 | 0 | if ((res = mp_int_sub_value(0 TEMP0 (4), 1, TEMP0 (4))) != MP_OK) |
1596 | 0 | goto CLEANUP; |
1597 | 0 | } |
1598 | 0 | if (0 (res = mp_int_copy(0 TEMP0 (4), TEMP0 (1))) != MP_OK) |
1599 | 0 | goto CLEANUP; |
1600 | 0 | } |
1601 | 0 |
|
1602 | 0 | if (0 (res = mp_int_copy(0 TEMP0 (1), c)) != MP_OK) |
1603 | 0 | goto CLEANUP; |
1604 | 0 |
|
1605 | 0 | /* If the original value of a was negative, flip the output sign. */ |
1606 | 0 | if (0 flips0 ) |
1607 | 0 | (void) mp_int_neg(c, c); /* cannot fail */ |
1608 | 0 |
|
1609 | 0 | CLEANUP_TEMP0 (); |
1610 | 0 | return res; |
1611 | 0 | } |
1612 | | |
1613 | | mp_result mp_int_to_int(mp_int z, mp_small *out) |
1614 | 44.6M | { |
1615 | 44.6M | mp_usmall uv = 0; |
1616 | 44.6M | mp_size uz; |
1617 | 44.6M | mp_digit *dz; |
1618 | 44.6M | mp_sign sz; |
1619 | 44.6M | |
1620 | 44.6M | CHECK(z != NULL); |
1621 | 44.6M | |
1622 | 44.6M | /* Make sure the value is representable as a small integer */ |
1623 | 44.6M | sz = MP_SIGN(z); |
1624 | 44.6M | if ((sz == MP_ZPOS && 44.6M mp_int_compare_value(z, 36.5M MP_SMALL_MAX36.5M ) > 0) || |
1625 | 38.7M | mp_int_compare_value(z, 38.7M MP_SMALL_MIN38.7M ) < 0) |
1626 | 8.32M | return MP_RANGE; |
1627 | 36.2M | |
1628 | 36.2M | uz = 36.2M MP_USED36.2M (z); |
1629 | 36.2M | dz = MP_DIGITS(z) + uz - 1; |
1630 | 36.2M | |
1631 | 80.7M | while (uz > 080.7M ) { |
1632 | 44.4M | uv <<= MP_DIGIT_BIT/2; |
1633 | 44.4M | uv = (uv << (MP_DIGIT_BIT/2)) | *dz--; |
1634 | 44.4M | --uz; |
1635 | 44.4M | } |
1636 | 36.2M | |
1637 | 36.2M | if (out) |
1638 | 36.2M | *out = (mp_small)((sz == MP_NEG) ? 36.2M -uv5.62M : uv30.6M ); |
1639 | 44.6M | |
1640 | 44.6M | return MP_OK; |
1641 | 44.6M | } |
1642 | | |
1643 | | mp_result mp_int_to_uint(mp_int z, mp_usmall *out) |
1644 | 0 | { |
1645 | 0 | mp_usmall uv = 0; |
1646 | 0 | mp_size uz; |
1647 | 0 | mp_digit *dz; |
1648 | 0 | mp_sign sz; |
1649 | 0 | |
1650 | 0 | CHECK(z != NULL); |
1651 | 0 |
|
1652 | 0 | /* Make sure the value is representable as an unsigned small integer */ |
1653 | 0 | sz = MP_SIGN(z); |
1654 | 0 | if (sz == MP_NEG || 0 mp_int_compare_uvalue(z, 0 MP_USMALL_MAX0 ) > 0) |
1655 | 0 | return MP_RANGE; |
1656 | 0 | |
1657 | 0 | uz = 0 MP_USED0 (z); |
1658 | 0 | dz = MP_DIGITS(z) + uz - 1; |
1659 | 0 | |
1660 | 0 | while (uz > 00 ) { |
1661 | 0 | uv <<= MP_DIGIT_BIT/2; |
1662 | 0 | uv = (uv << (MP_DIGIT_BIT/2)) | *dz--; |
1663 | 0 | --uz; |
1664 | 0 | } |
1665 | 0 | |
1666 | 0 | if (out) |
1667 | 0 | *out = uv; |
1668 | 0 | |
1669 | 0 | return MP_OK; |
1670 | 0 | } |
1671 | | |
1672 | | mp_result mp_int_to_string(mp_int z, mp_size radix, |
1673 | | char *str, int limit) |
1674 | 1.43k | { |
1675 | 1.43k | mp_result res; |
1676 | 1.43k | int cmp = 0; |
1677 | 1.43k | |
1678 | 1.43k | CHECK(z != NULL && str != NULL && limit >= 2); |
1679 | 1.43k | |
1680 | 1.43k | if (radix < 1.43k MP_MIN_RADIX1.43k || radix > 1.43k MP_MAX_RADIX1.43k ) |
1681 | 0 | return MP_RANGE; |
1682 | 1.43k | |
1683 | 1.43k | if (1.43k CMPZ1.43k (z) == 01.43k ) { |
1684 | 0 | *str++ = s_val2ch(0, 1); |
1685 | 0 | } |
1686 | 1.43k | else { |
1687 | 1.43k | mpz_t tmp; |
1688 | 1.43k | char *h, *t; |
1689 | 1.43k | |
1690 | 1.43k | if ((res = mp_int_init_copy(&tmp, z)) != MP_OK) |
1691 | 0 | return res; |
1692 | 1.43k | |
1693 | 1.43k | if (1.43k MP_SIGN1.43k (z) == MP_NEG1.43k ) { |
1694 | 741 | *str++ = '-'; |
1695 | 741 | --limit; |
1696 | 741 | } |
1697 | 1.43k | h = str; |
1698 | 1.43k | |
1699 | 1.43k | /* Generate digits in reverse order until finished or limit reached */ |
1700 | 27.5k | for (/* */; limit > 027.5k ; --limit26.0k ) { |
1701 | 27.5k | mp_digit d; |
1702 | 27.5k | |
1703 | 27.5k | if ((cmp = 27.5k CMPZ27.5k (&tmp)) == 0) |
1704 | 1.43k | break; |
1705 | 26.0k | |
1706 | 26.0k | d = s_ddiv(&tmp, (mp_digit)radix); |
1707 | 26.0k | *str++ = s_val2ch(d, 1); |
1708 | 26.0k | } |
1709 | 1.43k | t = str - 1; |
1710 | 1.43k | |
1711 | 1.43k | /* Put digits back in correct output order */ |
1712 | 13.9k | while (h < t13.9k ) { |
1713 | 12.4k | char tc = *h; |
1714 | 12.4k | *h++ = *t; |
1715 | 12.4k | *t-- = tc; |
1716 | 12.4k | } |
1717 | 1.43k | |
1718 | 1.43k | mp_int_clear(&tmp); |
1719 | 1.43k | } |
1720 | 1.43k | |
1721 | 1.43k | *str = '\0'; |
1722 | 1.43k | if (cmp == 0) |
1723 | 1.43k | return MP_OK; |
1724 | 1.43k | else |
1725 | 0 | return MP_TRUNC; |
1726 | 0 | } |
1727 | | |
1728 | | mp_result mp_int_string_len(mp_int z, mp_size radix) |
1729 | 1.54k | { |
1730 | 1.54k | int len; |
1731 | 1.54k | |
1732 | 1.54k | CHECK(z != NULL); |
1733 | 1.54k | |
1734 | 1.54k | if (radix < 1.54k MP_MIN_RADIX1.54k || radix > 1.54k MP_MAX_RADIX1.54k ) |
1735 | 0 | return MP_RANGE; |
1736 | 1.54k | |
1737 | 1.54k | len = s_outlen(z, radix) + 1; /* for terminator */ |
1738 | 1.54k | |
1739 | 1.54k | /* Allow for sign marker on negatives */ |
1740 | 1.54k | if (MP_SIGN1.54k (z) == MP_NEG1.54k ) |
1741 | 771 | len += 1; |
1742 | 1.54k | |
1743 | 1.54k | return len; |
1744 | 1.54k | } |
1745 | | |
1746 | | /* Read zero-terminated string into z */ |
1747 | | mp_result mp_int_read_string(mp_int z, mp_size radix, const char *str) |
1748 | 150 | { |
1749 | 150 | return mp_int_read_cstring(z, radix, str, NULL); |
1750 | 150 | } |
1751 | | |
1752 | | mp_result mp_int_read_cstring(mp_int z, mp_size radix, const char *str, char **end) |
1753 | 150 | { |
1754 | 150 | int ch; |
1755 | 150 | |
1756 | 150 | CHECK(z != NULL && str != NULL); |
1757 | 150 | |
1758 | 150 | if (radix < 150 MP_MIN_RADIX150 || radix > 150 MP_MAX_RADIX150 ) |
1759 | 0 | return MP_RANGE; |
1760 | 150 | |
1761 | 150 | /* Skip leading whitespace */ |
1762 | 150 | while (150 isspace((int)*str)150 ) |
1763 | 0 | ++str; |
1764 | 150 | |
1765 | 150 | /* Handle leading sign tag (+/-, positive default) */ |
1766 | 150 | switch (*str) { |
1767 | 55 | case '-': |
1768 | 55 | MP_SIGN(z) = MP_NEG; |
1769 | 55 | ++str; |
1770 | 55 | break; |
1771 | 0 | case '+': |
1772 | 0 | ++str; /* fallthrough */ |
1773 | 95 | default: |
1774 | 95 | MP_SIGN(z) = MP_ZPOS; |
1775 | 0 | break; |
1776 | 150 | } |
1777 | 150 | |
1778 | 150 | /* Skip leading zeroes */ |
1779 | 150 | while (150 (ch = s_ch2val(*str, radix)) == 0150 ) |
1780 | 0 | ++str; |
1781 | 150 | |
1782 | 150 | /* Make sure there is enough space for the value */ |
1783 | 150 | if (!s_pad(z, s_inlen(strlen(str), radix))) |
1784 | 0 | return MP_MEMORY; |
1785 | 150 | |
1786 | 150 | MP_USED150 (z) = 1; z->digits[0] = 0; |
1787 | 150 | |
1788 | 1.83k | while (*str != '\0' && 1.83k ((ch = s_ch2val(*str, radix)) >= 0)1.68k ) { |
1789 | 1.68k | s_dmul(z, (mp_digit)radix); |
1790 | 1.68k | s_dadd(z, (mp_digit)ch); |
1791 | 1.68k | ++str; |
1792 | 1.68k | } |
1793 | 150 | |
1794 | 150 | CLAMP(z); |
1795 | 150 | |
1796 | 150 | /* Override sign for zero, even if negative specified. */ |
1797 | 150 | if (CMPZ150 (z) == 0150 ) |
1798 | 0 | MP_SIGN0 (z) = MP_ZPOS0 ; |
1799 | 150 | |
1800 | 150 | if (end != NULL) |
1801 | 0 | *end = (char *)str; |
1802 | 150 | |
1803 | 150 | /* Return a truncation error if the string has unprocessed characters |
1804 | 150 | remaining, so the caller can tell if the whole string was done */ |
1805 | 150 | if (*str != '\0') |
1806 | 0 | return MP_TRUNC; |
1807 | 150 | else |
1808 | 150 | return MP_OK; |
1809 | 0 | } |
1810 | | |
1811 | | mp_result mp_int_count_bits(mp_int z) |
1812 | 3.97k | { |
1813 | 3.97k | mp_size nbits = 0, uz; |
1814 | 3.97k | mp_digit d; |
1815 | 3.97k | |
1816 | 3.97k | CHECK(z != NULL); |
1817 | 3.97k | |
1818 | 3.97k | uz = MP_USED(z); |
1819 | 3.97k | if (uz == 1 && 3.97k z->digits[0] == 02.63k ) |
1820 | 0 | return 1; |
1821 | 3.97k | |
1822 | 3.97k | --uz; |
1823 | 3.97k | nbits = uz * MP_DIGIT_BIT; |
1824 | 3.97k | d = z->digits[uz]; |
1825 | 3.97k | |
1826 | 62.7k | while (d != 062.7k ) { |
1827 | 58.7k | d >>= 1; |
1828 | 58.7k | ++nbits; |
1829 | 58.7k | } |
1830 | 3.97k | |
1831 | 3.97k | return nbits; |
1832 | 3.97k | } |
1833 | | |
1834 | | mp_result mp_int_to_binary(mp_int z, unsigned char *buf, int limit) |
1835 | 0 | { |
1836 | 0 | static const int PAD_FOR_2C = 1; |
1837 | 0 |
|
1838 | 0 | mp_result res; |
1839 | 0 | int limpos = limit; |
1840 | 0 |
|
1841 | 0 | CHECK(z != NULL && buf != NULL); |
1842 | 0 |
|
1843 | 0 | res = s_tobin(z, buf, &limpos, PAD_FOR_2C); |
1844 | 0 |
|
1845 | 0 | if (MP_SIGN0 (z) == MP_NEG0 ) |
1846 | 0 | s_2comp(buf, limpos); |
1847 | 0 |
|
1848 | 0 | return res; |
1849 | 0 | } |
1850 | | |
1851 | | mp_result mp_int_read_binary(mp_int z, unsigned char *buf, int len) |
1852 | 0 | { |
1853 | 0 | mp_size need, i; |
1854 | 0 | unsigned char *tmp; |
1855 | 0 | mp_digit *dz; |
1856 | 0 |
|
1857 | 0 | CHECK(z != NULL && buf != NULL && len > 0); |
1858 | 0 |
|
1859 | 0 | /* Figure out how many digits are needed to represent this value */ |
1860 | 0 | need = ((len * CHAR_BIT) + (MP_DIGIT_BIT0 - 1)) / MP_DIGIT_BIT0 ; |
1861 | 0 | if (!s_pad(z, need)) |
1862 | 0 | return MP_MEMORY; |
1863 | 0 |
|
1864 | 0 | mp_int_zero(z); |
1865 | 0 |
|
1866 | 0 | /* If the high-order bit is set, take the 2's complement before reading the |
1867 | 0 | value (it will be restored afterward) */ |
1868 | 0 | if (buf[0] >> (CHAR_BIT - 1)0 ) { |
1869 | 0 | MP_SIGN(z) = MP_NEG; |
1870 | 0 | s_2comp(buf, len); |
1871 | 0 | } |
1872 | 0 |
|
1873 | 0 | dz = MP_DIGITS(z); |
1874 | 0 | for (tmp = buf, i = len; i > 00 ; --i, ++tmp0 ) { |
1875 | 0 | s_qmul(z, (mp_size) CHAR_BIT); |
1876 | 0 | *dz |= *tmp; |
1877 | 0 | } |
1878 | 0 |
|
1879 | 0 | /* Restore 2's complement if we took it before */ |
1880 | 0 | if (MP_SIGN0 (z) == MP_NEG0 ) |
1881 | 0 | s_2comp(buf, len); |
1882 | 0 |
|
1883 | 0 | return MP_OK; |
1884 | 0 | } |
1885 | | |
1886 | | mp_result mp_int_binary_len(mp_int z) |
1887 | 0 | { |
1888 | 0 | mp_result res = mp_int_count_bits(z); |
1889 | 0 | int bytes = mp_int_unsigned_len(z); |
1890 | 0 |
|
1891 | 0 | if (res <= 0) |
1892 | 0 | return res; |
1893 | 0 |
|
1894 | 0 | bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT; |
1895 | 0 |
|
1896 | 0 | /* If the highest-order bit falls exactly on a byte boundary, we need to pad |
1897 | 0 | with an extra byte so that the sign will be read correctly when reading it |
1898 | 0 | back in. */ |
1899 | 0 | if (bytes * CHAR_BIT == res) |
1900 | 0 | ++bytes; |
1901 | 0 |
|
1902 | 0 | return bytes; |
1903 | 0 | } |
1904 | | |
1905 | | mp_result mp_int_to_unsigned(mp_int z, unsigned char *buf, int limit) |
1906 | 0 | { |
1907 | 0 | static const int NO_PADDING = 0; |
1908 | 0 |
|
1909 | 0 | CHECK(z != NULL && buf != NULL); |
1910 | 0 |
|
1911 | 0 | return s_tobin(z, buf, &limit, NO_PADDING); |
1912 | 0 | } |
1913 | | |
1914 | | mp_result mp_int_read_unsigned(mp_int z, unsigned char *buf, int len) |
1915 | 0 | { |
1916 | 0 | mp_size need, i; |
1917 | 0 | unsigned char *tmp; |
1918 | 0 |
|
1919 | 0 | CHECK(z != NULL && buf != NULL && len > 0); |
1920 | 0 |
|
1921 | 0 | /* Figure out how many digits are needed to represent this value */ |
1922 | 0 | need = ((len * CHAR_BIT) + (MP_DIGIT_BIT0 - 1)) / MP_DIGIT_BIT0 ; |
1923 | 0 | if (!s_pad(z, need)) |
1924 | 0 | return MP_MEMORY; |
1925 | 0 |
|
1926 | 0 | mp_int_zero(z); |
1927 | 0 |
|
1928 | 0 | for (tmp = buf, i = len; i > 00 ; --i, ++tmp0 ) { |
1929 | 0 | (void) s_qmul(z, CHAR_BIT); |
1930 | 0 | *MP_DIGITS(z) |= *tmp; |
1931 | 0 | } |
1932 | 0 |
|
1933 | 0 | return MP_OK; |
1934 | 0 | } |
1935 | | |
1936 | | mp_result mp_int_unsigned_len(mp_int z) |
1937 | 2.43k | { |
1938 | 2.43k | mp_result res = mp_int_count_bits(z); |
1939 | 2.43k | int bytes; |
1940 | 2.43k | |
1941 | 2.43k | if (res <= 0) |
1942 | 0 | return res; |
1943 | 2.43k | |
1944 | 2.43k | bytes = (res + (CHAR_BIT - 1)) / CHAR_BIT; |
1945 | 2.43k | |
1946 | 2.43k | return bytes; |
1947 | 2.43k | } |
1948 | | |
1949 | | const char *mp_error_string(mp_result res) |
1950 | 0 | { |
1951 | 0 | int ix; |
1952 | 0 | if (res > 0) |
1953 | 0 | return s_unknown_err; |
1954 | 0 |
|
1955 | 0 | res = -res; |
1956 | 0 | for (ix = 0; ix < res && 0 s_error_msg[ix] != NULL0 ; ++ix0 ) |
1957 | 0 | ; |
1958 | 0 |
|
1959 | 0 | if (s_error_msg[ix] != NULL) |
1960 | 0 | return s_error_msg[ix]; |
1961 | 0 | else |
1962 | 0 | return s_unknown_err; |
1963 | 0 | } |
1964 | | |
1965 | | /*------------------------------------------------------------------------*/ |
1966 | | /* Private functions for internal use. These make assumptions. */ |
1967 | | |
1968 | | STATIC mp_digit *s_alloc(mp_size num) |
1969 | 21.1M | { |
1970 | 21.1M | mp_digit *out = malloc(num * sizeof(mp_digit)); |
1971 | 21.1M | |
1972 | 21.1M | assert(out != NULL); /* for debugging */ |
1973 | | #if DEBUG > 1 |
1974 | | { |
1975 | | mp_digit v = (mp_digit) 0xdeadbeef; |
1976 | | int ix; |
1977 | | |
1978 | | for (ix = 0; ix < num; ++ix) |
1979 | | out[ix] = v; |
1980 | | } |
1981 | | #endif |
1982 | | |
1983 | 21.1M | return out; |
1984 | 21.1M | } |
1985 | | |
1986 | | STATIC mp_digit *s_realloc(mp_digit *old, mp_size osize, mp_size nsize) |
1987 | 645k | { |
1988 | | #if DEBUG > 1 |
1989 | | mp_digit *new = s_alloc(nsize); |
1990 | | int ix; |
1991 | | |
1992 | | for (ix = 0; ix < nsize; ++ix) |
1993 | | new[ix] = (mp_digit) 0xdeadbeef; |
1994 | | |
1995 | | memcpy(new, old, osize * sizeof(mp_digit)); |
1996 | | #else |
1997 | | mp_digit *new = realloc(old, nsize * sizeof(mp_digit)); |
1998 | 645k | |
1999 | 645k | assert(new != NULL); /* for debugging */ |
2000 | 645k | #endif |
2001 | 645k | return new; |
2002 | 645k | } |
2003 | | |
2004 | | STATIC void s_free(void *ptr) |
2005 | 21.1M | { |
2006 | 21.1M | free(ptr); |
2007 | 21.1M | } |
2008 | | |
2009 | | STATIC int s_pad(mp_int z, mp_size min) |
2010 | 200M | { |
2011 | 200M | if (MP_ALLOC200M (z) < min200M ) { |
2012 | 14.2M | mp_size nsize = ROUND_PREC(min); |
2013 | 14.2M | mp_digit *tmp; |
2014 | 14.2M | |
2015 | 14.2M | if ((void *)z->digits == (void *)z14.2M ) { |
2016 | 13.6M | if ((tmp = s_alloc(nsize)) == NULL) |
2017 | 0 | return 0; |
2018 | 13.6M | |
2019 | 13.6M | COPY13.6M (MP_DIGITS(z), tmp, MP_USED(z)); |
2020 | 13.6M | } |
2021 | 645k | else if (645k (tmp = s_realloc(645k MP_DIGITS645k (z), MP_ALLOC645k (z), nsize)) == NULL) |
2022 | 0 | return 0; |
2023 | 14.2M | |
2024 | 14.2M | MP_DIGITS14.2M (z) = tmp; |
2025 | 14.2M | MP_ALLOC(z) = nsize; |
2026 | 14.2M | } |
2027 | 200M | |
2028 | 200M | return 1; |
2029 | 200M | } |
2030 | | |
2031 | | /* Note: This will not work correctly when value == MP_SMALL_MIN */ |
2032 | | STATIC void s_fake(mp_int z, mp_small value, mp_digit vbuf[]) |
2033 | 1.03M | { |
2034 | 1.03M | mp_usmall uv = (mp_usmall) (value < 0) ? -value6.44k : value1.02M ; |
2035 | 1.03M | s_ufake(z, uv, vbuf); |
2036 | 1.03M | if (value < 0) |
2037 | 6.44k | z->sign = MP_NEG; |
2038 | 1.03M | } |
2039 | | |
2040 | | STATIC void s_ufake(mp_int z, mp_usmall value, mp_digit vbuf[]) |
2041 | 49.2M | { |
2042 | 49.2M | mp_size ndig = (mp_size) s_uvpack(value, vbuf); |
2043 | 49.2M | |
2044 | 49.2M | z->used = ndig; |
2045 | 49.2M | z->alloc = MP_VALUE_DIGITS(value); |
2046 | 49.2M | z->sign = MP_ZPOS; |
2047 | 49.2M | z->digits = vbuf; |
2048 | 49.2M | } |
2049 | | |
2050 | | STATIC int s_cdig(mp_digit *da, mp_digit *db, mp_size len) |
2051 | 57.5M | { |
2052 | 57.5M | mp_digit *dat = da + len - 1, *dbt = db + len - 1; |
2053 | 57.5M | |
2054 | 76.9M | for (/* */; len != 076.9M ; --len, --dat, --dbt19.3M ) { |
2055 | 68.7M | if (*dat > *dbt) |
2056 | 13.3M | return 1; |
2057 | 55.4M | else if (55.4M *dat < *dbt55.4M ) |
2058 | 36.0M | return -1; |
2059 | 68.7M | } |
2060 | 57.5M | |
2061 | 8.16M | return 0; |
2062 | 57.5M | } |
2063 | | |
2064 | | STATIC int s_uvpack(mp_usmall uv, mp_digit t[]) |
2065 | 49.2M | { |
2066 | 49.2M | int ndig = 0; |
2067 | 49.2M | |
2068 | 49.2M | if (uv == 0) |
2069 | 856k | t[ndig++] = 0; |
2070 | 48.4M | else { |
2071 | 141M | while (uv != 0141M ) { |
2072 | 93.0M | t[ndig++] = (mp_digit) uv; |
2073 | 93.0M | uv >>= MP_DIGIT_BIT/2; |
2074 | 93.0M | uv >>= MP_DIGIT_BIT/2; |
2075 | 93.0M | } |
2076 | 48.4M | } |
2077 | 49.2M | |
2078 | 49.2M | return ndig; |
2079 | 49.2M | } |
2080 | | |
2081 | | STATIC int s_ucmp(mp_int a, mp_int b) |
2082 | 142M | { |
2083 | 142M | mp_size ua = MP_USED142M (a), ub = MP_USED142M (b); |
2084 | 142M | |
2085 | 142M | if (ua > ub) |
2086 | 17.1M | return 1; |
2087 | 124M | else if (124M ub > ua124M ) |
2088 | 67.3M | return -1; |
2089 | 124M | else |
2090 | 124M | return s_cdig(57.5M MP_DIGITS57.5M (a), MP_DIGITS57.5M (b), ua); |
2091 | 0 | } |
2092 | | |
2093 | | STATIC int s_vcmp(mp_int a, mp_small v) |
2094 | 47.9M | { |
2095 | 47.9M | mp_usmall uv = (v < 0) ? -(mp_usmall) v8.06M : (mp_usmall) v39.8M ; |
2096 | 47.9M | return s_uvcmp(a, uv); |
2097 | 47.9M | } |
2098 | | |
2099 | | STATIC int s_uvcmp(mp_int a, mp_usmall uv) |
2100 | 47.9M | { |
2101 | 47.9M | mpz_t vtmp; |
2102 | 47.9M | mp_digit vdig[MP_VALUE_DIGITS(uv)]; |
2103 | 47.9M | |
2104 | 47.9M | s_ufake(&vtmp, uv, vdig); |
2105 | 47.9M | return s_ucmp(a, &vtmp); |
2106 | 47.9M | } |
2107 | | |
2108 | | STATIC mp_digit s_uadd(mp_digit *da, mp_digit *db, mp_digit *dc, |
2109 | | mp_size size_a, mp_size size_b) |
2110 | 8.11M | { |
2111 | 8.11M | mp_size pos; |
2112 | 8.11M | mp_word w = 0; |
2113 | 8.11M | |
2114 | 8.11M | /* Insure that da is the longer of the two to simplify later code */ |
2115 | 8.11M | if (size_b > size_a8.11M ) { |
2116 | 823k | SWAP(mp_digit *, da, db); |
2117 | 823k | SWAP(mp_size, size_a, size_b); |
2118 | 823k | } |
2119 | 8.11M | |
2120 | 8.11M | /* Add corresponding digits until the shorter number runs out */ |
2121 | 17.5M | for (pos = 0; pos < size_b17.5M ; ++pos, ++da, ++db, ++dc9.47M ) { |
2122 | 9.47M | w = w + (mp_word) *da + (mp_word) *db; |
2123 | 9.47M | *dc = LOWER_HALF(w); |
2124 | 9.47M | w = UPPER_HALF(w); |
2125 | 9.47M | } |
2126 | 8.11M | |
2127 | 8.11M | /* Propagate carries as far as necessary */ |
2128 | 13.6M | for (/* */; pos < size_a13.6M ; ++pos, ++da, ++dc5.51M ) { |
2129 | 5.51M | w = w + *da; |
2130 | 5.51M | |
2131 | 5.51M | *dc = LOWER_HALF(w); |
2132 | 5.51M | w = UPPER_HALF(w); |
2133 | 5.51M | } |
2134 | 8.11M | |
2135 | 8.11M | /* Return carry out */ |
2136 | 8.11M | return (mp_digit)w; |
2137 | 8.11M | } |
2138 | | |
2139 | | STATIC void s_usub(mp_digit *da, mp_digit *db, mp_digit *dc, |
2140 | | mp_size size_a, mp_size size_b) |
2141 | 85.6M | { |
2142 | 85.6M | mp_size pos; |
2143 | 85.6M | mp_word w = 0; |
2144 | 85.6M | |
2145 | 85.6M | /* We assume that |a| >= |b| so this should definitely hold */ |
2146 | 85.6M | assert(size_a >= size_b); |
2147 | 85.6M | |
2148 | 85.6M | /* Subtract corresponding digits and propagate borrow */ |
2149 | 226M | for (pos = 0; pos < size_b226M ; ++pos, ++da, ++db, ++dc140M ) { |
2150 | 140M | w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */ |
2151 | 140M | (mp_word)*da) - w - (mp_word)*db; |
2152 | 140M | |
2153 | 140M | *dc = LOWER_HALF(w); |
2154 | 140M | w = (UPPER_HALF(w) == 0); |
2155 | 140M | } |
2156 | 85.6M | |
2157 | 85.6M | /* Finish the subtraction for remaining upper digits of da */ |
2158 | 170M | for (/* */; pos < size_a170M ; ++pos, ++da, ++dc85.1M ) { |
2159 | 85.1M | w = ((mp_word)MP_DIGIT_MAX + 1 + /* MP_RADIX */ |
2160 | 85.1M | (mp_word)*da) - w; |
2161 | 85.1M | |
2162 | 85.1M | *dc = LOWER_HALF(w); |
2163 | 85.1M | w = (UPPER_HALF(w) == 0); |
2164 | 85.1M | } |
2165 | 85.6M | |
2166 | 85.6M | /* If there is a borrow out at the end, it violates the precondition */ |
2167 | 85.6M | assert(w == 0); |
2168 | 85.6M | } |
2169 | | |
2170 | | STATIC int s_kmul(mp_digit *da, mp_digit *db, mp_digit *dc, |
2171 | | mp_size size_a, mp_size size_b) |
2172 | 7.48M | { |
2173 | 7.48M | mp_size bot_size; |
2174 | 7.48M | |
2175 | 7.48M | /* Make sure b is the smaller of the two input values */ |
2176 | 7.48M | if (size_b > size_a7.48M ) { |
2177 | 3.03M | SWAP(mp_digit *, da, db); |
2178 | 3.03M | SWAP(mp_size, size_a, size_b); |
2179 | 3.03M | } |
2180 | 7.48M | |
2181 | 7.48M | /* Insure that the bottom is the larger half in an odd-length split; the code |
2182 | 7.48M | below relies on this being true. |
2183 | 7.48M | */ |
2184 | 7.48M | bot_size = (size_a + 1) / 2; |
2185 | 7.48M | |
2186 | 7.48M | /* If the values are big enough to bother with recursion, use the Karatsuba |
2187 | 7.48M | algorithm to compute the product; otherwise use the normal multiplication |
2188 | 7.48M | algorithm |
2189 | 7.48M | */ |
2190 | 7.48M | if (multiply_threshold && |
2191 | 7.48M | size_a >= multiply_threshold && |
2192 | 7.48M | size_b > bot_size2.41k ) { |
2193 | 36 | |
2194 | 36 | mp_digit *t1, *t2, *t3, carry; |
2195 | 36 | |
2196 | 36 | mp_digit *a_top = da + bot_size; |
2197 | 36 | mp_digit *b_top = db + bot_size; |
2198 | 36 | |
2199 | 36 | mp_size at_size = size_a - bot_size; |
2200 | 36 | mp_size bt_size = size_b - bot_size; |
2201 | 36 | mp_size buf_size = 2 * bot_size; |
2202 | 36 | |
2203 | 36 | /* Do a single allocation for all three temporary buffers needed; each |
2204 | 36 | buffer must be big enough to hold the product of two bottom halves, and |
2205 | 36 | one buffer needs space for the completed product; twice the space is |
2206 | 36 | plenty. |
2207 | 36 | */ |
2208 | 36 | if ((t1 = s_alloc(4 * buf_size)) == NULL36 ) return 00 ; |
2209 | 36 | t2 = t1 + buf_size; |
2210 | 36 | t3 = t2 + buf_size; |
2211 | 36 | ZERO(t1, 4 * buf_size); |
2212 | 36 | |
2213 | 36 | /* t1 and t2 are initially used as temporaries to compute the inner product |
2214 | 36 | (a1 + a0)(b1 + b0) = a1b1 + a1b0 + a0b1 + a0b0 |
2215 | 36 | */ |
2216 | 36 | carry = s_uadd(da, a_top, t1, bot_size, at_size); /* t1 = a1 + a0 */ |
2217 | 36 | t1[bot_size] = carry; |
2218 | 36 | |
2219 | 36 | carry = s_uadd(db, b_top, t2, bot_size, bt_size); /* t2 = b1 + b0 */ |
2220 | 36 | t2[bot_size] = carry; |
2221 | 36 | |
2222 | 36 | (void) s_kmul(t1, t2, t3, bot_size + 1, bot_size + 1); /* t3 = t1 * t2 */ |
2223 | 36 | |
2224 | 36 | /* Now we'll get t1 = a0b0 and t2 = a1b1, and subtract them out so that |
2225 | 36 | we're left with only the pieces we want: t3 = a1b0 + a0b1 |
2226 | 36 | */ |
2227 | 36 | ZERO(t1, buf_size); |
2228 | 36 | ZERO(t2, buf_size); |
2229 | 36 | (void) s_kmul(da, db, t1, bot_size, bot_size); /* t1 = a0 * b0 */ |
2230 | 36 | (void) s_kmul(a_top, b_top, t2, at_size, bt_size); /* t2 = a1 * b1 */ |
2231 | 36 | |
2232 | 36 | /* Subtract out t1 and t2 to get the inner product */ |
2233 | 36 | s_usub(t3, t1, t3, buf_size + 2, buf_size); |
2234 | 36 | s_usub(t3, t2, t3, buf_size + 2, buf_size); |
2235 | 36 | |
2236 | 36 | /* Assemble the output value */ |
2237 | 36 | COPY(t1, dc, buf_size); |
2238 | 36 | carry = s_uadd(t3, dc + bot_size, dc + bot_size, |
2239 | 36 | buf_size + 1, buf_size); |
2240 | 36 | assert(carry == 0); |
2241 | 36 | |
2242 | 36 | carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size, |
2243 | 36 | buf_size, buf_size); |
2244 | 36 | assert(carry == 0); |
2245 | 36 | |
2246 | 36 | s_free(t1); /* note t2 and t3 are just internal pointers to t1 */ |
2247 | 36 | } |
2248 | 7.48M | else { |
2249 | 7.48M | s_umul(da, db, dc, size_a, size_b); |
2250 | 7.48M | } |
2251 | 7.48M | |
2252 | 7.48M | return 1; |
2253 | 7.48M | } |
2254 | | |
2255 | | STATIC void s_umul(mp_digit *da, mp_digit *db, mp_digit *dc, |
2256 | | mp_size size_a, mp_size size_b) |
2257 | 7.48M | { |
2258 | 7.48M | mp_size a, b; |
2259 | 7.48M | mp_word w; |
2260 | 7.48M | |
2261 | 25.7M | for (a = 0; a < size_a25.7M ; ++a, ++dc, ++da18.2M ) { |
2262 | 18.2M | mp_digit *dct = dc; |
2263 | 18.2M | mp_digit *dbt = db; |
2264 | 18.2M | |
2265 | 18.2M | if (*da == 0) |
2266 | 1.31M | continue; |
2267 | 16.9M | |
2268 | 16.9M | w = 0; |
2269 | 40.2M | for (b = 0; b < size_b40.2M ; ++b, ++dbt, ++dct23.3M ) { |
2270 | 23.3M | w = (mp_word)*da * (mp_word)*dbt + w + (mp_word)*dct; |
2271 | 23.3M | |
2272 | 23.3M | *dct = LOWER_HALF(w); |
2273 | 23.3M | w = UPPER_HALF(w); |
2274 | 23.3M | } |
2275 | 18.2M | |
2276 | 18.2M | *dct = (mp_digit)w; |
2277 | 18.2M | } |
2278 | 7.48M | } |
2279 | | |
2280 | | STATIC int s_ksqr(mp_digit *da, mp_digit *dc, mp_size size_a) |
2281 | 0 | { |
2282 | 0 | if (multiply_threshold && 0 size_a > multiply_threshold0 ) { |
2283 | 0 | mp_size bot_size = (size_a + 1) / 2; |
2284 | 0 | mp_digit *a_top = da + bot_size; |
2285 | 0 | mp_digit *t1, *t2, *t3, carry; |
2286 | 0 | mp_size at_size = size_a - bot_size; |
2287 | 0 | mp_size buf_size = 2 * bot_size; |
2288 | 0 |
|
2289 | 0 | if ((t1 = s_alloc(4 * buf_size)) == NULL0 ) return 00 ; |
2290 | 0 | t2 = t1 + buf_size; |
2291 | 0 | t3 = t2 + buf_size; |
2292 | 0 | ZERO(t1, 4 * buf_size); |
2293 | 0 |
|
2294 | 0 | (void) s_ksqr(da, t1, bot_size); /* t1 = a0 ^ 2 */ |
2295 | 0 | (void) s_ksqr(a_top, t2, at_size); /* t2 = a1 ^ 2 */ |
2296 | 0 |
|
2297 | 0 | (void) s_kmul(da, a_top, t3, bot_size, at_size); /* t3 = a0 * a1 */ |
2298 | 0 |
|
2299 | 0 | /* Quick multiply t3 by 2, shifting left (can't overflow) */ |
2300 | 0 | { |
2301 | 0 | int i, top = bot_size + at_size; |
2302 | 0 | mp_word w, save = 0; |
2303 | 0 |
|
2304 | 0 | for (i = 0; i < top0 ; ++i0 ) { |
2305 | 0 | w = t3[i]; |
2306 | 0 | w = (w << 1) | save; |
2307 | 0 | t3[i] = LOWER_HALF(w); |
2308 | 0 | save = UPPER_HALF(w); |
2309 | 0 | } |
2310 | 0 | t3[i] = LOWER_HALF(save); |
2311 | 0 | } |
2312 | 0 |
|
2313 | 0 | /* Assemble the output value */ |
2314 | 0 | COPY(t1, dc, 2 * bot_size); |
2315 | 0 | carry = s_uadd(t3, dc + bot_size, dc + bot_size, |
2316 | 0 | buf_size + 1, buf_size); |
2317 | 0 | assert(carry == 0); |
2318 | 0 |
|
2319 | 0 | carry = s_uadd(t2, dc + 2*bot_size, dc + 2*bot_size, |
2320 | 0 | buf_size, buf_size); |
2321 | 0 | assert(carry == 0); |
2322 | 0 |
|
2323 | 0 | s_free(t1); /* note that t2 and t2 are internal pointers only */ |
2324 | 0 |
|
2325 | 0 | } |
2326 | 0 | else { |
2327 | 0 | s_usqr(da, dc, size_a); |
2328 | 0 | } |
2329 | 0 |
|
2330 | 0 | return 1; |
2331 | 0 | } |
2332 | | |
2333 | | STATIC void s_usqr(mp_digit *da, mp_digit *dc, mp_size size_a) |
2334 | 0 | { |
2335 | 0 | mp_size i, j; |
2336 | 0 | mp_word w; |
2337 | 0 |
|
2338 | 0 | for (i = 0; i < size_a0 ; ++i, dc += 2, ++da0 ) { |
2339 | 0 | mp_digit *dct = dc, *dat = da; |
2340 | 0 |
|
2341 | 0 | if (*da == 0) |
2342 | 0 | continue; |
2343 | 0 |
|
2344 | 0 | /* Take care of the first digit, no rollover */ |
2345 | 0 | w = (mp_word)*dat * (mp_word)*dat + (mp_word)*dct; |
2346 | 0 | *dct = LOWER_HALF(w); |
2347 | 0 | w = UPPER_HALF(w); |
2348 | 0 | ++dat; ++dct; |
2349 | 0 |
|
2350 | 0 | for (j = i + 1; j < size_a0 ; ++j, ++dat, ++dct0 ) { |
2351 | 0 | mp_word t = (mp_word)*da * (mp_word)*dat; |
2352 | 0 | mp_word u = w + (mp_word)*dct, ov = 0; |
2353 | 0 |
|
2354 | 0 | /* Check if doubling t will overflow a word */ |
2355 | 0 | if (HIGH_BIT_SET(t)) |
2356 | 0 | ov = 1; |
2357 | 0 |
|
2358 | 0 | w = t + t; |
2359 | 0 |
|
2360 | 0 | /* Check if adding u to w will overflow a word */ |
2361 | 0 | if (ADD_WILL_OVERFLOW(w, u)) |
2362 | 0 | ov = 1; |
2363 | 0 |
|
2364 | 0 | w += u; |
2365 | 0 |
|
2366 | 0 | *dct = LOWER_HALF(w); |
2367 | 0 | w = UPPER_HALF(w); |
2368 | 0 | if (ov0 ) { |
2369 | 0 | w += MP_DIGIT_MAX; /* MP_RADIX */ |
2370 | 0 | ++w; |
2371 | 0 | } |
2372 | 0 | } |
2373 | 0 |
|
2374 | 0 | w = w + *dct; |
2375 | 0 | *dct = (mp_digit)w; |
2376 | 0 | while ((w = 0 UPPER_HALF0 (w)) != 0) { |
2377 | 0 | ++dct; w = w + *dct; |
2378 | 0 | *dct = LOWER_HALF(w); |
2379 | 0 | } |
2380 | 0 |
|
2381 | 0 | assert(w == 0); |
2382 | 0 | } |
2383 | 0 | } |
2384 | | |
2385 | | STATIC void s_dadd(mp_int a, mp_digit b) |
2386 | 1.68k | { |
2387 | 1.68k | mp_word w = 0; |
2388 | 1.68k | mp_digit *da = MP_DIGITS(a); |
2389 | 1.68k | mp_size ua = MP_USED(a); |
2390 | 1.68k | |
2391 | 1.68k | w = (mp_word)*da + b; |
2392 | 1.68k | *da++ = LOWER_HALF(w); |
2393 | 1.68k | w = UPPER_HALF(w); |
2394 | 1.68k | |
2395 | 1.88k | for (ua -= 1; ua > 01.88k ; --ua, ++da200 ) { |
2396 | 200 | w = (mp_word)*da + w; |
2397 | 200 | |
2398 | 200 | *da = LOWER_HALF(w); |
2399 | 200 | w = UPPER_HALF(w); |
2400 | 200 | } |
2401 | 1.68k | |
2402 | 1.68k | if (w1.68k ) { |
2403 | 2 | *da = (mp_digit)w; |
2404 | 2 | MP_USED(a) += 1; |
2405 | 2 | } |
2406 | 1.68k | } |
2407 | | |
2408 | | STATIC void s_dmul(mp_int a, mp_digit b) |
2409 | 1.68k | { |
2410 | 1.68k | mp_word w = 0; |
2411 | 1.68k | mp_digit *da = MP_DIGITS(a); |
2412 | 1.68k | mp_size ua = MP_USED(a); |
2413 | 1.68k | |
2414 | 3.54k | while (ua > 03.54k ) { |
2415 | 1.86k | w = (mp_word)*da * b + w; |
2416 | 1.86k | *da++ = LOWER_HALF(w); |
2417 | 1.86k | w = UPPER_HALF(w); |
2418 | 1.86k | --ua; |
2419 | 1.86k | } |
2420 | 1.68k | |
2421 | 1.68k | if (w1.68k ) { |
2422 | 20 | *da = (mp_digit)w; |
2423 | 20 | MP_USED(a) += 1; |
2424 | 20 | } |
2425 | 1.68k | } |
2426 | | |
2427 | | STATIC void s_dbmul(mp_digit *da, mp_digit b, mp_digit *dc, mp_size size_a) |
2428 | 2.08M | { |
2429 | 2.08M | mp_word w = 0; |
2430 | 2.08M | |
2431 | 10.6M | while (size_a > 010.6M ) { |
2432 | 8.54M | w = (mp_word)*da++ * (mp_word)b + w; |
2433 | 8.54M | |
2434 | 8.54M | *dc++ = LOWER_HALF(w); |
2435 | 8.54M | w = UPPER_HALF(w); |
2436 | 8.54M | --size_a; |
2437 | 8.54M | } |
2438 | 2.08M | |
2439 | 2.08M | if (w) |
2440 | 0 | *dc = 0 LOWER_HALF0 (w); |
2441 | 2.08M | } |
2442 | | |
2443 | | STATIC mp_digit s_ddiv(mp_int a, mp_digit b) |
2444 | 899k | { |
2445 | 899k | mp_word w = 0, qdigit; |
2446 | 899k | mp_size ua = MP_USED(a); |
2447 | 899k | mp_digit *da = MP_DIGITS(a) + ua - 1; |
2448 | 899k | |
2449 | 3.05M | for (/* */; ua > 03.05M ; --ua, --da2.15M ) { |
2450 | 2.15M | w = (w << MP_DIGIT_BIT) | *da; |
2451 | 2.15M | |
2452 | 2.15M | if (w >= b2.15M ) { |
2453 | 1.62M | qdigit = w / b; |
2454 | 1.62M | w = w % b; |
2455 | 1.62M | } |
2456 | 531k | else { |
2457 | 531k | qdigit = 0; |
2458 | 531k | } |
2459 | 2.15M | |
2460 | 2.15M | *da = (mp_digit)qdigit; |
2461 | 2.15M | } |
2462 | 899k | |
2463 | 899k | CLAMP(a); |
2464 | 899k | return (mp_digit)w; |
2465 | 899k | } |
2466 | | |
2467 | | STATIC void s_qdiv(mp_int z, mp_size p2) |
2468 | 86.0M | { |
2469 | 86.0M | mp_size ndig = p2 / MP_DIGIT_BIT86.0M , nbits = p2 % MP_DIGIT_BIT86.0M ; |
2470 | 86.0M | mp_size uz = MP_USED(z); |
2471 | 86.0M | |
2472 | 86.0M | if (ndig86.0M ) { |
2473 | 262k | mp_size mark; |
2474 | 262k | mp_digit *to, *from; |
2475 | 262k | |
2476 | 262k | if (ndig >= uz262k ) { |
2477 | 0 | mp_int_zero(z); |
2478 | 0 | return; |
2479 | 0 | } |
2480 | 262k | |
2481 | 262k | to = 262k MP_DIGITS262k (z); from = to + ndig; |
2482 | 262k | |
2483 | 608k | for (mark = ndig; mark < uz608k ; ++mark345k ) |
2484 | 345k | *to++ = *from++; |
2485 | 262k | |
2486 | 262k | MP_USED(z) = uz - ndig; |
2487 | 262k | } |
2488 | 86.0M | |
2489 | 86.0M | if (86.0M nbits86.0M ) { |
2490 | 81.3M | mp_digit d = 0, *dz, save; |
2491 | 81.3M | mp_size up = MP_DIGIT_BIT - nbits; |
2492 | 81.3M | |
2493 | 81.3M | uz = MP_USED(z); |
2494 | 81.3M | dz = MP_DIGITS(z) + uz - 1; |
2495 | 81.3M | |
2496 | 295M | for (/* */; uz > 0295M ; --uz, --dz214M ) { |
2497 | 214M | save = *dz; |
2498 | 214M | |
2499 | 214M | *dz = (*dz >> nbits) | (d << up); |
2500 | 214M | d = save; |
2501 | 214M | } |
2502 | 81.3M | |
2503 | 81.3M | CLAMP(z); |
2504 | 81.3M | } |
2505 | 86.0M | |
2506 | 86.0M | if (MP_USED86.0M (z) == 1 && 86.0M z->digits[0] == 033.3M ) |
2507 | 779k | MP_SIGN779k (z) = MP_ZPOS779k ; |
2508 | 86.0M | } |
2509 | | |
2510 | | STATIC void s_qmod(mp_int z, mp_size p2) |
2511 | 138k | { |
2512 | 138k | mp_size start = p2 / MP_DIGIT_BIT138k + 1, rest = p2 % MP_DIGIT_BIT138k ; |
2513 | 138k | mp_size uz = MP_USED(z); |
2514 | 138k | mp_digit mask = (1u << rest) - 1; |
2515 | 138k | |
2516 | 138k | if (start <= uz138k ) { |
2517 | 138k | MP_USED(z) = start; |
2518 | 138k | z->digits[start - 1] &= mask; |
2519 | 138k | CLAMP(z); |
2520 | 138k | } |
2521 | 138k | } |
2522 | | |
2523 | | STATIC int s_qmul(mp_int z, mp_size p2) |
2524 | 4.42M | { |
2525 | 4.42M | mp_size uz, need, rest, extra, i; |
2526 | 4.42M | mp_digit *from, *to, d; |
2527 | 4.42M | |
2528 | 4.42M | if (p2 == 0) |
2529 | 1.46M | return 1; |
2530 | 2.95M | |
2531 | 2.95M | uz = 2.95M MP_USED2.95M (z); |
2532 | 2.95M | need = p2 / MP_DIGIT_BIT2.95M ; rest = p2 % MP_DIGIT_BIT2.95M ; |
2533 | 2.95M | |
2534 | 2.95M | /* Figure out if we need an extra digit at the top end; this occurs if the |
2535 | 2.95M | topmost `rest' bits of the high-order digit of z are not zero, meaning |
2536 | 2.95M | they will be shifted off the end if not preserved */ |
2537 | 2.95M | extra = 0; |
2538 | 2.95M | if (rest != 02.95M ) { |
2539 | 2.89M | mp_digit *dz = MP_DIGITS(z) + uz - 1; |
2540 | 2.89M | |
2541 | 2.89M | if ((*dz >> (2.89M MP_DIGIT_BIT2.89M - rest)) != 0) |
2542 | 582k | extra = 1; |
2543 | 2.89M | } |
2544 | 2.95M | |
2545 | 2.95M | if (!s_pad(z, uz + need + extra)) |
2546 | 0 | return 0; |
2547 | 2.95M | |
2548 | 2.95M | /* If we need to shift by whole digits, do that in one pass, then |
2549 | 2.95M | to back and shift by partial digits. |
2550 | 2.95M | */ |
2551 | 2.95M | if (2.95M need > 02.95M ) { |
2552 | 87.5k | from = MP_DIGITS(z) + uz - 1; |
2553 | 87.5k | to = from + need; |
2554 | 87.5k | |
2555 | 175k | for (i = 0; i < uz175k ; ++i87.7k ) |
2556 | 87.7k | *to-- = *from--; |
2557 | 87.5k | |
2558 | 87.5k | ZERO(MP_DIGITS(z), need); |
2559 | 87.5k | uz += need; |
2560 | 87.5k | } |
2561 | 2.95M | |
2562 | 2.95M | if (rest2.95M ) { |
2563 | 2.89M | d = 0; |
2564 | 10.4M | for (i = need, from = MP_DIGITS2.89M (z) + need; i < uz10.4M ; ++i, ++from7.55M ) { |
2565 | 7.55M | mp_digit save = *from; |
2566 | 7.55M | |
2567 | 7.55M | *from = (*from << rest) | (d >> (MP_DIGIT_BIT - rest)); |
2568 | 7.55M | d = save; |
2569 | 7.55M | } |
2570 | 2.89M | |
2571 | 2.89M | d >>= (MP_DIGIT_BIT - rest); |
2572 | 2.89M | if (d != 02.89M ) { |
2573 | 582k | *from = d; |
2574 | 582k | uz += extra; |
2575 | 582k | } |
2576 | 2.89M | } |
2577 | 2.95M | |
2578 | 2.95M | MP_USED(z) = uz; |
2579 | 2.95M | CLAMP(z); |
2580 | 4.42M | |
2581 | 4.42M | return 1; |
2582 | 4.42M | } |
2583 | | |
2584 | | /* Compute z = 2^p2 - |z|; requires that 2^p2 >= |z| |
2585 | | The sign of the result is always zero/positive. |
2586 | | */ |
2587 | | STATIC int s_qsub(mp_int z, mp_size p2) |
2588 | 0 | { |
2589 | 0 | mp_digit hi = (1 << (p2 % MP_DIGIT_BIT)), *zp; |
2590 | 0 | mp_size tdig = (p2 / MP_DIGIT_BIT), pos; |
2591 | 0 | mp_word w = 0; |
2592 | 0 |
|
2593 | 0 | if (!s_pad(z, tdig + 1)) |
2594 | 0 | return 0; |
2595 | 0 |
|
2596 | 0 | for (pos = 0, zp = 0 MP_DIGITS0 (z); pos < tdig0 ; ++pos, ++zp0 ) { |
2597 | 0 | w = ((mp_word) MP_DIGIT_MAX + 1) - w - (mp_word)*zp; |
2598 | 0 |
|
2599 | 0 | *zp = LOWER_HALF(w); |
2600 | 0 | w = UPPER_HALF0 (w) ? 00 : 10 ; |
2601 | 0 | } |
2602 | 0 |
|
2603 | 0 | w = ((mp_word) MP_DIGIT_MAX + 1 + hi) - w - (mp_word)*zp; |
2604 | 0 | *zp = LOWER_HALF(w); |
2605 | 0 |
|
2606 | 0 | assert(UPPER_HALF(w) != 0); /* no borrow out should be possible */ |
2607 | 0 |
|
2608 | 0 | MP_SIGN(z) = MP_ZPOS; |
2609 | 0 | CLAMP(z); |
2610 | 0 |
|
2611 | 0 | return 1; |
2612 | 0 | } |
2613 | | |
2614 | | STATIC int s_dp2k(mp_int z) |
2615 | 84.2M | { |
2616 | 84.2M | int k = 0; |
2617 | 84.2M | mp_digit *dp = MP_DIGITS(z), d; |
2618 | 84.2M | |
2619 | 84.2M | if (MP_USED84.2M (z) == 1 && 84.2M *dp == 030.2M ) |
2620 | 0 | return 1; |
2621 | 84.2M | |
2622 | 84.7M | while (84.2M *dp == 084.7M ) { |
2623 | 434k | k += MP_DIGIT_BIT; |
2624 | 434k | ++dp; |
2625 | 434k | } |
2626 | 84.2M | |
2627 | 84.2M | d = *dp; |
2628 | 245M | while ((d & 1) == 0245M ) { |
2629 | 161M | d >>= 1; |
2630 | 161M | ++k; |
2631 | 161M | } |
2632 | 84.2M | |
2633 | 84.2M | return k; |
2634 | 84.2M | } |
2635 | | |
2636 | | STATIC int s_isp2(mp_int z) |
2637 | 2.68M | { |
2638 | 2.68M | mp_size uz = MP_USED(z), k = 0; |
2639 | 2.68M | mp_digit *dz = MP_DIGITS(z), d; |
2640 | 2.68M | |
2641 | 2.73M | while (uz > 12.73M ) { |
2642 | 898k | if (*dz++ != 0) |
2643 | 840k | return -1; |
2644 | 57.1k | k += 57.1k MP_DIGIT_BIT57.1k ; |
2645 | 898k | --uz; |
2646 | 898k | } |
2647 | 2.68M | |
2648 | 1.84M | d = *dz; |
2649 | 5.24M | while (d > 15.24M ) { |
2650 | 4.27M | if (d & 1) |
2651 | 875k | return -1; |
2652 | 3.39M | ++k; d >>= 1; |
2653 | 3.39M | } |
2654 | 1.84M | |
2655 | 966k | return (int) k; |
2656 | 2.68M | } |
2657 | | |
2658 | | STATIC int s_2expt(mp_int z, mp_small k) |
2659 | 0 | { |
2660 | 0 | mp_size ndig, rest; |
2661 | 0 | mp_digit *dz; |
2662 | 0 |
|
2663 | 0 | ndig = (k + MP_DIGIT_BIT0 ) / MP_DIGIT_BIT0 ; |
2664 | 0 | rest = k % MP_DIGIT_BIT; |
2665 | 0 |
|
2666 | 0 | if (!s_pad(z, ndig)) |
2667 | 0 | return 0; |
2668 | 0 |
|
2669 | 0 | dz = 0 MP_DIGITS0 (z); |
2670 | 0 | ZERO(dz, ndig); |
2671 | 0 | *(dz + ndig - 1) = (1 << rest); |
2672 | 0 | MP_USED(z) = ndig; |
2673 | 0 |
|
2674 | 0 | return 1; |
2675 | 0 | } |
2676 | | |
2677 | | STATIC int s_norm(mp_int a, mp_int b) |
2678 | 842k | { |
2679 | 842k | mp_digit d = b->digits[MP_USED(b) - 1]; |
2680 | 842k | int k = 0; |
2681 | 842k | |
2682 | 15.0M | while (d < (1u << (mp_digit)(15.0M MP_DIGIT_BIT15.0M - 1))) { /* d < (MP_RADIX / 2) */ |
2683 | 14.2M | d <<= 1; |
2684 | 14.2M | ++k; |
2685 | 14.2M | } |
2686 | 842k | |
2687 | 842k | /* These multiplications can't fail */ |
2688 | 842k | if (k != 0842k ) { |
2689 | 827k | (void) s_qmul(a, (mp_size) k); |
2690 | 827k | (void) s_qmul(b, (mp_size) k); |
2691 | 827k | } |
2692 | 842k | |
2693 | 842k | return k; |
2694 | 842k | } |
2695 | | |
2696 | | STATIC mp_result s_brmu(mp_int z, mp_int m) |
2697 | 0 | { |
2698 | 0 | mp_size um = MP_USED(m) * 2; |
2699 | 0 |
|
2700 | 0 | if (!s_pad(z, um)) |
2701 | 0 | return MP_MEMORY; |
2702 | 0 |
|
2703 | 0 | s_2expt(z, 0 MP_DIGIT_BIT0 * um); |
2704 | 0 | return mp_int_div(z, m, z, NULL); |
2705 | 0 | } |
2706 | | |
2707 | | STATIC int s_reduce(mp_int x, mp_int m, mp_int mu, mp_int q1, mp_int q2) |
2708 | 0 | { |
2709 | 0 | mp_size um = MP_USED(m), umb_p1, umb_m1; |
2710 | 0 |
|
2711 | 0 | umb_p1 = (um + 1) * MP_DIGIT_BIT; |
2712 | 0 | umb_m1 = (um - 1) * MP_DIGIT_BIT; |
2713 | 0 |
|
2714 | 0 | if (mp_int_copy(x, q1) != MP_OK) |
2715 | 0 | return 0; |
2716 | 0 |
|
2717 | 0 | /* Compute q2 = floor((floor(x / b^(k-1)) * mu) / b^(k+1)) */ |
2718 | 0 | s_qdiv(q1, umb_m1); |
2719 | 0 | UMUL(q1, mu, q2); |
2720 | 0 | s_qdiv(q2, umb_p1); |
2721 | 0 |
|
2722 | 0 | /* Set x = x mod b^(k+1) */ |
2723 | 0 | s_qmod(x, umb_p1); |
2724 | 0 |
|
2725 | 0 | /* Now, q is a guess for the quotient a / m. |
2726 | 0 | Compute x - q * m mod b^(k+1), replacing x. This may be off |
2727 | 0 | by a factor of 2m, but no more than that. |
2728 | 0 | */ |
2729 | 0 | UMUL(q2, m, q1); |
2730 | 0 | s_qmod(q1, umb_p1); |
2731 | 0 | (void) mp_int_sub(x, q1, x); /* can't fail */ |
2732 | 0 |
|
2733 | 0 | /* The result may be < 0; if it is, add b^(k+1) to pin it in the proper |
2734 | 0 | range. */ |
2735 | 0 | if ((0 CMPZ0 (x) < 0) && !s_qsub(x, umb_p1)0 ) |
2736 | 0 | return 0; |
2737 | 0 |
|
2738 | 0 | /* If x > m, we need to back it off until it is in range. This will be |
2739 | 0 | required at most twice. */ |
2740 | 0 | if (0 mp_int_compare(x, m) >= 00 ) { |
2741 | 0 | (void) mp_int_sub(x, m, x); |
2742 | 0 | if (mp_int_compare(x, m) >= 0) |
2743 | 0 | (void) mp_int_sub(x, m, x); |
2744 | 0 | } |
2745 | 0 |
|
2746 | 0 | /* At this point, x has been properly reduced. */ |
2747 | 0 | return 1; |
2748 | 0 | } |
2749 | | |
2750 | | /* Perform modular exponentiation using Barrett's method, where mu is the |
2751 | | reduction constant for m. Assumes a < m, b > 0. */ |
2752 | | STATIC mp_result s_embar(mp_int a, mp_int b, mp_int m, mp_int mu, mp_int c) |
2753 | 0 | { |
2754 | 0 | mp_digit *db, *dbt, umu, d; |
2755 | 0 | mp_result res; |
2756 | 0 | DECLARE_TEMP(3); |
2757 | 0 |
|
2758 | 0 | umu = MP_USED0 (mu); db = MP_DIGITS0 (b); dbt = db + MP_USED0 (b) - 1; |
2759 | 0 |
|
2760 | 0 | while (last__ < 30 ) { |
2761 | 0 | SETUP(mp_int_init_size(LAST_TEMP(), 4 * umu)); |
2762 | 0 | ZERO0 (MP_DIGITS(TEMP(last__ - 1)), MP_ALLOC(TEMP(last__ - 1))); |
2763 | 0 | } |
2764 | 0 |
|
2765 | 0 | (void) mp_int_set_value(c, 1); |
2766 | 0 |
|
2767 | 0 | /* Take care of low-order digits */ |
2768 | 0 | while (db < dbt0 ) { |
2769 | 0 | int i; |
2770 | 0 |
|
2771 | 0 | for (d = *db, i = MP_DIGIT_BIT0 ; i > 00 ; --i, d >>= 10 ) { |
2772 | 0 | if (d & 10 ) { |
2773 | 0 | /* The use of a second temporary avoids allocation */ |
2774 | 0 | UMUL(c, a, TEMP(0)); |
2775 | 0 | if (!s_reduce(0 TEMP0 (0), m, mu, TEMP0 (1), TEMP0 (2))) { |
2776 | 0 | res = MP_MEMORY; goto CLEANUP; |
2777 | 0 | } |
2778 | 0 | mp_int_copy(0 TEMP0 (0), c); |
2779 | 0 | } |
2780 | 0 |
|
2781 | 0 |
|
2782 | 0 | USQR0 (a, TEMP(0)); |
2783 | 0 | assert(MP_SIGN(TEMP(0)) == MP_ZPOS); |
2784 | 0 | if (!s_reduce(0 TEMP0 (0), m, mu, TEMP0 (1), TEMP0 (2))) { |
2785 | 0 | res = MP_MEMORY; goto CLEANUP; |
2786 | 0 | } |
2787 | 0 | assert(MP_SIGN(TEMP(0)) == MP_ZPOS); |
2788 | 0 | mp_int_copy(TEMP(0), a); |
2789 | 0 | } |
2790 | 0 |
|
2791 | 0 | ++db; |
2792 | 0 | } |
2793 | 0 |
|
2794 | 0 | /* Take care of highest-order digit */ |
2795 | 0 | d = *dbt; |
2796 | 0 | for (;;) { |
2797 | 0 | if (d & 10 ) { |
2798 | 0 | UMUL(c, a, TEMP(0)); |
2799 | 0 | if (!s_reduce(0 TEMP0 (0), m, mu, TEMP0 (1), TEMP0 (2))) { |
2800 | 0 | res = MP_MEMORY; goto CLEANUP; |
2801 | 0 | } |
2802 | 0 | mp_int_copy(0 TEMP0 (0), c); |
2803 | 0 | } |
2804 | 0 |
|
2805 | 0 | d >>= 1; |
2806 | 0 | if (!d0 ) break0 ; |
2807 | 0 |
|
2808 | 0 | USQR0 (a, TEMP(0)); |
2809 | 0 | if (!s_reduce(0 TEMP0 (0), m, mu, TEMP0 (1), TEMP0 (2))) { |
2810 | 0 | res = MP_MEMORY; goto CLEANUP; |
2811 | 0 | } |
2812 | 0 | (void) mp_int_copy(0 TEMP0 (0), a); |
2813 | 0 | } |
2814 | 0 |
|
2815 | 0 | CLEANUP_TEMP0 (); |
2816 | 0 | return res; |
2817 | 0 | } |
2818 | | |
2819 | | /* Division of nonnegative integers |
2820 | | |
2821 | | This function implements division algorithm for unsigned multi-precision |
2822 | | integers. The algorithm is based on Algorithm D from Knuth's "The Art of |
2823 | | Computer Programming", 3rd ed. 1998, pg 272-273. |
2824 | | |
2825 | | We diverge from Knuth's algorithm in that we do not perform the subtraction |
2826 | | from the remainder until we have determined that we have the correct |
2827 | | quotient digit. This makes our algorithm less efficient that Knuth because |
2828 | | we might have to perform multiple multiplication and comparison steps before |
2829 | | the subtraction. The advantage is that it is easy to implement and ensure |
2830 | | correctness without worrying about underflow from the subtraction. |
2831 | | |
2832 | | inputs: u a n+m digit integer in base b (b is 2^MP_DIGIT_BIT) |
2833 | | v a n digit integer in base b (b is 2^MP_DIGIT_BIT) |
2834 | | n >= 1 |
2835 | | m >= 0 |
2836 | | outputs: u / v stored in u |
2837 | | u % v stored in v |
2838 | | */ |
2839 | 1.71M | STATIC mp_result s_udiv_knuth(mp_int u, mp_int v) { |
2840 | 1.71M | mpz_t q, r, t; |
2841 | 1.71M | mp_result |
2842 | 1.71M | res = MP_OK; |
2843 | 1.71M | int k,j; |
2844 | 1.71M | mp_size m,n; |
2845 | 1.71M | |
2846 | 1.71M | /* Force signs to positive */ |
2847 | 1.71M | MP_SIGN(u) = MP_ZPOS; |
2848 | 1.71M | MP_SIGN(v) = MP_ZPOS; |
2849 | 1.71M | |
2850 | 1.71M | /* Use simple division algorithm when v is only one digit long */ |
2851 | 1.71M | if (MP_USED1.71M (v) == 11.71M ) { |
2852 | 873k | mp_digit d, rem; |
2853 | 873k | d = v->digits[0]; |
2854 | 873k | rem = s_ddiv(u, d); |
2855 | 873k | mp_int_set_value(v, rem); |
2856 | 873k | return MP_OK; |
2857 | 873k | } |
2858 | 842k | |
2859 | 842k | /* Algorithm D |
2860 | 842k | |
2861 | 842k | The n and m variables are defined as used by Knuth. |
2862 | 842k | u is an n digit number with digits u_{n-1}..u_0. |
2863 | 842k | v is an n+m digit number with digits from v_{m+n-1}..v_0. |
2864 | 842k | We require that n > 1 and m >= 0 |
2865 | 842k | */ |
2866 | 842k | n = 842k MP_USED842k (v); |
2867 | 842k | m = MP_USED(u) - n; |
2868 | 842k | assert(n > 1); |
2869 | 842k | assert(m >= 0); |
2870 | 842k | |
2871 | 842k | /* D1: Normalize. |
2872 | 842k | The normalization step provides the necessary condition for Theorem B, |
2873 | 842k | which states that the quotient estimate for q_j, call it qhat |
2874 | 842k | |
2875 | 842k | qhat = u_{j+n}u_{j+n-1} / v_{n-1} |
2876 | 842k | |
2877 | 842k | is bounded by |
2878 | 842k | |
2879 | 842k | qhat - 2 <= q_j <= qhat. |
2880 | 842k | |
2881 | 842k | That is, qhat is always greater than the actual quotient digit q, |
2882 | 842k | and it is never more than two larger than the actual quotient digit. |
2883 | 842k | */ |
2884 | 842k | k = s_norm(u, v); |
2885 | 842k | |
2886 | 842k | /* Extend size of u by one if needed. |
2887 | 842k | |
2888 | 842k | The algorithm begins with a value of u that has one more digit of input. |
2889 | 842k | The normalization step sets u_{m+n}..u_0 = 2^k * u_{m+n-1}..u_0. If the |
2890 | 842k | multiplication did not increase the number of digits of u, we need to add |
2891 | 842k | a leading zero here. |
2892 | 842k | */ |
2893 | 842k | if (k == 0 || 842k MP_USED827k (u) != m + n + 1827k ) { |
2894 | 318k | if (!s_pad(u, m+n+1)) |
2895 | 0 | return MP_MEMORY; |
2896 | 318k | u->digits[m+n] = 0; |
2897 | 318k | u->used = m+n+1; |
2898 | 318k | } |
2899 | 842k | |
2900 | 842k | /* Add a leading 0 to v. |
2901 | 842k | |
2902 | 842k | The multiplication in step D4 multiplies qhat * 0v_{n-1}..v_0. We need to |
2903 | 842k | add the leading zero to v here to ensure that the multiplication will |
2904 | 842k | produce the full n+1 digit result. |
2905 | 842k | */ |
2906 | 842k | if (842k !s_pad(v, n+1)842k ) return MP_MEMORY0 ; v->digits[n] = 0; |
2907 | 842k | |
2908 | 842k | /* Initialize temporary variables q and t. |
2909 | 842k | q allocates space for m+1 digits to store the quotient digits |
2910 | 842k | t allocates space for n+1 digits to hold the result of q_j*v |
2911 | 842k | */ |
2912 | 842k | if ((res = mp_int_init_size(&q, m + 1)) != MP_OK842k ) return res0 ; |
2913 | 842k | if (842k (res = mp_int_init_size(&t, n + 1)) != MP_OK842k ) goto CLEANUP0 ; |
2914 | 842k | |
2915 | 842k | /* D2: Initialize j */ |
2916 | 842k | j = m; |
2917 | 842k | r.digits = MP_DIGITS(u) + j; /* The contents of r are shared with u */ |
2918 | 842k | r.used = n + 1; |
2919 | 842k | r.sign = MP_ZPOS; |
2920 | 842k | r.alloc = MP_ALLOC(u); |
2921 | 842k | ZERO(t.digits, t.alloc); |
2922 | 842k | |
2923 | 842k | /* Calculate the m+1 digits of the quotient result */ |
2924 | 2.76M | for (; j >= 02.76M ; j--1.92M ) { |
2925 | 1.92M | /* D3: Calculate q' */ |
2926 | 1.92M | /* r->digits is aligned to position j of the number u */ |
2927 | 1.92M | mp_word pfx, qhat; |
2928 | 1.92M | pfx = r.digits[n]; |
2929 | 1.92M | pfx <<= MP_DIGIT_BIT / 2; |
2930 | 1.92M | pfx <<= MP_DIGIT_BIT / 2; |
2931 | 1.92M | pfx |= r.digits[n-1]; /* pfx = u_{j+n}{j+n-1} */ |
2932 | 1.92M | |
2933 | 1.92M | qhat = pfx / v->digits[n-1]; |
2934 | 1.92M | /* Check to see if qhat > b, and decrease qhat if so. |
2935 | 1.92M | Theorem B guarantess that qhat is at most 2 larger than the |
2936 | 1.92M | actual value, so it is possible that qhat is greater than |
2937 | 1.92M | the maximum value that will fit in a digit */ |
2938 | 1.92M | if (qhat > 1.92M MP_DIGIT_MAX1.92M ) |
2939 | 86 | qhat = 86 MP_DIGIT_MAX86 ; |
2940 | 1.92M | |
2941 | 1.92M | /* D4,D5,D6: Multiply qhat * v and test for a correct value of q |
2942 | 1.92M | |
2943 | 1.92M | We proceed a bit different than the way described by Knuth. This way is |
2944 | 1.92M | simpler but less efficent. Instead of doing the multiply and subtract |
2945 | 1.92M | then checking for underflow, we first do the multiply of qhat * v and |
2946 | 1.92M | see if it is larger than the current remainder r. If it is larger, we |
2947 | 1.92M | decrease qhat by one and try again. We may need to decrease qhat one |
2948 | 1.92M | more time before we get a value that is smaller than r. |
2949 | 1.92M | |
2950 | 1.92M | This way is less efficent than Knuth becuase we do more multiplies, but |
2951 | 1.92M | we do not need to worry about underflow this way. |
2952 | 1.92M | */ |
2953 | 1.92M | /* t = qhat * v */ |
2954 | 1.92M | s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1); t.used = n + 1; |
2955 | 1.92M | CLAMP(&t); |
2956 | 1.92M | |
2957 | 1.92M | /* Clamp r for the comparison. Comparisons do not like leading zeros. */ |
2958 | 1.92M | CLAMP(&r); |
2959 | 1.92M | if (s_ucmp(&t, &r) > 01.92M ) { /* would the remainder be negative? */ |
2960 | 157k | qhat -= 1; /* try a smaller q */ |
2961 | 157k | s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1); |
2962 | 157k | t.used = n + 1; CLAMP(&t); |
2963 | 157k | if (s_ucmp(&t, &r) > 0157k ) { /* would the remainder be negative? */ |
2964 | 4.39k | assert(qhat > 0); |
2965 | 4.39k | qhat -= 1; /* try a smaller q */ |
2966 | 4.39k | s_dbmul(MP_DIGITS(v), (mp_digit) qhat, t.digits, n+1); |
2967 | 4.39k | t.used = n + 1; CLAMP(&t); |
2968 | 4.39k | } |
2969 | 157k | assert(s_ucmp(&t, &r) <= 0 && "The mathematics failed us."); |
2970 | 157k | } |
2971 | 1.92M | /* Unclamp r. The D algorithm expects r = u_{j+n}..u_j to always be n+1 |
2972 | 1.92M | digits long. */ |
2973 | 1.92M | r.used = n + 1; |
2974 | 1.92M | |
2975 | 1.92M | /* D4: Multiply and subtract |
2976 | 1.92M | |
2977 | 1.92M | Note: The multiply was completed above so we only need to subtract here. |
2978 | 1.92M | */ |
2979 | 1.92M | s_usub(r.digits, t.digits, r.digits, r.used, t.used); |
2980 | 1.92M | |
2981 | 1.92M | /* D5: Test remainder |
2982 | 1.92M | |
2983 | 1.92M | Note: Not needed because we always check that qhat is the correct value |
2984 | 1.92M | before performing the subtract. Value cast to mp_digit to prevent |
2985 | 1.92M | warning, qhat has been clamped to MP_DIGIT_MAX |
2986 | 1.92M | */ |
2987 | 1.92M | q.digits[j] = (mp_digit)qhat; |
2988 | 1.92M | |
2989 | 1.92M | /* D6: Add back |
2990 | 1.92M | Note: Not needed because we always check that qhat is the correct value |
2991 | 1.92M | before performing the subtract. |
2992 | 1.92M | */ |
2993 | 1.92M | |
2994 | 1.92M | /* D7: Loop on j */ |
2995 | 1.92M | r.digits--; |
2996 | 1.92M | ZERO(t.digits, t.alloc); |
2997 | 1.92M | } |
2998 | 842k | |
2999 | 842k | /* Get rid of leading zeros in q */ |
3000 | 842k | q.used = m + 1; |
3001 | 842k | CLAMP(&q); |
3002 | 842k | |
3003 | 842k | /* Denormalize the remainder */ |
3004 | 842k | CLAMP(u); /* use u here because the r.digits pointer is off-by-one */ |
3005 | 842k | if (k != 0) |
3006 | 827k | s_qdiv(u, k); |
3007 | 842k | |
3008 | 842k | mp_int_copy(u, v); /* ok: 0 <= r < v */ |
3009 | 842k | mp_int_copy(&q, u); /* ok: q <= u */ |
3010 | 842k | |
3011 | 842k | mp_int_clear(&t); |
3012 | 842k | CLEANUP: |
3013 | 842k | mp_int_clear(&q); |
3014 | 842k | return res; |
3015 | 1.71M | } |
3016 | | |
3017 | | STATIC int s_outlen(mp_int z, mp_size r) |
3018 | 1.54k | { |
3019 | 1.54k | mp_result bits; |
3020 | 1.54k | double raw; |
3021 | 1.54k | |
3022 | 1.54k | assert(r >= MP_MIN_RADIX && r <= MP_MAX_RADIX); |
3023 | 1.54k | |
3024 | 1.54k | bits = mp_int_count_bits(z); |
3025 | 1.54k | raw = (double)bits * s_log2[r]; |
3026 | 1.54k | |
3027 | 1.54k | return (int)(raw + 0.999999); |
3028 | 1.54k | } |
3029 | | |
3030 | | STATIC mp_size s_inlen(int len, mp_size r) |
3031 | 150 | { |
3032 | 150 | double raw = (double)len / s_log2[r]; |
3033 | 150 | mp_size bits = (mp_size)(raw + 0.5); |
3034 | 150 | |
3035 | 150 | return (mp_size)((bits + (MP_DIGIT_BIT150 - 1)) / MP_DIGIT_BIT150 ) + 1; |
3036 | 150 | } |
3037 | | |
3038 | | STATIC int s_ch2val(char c, int r) |
3039 | 1.83k | { |
3040 | 1.83k | int out; |
3041 | 1.83k | |
3042 | 1.83k | if (isdigit((unsigned char) c)) |
3043 | 1.83k | out = c - '0'; |
3044 | 0 | else if (0 r > 10 && 0 isalpha((unsigned char) c)0 ) |
3045 | 0 | out = toupper(c) - 'A' + 10; |
3046 | 0 | else |
3047 | 0 | return -1; |
3048 | 1.83k | |
3049 | 1.83k | return (out >= r) ? 1.83k -10 : out1.83k ; |
3050 | 1.83k | } |
3051 | | |
3052 | | STATIC char s_val2ch(int v, int caps) |
3053 | 26.0k | { |
3054 | 26.0k | assert(v >= 0); |
3055 | 26.0k | |
3056 | 26.0k | if (v < 10) |
3057 | 26.0k | return v + '0'; |
3058 | 0 | else { |
3059 | 0 | char out = (v - 10) + 'a'; |
3060 | 0 |
|
3061 | 0 | if (caps) |
3062 | 0 | return toupper(out); |
3063 | 0 | else |
3064 | 0 | return out; |
3065 | 0 | } |
3066 | 26.0k | } |
3067 | | |
3068 | | STATIC void s_2comp(unsigned char *buf, int len) |
3069 | 0 | { |
3070 | 0 | int i; |
3071 | 0 | unsigned short s = 1; |
3072 | 0 |
|
3073 | 0 | for (i = len - 1; i >= 00 ; --i0 ) { |
3074 | 0 | unsigned char c = ~buf[i]; |
3075 | 0 |
|
3076 | 0 | s = c + s; |
3077 | 0 | c = s & UCHAR_MAX; |
3078 | 0 | s >>= CHAR_BIT; |
3079 | 0 |
|
3080 | 0 | buf[i] = c; |
3081 | 0 | } |
3082 | 0 |
|
3083 | 0 | /* last carry out is ignored */ |
3084 | 0 | } |
3085 | | |
3086 | | STATIC mp_result s_tobin(mp_int z, unsigned char *buf, int *limpos, int pad) |
3087 | 0 | { |
3088 | 0 | mp_size uz; |
3089 | 0 | mp_digit *dz; |
3090 | 0 | int pos = 0, limit = *limpos; |
3091 | 0 |
|
3092 | 0 | uz = MP_USED0 (z); dz = MP_DIGITS0 (z); |
3093 | 0 | while (uz > 0 && 0 pos < limit0 ) { |
3094 | 0 | mp_digit d = *dz++; |
3095 | 0 | int i; |
3096 | 0 |
|
3097 | 0 | for (i = sizeof(mp_digit); i > 0 && 0 pos < limit0 ; --i0 ) { |
3098 | 0 | buf[pos++] = (unsigned char)d; |
3099 | 0 | d >>= CHAR_BIT; |
3100 | 0 |
|
3101 | 0 | /* Don't write leading zeroes */ |
3102 | 0 | if (d == 0 && 0 uz == 10 ) |
3103 | 0 | i = 0; /* exit loop without signaling truncation */ |
3104 | 0 | } |
3105 | 0 |
|
3106 | 0 | /* Detect truncation (loop exited with pos >= limit) */ |
3107 | 0 | if (i > 00 ) break0 ; |
3108 | 0 |
|
3109 | 0 | --uz; |
3110 | 0 | } |
3111 | 0 |
|
3112 | 0 | if (pad != 0 && 0 (buf[pos - 1] >> (CHAR_BIT - 1))0 ) { |
3113 | 0 | if (pos < limit) |
3114 | 0 | buf[pos++] = 0; |
3115 | 0 | else |
3116 | 0 | uz = 1; |
3117 | 0 | } |
3118 | 0 |
|
3119 | 0 | /* Digits are in reverse order, fix that */ |
3120 | 0 | REV(unsigned char, buf, pos); |
3121 | 0 |
|
3122 | 0 | /* Return the number of bytes actually written */ |
3123 | 0 | *limpos = pos; |
3124 | 0 |
|
3125 | 0 | return (uz == 0) ? MP_OK0 : MP_TRUNC0 ; |
3126 | 0 | } |
3127 | | |
3128 | | #if DEBUG |
3129 | | void s_print(char *tag, mp_int z) |
3130 | | { |
3131 | | int i; |
3132 | | |
3133 | | fprintf(stderr, "%s: %c ", tag, |
3134 | | (MP_SIGN(z) == MP_NEG) ? '-' : '+'); |
3135 | | |
3136 | | for (i = MP_USED(z) - 1; i >= 0; --i) |
3137 | | fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), z->digits[i]); |
3138 | | |
3139 | | fputc('\n', stderr); |
3140 | | |
3141 | | } |
3142 | | |
3143 | | void s_print_buf(char *tag, mp_digit *buf, mp_size num) |
3144 | | { |
3145 | | int i; |
3146 | | |
3147 | | fprintf(stderr, "%s: ", tag); |
3148 | | |
3149 | | for (i = num - 1; i >= 0; --i) |
3150 | | fprintf(stderr, "%0*X", (int)(MP_DIGIT_BIT / 4), buf[i]); |
3151 | | |
3152 | | fputc('\n', stderr); |
3153 | | } |
3154 | | #endif |
3155 | | |
3156 | | /* Here there be dragons */ |