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