longobject.c revision 8b27d929a606b9723efb6cf3045a5a4b5418acb3
1/***********************************************************
2Copyright 1991 by Stichting Mathematisch Centrum, Amsterdam, The
3Netherlands.
4
5                        All Rights Reserved
6
7Permission to use, copy, modify, and distribute this software and its
8documentation for any purpose and without fee is hereby granted,
9provided that the above copyright notice appear in all copies and that
10both that copyright notice and this permission notice appear in
11supporting documentation, and that the names of Stichting Mathematisch
12Centrum or CWI not be used in advertising or publicity pertaining to
13distribution of the software without specific, written prior permission.
14
15STICHTING MATHEMATISCH CENTRUM DISCLAIMS ALL WARRANTIES WITH REGARD TO
16THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
17FITNESS, IN NO EVENT SHALL STICHTING MATHEMATISCH CENTRUM BE LIABLE
18FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
19WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
20ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT
21OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
22
23******************************************************************/
24
25/* Long (arbitrary precision) integer object implementation */
26
27/* XXX The functional organization of this file is terrible */
28
29#include "allobjects.h"
30#include "intrcheck.h"
31#include "longintrepr.h"
32#include <math.h>
33#include <assert.h>
34
35#define ABS(x) ((x) < 0 ? -(x) : (x))
36
37/* Forward */
38static longobject *long_normalize PROTO((longobject *));
39static longobject *mul1 PROTO((longobject *, wdigit));
40static longobject *muladd1 PROTO((longobject *, wdigit, wdigit));
41static longobject *divrem1 PROTO((longobject *, wdigit, digit *));
42
43static int ticker;	/* XXX Could be shared with ceval? */
44
45#define INTRCHECK(block) \
46	if (--ticker < 0) { \
47		ticker = 100; \
48		if (intrcheck()) { block; } \
49	}
50
51/* Normalize (remove leading zeros from) a long int object.
52   Doesn't attempt to free the storage--in most cases, due to the nature
53   of the algorithms used, this could save at most be one word anyway. */
54
55static longobject *
56long_normalize(v)
57	register longobject *v;
58{
59	int j = ABS(v->ob_size);
60	register int i = j;
61
62	while (i > 0 && v->ob_digit[i-1] == 0)
63		--i;
64	if (i != j)
65		v->ob_size = (v->ob_size < 0) ? -(i) : i;
66	return v;
67}
68
69/* Allocate a new long int object with size digits.
70   Return NULL and set exception if we run out of memory. */
71
72longobject *
73alloclongobject(size)
74	int size;
75{
76	return NEWVAROBJ(longobject, &Longtype, size);
77}
78
79/* Create a new long int object from a C long int */
80
81object *
82newlongobject(ival)
83	long ival;
84{
85	/* Assume a C long fits in at most 3 'digits' */
86	longobject *v = alloclongobject(3);
87	if (v != NULL) {
88		if (ival < 0) {
89			ival = -ival;
90			v->ob_size = -(v->ob_size);
91		}
92		v->ob_digit[0] = ival & MASK;
93		v->ob_digit[1] = (ival >> SHIFT) & MASK;
94		v->ob_digit[2] = (ival >> (2*SHIFT)) & MASK;
95		v = long_normalize(v);
96	}
97	return (object *)v;
98}
99
100/* Create a new long int object from a C double */
101
102object *
103dnewlongobject(dval)
104	double dval;
105{
106	longobject *v;
107	double frac;
108	int i, ndig, expo, neg;
109	neg = 0;
110	if (dval < 0.0) {
111		neg = 1;
112		dval = -dval;
113	}
114	frac = frexp(dval, &expo); /* dval = frac*2**expo; 0.0 <= frac < 1.0 */
115	if (expo <= 0)
116		return newlongobject(0L);
117	ndig = (expo-1) / SHIFT + 1; /* Number of 'digits' in result */
118	v = alloclongobject(ndig);
119	if (v == NULL)
120		return NULL;
121	frac = ldexp(frac, (expo-1) % SHIFT + 1);
122	for (i = ndig; --i >= 0; ) {
123		long bits = (long)frac;
124		v->ob_digit[i] = bits;
125		frac = frac - (double)bits;
126		frac = ldexp(frac, SHIFT);
127	}
128	if (neg)
129		v->ob_size = -(v->ob_size);
130	return (object *)v;
131}
132
133/* Get a C long int from a long int object.
134   Returns -1 and sets an error condition if overflow occurs. */
135
136long
137getlongvalue(vv)
138	object *vv;
139{
140	register longobject *v;
141	long x, prev;
142	int i, sign;
143
144	if (vv == NULL || !is_longobject(vv)) {
145		err_badcall();
146		return -1;
147	}
148	v = (longobject *)vv;
149	i = v->ob_size;
150	sign = 1;
151	x = 0;
152	if (i < 0) {
153		sign = -1;
154		i = -(i);
155	}
156	while (--i >= 0) {
157		prev = x;
158		x = (x << SHIFT) + v->ob_digit[i];
159		if ((x >> SHIFT) != prev) {
160			err_setstr(ValueError,
161				"long int too long to convert");
162			return -1;
163		}
164	}
165	return x * sign;
166}
167
168/* Get a C double from a long int object.  No overflow check. */
169
170double
171dgetlongvalue(vv)
172	object *vv;
173{
174	register longobject *v;
175	double x;
176	double multiplier = (double) (1L << SHIFT);
177	int i, sign;
178
179	if (vv == NULL || !is_longobject(vv)) {
180		err_badcall();
181		return -1;
182	}
183	v = (longobject *)vv;
184	i = v->ob_size;
185	sign = 1;
186	x = 0.0;
187	if (i < 0) {
188		sign = -1;
189		i = -(i);
190	}
191	while (--i >= 0) {
192		x = x*multiplier + v->ob_digit[i];
193	}
194	return x * sign;
195}
196
197/* Multiply by a single digit, ignoring the sign. */
198
199static longobject *
200mul1(a, n)
201	longobject *a;
202	wdigit n;
203{
204	return muladd1(a, n, (digit)0);
205}
206
207/* Multiply by a single digit and add a single digit, ignoring the sign. */
208
209static longobject *
210muladd1(a, n, extra)
211	longobject *a;
212	wdigit n;
213	wdigit extra;
214{
215	int size_a = ABS(a->ob_size);
216	longobject *z = alloclongobject(size_a+1);
217	twodigits carry = extra;
218	int i;
219
220	if (z == NULL)
221		return NULL;
222	for (i = 0; i < size_a; ++i) {
223		carry += (twodigits)a->ob_digit[i] * n;
224		z->ob_digit[i] = carry & MASK;
225		carry >>= SHIFT;
226	}
227	z->ob_digit[i] = carry;
228	return long_normalize(z);
229}
230
231/* Divide a long integer by a digit, returning both the quotient
232   (as function result) and the remainder (through *prem).
233   The sign of a is ignored; n should not be zero. */
234
235static longobject *
236divrem1(a, n, prem)
237	longobject *a;
238	wdigit n;
239	digit *prem;
240{
241	int size = ABS(a->ob_size);
242	longobject *z;
243	int i;
244	twodigits rem = 0;
245
246	assert(n > 0 && n <= MASK);
247	z = alloclongobject(size);
248	if (z == NULL)
249		return NULL;
250	for (i = size; --i >= 0; ) {
251		rem = (rem << SHIFT) + a->ob_digit[i];
252		z->ob_digit[i] = rem/n;
253		rem %= n;
254	}
255	*prem = rem;
256	return long_normalize(z);
257}
258
259/* Convert a long int object to a string, using a given conversion base.
260   Return a string object.
261   If base is 8 or 16, add the proper prefix '0' or '0x'.
262   External linkage: used in bltinmodule.c by hex() and oct(). */
263
264object *
265long_format(aa, base)
266	object *aa;
267	int base;
268{
269	register longobject *a = (longobject *)aa;
270	stringobject *str;
271	int i;
272	int size_a = ABS(a->ob_size);
273	char *p;
274	int bits;
275	char sign = '\0';
276
277	if (a == NULL || !is_longobject(a)) {
278		err_badcall();
279		return NULL;
280	}
281	assert(base >= 2 && base <= 36);
282
283	/* Compute a rough upper bound for the length of the string */
284	i = base;
285	bits = 0;
286	while (i > 1) {
287		++bits;
288		i >>= 1;
289	}
290	i = 6 + (size_a*SHIFT + bits-1) / bits;
291	str = (stringobject *) newsizedstringobject((char *)0, i);
292	if (str == NULL)
293		return NULL;
294	p = GETSTRINGVALUE(str) + i;
295	*p = '\0';
296	*--p = 'L';
297	if (a->ob_size < 0)
298		sign = '-';
299
300	INCREF(a);
301	do {
302		digit rem;
303		longobject *temp = divrem1(a, (digit)base, &rem);
304		if (temp == NULL) {
305			DECREF(a);
306			DECREF(str);
307			return NULL;
308		}
309		if (rem < 10)
310			rem += '0';
311		else
312			rem += 'A'-10;
313		assert(p > GETSTRINGVALUE(str));
314		*--p = rem;
315		DECREF(a);
316		a = temp;
317		INTRCHECK({
318			DECREF(a);
319			DECREF(str);
320			err_set(KeyboardInterrupt);
321			return NULL;
322		})
323	} while (ABS(a->ob_size) != 0);
324	DECREF(a);
325	if (base == 8)
326		*--p = '0';
327	else if (base == 16) {
328		*--p = 'x';
329		*--p = '0';
330	}
331	else if (base != 10) {
332		*--p = '#';
333		*--p = '0' + base%10;
334		if (base > 10)
335			*--p = '0' + base/10;
336	}
337	if (sign)
338		*--p = sign;
339	if (p != GETSTRINGVALUE(str)) {
340		char *q = GETSTRINGVALUE(str);
341		assert(p > q);
342		do {
343		} while ((*q++ = *p++) != '\0');
344		q--;
345		resizestring((object **)&str, (int) (q - GETSTRINGVALUE(str)));
346	}
347	return (object *)str;
348}
349
350/* Convert a string to a long int object, in a given base.
351   Base zero implies a default depending on the number.
352   External linkage: used in compile.c for literals. */
353
354object *
355long_scan(str, base)
356	char *str;
357	int base;
358{
359	int sign = 1;
360	longobject *z;
361
362	assert(base == 0 || base >= 2 && base <= 36);
363	if (*str == '+')
364		++str;
365	else if (*str == '-') {
366		++str;
367		sign = -1;
368	}
369	if (base == 0) {
370		if (str[0] != '0')
371			base = 10;
372		else if (str[1] == 'x' || str[1] == 'X')
373			base = 16;
374		else
375			base = 8;
376	}
377	if (base == 16 && str[0] == '0' && (str[1] == 'x' || str[1] == 'X'))
378		str += 2;
379	z = alloclongobject(0);
380	for ( ; z != NULL; ++str) {
381		int k = -1;
382		longobject *temp;
383
384		if (*str <= '9')
385			k = *str - '0';
386		else if (*str >= 'a')
387			k = *str - 'a' + 10;
388		else if (*str >= 'A')
389			k = *str - 'A' + 10;
390		if (k < 0 || k >= base)
391			break;
392		temp = muladd1(z, (digit)base, (digit)k);
393		DECREF(z);
394		z = temp;
395	}
396	if (sign < 0 && z != NULL && z->ob_size != 0)
397		z->ob_size = -(z->ob_size);
398	return (object *) z;
399}
400
401static longobject *x_divrem PROTO((longobject *, longobject *, longobject **));
402static object *long_pos PROTO((longobject *));
403static long_divrem PROTO((longobject *, longobject *,
404	longobject **, longobject **));
405
406/* Long division with remainder, top-level routine */
407
408static int
409long_divrem(a, b, pdiv, prem)
410	longobject *a, *b;
411	longobject **pdiv;
412	longobject **prem;
413{
414	int size_a = ABS(a->ob_size), size_b = ABS(b->ob_size);
415	longobject *z;
416
417	if (size_b == 0) {
418		err_setstr(ZeroDivisionError, "long division or modulo");
419		return -1;
420	}
421	if (size_a < size_b ||
422			size_a == size_b &&
423			a->ob_digit[size_a-1] < b->ob_digit[size_b-1]) {
424		/* |a| < |b|. */
425		*pdiv = alloclongobject(0);
426		INCREF(a);
427		*prem = (longobject *) a;
428		return 0;
429	}
430	if (size_b == 1) {
431		digit rem = 0;
432		z = divrem1(a, b->ob_digit[0], &rem);
433		if (z == NULL)
434			return -1;
435		*prem = (longobject *) newlongobject((long)rem);
436	}
437	else {
438		z = x_divrem(a, b, prem);
439		if (z == NULL)
440			return -1;
441	}
442	/* Set the signs.
443	   The quotient z has the sign of a*b;
444	   the remainder r has the sign of a,
445	   so a = b*z + r. */
446	if ((a->ob_size < 0) != (b->ob_size < 0))
447		z->ob_size = -(z->ob_size);
448	if (a->ob_size < 0 && (*prem)->ob_size != 0)
449		(*prem)->ob_size = -((*prem)->ob_size);
450	*pdiv = z;
451	return 0;
452}
453
454/* Unsigned long division with remainder -- the algorithm */
455
456static longobject *
457x_divrem(v1, w1, prem)
458	longobject *v1, *w1;
459	longobject **prem;
460{
461	int size_v = ABS(v1->ob_size), size_w = ABS(w1->ob_size);
462	digit d = (twodigits)BASE / (w1->ob_digit[size_w-1] + 1);
463	longobject *v = mul1(v1, d);
464	longobject *w = mul1(w1, d);
465	longobject *a;
466	int j, k;
467
468	if (v == NULL || w == NULL) {
469		XDECREF(v);
470		XDECREF(w);
471		return NULL;
472	}
473
474	assert(size_v >= size_w && size_w > 1); /* Assert checks by div() */
475	assert(v->ob_refcnt == 1); /* Since v will be used as accumulator! */
476	assert(size_w == ABS(w->ob_size)); /* That's how d was calculated */
477
478	size_v = ABS(v->ob_size);
479	a = alloclongobject(size_v - size_w + 1);
480
481	for (j = size_v, k = a->ob_size-1; a != NULL && k >= 0; --j, --k) {
482		digit vj = (j >= size_v) ? 0 : v->ob_digit[j];
483		twodigits q;
484		stwodigits carry = 0;
485		int i;
486
487		INTRCHECK({
488			DECREF(a);
489			a = NULL;
490			err_set(KeyboardInterrupt);
491			break;
492		})
493		if (vj == w->ob_digit[size_w-1])
494			q = MASK;
495		else
496			q = (((twodigits)vj << SHIFT) + v->ob_digit[j-1]) /
497				w->ob_digit[size_w-1];
498
499		while (w->ob_digit[size_w-2]*q >
500				((
501					((twodigits)vj << SHIFT)
502					+ v->ob_digit[j-1]
503					- q*w->ob_digit[size_w-1]
504								) << SHIFT)
505				+ v->ob_digit[j-2])
506			--q;
507
508		for (i = 0; i < size_w && i+k < size_v; ++i) {
509			twodigits z = w->ob_digit[i] * q;
510			digit zz = z >> SHIFT;
511			carry += v->ob_digit[i+k] - z + ((twodigits)zz << SHIFT);
512			v->ob_digit[i+k] = carry & MASK;
513			carry = (carry >> SHIFT) - zz;
514		}
515
516		if (i+k < size_v) {
517			carry += v->ob_digit[i+k];
518			v->ob_digit[i+k] = 0;
519		}
520
521		if (carry == 0)
522			a->ob_digit[k] = q;
523		else {
524			assert(carry == -1);
525			a->ob_digit[k] = q-1;
526			carry = 0;
527			for (i = 0; i < size_w && i+k < size_v; ++i) {
528				carry += v->ob_digit[i+k] + w->ob_digit[i];
529				v->ob_digit[i+k] = carry & MASK;
530				carry >>= SHIFT;
531			}
532		}
533	} /* for j, k */
534
535	if (a != NULL) {
536		a = long_normalize(a);
537		*prem = divrem1(v, d, &d);
538		/* d receives the (unused) remainder */
539		if (*prem == NULL) {
540			DECREF(a);
541			a = NULL;
542		}
543	}
544	DECREF(v);
545	DECREF(w);
546	return a;
547}
548
549/* Methods */
550
551/* Forward */
552static void long_dealloc PROTO((object *));
553static int long_print PROTO((object *, FILE *, int));
554static object *long_repr PROTO((object *));
555static int long_compare PROTO((longobject *, longobject *));
556
557static object *long_add PROTO((longobject *, longobject *));
558static object *long_sub PROTO((longobject *, longobject *));
559static object *long_mul PROTO((longobject *, longobject *));
560static object *long_div PROTO((longobject *, longobject *));
561static object *long_mod PROTO((longobject *, longobject *));
562static object *long_divmod PROTO((longobject *, longobject *));
563static object *long_pow PROTO((longobject *, longobject *));
564static object *long_neg PROTO((longobject *));
565static object *long_pos PROTO((longobject *));
566static object *long_abs PROTO((longobject *));
567static int long_nonzero PROTO((longobject *));
568static object *long_invert PROTO((longobject *));
569static object *long_lshift PROTO((longobject *, longobject *));
570static object *long_rshift PROTO((longobject *, longobject *));
571static object *long_and PROTO((longobject *, longobject *));
572static object *long_xor PROTO((longobject *, longobject *));
573static object *long_or PROTO((longobject *, longobject *));
574
575static void
576long_dealloc(v)
577	object *v;
578{
579	DEL(v);
580}
581
582/* ARGSUSED */
583static int
584long_print(v, fp, flags)
585	object *v;
586	FILE *fp;
587	int flags; /* Not used but required by interface */
588{
589	stringobject *str = (stringobject *) long_format(v, 10);
590	if (str == NULL)
591		return -1;
592	fprintf(fp, "%s", GETSTRINGVALUE(str));
593	DECREF(str);
594	return 0;
595}
596
597static object *
598long_repr(v)
599	object *v;
600{
601	return long_format(v, 10);
602}
603
604static int
605long_compare(a, b)
606	longobject *a, *b;
607{
608	int sign;
609
610	if (a->ob_size != b->ob_size) {
611		if (ABS(a->ob_size) == 0 && ABS(b->ob_size) == 0)
612			sign = 0;
613		else
614			sign = a->ob_size - b->ob_size;
615	}
616	else {
617		int i = ABS(a->ob_size);
618		while (--i >= 0 && a->ob_digit[i] == b->ob_digit[i])
619			;
620		if (i < 0)
621			sign = 0;
622		else
623			sign = (int)a->ob_digit[i] - (int)b->ob_digit[i];
624	}
625	return sign < 0 ? -1 : sign > 0 ? 1 : 0;
626}
627
628/* Add the absolute values of two long integers. */
629
630static longobject *x_add PROTO((longobject *, longobject *));
631static longobject *
632x_add(a, b)
633	longobject *a, *b;
634{
635	int size_a = ABS(a->ob_size), size_b = ABS(b->ob_size);
636	longobject *z;
637	int i;
638	digit carry = 0;
639
640	/* Ensure a is the larger of the two: */
641	if (size_a < size_b) {
642		{ longobject *temp = a; a = b; b = temp; }
643		{ int size_temp = size_a; size_a = size_b; size_b = size_temp; }
644	}
645	z = alloclongobject(size_a+1);
646	if (z == NULL)
647		return NULL;
648	for (i = 0; i < size_b; ++i) {
649		carry += a->ob_digit[i] + b->ob_digit[i];
650		z->ob_digit[i] = carry & MASK;
651		/* The following assumes unsigned shifts don't
652		   propagate the sign bit. */
653		carry >>= SHIFT;
654	}
655	for (; i < size_a; ++i) {
656		carry += a->ob_digit[i];
657		z->ob_digit[i] = carry & MASK;
658		carry >>= SHIFT;
659	}
660	z->ob_digit[i] = carry;
661	return long_normalize(z);
662}
663
664/* Subtract the absolute values of two integers. */
665
666static longobject *x_sub PROTO((longobject *, longobject *));
667static longobject *
668x_sub(a, b)
669	longobject *a, *b;
670{
671	int size_a = ABS(a->ob_size), size_b = ABS(b->ob_size);
672	longobject *z;
673	int i;
674	int sign = 1;
675	digit borrow = 0;
676
677	/* Ensure a is the larger of the two: */
678	if (size_a < size_b) {
679		sign = -1;
680		{ longobject *temp = a; a = b; b = temp; }
681		{ int size_temp = size_a; size_a = size_b; size_b = size_temp; }
682	}
683	else if (size_a == size_b) {
684		/* Find highest digit where a and b differ: */
685		i = size_a;
686		while (--i >= 0 && a->ob_digit[i] == b->ob_digit[i])
687			;
688		if (i < 0)
689			return alloclongobject(0);
690		if (a->ob_digit[i] < b->ob_digit[i]) {
691			sign = -1;
692			{ longobject *temp = a; a = b; b = temp; }
693		}
694		size_a = size_b = i+1;
695	}
696	z = alloclongobject(size_a);
697	if (z == NULL)
698		return NULL;
699	for (i = 0; i < size_b; ++i) {
700		/* The following assumes unsigned arithmetic
701		   works module 2**N for some N>SHIFT. */
702		borrow = a->ob_digit[i] - b->ob_digit[i] - borrow;
703		z->ob_digit[i] = borrow & MASK;
704		borrow >>= SHIFT;
705		borrow &= 1; /* Keep only one sign bit */
706	}
707	for (; i < size_a; ++i) {
708		borrow = a->ob_digit[i] - borrow;
709		z->ob_digit[i] = borrow & MASK;
710		borrow >>= SHIFT;
711	}
712	assert(borrow == 0);
713	if (sign < 0)
714		z->ob_size = -(z->ob_size);
715	return long_normalize(z);
716}
717
718static object *
719long_add(a, b)
720	longobject *a;
721	longobject *b;
722{
723	longobject *z;
724
725	if (a->ob_size < 0) {
726		if (b->ob_size < 0) {
727			z = x_add(a, b);
728			if (z != NULL && z->ob_size != 0)
729				z->ob_size = -(z->ob_size);
730		}
731		else
732			z = x_sub(b, a);
733	}
734	else {
735		if (b->ob_size < 0)
736			z = x_sub(a, b);
737		else
738			z = x_add(a, b);
739	}
740	return (object *)z;
741}
742
743static object *
744long_sub(a, b)
745	longobject *a;
746	longobject *b;
747{
748	longobject *z;
749
750	if (a->ob_size < 0) {
751		if (b->ob_size < 0)
752			z = x_sub(a, b);
753		else
754			z = x_add(a, b);
755		if (z != NULL && z->ob_size != 0)
756			z->ob_size = -(z->ob_size);
757	}
758	else {
759		if (b->ob_size < 0)
760			z = x_add(a, b);
761		else
762			z = x_sub(a, b);
763	}
764	return (object *)z;
765}
766
767static object *
768long_mul(a, b)
769	longobject *a;
770	longobject *b;
771{
772	int size_a;
773	int size_b;
774	longobject *z;
775	int i;
776
777	size_a = ABS(a->ob_size);
778	size_b = ABS(b->ob_size);
779	z = alloclongobject(size_a + size_b);
780	if (z == NULL)
781		return NULL;
782	for (i = 0; i < z->ob_size; ++i)
783		z->ob_digit[i] = 0;
784	for (i = 0; i < size_a; ++i) {
785		twodigits carry = 0;
786		twodigits f = a->ob_digit[i];
787		int j;
788
789		INTRCHECK({
790			DECREF(z);
791			err_set(KeyboardInterrupt);
792			return NULL;
793		})
794		for (j = 0; j < size_b; ++j) {
795			carry += z->ob_digit[i+j] + b->ob_digit[j] * f;
796			z->ob_digit[i+j] = carry & MASK;
797			carry >>= SHIFT;
798		}
799		for (; carry != 0; ++j) {
800			assert(i+j < z->ob_size);
801			carry += z->ob_digit[i+j];
802			z->ob_digit[i+j] = carry & MASK;
803			carry >>= SHIFT;
804		}
805	}
806	if (a->ob_size < 0)
807		z->ob_size = -(z->ob_size);
808	if (b->ob_size < 0)
809		z->ob_size = -(z->ob_size);
810	return (object *) long_normalize(z);
811}
812
813/* The / and % operators are now defined in terms of divmod().
814   The expression a mod b has the value a - b*floor(a/b).
815   The long_divrem function gives the remainder after division of
816   |a| by |b|, with the sign of a.  This is also expressed
817   as a - b*trunc(a/b), if trunc truncates towards zero.
818   Some examples:
819   	 a	 b	a rem b		a mod b
820   	 13	 10	 3		 3
821   	-13	 10	-3		 7
822   	 13	-10	 3		-7
823   	-13	-10	-3		-3
824   So, to get from rem to mod, we have to add b if a and b
825   have different signs.  We then subtract one from the 'div'
826   part of the outcome to keep the invariant intact. */
827
828static int l_divmod PROTO((longobject *, longobject *,
829	longobject **, longobject **));
830static int
831l_divmod(v, w, pdiv, pmod)
832	longobject *v;
833	longobject *w;
834	longobject **pdiv;
835	longobject **pmod;
836{
837	longobject *div, *mod;
838
839	if (long_divrem(v, w, &div, &mod) < 0)
840		return -1;
841	if (mod->ob_size < 0 && w->ob_size > 0 ||
842				mod->ob_size > 0 && w->ob_size < 0) {
843		longobject *temp;
844		longobject *one;
845		temp = (longobject *) long_add(mod, w);
846		DECREF(mod);
847		mod = temp;
848		if (mod == NULL) {
849			DECREF(div);
850			return -1;
851		}
852		one = (longobject *) newlongobject(1L);
853		if (one == NULL ||
854		    (temp = (longobject *) long_sub(div, one)) == NULL) {
855			DECREF(mod);
856			DECREF(div);
857			return -1;
858		}
859		DECREF(div);
860		div = temp;
861	}
862	*pdiv = div;
863	*pmod = mod;
864	return 0;
865}
866
867static object *
868long_div(v, w)
869	longobject *v;
870	longobject *w;
871{
872	longobject *div, *mod;
873	if (l_divmod(v, w, &div, &mod) < 0)
874		return NULL;
875	DECREF(mod);
876	return (object *)div;
877}
878
879static object *
880long_mod(v, w)
881	longobject *v;
882	longobject *w;
883{
884	longobject *div, *mod;
885	if (l_divmod(v, w, &div, &mod) < 0)
886		return NULL;
887	DECREF(div);
888	return (object *)mod;
889}
890
891static object *
892long_divmod(v, w)
893	longobject *v;
894	longobject *w;
895{
896	object *z;
897	longobject *div, *mod;
898	if (l_divmod(v, w, &div, &mod) < 0)
899		return NULL;
900	z = newtupleobject(2);
901	if (z != NULL) {
902		settupleitem(z, 0, (object *) div);
903		settupleitem(z, 1, (object *) mod);
904	}
905	else {
906		DECREF(div);
907		DECREF(mod);
908	}
909	return z;
910}
911
912static object *
913long_pow(a, b)
914	longobject *a;
915	longobject *b;
916{
917	longobject *z;
918	int size_b, i;
919
920	size_b = b->ob_size;
921	if (size_b < 0) {
922		err_setstr(ValueError, "long integer to the negative power");
923		return NULL;
924	}
925
926	z = (longobject *)newlongobject(1L);
927
928	INCREF(a);
929	for (i = 0; i < size_b; ++i) {
930		digit bi = b->ob_digit[i];
931		int j;
932
933		for (j = 0; j < SHIFT; ++j) {
934			longobject *temp;
935
936			if (bi & 1) {
937				temp = (longobject *)long_mul(z, a);
938				DECREF(z);
939				z = temp;
940				if (z == NULL)
941					break;
942			}
943			bi >>= 1;
944			if (bi == 0 && i+1 == size_b)
945				break;
946			temp = (longobject *)long_mul(a, a);
947			DECREF(a);
948			a = temp;
949			if (a == NULL) {
950				DECREF(z);
951				z = NULL;
952				break;
953			}
954		}
955		if (a == NULL)
956			break;
957	}
958	XDECREF(a);
959	return (object *)z;
960}
961
962static object *
963long_invert(v)
964	longobject *v;
965{
966	/* Implement ~x as -(x+1) */
967	longobject *x;
968	longobject *w;
969	w = (longobject *)newlongobject(1L);
970	if (w == NULL)
971		return NULL;
972	x = (longobject *) long_add(v, w);
973	DECREF(w);
974	if (x == NULL)
975		return NULL;
976	if (x->ob_size != 0)
977		x->ob_size = -(x->ob_size);
978	return (object *)x;
979}
980
981static object *
982long_pos(v)
983	longobject *v;
984{
985	INCREF(v);
986	return (object *)v;
987}
988
989static object *
990long_neg(v)
991	longobject *v;
992{
993	longobject *z;
994	int i, n;
995	n = ABS(v->ob_size);
996	if (n == 0) {
997		/* -0 == 0 */
998		INCREF(v);
999		return (object *) v;
1000	}
1001	z = alloclongobject(ABS(n));
1002	if (z == NULL)
1003		return NULL;
1004	for (i = 0; i < n; i++)
1005		z->ob_digit[i] = v->ob_digit[i];
1006	z->ob_size = -(v->ob_size);
1007	return (object *)z;
1008}
1009
1010static object *
1011long_abs(v)
1012	longobject *v;
1013{
1014	if (v->ob_size < 0)
1015		return long_neg(v);
1016	else {
1017		INCREF(v);
1018		return (object *)v;
1019	}
1020}
1021
1022static int
1023long_nonzero(v)
1024	longobject *v;
1025{
1026	return ABS(v->ob_size) != 0;
1027}
1028
1029static object *
1030long_rshift(a, b)
1031	longobject *a;
1032	longobject *b;
1033{
1034	longobject *z;
1035	long shiftby;
1036	int newsize, wordshift, loshift, hishift, i, j;
1037	digit lomask, himask;
1038
1039	if (a->ob_size < 0) {
1040		/* Right shifting negative numbers is harder */
1041		longobject *a1, *a2, *a3;
1042		a1 = (longobject *) long_invert(a);
1043		if (a1 == NULL) return NULL;
1044		a2 = (longobject *) long_rshift(a1, b);
1045		DECREF(a1);
1046		if (a2 == NULL) return NULL;
1047		a3 = (longobject *) long_invert(a2);
1048		DECREF(a2);
1049		return (object *) a3;
1050	}
1051
1052	shiftby = getlongvalue((object *)b);
1053	if (shiftby == -1L && err_occurred())
1054		return NULL;
1055	if (shiftby < 0) {
1056		err_setstr(ValueError, "negative shift count");
1057		return NULL;
1058	}
1059	wordshift = shiftby / SHIFT;
1060	newsize = ABS(a->ob_size) - wordshift;
1061	if (newsize <= 0) {
1062		z = alloclongobject(0);
1063		return (object *)z;
1064	}
1065	loshift = shiftby % SHIFT;
1066	hishift = SHIFT - loshift;
1067	lomask = ((digit)1 << hishift) - 1;
1068	himask = MASK ^ lomask;
1069	z = alloclongobject(newsize);
1070	if (z == NULL)
1071		return NULL;
1072	if (a->ob_size < 0)
1073		z->ob_size = -(z->ob_size);
1074	for (i = 0, j = wordshift; i < newsize; i++, j++) {
1075		z->ob_digit[i] = (a->ob_digit[j] >> loshift) & lomask;
1076		if (i+1 < newsize)
1077			z->ob_digit[i] |=
1078			  (a->ob_digit[j+1] << hishift) & himask;
1079	}
1080	return (object *) long_normalize(z);
1081}
1082
1083static object *
1084long_lshift(a, b)
1085	longobject *a;
1086	longobject *b;
1087{
1088	longobject *z;
1089	long shiftby;
1090	int newsize, wordshift, loshift, hishift, i, j;
1091	digit lomask, himask;
1092
1093	shiftby = getlongvalue((object *)b);
1094	if (shiftby == -1L && err_occurred())
1095		return NULL;
1096	if (shiftby < 0) {
1097		err_setstr(ValueError, "negative shift count");
1098		return NULL;
1099	}
1100	if (shiftby > MASK) {
1101		err_setstr(ValueError, "outrageous left shift count");
1102		return NULL;
1103	}
1104	if (shiftby % SHIFT == 0) {
1105		wordshift = shiftby / SHIFT;
1106		loshift = 0;
1107		hishift = SHIFT;
1108		newsize = ABS(a->ob_size) + wordshift;
1109		lomask = MASK;
1110		himask = 0;
1111	}
1112	else {
1113		wordshift = shiftby / SHIFT + 1;
1114		loshift = SHIFT - shiftby%SHIFT;
1115		hishift = shiftby % SHIFT;
1116		newsize = ABS(a->ob_size) + wordshift;
1117		lomask = ((digit)1 << hishift) - 1;
1118		himask = MASK ^ lomask;
1119	}
1120	z = alloclongobject(newsize);
1121	if (z == NULL)
1122		return NULL;
1123	if (a->ob_size < 0)
1124		z->ob_size = -(z->ob_size);
1125	for (i = 0; i < wordshift; i++)
1126		z->ob_digit[i] = 0;
1127	for (i = wordshift, j = 0; i < newsize; i++, j++) {
1128		if (i > 0)
1129			z->ob_digit[i-1] |=
1130				(a->ob_digit[j] << hishift) & himask;
1131		z->ob_digit[i] =
1132			(a->ob_digit[j] >> loshift) & lomask;
1133	}
1134	return (object *) long_normalize(z);
1135}
1136
1137
1138/* Bitwise and/xor/or operations */
1139
1140#define MAX(x, y) ((x) < (y) ? (y) : (x))
1141#define MIN(x, y) ((x) > (y) ? (y) : (x))
1142
1143static object *long_bitwise PROTO((longobject *, int, longobject *));
1144static object *
1145long_bitwise(a, op, b)
1146	longobject *a;
1147	int op; /* '&', '|', '^' */
1148	longobject *b;
1149{
1150	digit maska, maskb; /* 0 or MASK */
1151	int negz;
1152	int size_a, size_b, size_z;
1153	longobject *z;
1154	int i;
1155	digit diga, digb;
1156	object *v;
1157
1158	if (a->ob_size < 0) {
1159		a = (longobject *) long_invert(a);
1160		maska = MASK;
1161	}
1162	else {
1163		INCREF(a);
1164		maska = 0;
1165	}
1166	if (b->ob_size < 0) {
1167		b = (longobject *) long_invert(b);
1168		maskb = MASK;
1169	}
1170	else {
1171		INCREF(b);
1172		maskb = 0;
1173	}
1174
1175	size_a = a->ob_size;
1176	size_b = b->ob_size;
1177	size_z = MAX(size_a, size_b);
1178	z = alloclongobject(size_z);
1179	if (a == NULL || b == NULL || z == NULL) {
1180		XDECREF(a);
1181		XDECREF(b);
1182		XDECREF(z);
1183		return NULL;
1184	}
1185
1186	negz = 0;
1187	switch (op) {
1188	case '^':
1189		if (maska != maskb) {
1190			maska ^= MASK;
1191			negz = -1;
1192		}
1193		break;
1194	case '&':
1195		if (maska && maskb) {
1196			op = '|';
1197			maska ^= MASK;
1198			maskb ^= MASK;
1199			negz = -1;
1200		}
1201		break;
1202	case '|':
1203		if (maska || maskb) {
1204			op = '&';
1205			maska ^= MASK;
1206			maskb ^= MASK;
1207			negz = -1;
1208		}
1209		break;
1210	}
1211
1212	for (i = 0; i < size_z; ++i) {
1213		diga = (i < size_a ? a->ob_digit[i] : 0) ^ maska;
1214		digb = (i < size_b ? b->ob_digit[i] : 0) ^ maskb;
1215		switch (op) {
1216		case '&': z->ob_digit[i] = diga & digb; break;
1217		case '|': z->ob_digit[i] = diga | digb; break;
1218		case '^': z->ob_digit[i] = diga ^ digb; break;
1219		}
1220	}
1221
1222	DECREF(a);
1223	DECREF(b);
1224	z = long_normalize(z);
1225	if (negz == 0)
1226		return (object *) z;
1227	v = long_invert(z);
1228	DECREF(z);
1229	return v;
1230}
1231
1232static object *
1233long_and(a, b)
1234	longobject *a;
1235	longobject *b;
1236{
1237	return long_bitwise(a, '&', b);
1238}
1239
1240static object *
1241long_xor(a, b)
1242	longobject *a;
1243	longobject *b;
1244{
1245	return long_bitwise(a, '^', b);
1246}
1247
1248static object *
1249long_or(a, b)
1250	longobject *a;
1251	longobject *b;
1252{
1253	return long_bitwise(a, '|', b);
1254}
1255
1256#define UF (object* (*) FPROTO((object *))) /* Unary function */
1257#define BF (object* (*) FPROTO((object *, object *))) /* Binary function */
1258#define IF (int (*) FPROTO((object *))) /* Int function */
1259
1260static number_methods long_as_number = {
1261	BF long_add,	/*nb_add*/
1262	BF long_sub,	/*nb_subtract*/
1263	BF long_mul,	/*nb_multiply*/
1264	BF long_div,	/*nb_divide*/
1265	BF long_mod,	/*nb_remainder*/
1266	BF long_divmod,	/*nb_divmod*/
1267	BF long_pow,	/*nb_power*/
1268	UF long_neg,	/*nb_negative*/
1269	UF long_pos,	/*tp_positive*/
1270	UF long_abs,	/*tp_absolute*/
1271	IF long_nonzero,/*tp_nonzero*/
1272	UF long_invert,	/*nb_invert*/
1273	BF long_lshift,	/*nb_lshift*/
1274	BF long_rshift,	/*nb_rshift*/
1275	BF long_and,	/*nb_and*/
1276	BF long_xor,	/*nb_xor*/
1277	BF long_or,	/*nb_or*/
1278};
1279
1280typeobject Longtype = {
1281	OB_HEAD_INIT(&Typetype)
1282	0,
1283	"long int",
1284	sizeof(longobject) - sizeof(digit),
1285	sizeof(digit),
1286	long_dealloc,	/*tp_dealloc*/
1287	long_print,	/*tp_print*/
1288	0,		/*tp_getattr*/
1289	0,		/*tp_setattr*/
1290	(int (*) FPROTO((object *, object *)))
1291	long_compare,	/*tp_compare*/
1292	long_repr,	/*tp_repr*/
1293	&long_as_number,/*tp_as_number*/
1294	0,		/*tp_as_sequence*/
1295	0,		/*tp_as_mapping*/
1296};
1297