?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
improve treatment of multi-line replies, ignore empty lines
[python/dscho.git] / Modules / mpzmodule.c
blobce31a4bc4736261a435c394d42468ac67818b28e
1 /***********************************************************
2 Copyright 1991-1995 by Stichting Mathematisch Centrum, Amsterdam,
3 The Netherlands.
4
5 All Rights Reserved
6
7 Permission to use, copy, modify, and distribute this software and its
8 documentation for any purpose and without fee is hereby granted,
9 provided that the above copyright notice appear in all copies and that
10 both that copyright notice and this permission notice appear in
11 supporting documentation, and that the names of Stichting Mathematisch
12 Centrum or CWI not be used in advertising or publicity pertaining to
13 distribution of the software without specific, written prior permission.
14
15 STICHTING MATHEMATISCH CENTRUM DISCLAIMS ALL WARRANTIES WITH REGARD TO
16 THIS SOFTWARE, INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND
17 FITNESS, IN NO EVENT SHALL STICHTING MATHEMATISCH CENTRUM BE LIABLE
18 FOR ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
19 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
20 ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT
21 OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.
22
23 ******************************************************************/
24
25 /* MPZ module */
26
27 /* This module provides an interface to an alternate Multi-Precision
28 library, GNU MP in this case */
29
30 /* XXX note: everywhere where mpz_size is called,
31 sizeof (limb) == sizeof (long) has been assumed. */
32
33
34 /* MPZ objects */
35
36 #include "allobjects.h"
37
38 #include <assert.h>
39 #include <sys/types.h> /* For size_t */
40
41 /*
42 ** These are the cpp-flags used in this file...
43 **
44 **
45 ** MPZ_MDIV_BUG works around the mpz_m{div,mod,...} routines.
46 ** This bug has been fixed in a later release of
47 ** GMP.
48 **
49 ** MPZ_GET_STR_BUG mpz_get_str corrupts memory, seems to be fixed
50 ** in a later release
51 **
52 ** MPZ_DEBUG generates a bunch of diagnostic messages
53 **
54 ** MPZ_SPARE_MALLOC if set, results in extra code that tries to
55 ** minimize the creation of extra objects.
56 **
57 ** MPZ_TEST_DIV extra diagnostic output on stderr, when division
58 ** routines are involved
59 **
60 ** MPZ_LIB_DOES_CHECKING if set, assumes that mpz library doesn't call
61 ** alloca with arg < 0 (when casted to a signed
62 ** integral type).
63 **
64 ** MPZ_CONVERSIONS_AS_METHODS if set, presents the conversions as
65 ** methods. e.g., `mpz(5).long() == 5L'
66 ** Later, Guido provided an interface to the
67 ** standard functions. So this flag has no been
68 ** cleared, and `long(mpz(5)) == 5L'
69 **
70 ** MP_TEST_ALLOC If set, you would discover why MPZ_GET_STR_BUG
71 ** is needed
72 **
73 ** MAKEDUMMYINT Must be set if dynamic linking will be used
74 */
75
76
77 /*
78 ** IMHO, mpz_m{div,mod,divmod}() do the wrong things when the denominator < 0
79 ** I assume that this will be fixed in a future release
80 */
81 /*#define MPZ_MDIV_BUG fixed the (for me) nexessary parts in libgmp.a */
82 /*
83 ** IMO, mpz_get_str() assumes a bit too large target space, if he doesn't
84 ** allocate it himself
85 */
86 #define MPZ_GET_STR_BUG
87
88 #include "gmp.h"
89 typedef struct {
90 OB_HEAD
91 MP_INT mpz; /* the actual number */
92 } mpzobject;
93
94 staticforward typeobject MPZtype;
95
96 #define is_mpzobject(v) ((v)->ob_type == &MPZtype)
97
98 static const char initialiser_name[] = "mpz";
99
100 /* #define MPZ_DEBUG */
101
102 static mpzobject *
103 newmpzobject()
104 {
105 mpzobject *mpzp;
106
107
108 #ifdef MPZ_DEBUG
109 fputs( "mpz_object() called...\n", stderr );
110 #endif /* def MPZ_DEBUG */
111 mpzp = NEWOBJ(mpzobject, &MPZtype);
112 if (mpzp == NULL)
113 return NULL;
114
115 mpz_init(&mpzp->mpz); /* actual initialisation */
116 return mpzp;
117 } /* newmpzobject() */
118
119 #ifdef MPZ_GET_STR_BUG
120 #include "gmp-impl.h"
121 #include "longlong.h"
122 #endif /* def MPZ_GET_STR_BUG */
123
124 static object *
125 mpz_format(objp, base, withname)
126 object *objp;
127 int base;
128 unsigned char withname;
129 {
130 mpzobject *mpzp = (mpzobject *)objp;
131 stringobject *strobjp;
132 int i;
133 int cmpres;
134 int taglong;
135 char *cp;
136 char prefix[5], *tcp;
137
138
139 tcp = &prefix[0];
140
141 if (mpzp == NULL || !is_mpzobject(mpzp)) {
142 err_badcall();
143 return NULL;
144 }
145
146 assert(base >= 2 && base <= 36);
147
148 if (withname)
149 i = strlen(initialiser_name) + 2; /* e.g. 'mpz(' + ')' */
150 else
151 i = 0;
152
153 if ((cmpres = mpz_cmp_si(&mpzp->mpz, 0L)) == 0)
154 base = 10; /* '0' in every base, right */
155 else if (cmpres < 0) {
156 *tcp++ = '-';
157 i += 1; /* space to hold '-' */
158 }
159
160 #ifdef MPZ_DEBUG
161 fprintf(stderr, "mpz_format: mpz_sizeinbase %d\n",
162 (int)mpz_sizeinbase(&mpzp->mpz, base));
163 #endif /* def MPZ_DEBUG */
164 #ifdef MPZ_GET_STR_BUG
165 i += ((size_t) abs(mpzp->mpz.size) * BITS_PER_MP_LIMB
166 * __mp_bases[base].chars_per_bit_exactly) + 1;
167 #else /* def MPZ_GET_STR_BUG */
168 i += (int)mpz_sizeinbase(&mpzp->mpz, base);
169 #endif /* def MPZ_GET_STR_BUG else */
170
171 if (base == 16) {
172 *tcp++ = '0';
173 *tcp++ = 'x';
174 i += 2; /* space to hold '0x' */
175 }
176 else if (base == 8) {
177 *tcp++ = '0';
178 i += 1; /* space to hold the extra '0' */
179 }
180 else if (base > 10) {
181 *tcp++ = '0' + base / 10;
182 *tcp++ = '0' + base % 10;
183 *tcp++ = '#';
184 i += 3; /* space to hold e.g. '12#' */
185 }
186 else if (base < 10) {
187 *tcp++ = '0' + base;
188 *tcp++ = '#';
189 i += 2; /* space to hold e.g. '6#' */
190 }
191
192 /*
193 ** the following code looks if we need a 'L' attached to the number
194 ** it will also attach an 'L' to the value -0x80000000
195 */
196 taglong = 0;
197 if (mpz_size(&mpzp->mpz) > 1
198 || (long)mpz_get_ui(&mpzp->mpz) < 0L) {
199 taglong = 1;
200 i += 1; /* space to hold 'L' */
201 }
202
203 #ifdef MPZ_DEBUG
204 fprintf(stderr, "mpz_format: requesting string size %d\n", i);
205 #endif /* def MPZ_DEBUG */
206 if ((strobjp = (stringobject *)newsizedstringobject((char *)0, i))
207 == NULL)
208 return NULL;
209
210 /* get the beginning of the string memory and start copying things */
211 cp = GETSTRINGVALUE(strobjp);
212 if (withname) {
213 strcpy(cp, initialiser_name);
214 cp += strlen(initialiser_name);
215 *cp++ = '('; /*')'*/
216 }
217
218 /* copy the already prepared prefix; e.g. sign and base indicator */
219 *tcp = '\0';
220 strcpy(cp, prefix);
221 cp += tcp - prefix;
222
223 /* since' we have the sign already, let the lib think it's a positive
224 number */
225 if (cmpres < 0)
226 mpz_neg(&mpzp->mpz,&mpzp->mpz); /* hack Hack HAck HACk HACK */
227 (void)mpz_get_str(cp, base, &mpzp->mpz);
228 if (cmpres < 0)
229 mpz_neg(&mpzp->mpz,&mpzp->mpz); /* hack Hack HAck HACk HACK */
230 #ifdef MPZ_DEBUG
231 fprintf(stderr, "mpz_format: base (ultim) %d, mpz_get_str: %s\n",
232 base, cp);
233 #endif /* def MPZ_DEBUG */
234 cp += strlen(cp);
235
236 if (taglong)
237 *cp++ = 'L';
238 if (withname)
239 *cp++ = /*'('*/ ')';
240
241 *cp = '\0';
242
243 #ifdef MPZ_DEBUG
244 fprintf(stderr,
245 "mpz_format: cp (str end) 0x%x, begin 0x%x, diff %d, i %d\n",
246 cp, GETSTRINGVALUE(strobjp), cp - GETSTRINGVALUE(strobjp), i);
247 #endif /* def MPZ_DEBUG */
248 assert(cp - GETSTRINGVALUE(strobjp) <= i);
249
250 if (cp - GETSTRINGVALUE(strobjp) != i) {
251 strobjp->ob_size -= i - (cp - GETSTRINGVALUE(strobjp));
252 }
253
254 return (object *)strobjp;
255 } /* mpz_format() */
256
257 /* MPZ methods */
258
259 static void
260 mpz_dealloc(mpzp)
261 mpzobject *mpzp;
262 {
263 #ifdef MPZ_DEBUG
264 fputs( "mpz_dealloc() called...\n", stderr );
265 #endif /* def MPZ_DEBUG */
266 mpz_clear(&mpzp->mpz);
267 DEL(mpzp);
268 } /* mpz_dealloc() */
269
270
271 /* pointers to frequently used values 0, 1 and -1 */
272 static mpzobject *mpz_value_zero, *mpz_value_one, *mpz_value_mone;
273
274 static int
275 mpz_compare(a, b)
276 mpzobject *a, *b;
277 {
278 int cmpres;
279
280
281 /* guido sez it's better to return -1, 0 or 1 */
282 return (cmpres = mpz_cmp( &a->mpz, &b->mpz )) == 0 ? 0
283 : cmpres > 0 ? 1 : -1;
284 } /* mpz_compare() */
285
286 static object *
287 mpz_addition(a, b)
288 mpzobject *a;
289 mpzobject *b;
290 {
291 mpzobject *z;
292
293
294 #ifdef MPZ_SPARE_MALLOC
295 if (mpz_cmp_ui(&a->mpz, (unsigned long int)0) == 0) {
296 INCREF(b);
297 return (object *)b;
298 }
299
300 if (mpz_cmp_ui(&b->mpz, (unsigned long int)0) == 0) {
301 INCREF(a);
302 return (object *)a;
303 }
304 #endif /* def MPZ_SPARE_MALLOC */
305
306 if ((z = newmpzobject()) == NULL)
307 return NULL;
308
309 mpz_add(&z->mpz, &a->mpz, &b->mpz);
310 return (object *)z;
311 } /* mpz_addition() */
312
313 static object *
314 mpz_substract(a, b)
315 mpzobject *a;
316 mpzobject *b;
317 {
318 mpzobject *z;
319
320
321 #ifdef MPZ_SPARE_MALLOC
322 if (mpz_cmp_ui(&b->mpz, (unsigned long int)0) == 0) {
323 INCREF(a);
324 return (object *)a;
325 }
326 #endif /* MPZ_SPARE_MALLOC */
327
328 if ((z = newmpzobject()) == NULL)
329 return NULL;
330
331 mpz_sub(&z->mpz, &a->mpz, &b->mpz);
332 return (object *)z;
333 } /* mpz_substract() */
334
335 static object *
336 mpz_multiply(a, b)
337 mpzobject *a;
338 mpzobject *b;
339 {
340 #ifdef MPZ_SPARE_MALLOC
341 int cmpres;
342 #endif /* def MPZ_SPARE_MALLOC */
343 mpzobject *z;
344
345
346 #ifdef MPZ_SPARE_MALLOC
347 if ((cmpres = mpz_cmp_ui(&a->mpz, (unsigned long int)0)) == 0) {
348 INCREF(mpz_value_zero);
349 return (object *)mpz_value_zero;
350 }
351 if (cmpres > 0 && mpz_cmp_ui(&a->mpz, (unsigned long int)1) == 0) {
352 INCREF(b);
353 return (object *)b;
354 }
355
356 if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long_int)0)) == 0) {
357 INCREF(mpz_value_zero);
358 return (object *)mpz_value_zero;
359 }
360 if (cmpres > 0 && mpz_cmp_ui(&b->mpz, (unsigned long int)1) == 0) {
361 INCREF(a);
362 return (object *)a;
363 }
364 #endif /* MPZ_SPARE_MALLOC */
365
366 if ((z = newmpzobject()) == NULL)
367 return NULL;
368
369 mpz_mul( &z->mpz, &a->mpz, &b->mpz );
370 return (object *)z;
371
372 } /* mpz_multiply() */
373
374 static object *
375 mpz_divide(a, b)
376 mpzobject *a;
377 mpzobject *b;
378 {
379 #ifdef MPZ_SPARE_MALLOC
380 int cmpres;
381 #endif /* def MPZ_SPARE_MALLOC */
382 mpzobject *z;
383
384
385 if ((
386 #ifdef MPZ_SPARE_MALLOC
387 cmpres =
388 #endif /* def MPZ_SPARE_MALLOC */
389 mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
390 err_setstr(ZeroDivisionError, "mpz./ by zero");
391 return NULL;
392 }
393 #ifdef MPZ_SPARE_MALLOC
394 if (cmpres > 0 && mpz_cmp_ui(&b->mpz(unsigned long int)1) == 0) {
395 INCREF(a);
396 return (object *)a;
397 }
398 #endif /* def MPZ_SPARE_MALLOC */
399
400 if ((z = newmpzobject()) == NULL)
401 return NULL;
402
403 #ifdef MPZ_TEST_DIV
404 fputs("mpz_divide: div result", stderr);
405 mpz_div(&z->mpz, &a->mpz, &b->mpz);
406 mpz_out_str(stderr, 10, &z->mpz);
407 putc('\n', stderr);
408 #endif /* def MPZ_TEST_DIV */
409 #ifdef MPZ_MDIV_BUG
410 if ((mpz_cmp_ui(&a->mpz, (unsigned long int)0) < 0)
411 != (mpz_cmp_ui(&b->mpz, (unsigned long int)0) < 0)) {
412 /*
413 ** numerator has other sign than denominator: we have
414 ** to look at the remainder for a correction, since mpz_mdiv
415 ** also calls mpz_divmod, I can as well do it myself
416 */
417 MP_INT tmpmpz;
418
419
420 mpz_init(&tmpmpz);
421 mpz_divmod(&z->mpz, &tmpmpz, &a->mpz, &b->mpz);
422
423 if (mpz_cmp_ui(&tmpmpz, (unsigned long int)0) != 0)
424 mpz_sub_ui(&z->mpz, &z->mpz, (unsigned long int)1);
425
426 mpz_clear(&tmpmpz);
427 }
428 else
429 mpz_div(&z->mpz, &a->mpz, &b->mpz);
430 /* the ``naive'' implementation does it right for operands
431 having the same sign */
432
433 #else /* def MPZ_MDIV_BUG */
434 mpz_mdiv(&z->mpz, &a->mpz, &b->mpz);
435 #endif /* def MPZ_MDIV_BUG else */
436 #ifdef MPZ_TEST_DIV
437 fputs("mpz_divide: mdiv result", stderr);
438 mpz_out_str(stderr, 10, &z->mpz);
439 putc('\n', stderr);
440 #endif /* def MPZ_TEST_DIV */
441 return (object *)z;
442
443 } /* mpz_divide() */
444
445 static object *
446 mpz_remainder(a, b)
447 mpzobject *a;
448 mpzobject *b;
449 {
450 #ifdef MPZ_SPARE_MALLOC
451 int cmpres;
452 #endif /* def MPZ_SPARE_MALLOC */
453 mpzobject *z;
454
455
456 if ((
457 #ifdef MPZ_SPARE_MALLOC
458 cmpres =
459 #endif /* def MPZ_SPARE_MALLOC */
460 mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
461 err_setstr(ZeroDivisionError, "mpz.% by zero");
462 return NULL;
463 }
464 #ifdef MPZ_SPARE_MALLOC
465 if (cmpres > 0) {
466 if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)2)) == 0) {
467 INCREF(mpz_value_one);
468 return (object *)mpz_value_one;
469 }
470 if (cmpres < 0) {
471 /* b must be 1 now */
472 INCREF(mpz_value_zero);
473 return (object *)mpz_value_zero;
474 }
475 }
476 #endif /* def MPZ_SPARE_MALLOC */
477
478 if ((z = newmpzobject()) == NULL)
479 return NULL;
480
481 #ifdef MPZ_TEST_DIV
482 fputs("mpz_remain: mod result", stderr);
483 mpz_mod(&z->mpz, &a->mpz, &b->mpz);
484 mpz_out_str(stderr, 10, &z->mpz);
485 putc('\n', stderr);
486 #endif /* def MPZ_TEST_DIV */
487 #ifdef MPZ_MDIV_BUG
488
489 /* the ``naive'' implementation does it right for operands
490 having the same sign */
491 mpz_mod(&z->mpz, &a->mpz, &b->mpz);
492
493 /* assumption: z, a and b all point to different locations */
494 if ((mpz_cmp_ui(&a->mpz, (unsigned long int)0) < 0)
495 != (mpz_cmp_ui(&b->mpz, (unsigned long int)0) < 0)
496 && mpz_cmp_ui(&z->mpz, (unsigned long int)0) != 0)
497 mpz_add(&z->mpz, &z->mpz, &b->mpz);
498 /*
499 ** numerator has other sign than denominator: we have
500 ** to look at the remainder for a correction, since mpz_mdiv
501 ** also calls mpz_divmod, I can as well do it myself
502 */
503 #else /* def MPZ_MDIV_BUG */
504 mpz_mmod(&z->mpz, &a->mpz, &b->mpz);
505 #endif /* def MPZ_MDIV_BUG else */
506 #ifdef MPZ_TEST_DIV
507 fputs("mpz_remain: mmod result", stderr);
508 mpz_out_str(stderr, 10, &z->mpz);
509 putc('\n', stderr);
510 #endif /* def MPZ_TEST_DIV */
511 return (object *)z;
512
513 } /* mpz_remainder() */
514
515 static object *
516 mpz_div_and_mod(a, b)
517 mpzobject *a;
518 mpzobject *b;
519 {
520 object *z = NULL;
521 mpzobject *x = NULL, *y = NULL;
522
523
524 if (mpz_cmp_ui(&b->mpz, (unsigned long int)0) == 0) {
525 err_setstr(ZeroDivisionError, "mpz.divmod by zero");
526 return NULL;
527 }
528
529 if ((z = newtupleobject(2)) == NULL
530 || (x = newmpzobject()) == NULL
531 || (y = newmpzobject()) == NULL) {
532 XDECREF(z);
533 XDECREF(x);
534 XDECREF(y);
535 return NULL;
536 }
537
538 #ifdef MPZ_TEST_DIV
539 fputs("mpz_divmod: dm result", stderr);
540 mpz_divmod(&x->mpz, &y->mpz, &a->mpz, &b->mpz);
541 mpz_out_str(stderr, 10, &x->mpz);
542 putc('\n', stderr);
543 mpz_out_str(stderr, 10, &y->mpz);
544 putc('\n', stderr);
545 #endif /* def MPZ_TEST_DIV */
546 #ifdef MPZ_MDIV_BUG
547 mpz_divmod(&x->mpz, &y->mpz, &a->mpz, &b->mpz);
548 if ((mpz_cmp_ui(&a->mpz, (unsigned long int)0) < 0)
549 != (mpz_cmp_ui(&b->mpz, (unsigned long int)0) < 0)
550 && mpz_cmp_ui(&y->mpz, (unsigned long int)0) != 0) {
551 /*
552 ** numerator has other sign than denominator: we have
553 ** to look at the remainder for a correction.
554 */
555 mpz_add(&y->mpz, &y->mpz, &b->mpz);
556 mpz_sub_ui(&x->mpz, &x->mpz, (unsigned long int)1);
557 }
558 #else /* def MPZ_MDIV_BUG */
559 mpz_mdivmod( &x->mpz, &y->mpz, &a->mpz, &b->mpz );
560 #endif /* def MPZ_MDIV_BUG else */
561 #ifdef MPZ_TEST_DIV
562 fputs("mpz_divmod: mdm result", stderr);
563 mpz_out_str(stderr, 10, &x->mpz);
564 putc('\n', stderr);
565 mpz_out_str(stderr, 10, &y->mpz);
566 putc('\n', stderr);
567 #endif /* def MPZ_TEST_DIV */
568
569 (void)settupleitem(z, 0, (object *)x);
570 (void)settupleitem(z, 1, (object *)y);
571
572 return z;
573 } /* mpz_div_and_mod() */
574
575 static object *
576 mpz_power(a, b, m)
577 mpzobject *a;
578 mpzobject *b;
579 mpzobject *m;
580 {
581 mpzobject *z;
582 int cmpres;
583 long int longtmp1, longtmp2;
584
585 if ((object *)m!=Py_None)
586 {
587 mpzobject *z2;
588 INCREF(Py_None);
589 z=mpz_power(a, b, (mpzobject *)Py_None);
590 DECREF(Py_None);
591 if (z==NULL) return(z);
592 z2=mpz_remainder(z, m);
593 DECREF(z);
594 return((object *)z2);
595 }
596
597 if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
598 /* the gnu-mp lib sets pow(0,0) to 0, we to 1 */
599
600 INCREF(mpz_value_one);
601 return (object *)mpz_value_one;
602 }
603
604 if (cmpres < 0) {
605 err_setstr(ValueError, "mpz.pow to negative exponent");
606 return NULL;
607 }
608
609 if ((cmpres = mpz_cmp_ui(&a->mpz, (unsigned long int)0)) == 0) {
610 /* the base is 0 */
611
612 INCREF(mpz_value_zero);
613 return (object *)mpz_value_zero;
614 }
615 else if (cmpres > 0
616 && mpz_cmp_ui(&a->mpz, (unsigned long int)1) == 0) {
617 /* the base is 1 */
618
619 INCREF(mpz_value_one);
620 return (object *)mpz_value_one;
621 }
622 else if (cmpres < 0
623 && mpz_cmp_si(&a->mpz, (long int)-1) == 0) {
624
625 MP_INT tmpmpz;
626 /* the base is -1: pow(-1, any) == 1,-1 for even,uneven b */
627 /* XXX this code needs to be optimized: what's better?
628 mpz_mmod_ui or mpz_mod_2exp, I choose for the latter
629 for *un*obvious reasons */
630
631 /* is the exponent even? */
632 mpz_init(&tmpmpz);
633
634 /* look to the remainder after a division by (1 << 1) */
635 mpz_mod_2exp(&tmpmpz, &b->mpz, (unsigned long int)1);
636
637 if (mpz_cmp_ui(&tmpmpz, (unsigned int)0) == 0) {
638 mpz_clear(&tmpmpz);
639 INCREF(mpz_value_one);
640 return (object *)mpz_value_one;
641 }
642 mpz_clear(&tmpmpz);
643 INCREF(mpz_value_mone);
644 return (object *)mpz_value_mone;
645 }
646
647 #ifdef MPZ_LIB_DOES_CHECKING
648 /* check if it's doable: sizeof(exp) > sizeof(long) &&
649 abs(base) > 1 ?? --> No Way */
650 if (mpz_size(&b->mpz) > 1)
651 return (object *)err_nomem();
652 #else /* def MPZ_LIB_DOES_CHECKING */
653 /* wet finger method */
654 if (mpz_cmp_ui(&b->mpz, (unsigned long int)0x10000) >= 0) {
655 err_setstr(ValueError, "mpz.pow outrageous exponent");
656 return NULL;
657 }
658 #endif /* def MPZ_LIB_DOES_CHECKING else */
659
660 if ((z = newmpzobject()) == NULL)
661 return NULL;
662
663 mpz_pow_ui(&z->mpz, &a->mpz, mpz_get_ui(&b->mpz));
664
665 return (object *)z;
666 } /* mpz_power() */
667
668
669 static object *
670 mpz_negative(v)
671 mpzobject *v;
672 {
673 mpzobject *z;
674
675
676 #ifdef MPZ_SPARE_MALLOC
677 if (mpz_cmp_ui(&v->mpz, (unsigned long int)0) == 0) {
678 /* -0 == 0 */
679 INCREF(v);
680 return (object *)v;
681 }
682 #endif /* def MPZ_SPARE_MALLOC */
683
684 if ((z = newmpzobject()) == NULL)
685 return NULL;
686
687 mpz_neg(&z->mpz, &v->mpz);
688 return (object *)z;
689 } /* mpz_negative() */
690
691
692 static object *
693 mpz_positive(v)
694 mpzobject *v;
695 {
696 INCREF(v);
697 return (object *)v;
698 } /* mpz_positive() */
699
700
701 static object *
702 mpz_absolute(v)
703 mpzobject *v;
704 {
705 mpzobject *z;
706
707
708 if (mpz_cmp_ui(&v->mpz, (unsigned long int)0) >= 0) {
709 INCREF(v);
710 return (object *)v;
711 }
712
713 if ((z = newmpzobject()) == NULL)
714 return NULL;
715
716 mpz_neg(&z->mpz, &v->mpz);
717 return (object *)z;
718 } /* mpz_absolute() */
719
720 static int
721 mpz_nonzero(v)
722 mpzobject *v;
723 {
724 return mpz_cmp_ui(&v->mpz, (unsigned long int)0) != 0;
725 } /* mpz_nonzero() */
726
727 static object *
728 mpz_invert(v)
729 mpzobject *v;
730 {
731 mpzobject *z;
732
733
734 /* I think mpz_com does exactly what needed */
735 if ((z = newmpzobject()) == NULL)
736 return NULL;
737
738 mpz_com(&z->mpz, &v->mpz);
739 return (object *)z;
740 } /* mpz_invert() */
741
742 static object *
743 mpz_lshift(a, b)
744 mpzobject *a;
745 mpzobject *b;
746 {
747 int cmpres;
748 mpzobject *z;
749
750
751 if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
752 /* a << 0 == a */
753 INCREF(a);
754 return (object *)a;
755 }
756
757 if (cmpres < 0) {
758 err_setstr(ValueError, "mpz.<< negative shift count");
759 return NULL;
760 }
761
762 #ifdef MPZ_LIB_DOES_CHECKING
763 if (mpz_size(&b->mpz) > 1)
764 return (object *)err_nomem();
765 #else /* def MPZ_LIB_DOES_CHECKING */
766 /* wet finger method */
767 if (mpz_cmp_ui(&b->mpz, (unsigned long int)0x10000) >= 0) {
768 err_setstr(ValueError, "mpz.<< outrageous shift count");
769 return NULL;
770 }
771 #endif /* def MPZ_LIB_DOES_CHECKING else */
772
773 if ((z = newmpzobject()) == NULL)
774 return NULL;
775
776 mpz_mul_2exp(&z->mpz, &a->mpz, mpz_get_ui(&b->mpz));
777 return (object *)z;
778 } /* mpz_lshift() */
779
780 static object *
781 mpz_rshift(a, b)
782 mpzobject *a;
783 mpzobject *b;
784 {
785 int cmpres;
786 mpzobject *z;
787
788
789 if ((cmpres = mpz_cmp_ui(&b->mpz, (unsigned long int)0)) == 0) {
790 /* a >> 0 == a */
791 INCREF(a);
792 return (object *)a;
793 }
794
795 if (cmpres < 0) {
796 err_setstr(ValueError, "mpz.>> negative shift count");
797 return NULL;
798 }
799
800 if (mpz_size(&b->mpz) > 1)
801 return (object *)err_nomem();
802
803 if ((z = newmpzobject()) == NULL)
804 return NULL;
805
806 mpz_div_2exp(&z->mpz, &a->mpz, mpz_get_ui(&b->mpz));
807 return (object *)z;
808 } /* mpz_rshift() */
809
810 static object *
811 mpz_andfunc(a, b)
812 mpzobject *a;
813 mpzobject *b;
814 {
815 mpzobject *z;
816
817
818 if ((z = newmpzobject()) == NULL)
819 return NULL;
820
821 mpz_and(&z->mpz, &a->mpz, &b->mpz);
822 return (object *)z;
823 } /* mpz_andfunc() */
824
825 /* hack Hack HAck HACk HACK, XXX this code is dead slow */
826 void
827 mpz_xor(res, op1, op2)
828 MP_INT *res;
829 const MP_INT *op1;
830 const MP_INT *op2;
831 {
832 MP_INT tmpmpz;
833
834 mpz_init(&tmpmpz);
835
836 mpz_and(res, op1, op2);
837 mpz_com(&tmpmpz, res);
838 mpz_ior(res, op1, op2);
839 mpz_and(res, res, &tmpmpz);
840
841 mpz_clear(&tmpmpz);
842 } /* mpz_xor() HACK */
843
844 static object *
845 mpz_xorfunc(a, b)
846 mpzobject *a;
847 mpzobject *b;
848 {
849 mpzobject *z;
850
851
852 if ((z = newmpzobject()) == NULL)
853 return NULL;
854
855 mpz_xor(&z->mpz, &a->mpz, &b->mpz);
856 return (object *)z;
857 } /* mpz_xorfunc() */
858
859 static object *
860 mpz_orfunc(a, b)
861 mpzobject *a;
862 mpzobject *b;
863 {
864 mpzobject *z;
865
866
867 if ((z = newmpzobject()) == NULL)
868 return NULL;
869
870 mpz_ior(&z->mpz, &a->mpz, &b->mpz);
871 return (object *)z;
872 } /* mpz_orfunc() */
873
874 /* MPZ initialisation */
875
876 #include "longintrepr.h"
877
878 static object *
879 MPZ_mpz(self, args)
880 object *self;
881 object *args;
882 {
883 mpzobject *mpzp;
884 object *objp;
885
886
887 #ifdef MPZ_DEBUG
888 fputs("MPZ_mpz() called...\n", stderr);
889 #endif /* def MPZ_DEBUG */
890
891 if (!getargs(args, "O", &objp))
892 return NULL;
893
894 /* at least we know it's some object */
895 /* note DON't DECREF args NEITHER objp */
896
897 if (is_intobject(objp)) {
898 long lval;
899
900 if (!getargs(objp, "l", &lval))
901 return NULL;
902
903 if (lval == (long)0) {
904 INCREF(mpz_value_zero);
905 mpzp = mpz_value_zero;
906 }
907 else if (lval == (long)1) {
908 INCREF(mpz_value_one);
909 mpzp = mpz_value_one;
910 }
911 else if ((mpzp = newmpzobject()) == NULL)
912 return NULL;
913 else mpz_set_si(&mpzp->mpz, lval);
914 }
915 else if (is_longobject(objp)) {
916 MP_INT mplongdigit;
917 int i;
918 unsigned char isnegative;
919
920
921 if ((mpzp = newmpzobject()) == NULL)
922 return NULL;
923
924 mpz_set_si(&mpzp->mpz, 0L);
925 mpz_init(&mplongdigit);
926
927 /* how we're gonna handle this? */
928 if (isnegative = ((i = ((longobject *)objp)->ob_size) < 0) )
929 i = -i;
930
931 while (i--) {
932 mpz_set_ui(&mplongdigit,
933 (unsigned long)
934 ((longobject *)objp)->ob_digit[i]);
935 mpz_mul_2exp(&mplongdigit,&mplongdigit,
936 (unsigned long int)i * SHIFT);
937 mpz_ior(&mpzp->mpz, &mpzp->mpz, &mplongdigit);
938 }
939
940 if (isnegative)
941 mpz_neg(&mpzp->mpz, &mpzp->mpz);
942
943 /* get rid of allocation for tmp variable */
944 mpz_clear(&mplongdigit);
945 }
946 else if (is_stringobject(objp)) {
947 char *cp;
948 int len;
949 MP_INT mplongdigit;
950
951 if (!getargs(objp, "s#", &cp, &len))
952 return NULL;
953
954 if ((mpzp = newmpzobject()) == NULL)
955 return NULL;
956
957 mpz_set_si(&mpzp->mpz, 0L);
958 mpz_init(&mplongdigit);
959
960 /* let's do it the same way as with the long conversion:
961 without thinking how it can be faster (-: :-) */
962
963 cp += len;
964 while (len--) {
965 mpz_set_ui(&mplongdigit, (unsigned long)*--cp );
966 mpz_mul_2exp(&mplongdigit,&mplongdigit,
967 (unsigned long int)len * 8);
968 mpz_ior(&mpzp->mpz, &mpzp->mpz, &mplongdigit);
969 }
970
971 /* get rid of allocation for tmp variable */
972 mpz_clear(&mplongdigit);
973 }
974 else if (is_mpzobject(objp)) {
975 INCREF(objp);
976 mpzp = (mpzobject *)objp;
977 }
978 else {
979 err_setstr(TypeError,
980 "mpz.mpz() expects integer, long, string or mpz object argument");
981 return NULL;
982 }
983
984
985 #ifdef MPZ_DEBUG
986 fputs("MPZ_mpz: created mpz=", stderr);
987 mpz_out_str(stderr, 10, &mpzp->mpz);
988 putc('\n', stderr);
989 #endif /* def MPZ_DEBUG */
990 return (object *)mpzp;
991 } /* MPZ_mpz() */
992
993 static mpzobject *
994 mpz_mpzcoerce(z)
995 object *z;
996 {
997 /* shortcut: 9 out of 10 times the type is already ok */
998 if (is_mpzobject(z)) {
999 INCREF(z);
1000 return (mpzobject *)z; /* coercion succeeded */
1001 }
1002
1003 /* what types do we accept?: intobjects and longobjects */
1004 if (is_intobject(z) || is_longobject(z))
1005 return (mpzobject *)MPZ_mpz((object *)NULL, z);
1006
1007 err_setstr(TypeError, "number coercion (to mpzobject) failed");
1008 return NULL;
1009 } /* mpz_mpzcoerce() */
1010
1011 /* Forward */
1012 static void mpz_divm PROTO((MP_INT *res, const MP_INT *num,
1013 const MP_INT *den, const MP_INT *mod));
1014
1015 static object *
1016 MPZ_powm(self, args)
1017 object *self;
1018 object *args;
1019 {
1020 object *base, *exp, *mod;
1021 mpzobject *mpzbase = NULL, *mpzexp = NULL, *mpzmod = NULL;
1022 mpzobject *z;
1023 int tstres;
1024
1025
1026 if (!getargs(args, "(OOO)", &base, &exp, &mod))
1027 return NULL;
1028
1029 if ((mpzbase = mpz_mpzcoerce(base)) == NULL
1030 || (mpzexp = mpz_mpzcoerce(exp)) == NULL
1031 || (mpzmod = mpz_mpzcoerce(mod)) == NULL
1032 || (z = newmpzobject()) == NULL) {
1033 XDECREF(mpzbase);
1034 XDECREF(mpzexp);
1035 XDECREF(mpzmod);
1036 return NULL;
1037 }
1038
1039 if ((tstres=mpz_cmp_ui(&mpzexp->mpz, (unsigned long int)0)) == 0) {
1040 INCREF(mpz_value_one);
1041 return (object *)mpz_value_one;
1042 }
1043
1044 if (tstres < 0) {
1045 MP_INT absexp;
1046 /* negative exp */
1047
1048 mpz_init_set(&absexp, &mpzexp->mpz);
1049 mpz_abs(&absexp, &absexp);
1050 mpz_powm(&z->mpz, &mpzbase->mpz, &absexp, &mpzmod->mpz);
1051
1052 mpz_divm(&z->mpz, &mpz_value_one->mpz, &z->mpz, &mpzmod->mpz);
1053
1054 mpz_clear(&absexp);
1055 }
1056 else {
1057 mpz_powm(&z->mpz, &mpzbase->mpz, &mpzexp->mpz, &mpzmod->mpz);
1058 }
1059
1060 DECREF(mpzbase);
1061 DECREF(mpzexp);
1062 DECREF(mpzmod);
1063
1064 return (object *)z;
1065 } /* MPZ_powm() */
1066
1067
1068 static object *
1069 MPZ_gcd(self, args)
1070 object *self;
1071 object *args;
1072 {
1073 object *op1, *op2;
1074 mpzobject *mpzop1 = NULL, *mpzop2 = NULL;
1075 mpzobject *z;
1076
1077
1078 if (!getargs(args, "(OO)", &op1, &op2))
1079 return NULL;
1080
1081 if ((mpzop1 = mpz_mpzcoerce(op1)) == NULL
1082 || (mpzop2 = mpz_mpzcoerce(op2)) == NULL
1083 || (z = newmpzobject()) == NULL) {
1084 XDECREF(mpzop1);
1085 XDECREF(mpzop2);
1086 return NULL;
1087 }
1088
1089 /* ok, we have three mpzobjects, and an initialised result holder */
1090 mpz_gcd(&z->mpz, &mpzop1->mpz, &mpzop2->mpz);
1091
1092 DECREF(mpzop1);
1093 DECREF(mpzop2);
1094
1095 return (object *)z;
1096 } /* MPZ_gcd() */
1097
1098
1099 static object *
1100 MPZ_gcdext(self, args)
1101 object *self;
1102 object *args;
1103 {
1104 object *op1, *op2, *z = NULL;
1105 mpzobject *mpzop1 = NULL, *mpzop2 = NULL;
1106 mpzobject *g = NULL, *s = NULL, *t = NULL;
1107
1108
1109 if (!getargs(args, "(OO)", &op1, &op2))
1110 return NULL;
1111
1112 if ((mpzop1 = mpz_mpzcoerce(op1)) == NULL
1113 || (mpzop2 = mpz_mpzcoerce(op2)) == NULL
1114 || (z = newtupleobject(3)) == NULL
1115 || (g = newmpzobject()) == NULL
1116 || (s = newmpzobject()) == NULL
1117 || (t = newmpzobject()) == NULL) {
1118 XDECREF(mpzop1);
1119 XDECREF(mpzop2);
1120 XDECREF(z);
1121 XDECREF(g);
1122 XDECREF(s);
1123 /*XDECREF(t);*/
1124 return NULL;
1125 }
1126
1127 mpz_gcdext(&g->mpz, &s->mpz, &t->mpz, &mpzop1->mpz, &mpzop2->mpz);
1128
1129 DECREF(mpzop1);
1130 DECREF(mpzop2);
1131
1132 (void)settupleitem(z, 0, (object *)g);
1133 (void)settupleitem(z, 1, (object *)s);
1134 (void)settupleitem(z, 2, (object *)t);
1135
1136 return (object *)z;
1137 } /* MPZ_gcdext() */
1138
1139
1140 static object *
1141 MPZ_sqrt(self, args)
1142 object *self;
1143 object *args;
1144 {
1145 object *op;
1146 mpzobject *mpzop = NULL;
1147 mpzobject *z;
1148
1149
1150 if (!getargs(args, "O", &op))
1151 return NULL;
1152
1153 if ((mpzop = mpz_mpzcoerce(op)) == NULL
1154 || (z = newmpzobject()) == NULL) {
1155 XDECREF(mpzop);
1156 return NULL;
1157 }
1158
1159 mpz_sqrt(&z->mpz, &mpzop->mpz);
1160
1161 DECREF(mpzop);
1162
1163 return (object *)z;
1164 } /* MPZ_sqrt() */
1165
1166
1167 static object *
1168 MPZ_sqrtrem(self, args)
1169 object *self;
1170 object *args;
1171 {
1172 object *op, *z = NULL;
1173 mpzobject *mpzop = NULL;
1174 mpzobject *root = NULL, *rem = NULL;
1175
1176
1177 if (!getargs(args, "O", &op))
1178 return NULL;
1179
1180 if ((mpzop = mpz_mpzcoerce(op)) == NULL
1181 || (z = newtupleobject(2)) == NULL
1182 || (root = newmpzobject()) == NULL
1183 || (rem = newmpzobject()) == NULL) {
1184 XDECREF(mpzop);
1185 XDECREF(z);
1186 XDECREF(root);
1187 /*XDECREF(rem);*/
1188 return NULL;
1189 }
1190
1191 mpz_sqrtrem(&root->mpz, &rem->mpz, &mpzop->mpz);
1192
1193 DECREF(mpzop);
1194
1195 (void)settupleitem(z, 0, (object *)root);
1196 (void)settupleitem(z, 1, (object *)rem);
1197
1198 return (object *)z;
1199 } /* MPZ_sqrtrem() */
1200
1201
1202 static void
1203 #if __STDC__
1204 mpz_divm(MP_INT *res, const MP_INT *num, const MP_INT *den, const MP_INT *mod)
1205 #else
1206 mpz_divm(res, num, den, mod)
1207 MP_INT *res;
1208 const MP_INT *num;
1209 const MP_INT *den;
1210 const MP_INT *mod;
1211 #endif
1212 {
1213 MP_INT s0, s1, q, r, x, d0, d1;
1214
1215 mpz_init_set(&s0, num);
1216 mpz_init_set_ui(&s1, 0);
1217 mpz_init(&q);
1218 mpz_init(&r);
1219 mpz_init(&x);
1220 mpz_init_set(&d0, den);
1221 mpz_init_set(&d1, mod);
1222
1223 while (d1.size != 0) {
1224 mpz_divmod(&q, &r, &d0, &d1);
1225 mpz_set(&d0, &d1);
1226 mpz_set(&d1, &r);
1227
1228 mpz_mul(&x, &s1, &q);
1229 mpz_sub(&x, &s0, &x);
1230 mpz_set(&s0, &s1);
1231 mpz_set(&s1, &x);
1232 }
1233
1234 if (d0.size != 1 || d0.d[0] != 1)
1235 res->size = 0; /* trouble: the gcd != 1; set s to zero */
1236 else {
1237 #ifdef MPZ_MDIV_BUG
1238 /* watch out here! first check the signs, and then perform
1239 the mpz_mod() since mod could point to res */
1240 if ((s0.size < 0) != (mod->size < 0)) {
1241 mpz_mod(res, &s0, mod);
1242
1243 if (res->size)
1244 mpz_add(res, res, mod);
1245 }
1246 else
1247 mpz_mod(res, &s0, mod);
1248
1249 #else /* def MPZ_MDIV_BUG */
1250 mpz_mmod(res, &s0, mod);
1251 #endif /* def MPZ_MDIV_BUG else */
1252 }
1253
1254 mpz_clear(&s0);
1255 mpz_clear(&s1);
1256 mpz_clear(&q);
1257 mpz_clear(&r);
1258 mpz_clear(&x);
1259 mpz_clear(&d0);
1260 mpz_clear(&d1);
1261 } /* mpz_divm() */
1262
1263
1264 static object *
1265 MPZ_divm(self, args)
1266 object *self;
1267 object *args;
1268 {
1269 object *num, *den, *mod;
1270 mpzobject *mpznum, *mpzden, *mpzmod = NULL;
1271 mpzobject *z = NULL;
1272
1273
1274 if (!getargs(args, "(OOO)", &num, &den, &mod))
1275 return NULL;
1276
1277 if ((mpznum = mpz_mpzcoerce(num)) == NULL
1278 || (mpzden = mpz_mpzcoerce(den)) == NULL
1279 || (mpzmod = mpz_mpzcoerce(mod)) == NULL
1280 || (z = newmpzobject()) == NULL ) {
1281 XDECREF(mpznum);
1282 XDECREF(mpzden);
1283 XDECREF(mpzmod);
1284 return NULL;
1285 }
1286
1287 mpz_divm(&z->mpz, &mpznum->mpz, &mpzden->mpz, &mpzmod->mpz);
1288
1289 DECREF(mpznum);
1290 DECREF(mpzden);
1291 DECREF(mpzmod);
1292
1293 if (mpz_cmp_ui(&z->mpz, (unsigned long int)0) == 0) {
1294 DECREF(z);
1295 err_setstr(ValueError, "gcd(den, mod) != 1 or num == 0");
1296 return NULL;
1297 }
1298
1299 return (object *)z;
1300 } /* MPZ_divm() */
1301
1302
1303 /* MPZ methods-as-attributes */
1304 #ifdef MPZ_CONVERSIONS_AS_METHODS
1305 static object *
1306 mpz_int(self, args)
1307 mpzobject *self;
1308 object *args;
1309 #else /* def MPZ_CONVERSIONS_AS_METHODS */
1310 static object *
1311 mpz_int(self)
1312 mpzobject *self;
1313 #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
1314 {
1315 long sli;
1316
1317
1318 #ifdef MPZ_CONVERSIONS_AS_METHODS
1319 if (!getnoarg(args))
1320 return NULL;
1321 #endif /* def MPZ_CONVERSIONS_AS_METHODS */
1322
1323 if (mpz_size(&self->mpz) > 1
1324 || (sli = (long)mpz_get_ui(&self->mpz)) < (long)0 ) {
1325 err_setstr(ValueError, "mpz.int() arg too long to convert");
1326 return NULL;
1327 }
1328
1329 if (mpz_cmp_ui(&self->mpz, (unsigned long)0) < 0)
1330 sli = -sli;
1331
1332 return newintobject(sli);
1333 } /* mpz_int() */
1334
1335 static object *
1336 #ifdef MPZ_CONVERSIONS_AS_METHODS
1337 mpz_long(self, args)
1338 mpzobject *self;
1339 object *args;
1340 #else /* def MPZ_CONVERSIONS_AS_METHODS */
1341 mpz_long(self)
1342 mpzobject *self;
1343 #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
1344 {
1345 int i, isnegative;
1346 unsigned long int uli;
1347 longobject *longobjp;
1348 int ldcount;
1349 int bitpointer, newbitpointer;
1350 MP_INT mpzscratch;
1351
1352
1353 #ifdef MPZ_CONVERSIONS_AS_METHODS
1354 if (!getnoarg(args))
1355 return NULL;
1356 #endif /* def MPZ_CONVERSIONS_AS_METHODS */
1357
1358 /* determine length of python-long to be allocated */
1359 if ((longobjp = alloclongobject(i = (int)
1360 ((mpz_size(&self->mpz) * BITS_PER_MP_LIMB
1361 + SHIFT - 1) /
1362 SHIFT))) == NULL)
1363 return NULL;
1364
1365 /* determine sign, and copy self to scratch var */
1366 mpz_init_set(&mpzscratch, &self->mpz);
1367 if (isnegative = (mpz_cmp_ui(&self->mpz, (unsigned long int)0) < 0))
1368 mpz_neg(&mpzscratch, &mpzscratch);
1369
1370 /* let those bits come, let those bits go,
1371 e.g. dismantle mpzscratch, build longobject */
1372
1373 bitpointer = 0; /* the number of valid bits in stock */
1374 newbitpointer = 0;
1375 ldcount = 0; /* the python-long limb counter */
1376 uli = (unsigned long int)0;
1377 while (i--) {
1378 longobjp->ob_digit[ldcount] = uli & MASK;
1379
1380 /* check if we've had enough bits for this digit */
1381 if (bitpointer < SHIFT) {
1382 uli = mpz_get_ui(&mpzscratch);
1383 longobjp->ob_digit[ldcount] |=
1384 (uli << bitpointer) & MASK;
1385 uli >>= SHIFT-bitpointer;
1386 bitpointer += BITS_PER_MP_LIMB;
1387 mpz_div_2exp(&mpzscratch, &mpzscratch,
1388 BITS_PER_MP_LIMB);
1389 }
1390 else
1391 uli >>= SHIFT;
1392 bitpointer -= SHIFT;
1393 ldcount++;
1394 }
1395
1396 assert(mpz_cmp_ui(&mpzscratch, (unsigned long int)0) == 0);
1397 mpz_clear(&mpzscratch);
1398 assert(ldcount <= longobjp->ob_size);
1399
1400 /* long_normalize() is file-static */
1401 /* longobjp = long_normalize(longobjp); */
1402 while (ldcount > 0 && longobjp->ob_digit[ldcount-1] == 0)
1403 ldcount--;
1404 longobjp->ob_size = ldcount;
1405
1406
1407 if (isnegative)
1408 longobjp->ob_size = -longobjp->ob_size;
1409
1410 return (object *)longobjp;
1411
1412 } /* mpz_long() */
1413
1414
1415 /* I would have avoided pow() anyways, so ... */
1416 static const double multiplier = 256.0 * 256.0 * 256.0 * 256.0;
1417
1418 #ifdef MPZ_CONVERSIONS_AS_METHODS
1419 static object *
1420 mpz_float(self, args)
1421 mpzobject *self;
1422 object *args;
1423 #else /* def MPZ_CONVERSIONS_AS_METHODS */
1424 static object *
1425 mpz_float(self)
1426 mpzobject *self;
1427 #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
1428 {
1429 int i, isnegative;
1430 double x;
1431 double mulstate;
1432 MP_INT mpzscratch;
1433
1434
1435 #ifdef MPZ_CONVERSIONS_AS_METHODS
1436 if (!getnoarg(args))
1437 return NULL;
1438 #endif /* def MPZ_CONVERSIONS_AS_METHODS */
1439
1440 i = (int)mpz_size(&self->mpz);
1441
1442 /* determine sign, and copy abs(self) to scratch var */
1443 if (isnegative = (mpz_cmp_ui(&self->mpz, (unsigned long int)0) < 0)) {
1444 mpz_init(&mpzscratch);
1445 mpz_neg(&mpzscratch, &self->mpz);
1446 }
1447 else
1448 mpz_init_set(&mpzscratch, &self->mpz);
1449
1450 /* let those bits come, let those bits go,
1451 e.g. dismantle mpzscratch, build floatobject */
1452
1453 x = 0.0;
1454 mulstate = 1.0;
1455 while (i--) {
1456 x += mulstate * mpz_get_ui(&mpzscratch);
1457 mulstate *= multiplier;
1458 mpz_div_2exp(&mpzscratch, &mpzscratch, BITS_PER_MP_LIMB);
1459 }
1460
1461 assert(mpz_cmp_ui(&mpzscratch, (unsigned long int)0) == 0);
1462 mpz_clear(&mpzscratch);
1463
1464 if (isnegative)
1465 x = -x;
1466
1467 return newfloatobject(x);
1468
1469 } /* mpz_float() */
1470
1471 #ifdef MPZ_CONVERSIONS_AS_METHODS
1472 static object *
1473 mpz_hex(self, args)
1474 mpzobject *self;
1475 object *args;
1476 #else /* def MPZ_CONVERSIONS_AS_METHODS */
1477 static object *
1478 mpz_hex(self)
1479 mpzobject *self;
1480 #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
1481 {
1482 #ifdef MPZ_CONVERSIONS_AS_METHODS
1483 if (!getnoarg(args))
1484 return NULL;
1485 #endif /* def MPZ_CONVERSIONS_AS_METHODS */
1486
1487 return mpz_format(self, 16, (unsigned char)1);
1488 } /* mpz_hex() */
1489
1490 #ifdef MPZ_CONVERSIONS_AS_METHODS
1491 static object *
1492 mpz_oct(self, args)
1493 mpzobject *self;
1494 object *args;
1495 #else /* def MPZ_CONVERSIONS_AS_METHODS */
1496 static object *
1497 mpz_oct(self)
1498 mpzobject *self;
1499 #endif /* def MPZ_CONVERSIONS_AS_METHODS else */
1500 {
1501 #ifdef MPZ_CONVERSIONS_AS_METHODS
1502 if (!getnoarg(args))
1503 return NULL;
1504 #endif /* def MPZ_CONVERSIONS_AS_METHODS */
1505
1506 return mpz_format(self, 8, (unsigned char)1);
1507 } /* mpz_oct() */
1508
1509 static object *
1510 mpz_binary(self, args)
1511 mpzobject *self;
1512 object *args;
1513 {
1514 int size;
1515 stringobject *strobjp;
1516 char *cp;
1517 MP_INT mp;
1518 unsigned long ldigit;
1519
1520 if (!getnoarg(args))
1521 return NULL;
1522
1523 if (mpz_cmp_ui(&self->mpz, (unsigned long int)0) < 0) {
1524 err_setstr(ValueError, "mpz.binary() arg must be >= 0");
1525 return NULL;
1526 }
1527
1528 mpz_init_set(&mp, &self->mpz);
1529 size = (int)mpz_size(&mp);
1530
1531 if ((strobjp = (stringobject *)
1532 newsizedstringobject((char *)0,
1533 size * sizeof (unsigned long int))) == NULL)
1534 return NULL;
1535
1536 /* get the beginning of the string memory and start copying things */
1537 cp = GETSTRINGVALUE(strobjp);
1538
1539 /* this has been programmed using a (fairly) decent lib-i/f it could
1540 be must faster if we looked into the GMP lib */
1541 while (size--) {
1542 ldigit = mpz_get_ui(&mp);
1543 mpz_div_2exp(&mp, &mp, BITS_PER_MP_LIMB);
1544 *cp++ = (unsigned char)(ldigit & 0xFF);
1545 *cp++ = (unsigned char)((ldigit >>= 8) & 0xFF);
1546 *cp++ = (unsigned char)((ldigit >>= 8) & 0xFF);
1547 *cp++ = (unsigned char)((ldigit >>= 8) & 0xFF);
1548 }
1549
1550 while (strobjp->ob_size && !*--cp)
1551 strobjp->ob_size--;
1552
1553 return (object *)strobjp;
1554 } /* mpz_binary() */
1555
1556
1557 static struct methodlist mpz_methods[] = {
1558 #ifdef MPZ_CONVERSIONS_AS_METHODS
1559 {"int", mpz_int},
1560 {"long", mpz_long},
1561 {"float", mpz_float},
1562 {"hex", mpz_hex},
1563 {"oct", mpz_oct},
1564 #endif /* def MPZ_CONVERSIONS_AS_METHODS */
1565 {"binary", (object *(*)(object *, object *))mpz_binary},
1566 {NULL, NULL} /* sentinel */
1567 };
1568
1569 static object *
1570 mpz_getattr(self, name)
1571 mpzobject *self;
1572 char *name;
1573 {
1574 return findmethod(mpz_methods, (object *)self, name);
1575 } /* mpz_getattr() */
1576
1577
1578 static int
1579 mpz_coerce(pv, pw)
1580 object **pv;
1581 object **pw;
1582 {
1583 object *z;
1584
1585 #ifdef MPZ_DEBUG
1586 fputs("mpz_coerce() called...\n", stderr);
1587 #endif /* def MPZ_DEBUG */
1588
1589 assert(is_mpzobject(*pv));
1590
1591 /* always convert other arg to mpz value, except for floats */
1592 if (!is_floatobject(*pw)) {
1593 if ((z = (object *)mpz_mpzcoerce(*pw)) == NULL)
1594 return -1; /* -1: an error always has been set */
1595
1596 INCREF(*pv);
1597 *pw = z;
1598 }
1599 else {
1600 if ((z = mpz_float(*pv, NULL)) == NULL)
1601 return -1;
1602
1603 INCREF(*pw);
1604 *pv = z;
1605 }
1606 return 0; /* coercion succeeded */
1607
1608 } /* mpz_coerce() */
1609
1610
1611 static object *
1612 mpz_repr(v)
1613 object *v;
1614 {
1615 return mpz_format(v, 10, (unsigned char)1);
1616 } /* mpz_repr() */
1617
1618
1619
1620 #define UF (unaryfunc)
1621 #define BF (binaryfunc)
1622 #define TF (ternaryfunc)
1623 #define IF (inquiry)
1624 #define CF (coercion)
1625
1626 static number_methods mpz_as_number = {
1627 BF mpz_addition, /*nb_add*/
1628 BF mpz_substract, /*nb_subtract*/
1629 BF mpz_multiply, /*nb_multiply*/
1630 BF mpz_divide, /*nb_divide*/
1631 BF mpz_remainder, /*nb_remainder*/
1632 BF mpz_div_and_mod, /*nb_divmod*/
1633 TF mpz_power, /*nb_power*/
1634 UF mpz_negative, /*nb_negative*/
1635 UF mpz_positive, /*tp_positive*/
1636 UF mpz_absolute, /*tp_absolute*/
1637 IF mpz_nonzero, /*tp_nonzero*/
1638 UF mpz_invert, /*nb_invert*/
1639 BF mpz_lshift, /*nb_lshift*/
1640 BF mpz_rshift, /*nb_rshift*/
1641 BF mpz_andfunc, /*nb_and*/
1642 BF mpz_xorfunc, /*nb_xor*/
1643 BF mpz_orfunc, /*nb_or*/
1644 CF mpz_coerce, /*nb_coerce*/
1645 #ifndef MPZ_CONVERSIONS_AS_METHODS
1646 UF mpz_int, /*nb_int*/
1647 UF mpz_long, /*nb_long*/
1648 UF mpz_float, /*nb_float*/
1649 UF mpz_oct, /*nb_oct*/
1650 UF mpz_hex, /*nb_hex*/
1651 #endif /* ndef MPZ_CONVERSIONS_AS_METHODS */
1652 };
1653
1654 static typeobject MPZtype = {
1655 OB_HEAD_INIT(&Typetype)
1656 0, /*ob_size*/
1657 "mpz", /*tp_name*/
1658 sizeof(mpzobject), /*tp_size*/
1659 0, /*tp_itemsize*/
1660 /* methods */
1661 (destructor)mpz_dealloc, /*tp_dealloc*/
1662 0, /*tp_print*/
1663 (getattrfunc)mpz_getattr, /*tp_getattr*/
1664 0, /*tp_setattr*/
1665 (cmpfunc)mpz_compare, /*tp_compare*/
1666 (reprfunc)mpz_repr, /*tp_repr*/
1667 &mpz_as_number, /*tp_as_number*/
1668 };
1669
1670 /* List of functions exported by this module */
1671
1672 static struct methodlist mpz_functions[] = {
1673 #if 0
1674 {initialiser_name, MPZ_mpz},
1675 #else /* 0 */
1676 /* until guido ``fixes'' struct methodlist */
1677 {(char *)initialiser_name, MPZ_mpz},
1678 #endif /* 0 else */
1679 {"powm", MPZ_powm},
1680 {"gcd", MPZ_gcd},
1681 {"gcdext", MPZ_gcdext},
1682 {"sqrt", MPZ_sqrt},
1683 {"sqrtrem", MPZ_sqrtrem},
1684 {"divm", MPZ_divm},
1685 {NULL, NULL} /* Sentinel */
1686 };
1687
1688
1689 /* #define MP_TEST_ALLOC */
1690
1691 #ifdef MP_TEST_ALLOC
1692 #define MP_TEST_SIZE 4
1693 static const char mp_test_magic[MP_TEST_SIZE] = {'\xAA','\xAA','\xAA','\xAA'};
1694 static mp_test_error( location )
1695 int *location;
1696 {
1697 /* assumptions: *alloc returns address dividable by 4,
1698 mpz_* routines allocate in chunks dividable by four */
1699 fprintf(stderr, "MP_TEST_ERROR: location holds 0x%08d\n", *location );
1700 fatal("MP_TEST_ERROR");
1701 } /* static mp_test_error() */
1702 #define MP_EXTRA_ALLOC(size) ((size) + MP_TEST_SIZE)
1703 #define MP_SET_TEST(basep,size) (void)memcpy( ((char *)(basep))+(size), mp_test_magic, MP_TEST_SIZE)
1704 #define MP_DO_TEST(basep,size) if ( !memcmp( ((char *)(basep))+(size), mp_test_magic, MP_TEST_SIZE ) ) \
1705 ; \
1706 else \
1707 mp_test_error((int *)((char *)(basep) + size))
1708 #else /* def MP_TEST_ALLOC */
1709 #define MP_EXTRA_ALLOC(size) (size)
1710 #define MP_SET_TEST(basep,size)
1711 #define MP_DO_TEST(basep,size)
1712 #endif /* def MP_TEST_ALLOC else */
1713
1714 void *mp_allocate( alloc_size )
1715 size_t alloc_size;
1716 {
1717 void *res;
1718
1719 #ifdef MPZ_DEBUG
1720 fprintf(stderr, "mp_allocate : size %ld\n",
1721 alloc_size);
1722 #endif /* def MPZ_DEBUG */
1723
1724 if ( (res = malloc(MP_EXTRA_ALLOC(alloc_size))) == NULL )
1725 fatal("mp_allocate failure");
1726
1727 #ifdef MPZ_DEBUG
1728 fprintf(stderr, "mp_allocate : address 0x%08x\n", res);
1729 #endif /* def MPZ_DEBUG */
1730
1731 MP_SET_TEST(res,alloc_size);
1732
1733 return res;
1734 } /* mp_allocate() */
1735
1736
1737 void *mp_reallocate( ptr, old_size, new_size )
1738 void *ptr;
1739 size_t old_size;
1740 size_t new_size;
1741 {
1742 void *res;
1743
1744 #ifdef MPZ_DEBUG
1745 fprintf(stderr, "mp_reallocate: old address 0x%08x, old size %ld\n",
1746 ptr, old_size);
1747 #endif /* def MPZ_DEBUG */
1748
1749 MP_DO_TEST(ptr, old_size);
1750
1751 if ( (res = realloc(ptr, MP_EXTRA_ALLOC(new_size))) == NULL )
1752 fatal("mp_reallocate failure");
1753
1754 #ifdef MPZ_DEBUG
1755 fprintf(stderr, "mp_reallocate: new address 0x%08x, new size %ld\n",
1756 res, new_size);
1757 #endif /* def MPZ_DEBUG */
1758
1759 MP_SET_TEST(res, new_size);
1760
1761 return res;
1762 } /* mp_reallocate() */
1763
1764
1765 void mp_free( ptr, size )
1766 void *ptr;
1767 size_t size;
1768 {
1769
1770 #ifdef MPZ_DEBUG
1771 fprintf(stderr, "mp_free : old address 0x%08x, old size %ld\n",
1772 ptr, size);
1773 #endif /* def MPZ_DEBUG */
1774
1775 MP_DO_TEST(ptr, size);
1776 free(ptr);
1777 } /* mp_free() */
1778
1779
1780
1781 /* Initialize this module. */
1782
1783 void
1784 initmpz()
1785 {
1786 #ifdef MPZ_DEBUG
1787 fputs( "initmpz() called...\n", stderr );
1788 #endif /* def MPZ_DEBUG */
1789
1790 mp_set_memory_functions( mp_allocate, mp_reallocate, mp_free );
1791 (void)initmodule("mpz", mpz_functions);
1792
1793 /* create some frequently used constants */
1794 if ((mpz_value_zero = newmpzobject()) == NULL)
1795 fatal("initmpz: can't initialize mpz constants");
1796 mpz_set_ui(&mpz_value_zero->mpz, (unsigned long int)0);
1797
1798 if ((mpz_value_one = newmpzobject()) == NULL)
1799 fatal("initmpz: can't initialize mpz constants");
1800 mpz_set_ui(&mpz_value_one->mpz, (unsigned long int)1);
1801
1802 if ((mpz_value_mone = newmpzobject()) == NULL)
1803 fatal("initmpz: can't initialize mpz constants");
1804 mpz_set_si(&mpz_value_mone->mpz, (long)-1);
1805
1806 } /* initmpz() */
1807 #ifdef MAKEDUMMYINT
1808 int _mpz_dummy_int; /* XXX otherwise, we're .bss-less (DYNLOAD->Jack?) */
1809 #endif /* def MAKEDUMMYINT */