1 /*
2 * Copyright (c) 2002 by The XFree86 Project, Inc.
3 *
4 * Permission is hereby granted, free of charge, to any person obtaining a
5 * copy of this software and associated documentation files (the "Software"),
6 * to deal in the Software without restriction, including without limitation
7 * the rights to use, copy, modify, merge, publish, distribute, sublicense,
8 * and/or sell copies of the Software, and to permit persons to whom the
9 * Software is furnished to do so, subject to the following conditions:
10 *
11 * The above copyright notice and this permission notice shall be included in
12 * all copies or substantial portions of the Software.
13 *
14 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
15 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
16 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
17 * THE XFREE86 PROJECT BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
18 * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF
19 * OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
20 * SOFTWARE.
21 *
22 * Except as contained in this notice, the name of the XFree86 Project shall
23 * not be used in advertising or otherwise to promote the sale, use or other
24 * dealings in this Software without prior written authorization from the
25 * XFree86 Project.
26 *
27 * Author: Paulo César Pereira de Andrade
28 */
29
30 /* $XFree86: xc/programs/xedit/lisp/mp/mpi.c,v 1.12 2002/11/20 07:44:43 paulo Exp $ */
31
32 #include "mp.h"
33
34 #ifdef __UNIXOS2__
35 # define finite(x) isfinite(x)
36 #endif
37
38 /*
39 * Prototypes
40 */
41 /* do the hard work of mpi_add and mpi_sub */
42 static void mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub);
43
44 /* logical functions implementation */
45 static INLINE BNS mpi_logic(BNS op1, BNS op2, BNS op);
46 static void mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op);
47
48 /* internal mpi_seti, whithout memory allocation */
49 static void _mpi_seti(mpi *rop, long si);
50
51 /*
52 * Initialization
53 */
54 static BNS onedig[1] = { 1 };
55 static mpi mpone = { 1, 1, 0, (BNS*)&onedig };
56
57 /*
58 * Implementation
59 */
60 void
mpi_init(mpi * op)61 mpi_init(mpi *op)
62 {
63 op->sign = 0;
64 op->size = op->alloc = 1;
65 op->digs = mp_malloc(sizeof(BNS));
66 op->digs[0] = 0;
67 }
68
69 void
mpi_clear(mpi * op)70 mpi_clear(mpi *op)
71 {
72 op->sign = 0;
73 op->size = op->alloc = 0;
74 mp_free(op->digs);
75 }
76
77 void
mpi_set(mpi * rop,mpi * op)78 mpi_set(mpi *rop, mpi *op)
79 {
80 if (rop != op) {
81 if (rop->alloc < op->size) {
82 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size);
83 rop->alloc = op->size;
84 }
85 rop->size = op->size;
86 memcpy(rop->digs, op->digs, sizeof(BNS) * op->size);
87 rop->sign = op->sign;
88 }
89 }
90
91 void
mpi_seti(mpi * rop,long si)92 mpi_seti(mpi *rop, long si)
93 {
94 unsigned long ui;
95 int sign = si < 0;
96 int size;
97
98 if (si == MINSLONG) {
99 ui = MINSLONG;
100 size = 2;
101 }
102 else {
103 if (sign)
104 ui = -si;
105 else
106 ui = si;
107 if (ui < CARRY)
108 size = 1;
109 else
110 size = 2;
111 }
112
113 if (rop->alloc < size) {
114 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
115 rop->alloc = size;
116 }
117 rop->size = size;
118
119 /* store data in small mp integer */
120 rop->digs[0] = (BNS)ui;
121 if (size > 1)
122 rop->digs[1] = (BNS)(ui >> BNSBITS);
123 rop->size = size;
124
125 /* adjust result sign */
126 rop->sign = sign;
127 }
128
129 static void
_mpi_seti(mpi * rop,long si)130 _mpi_seti(mpi *rop, long si)
131 {
132 unsigned long ui;
133 int sign = si < 0;
134 int size;
135
136 if (si == MINSLONG) {
137 ui = MINSLONG;
138 size = 2;
139 }
140 else {
141 if (sign)
142 ui = -si;
143 else
144 ui = si;
145 if (ui < CARRY)
146 size = 1;
147 else
148 size = 2;
149 }
150
151 rop->digs[0] = (BNS)ui;
152 if (size > 1)
153 rop->digs[1] = (BNS)(ui >> BNSBITS);
154 rop->size = size;
155
156 rop->sign = sign;
157 }
158
159 void
mpi_setd(mpi * rop,double d)160 mpi_setd(mpi *rop, double d)
161 {
162 long i;
163 double mantissa;
164 int shift, exponent;
165 BNI size;
166
167 if (isnan(d))
168 d = 0.0;
169 else if (!finite(d))
170 d = copysign(1.0, d) * DBL_MAX;
171
172 /* check if number is larger than 1 */
173 if (fabs(d) < 1.0) {
174 rop->digs[0] = 0;
175 rop->size = 1;
176 rop->sign = d < 0.0;
177
178 return;
179 }
180
181 mantissa = frexp(d, &exponent);
182 if (mantissa < 0)
183 mantissa = -mantissa;
184
185 size = (exponent + (BNSBITS - 1)) / BNSBITS;
186 shift = BNSBITS - (exponent & (BNSBITS - 1));
187
188 /* adjust amount of memory */
189 if (rop->alloc < size) {
190 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
191 rop->alloc = size;
192 }
193 rop->size = size;
194
195 /* adjust the exponent */
196 if (shift < BNSBITS)
197 mantissa = ldexp(mantissa, -shift);
198
199 /* convert double */
200 for (i = size - 1; i >= 0 && mantissa != 0.0; i--) {
201 mantissa = ldexp(mantissa, BNSBITS);
202 rop->digs[i] = (BNS)mantissa;
203 mantissa -= rop->digs[i];
204 }
205 for (; i >= 0; i--)
206 rop->digs[i] = 0;
207
208 /* normalize */
209 if (size > 1 && rop->digs[size - 1] == 0)
210 --rop->size;
211
212 rop->sign = d < 0.0;
213 }
214
215 /* how many BNS in the given base, log(base) / log(CARRY) */
216 #ifdef LONG64
217 static double str_bases[37] = {
218 0.0000000000000000, 0.0000000000000000, 0.0312500000000000,
219 0.0495300781475362, 0.0625000000000000, 0.0725602529652301,
220 0.0807800781475362, 0.0877298413143002, 0.0937500000000000,
221 0.0990601562950723, 0.1038102529652301, 0.1081072380824156,
222 0.1120300781475362, 0.1156387411919092, 0.1189798413143002,
223 0.1220903311127662, 0.1250000000000000, 0.1277332137890731,
224 0.1303101562950723, 0.1327477347951120, 0.1350602529652300,
225 0.1372599194618363, 0.1393572380824156, 0.1413613111267817,
226 0.1432800781475362, 0.1451205059304602, 0.1468887411919092,
227 0.1485902344426084, 0.1502298413143002, 0.1518119060977367,
228 0.1533403311127662, 0.1548186346995899, 0.1562500000000000,
229 0.1576373162299517, 0.1589832137890731, 0.1602900942795302,
230 0.1615601562950723,
231 };
232 #else
233 static double str_bases[37] = {
234 0.0000000000000000, 0.0000000000000000, 0.0625000000000000,
235 0.0990601562950723, 0.1250000000000000, 0.1451205059304602,
236 0.1615601562950723, 0.1754596826286003, 0.1875000000000000,
237 0.1981203125901446, 0.2076205059304602, 0.2162144761648311,
238 0.2240601562950723, 0.2312774823838183, 0.2379596826286003,
239 0.2441806622255325, 0.2500000000000000, 0.2554664275781462,
240 0.2606203125901445, 0.2654954695902241, 0.2701205059304602,
241 0.2745198389236725, 0.2787144761648311, 0.2827226222535633,
242 0.2865601562950723, 0.2902410118609203, 0.2937774823838183,
243 0.2971804688852168, 0.3004596826286003, 0.3036238121954733,
244 0.3066806622255324, 0.3096372693991797, 0.3125000000000000,
245 0.3152746324599034, 0.3179664275781462, 0.3205801885590604,
246 0.3231203125901446,
247 };
248 #endif
249
250 void
mpi_setstr(mpi * rop,char * str,int base)251 mpi_setstr(mpi *rop, char *str, int base)
252 {
253 long i; /* counter */
254 int sign; /* result sign */
255 BNI carry; /* carry value */
256 BNI value; /* temporary value */
257 BNI size; /* size of result */
258 char *ptr; /* end of valid input */
259
260 /* initialization */
261 sign = 0;
262 carry = 0;
263
264 /* skip leading spaces */
265 while (isspace(*str))
266 ++str;
267
268 /* check if sign supplied */
269 if (*str == '-') {
270 sign = 1;
271 ++str;
272 }
273 else if (*str == '+')
274 ++str;
275
276 /* skip leading zeros */
277 while (*str == '0')
278 ++str;
279
280 ptr = str;
281 while (*ptr) {
282 if (*ptr >= '0' && *ptr <= '9') {
283 if (*ptr - '0' >= base)
284 break;
285 }
286 else if (*ptr >= 'A' && *ptr <= 'Z') {
287 if (*ptr - 'A' + 10 >= base)
288 break;
289 }
290 else if (*ptr >= 'a' && *ptr <= 'z') {
291 if (*ptr - 'a' + 10 >= base)
292 break;
293 }
294 else
295 break;
296 ++ptr;
297 }
298
299 /* resulting size */
300 size = (ptr - str) * str_bases[base] + 1;
301
302 /* make sure rop has enough storage */
303 if (rop->alloc < size) {
304 rop->digs = mp_realloc(rop->digs, size * sizeof(BNS));
305 rop->alloc = size;
306 }
307 rop->size = size;
308
309 /* initialize rop to zero */
310 memset(rop->digs, '\0', size * sizeof(BNS));
311
312 /* set result sign */
313 rop->sign = sign;
314
315 /* convert string */
316 for (; str < ptr; str++) {
317 value = *str;
318 if (islower(value))
319 value = toupper(value);
320 value = value > '9' ? value - 'A' + 10 : value - '0';
321 value += (BNI)rop->digs[0] * base;
322 carry = value >> BNSBITS;
323 rop->digs[0] = (BNS)value;
324 for (i = 1; i < size; i++) {
325 value = (BNI)rop->digs[i] * base + carry;
326 carry = value >> BNSBITS;
327 rop->digs[i] = (BNS)value;
328 }
329 }
330
331 /* normalize */
332 if (rop->size > 1 && rop->digs[rop->size - 1] == 0)
333 --rop->size;
334 }
335
336 void
mpi_add(mpi * rop,mpi * op1,mpi * op2)337 mpi_add(mpi *rop, mpi *op1, mpi *op2)
338 {
339 mpi_addsub(rop, op1, op2, 0);
340 }
341
342 void
mpi_addi(mpi * rop,mpi * op1,long op2)343 mpi_addi(mpi *rop, mpi *op1, long op2)
344 {
345 BNS digs[2];
346 mpi op;
347
348 op.digs = (BNS*)digs;
349 _mpi_seti(&op, op2);
350
351 mpi_addsub(rop, op1, &op, 0);
352 }
353
354 void
mpi_sub(mpi * rop,mpi * op1,mpi * op2)355 mpi_sub(mpi *rop, mpi *op1, mpi *op2)
356 {
357 mpi_addsub(rop, op1, op2, 1);
358 }
359
360 void
mpi_subi(mpi * rop,mpi * op1,long op2)361 mpi_subi(mpi *rop, mpi *op1, long op2)
362 {
363 BNS digs[2];
364 mpi op;
365
366 op.digs = (BNS*)digs;
367 _mpi_seti(&op, op2);
368
369 mpi_addsub(rop, op1, &op, 1);
370 }
371
372 static void
mpi_addsub(mpi * rop,mpi * op1,mpi * op2,int sub)373 mpi_addsub(mpi *rop, mpi *op1, mpi *op2, int sub)
374 {
375 long xlen; /* maximum result size */
376
377 if (sub ^ (op1->sign == op2->sign)) {
378 /* plus one for possible carry */
379 xlen = MAX(op1->size, op2->size) + 1;
380 if (rop->alloc < xlen) {
381 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen);
382 rop->alloc = xlen;
383 }
384 rop->size = mp_add(rop->digs, op1->digs, op2->digs,
385 op1->size, op2->size);
386 rop->sign = op1->sign;
387 }
388 else {
389 long cmp; /* check for larger operator */
390
391 cmp = mpi_cmpabs(op1, op2);
392 if (cmp == 0) {
393 rop->digs[0] = 0;
394 rop->size = 1;
395 rop->sign = 0;
396 }
397 else {
398 xlen = MAX(op1->size, op2->size);
399 if (rop->alloc < xlen) {
400 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xlen);
401 rop->alloc = xlen;
402 }
403 if (cmp > 0) {
404 rop->size = mp_sub(rop->digs, op1->digs, op2->digs,
405 op1->size, op2->size);
406 rop->sign = op1->sign;
407 }
408 else {
409 rop->size = mp_sub(rop->digs, op2->digs, op1->digs,
410 op2->size, op1->size);
411 rop->sign = sub ^ op2->sign;
412 }
413 }
414 }
415 }
416
417 void
mpi_mul(mpi * rop,mpi * op1,mpi * op2)418 mpi_mul(mpi *rop, mpi *op1, mpi *op2)
419 {
420 int sign; /* sign flag */
421 BNS *digs; /* result data */
422 long xsize; /* result size */
423
424 /* get result sign */
425 sign = op1->sign ^ op2->sign;
426
427 /* check for special cases */
428 if (op1->size == 1) {
429 if (*op1->digs == 0) {
430 /* multiply by 0 */
431 mpi_seti(rop, 0);
432 return;
433 }
434 else if (*op1->digs == 1) {
435 /* multiply by +-1 */
436 if (rop->alloc < op2->size) {
437 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op2->size);
438 rop->alloc = op2->size;
439 }
440 rop->size = op2->size;
441 memmove(rop->digs, op2->digs, sizeof(BNS) * op2->size);
442 rop->sign = op2->size > 1 || *op2->digs ? sign : 0;
443
444 return;
445 }
446 }
447 else if (op2->size == 1) {
448 if (*op2->digs == 0) {
449 /* multiply by 0 */
450 mpi_seti(rop, 0);
451 return;
452 }
453 else if (*op2->digs == 1) {
454 /* multiply by +-1 */
455 if (rop->alloc < op1->size) {
456 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op1->size);
457 rop->alloc = op1->size;
458 }
459 rop->size = op1->size;
460 memmove(rop->digs, op1->digs, sizeof(BNS) * op1->size);
461 rop->sign = op1->size > 1 || *op1->digs ? sign : 0;
462
463 return;
464 }
465 }
466
467 /* allocate result data and set it to zero */
468 xsize = op1->size + op2->size;
469 if (rop->digs == op1->digs || rop->digs == op2->digs)
470 /* rop is also an operand */
471 digs = mp_calloc(1, sizeof(BNS) * xsize);
472 else {
473 if (rop->alloc < xsize) {
474 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize);
475 rop->alloc = xsize;
476 }
477 digs = rop->digs;
478 memset(digs, '\0', sizeof(BNS) * xsize);
479 }
480
481 /* multiply operands */
482 xsize = mp_mul(digs, op1->digs, op2->digs, op1->size, op2->size);
483
484 /* store result in rop */
485 if (digs != rop->digs) {
486 /* if rop was an operand, free old data */
487 mp_free(rop->digs);
488 rop->digs = digs;
489 }
490 rop->size = xsize;
491
492 /* set result sign */
493 rop->sign = sign;
494 }
495
496 void
mpi_muli(mpi * rop,mpi * op1,long op2)497 mpi_muli(mpi *rop, mpi *op1, long op2)
498 {
499 BNS digs[2];
500 mpi op;
501
502 op.digs = (BNS*)digs;
503 _mpi_seti(&op, op2);
504
505 mpi_mul(rop, op1, &op);
506 }
507
508 void
mpi_div(mpi * rop,mpi * num,mpi * den)509 mpi_div(mpi *rop, mpi *num, mpi *den)
510 {
511 mpi_divqr(rop, NULL, num, den);
512 }
513
514 void
mpi_rem(mpi * rop,mpi * num,mpi * den)515 mpi_rem(mpi *rop, mpi *num, mpi *den)
516 {
517 mpi_divqr(NULL, rop, num, den);
518 }
519
520 /*
521 * Could/should be changed to not allocate qdigs if qrop is NULL
522 * Performance wouldn't suffer too much with a test on every loop iteration.
523 */
524 void
mpi_divqr(mpi * qrop,mpi * rrop,mpi * num,mpi * den)525 mpi_divqr(mpi *qrop, mpi *rrop, mpi *num, mpi *den)
526 {
527 long i, j; /* counters */
528 int qsign; /* sign of quotient */
529 int rsign; /* sign of remainder */
530 BNI qsize; /* size of quotient */
531 BNI rsize; /* size of remainder */
532 BNS qest; /* estimative of quotient value */
533 BNS *qdigs, *rdigs; /* work copy or result */
534 BNS *ndigs, *ddigs; /* work copy or divisor and dividend */
535 BNI value; /* temporary result */
536 long svalue; /* signed temporary result (2's complement) */
537 BNS carry, scarry, denorm; /* carry and normalization */
538 BNI dpos, npos; /* offsets in data */
539
540 /* get signs */
541 rsign = num->sign;
542 qsign = rsign ^ den->sign;
543
544 /* check for special case */
545 if (num->size < den->size) {
546 /* quotient is zero and remainder is numerator */
547 if (rrop && rrop->digs != num->digs) {
548 if (rrop->alloc < num->size) {
549 rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * num->size);
550 rrop->alloc = num->size;
551 }
552 rrop->size = num->size;
553 memcpy(rrop->digs, num->digs, sizeof(BNS) * num->size);
554 rrop->sign = rsign;
555 }
556 if (qrop)
557 mpi_seti(qrop, 0);
558
559 return;
560 }
561
562 /* estimate result sizes */
563 rsize = den->size;
564 qsize = num->size - den->size + 1;
565
566 /* offsets */
567 npos = num->size - 1;
568 dpos = den->size - 1;
569
570 /* allocate space for quotient and remainder */
571 if (qrop == NULL || qrop->digs == num->digs || qrop->digs == den->digs)
572 qdigs = mp_calloc(1, sizeof(BNS) * qsize);
573 else {
574 if (qrop->alloc < qsize) {
575 qrop->digs = mp_realloc(qrop->digs, sizeof(BNS) * qsize);
576 qrop->alloc = qsize;
577 }
578 memset(qrop->digs, '\0', sizeof(BNS) * qsize);
579 qdigs = qrop->digs;
580 }
581 if (rrop) {
582 if (rrop->digs == num->digs || rrop->digs == den->digs)
583 rdigs = mp_calloc(1, sizeof(BNS) * rsize);
584 else {
585 if (rrop->alloc < rsize) {
586 rrop->digs = mp_realloc(rrop->digs, sizeof(BNS) * rsize);
587 rrop->alloc = rsize;
588 }
589 memset(rrop->digs, '\0', sizeof(BNS) * rsize);
590 rdigs = rrop->digs;
591 }
592 }
593 else
594 rdigs = NULL; /* fix gcc warning */
595
596 /* special case, only one word in divisor */
597 if (dpos == 0) {
598 for (carry = 0, i = npos; i >= 0; i--) {
599 value = ((BNI)carry << BNSBITS) + num->digs[i];
600 qdigs[i] = (BNS)(value / den->digs[0]);
601 carry = (BNS)(value % den->digs[0]);
602 }
603 if (rrop)
604 rdigs[0] = carry;
605
606 goto mpi_divqr_done;
607 }
608
609 /* make work copy of numerator */
610 ndigs = mp_malloc(sizeof(BNS) * (num->size + 1));
611 /* allocate one extra word an update offset */
612 memcpy(ndigs, num->digs, sizeof(BNS) * num->size);
613 ndigs[num->size] = 0;
614 ++npos;
615
616 /* normalize */
617 denorm = (BNS)((BNI)CARRY / ((BNI)(den->digs[dpos]) + 1));
618
619 if (denorm > 1) {
620 /* i <= num->size because ndigs has an extra word */
621 for (carry = 0, i = 0; i <= num->size; i++) {
622 value = ndigs[i] * (BNI)denorm + carry;
623 ndigs[i] = (BNS)value;
624 carry = (BNS)(value >> BNSBITS);
625 }
626 /* make work copy of denominator */
627 ddigs = mp_malloc(sizeof(BNS) * den->size);
628 memcpy(ddigs, den->digs, sizeof(BNS) * den->size);
629 for (carry = 0, i = 0; i < den->size; i++) {
630 value = ddigs[i] * (BNI)denorm + carry;
631 ddigs[i] = (BNS)value;
632 carry = (BNS)(value >> BNSBITS);
633 }
634 }
635 else
636 /* only allocate copy of denominator if going to change it */
637 ddigs = den->digs;
638
639 /* divide mp integers */
640 for (j = qsize - 1; j >= 0; j--, npos--) {
641 /* estimate quotient */
642 if (ndigs[npos] == ddigs[dpos])
643 qest = (BNS)SMASK;
644 else
645 qest = (BNS)((((BNI)(ndigs[npos]) << BNSBITS) + ndigs[npos - 1]) /
646 ddigs[dpos]);
647
648 while ((value = ((BNI)(ndigs[npos]) << BNSBITS) + ndigs[npos - 1] -
649 qest * (BNI)(ddigs[dpos])) < CARRY &&
650 ddigs[dpos - 1] * (BNI)qest >
651 (value << BNSBITS) + ndigs[npos - 2])
652 --qest;
653
654 /* multiply and subtract */
655 carry = scarry = 0;
656 for (i = 0; i < den->size; i++) {
657 value = qest * (BNI)ddigs[i] + carry;
658 carry = (BNS)(value >> BNSBITS);
659 svalue = (long)ndigs[npos - dpos + i - 1] - (long)(value & SMASK) -
660 (long)scarry;
661 ndigs[npos - dpos + i - 1] = (BNS)svalue;
662 scarry = svalue < 0;
663 }
664
665 svalue = (long)ndigs[npos] - (long)(carry & SMASK) - (long)scarry;
666 ndigs[npos] = (BNS)svalue;
667
668 if (svalue & LMASK) {
669 /* quotient too big */
670 --qest;
671 carry = 0;
672 for (i = 0; i < den->size; i++) {
673 value = ndigs[npos - dpos + i - 1] + (BNI)carry + (BNI)ddigs[i];
674 ndigs[npos - dpos + i - 1] = (BNS)value;
675 carry = (BNS)(value >> BNSBITS);
676 }
677 ndigs[npos] += carry;
678 }
679
680 qdigs[j] = qest;
681 }
682
683 /* calculate remainder */
684 if (rrop) {
685 for (carry = 0, j = dpos; j >= 0; j--) {
686 value = ((BNI)carry << BNSBITS) + ndigs[j];
687 rdigs[j] = (BNS)(value / denorm);
688 carry = (BNS)(value % denorm);
689 }
690 }
691
692 mp_free(ndigs);
693 if (ddigs != den->digs)
694 mp_free(ddigs);
695
696 mpi_divqr_done:
697 if (rrop) {
698 if (rrop->digs != rdigs)
699 mp_free(rrop->digs);
700 /* normalize remainder */
701 for (i = rsize - 1; i >= 0; i--)
702 if (rdigs[i] != 0)
703 break;
704 if (i != rsize - 1) {
705 if (i < 0) {
706 rsign = 0;
707 rsize = 1;
708 }
709 else
710 rsize = i + 1;
711 }
712 rrop->digs = rdigs;
713 rrop->sign = rsign;
714 rrop->size = rsize;
715 }
716
717 /* normalize quotient */
718 if (qrop) {
719 if (qrop->digs != qdigs)
720 mp_free(qrop->digs);
721 for (i = qsize - 1; i >= 0; i--)
722 if (qdigs[i] != 0)
723 break;
724 if (i != qsize - 1) {
725 if (i < 0) {
726 qsign = 0;
727 qsize = 1;
728 }
729 else
730 qsize = i + 1;
731 }
732 qrop->digs = qdigs;
733 qrop->sign = qsign;
734 qrop->size = qsize;
735 }
736 else
737 mp_free(qdigs);
738 }
739
740 long
mpi_divqri(mpi * qrop,mpi * num,long den)741 mpi_divqri(mpi *qrop, mpi *num, long den)
742 {
743 BNS ddigs[2];
744 mpi dop, rrop;
745 long remainder;
746
747 dop.digs = (BNS*)ddigs;
748 _mpi_seti(&dop, den);
749
750 memset(&rrop, '\0', sizeof(mpi));
751 mpi_init(&rrop);
752 mpi_divqr(qrop, &rrop, num, &dop);
753 remainder = rrop.digs[0];
754 if (rrop.size > 1)
755 remainder |= (BNI)(rrop.digs[1]) << BNSBITS;
756 if (rrop.sign)
757 remainder = -remainder;
758 mpi_clear(&rrop);
759
760 return (remainder);
761 }
762
763 void
mpi_divi(mpi * rop,mpi * num,long den)764 mpi_divi(mpi *rop, mpi *num, long den)
765 {
766 BNS ddigs[2];
767 mpi dop;
768
769 dop.digs = (BNS*)ddigs;
770 _mpi_seti(&dop, den);
771
772 mpi_divqr(rop, NULL, num, &dop);
773 }
774
775 long
mpi_remi(mpi * num,long den)776 mpi_remi(mpi *num, long den)
777 {
778 return (mpi_divqri(NULL, num, den));
779 }
780
781 void
mpi_mod(mpi * rop,mpi * num,mpi * den)782 mpi_mod(mpi *rop, mpi *num, mpi *den)
783 {
784 mpi_rem(rop, num, den);
785 if (num->sign ^ den->sign)
786 mpi_add(rop, rop, den);
787 }
788
789 long
mpi_modi(mpi * num,long den)790 mpi_modi(mpi *num, long den)
791 {
792 long remainder;
793
794 remainder = mpi_remi(num, den);
795 if (num->sign ^ (den < 0))
796 remainder += den;
797
798 return (remainder);
799 }
800
801 void
mpi_gcd(mpi * rop,mpi * num,mpi * den)802 mpi_gcd(mpi *rop, mpi *num, mpi *den)
803 {
804 long cmp;
805 mpi rem;
806
807 /* check if result already given */
808 cmp = mpi_cmpabs(num, den);
809
810 /* check if num is equal to den or if num is zero */
811 if (cmp == 0 || (num->size == 1 && num->digs[0] == 0)) {
812 mpi_set(rop, den);
813 rop->sign = 0;
814 return;
815 }
816 /* check if den is not zero */
817 if (den->size == 1 && den->digs[0] == 0) {
818 mpi_set(rop, num);
819 rop->sign = 0;
820 return;
821 }
822
823 /* don't call mpi_init, relies on realloc(0, size) == malloc(size) */
824 memset(&rem, '\0', sizeof(mpi));
825
826 /* if num larger than den */
827 if (cmp > 0) {
828 mpi_rem(&rem, num, den);
829 if (rem.size == 1 && rem.digs[0] == 0) {
830 /* exact division */
831 mpi_set(rop, den);
832 rop->sign = 0;
833 mpi_clear(&rem);
834 return;
835 }
836 mpi_set(rop, den);
837 }
838 else {
839 mpi_rem(&rem, den, num);
840 if (rem.size == 1 && rem.digs[0] == 0) {
841 /* exact division */
842 mpi_set(rop, num);
843 rop->sign = 0;
844 mpi_clear(&rem);
845 return;
846 }
847 mpi_set(rop, num);
848 }
849
850 /* loop using positive values */
851 rop->sign = rem.sign = 0;
852
853 /* cannot optimize this inverting rem/rop assignment earlier
854 * because rop mais be an operand */
855 mpi_swap(rop, &rem);
856
857 /* Euclides algorithm */
858 for (;;) {
859 mpi_rem(&rem, &rem, rop);
860 if (rem.size == 1 && rem.digs[0] == 0)
861 break;
862 mpi_swap(rop, &rem);
863 }
864 mpi_clear(&rem);
865 }
866
867 void
mpi_lcm(mpi * rop,mpi * num,mpi * den)868 mpi_lcm(mpi *rop, mpi *num, mpi *den)
869 {
870 mpi gcd;
871
872 /* check for zero operand */
873 if ((num->size == 1 && num->digs[0] == 0) ||
874 (den->size == 1 && den->digs[0] == 0)) {
875 rop->digs[0] = 0;
876 rop->sign = 0;
877 return;
878 }
879
880 /* don't call mpi_init, relies on realloc(0, size) == malloc(size) */
881 memset(&gcd, '\0', sizeof(mpi));
882
883 mpi_gcd(&gcd, num, den);
884 mpi_div(&gcd, den, &gcd);
885 mpi_mul(rop, &gcd, num);
886 rop->sign = 0;
887
888 mpi_clear(&gcd);
889 }
890
891 void
mpi_pow(mpi * rop,mpi * op,unsigned long exp)892 mpi_pow(mpi *rop, mpi *op, unsigned long exp)
893 {
894 mpi zop, top;
895
896 if (exp == 2) {
897 mpi_mul(rop, op, op);
898 return;
899 }
900 /* check for op**0 */
901 else if (exp == 0) {
902 rop->digs[0] = 1;
903 rop->size = 1;
904 rop->sign = 0;
905 return;
906 }
907 else if (exp == 1) {
908 mpi_set(rop, op);
909 return;
910 }
911 else if (op->size == 1) {
912 if (op->digs[0] == 0) {
913 mpi_seti(rop, 0);
914 return;
915 }
916 else if (op->digs[0] == 1) {
917 mpi_seti(rop, op->sign && (exp & 1) ? -1 : 1);
918 return;
919 }
920 }
921
922 memset(&zop, '\0', sizeof(mpi));
923 memset(&top, '\0', sizeof(mpi));
924 mpi_set(&zop, op);
925 mpi_set(&top, op);
926 for (--exp; exp; exp >>= 1) {
927 if (exp & 1)
928 mpi_mul(&zop, &top, &zop);
929 mpi_mul(&top, &top, &top);
930 }
931
932 mpi_clear(&top);
933 rop->sign = zop.sign;
934 mp_free(rop->digs);
935 rop->digs = zop.digs;
936 rop->size = zop.size;
937 }
938
939 /* Find integer root of given number using the iteration
940 * x{n+1} = ((K-1) * x{n} + N / x{n}^(K-1)) / K
941 */
942 int
mpi_root(mpi * rop,mpi * op,unsigned long nth)943 mpi_root(mpi *rop, mpi *op, unsigned long nth)
944 {
945 long bits, cmp;
946 int exact;
947 int sign;
948 mpi *r, t, temp, quot, old, rem;
949
950 sign = op->sign;
951
952 /* divide by zero op**1/0 */
953 if (nth == 0) {
954 int one = 1, zero = 0;
955 one = one / zero;
956 }
957 /* result is complex */
958 if (sign && !(nth & 1)) {
959 int one = 1, zero = 0;
960 one = one / zero;
961 }
962
963 /* special case op**1/1 = op */
964 if (nth == 1) {
965 mpi_set(rop, op);
966 return (1);
967 }
968
969 bits = mpi_getsize(op, 2) - 2;
970
971 if (bits < 0 || bits / nth == 0) {
972 /* integral root is surely less than 2 */
973 exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0);
974 mpi_seti(rop, sign ? -1 : op->digs[0] == 0 ? 0 : 1);
975
976 return (exact == 1);
977 }
978
979 /* initialize */
980 if (rop != op)
981 r = rop;
982 else {
983 r = &t;
984 memset(r, '\0', sizeof(mpi));
985 }
986 memset(&temp, '\0', sizeof(mpi));
987 memset(", '\0', sizeof(mpi));
988 memset(&old, '\0', sizeof(mpi));
989 memset(&rem, '\0', sizeof(mpi));
990
991 if (sign)
992 r->sign = 0;
993
994 /* root aproximation */
995 mpi_ash(r, op, -(bits - (bits / nth)));
996
997 for (;;) {
998 mpi_pow(&temp, r, nth - 1);
999 mpi_divqr(", &rem, op, &temp);
1000 cmp = mpi_cmpabs(r, ");
1001 if (cmp == 0) {
1002 exact = mpi_cmpi(&rem, 0) == 0;
1003 break;
1004 }
1005 else if (cmp < 0) {
1006 if (mpi_cmpabs(r, &old) == 0) {
1007 exact = 0;
1008 break;
1009 }
1010 mpi_set(&old, r);
1011 }
1012 mpi_muli(&temp, r, nth - 1);
1013 mpi_add(", ", &temp);
1014 mpi_divi(r, ", nth);
1015 }
1016
1017 mpi_clear(&temp);
1018 mpi_clear(");
1019 mpi_clear(&old);
1020 mpi_clear(&rem);
1021 if (r != rop) {
1022 mpi_set(rop, r);
1023 mpi_clear(r);
1024 }
1025 rop->sign = sign;
1026
1027 return (exact);
1028 }
1029
1030 /*
1031 * Find square root using the iteration:
1032 * x{n+1} = (x{n}+N/x{n})/2
1033 */
1034 int
mpi_sqrt(mpi * rop,mpi * op)1035 mpi_sqrt(mpi *rop, mpi *op)
1036 {
1037 long bits, cmp;
1038 int exact;
1039 mpi *r, t, quot, rem, old;
1040
1041 /* result is complex */
1042 if (op->sign) {
1043 int one = 1, zero = 0;
1044 one = one / zero;
1045 }
1046
1047 bits = mpi_getsize(op, 2) - 2;
1048
1049 if (bits < 2) {
1050 /* integral root is surely less than 2 */
1051 exact = op->size == 1 && (op->digs[0] == 1 || op->digs[0] == 0);
1052 mpi_seti(rop, op->digs[0] == 0 ? 0 : 1);
1053
1054 return (exact == 1);
1055 }
1056
1057 /* initialize */
1058 if (rop != op)
1059 r = rop;
1060 else {
1061 r = &t;
1062 memset(r, '\0', sizeof(mpi));
1063 }
1064 memset(", '\0', sizeof(mpi));
1065 memset(&rem, '\0', sizeof(mpi));
1066 memset(&old, '\0', sizeof(mpi));
1067
1068 /* root aproximation */
1069 mpi_ash(r, op, -(bits - (bits / 2)));
1070
1071 for (;;) {
1072 if (mpi_cmpabs(r, &old) == 0) {
1073 exact = 0;
1074 break;
1075 }
1076 mpi_divqr(", &rem, op, r);
1077 cmp = mpi_cmpabs(", r);
1078 if (cmp == 0) {
1079 exact = mpi_cmpi(&rem, 0) == 0;
1080 break;
1081 }
1082 else if (cmp > 0 && rem.size == 1 && rem.digs[0] == 0)
1083 mpi_subi(", ", 1);
1084 mpi_set(&old, r);
1085 mpi_add(r, r, ");
1086 mpi_ash(r, r, -1);
1087 }
1088 mpi_clear(");
1089 mpi_clear(&rem);
1090 mpi_clear(&old);
1091 if (r != rop) {
1092 mpi_set(rop, r);
1093 mpi_clear(r);
1094 }
1095
1096 return (exact);
1097 }
1098
1099 void
mpi_ash(mpi * rop,mpi * op,long shift)1100 mpi_ash(mpi *rop, mpi *op, long shift)
1101 {
1102 long i; /* counter */
1103 long xsize; /* maximum result size */
1104 BNS *digs;
1105
1106 /* check for 0 shift, multiply/divide by 1 */
1107 if (shift == 0) {
1108 if (rop != op) {
1109 if (rop->alloc < op->size) {
1110 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * op->size);
1111 rop->alloc = op->size;
1112 }
1113 rop->size = op->size;
1114 memcpy(rop->digs, op->digs, sizeof(BNS) * op->size);
1115 }
1116
1117 return;
1118 }
1119 else if (op->size == 1 && op->digs[0] == 0) {
1120 rop->sign = 0;
1121 rop->size = 1;
1122 rop->digs[0] = 0;
1123
1124 return;
1125 }
1126
1127 /* check shift and initialize */
1128 if (shift > 0)
1129 xsize = op->size + (shift / BNSBITS) + 1;
1130 else {
1131 xsize = op->size - ((-shift) / BNSBITS);
1132 if (xsize <= 0) {
1133 rop->size = 1;
1134 rop->sign = op->sign;
1135 rop->digs[0] = op->sign ? 1 : 0;
1136
1137 return;
1138 }
1139 }
1140
1141 /* allocate/adjust memory for result */
1142 if (rop == op)
1143 digs = mp_malloc(sizeof(BNS) * xsize);
1144 else {
1145 if (rop->alloc < xsize) {
1146 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * xsize);
1147 rop->alloc = xsize;
1148 }
1149 digs = rop->digs;
1150 }
1151
1152 /* left shift, multiply by power of two */
1153 if (shift > 0)
1154 rop->size = mp_lshift(digs, op->digs, op->size, shift);
1155 /* right shift, divide by power of two */
1156 else {
1157 long carry = 0;
1158
1159 if (op->sign) {
1160 BNI words, bits;
1161
1162 words = -shift / BNSBITS;
1163 bits = -shift % BNSBITS;
1164 for (i = 0; i < words; i++)
1165 carry |= op->digs[xsize + i];
1166 if (!carry) {
1167 for (i = 0; i < bits; i++)
1168 if (op->digs[op->size - xsize] & (1 << i)) {
1169 carry = 1;
1170 break;
1171 }
1172 }
1173 }
1174 rop->size = mp_rshift(digs, op->digs, op->size, -shift);
1175
1176 if (carry)
1177 /* emulates two's complement subtracting 1 from the result */
1178 rop->size = mp_add(digs, digs, mpone.digs, rop->size, 1);
1179 }
1180
1181 if (rop->digs != digs) {
1182 mp_free(rop->digs);
1183 rop->alloc = rop->size;
1184 rop->digs = digs;
1185 }
1186 rop->sign = op->sign;
1187 }
1188
1189 static INLINE BNS
mpi_logic(BNS op1,BNS op2,BNS op)1190 mpi_logic(BNS op1, BNS op2, BNS op)
1191 {
1192 switch (op) {
1193 case '&':
1194 return (op1 & op2);
1195 case '|':
1196 return (op1 | op2);
1197 case '^':
1198 return (op1 ^ op2);
1199 }
1200
1201 return (SMASK);
1202 }
1203
1204 static void
mpi_log(mpi * rop,mpi * op1,mpi * op2,BNS op)1205 mpi_log(mpi *rop, mpi *op1, mpi *op2, BNS op)
1206 {
1207 long i; /* counter */
1208 long c, c1, c2; /* carry */
1209 BNS *digs, *digs1, *digs2; /* pointers to mp data */
1210 BNI size, size1, size2;
1211 BNS sign, sign1, sign2;
1212 BNS n, n1, n2; /* logical operands */
1213 BNI sum;
1214
1215 /* initialize */
1216 size1 = op1->size;
1217 size2 = op2->size;
1218
1219 sign1 = op1->sign ? SMASK : 0;
1220 sign2 = op2->sign ? SMASK : 0;
1221
1222 sign = mpi_logic(sign1, sign2, op);
1223
1224 size = MAX(size1, size2);
1225 if (sign)
1226 ++size;
1227 if (rop->alloc < size) {
1228 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * size);
1229 rop->alloc = size;
1230 }
1231
1232 digs = rop->digs;
1233 digs1 = op1->digs;
1234 digs2 = op2->digs;
1235
1236 c = c1 = c2 = 1;
1237
1238 /* apply logical operation */
1239 for (i = 0; i < size; i++) {
1240 if (i >= size1)
1241 n1 = sign1;
1242 else if (sign1) {
1243 sum = (BNI)(BNS)(~digs1[i]) + c1;
1244 c1 = (long)(sum >> BNSBITS);
1245 n1 = (BNS)sum;
1246 }
1247 else
1248 n1 = digs1[i];
1249
1250 if (i >= size2)
1251 n2 = sign2;
1252 else if (sign2) {
1253 sum = (BNI)(BNS)(~digs2[i]) + c2;
1254 c2 = (long)(sum >> BNSBITS);
1255 n2 = (BNS)sum;
1256 }
1257 else
1258 n2 = digs2[i];
1259
1260 n = mpi_logic(n1, n2, op);
1261 if (sign) {
1262 sum = (BNI)(BNS)(~n) + c;
1263 c = (long)(sum >> BNSBITS);
1264 digs[i] = (BNS)sum;
1265 }
1266 else
1267 digs[i] = n;
1268 }
1269
1270 /* normalize */
1271 for (i = size - 1; i >= 0; i--)
1272 if (digs[i] != 0)
1273 break;
1274 if (i != size - 1) {
1275 if (i < 0) {
1276 sign = 0;
1277 size = 1;
1278 }
1279 else
1280 size = i + 1;
1281 }
1282
1283 rop->sign = sign != 0;
1284 rop->size = size;
1285 }
1286
1287 void
mpi_and(mpi * rop,mpi * op1,mpi * op2)1288 mpi_and(mpi *rop, mpi *op1, mpi *op2)
1289 {
1290 mpi_log(rop, op1, op2, '&');
1291 }
1292
1293 void
mpi_ior(mpi * rop,mpi * op1,mpi * op2)1294 mpi_ior(mpi *rop, mpi *op1, mpi *op2)
1295 {
1296 mpi_log(rop, op1, op2, '|');
1297 }
1298
1299 void
mpi_xor(mpi * rop,mpi * op1,mpi * op2)1300 mpi_xor(mpi *rop, mpi *op1, mpi *op2)
1301 {
1302 mpi_log(rop, op1, op2, '^');
1303 }
1304
1305 void
mpi_com(mpi * rop,mpi * op)1306 mpi_com(mpi *rop, mpi *op)
1307 {
1308 static BNS digs[1] = { 1 };
1309 static mpi one = { 0, 1, 1, (BNS*)&digs };
1310
1311 mpi_log(rop, rop, &one, '^');
1312 }
1313
1314 void
mpi_neg(mpi * rop,mpi * op)1315 mpi_neg(mpi *rop, mpi *op)
1316 {
1317 int sign = op->sign ^ 1;
1318
1319 if (rop->digs != op->digs) {
1320 if (rop->alloc < op->size) {
1321 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size);
1322 rop->alloc = op->size;
1323 }
1324 rop->size = op->size;
1325 memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size);
1326 }
1327
1328 rop->sign = sign;
1329 }
1330
1331 void
mpi_abs(mpi * rop,mpi * op)1332 mpi_abs(mpi *rop, mpi *op)
1333 {
1334 if (rop->digs != op->digs) {
1335 if (rop->alloc < op->size) {
1336 rop->digs = mp_realloc(rop->digs, sizeof(BNS) * rop->size);
1337 rop->alloc = op->size;
1338 }
1339 rop->size = op->size;
1340 memcpy(rop->digs, op->digs, sizeof(BNS) * rop->size);
1341 }
1342
1343 rop->sign = 0;
1344 }
1345
1346 int
mpi_cmp(mpi * op1,mpi * op2)1347 mpi_cmp(mpi *op1, mpi *op2)
1348 {
1349 if (op1->sign ^ op2->sign)
1350 return (op1->sign ? -1 : 1);
1351
1352 if (op1->size == op2->size) {
1353 long i, cmp = 0;
1354
1355 for (i = op1->size - 1; i >= 0; i--)
1356 if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0)
1357 break;
1358
1359 return (cmp == 0 ? 0 : (cmp < 0) ^ op1->sign ? -1 : 1);
1360 }
1361
1362 return ((op1->size < op2->size) ^ op1->sign ? -1 : 1);
1363 }
1364
1365 int
mpi_cmpi(mpi * op1,long op2)1366 mpi_cmpi(mpi *op1, long op2)
1367 {
1368 long cmp;
1369
1370 if (op1->size > 2)
1371 return (op1->sign ? -1 : 1);
1372
1373 cmp = op1->digs[0];
1374 if (op1->size == 2) {
1375 cmp |= (long)op1->digs[1] << BNSBITS;
1376 if (cmp == MINSLONG)
1377 return (op2 == cmp && op1->sign ? 0 : op1->sign ? -1 : 1);
1378 }
1379 if (op1->sign)
1380 cmp = -cmp;
1381
1382 return (cmp - op2);
1383 }
1384
1385 int
mpi_cmpabs(mpi * op1,mpi * op2)1386 mpi_cmpabs(mpi *op1, mpi *op2)
1387 {
1388 if (op1->size == op2->size) {
1389 long i, cmp = 0;
1390
1391 for (i = op1->size - 1; i >= 0; i--)
1392 if ((cmp = (long)op1->digs[i] - (long)op2->digs[i]) != 0)
1393 break;
1394
1395 return (cmp);
1396 }
1397
1398 return ((op1->size < op2->size) ? -1 : 1);
1399 }
1400
1401 int
mpi_cmpabsi(mpi * op1,long op2)1402 mpi_cmpabsi(mpi *op1, long op2)
1403 {
1404 unsigned long cmp;
1405
1406 if (op1->size > 2)
1407 return (1);
1408
1409 cmp = op1->digs[0];
1410 if (op1->size == 2)
1411 cmp |= (unsigned long)op1->digs[1] << BNSBITS;
1412
1413 return (cmp > op2 ? 1 : cmp == op2 ? 0 : -1);
1414 }
1415
1416 int
mpi_sgn(mpi * op)1417 mpi_sgn(mpi *op)
1418 {
1419 return (op->sign ? -1 : op->size > 1 || op->digs[0] ? 1 : 0);
1420 }
1421
1422 void
mpi_swap(mpi * op1,mpi * op2)1423 mpi_swap(mpi *op1, mpi *op2)
1424 {
1425 if (op1 != op2) {
1426 mpi swap;
1427
1428 memcpy(&swap, op1, sizeof(mpi));
1429 memcpy(op1, op2, sizeof(mpi));
1430 memcpy(op2, &swap, sizeof(mpi));
1431 }
1432 }
1433
1434 int
mpi_fiti(mpi * op)1435 mpi_fiti(mpi *op)
1436 {
1437 if (op->size == 1)
1438 return (1);
1439 else if (op->size == 2) {
1440 unsigned long value = ((BNI)(op->digs[1]) << BNSBITS) | op->digs[0];
1441
1442 if (value & MINSLONG)
1443 return (op->sign && value == MINSLONG) ? 1 : 0;
1444
1445 return (1);
1446 }
1447
1448 return (0);
1449 }
1450
1451 long
mpi_geti(mpi * op)1452 mpi_geti(mpi *op)
1453 {
1454 long value;
1455
1456 value = op->digs[0];
1457 if (op->size > 1)
1458 value |= (BNI)(op->digs[1]) << BNSBITS;
1459
1460 return (op->sign && value != MINSLONG ? -value : value);
1461 }
1462
1463 double
mpi_getd(mpi * op)1464 mpi_getd(mpi *op)
1465 {
1466 long i, len;
1467 double d = 0.0;
1468 int exponent;
1469
1470 #define FLOATDIGS sizeof(double) / sizeof(BNS)
1471
1472 switch (op->size) {
1473 case 2:
1474 d = (BNI)(op->digs[1]) << BNSBITS;
1475 case 1:
1476 d += op->digs[0];
1477 return (op->sign ? -d : d);
1478 default:
1479 break;
1480 }
1481
1482 for (i = 0, len = op->size; len > 0 && i < FLOATDIGS; i++)
1483 d = ldexp(d, BNSBITS) + op->digs[--len];
1484 d = frexp(d, &exponent);
1485 if (len > 0)
1486 exponent += len * BNSBITS;
1487
1488 if (d == 0.0)
1489 return (d);
1490
1491 d = ldexp(d, exponent);
1492
1493 return (op->sign ? -d : d);
1494 }
1495
1496 /* how many digits in a given base, floor(log(CARRY) / log(base)) */
1497 #ifdef LONG64
1498 static char dig_bases[37] = {
1499 0, 0, 32, 20, 16, 13, 12, 11, 10, 10, 9, 9, 8, 8, 8, 8,
1500 8, 7, 7, 7, 7, 7, 7, 7, 6, 6, 6, 6, 6, 6, 6, 6,
1501 6, 6, 6, 6, 6,
1502 };
1503 #else
1504 static char dig_bases[37] = {
1505 0, 0, 16, 10, 8, 6, 6, 5, 5, 5, 4, 4, 4, 4, 4, 4,
1506 4, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
1507 3, 3, 3, 3, 3,
1508 };
1509 #endif
1510
1511 /* how many digits per bit in a given base, log(2) / log(base) */
1512 static double bit_bases[37] = {
1513 0.0000000000000000, 0.0000000000000000, 1.0000000000000000,
1514 0.6309297535714575, 0.5000000000000000, 0.4306765580733931,
1515 0.3868528072345416, 0.3562071871080222, 0.3333333333333334,
1516 0.3154648767857287, 0.3010299956639811, 0.2890648263178878,
1517 0.2789429456511298, 0.2702381544273197, 0.2626495350371936,
1518 0.2559580248098155, 0.2500000000000000, 0.2446505421182260,
1519 0.2398124665681315, 0.2354089133666382, 0.2313782131597592,
1520 0.2276702486969530, 0.2242438242175754, 0.2210647294575037,
1521 0.2181042919855316, 0.2153382790366965, 0.2127460535533632,
1522 0.2103099178571525, 0.2080145976765095, 0.2058468324604344,
1523 0.2037950470905062, 0.2018490865820999, 0.2000000000000000,
1524 0.1982398631705605, 0.1965616322328226, 0.1949590218937863,
1525 0.1934264036172708,
1526 };
1527
1528 /* normalization base for string conversion, pow(base, dig_bases[base]) & ~CARRY */
1529 #ifdef LONG64
1530 static BNS big_bases[37] = {
1531 0x00000001, 0x00000001, 0x00000000, 0xCFD41B91, 0x00000000, 0x48C27395,
1532 0x81BF1000, 0x75DB9C97, 0x40000000, 0xCFD41B91, 0x3B9ACA00, 0x8C8B6D2B,
1533 0x19A10000, 0x309F1021, 0x57F6C100, 0x98C29B81, 0x00000000, 0x18754571,
1534 0x247DBC80, 0x3547667B, 0x4C4B4000, 0x6B5A6E1D, 0x94ACE180, 0xCAF18367,
1535 0x0B640000, 0x0E8D4A51, 0x1269AE40, 0x17179149, 0x1CB91000, 0x23744899,
1536 0x2B73A840, 0x34E63B41, 0x40000000, 0x4CFA3CC1, 0x5C13D840, 0x6D91B519,
1537 0x81BF1000,
1538 };
1539 #else
1540 static BNS big_bases[37] = {
1541 0x0001, 0x0001, 0x0000, 0xE6A9, 0x0000, 0x3D09, 0xB640, 0x41A7, 0x8000,
1542 0xE6A9, 0x2710, 0x3931, 0x5100, 0x6F91, 0x9610, 0xC5C1, 0x0000, 0x1331,
1543 0x16C8, 0x1ACB, 0x1F40, 0x242D, 0x2998, 0x2F87, 0x3600, 0x3D09, 0x44A8,
1544 0x4CE3, 0x55C0, 0x5F45, 0x6978, 0x745F, 0x8000, 0x8C61, 0x9988, 0xA77B,
1545 0xb640,
1546 };
1547 #endif
1548
1549 unsigned long
mpi_getsize(mpi * op,int base)1550 mpi_getsize(mpi *op, int base)
1551 {
1552 unsigned long value, bits;
1553
1554 value = op->digs[op->size - 1];
1555
1556 /* count leading bits */
1557 if (value) {
1558 unsigned long count, carry;
1559
1560 for (count = 0, carry = CARRY >> 1; carry; count++, carry >>= 1)
1561 if (value & carry)
1562 break;
1563
1564 bits = BNSBITS - count;
1565 }
1566 else
1567 bits = 0;
1568
1569 return ((bits + (op->size - 1) * BNSBITS) * bit_bases[base] + 1);
1570 }
1571
1572 char *
mpi_getstr(char * str,mpi * op,int base)1573 mpi_getstr(char *str, mpi *op, int base)
1574 {
1575 long i; /* counter */
1576 BNS *digs, *xdigs; /* copy of op data */
1577 BNI size; /* size of op */
1578 BNI digits; /* digits per word in given base */
1579 BNI bigbase; /* big base of given base */
1580 BNI strsize; /* size of resulting string */
1581 char *cp; /* pointer in str for conversion */
1582
1583 /* initialize */
1584 size = op->size;
1585 strsize = mpi_getsize(op, base) + op->sign + 1;
1586
1587 if (str == NULL)
1588 str = mp_malloc(strsize);
1589
1590 /* check for zero */
1591 if (size == 1 && op->digs[0] == 0) {
1592 str[0] = '0';
1593 str[1] = '\0';
1594
1595 return (str);
1596 }
1597
1598 digits = dig_bases[base];
1599 bigbase = big_bases[base];
1600
1601 cp = str + strsize;
1602 *--cp = '\0';
1603
1604 /* make copy of op data and adjust digs */
1605 xdigs = mp_malloc(size * sizeof(BNS));
1606 memcpy(xdigs, op->digs, size * sizeof(BNS));
1607 digs = xdigs + size - 1;
1608
1609 /* convert to string */
1610 for (;;) {
1611 long count = -1;
1612 BNI value;
1613 BNS quotient, remainder = 0;
1614
1615 /* if power of two base */
1616 if ((base & (base - 1)) == 0) {
1617 for (i = 0; i < size; i++) {
1618 quotient = remainder;
1619 remainder = digs[-i];
1620 digs[-i] = quotient;
1621 if (count < 0 && quotient)
1622 count = i;
1623 }
1624 }
1625 else {
1626 for (i = 0; i < size; i++) {
1627 value = digs[-i] + ((BNI)remainder << BNSBITS);
1628 quotient = (BNS)(value / bigbase);
1629 remainder = (BNS)(value % bigbase);
1630 digs[-i] = quotient;
1631 if (count < 0 && quotient)
1632 count = i;
1633 }
1634 }
1635 quotient = remainder;
1636 for (i = 0; i < digits; i++) {
1637 if (quotient == 0 && count < 0)
1638 break;
1639 remainder = quotient % base;
1640 quotient /= base;
1641 *--cp = remainder < 10 ? remainder + '0' : remainder - 10 + 'A';
1642 }
1643 if (count < 0)
1644 break;
1645 digs -= count;
1646 size -= count;
1647 }
1648
1649 /* adjust sign */
1650 if (op->sign)
1651 *--cp = '-';
1652
1653 /* remove any extra characters */
1654 if (cp > str)
1655 strcpy(str, cp);
1656
1657 mp_free(xdigs);
1658
1659 return (str);
1660 }
1661