1 /*
2 ** huge.c
3 **
4 ** Arbitrary precision integer library from Python sources.
5 **
6 ** This is a minor modification of the file "huge-number.c" taken from
7 ** mirrordir-0.10.49 which in turn contains these copyrights ...
8 **
9 ** $Id: huge.c,v 1.1.1.1 2001/04/12 18:08:01 ndwinton Exp $
10 */
11
12 /* huge-number.c: arbitrary precision integer library from Python sources
13 This has nothing to do with cryptography.
14 Copyright (C) 1998 Paul Sheer
15
16 This program is free software; you can redistribute it and/or modify
17 it under the terms of the GNU General Public License as published by
18 the Free Software Foundation; either version 2 of the License, or
19 (at your option) any later version.
20
21 This program is distributed in the hope that it will be useful,
22 but WITHOUT ANY WARRANTY; without even the implied warranty of
23 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 GNU General Public License for more details.
25
26 You should have received a copy of the GNU General Public License
27 along with this program; if not, write to the Free Software
28 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
29 */
30
31 /* This file was taken from the Python source for `long' type
32 integers. I have changed it to compile independently of the
33 Python source, and added the optimisation that GNU C can
34 use 31 bit digits instead of Python's 15 bit. You can download
35 the original from www.python.org. This file bears little
36 resemblance to the original though - paul */
37
38 /***********************************************************
39 Copyright 1991-1995 by Stichting Mathematisch Centrum, Amsterdam,
40 The Netherlands.
41
42 All Rights Reserved
43
44 Permission to use, copy, modify, and distribute this software and its
45 documentation for any purpose and without fee is hereby granted,
46 provided that the above copyright notice appear in all copies and that
47 both that copyright notice and this permission notice appear in
48 supporting documentation, and that the names of Stichting Mathematisch
49 Centrum or CWI or Corporation for National Research Initiatives or
50 CNRI not be used in advertising or publicity pertaining to
51 distribution of the software without specific, written prior
52 permission.
53
54 While CWI is the initial source for this software, a modified version
55 is made available by the Corporation for National Research Initiatives
56 (CNRI) at the Internet address ftp://ftp.python.org.
57
58 STICHTING MATHEMATISCH CENTRUM AND CNRI DISCLAIM ALL WARRANTIES WITH
59 REGARD TO THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF
60 MERCHANTABILITY AND FITNESS, IN NO EVENT SHALL STICHTING MATHEMATISCH
61 CENTRUM OR CNRI BE LIABLE FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL
62 DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR
63 PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER
64 TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
65 PERFORMANCE OF THIS SOFTWARE.
66
67 ******************************************************************/
68
69 #include <stdlib.h>
70 #include <string.h>
71 #include <stdio.h>
72 #include "huge.h"
73
74 #undef ABS
75 #undef MAX
76 #undef MIN
77 #define ABS(x) ((x) >= 0 ? (x) : -(x))
78 #define MAX(x, y) ((x) < (y) ? (y) : (x))
79 #define MIN(x, y) ((x) > (y) ? (y) : (x))
80
81 #define ob_size size
82 #define ob_digit d
83
84 #ifdef __GNUC__
85 #define huge_assert(x) { if (!(x)) { fprintf (stderr, "huge: assertion failed, %s:%d\n", __FILE__, __LINE__); abort(); } }
86 #else
87 #define huge_assert(x) { if (!(x)) abort(); }
88 #endif
89
90 static Huge *huge_normalize (Huge *);
91 static Huge *mul1 (Huge *, wdigit);
92 static Huge *muladd1 (Huge *, wdigit, wdigit);
93 static Huge *divrem1 (Huge *, wdigit, digit *);
94 static Huge *x_divrem (Huge * v1, Huge * w1, Huge ** prem);
95
96 #define huge_error(x) fprintf (stderr, "huge_%s\n", x)
97
98 /* Normalize (remove leading zeros from) a long int object.
99 Doesn't attempt to free the storage--in most cases, due to the nature
100 of the algorithms used, this could save at most be one word anyway. */
101
huge_normalize(Huge * v)102 static Huge *huge_normalize (Huge * v)
103 {
104 int j = ABS (v->ob_size);
105 int i = j;
106
107 while (i > 0 && v->ob_digit[i - 1] == 0)
108 --i;
109 if (i != j)
110 v->ob_size = (v->ob_size < 0) ? -(i) : i;
111 return v;
112 }
113
huge_new(int size)114 Huge *huge_new (int size)
115 {
116 Huge *h;
117 char *x;
118 h = malloc (sizeof (Huge) + size * sizeof (digit));
119 x = (char *) h;
120 x += sizeof (Huge);
121 h->d = (digit *) x;
122 h->size = size;
123 memset (h->d, 0, size * sizeof (digit));
124 return h;
125 }
126
huge_copy(Huge * a,Huge * b)127 void huge_copy (Huge * a, Huge * b)
128 {
129 int i;
130 for (i = 0; i < ABS (b->size); i++)
131 a->d[i] = b->d[i];
132 a->size = b->size;
133 }
134
huge_dup(Huge * a)135 Huge *huge_dup (Huge * a)
136 {
137 Huge *b;
138 if (!a)
139 return 0;
140 b = huge_new (ABS (a->ob_size));
141 huge_copy (b, a);
142 return b;
143 }
144
145 /* Create a new long int object from a C long int */
146
huge_from_long(long ival)147 Huge *huge_from_long (long ival)
148 {
149 /* Assume a C long fits in at most 5 'digits' */
150 /* Works on both 32- and 64-bit machines */
151 Huge *v = huge_new (5);
152 unsigned long t = ival;
153 int i;
154 if (ival < 0) {
155 t = -ival;
156 v->ob_size = -(v->ob_size);
157 }
158 for (i = 0; i < 5; i++) {
159 v->ob_digit[i] = (digit) (t & MASK);
160 t >>= SHIFT;
161 }
162 return huge_normalize (v);
163 }
164
huge_set_bit(Huge * v,int i)165 Huge *huge_set_bit (Huge * v, int i)
166 {
167 Huge *w;
168 w = huge_new (MAX (ABS (v->ob_size), i / SHIFT + 1));
169 huge_copy (w, v);
170 w->d[i / SHIFT] |= (1 << (i % SHIFT));
171 return w;
172 }
173
huge_clear_bit(Huge * v,int i)174 void huge_clear_bit (Huge * v, int i)
175 {
176 if (i / SHIFT < ABS (v->ob_size))
177 v->d[i / SHIFT] &= ~(1 << (i % SHIFT));
178 huge_normalize (v);
179 }
180
_huge_get_char(Huge * a,int j)181 static inline unsigned char _huge_get_char (Huge * a, int j)
182 {
183 twodigits r = 0;
184 int i;
185 i = j * 8 / SHIFT;
186 if (i < a->size) {
187 r = a->d[i];
188 if (++i < a->size)
189 r |= (twodigits) a->d[i] << SHIFT;
190 }
191 r >>= ((j * 8) % SHIFT);
192 return r & 0xFF;
193 }
194
195 /* result must be free'd */
huge_as_binary(Huge * a,int * l)196 char *huge_as_binary (Huge * a, int *l)
197 {
198 char *s;
199 int i;
200 *l = (a->size * SHIFT) / 8 + 1;
201 s = malloc (*l + 1);
202 for (i = 0; i < *l; i++)
203 s[i] = _huge_get_char (a, i);
204 while (*l > 0 && !s[*l - 1])
205 (*l)--;
206 return s;
207 }
208
209 /* result must be free'd */
huge_from_binary(unsigned char * s,int n)210 Huge *huge_from_binary (unsigned char *s, int n)
211 {
212 Huge *z;
213 int i, size;
214 digit *d;
215 size = n * 8 / SHIFT;
216 z = huge_new (size + 1);
217 d = z->d;
218 for (i = 0; i < size + 1; i++) {
219 int j;
220 twodigits t = 0;
221 int r;
222 r = i * SHIFT / 8;
223 for (j = 0; j < SHIFT / 8 + 3 && r < n; j++, r++)
224 t |= (twodigits) s[r] << (j * 8);
225 t >>= ((i * SHIFT) % 8);
226 *d++ = (digit) t & MASK;
227 }
228 return huge_normalize (z);
229 }
230
231 /* Create a new long int object from a C unsigned long int */
232
huge_from_unsigned_long(unsigned long ival)233 Huge *huge_from_unsigned_long (unsigned long ival)
234 {
235 unsigned long t = ival;
236 int i;
237 /* Assume a C long fits in at most 5 'digits' */
238 /* Works on both 32- and 64-bit machines */
239 Huge *v = huge_new (5);
240 for (i = 0; i < 5; i++) {
241 v->ob_digit[i] = (digit) (t & MASK);
242 t >>= SHIFT;
243 }
244 return huge_normalize (v);
245 }
246
247 /* Get a C long int from a long int object.
248 Returns -1 and sets an error condition if overflow occurs. */
249
huge_as_long(Huge * v)250 long huge_as_long (Huge * v)
251 {
252 long x, prev;
253 int i, sign;
254
255 i = v->ob_size;
256 sign = 1;
257 x = 0;
258 if (i < 0) {
259 sign = -1;
260 i = -(i);
261 }
262 while (--i >= 0) {
263 prev = x;
264 x = (x << SHIFT) + v->ob_digit[i];
265 if ((x >> SHIFT) != prev) {
266 huge_error ("as_long(): long int too long to convert");
267 return -1;
268 }
269 }
270 return x * sign;
271 }
272
273 /* Get a C long int from a long int object.
274 Returns -1 and sets an error condition if overflow occurs. */
275
huge_as_unsigned_long(Huge * v)276 unsigned long huge_as_unsigned_long (Huge * v)
277 {
278 unsigned long x, prev;
279 int i;
280
281 i = v->ob_size;
282 x = 0;
283 if (i < 0) {
284 huge_error ("as_unsigned_long(): can't convert negative value to unsigned long");
285 return (unsigned long) -1;
286 }
287 while (--i >= 0) {
288 prev = x;
289 x = (x << SHIFT) + v->ob_digit[i];
290 if ((x >> SHIFT) != prev) {
291 huge_error ("as_unsigned_long(): long int too long to convert");
292 return (unsigned long) -1;
293 }
294 }
295 return x;
296 }
297
298 /* Get a C double from a long int object. */
299
300
301 /* Multiply by a single digit, ignoring the sign. */
302
mul1(Huge * a,wdigit n)303 static Huge *mul1 (Huge * a, wdigit n)
304 {
305 return muladd1 (a, n, (digit) 0);
306 }
307
308 /*
309 * gcc knows about 64bit product, so no optimisation needed:
310 *
311 * pushl -8(%ebp)
312 * pushl $.LC2
313 * call printf
314 *.stabn 68,0,47,.LM64-huge_from_long
315 *.LM64:
316 * pushl %edi
317 * pushl $.LC2
318 * call printf
319 *.stabn 68,0,48,.LM65-huge_from_long
320 *.LM65:
321 * movl -8(%ebp),%eax
322 * imull %edi
323 * movl %eax,-16(%ebp)
324 * movl %edx,-12(%ebp)
325 *.stabn 68,0,49,.LM66-huge_from_long
326 *.LM66:
327 * pushl -12(%ebp)
328 * pushl -16(%ebp)
329 * pushl $.LC2
330 * call printf
331 */
332
muladd1(Huge * a,wdigit n,wdigit extra)333 static Huge *muladd1 (Huge * a, wdigit n, wdigit extra)
334 {
335 int size_a = ABS (a->ob_size);
336 Huge *z = huge_new (size_a + 1);
337 twodigits carry = extra;
338 int i;
339 for (i = 0; i < size_a; ++i) {
340 carry += (twodigits) a->ob_digit[i] * n;
341 z->ob_digit[i] = (digit) (carry & MASK);
342 carry >>= SHIFT;
343 }
344 z->ob_digit[i] = (digit) carry;
345 return huge_normalize (z);
346 }
347
348 /* Divide a long integer by a digit, returning both the quotient
349 (as function result) and the remainder (through *prem).
350 The sign of a is ignored; n should not be zero. */
351
divrem1(Huge * a,wdigit n,digit * prem)352 static Huge *divrem1 (Huge * a, wdigit n, digit * prem)
353 {
354 int size = ABS (a->ob_size);
355 Huge *z;
356 int i;
357 twodigits rem = 0;
358
359 huge_assert (n > 0 && n <= MASK);
360 z = huge_new (size);
361 for (i = size; --i >= 0;) {
362 rem = (rem << SHIFT) + a->ob_digit[i];
363 z->ob_digit[i] = (digit) (rem / n);
364 rem %= n;
365 }
366 *prem = (digit) rem;
367 return huge_normalize (z);
368 }
369
370 /* Convert a long int object to a string, using a given conversion base.
371 Return a string object.
372
373 NDW: The following does not apply here ....
374 If base is 8 or 16, add the proper prefix '0' or '0x'.
375 External linkage: used in bltinmodule.c by hex() and oct(). */
376
huge_format(Huge * a,int base)377 char *huge_format (Huge * a, int base)
378 {
379 char *str;
380 int i;
381 int size_a = ABS (a->ob_size);
382 char *p;
383 int bits;
384 char sign = '\0';
385
386 a = huge_dup (a);
387 huge_assert (base >= 2 && base <= 36);
388
389 /* Compute a rough upper bound for the length of the string */
390 i = base;
391 bits = 0;
392 while (i > 1) {
393 ++bits;
394 i >>= 1;
395 }
396 i = 6 + (size_a * SHIFT + bits - 1) / bits;
397 str = malloc (i + 1);
398 p = str + i;
399 *p = '\0';
400 #ifdef ORIGINAL_BEHAVIOUR
401 *--p = 'L';
402 #endif
403 if (a->ob_size < 0) {
404 sign = '-';
405 a->ob_size = ABS (a->ob_size);
406 }
407
408 do {
409 digit rem;
410 Huge *temp = divrem1 (a, (digit) base, &rem);
411 if (temp == 0) {
412 huge_free (a);
413 free (str);
414 return 0;
415 }
416 if (rem < 10)
417 rem += '0';
418 else
419 rem += 'a' - 10;
420 huge_assert (p > str);
421 *--p = (char) rem;
422 huge_free (a);
423 a = temp;
424 } while (ABS (a->ob_size) != 0);
425 huge_free (a);
426 #ifdef ORIGINAL_BEHAVIOUR
427 /* NDW -- removed this for GMP compatibility */
428 if (base == 8) {
429 if (size_a != 0)
430 *--p = '0';
431 } else if (base == 16) {
432 *--p = 'x';
433 *--p = '0';
434 } else if (base != 10) {
435 *--p = '#';
436 *--p = '0' + base % 10;
437 if (base > 10)
438 *--p = '0' + base / 10;
439 }
440 #endif
441 if (sign)
442 *--p = sign;
443 if (p != str) {
444 char *q = str;
445 huge_assert (p > q);
446 do {
447 } while ((*q++ = *p++) != '\0');
448 q--;
449 }
450 return str;
451 }
452
huge_from_string(char * str,char ** pend,int base)453 Huge *huge_from_string (char *str, char **pend, int base)
454 {
455 int sign = 1;
456 Huge *z;
457
458 while (*str != '\0' && strchr ("\t\n ", *str))
459 str++;
460 if (*str == '+')
461 ++str;
462 else if (*str == '-') {
463 ++str;
464 sign = -1;
465 }
466 while (*str != '\0' && strchr ("\t\n ", *str))
467 str++;
468 if (base == 0) {
469 if (str[0] != '0')
470 base = 10;
471 else if (str[1] == 'x' || str[1] == 'X')
472 base = 16;
473 else
474 base = 8;
475 }
476 if (base == 16 && str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
477 str += 2;
478 z = huge_new (0);
479 for (; z != 0; ++str) {
480 int k = -1;
481 Huge *temp;
482
483 if (*str <= '9')
484 k = *str - '0';
485 else if (*str >= 'a')
486 k = *str - 'a' + 10;
487 else if (*str >= 'A')
488 k = *str - 'A' + 10;
489 if (k < 0 || k >= base)
490 break;
491 temp = muladd1 (z, (digit) base, (digit) k);
492 huge_free (z);
493 z = temp;
494 }
495 if (sign < 0 && z != 0 && z->ob_size != 0)
496 z->ob_size = -(z->ob_size);
497 if (pend)
498 *pend = str;
499 return huge_normalize (z);
500 }
501
502 /* Long division with remainder, top-level routine */
503
_huge_divrem(Huge * a,Huge * b,Huge ** pdiv,Huge ** prem)504 static int _huge_divrem (Huge * a, Huge * b, Huge ** pdiv, Huge ** prem)
505 {
506 int size_a = ABS (a->ob_size), size_b = ABS (b->ob_size);
507 Huge *z;
508
509 if (!size_b)
510 huge_error ("divrem(): divide by zero");
511 if (size_a < size_b ||
512 (size_a == size_b &&
513 a->ob_digit[size_a - 1] < b->ob_digit[size_b - 1])) {
514 /* |a| < |b|. */
515 *pdiv = huge_new (0);
516 *prem = huge_dup (a);
517 return 0;
518 }
519 if (size_b == 1) {
520 digit rem = 0;
521 z = divrem1 (a, b->ob_digit[0], &rem);
522 if (z == 0)
523 return -1;
524 *prem = huge_from_long ((long) rem);
525 } else {
526 z = x_divrem (a, b, prem);
527 if (z == 0)
528 return -1;
529 }
530 /* Set the signs.
531 The quotient z has the sign of a*b;
532 the remainder r has the sign of a,
533 so a = b*z + r. */
534 if ((a->ob_size < 0) != (b->ob_size < 0))
535 z->ob_size = -(z->ob_size);
536 if (a->ob_size < 0 && (*prem)->ob_size != 0)
537 (*prem)->ob_size = -((*prem)->ob_size);
538 *pdiv = z;
539 return 0;
540 }
541
542 /* Unsigned long division with remainder -- the algorithm */
543
x_divrem(Huge * v1,Huge * w1,Huge ** prem)544 static Huge *x_divrem (Huge * v1, Huge * w1, Huge ** prem)
545 {
546 int size_v = ABS (v1->ob_size), size_w = ABS (w1->ob_size);
547 digit d = (digit) ((twodigits) BASE / (w1->ob_digit[size_w - 1] + 1));
548 Huge *v = mul1 (v1, d);
549 Huge *w = mul1 (w1, d);
550 Huge *a;
551 int j, k;
552
553 if (v == 0 || w == 0) {
554 huge_free (v);
555 huge_free (w);
556 return 0;
557 }
558 huge_assert (size_v >= size_w && size_w > 1); /* Assert checks by div() */
559 huge_assert (size_w == ABS (w->ob_size)); /* That's how d was calculated */
560
561 size_v = ABS (v->ob_size);
562 a = huge_new (size_v - size_w + 1);
563
564 for (j = size_v, k = a->ob_size - 1; a != 0 && k >= 0; --j, --k) {
565 digit vj = (j >= size_v) ? 0 : v->ob_digit[j];
566 twodigits q;
567 stwodigits carry = 0;
568 int i;
569
570 if (vj == w->ob_digit[size_w - 1])
571 q = MASK;
572 else
573 q = (((twodigits) vj << SHIFT) + v->ob_digit[j - 1]) /
574 w->ob_digit[size_w - 1];
575
576 while (w->ob_digit[size_w - 2] * q >
577 ((
578 ((twodigits) vj << SHIFT)
579 + v->ob_digit[j - 1]
580 - q * w->ob_digit[size_w - 1]
581 ) << SHIFT)
582 + v->ob_digit[j - 2])
583 --q;
584
585 for (i = 0; i < size_w && i + k < size_v; ++i) {
586 twodigits z = w->ob_digit[i] * q;
587 digit zz = (digit) (z >> SHIFT);
588 carry += v->ob_digit[i + k] - z
589 + ((twodigits) zz << SHIFT);
590 v->ob_digit[i + k] = carry & MASK;
591 carry = (carry >> SHIFT) - zz;
592 }
593
594 if (i + k < size_v) {
595 carry += v->ob_digit[i + k];
596 v->ob_digit[i + k] = 0;
597 }
598 if (carry == 0)
599 a->ob_digit[k] = (digit) q;
600 else {
601 huge_assert (carry == -1);
602 a->ob_digit[k] = (digit) q - 1;
603 carry = 0;
604 for (i = 0; i < size_w && i + k < size_v; ++i) {
605 carry += v->ob_digit[i + k] + w->ob_digit[i];
606 v->ob_digit[i + k] = carry & MASK;
607 carry >>= SHIFT;
608 }
609 }
610 } /* for j, k */
611
612 if (a == 0)
613 *prem = 0;
614 else {
615 a = huge_normalize (a);
616 *prem = divrem1 (v, d, &d);
617 /* d receives the (unused) remainder */
618 if (*prem == 0) {
619 huge_free (a);
620 a = 0;
621 }
622 }
623 huge_free (v);
624 huge_free (w);
625 return a;
626 }
627
huge_compare(Huge * a,Huge * b)628 int huge_compare (Huge * a, Huge * b)
629 {
630 int sign;
631
632 if (a->ob_size != b->ob_size) {
633 if (ABS (a->ob_size) == 0 && ABS (b->ob_size) == 0)
634 sign = 0;
635 else
636 sign = a->ob_size - b->ob_size;
637 } else {
638 int i = ABS (a->ob_size);
639 while (--i >= 0 && a->ob_digit[i] == b->ob_digit[i]);
640 if (i < 0)
641 sign = 0;
642 else {
643 sign = (int) a->ob_digit[i] - (int) b->ob_digit[i];
644 if (a->ob_size < 0)
645 sign = -sign;
646 }
647 }
648 return sign < 0 ? -1 : sign > 0 ? 1 : 0;
649 }
650
651 /* Add the absolute values of two long integers. */
652
x_add(Huge * a,Huge * b)653 static Huge *x_add (Huge * a, Huge * b)
654 {
655 int size_a = ABS (a->ob_size), size_b = ABS (b->ob_size);
656 Huge *z;
657 int i;
658 digit carry = 0;
659
660 /* Ensure a is the larger of the two: */
661 if (size_a < size_b) {
662 {
663 Huge *temp = a;
664 a = b;
665 b = temp;
666 }
667 {
668 int size_temp = size_a;
669 size_a = size_b;
670 size_b = size_temp;
671 }
672 }
673 z = huge_new (size_a + 1);
674 for (i = 0; i < size_b; ++i) {
675 carry += a->ob_digit[i] + b->ob_digit[i];
676 z->ob_digit[i] = carry & MASK;
677 /* The following assumes unsigned shifts don't
678 propagate the sign bit. */
679 carry >>= SHIFT;
680 }
681 for (; i < size_a; ++i) {
682 carry += a->ob_digit[i];
683 z->ob_digit[i] = carry & MASK;
684 carry >>= SHIFT;
685 }
686 z->ob_digit[i] = carry;
687 return huge_normalize (z);
688 }
689
690 /* Subtract the absolute values of two integers. */
691
x_sub(Huge * a,Huge * b)692 static Huge *x_sub (Huge * a, Huge * b)
693 {
694 int size_a = ABS (a->ob_size), size_b = ABS (b->ob_size);
695 Huge *z;
696 int i;
697 int sign = 1;
698 digit borrow = 0;
699
700 /* Ensure a is the larger of the two: */
701 if (size_a < size_b) {
702 sign = -1;
703 {
704 Huge *temp = a;
705 a = b;
706 b = temp;
707 }
708 {
709 int size_temp = size_a;
710 size_a = size_b;
711 size_b = size_temp;
712 }
713 } else if (size_a == size_b) {
714 /* Find highest digit where a and b differ: */
715 i = size_a;
716 while (--i >= 0 && a->ob_digit[i] == b->ob_digit[i]);
717 if (i < 0)
718 return huge_new (0);
719 if (a->ob_digit[i] < b->ob_digit[i]) {
720 sign = -1;
721 {
722 Huge *temp = a;
723 a = b;
724 b = temp;
725 }
726 }
727 size_a = size_b = i + 1;
728 }
729 z = huge_new (size_a);
730 for (i = 0; i < size_b; ++i) {
731 /* The following assumes unsigned arithmetic
732 works module 2**N for some N>SHIFT. */
733 borrow = a->ob_digit[i] - b->ob_digit[i] - borrow;
734 z->ob_digit[i] = borrow & MASK;
735 borrow >>= SHIFT;
736 borrow &= 1; /* Keep only one sign bit */
737 }
738 for (; i < size_a; ++i) {
739 borrow = a->ob_digit[i] - borrow;
740 z->ob_digit[i] = borrow & MASK;
741 borrow >>= SHIFT;
742 }
743 huge_assert (borrow == 0);
744 if (sign < 0)
745 z->ob_size = -(z->ob_size);
746 return huge_normalize (z);
747 }
748
huge_add(Huge * a,Huge * b)749 Huge *huge_add (Huge * a, Huge * b)
750 {
751 Huge *z;
752
753 if (a->ob_size < 0) {
754 if (b->ob_size < 0) {
755 z = x_add (a, b);
756 if (z != 0 && z->ob_size != 0)
757 z->ob_size = -(z->ob_size);
758 } else
759 z = x_sub (b, a);
760 } else {
761 if (b->ob_size < 0)
762 z = x_sub (a, b);
763 else
764 z = x_add (a, b);
765 }
766 return (Huge *) z;
767 }
768
huge_sub(Huge * a,Huge * b)769 Huge *huge_sub (Huge * a, Huge * b)
770 {
771 Huge *z;
772
773 if (a->ob_size < 0) {
774 if (b->ob_size < 0)
775 z = x_sub (a, b);
776 else
777 z = x_add (a, b);
778 if (z != 0 && z->ob_size != 0)
779 z->ob_size = -(z->ob_size);
780 } else {
781 if (b->ob_size < 0)
782 z = x_add (a, b);
783 else
784 z = x_sub (a, b);
785 }
786 return (Huge *) z;
787 }
788
huge_mul(Huge * a,Huge * b)789 Huge *huge_mul (Huge * a, Huge * b)
790 {
791 int size_a;
792 int size_b;
793 Huge *z;
794 int i;
795
796 size_a = ABS (a->ob_size);
797 size_b = ABS (b->ob_size);
798 z = huge_new (size_a + size_b);
799 for (i = 0; i < z->ob_size; ++i)
800 z->ob_digit[i] = 0;
801 for (i = 0; i < size_a; ++i) {
802 twodigits carry = 0;
803 twodigits f = a->ob_digit[i];
804 int j;
805 for (j = 0; j < size_b; ++j) {
806 carry += z->ob_digit[i + j] + b->ob_digit[j] * f;
807 z->ob_digit[i + j] = (digit) (carry & MASK);
808 carry >>= SHIFT;
809 }
810 for (; carry != 0; ++j) {
811 huge_assert (i + j < z->ob_size);
812 carry += z->ob_digit[i + j];
813 z->ob_digit[i + j] = (digit) (carry & MASK);
814 carry >>= SHIFT;
815 }
816 }
817 if (a->ob_size < 0)
818 z->ob_size = -(z->ob_size);
819 if (b->ob_size < 0)
820 z->ob_size = -(z->ob_size);
821 return (Huge *) huge_normalize (z);
822 }
823
824 /* The / and % operators are now defined in terms of divmod().
825 The expression a mod b has the value a - b*floor(a/b).
826 The huge_divrem function gives the remainder after division of
827 |a| by |b|, with the sign of a. This is also expressed
828 as a - b*trunc(a/b), if trunc truncates towards zero.
829 Some examples:
830 a b a rem b a mod b
831 13 10 3 3
832 -13 10 -3 7
833 13 -10 3 -7
834 -13 -10 -3 -3
835 So, to get from rem to mod, we have to add b if a and b
836 have different signs. We then subtract one from the 'divisor'
837 part of the outcome to keep the invariant intact. */
838
l_divmod(Huge * v,Huge * w,Huge ** pdiv,Huge ** pmod)839 static int l_divmod (Huge * v, Huge * w, Huge ** pdiv, Huge ** pmod)
840 {
841 Huge *divisor, *mod;
842
843 if (_huge_divrem (v, w, &divisor, &mod) < 0)
844 return -1;
845 if ((mod->ob_size < 0 && w->ob_size > 0) ||
846 (mod->ob_size > 0 && w->ob_size < 0)) {
847 Huge *temp;
848 Huge *one;
849 temp = (Huge *) huge_add (mod, w);
850 huge_free (mod);
851 mod = temp;
852 if (mod == 0) {
853 huge_free (divisor);
854 return -1;
855 }
856 one = huge_from_long (1L);
857 if ((temp = (Huge *) huge_sub (divisor, one)) == 0) {
858 huge_free (mod);
859 huge_free (divisor);
860 huge_free (one);
861 return -1;
862 }
863 huge_free (one);
864 huge_free (divisor);
865 divisor = temp;
866 }
867 *pdiv = divisor;
868 *pmod = mod;
869 return 0;
870 }
871
huge_div(Huge * v,Huge * w)872 Huge *huge_div (Huge * v, Huge * w)
873 {
874 Huge *divisor, *mod;
875 if (l_divmod (v, w, &divisor, &mod) < 0)
876 return 0;
877 huge_free (mod);
878 return (Huge *) divisor;
879 }
880
huge_mod(Huge * v,Huge * w)881 Huge *huge_mod (Huge * v, Huge * w)
882 {
883 Huge *divisor, *mod;
884 if (l_divmod (v, w, &divisor, &mod) < 0)
885 return 0;
886 huge_free (divisor);
887 return (Huge *) mod;
888 }
889
huge_divmod(Huge * v,Huge * w,Huge ** remainder)890 Huge *huge_divmod (Huge * v, Huge * w, Huge ** remainder)
891 {
892 Huge *divisor, *mod;
893 if (l_divmod (v, w, &divisor, &mod) < 0)
894 return 0;
895 if (remainder)
896 *remainder = mod;
897 return divisor;
898 }
899
huge_powmod(Huge * a,Huge * b,Huge * c)900 Huge *huge_powmod (Huge * a, Huge * b, Huge * c)
901 {
902 Huge *z = 0, *divisor = 0, *mod = 0;
903 int size_b, i;
904
905 a = huge_dup (a);
906 size_b = b->ob_size;
907 if (size_b < 0) {
908 huge_error ("pow(): long integer to the negative power");
909 return 0;
910 }
911 z = (Huge *) huge_from_long (1L);
912 for (i = 0; i < size_b; ++i) {
913 digit bi = b->ob_digit[i];
914 int j;
915
916 for (j = 0; j < SHIFT; ++j) {
917 Huge *temp = 0;
918
919 if (bi & 1) {
920 temp = (Huge *) huge_mul (z, a);
921 huge_free (z);
922 if (c != 0 && temp != 0) {
923 l_divmod (temp, c, &divisor, &mod);
924 huge_free (divisor);
925 huge_free (temp);
926 temp = mod;
927 }
928 z = temp;
929 if (z == 0)
930 break;
931 }
932 bi >>= 1;
933 if (bi == 0 && i + 1 == size_b)
934 break;
935 temp = (Huge *) huge_mul (a, a);
936 huge_free (a);
937 if ((Huge *) c != 0 && temp != 0) {
938 l_divmod (temp, c, &divisor, &mod);
939 huge_free (divisor);
940 huge_free (temp);
941 temp = mod;
942 }
943 a = temp;
944 if (a == 0) {
945 huge_free (z);
946 z = 0;
947 break;
948 }
949 }
950 if (a == 0 || z == 0)
951 break;
952 }
953 huge_free (a);
954 if ((Huge *) c != 0 && z != 0) {
955 l_divmod (z, c, &divisor, &mod);
956 huge_free (divisor);
957 huge_free (z);
958 z = mod;
959 }
960 return (Huge *) z;
961 }
962
huge_pow(Huge * a,Huge * b)963 Huge *huge_pow (Huge * a, Huge * b)
964 {
965 return huge_powmod (a, b, 0);
966 }
967
huge_invert(Huge * v)968 Huge *huge_invert (Huge * v)
969 {
970 /* Implement ~x as -(x+1) */
971 Huge *x;
972 Huge *w;
973 w = (Huge *) huge_from_long (1L);
974 if (w == 0)
975 return 0;
976 x = (Huge *) huge_add (v, w);
977 huge_free (w);
978 if (x == 0)
979 return 0;
980 if (x->ob_size != 0)
981 x->ob_size = -(x->ob_size);
982 return (Huge *) x;
983 }
984
huge_neg(Huge * v)985 Huge *huge_neg (Huge * v)
986 {
987 Huge *z;
988 int i, n;
989 n = ABS (v->ob_size);
990 /* -0 == 0 */
991 if (!n)
992 return huge_dup (v);
993 z = huge_new (n);
994 for (i = 0; i < n; i++)
995 z->ob_digit[i] = v->ob_digit[i];
996 z->ob_size = -(v->ob_size);
997 return (Huge *) z;
998 }
999
huge_abs(Huge * v)1000 Huge *huge_abs (Huge * v)
1001 {
1002 if (v->ob_size < 0)
1003 return huge_neg (v);
1004 else
1005 return huge_dup (v);
1006 }
1007
huge_nonzero(Huge * v)1008 int huge_nonzero (Huge * v)
1009 {
1010 if (!v)
1011 return 0;
1012 return v->ob_size != 0;
1013 }
1014
huge_rshift(Huge * a,int shiftby)1015 Huge *huge_rshift (Huge * a, int shiftby)
1016 {
1017 Huge *z;
1018 int newsize, wordshift, loshift, hishift, i, j;
1019 digit lomask, himask;
1020
1021 if (a->ob_size < 0) {
1022 /* Right shifting negative numbers is harder */
1023 Huge *a1, *a2, *a3;
1024 a1 = (Huge *) huge_invert (a);
1025 if (a1 == 0)
1026 return 0;
1027 a2 = (Huge *) huge_rshift (a1, shiftby);
1028 huge_free (a1);
1029 if (a2 == 0)
1030 return 0;
1031 a3 = (Huge *) huge_invert (a2);
1032 huge_free (a2);
1033 return (Huge *) a3;
1034 }
1035 if (shiftby < 0) {
1036 huge_error ("rshift(): negative shift count");
1037 return 0;
1038 }
1039 wordshift = shiftby / SHIFT;
1040 newsize = ABS (a->ob_size) - wordshift;
1041 if (newsize <= 0) {
1042 z = huge_new (0);
1043 return (Huge *) z;
1044 }
1045 loshift = shiftby % SHIFT;
1046 hishift = SHIFT - loshift;
1047 lomask = ((digit) 1 << hishift) - 1;
1048 himask = MASK ^ lomask;
1049 z = huge_new (newsize);
1050 if (a->ob_size < 0)
1051 z->ob_size = -(z->ob_size);
1052 for (i = 0, j = wordshift; i < newsize; i++, j++) {
1053 z->ob_digit[i] = (a->ob_digit[j] >> loshift) & lomask;
1054 if (i + 1 < newsize)
1055 z->ob_digit[i] |=
1056 (a->ob_digit[j + 1] << hishift) & himask;
1057 }
1058 return (Huge *) huge_normalize (z);
1059 }
1060
huge_lshift(Huge * a,int shiftby)1061 Huge *huge_lshift (Huge * a, int shiftby)
1062 {
1063 /* This version due to Tim Peters */
1064 Huge *z;
1065 int oldsize, newsize, wordshift, remshift, i, j;
1066 twodigits accum;
1067
1068 if (shiftby < 0) {
1069 huge_error ("lshift(): negative shift count");
1070 return 0;
1071 }
1072 /* wordshift, remshift = divmod(shiftby, SHIFT) */
1073 wordshift = (int) shiftby / SHIFT;
1074 remshift = (int) shiftby - wordshift * SHIFT;
1075
1076 oldsize = ABS (a->ob_size);
1077 newsize = oldsize + wordshift;
1078 if (remshift)
1079 ++newsize;
1080 z = huge_new (newsize);
1081 if (a->ob_size < 0)
1082 z->ob_size = -(z->ob_size);
1083 for (i = 0; i < wordshift; i++)
1084 z->ob_digit[i] = 0;
1085 accum = 0;
1086 for (i = wordshift, j = 0; j < oldsize; i++, j++) {
1087 accum |= a->ob_digit[j] << remshift;
1088 z->ob_digit[i] = (digit) (accum & MASK);
1089 accum >>= SHIFT;
1090 }
1091 if (remshift)
1092 z->ob_digit[newsize - 1] = (digit) accum;
1093 else
1094 huge_assert (!accum);
1095 return (Huge *) huge_normalize (z);
1096 }
1097
1098
1099 /* Bitwise and/xor/or operations */
1100
1101 /* op = '&', '|', '^' */
huge_bitwise(Huge * a,int op,Huge * b)1102 static Huge *huge_bitwise (Huge * a, int op, Huge * b)
1103 {
1104 digit maska, maskb; /* 0 or MASK */
1105 int negz;
1106 int size_a, size_b, size_z;
1107 Huge *z;
1108 int i;
1109 digit diga, digb;
1110 Huge *v;
1111
1112 if (a->ob_size < 0) {
1113 a = (Huge *) huge_invert (a);
1114 maska = MASK;
1115 } else {
1116 a = huge_dup (a);
1117 maska = 0;
1118 }
1119 if (b->ob_size < 0) {
1120 b = (Huge *) huge_invert (b);
1121 maskb = MASK;
1122 } else {
1123 b = huge_dup (b);
1124 maskb = 0;
1125 }
1126
1127 size_a = a->ob_size;
1128 size_b = b->ob_size;
1129 size_z = MAX (size_a, size_b);
1130 z = huge_new (size_z);
1131 if (a == 0 || b == 0) {
1132 huge_free (a);
1133 huge_free (b);
1134 huge_free (z);
1135 return 0;
1136 }
1137 negz = 0;
1138 switch (op) {
1139 case '^':
1140 if (maska != maskb) {
1141 maska ^= MASK;
1142 negz = -1;
1143 }
1144 break;
1145 case '&':
1146 if (maska && maskb) {
1147 op = '|';
1148 maska ^= MASK;
1149 maskb ^= MASK;
1150 negz = -1;
1151 }
1152 break;
1153 case '|':
1154 if (maska || maskb) {
1155 op = '&';
1156 maska ^= MASK;
1157 maskb ^= MASK;
1158 negz = -1;
1159 }
1160 break;
1161 }
1162
1163 for (i = 0; i < size_z; ++i) {
1164 diga = (i < size_a ? a->ob_digit[i] : 0) ^ maska;
1165 digb = (i < size_b ? b->ob_digit[i] : 0) ^ maskb;
1166 switch (op) {
1167 case '&':
1168 z->ob_digit[i] = diga & digb;
1169 break;
1170 case '|':
1171 z->ob_digit[i] = diga | digb;
1172 break;
1173 case '^':
1174 z->ob_digit[i] = diga ^ digb;
1175 break;
1176 }
1177 }
1178
1179 huge_free (a);
1180 huge_free (b);
1181 z = huge_normalize (z);
1182 if (negz == 0)
1183 return (Huge *) z;
1184 v = huge_invert (z);
1185 huge_free (z);
1186 return v;
1187 }
1188
huge_and(Huge * a,Huge * b)1189 Huge *huge_and (Huge * a, Huge * b)
1190 {
1191 return huge_bitwise (a, '&', b);
1192 }
1193
huge_xor(Huge * a,Huge * b)1194 Huge *huge_xor (Huge * a, Huge * b)
1195 {
1196 return huge_bitwise (a, '^', b);
1197 }
1198
huge_or(Huge * a,Huge * b)1199 Huge *huge_or (Huge * a, Huge * b)
1200 {
1201 return huge_bitwise (a, '|', b);
1202 }
1203
huge_oct(Huge * v)1204 char *huge_oct (Huge * v)
1205 {
1206 return huge_format (v, 8);
1207 }
1208
huge_hex(Huge * v)1209 char *huge_hex (Huge * v)
1210 {
1211 return huge_format (v, 16);
1212 }
1213
huge_dec(Huge * v)1214 char *huge_dec (Huge * v)
1215 {
1216 return huge_format (v, 10);
1217 }
1218
xhuge_log(Huge * h,char * msg,char * file,int line)1219 void xhuge_log(Huge *h, char *msg, char *file, int line)
1220 {
1221 static FILE *f = 0;
1222 char *p = 0;
1223 if (!f)
1224 f = fopen ("huge.log", "w");
1225 fprintf (f, "%s: %d:\n%s: %s\n", file, line, msg, p = huge_hex(h));
1226 fflush (f);
1227 if (p)
1228 free (p);
1229 }
1230
1231
1232
1233