chiark - git - mdw - catacomb-python/blob - mp.c

   1 /* -*-c-*-
   2  *
   3  * Multiprecision arithmetic
   4  *
   5  * (c) 2004 Straylight/Edgeware
   6  */
   7
   8 /*----- Licensing notice --------------------------------------------------*
   9  *
  10  * This file is part of the Python interface to Catacomb.
  11  *
  12  * Catacomb/Python is free software; you can redistribute it and/or modify
  13  * it under the terms of the GNU General Public License as published by
  14  * the Free Software Foundation; either version 2 of the License, or
  15  * (at your option) any later version.
  16  *
  17  * Catacomb/Python is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  20  * GNU General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU General Public License
  23  * along with Catacomb/Python; if not, write to the Free Software Foundation,
  24  * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
  25  */
  26
  27 /*----- Header files ------------------------------------------------------*/
  28
  29 #include "catacomb-python.h"
  30
  31 /*----- General utilities -------------------------------------------------*/
  32
  33 PyTypeObject *mp_pytype = 0;
  34 PyTypeObject *gf_pytype = 0;
  35
  36 mp *mp_frompylong(PyObject *obj)
  37 {
  38   unsigned long bits;
  39   PyLongObject *l = (PyLongObject *)obj;
  40   int sz;
  41   size_t w;
  42   mpd r = 0;
  43   int b = 0;
  44   int i;
  45   mp *x;
  46   mpw *p;
  47
  48   sz = l->ob_size;
  49   if (sz < 0) sz = -sz;
  50   assert(MPW_BITS >= SHIFT);
  51   bits = (unsigned long)sz * SHIFT;
  52   w = (bits + MPW_BITS - 1)/MPW_BITS;
  53   x = mp_new(w, l->ob_size < 0 ? MP_NEG : 0);
  54   p = x->v;
  55   for (i = 0; i < sz; i++) {
  56     r |= (mpd)l->ob_digit[i] << b;
  57     b += SHIFT;
  58     while (b >= MPW_BITS) {
  59       *p++ = MPW(r);
  60       r >>= MPW_BITS;
  61       b -= MPW_BITS;
  62     }
  63   }
  64   while (r) {
  65     *p++ = MPW(r);
  66     r >>= MPW_BITS;
  67   }
  68   x->vl = p;
  69   MP_SHRINK(x);
  70   return (x);
  71 }
  72
  73 PyObject *mp_topylong(mp *x)
  74 {
  75   unsigned long bits = mp_bits(x);
  76   int sz = (bits + SHIFT - 1)/SHIFT;
  77   PyLongObject *l = _PyLong_New(sz);
  78   mpd r = 0;
  79   int b = 0;
  80   mpw *p = x->v;
  81   int i = 0;
  82
  83   assert(MPW_BITS >= SHIFT);
  84   while (i < sz && p < x->vl) {
  85     r |= (mpd)*p++ << b;
  86     b += MPW_BITS;
  87     while (i < sz && b >= SHIFT) {
  88       l->ob_digit[i++] = r & MASK;
  89       r >>= SHIFT;
  90       b -= SHIFT;
  91     }
  92   }
  93   while (i < sz && r) {
  94     l->ob_digit[i++] = r & MASK;
  95     r >>= SHIFT;
  96   }
  97   l->ob_size = (x->f & MP_NEG) ? -sz : sz;
  98   return ((PyObject *)l);
  99 }
 100
 101 mp *mp_frompyobject(PyObject *o, int radix)
 102 {
 103   mp *x;
 104
 105   if (PyString_Check(o)) {
 106     mptext_stringctx sc;
 107     mp *x;
 108     sc.buf = PyString_AS_STRING(o);
 109     sc.lim = sc.buf + PyString_GET_SIZE(o);
 110     x = mp_read(MP_NEW, radix, &mptext_stringops, &sc);
 111     if (!x) return (0);
 112     if (sc.buf < sc.lim) { MP_DROP(x); return (0); }
 113     return (x);
 114   }
 115   if ((x = tomp(o)) != 0)
 116     return (x);
 117   return (0);
 118 }
 119
 120 PyObject *mp_topystring(mp *x, int radix, const char *xpre,
 121                         const char *pre, const char *post)
 122 {
 123   int len = mptext_len(x, radix) + 1;
 124   mptext_stringctx sc;
 125   PyObject *o;
 126   size_t xprelen = xpre ? strlen(xpre) : 0;
 127   size_t prelen = pre ? strlen(pre) : 0;
 128   size_t postlen = post ? strlen(post) : 0;
 129   char *p;
 130   MP_COPY(x);
 131   o = PyString_FromStringAndSize(0, len + 1 + xprelen + prelen + postlen);
 132   p = PyString_AS_STRING(o);
 133   sc.buf = p;
 134   if (xpre) { memcpy(sc.buf, xpre, xprelen); sc.buf += xprelen; }
 135   if (MP_NEGP(x)) { *sc.buf++ = '-'; x = mp_neg(x, x); }
 136   if (pre) { memcpy(sc.buf, pre, prelen); sc.buf += prelen; }
 137   sc.lim = sc.buf + len;
 138   mp_write(x, radix, &mptext_stringops, &sc);
 139   if (post) { memcpy(sc.buf, post, postlen); sc.buf += postlen; }
 140   MP_DROP(x);
 141   _PyString_Resize(&o, sc.buf - p);
 142   return (o);
 143 }
 144
 145 static int good_radix_p(int r, int readp)
 146 {
 147   return ((r >= -255 && r <= -2) ||
 148           (readp && r == 0) ||
 149           (r >= 2 && r <= 62));
 150 }
 151
 152 PyObject *mp_pywrap(mp *x)
 153 {
 154   mp_pyobj *z = PyObject_New(mp_pyobj, mp_pytype);
 155   z->x = x;
 156   return ((PyObject *)z);
 157 }
 158
 159 PyObject *gf_pywrap(mp *x)
 160 {
 161   mp_pyobj *z = PyObject_New(mp_pyobj, gf_pytype);
 162   z->x = x;
 163   return ((PyObject *)z);
 164 }
 165
 166 int mp_tolong_checked(mp *x, long *l, int must)
 167 {
 168   static mp *longmin = 0, *longmax = 0;
 169   int rc = -1;
 170
 171   if (!longmax) {
 172     longmin = mp_fromlong(MP_NEW, LONG_MIN);
 173     longmax = mp_fromlong(MP_NEW, LONG_MAX);
 174   }
 175   if (MP_CMP(x, <, longmin) || MP_CMP(x, >, longmax)) {
 176     if (must) VALERR("mp out of range for int");
 177     else goto end;
 178   }
 179   *l = mp_tolong(x);
 180   rc = 0;
 181 end:
 182   return (rc);
 183 }
 184
 185 /*----- Python interface --------------------------------------------------*/
 186
 187 static void mp_pydealloc(PyObject *o)
 188 {
 189   MP_DROP(MP_X(o));
 190   FREEOBJ(o);
 191 }
 192
 193 static PyObject *mp_pyrepr(PyObject *o)
 194   { return mp_topystring(MP_X(o), 10, "MP(", 0, "L)"); }
 195
 196 static PyObject *mp_pystr(PyObject *o)
 197   { return mp_topystring(MP_X(o), 10, 0, 0, 0); }
 198
 199 mp *tomp(PyObject *o)
 200 {
 201   PyObject *l;
 202   mp *x;
 203
 204   if (!o)
 205     return (0);
 206   else if (MP_PYCHECK(o) || GF_PYCHECK(o))
 207     return (MP_COPY(MP_X(o)));
 208   else if (FE_PYCHECK(o))
 209     return (F_OUT(FE_F(o), MP_NEW, FE_X(o)));
 210   else if (PFILT_PYCHECK(o))
 211     return (MP_COPY(PFILT_F(o)->m));
 212   else if (ECPT_PYCHECK(o)) {
 213     ec p = EC_INIT;
 214     getecptout(&p, o);
 215     x = MP_COPY(p.x);
 216     EC_DESTROY(&p);
 217     return (x);
 218   } else if (GE_PYCHECK(o)) {
 219     if ((x = G_TOINT(GE_G(o), MP_NEW, GE_X(o))) == 0)
 220       return (0);
 221     return (x);
 222   } else if (PyInt_Check(o))
 223     return (mp_fromlong(MP_NEW, PyInt_AS_LONG(o)));
 224   else if ((l = PyNumber_Long(o)) != 0) {
 225     x = mp_frompylong(l);
 226     Py_DECREF(l);
 227     return (x);
 228   } else {
 229     PyErr_Clear();
 230     return (0);
 231   }
 232 }
 233
 234 mp *getmp(PyObject *o)
 235 {
 236   mp *x = 0;
 237   if (!o) return (0);
 238   if ((x = tomp(o)) == 0) {
 239     PyErr_Format(PyExc_TypeError, "can't convert %.100s to mp",
 240                  o->ob_type->tp_name);
 241   }
 242   return (x);
 243 }
 244
 245 int convmp(PyObject *o, void *p)
 246 {
 247   mp *x;
 248   if ((x = getmp(o)) == 0) return (0);
 249   *(mp **)p = x;
 250   return (1);
 251 }
 252
 253 mp *getgf(PyObject *o)
 254 {
 255   mp *x = 0;
 256   if (!o) return (0);
 257   if ((x = tomp(o)) == 0) {
 258     PyErr_Format(PyExc_TypeError, "can't convert %.100s to gf",
 259                  o->ob_type->tp_name);
 260   }
 261   return (x);
 262 }
 263
 264 int convgf(PyObject *o, void *p)
 265 {
 266   mp *x;
 267   if ((x = getgf(o)) == 0) return (0);
 268   *(mp **)p = x;
 269   return (1);
 270 }
 271
 272 static mp *implicitmp(PyObject *o)
 273 {
 274   if (!o ||
 275       GF_PYCHECK(o) ||
 276       ECPT_PYCHECK(o) ||
 277       FE_PYCHECK(o) ||
 278       GE_PYCHECK(o))
 279     return (0);
 280   return (tomp(o));
 281 }
 282
 283 static mp *implicitgf(PyObject *o)
 284 {
 285   if (!o ||
 286       MP_PYCHECK(o) ||
 287       ECPT_PYCHECK(o) ||
 288       FE_PYCHECK(o) ||
 289       GE_PYCHECK(o))
 290     return (0);
 291   return (tomp(o));
 292 }
 293
 294 static int mpbinop(PyObject *x, PyObject *y, mp **xx, mp **yy)
 295 {
 296   if ((*xx = implicitmp(x)) == 0)
 297     return (-1);
 298   if ((*yy = implicitmp(y)) == 0) {
 299     MP_DROP(*xx);
 300     return (-1);
 301   }
 302   return (0);
 303 }
 304
 305 static int gfbinop(PyObject *x, PyObject *y, mp **xx, mp **yy)
 306 {
 307   if ((*xx = implicitgf(x)) == 0)
 308     return (-1);
 309   if ((*yy = implicitgf(y)) == 0) {
 310     MP_DROP(*xx);
 311     return (-1);
 312   }
 313   return (0);
 314 }
 315
 316 #define gf_and mp_and
 317 #define gf_or mp_or
 318 #define gf_xor mp_xor
 319 #define BINOP(pre, name)                                                \
 320   static PyObject *pre##_py##name(PyObject *x, PyObject *y) {           \
 321     mp *xx, *yy, *zz;                                                   \
 322     if (pre##binop(x, y, &xx, &yy)) RETURN_NOTIMPL;                     \
 323     zz = pre##_##name(MP_NEW, xx, yy);                                  \
 324     MP_DROP(xx); MP_DROP(yy);                                           \
 325     return (pre##_pywrap(zz));                                          \
 326   }
 327 BINOP(mp, add)
 328 BINOP(mp, sub)
 329 BINOP(mp, mul)
 330 BINOP(mp, and2c)
 331 BINOP(mp, or2c)
 332 BINOP(mp, xor2c)
 333 BINOP(gf, add)
 334 BINOP(gf, sub)
 335 BINOP(gf, mul)
 336 BINOP(gf, and)
 337 BINOP(gf, or)
 338 BINOP(gf, xor)
 339 #undef BINOP
 340
 341 static mp *mp_abs(mp *d, mp *x)
 342 {
 343   if (MP_NEGP(x))
 344     return (mp_neg(d, x));
 345   MP_COPY(x);
 346   if (d) MP_DROP(d);
 347   return (x);
 348 }
 349
 350 #define UNOP(pre, name)                                                 \
 351   static PyObject *pre##_py##name(PyObject *x)                          \
 352     { return mp_pywrap(pre##_##name(MP_NEW, MP_X(x))); }
 353 UNOP(mp, neg)
 354 UNOP(mp, abs)
 355 UNOP(mp, not2c)
 356 #undef UNOP
 357
 358 static PyObject *mp_pyid(PyObject *x) { RETURN_OBJ(x); }
 359
 360 #define gf_lsr mp_lsr
 361 #define SHIFTOP(pre, name, rname)                                       \
 362   static PyObject *pre##_py##name(PyObject *x, PyObject *y) {           \
 363     mp *xx, *yy;                                                        \
 364     PyObject *z = 0;                                                    \
 365     long n;                                                             \
 366     if (pre##binop(x, y, &xx, &yy)) RETURN_NOTIMPL;                     \
 367     if (mp_tolong_checked(yy, &n, 1)) goto end;                         \
 368     if (n < 0)                                                          \
 369       z = pre##_pywrap(mp_##rname(MP_NEW, xx, -n));                     \
 370     else                                                                \
 371       z = pre##_pywrap(mp_##name(MP_NEW, xx, n));                       \
 372   end:                                                                  \
 373     MP_DROP(xx); MP_DROP(yy);                                           \
 374     return (z);                                                         \
 375   }
 376 SHIFTOP(mp, lsl2c, lsr2c)
 377 SHIFTOP(mp, lsr2c, lsl2c)
 378 SHIFTOP(gf, lsl, lsr)
 379 SHIFTOP(gf, lsr, lsl)
 380 #undef SHIFTOP
 381
 382 #define DIVOP(pre, name, qq, rr, gather)                                \
 383   static PyObject *pre##_py##name(PyObject *x, PyObject *y) {           \
 384     mp *xx, *yy;                                                        \
 385     PyObject *z = 0;                                                    \
 386     INIT_##qq(q) INIT_##rr(r)                                           \
 387     if (pre##binop(x, y, &xx, &yy)) RETURN_NOTIMPL;                     \
 388     if (MP_ZEROP(yy))                                                   \
 389       ZDIVERR("division by zero");                                      \
 390     pre##_div(ARG_##qq(q), ARG_##rr(r), xx, yy);                        \
 391     z = gather;                                                         \
 392   end:                                                                  \
 393     MP_DROP(xx); MP_DROP(yy);                                           \
 394     return (z);                                                         \
 395   }
 396 #define INIT_YES(p) mp *p = MP_NEW;
 397 #define INIT_NO(p)
 398 #define ARG_YES(p) &p
 399 #define ARG_NO(p) 0
 400 DIVOP(mp, divmod, YES, YES,
 401       Py_BuildValue("(NN)", mp_pywrap(q), mp_pywrap(r)))
 402 DIVOP(mp, div, YES, NO, mp_pywrap(q))
 403 DIVOP(mp, mod, NO, YES, mp_pywrap(r))
 404 DIVOP(gf, divmod, YES, YES,
 405       Py_BuildValue("(NN)", gf_pywrap(q), gf_pywrap(r)))
 406 DIVOP(gf, div, YES, NO, gf_pywrap(q))
 407 DIVOP(gf, mod, NO, YES, gf_pywrap(r))
 408 #undef INIT_YES
 409 #undef INIT_NO
 410 #undef ARG_YES
 411 #undef ARG_NO
 412 #undef DIVOP
 413
 414 static mp *mp_modinv_checked(mp *d, mp *x, mp *p)
 415 {
 416   mp *g = MP_NEW;
 417   mp_gcd(&g, 0, &d, p, x);
 418   if (!MP_EQ(g, MP_ONE)) {
 419     MP_DROP(g); MP_DROP(d);
 420     PyErr_SetString(PyExc_ZeroDivisionError, "no modular inverse");
 421     return (0);
 422   }
 423   MP_DROP(g); return (d);
 424 }
 425
 426 static PyObject *mp_pyexp(PyObject *x, PyObject *y, PyObject *z)
 427 {
 428   mp *xx = 0, *yy = 0, *zz = 0;
 429   mp *r = 0;
 430   PyObject *rc = 0;
 431
 432   if ((xx = implicitmp(x)) == 0 || (yy = implicitmp(y)) == 0 ||
 433       (z && z != Py_None && (zz = tomp(z)) == 0)) {
 434     mp_drop(xx); mp_drop(yy); mp_drop(zz);
 435     RETURN_NOTIMPL;
 436   }
 437   if (!z || z == Py_None) {
 438     if (MP_NEGP(yy)) VALERR("negative exponent");
 439     r = mp_exp(MP_NEW, xx, yy);
 440   } else {
 441     if (!MP_POSP(zz)) VALERR("modulus must be positive");
 442     if (MP_NEGP(yy)) {
 443       if ((xx = mp_modinv_checked(xx, xx, zz)) == 0) goto end;
 444       yy = mp_neg(yy, yy);
 445     }
 446     if (MP_ODDP(zz)) {
 447       mpmont mm;
 448       mpmont_create(&mm, zz);
 449       r = mpmont_exp(&mm, MP_NEW, xx, yy);
 450       mpmont_destroy(&mm);
 451     } else {
 452       mpbarrett mb;
 453       mpbarrett_create(&mb, zz);
 454       r = mpbarrett_exp(&mb, MP_NEW, xx, yy);
 455       mpbarrett_destroy(&mb);
 456     }
 457   }
 458   rc = mp_pywrap(r);
 459 end:
 460   mp_drop(xx); mp_drop(yy); mp_drop(zz);
 461   return (rc);
 462 }
 463
 464 static mp *gf_modinv_checked(mp *d, mp *x, mp *p)
 465 {
 466   mp *g = MP_NEW;
 467   gf_gcd(&g, 0, &d, p, x);
 468   if (!MP_EQ(g, MP_ONE)) {
 469     MP_DROP(g); MP_DROP(d);
 470     PyErr_SetString(PyExc_ZeroDivisionError, "no modular inverse");
 471     return (0);
 472   }
 473   MP_DROP(g); return (d);
 474 }
 475
 476 #define BASEOP(name, radix, pre)                                        \
 477   static PyObject *mp_py##name(PyObject *x)                             \
 478     { return mp_topystring(MP_X(x), radix, 0, pre, 0); }
 479 BASEOP(oct, 8, "0");
 480 BASEOP(hex, 16, "0x");
 481 #undef BASEOP
 482
 483 static int mp_pynonzerop(PyObject *x) { return !MP_ZEROP(MP_X(x)); }
 484
 485 static PyObject *mp_pyint(PyObject *x)
 486 {
 487   long l;
 488   if (!mp_tolong_checked(MP_X(x), &l, 0)) return (PyInt_FromLong(l));
 489   else return mp_topylong(MP_X(x));
 490 }
 491 static PyObject *mp_pylong(PyObject *x)
 492   { return (mp_topylong(MP_X(x))); }
 493 static PyObject *mp_pyfloat(PyObject *x)
 494 {
 495   PyObject *l = mp_topylong(MP_X(x));
 496   double f = PyLong_AsDouble(l);
 497   Py_DECREF(l);
 498   return (PyFloat_FromDouble(f));
 499 }
 500
 501 #define COERCE(pre, PRE)                                                \
 502   static int pre##_pycoerce(PyObject **x, PyObject **y)                 \
 503   {                                                                     \
 504     mp *z;                                                              \
 505                                                                         \
 506     if (PRE##_PYCHECK(*y)) {                                            \
 507       Py_INCREF(*x); Py_INCREF(*y);                                     \
 508       return (0);                                                       \
 509     }                                                                   \
 510     if ((z = tomp(*y)) != 0) {                                          \
 511       Py_INCREF(*x);                                                    \
 512       *y = pre##_pywrap(z);                                             \
 513       return (0);                                                       \
 514     }                                                                   \
 515     return (1);                                                         \
 516   }
 517 COERCE(mp, MP)
 518 COERCE(gf, GF)
 519 #undef COERCE
 520
 521 static int mp_pycompare(PyObject *x, PyObject *y)
 522   { return mp_cmp(MP_X(x), MP_X(y)); }
 523
 524 static PyObject *mp_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
 525 {
 526   PyObject *x;
 527   mp *z;
 528   mp_pyobj *zz = 0;
 529   int radix = 0;
 530   char *kwlist[] = { "x", "radix", 0 };
 531
 532   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O|i:new", kwlist, &x, &radix))
 533     goto end;
 534   if (MP_PYCHECK(x)) RETURN_OBJ(x);
 535   if (!good_radix_p(radix, 1)) VALERR("bad radix");
 536   if ((z = mp_frompyobject(x, radix)) == 0) {
 537     PyErr_Format(PyExc_TypeError, "can't convert %.100s to mp",
 538                  x->ob_type->tp_name);
 539     goto end;
 540   }
 541   zz = (mp_pyobj *)ty->tp_alloc(ty, 0);
 542   zz->x = z;
 543 end:
 544   return ((PyObject *)zz);
 545 }
 546
 547 long mphash(mp *x)
 548 {
 549   PyObject *l = mp_topylong(x);
 550   long h = PyObject_Hash(l);
 551   Py_DECREF(l); return (h);
 552 }
 553
 554 static long mp_pyhash(PyObject *me) { return (mphash(MP_X(me))); }
 555
 556 static PyObject *mpmeth_jacobi(PyObject *me, PyObject *arg)
 557 {
 558   mp *y = 0;
 559   PyObject *z = 0;
 560
 561   if (!PyArg_ParseTuple(arg, "O&:jacobi", convmp, &y)) goto end;
 562   z = PyInt_FromLong(mp_jacobi(y, MP_X(me)));
 563 end:
 564   if (y) MP_DROP(y);
 565   return (z);
 566 }
 567
 568 #define BITOP(pre, name, c)                                             \
 569   static PyObject *pre##meth_##name(PyObject *me, PyObject *arg)        \
 570   {                                                                     \
 571     unsigned long i;                                                    \
 572     if (!PyArg_ParseTuple(arg, "O&:" #name, convulong, &i)) return (0); \
 573     return (pre##_pywrap(mp_##name##c(MP_NEW, MP_X(me), i)));           \
 574   }
 575 BITOP(mp, setbit, 2c);
 576 BITOP(mp, clearbit, 2c);
 577 BITOP(gf, setbit, );
 578 BITOP(gf, clearbit, );
 579 #undef BITOP
 580
 581 static PyObject *mpmeth_testbit(PyObject *me, PyObject *arg)
 582 {
 583   unsigned long i;
 584   if (!PyArg_ParseTuple(arg, "O&:testbit", convulong, &i)) return (0);
 585   return (getbool(mp_testbit2c(MP_X(me), i)));
 586 }
 587
 588 static PyObject *gfmeth_testbit(PyObject *me, PyObject *arg)
 589 {
 590   unsigned long i;
 591   if (!PyArg_ParseTuple(arg, "O&:testbit", convulong, &i)) return (0);
 592   return (getbool(mp_testbit(MP_X(me), i)));
 593 }
 594
 595 static PyObject *mpmeth_odd(PyObject *me, PyObject *arg)
 596 {
 597   mp *t;
 598   size_t s;
 599
 600   if (!PyArg_ParseTuple(arg, ":odd")) return (0);
 601   t = mp_odd(MP_NEW, MP_X(me), &s);
 602   return (Py_BuildValue("(lN)", (long)s, mp_pywrap(t)));
 603 }
 604
 605 static PyObject *mpmeth_sqr(PyObject *me, PyObject *arg)
 606 {
 607   if (!PyArg_ParseTuple(arg, ":sqr")) return (0);
 608   return (mp_pywrap(mp_sqr(MP_NEW, MP_X(me))));
 609 }
 610
 611 static PyObject *mpmeth_sqrt(PyObject *me, PyObject *arg)
 612 {
 613   if (!PyArg_ParseTuple(arg, ":sqrt")) return (0);
 614   if (MP_NEGP(MP_X(me))) VALERR("negative root");
 615   return (mp_pywrap(mp_sqrt(MP_NEW, MP_X(me))));
 616 end:
 617   return (0);
 618 }
 619
 620 static PyObject *mpmeth_gcd(PyObject *me, PyObject *arg)
 621 {
 622   mp *y = 0, *zz = MP_NEW;
 623   PyObject *z = 0;
 624
 625   if (!PyArg_ParseTuple(arg, "O&:gcd", convmp, &y)) goto end;
 626   mp_gcd(&zz, 0, 0, MP_X(me), y);
 627   z = mp_pywrap(zz);
 628 end:
 629   if (y) MP_DROP(y);
 630   return (z);
 631 }
 632
 633 static PyObject *mpmeth_gcdx(PyObject *me, PyObject *arg)
 634 {
 635   PyObject *z = 0;
 636   mp *yy = 0, *zz = MP_NEW, *uu = MP_NEW, *vv = MP_NEW;
 637
 638   if (!PyArg_ParseTuple(arg, "O&:gcdx", convmp, &yy)) goto end;
 639   mp_gcd(&zz, &uu, &vv, MP_X(me), yy);
 640   z = Py_BuildValue("(NNN)",
 641                     mp_pywrap(zz), mp_pywrap(uu), mp_pywrap(vv));
 642 end:
 643   if (yy) MP_DROP(yy);
 644   return (z);
 645 }
 646
 647 static PyObject *mpmeth_modinv(PyObject *me, PyObject *arg)
 648 {
 649   PyObject *z = 0;
 650   mp *yy = 0, *zz = MP_NEW;
 651
 652   if (!PyArg_ParseTuple(arg, "O&:modinv", convmp, &yy) ||
 653       (zz = mp_modinv_checked(MP_NEW, yy, MP_X(me))) == 0)
 654     goto end;
 655   z = mp_pywrap(zz);
 656 end:
 657   if (yy) MP_DROP(yy);
 658   return (z);
 659 }
 660
 661 static PyObject *mpmeth_tostring(PyObject *me, PyObject *arg, PyObject *kw)
 662 {
 663   int radix = 10;
 664   char *kwlist[] = { "radix", 0 };
 665   if (!PyArg_ParseTupleAndKeywords(arg, kw, "|i:tostring", kwlist, &radix))
 666     goto end;
 667   if (!good_radix_p(radix, 0)) VALERR("bad radix");
 668   return (mp_topystring(MP_X(me), radix, 0, 0, 0));
 669 end:
 670   return (0);
 671 }
 672
 673 static PyObject *mpmeth_modsqrt(PyObject *me, PyObject *arg)
 674 {
 675   PyObject *z = 0;
 676   mp *yy = 0, *zz = MP_NEW;
 677
 678   if (!PyArg_ParseTuple(arg, "O&:modsqrt", convmp, &yy)) goto end;
 679   if ((zz = mp_modsqrt(MP_NEW, yy, MP_X(me))) == 0)
 680     VALERR("no modular square root");
 681   z = mp_pywrap(zz);
 682 end:
 683   if (yy) MP_DROP(yy);
 684   return (z);
 685 }
 686
 687 static PyObject *mpmeth_leastcongruent(PyObject *me, PyObject *arg)
 688 {
 689   mp *z, *b, *m;
 690   PyObject *rc = 0;
 691
 692   if (!PyArg_ParseTuple(arg, "O&O&:leastcongruent", convmp, &b, convmp, &m))
 693     goto end;
 694   z = mp_leastcongruent(MP_NEW, b, MP_X(me), m);
 695   rc = mp_pywrap(z);
 696 end:
 697   return (rc);
 698 }
 699
 700 #define STOREOP(name, c)                                                \
 701   static PyObject *mpmeth_##name(PyObject *me,                          \
 702                                  PyObject *arg, PyObject *kw)           \
 703   {                                                                     \
 704     long len = -1;                                                      \
 705     char *kwlist[] = { "len", 0 };                                      \
 706     PyObject *rc = 0;                                                   \
 707                                                                         \
 708     if (!PyArg_ParseTupleAndKeywords(arg, kw, "|l:" #name,              \
 709                                     kwlist, &len))                      \
 710       goto end;                                                         \
 711     if (len < 0) {                                                      \
 712       len = mp_octets##c(MP_X(me));                                     \
 713       if (!len) len = 1;                                                \
 714     }                                                                   \
 715     rc = bytestring_pywrap(0, len);                                     \
 716     mp_##name(MP_X(me), PyString_AS_STRING(rc), len);                   \
 717   end:                                                                  \
 718     return (rc);                                                        \
 719   }
 720 STOREOP(storel, )
 721 STOREOP(storeb, )
 722 STOREOP(storel2c, 2c)
 723 STOREOP(storeb2c, 2c)
 724 #undef STOREOP
 725
 726 #define BUFOP(ty, pyty)                                                 \
 727   static PyObject *meth__##pyty##_frombuf(PyObject *me, PyObject *arg)  \
 728   {                                                                     \
 729     buf b;                                                              \
 730     char *p;                                                            \
 731     Py_ssize_t sz;                                                      \
 732     PyObject *rc = 0;                                                   \
 733     mp *x;                                                              \
 734                                                                         \
 735     if (!PyArg_ParseTuple(arg, "Os#:frombuf", &me, &p, &sz)) goto end;  \
 736     buf_init(&b, p, sz);                                                \
 737     if ((x = buf_getmp(&b)) == 0) VALERR("malformed data");             \
 738     rc = Py_BuildValue("(NN)", ty##_pywrap(x),                          \
 739                        bytestring_pywrapbuf(&b));                       \
 740   end:                                                                  \
 741     return (rc);                                                        \
 742   }
 743 BUFOP(mp, MP)
 744 BUFOP(gf, GF)
 745 #undef BUFOP
 746
 747 static PyObject *mpmeth_tobuf(PyObject *me, PyObject *arg)
 748 {
 749   buf b;
 750   PyObject *rc;
 751   mp *x;
 752   size_t n;
 753
 754   if (!PyArg_ParseTuple(arg, ":tobuf")) return (0);
 755   x = MP_X(me);
 756   n = mp_octets(x) + 3;
 757   rc = bytestring_pywrap(0, n);
 758   buf_init(&b, PyString_AS_STRING(rc), n);
 759   buf_putmp(&b, x);
 760   assert(BOK(&b));
 761   _PyString_Resize(&rc, BLEN(&b));
 762   return (rc);
 763 }
 764
 765 static PyObject *mpmeth_primep(PyObject *me, PyObject *arg, PyObject *kw)
 766 {
 767   grand *r = &rand_global;
 768   char *kwlist[] = { "rng", 0 };
 769   PyObject *rc = 0;
 770
 771   if (!PyArg_ParseTupleAndKeywords(arg, kw, "|O&", kwlist, convgrand, &r))
 772     goto end;
 773   rc = getbool(pgen_primep(MP_X(me), r));
 774 end:
 775   return (rc);
 776 }
 777
 778 static PyObject *mpget_nbits(PyObject *me, void *hunoz)
 779   { return (PyInt_FromLong(mp_bits(MP_X(me)))); }
 780
 781 static PyObject *mpget_noctets(PyObject *me, void *hunoz)
 782   { return (PyInt_FromLong(mp_octets(MP_X(me)))); }
 783
 784 static PyObject *mpget_noctets2c(PyObject *me, void *hunoz)
 785   { return (PyInt_FromLong(mp_octets2c(MP_X(me)))); }
 786
 787 static PyGetSetDef mp_pygetset[] = {
 788 #define GETSETNAME(op, func) mp##op##_##func
 789   GET   (nbits,         "X.nbits -> bit length of X")
 790   GET   (noctets,       "X.noctets -> octet length of X")
 791   GET   (noctets2c,     "X.noctets2c -> two's complement octet length of X")
 792 #undef GETSETNAME
 793   { 0 }
 794 };
 795
 796 static PyMethodDef mp_pymethods[] = {
 797 #define METHNAME(func) mpmeth_##func
 798   METH  (jacobi,        "X.jacobi(Y) -> Jacobi symbol (Y/X) (NB inversion!)")
 799   METH  (setbit,        "X.setbit(N) -> X with bit N set")
 800   METH  (clearbit,      "X.clearbit(N) -> X with bit N clear")
 801   METH  (testbit,       "X.testbit(N) -> true/false if bit N set/clear in X")
 802   METH  (odd,           "X.odd() -> S, T where X = 2^S T with T odd")
 803   METH  (sqr,           "X.sqr() -> X^2")
 804   METH  (sqrt,          "X.sqrt() -> largest integer <= sqrt(X)")
 805   METH  (gcd,           "X.gcd(Y) -> gcd(X, Y)")
 806   METH  (gcdx,
 807          "X.gcdx(Y) -> (gcd(X, Y), U, V) with X U + Y V = gcd(X, Y)")
 808   METH  (modinv,        "X.modinv(Y) -> multiplicative inverse of Y mod X")
 809   METH  (modsqrt,       "X.modsqrt(Y) -> square root of Y mod X, if X prime")
 810   METH  (leastcongruent,
 811          "X.leastcongruent(B, M) -> smallest Z >= B with Z == X (mod M)")
 812   KWMETH(primep,        "X.primep(rng = rand) -> true/false if X is prime")
 813   KWMETH(tostring,      "X.tostring(radix = 10) -> STR")
 814   KWMETH(storel,        "X.storel(len = -1) -> little-endian bytes")
 815   KWMETH(storeb,        "X.storeb(len = -1) -> big-endian bytes")
 816   KWMETH(storel2c,
 817          "X.storel2c(len = -1) -> little-endian bytes, two's complement")
 818   KWMETH(storeb2c,
 819          "X.storeb2c(len = -1) -> big-endian bytes, two's complement")
 820   METH  (tobuf,         "X.tobuf() -> buffer format")
 821 #undef METHNAME
 822   { 0 }
 823 };
 824
 825 static PyNumberMethods mp_pynumber = {
 826   mp_pyadd,                             /* @nb_add@ */
 827   mp_pysub,                             /* @nb_subtract@ */
 828   mp_pymul,                             /* @nb_multiply@ */
 829   0,                                    /* @nb_divide@ */
 830   mp_pymod,                             /* @nb_remainder@ */
 831   mp_pydivmod,                          /* @nb_divmod@ */
 832   mp_pyexp,                             /* @nb_power@ */
 833   mp_pyneg,                             /* @nb_negative@ */
 834   mp_pyid,                              /* @nb_positive@ */
 835   mp_pyabs,                             /* @nb_absolute@ */
 836   mp_pynonzerop,                        /* @nb_nonzero@ */
 837   mp_pynot2c,                           /* @nb_invert@ */
 838   mp_pylsl2c,                           /* @nb_lshift@ */
 839   mp_pylsr2c,                           /* @nb_rshift@ */
 840   mp_pyand2c,                           /* @nb_and@ */
 841   mp_pyxor2c,                           /* @nb_xor@ */
 842   mp_pyor2c,                            /* @nb_or@ */
 843   mp_pycoerce,                          /* @nb_coerce@ */
 844   mp_pyint,                             /* @nb_int@ */
 845   mp_pylong,                            /* @nb_long@ */
 846   mp_pyfloat,                           /* @nb_float@ */
 847   mp_pyoct,                             /* @nb_oct@ */
 848   mp_pyhex,                             /* @nb_hex@ */
 849
 850   0,                                    /* @nb_inplace_add@ */
 851   0,                                    /* @nb_inplace_subtract@ */
 852   0,                                    /* @nb_inplace_multiply@ */
 853   0,                                    /* @nb_inplace_divide@ */
 854   0,                                    /* @nb_inplace_remainder@ */
 855   0,                                    /* @nb_inplace_power@ */
 856   0,                                    /* @nb_inplace_lshift@ */
 857   0,                                    /* @nb_inplace_rshift@ */
 858   0,                                    /* @nb_inplace_and@ */
 859   0,                                    /* @nb_inplace_xor@ */
 860   0,                                    /* @nb_inplace_or@ */
 861
 862   mp_pydiv,                             /* @nb_floor_divide@ */
 863   0,                                    /* @nb_true_divide@ */
 864   0,                                    /* @nb_inplace_floor_divide@ */
 865   0,                                    /* @nb_inplace_true_divide@ */
 866 };
 867
 868 static PyTypeObject mp_pytype_skel = {
 869   PyObject_HEAD_INIT(0) 0,              /* Header */
 870   "MP",                                 /* @tp_name@ */
 871   sizeof(mp_pyobj),                     /* @tp_basicsize@ */
 872   0,                                    /* @tp_itemsize@ */
 873
 874   mp_pydealloc,                         /* @tp_dealloc@ */
 875   0,                                    /* @tp_print@ */
 876   0,                                    /* @tp_getattr@ */
 877   0,                                    /* @tp_setattr@ */
 878   mp_pycompare,                         /* @tp_compare@ */
 879   mp_pyrepr,                            /* @tp_repr@ */
 880   &mp_pynumber,                         /* @tp_as_number@ */
 881   0,                                    /* @tp_as_sequence@ */
 882   0,                                    /* @tp_as_mapping@ */
 883   mp_pyhash,                            /* @tp_hash@ */
 884   0,                                    /* @tp_call@ */
 885   mp_pystr,                             /* @tp_str@ */
 886   0,                                    /* @tp_getattro@ */
 887   0,                                    /* @tp_setattro@ */
 888   0,                                    /* @tp_as_buffer@ */
 889   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
 890     Py_TPFLAGS_CHECKTYPES |
 891     Py_TPFLAGS_BASETYPE,
 892
 893   /* @tp_doc@ */
 894 "Multiprecision integers, similar to `long' but more efficient and\n\
 895 versatile.  Support all the standard arithmetic operations.\n\
 896 \n\
 897 Constructor mp(X, radix = R) attempts to convert X to an `mp'.  If\n\
 898 X is a string, it's read in radix-R form, or we look for a prefix\n\
 899 if R = 0.  Other acceptable things are ints and longs.\n\
 900 \n\
 901 Notes:\n\
 902 \n\
 903   * Use `//' for division.  MPs don't have `/' division.",
 904
 905   0,                                    /* @tp_traverse@ */
 906   0,                                    /* @tp_clear@ */
 907   0,                                    /* @tp_richcompare@ */
 908   0,                                    /* @tp_weaklistoffset@ */
 909   0,                                    /* @tp_iter@ */
 910   0,                                    /* @tp_iternext@ */
 911   mp_pymethods,                         /* @tp_methods@ */
 912   0,                                    /* @tp_members@ */
 913   mp_pygetset,                          /* @tp_getset@ */
 914   0,                                    /* @tp_base@ */
 915   0,                                    /* @tp_dict@ */
 916   0,                                    /* @tp_descr_get@ */
 917   0,                                    /* @tp_descr_set@ */
 918   0,                                    /* @tp_dictoffset@ */
 919   0,                                    /* @tp_init@ */
 920   PyType_GenericAlloc,                  /* @tp_alloc@ */
 921   mp_pynew,                             /* @tp_new@ */
 922   0,                                    /* @tp_free@ */
 923   0                                     /* @tp_is_gc@ */
 924 };
 925
 926 static PyObject *meth__MP_fromstring(PyObject *me,
 927                                      PyObject *arg, PyObject *kw)
 928 {
 929   int r = 0;
 930   char *p;
 931   Py_ssize_t len;
 932   PyObject *z = 0;
 933   mp *zz;
 934   mptext_stringctx sc;
 935   char *kwlist[] = { "class", "x", "radix", 0 };
 936
 937   if (!PyArg_ParseTupleAndKeywords(arg, kw, "Os#|i:fromstring",
 938                                    kwlist, &me, &p, &len, &r))
 939     goto end;
 940   if (!good_radix_p(r, 1)) VALERR("bad radix");
 941   sc.buf = p; sc.lim = p + len;
 942   if ((zz = mp_read(MP_NEW, r, &mptext_stringops, &sc)) == 0)
 943     VALERR("bad integer");
 944   z = Py_BuildValue("(Ns#)", mp_pywrap(zz),
 945                     sc.buf, (Py_ssize_t)(sc.lim - sc.buf));
 946 end:
 947   return (z);
 948 }
 949
 950 static PyObject *meth__MP_factorial(PyObject *me, PyObject *arg)
 951 {
 952   unsigned long i;
 953   mp *x;
 954   if (!PyArg_ParseTuple(arg, "OO&:factorial", &me, convulong, &i))
 955     return (0);
 956   x = mp_factorial(i);
 957   return mp_pywrap(x);
 958 }
 959
 960 static PyObject *meth__MP_fibonacci(PyObject *me, PyObject *arg)
 961 {
 962   long i;
 963   mp *x;
 964   if (!PyArg_ParseTuple(arg, "Ol:fibonacci", &me, &i)) return (0);
 965   x = mp_fibonacci(i);
 966   return mp_pywrap(x);
 967 }
 968
 969 #define LOADOP(pre, py, name)                                           \
 970   static PyObject *meth__##py##_##name(PyObject *me, PyObject *arg)     \
 971   {                                                                     \
 972     char *p;                                                            \
 973     Py_ssize_t len;                                                     \
 974     if (!PyArg_ParseTuple(arg, "Os#:" #name, &me, &p, &len)) return (0); \
 975     return (pre##_pywrap(mp_##name(MP_NEW, p, len)));                   \
 976   }
 977 LOADOP(mp, MP, loadl)
 978 LOADOP(mp, MP, loadb)
 979 LOADOP(mp, MP, loadl2c)
 980 LOADOP(mp, MP, loadb2c)
 981 LOADOP(gf, GF, loadl)
 982 LOADOP(gf, GF, loadb)
 983 #undef LOADOP
 984
 985 /*----- Products of small integers ----------------------------------------*/
 986
 987 typedef struct mpmul_pyobj {
 988   PyObject_HEAD
 989   int livep;
 990   mpmul mm;
 991 } mpmul_pyobj;
 992
 993 #define MPMUL_LIVEP(o) (((mpmul_pyobj *)(o))->livep)
 994 #define MPMUL_PY(o) (&((mpmul_pyobj *)(o))->mm)
 995
 996 static void mpmul_pydealloc(PyObject *me)
 997 {
 998   if (MPMUL_LIVEP(me))
 999     mp_drop(mpmul_done(MPMUL_PY(me)));
1000   FREEOBJ(me);
1001 }
1002
1003 static PyObject *mmmeth_factor(PyObject *me, PyObject *arg)
1004 {
1005   PyObject *q, *i;
1006   mp *x;
1007
1008   if (!MPMUL_LIVEP(me)) VALERR("MPMul object invalid");
1009   if (PyTuple_Size(arg) != 1)
1010     i = PyObject_GetIter(arg);
1011   else {
1012     if ((q = PyTuple_GetItem(arg, 0)) == 0) goto end;
1013     if ((i = PyObject_GetIter(q)) == 0) {
1014       PyErr_Clear(); /* that's ok */
1015       i = PyObject_GetIter(arg);
1016     }
1017   }
1018   if (!i) goto end;
1019   while ((q = PyIter_Next(i)) != 0) {
1020     x = getmp(q); Py_DECREF(q); if (!x) {
1021       Py_DECREF(i);
1022       goto end;
1023     }
1024     mpmul_add(MPMUL_PY(me), x);
1025     MP_DROP(x);
1026   }
1027   Py_DECREF(i);
1028   RETURN_ME;
1029 end:
1030   return (0);
1031 }
1032
1033 static PyObject *mmmeth_done(PyObject *me, PyObject *arg)
1034 {
1035   mp *x;
1036
1037   if (!PyArg_ParseTuple(arg, ":done")) goto end;
1038   if (!MPMUL_LIVEP(me)) VALERR("MPMul object invalid");
1039   x = mpmul_done(MPMUL_PY(me));
1040   MPMUL_LIVEP(me) = 0;
1041   return (mp_pywrap(x));
1042 end:
1043   return (0);
1044 }
1045
1046 static PyObject *mpmul_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1047 {
1048   mpmul_pyobj *mm;
1049
1050   if (kw) TYERR("keyword arguments not allowed here");
1051   mm = (mpmul_pyobj *)ty->tp_alloc(ty, 0);
1052   mpmul_init(&mm->mm);
1053   mm->livep = 1;
1054   if (mmmeth_factor((PyObject *)mm, arg) == 0) {
1055     Py_DECREF(mm);
1056     goto end;
1057   }
1058   return ((PyObject *)mm);
1059 end:
1060   return (0);
1061 }
1062
1063 static PyObject *mmget_livep(PyObject *me, void *hunoz)
1064   { return (getbool(MPMUL_LIVEP(me))); }
1065
1066 static PyGetSetDef mpmul_pygetset[] = {
1067 #define GETSETNAME(op, name) mm##op##_##name
1068   GET   (livep,                 "MM.livep -> flag: object still valid?")
1069 #undef GETSETNAME
1070   { 0 }
1071 };
1072
1073 static PyMethodDef mpmul_pymethods[] = {
1074 #define METHNAME(name) mmmeth_##name
1075   METH  (factor,                "MM.factor(ITERABLE) or MM.factor(I, ...)")
1076   METH  (done,                  "MM.done() -> PRODUCT")
1077 #undef METHNAME
1078   { 0 }
1079 };
1080
1081 static PyTypeObject *mpmul_pytype, mpmul_pytype_skel = {
1082   PyObject_HEAD_INIT(0) 0,              /* Header */
1083   "MPMul",                              /* @tp_name@ */
1084   sizeof(mpmul_pyobj),                  /* @tp_basicsize@ */
1085   0,                                    /* @tp_itemsize@ */
1086
1087   mpmul_pydealloc,                      /* @tp_dealloc@ */
1088   0,                                    /* @tp_print@ */
1089   0,                                    /* @tp_getattr@ */
1090   0,                                    /* @tp_setattr@ */
1091   0,                                    /* @tp_compare@ */
1092   0,                                    /* @tp_repr@ */
1093   0,                                    /* @tp_as_number@ */
1094   0,                                    /* @tp_as_sequence@ */
1095   0,                                    /* @tp_as_mapping@ */
1096   0,                                    /* @tp_hash@ */
1097   0,                                    /* @tp_call@ */
1098   0,                                    /* @tp_str@ */
1099   0,                                    /* @tp_getattro@ */
1100   0,                                    /* @tp_setattro@ */
1101   0,                                    /* @tp_as_buffer@ */
1102   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1103     Py_TPFLAGS_BASETYPE,
1104
1105   /* @tp_doc@ */
1106 "An object for multiplying many small integers.",
1107
1108   0,                                    /* @tp_traverse@ */
1109   0,                                    /* @tp_clear@ */
1110   0,                                    /* @tp_richcompare@ */
1111   0,                                    /* @tp_weaklistoffset@ */
1112   0,                                    /* @tp_iter@ */
1113   0,                                    /* @tp_iternext@ */
1114   mpmul_pymethods,                      /* @tp_methods@ */
1115   0,                                    /* @tp_members@ */
1116   mpmul_pygetset,                       /* @tp_getset@ */
1117   0,                                    /* @tp_base@ */
1118   0,                                    /* @tp_dict@ */
1119   0,                                    /* @tp_descr_get@ */
1120   0,                                    /* @tp_descr_set@ */
1121   0,                                    /* @tp_dictoffset@ */
1122   0,                                    /* @tp_init@ */
1123   PyType_GenericAlloc,                  /* @tp_alloc@ */
1124   mpmul_pynew,                          /* @tp_new@ */
1125   0,                                    /* @tp_free@ */
1126   0                                     /* @tp_is_gc@ */
1127 };
1128
1129 /*----- Montgomery reduction ----------------------------------------------*/
1130
1131 typedef struct mpmont_pyobj {
1132   PyObject_HEAD
1133   mpmont mm;
1134 } mpmont_pyobj;
1135
1136 #define MPMONT_PY(o) (&((mpmont_pyobj *)(o))->mm)
1137
1138 static PyObject *mmmeth_int(PyObject *me, PyObject *arg)
1139 {
1140   PyObject *z = 0;
1141   mp *yy = 0;
1142   mpmont *mm = MPMONT_PY(me);
1143
1144   if (!PyArg_ParseTuple(arg, "O&:in", convmp, &yy))
1145     goto end;
1146   mp_div(0, &yy, yy, mm->m);
1147   z = mp_pywrap(mpmont_mul(mm, MP_NEW, yy, mm->r2));
1148 end:
1149   if (yy) MP_DROP(yy);
1150   return (z);
1151 }
1152
1153 static PyObject *mmmeth_mul(PyObject *me, PyObject *arg)
1154 {
1155   PyObject *rc = 0;
1156   mp *yy = 0, *zz = 0;
1157
1158   if (!PyArg_ParseTuple(arg, "O&O&:mul", convmp, &yy, convmp, &zz))
1159     goto end;
1160   rc = mp_pywrap(mpmont_mul(MPMONT_PY(me), MP_NEW, yy, zz));
1161 end:
1162   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1163   return (rc);
1164 }
1165
1166 static PyObject *mmmeth_exp(PyObject *me, PyObject *arg)
1167 {
1168   PyObject *rc = 0;
1169   mp *yy = 0, *zz = 0;
1170
1171   if (!PyArg_ParseTuple(arg, "O&O&:exp", convmp, &yy, convmp, &zz))
1172     goto end;
1173   if (MP_NEGP(zz)) {
1174     if ((yy = mp_modinv_checked(yy, yy, MPMONT_PY(me)->m)) == 0) goto end;
1175     zz = mp_neg(zz, zz);
1176   }
1177   rc = mp_pywrap(mpmont_exp(MPMONT_PY(me), MP_NEW, yy, zz));
1178 end:
1179   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1180   return (rc);
1181 }
1182
1183 static PyObject *mmmeth_expr(PyObject *me, PyObject *arg)
1184 {
1185   PyObject *rc = 0;
1186   mp *yy = 0, *zz = 0;
1187
1188   if (!PyArg_ParseTuple(arg, "O&O&:expr", convmp, &yy, convmp, &zz))
1189     goto end;
1190   if (MP_NEGP(zz)) {
1191     yy = mpmont_reduce(MPMONT_PY(me), yy, yy);
1192     if ((yy = mp_modinv_checked(yy, yy, MPMONT_PY(me)->m)) == 0) goto end;
1193     yy = mpmont_mul(MPMONT_PY(me), yy, yy, MPMONT_PY(me)->r2);
1194     zz = mp_neg(zz, zz);
1195   }
1196   rc = mp_pywrap(mpmont_expr(MPMONT_PY(me), MP_NEW, yy, zz));
1197 end:
1198   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1199   return (rc);
1200 }
1201
1202 static PyObject *mm_mexpr_id(PyObject *me)
1203   { return mp_pywrap(MP_COPY(MPMONT_PY(me)->r)); }
1204
1205 static int mm_mexpr_fill(void *p, PyObject *me, PyObject *x, PyObject *y)
1206 {
1207   mp *xx = 0, *yy = 0;
1208   mp_expfactor *f = p;
1209   mpmont *mm = MPMONT_PY(me);
1210
1211   if ((xx = getmp(x)) == 0 || (yy = getmp(y)) == 0)
1212     goto fail;
1213   if (MP_NEGP(yy)) {
1214     xx = mpmont_reduce(mm, xx, xx);
1215     if ((xx = mp_modinv_checked(xx, xx, yy)) == 0)
1216       goto fail;
1217     xx = mpmont_mul(mm, xx, xx, mm->r2);
1218     yy = mp_neg(yy, yy);
1219   }
1220   f->base = xx;
1221   f->exp = yy;
1222   return (0);
1223
1224 fail:
1225   mp_drop(xx); mp_drop(yy);
1226   return (-1);
1227 }
1228
1229 static PyObject *mm_mexpr(PyObject *me, void *v, int n)
1230   { return mp_pywrap(mpmont_mexpr(MPMONT_PY(me), MP_NEW, v, n)); }
1231
1232 static void mp_mexp_drop(void *p)
1233 {
1234   mp_expfactor *f = p;
1235   mp_drop(f->base);
1236   mp_drop(f->exp);
1237 }
1238
1239 static PyObject *mmmeth_mexpr(PyObject *me, PyObject *arg)
1240 {
1241   return mexp_common(me, arg, sizeof(mp_expfactor),
1242                      mm_mexpr_id, mm_mexpr_fill, mm_mexpr, mp_mexp_drop);
1243 }
1244
1245 static PyObject *mp_mexp_id(PyObject *me)
1246   { return mp_pywrap(MP_ONE); }
1247
1248 static int mp_mexp_fill(void *p, PyObject *me, PyObject *x, PyObject *y)
1249 {
1250   mp *xx = 0, *yy = 0;
1251   mp_expfactor *f = p;
1252
1253   if ((xx = getmp(x)) == 0 || (yy = getmp(y)) == 0)
1254     goto fail;
1255   if (MP_NEGP(yy)) {
1256     if ((xx = mp_modinv_checked(xx, xx, yy)) == 0)
1257       goto fail;
1258     yy = mp_neg(yy, yy);
1259   }
1260   f->base = xx;
1261   f->exp = yy;
1262   return (0);
1263
1264 fail:
1265   mp_drop(xx); mp_drop(yy);
1266   return (-1);
1267 }
1268
1269 static PyObject *mm_mexp(PyObject *me, void *v, int n)
1270   { return mp_pywrap(mpmont_mexp(MPMONT_PY(me), MP_NEW, v, n)); }
1271
1272 static PyObject *mmmeth_mexp(PyObject *me, PyObject *arg)
1273 {
1274   return mexp_common(me, arg, sizeof(mp_expfactor),
1275                      mp_mexp_id, mp_mexp_fill, mm_mexp, mp_mexp_drop);
1276 }
1277
1278 #define mmmeth_ext mmmeth_reduce
1279 static PyObject *mmmeth_reduce(PyObject *me, PyObject *arg)
1280 {
1281   PyObject *z = 0;
1282   mp *yy = 0;
1283
1284   if (!PyArg_ParseTuple(arg, "O&", convmp, &yy)) goto end;
1285   z = mp_pywrap(mpmont_reduce(MPMONT_PY(me), MP_NEW, yy));
1286 end:
1287   return (z);
1288 }
1289
1290 static void mpmont_pydealloc(PyObject *me)
1291 {
1292   mpmont_destroy(MPMONT_PY(me));
1293   FREEOBJ(me);
1294 }
1295
1296 static PyObject *mpmont_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1297 {
1298   mpmont_pyobj *mm = 0;
1299   char *kwlist[] = { "m", 0 };
1300   mp *xx = 0;
1301
1302   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convmp, &xx))
1303     goto end;
1304   if (!MP_POSP(xx) || !MP_ODDP(xx)) VALERR("m must be positive and odd");
1305   mm = (mpmont_pyobj *)ty->tp_alloc(ty, 0);
1306   mpmont_create(&mm->mm, xx);
1307 end:
1308   if (xx) MP_DROP(xx);
1309   return ((PyObject *)mm);
1310 }
1311
1312 static PyObject *mmget_m(PyObject *me, void *hunoz)
1313   { return (mp_pywrap(MP_COPY(MPMONT_PY(me)->m))); }
1314
1315 static PyObject *mmget_r(PyObject *me, void *hunoz)
1316   { return (mp_pywrap(MP_COPY(MPMONT_PY(me)->r))); }
1317
1318 static PyObject *mmget_r2(PyObject *me, void *hunoz)
1319   { return (mp_pywrap(MP_COPY(MPMONT_PY(me)->r2))); }
1320
1321 static PyGetSetDef mpmont_pygetset[] = {
1322 #define GETSETNAME(op, name) mm##op##_##name
1323   GET   (m,             "M.m -> modulus for reduction")
1324   GET   (r,             "M.r -> multiplicative identity")
1325   GET   (r2,            "M.r2 -> M.r^2, Montgomerization factor")
1326 #undef GETSETNAME
1327   { 0 }
1328 };
1329
1330 static PyMethodDef mpmont_pymethods[] = {
1331 #define METHNAME(name) mmmeth_##name
1332   METH  (int,           "M.int(X) -> XR")
1333   METH  (mul,           "M.mul(XR, YR) -> ZR where Z = X Y")
1334   METH  (expr,          "M.expr(XR, N) -> ZR where Z = X^N mod M.m")
1335   METH  (mexpr,         "\
1336 M.mexpr([(XR0, N0), (XR1, N1), ...]) = ZR where Z = X0^N0 X1^N1 ... mod M.m\n\
1337 \t(the list may be flattened if this more convenient.)")
1338   METH  (reduce,        "M.reduce(XR) -> X")
1339   METH  (ext,           "M.ext(XR) -> X")
1340   METH  (exp,           "M.exp(X, N) -> X^N mod M.m")
1341   METH  (mexp,          "\
1342 M.mexp([(X0, N0), (X1, N1), ...]) = X0^N0 X1^N1 ... mod M.m\n\
1343 \t(the list may be flattened if this more convenient.)")
1344 #undef METHNAME
1345   { 0 }
1346 };
1347
1348 static PyTypeObject *mpmont_pytype, mpmont_pytype_skel = {
1349   PyObject_HEAD_INIT(0) 0,              /* Header */
1350   "MPMont",                             /* @tp_name@ */
1351   sizeof(mpmont_pyobj),                 /* @tp_basicsize@ */
1352   0,                                    /* @tp_itemsize@ */
1353
1354   mpmont_pydealloc,                     /* @tp_dealloc@ */
1355   0,                                    /* @tp_print@ */
1356   0,                                    /* @tp_getattr@ */
1357   0,                                    /* @tp_setattr@ */
1358   0,                                    /* @tp_compare@ */
1359   0,                                    /* @tp_repr@ */
1360   0,                                    /* @tp_as_number@ */
1361   0,                                    /* @tp_as_sequence@ */
1362   0,                                    /* @tp_as_mapping@ */
1363   0,                                    /* @tp_hash@ */
1364   0,                                    /* @tp_call@ */
1365   0,                                    /* @tp_str@ */
1366   0,                                    /* @tp_getattro@ */
1367   0,                                    /* @tp_setattro@ */
1368   0,                                    /* @tp_as_buffer@ */
1369   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1370     Py_TPFLAGS_BASETYPE,
1371
1372   /* @tp_doc@ */
1373 "A Montgomery reduction context.",
1374
1375   0,                                    /* @tp_traverse@ */
1376   0,                                    /* @tp_clear@ */
1377   0,                                    /* @tp_richcompare@ */
1378   0,                                    /* @tp_weaklistoffset@ */
1379   0,                                    /* @tp_iter@ */
1380   0,                                    /* @tp_iternext@ */
1381   mpmont_pymethods,                     /* @tp_methods@ */
1382   0,                                    /* @tp_members@ */
1383   mpmont_pygetset,                      /* @tp_getset@ */
1384   0,                                    /* @tp_base@ */
1385   0,                                    /* @tp_dict@ */
1386   0,                                    /* @tp_descr_get@ */
1387   0,                                    /* @tp_descr_set@ */
1388   0,                                    /* @tp_dictoffset@ */
1389   0,                                    /* @tp_init@ */
1390   PyType_GenericAlloc,                  /* @tp_alloc@ */
1391   mpmont_pynew,                         /* @tp_new@ */
1392   0,                                    /* @tp_free@ */
1393   0                                     /* @tp_is_gc@ */
1394 };
1395
1396 /*----- Barrett reduction -------------------------------------------------*/
1397
1398 typedef struct mpbarrett_pyobj {
1399   PyObject_HEAD
1400   mpbarrett mb;
1401 } mpbarrett_pyobj;
1402
1403 #define MPBARRETT_PY(o) (&((mpbarrett_pyobj *)(o))->mb)
1404
1405 static PyObject *mbmeth_exp(PyObject *me, PyObject *arg)
1406 {
1407   PyObject *rc = 0;
1408   mp *yy = 0, *zz = 0;
1409
1410   if (!PyArg_ParseTuple(arg, "O&O&:exp", convmp, &yy, convmp, &zz))
1411     goto end;
1412   if (MP_NEGP(zz)) {
1413     if ((yy = mp_modinv_checked(yy, yy, MPBARRETT_PY(me)->m)) == 0) goto end;
1414     zz = mp_neg(zz, zz);
1415   }
1416   rc = mp_pywrap(mpbarrett_exp(MPBARRETT_PY(me), MP_NEW, yy, zz));
1417 end:
1418   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1419   return (rc);
1420 }
1421
1422 static PyObject *mb_mexp(PyObject *me, void *v, int n)
1423   { return mp_pywrap(mpbarrett_mexp(MPBARRETT_PY(me), MP_NEW, v, n)); }
1424
1425 static PyObject *mbmeth_mexp(PyObject *me, PyObject *arg)
1426 {
1427   return mexp_common(me, arg, sizeof(mp_expfactor),
1428                      mp_mexp_id, mp_mexp_fill, mb_mexp, mp_mexp_drop);
1429 }
1430
1431 static PyObject *mbmeth_reduce(PyObject *me, PyObject *arg)
1432 {
1433   PyObject *z = 0;
1434   mp *yy = 0;
1435
1436   if (!PyArg_ParseTuple(arg, "O&:reduce", convmp, &yy))
1437     goto end;
1438   z = mp_pywrap(mpbarrett_reduce(MPBARRETT_PY(me), MP_NEW, yy));
1439 end:
1440   return (z);
1441 }
1442
1443 static void mpbarrett_pydealloc(PyObject *me)
1444 {
1445   mpbarrett_destroy(MPBARRETT_PY(me));
1446   FREEOBJ(me);
1447 }
1448
1449 static PyObject *mpbarrett_pynew(PyTypeObject *ty,
1450                                  PyObject *arg, PyObject *kw)
1451 {
1452   mpbarrett_pyobj *mb = 0;
1453   char *kwlist[] = { "m", 0 };
1454   mp *xx = 0;
1455
1456   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convmp, &xx))
1457     goto end;
1458   if (!MP_POSP(xx)) VALERR("m must be positive");
1459   mb = (mpbarrett_pyobj *)ty->tp_alloc(ty, 0);
1460   mpbarrett_create(&mb->mb, xx);
1461 end:
1462   if (xx) MP_DROP(xx);
1463   return ((PyObject *)mb);
1464 }
1465
1466 static PyObject *mbget_m(PyObject *me, void *hunoz)
1467   { return (mp_pywrap(MP_COPY(MPBARRETT_PY(me)->m))); }
1468
1469 static PyGetSetDef mpbarrett_pygetset[] = {
1470 #define GETSETNAME(op, name) mb##op##_##name
1471   GET   (m,             "B.m -> modulus for reduction")
1472 #undef GETSETNAME
1473   { 0 }
1474 };
1475
1476 static PyMethodDef mpbarrett_pymethods[] = {
1477 #define METHNAME(name) mbmeth_##name
1478   METH  (reduce,        "B.reduce(X) -> X mod B.m")
1479   METH  (exp,           "B.exp(X, N) -> X^N mod B.m")
1480   METH  (mexp,          "\
1481 B.mexp([(X0, N0), (X1, N1), ...]) = X0^N0 X1^N1 ... mod B.m\n\
1482 \t(the list may be flattened if this more convenient.)")
1483 #undef METHNAME
1484   { 0 }
1485 };
1486
1487 static PyTypeObject *mpbarrett_pytype, mpbarrett_pytype_skel = {
1488   PyObject_HEAD_INIT(0) 0,              /* Header */
1489   "MPBarrett",                          /* @tp_name@ */
1490   sizeof(mpbarrett_pyobj),              /* @tp_basicsize@ */
1491   0,                                    /* @tp_itemsize@ */
1492
1493   mpbarrett_pydealloc,                  /* @tp_dealloc@ */
1494   0,                                    /* @tp_print@ */
1495   0,                                    /* @tp_getattr@ */
1496   0,                                    /* @tp_setattr@ */
1497   0,                                    /* @tp_compare@ */
1498   0,                                    /* @tp_repr@ */
1499   0,                                    /* @tp_as_number@ */
1500   0,                                    /* @tp_as_sequence@ */
1501   0,                                    /* @tp_as_mapping@ */
1502   0,                                    /* @tp_hash@ */
1503   0,                                    /* @tp_call@ */
1504   0,                                    /* @tp_str@ */
1505   0,                                    /* @tp_getattro@ */
1506   0,                                    /* @tp_setattro@ */
1507   0,                                    /* @tp_as_buffer@ */
1508   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1509     Py_TPFLAGS_BASETYPE,
1510
1511   /* @tp_doc@ */
1512 "A Barrett reduction context.",
1513
1514   0,                                    /* @tp_traverse@ */
1515   0,                                    /* @tp_clear@ */
1516   0,                                    /* @tp_richcompare@ */
1517   0,                                    /* @tp_weaklistoffset@ */
1518   0,                                    /* @tp_iter@ */
1519   0,                                    /* @tp_iternext@ */
1520   mpbarrett_pymethods,                  /* @tp_methods@ */
1521   0,                                    /* @tp_members@ */
1522   mpbarrett_pygetset,                   /* @tp_getset@ */
1523   0,                                    /* @tp_base@ */
1524   0,                                    /* @tp_dict@ */
1525   0,                                    /* @tp_descr_get@ */
1526   0,                                    /* @tp_descr_set@ */
1527   0,                                    /* @tp_dictoffset@ */
1528   0,                                    /* @tp_init@ */
1529   PyType_GenericAlloc,                  /* @tp_alloc@ */
1530   mpbarrett_pynew,                      /* @tp_new@ */
1531   0,                                    /* @tp_free@ */
1532   0                                     /* @tp_is_gc@ */
1533 };
1534
1535 /*----- Nice prime reduction ----------------------------------------------*/
1536
1537 typedef struct mpreduce_pyobj {
1538   PyObject_HEAD
1539   mpreduce mr;
1540 } mpreduce_pyobj;
1541
1542 #define MPREDUCE_PY(o) (&((mpreduce_pyobj *)(o))->mr)
1543
1544 static PyObject *mrmeth_exp(PyObject *me, PyObject *arg)
1545 {
1546   PyObject *rc = 0;
1547   mp *yy = 0, *zz = 0;
1548
1549   if (!PyArg_ParseTuple(arg, "O&O&:exp", convmp, &yy, convmp, &zz))
1550     goto end;
1551   if (MP_NEGP(zz)) {
1552     if ((yy = mp_modinv_checked(yy, yy, MPREDUCE_PY(me)->p)) == 0) goto end;
1553     zz = mp_neg(zz, zz);
1554   }
1555   rc = mp_pywrap(mpreduce_exp(MPREDUCE_PY(me), MP_NEW, yy, zz));
1556 end:
1557   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1558   return (rc);
1559 }
1560
1561 static PyObject *mrmeth_reduce(PyObject *me, PyObject *arg)
1562 {
1563   PyObject *z = 0;
1564   mp *yy = 0;
1565
1566   if (!PyArg_ParseTuple(arg, "O&:reduce", convmp, &yy)) goto end;
1567   z = mp_pywrap(mpreduce_do(MPREDUCE_PY(me), MP_NEW, yy));
1568 end:
1569   return (z);
1570 }
1571
1572 static void mpreduce_pydealloc(PyObject *me)
1573 {
1574   mpreduce_destroy(MPREDUCE_PY(me));
1575   FREEOBJ(me);
1576 }
1577
1578 static PyObject *mpreduce_pynew(PyTypeObject *ty,
1579                                  PyObject *arg, PyObject *kw)
1580 {
1581   mpreduce_pyobj *mr = 0;
1582   mpreduce r;
1583   char *kwlist[] = { "m", 0 };
1584   mp *xx = 0;
1585
1586   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convmp, &xx))
1587     goto end;
1588   if (!MP_POSP(xx)) VALERR("m must be positive");
1589   if (mpreduce_create(&r, xx)) VALERR("bad modulus (must be 2^k - ...)");
1590   mr = (mpreduce_pyobj *)ty->tp_alloc(ty, 0);
1591   mr->mr = r;
1592 end:
1593   if (xx) MP_DROP(xx);
1594   return ((PyObject *)mr);
1595 }
1596
1597 static PyObject *mrget_m(PyObject *me, void *hunoz)
1598   { return (mp_pywrap(MP_COPY(MPREDUCE_PY(me)->p))); }
1599
1600 static PyGetSetDef mpreduce_pygetset[] = {
1601 #define GETSETNAME(op, name) mr##op##_##name
1602   GET   (m,             "R.m -> modulus for reduction")
1603 #undef GETSETNAME
1604   { 0 }
1605 };
1606
1607 static PyMethodDef mpreduce_pymethods[] = {
1608 #define METHNAME(name) mrmeth_##name
1609   METH  (reduce,        "R.reduce(X) -> X mod B.m")
1610   METH  (exp,           "R.exp(X, N) -> X^N mod B.m")
1611 #undef METHNAME
1612   { 0 }
1613 };
1614
1615 static PyTypeObject *mpreduce_pytype, mpreduce_pytype_skel = {
1616   PyObject_HEAD_INIT(0) 0,              /* Header */
1617   "MPReduce",                           /* @tp_name@ */
1618   sizeof(mpreduce_pyobj),               /* @tp_basicsize@ */
1619   0,                                    /* @tp_itemsize@ */
1620
1621   mpreduce_pydealloc,                   /* @tp_dealloc@ */
1622   0,                                    /* @tp_print@ */
1623   0,                                    /* @tp_getattr@ */
1624   0,                                    /* @tp_setattr@ */
1625   0,                                    /* @tp_compare@ */
1626   0,                                    /* @tp_repr@ */
1627   0,                                    /* @tp_as_number@ */
1628   0,                                    /* @tp_as_sequence@ */
1629   0,                                    /* @tp_as_mapping@ */
1630   0,                                    /* @tp_hash@ */
1631   0,                                    /* @tp_call@ */
1632   0,                                    /* @tp_str@ */
1633   0,                                    /* @tp_getattro@ */
1634   0,                                    /* @tp_setattro@ */
1635   0,                                    /* @tp_as_buffer@ */
1636   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1637     Py_TPFLAGS_BASETYPE,
1638
1639   /* @tp_doc@ */
1640 "A reduction context for reduction modulo primes of special form.",
1641
1642   0,                                    /* @tp_traverse@ */
1643   0,                                    /* @tp_clear@ */
1644   0,                                    /* @tp_richcompare@ */
1645   0,                                    /* @tp_weaklistoffset@ */
1646   0,                                    /* @tp_iter@ */
1647   0,                                    /* @tp_iternext@ */
1648   mpreduce_pymethods,                   /* @tp_methods@ */
1649   0,                                    /* @tp_members@ */
1650   mpreduce_pygetset,                    /* @tp_getset@ */
1651   0,                                    /* @tp_base@ */
1652   0,                                    /* @tp_dict@ */
1653   0,                                    /* @tp_descr_get@ */
1654   0,                                    /* @tp_descr_set@ */
1655   0,                                    /* @tp_dictoffset@ */
1656   0,                                    /* @tp_init@ */
1657   PyType_GenericAlloc,                  /* @tp_alloc@ */
1658   mpreduce_pynew,                       /* @tp_new@ */
1659   0,                                    /* @tp_free@ */
1660   0                                     /* @tp_is_gc@ */
1661 };
1662
1663 /*----- Chinese Remainder Theorem solution --------------------------------*/
1664
1665 typedef struct mpcrt_pyobj {
1666   PyObject_HEAD
1667   mpcrt c;
1668 } mpcrt_pyobj;
1669
1670 #define MPCRT_PY(o) (&((mpcrt_pyobj *)(o))->c)
1671
1672 static PyObject *mcmeth_solve(PyObject *me, PyObject *arg)
1673 {
1674   mpcrt *c = MPCRT_PY(me);
1675   PyObject *q = 0, *x, *z = 0;
1676   mp *xx;
1677   mp **v = 0;
1678   int i = 0, n = c->k;
1679
1680   Py_INCREF(me);
1681   if (PyTuple_Size(arg) == n)
1682     q = arg;
1683   else if (!PyArg_ParseTuple(arg, "O:solve", &q))
1684     goto end;
1685   Py_INCREF(q);
1686   if (!PySequence_Check(q)) TYERR("want a sequence of residues");
1687   if (PySequence_Size(q) != n) VALERR("residue count mismatch");
1688   v = xmalloc(n * sizeof(*v));
1689   for (i = 0; i < n; i++) {
1690     if ((x = PySequence_GetItem(q, i)) == 0) goto end;
1691     xx = getmp(x); Py_DECREF(x); if (!xx) goto end;
1692     v[i] = xx;
1693   }
1694   z = mp_pywrap(mpcrt_solve(c, MP_NEW, v));
1695 end:
1696   if (v) {
1697     n = i;
1698     for (i = 0; i < n; i++)
1699       MP_DROP(v[i]);
1700     xfree(v);
1701   }
1702   Py_DECREF(me);
1703   Py_XDECREF(q);
1704   return (z);
1705 }
1706
1707 static void mpcrt_pydealloc(PyObject *me)
1708 {
1709   mpcrt *c = MPCRT_PY(me);
1710   mpcrt_destroy(c);
1711   xfree(c->v);
1712 }
1713
1714 static PyObject *mpcrt_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1715 {
1716   mpcrt_mod *v = 0;
1717   int n, i = 0;
1718   char *kwlist[] = { "mv", 0 };
1719   PyObject *q = 0, *x;
1720   mp *xx;
1721   mpcrt_pyobj *c = 0;
1722
1723   if (PyTuple_Size(arg) > 1)
1724     q = arg;
1725   else if (!PyArg_ParseTupleAndKeywords(arg, kw, "O:new", kwlist, &q))
1726     goto end;
1727   Py_INCREF(q);
1728   if (!PySequence_Check(q)) TYERR("want a sequence of moduli");
1729   n = PySequence_Size(q);
1730   if (PyErr_Occurred()) goto end;
1731   if (!n) VALERR("want at least one modulus");
1732   v = xmalloc(n * sizeof(*v));
1733   for (i = 0; i < n; i++) {
1734     if ((x = PySequence_GetItem(q, i)) == 0) goto end;
1735     xx = getmp(x); Py_DECREF(x); if (!xx) goto end;
1736     v[i].m = xx; v[i].n = 0; v[i].ni = 0; v[i].nni = 0;
1737   }
1738   c = (mpcrt_pyobj *)ty->tp_alloc(ty, 0);
1739   mpcrt_create(&c->c, v, n, 0);
1740   Py_DECREF(q);
1741   return ((PyObject *)c);
1742
1743 end:
1744   if (v) {
1745     n = i;
1746     for (i = 0; i < n; i++)
1747       MP_DROP(v[i].m);
1748     xfree(v);
1749   }
1750   Py_XDECREF(q);
1751   return (0);
1752 }
1753
1754 static PyObject *mcget_product(PyObject *me, void *hunoz)
1755   { return (mp_pywrap(MP_COPY(MPCRT_PY(me)->mb.m))); }
1756
1757 static PyObject *mcget_moduli(PyObject *me, void *hunoz)
1758 {
1759   int i;
1760   PyObject *q;
1761   mpcrt *c = MPCRT_PY(me);
1762
1763   if ((q = PyList_New(c->k)) == 0) return (0);
1764   for (i = 0; i < c->k; i++)
1765     PyList_SetItem(q, i, mp_pywrap(c->v[i].m));
1766   return (q);
1767 }
1768
1769 static PyGetSetDef mpcrt_pygetset[] = {
1770 #define GETSETNAME(op, name) mc##op##_##name
1771   GET   (product,       "C.product -> product of moduli")
1772   GET   (moduli,        "C.moduli -> list of individual moduli")
1773 #undef GETSETNAME
1774   { 0 }
1775 };
1776
1777 static PyMethodDef mpcrt_pymethods[] = {
1778 #define METHNAME(name) mcmeth_##name
1779   METH  (solve,         "C.solve([R0, R1]) -> X mod C.product")
1780 #undef METHNAME
1781   { 0 }
1782 };
1783
1784 static PyTypeObject *mpcrt_pytype, mpcrt_pytype_skel = {
1785   PyObject_HEAD_INIT(0) 0,              /* Header */
1786   "MPCRT",                              /* @tp_name@ */
1787   sizeof(mpcrt_pyobj),                  /* @tp_basicsize@ */
1788   0,                                    /* @tp_itemsize@ */
1789
1790   mpcrt_pydealloc,                      /* @tp_dealloc@ */
1791   0,                                    /* @tp_print@ */
1792   0,                                    /* @tp_getattr@ */
1793   0,                                    /* @tp_setattr@ */
1794   0,                                    /* @tp_compare@ */
1795   0,                                    /* @tp_repr@ */
1796   0,                                    /* @tp_as_number@ */
1797   0,                                    /* @tp_as_sequence@ */
1798   0,                                    /* @tp_as_mapping@ */
1799   0,                                    /* @tp_hash@ */
1800   0,                                    /* @tp_call@ */
1801   0,                                    /* @tp_str@ */
1802   0,                                    /* @tp_getattro@ */
1803   0,                                    /* @tp_setattro@ */
1804   0,                                    /* @tp_as_buffer@ */
1805   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1806     Py_TPFLAGS_BASETYPE,
1807
1808   /* @tp_doc@ */
1809 "A context for the solution of Chinese Remainder Theorem problems.",
1810
1811   0,                                    /* @tp_traverse@ */
1812   0,                                    /* @tp_clear@ */
1813   0,                                    /* @tp_richcompare@ */
1814   0,                                    /* @tp_weaklistoffset@ */
1815   0,                                    /* @tp_iter@ */
1816   0,                                    /* @tp_iternext@ */
1817   mpcrt_pymethods,                      /* @tp_methods@ */
1818   0,                                    /* @tp_members@ */
1819   mpcrt_pygetset,                       /* @tp_getset@ */
1820   0,                                    /* @tp_base@ */
1821   0,                                    /* @tp_dict@ */
1822   0,                                    /* @tp_descr_get@ */
1823   0,                                    /* @tp_descr_set@ */
1824   0,                                    /* @tp_dictoffset@ */
1825   0,                                    /* @tp_init@ */
1826   PyType_GenericAlloc,                  /* @tp_alloc@ */
1827   mpcrt_pynew,                          /* @tp_new@ */
1828   0,                                    /* @tp_free@ */
1829   0                                     /* @tp_is_gc@ */
1830 };
1831
1832 /*----- Binary polynomials ------------------------------------------------*/
1833
1834 static PyObject *gf_pyrepr(PyObject *o)
1835   { return mp_topystring(MP_X(o), 16, "GF(", "0x", "L)"); }
1836
1837 static PyObject *gf_pyrichcompare(PyObject *x, PyObject *y, int op)
1838 {
1839   mp *xx, *yy;
1840   int xl, yl;
1841   int b;
1842
1843   if (gfbinop(x, y, &xx, &yy)) RETURN_NOTIMPL;
1844   switch (op) {
1845     case Py_EQ: b = MP_EQ(xx, yy); break;
1846     case Py_NE: b = !MP_EQ(xx, yy); break;
1847     default:
1848       xl = mp_bits(xx);
1849       yl = mp_bits(yy);
1850       switch (op) {
1851         case Py_LT: b = xl < yl; break;
1852         case Py_LE: b = xl <= yl; break;
1853         case Py_GT: b = xl > yl; break;
1854         case Py_GE: b = xl >= yl; break;
1855         default: abort();
1856       }
1857       break;
1858   }
1859   MP_DROP(xx); MP_DROP(yy);
1860   return (getbool(b));
1861 }
1862
1863 static PyObject *gf_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1864 {
1865   PyObject *x;
1866   mp *z;
1867   mp_pyobj *zz = 0;
1868   int radix = 0;
1869   char *kwlist[] = { "x", "radix", 0 };
1870
1871   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O|i:gf", kwlist, &x, &radix))
1872     goto end;
1873   if (GF_PYCHECK(x)) RETURN_OBJ(x);
1874   if (!good_radix_p(radix, 1)) VALERR("radix out of range");
1875   if ((z = mp_frompyobject(x, radix)) == 0) {
1876     PyErr_Format(PyExc_TypeError, "can't convert %.100s to gf",
1877                  x->ob_type->tp_name);
1878     goto end;
1879   }
1880   if (MP_NEGP(z)) {
1881     MP_DROP(z);
1882     VALERR("gf cannot be negative");
1883   }
1884   zz = (mp_pyobj *)ty->tp_alloc(ty, 0);
1885   zz->x = z;
1886 end:
1887   return ((PyObject *)zz);
1888 }
1889
1890 static long gf_pyhash(PyObject *me)
1891 {
1892   long i = mp_tolong(MP_X(me));
1893   i ^= 0xc7ecd67c; /* random perturbance */
1894   if (i == -1)
1895     i = -2;
1896   return (i);
1897 }
1898
1899 static PyObject *gf_pyexp(PyObject *x, PyObject *y, PyObject *z)
1900 {
1901   mp *xx = 0, *yy = 0, *zz = 0;
1902   mp *r = 0;
1903   PyObject *rc = 0;
1904
1905   if ((xx = tomp(x)) == 0 || (yy = tomp(y)) == 0 ||
1906       (z && z != Py_None && (zz = tomp(z)) == 0)) {
1907     mp_drop(xx); mp_drop(yy); mp_drop(zz);
1908     RETURN_NOTIMPL;
1909   }
1910   if (!z || z == Py_None) {
1911     if (MP_NEGP(yy)) VALERR("negative exponent");
1912     r = gf_exp(MP_NEW, xx, yy);
1913   } else {
1914     gfreduce gr;
1915     if (MP_ZEROP(zz)) ZDIVERR("zero modulus");
1916     if (MP_NEGP(yy)) {
1917       if ((xx = gf_modinv_checked(xx, xx, zz)) == 0) goto end;
1918       yy = mp_neg(yy, yy);
1919     }
1920     gfreduce_create(&gr, zz);
1921     r = gfreduce_exp(&gr, MP_NEW, xx, yy);
1922     gfreduce_destroy(&gr);
1923   }
1924   rc = gf_pywrap(r);
1925 end:
1926   mp_drop(xx); mp_drop(yy); mp_drop(zz);
1927   return (rc);
1928 }
1929
1930 static PyObject *gfmeth_sqr(PyObject *me, PyObject *arg)
1931 {
1932   if (!PyArg_ParseTuple(arg, ":sqr")) return (0);
1933   return (gf_pywrap(gf_sqr(MP_NEW, MP_X(me))));
1934 }
1935
1936 static PyObject *gfmeth_gcd(PyObject *me, PyObject *arg)
1937 {
1938   PyObject *z = 0;
1939   mp *yy = 0, *zz = MP_NEW;
1940
1941   if (!PyArg_ParseTuple(arg, "O&:gcd", convgf, &yy)) goto end;
1942   gf_gcd(&zz, 0, 0, MP_X(me), yy);
1943   z = gf_pywrap(zz);
1944 end:
1945   if (yy) MP_DROP(yy);
1946   return (z);
1947 }
1948
1949 static PyObject *gfmeth_gcdx(PyObject *me, PyObject *arg)
1950 {
1951   PyObject *z = 0;
1952   mp *yy = 0, *zz = MP_NEW, *uu = MP_NEW, *vv = MP_NEW;
1953
1954   if (!PyArg_ParseTuple(arg, "O&:gcdx", convgf, &yy))
1955     goto end;
1956   gf_gcd(&zz, &uu, &vv, MP_X(me), yy);
1957   z = Py_BuildValue("(NNN)",
1958                     gf_pywrap(zz), gf_pywrap(uu), gf_pywrap(vv));
1959 end:
1960   if (yy) MP_DROP(yy);
1961   return (z);
1962 }
1963
1964 static PyObject *gfmeth_modinv(PyObject *me, PyObject *arg)
1965 {
1966   PyObject *z = 0;
1967   mp *yy = 0, *zz = MP_NEW;
1968
1969   if (!PyArg_ParseTuple(arg, "O&:modinv", convgf, &yy) ||
1970       (zz = gf_modinv_checked(MP_NEW, yy, MP_X(me))) == 0)
1971     goto end;
1972   z = gf_pywrap(zz);
1973 end:
1974   if (yy) MP_DROP(yy);
1975   return (z);
1976 }
1977
1978 static PyObject *gfmeth_irreduciblep(PyObject *me, PyObject *arg)
1979 {
1980   if (!PyArg_ParseTuple(arg, ":irreduciblep")) return (0);
1981   return getbool(gf_irreduciblep(MP_X(me)));
1982 }
1983
1984 static PyObject *gfget_degree(PyObject *me, void *hunoz)
1985   { return (PyInt_FromLong(mp_bits(MP_X(me)) - 1)); }
1986
1987 static PyGetSetDef gf_pygetset[] = {
1988 #define GETSETNAME(op, name) gf##op##_##name
1989   GET   (degree,        "X.degree -> polynomial degree of X")
1990 #undef GETSETNAME
1991 #define GETSETNAME(op, name) mp##op##_##name
1992   GET   (nbits,         "X.nbits -> bit length of X")
1993   GET   (noctets,       "X.noctets -> octet length of X")
1994 #undef GETSETNAME
1995   { 0 }
1996 };
1997
1998 static PyMethodDef gf_pymethods[] = {
1999 #define METHNAME(func) gfmeth_##func
2000   METH  (setbit,        "X.setbit(N) -> X with bit N set")
2001   METH  (clearbit,      "X.clearbit(N) -> X with bit N clear")
2002   METH  (testbit,       "X.testbit(N) -> true/false if bit N set/clear in X")
2003   METH  (sqr,           "X.sqr() -> X^2")
2004   METH  (gcd,           "X.gcd(Y) -> gcd(X, Y)")
2005   METH  (gcdx,
2006          "X.gcdx(Y) -> (gcd(X, Y), U, V) with X U + Y V = gcd(X, Y)")
2007   METH  (modinv,        "X.modinv(Y) -> multiplicative inverse of Y mod X")
2008   METH  (irreduciblep,  "X.irreduciblep() -> true/false")
2009 #undef METHNAME
2010 #define METHNAME(func) mpmeth_##func
2011   KWMETH(tostring,      "X.tostring(radix = 10) -> STR")
2012   KWMETH(storel,        "X.storel(len = -1) -> little-endian bytes")
2013   KWMETH(storeb,        "X.storeb(len = -1) -> big-endian bytes")
2014   KWMETH(storel2c,
2015          "X.storel2c(len = -1) -> little-endian bytes, two's complement")
2016   KWMETH(storeb2c,
2017          "X.storeb2c(len = -1) -> big-endian bytes, two's complement")
2018   METH  (tobuf,         "X.tobuf() -> buffer format")
2019 #undef METHNAME
2020   { 0 }
2021 };
2022
2023 static PyNumberMethods gf_pynumber = {
2024   gf_pyadd,                             /* @nb_add@ */
2025   gf_pysub,                             /* @nb_subtract@ */
2026   gf_pymul,                             /* @nb_multiply@ */
2027   0,                                    /* @nb_divide@ */
2028   gf_pymod,                             /* @nb_remainder@ */
2029   gf_pydivmod,                          /* @nb_divmod@ */
2030   gf_pyexp,                             /* @nb_power@ */
2031   mp_pyid,                              /* @nb_negative@ */
2032   mp_pyid,                              /* @nb_positive@ */
2033   mp_pyid,                              /* @nb_absolute@ */
2034   mp_pynonzerop,                        /* @nb_nonzero@ */
2035   0 /* doesn't make any sense */,       /* @nb_invert@ */
2036   gf_pylsl,                             /* @nb_lshift@ */
2037   gf_pylsr,                             /* @nb_rshift@ */
2038   gf_pyand,                             /* @nb_and@ */
2039   gf_pyxor,                             /* @nb_xor@ */
2040   gf_pyor,                              /* @nb_or@ */
2041   gf_pycoerce,                          /* @nb_coerce@ */
2042   mp_pyint,                             /* @nb_int@ */
2043   mp_pylong,                            /* @nb_long@ */
2044   0 /* doesn't make any sense */,       /* @nb_float@ */
2045   mp_pyoct,                             /* @nb_oct@ */
2046   mp_pyhex,                             /* @nb_hex@ */
2047
2048   0,                                    /* @nb_inplace_add@ */
2049   0,                                    /* @nb_inplace_subtract@ */
2050   0,                                    /* @nb_inplace_multiply@ */
2051   0,                                    /* @nb_inplace_divide@ */
2052   0,                                    /* @nb_inplace_remainder@ */
2053   0,                                    /* @nb_inplace_power@ */
2054   0,                                    /* @nb_inplace_lshift@ */
2055   0,                                    /* @nb_inplace_rshift@ */
2056   0,                                    /* @nb_inplace_and@ */
2057   0,                                    /* @nb_inplace_xor@ */
2058   0,                                    /* @nb_inplace_or@ */
2059
2060   gf_pydiv,                             /* @nb_floor_divide@ */
2061   0,                                    /* @nb_true_divide@ */
2062   0,                                    /* @nb_inplace_floor_divide@ */
2063   0,                                    /* @nb_inplace_true_divide@ */
2064 };
2065
2066 static PyTypeObject gf_pytype_skel = {
2067   PyObject_HEAD_INIT(0) 0,              /* Header */
2068   "GF",                                 /* @tp_name@ */
2069   sizeof(mp_pyobj),                     /* @tp_basicsize@ */
2070   0,                                    /* @tp_itemsize@ */
2071
2072   mp_pydealloc,                         /* @tp_dealloc@ */
2073   0,                                    /* @tp_print@ */
2074   0,                                    /* @tp_getattr@ */
2075   0,                                    /* @tp_setattr@ */
2076   0,                                    /* @tp_compare@ */
2077   gf_pyrepr,                            /* @tp_repr@ */
2078   &gf_pynumber,                         /* @tp_as_number@ */
2079   0,                                    /* @tp_as_sequence@ */
2080   0,                                    /* @tp_as_mapping@ */
2081   gf_pyhash,                            /* @tp_hash@ */
2082   0,                                    /* @tp_call@ */
2083   mp_pyhex,                             /* @tp_str@ */
2084   0,                                    /* @tp_getattro@ */
2085   0,                                    /* @tp_setattro@ */
2086   0,                                    /* @tp_as_buffer@ */
2087   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
2088     Py_TPFLAGS_CHECKTYPES |
2089     Py_TPFLAGS_BASETYPE,
2090
2091   /* @tp_doc@ */
2092 "Binary polynomials.  Support almost all the standard arithmetic\n\
2093 operations.\n\
2094 \n\
2095 Constructor gf(X, radix = R) attempts to convert X to a `gf'.  If\n\
2096 X is a string, it's read in radix-R form, or we look for a prefix\n\
2097 if R = 0.  Other acceptable things are ints and longs.\n\
2098 \n\
2099 The name is hopelessly wrong from a technical point of view, but\n\
2100 but it's much easier to type than `p2' or `c2' or whatever.\n\
2101 \n\
2102 Notes:\n\
2103 \n\
2104   * Use `//' for division.  GFs don't have `/' division.",
2105
2106   0,                                    /* @tp_traverse@ */
2107   0,                                    /* @tp_clear@ */
2108   gf_pyrichcompare,                     /* @tp_richcompare@ */
2109   0,                                    /* @tp_weaklistoffset@ */
2110   0,                                    /* @tp_iter@ */
2111   0,                                    /* @tp_iternext@ */
2112   gf_pymethods,                         /* @tp_methods@ */
2113   0,                                    /* @tp_members@ */
2114   gf_pygetset,                          /* @tp_getset@ */
2115   0,                                    /* @tp_base@ */
2116   0,                                    /* @tp_dict@ */
2117   0,                                    /* @tp_descr_get@ */
2118   0,                                    /* @tp_descr_set@ */
2119   0,                                    /* @tp_dictoffset@ */
2120   0,                                    /* @tp_init@ */
2121   PyType_GenericAlloc,                  /* @tp_alloc@ */
2122   gf_pynew,                             /* @tp_new@ */
2123   0,                                    /* @tp_free@ */
2124   0                                     /* @tp_is_gc@ */
2125 };
2126
2127 static PyObject *meth__GF_fromstring(PyObject *me,
2128                                     PyObject *arg, PyObject *kw)
2129 {
2130   int r = 0;
2131   char *p;
2132   Py_ssize_t len;
2133   PyObject *z = 0;
2134   mp *zz;
2135   mptext_stringctx sc;
2136   char *kwlist[] = { "class", "x", "radix", 0 };
2137
2138   if (!PyArg_ParseTupleAndKeywords(arg, kw, "Os#|i:fromstring",
2139                                    kwlist, &me, &p, &len, &r))
2140     goto end;
2141   if (!good_radix_p(r, 1)) VALERR("bad radix");
2142   sc.buf = p; sc.lim = p + len;
2143   if ((zz = mp_read(MP_NEW, r, &mptext_stringops, &sc)) == 0 ||
2144       MP_NEGP(zz)) {
2145     if (zz) MP_DROP(zz);
2146     VALERR("bad binary polynomial");
2147   }
2148   z = Py_BuildValue("(Ns#)", gf_pywrap(zz),
2149                     sc.buf, (Py_ssize_t)(sc.lim - sc.buf));
2150 end:
2151   return (z);
2152 }
2153
2154 /*----- Sparse poly reduction ---------------------------------------------*/
2155
2156 typedef struct gfreduce_pyobj {
2157   PyObject_HEAD
2158   gfreduce mr;
2159 } gfreduce_pyobj;
2160
2161 #define GFREDUCE_PY(o) (&((gfreduce_pyobj *)(o))->mr)
2162
2163 static PyObject *grmeth_exp(PyObject *me, PyObject *arg)
2164 {
2165   PyObject *rc = 0;
2166   mp *yy = 0, *zz = 0;
2167
2168   if (!PyArg_ParseTuple(arg, "O&O&:exp", convgf, &yy, convgf, &zz))
2169     goto end;
2170   if (MP_NEGP(zz)) {
2171     if ((yy = gf_modinv_checked(yy, yy, GFREDUCE_PY(me)->p)) == 0) goto end;
2172     zz = mp_neg(zz, zz);
2173   }
2174   rc = gf_pywrap(gfreduce_exp(GFREDUCE_PY(me), MP_NEW, yy, zz));
2175 end:
2176   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
2177   return (rc);
2178 }
2179
2180 static PyObject *grmeth_trace(PyObject *me, PyObject *arg)
2181 {
2182   PyObject *rc = 0;
2183   mp *xx = 0;
2184
2185   if (!PyArg_ParseTuple(arg, "O&:trace", convgf, &xx)) goto end;
2186   rc = PyInt_FromLong(gfreduce_trace(GFREDUCE_PY(me), xx));
2187 end:
2188   if (xx) MP_DROP(xx);
2189   return (rc);
2190 }
2191
2192 static PyObject *grmeth_halftrace(PyObject *me, PyObject *arg)
2193 {
2194   PyObject *rc = 0;
2195   mp *xx = 0;
2196
2197   if (!PyArg_ParseTuple(arg, "O&:halftrace", convgf, &xx)) goto end;
2198   rc = gf_pywrap(gfreduce_halftrace(GFREDUCE_PY(me), MP_NEW, xx));
2199 end:
2200   if (xx) MP_DROP(xx);
2201   return (rc);
2202 }
2203
2204 static PyObject *grmeth_sqrt(PyObject *me, PyObject *arg)
2205 {
2206   PyObject *rc = 0;
2207   mp *xx = 0, *yy;
2208
2209   if (!PyArg_ParseTuple(arg, "O&:sqrt", convgf, &xx)) goto end;
2210   if ((yy = gfreduce_sqrt(GFREDUCE_PY(me), MP_NEW, xx)) == 0)
2211     VALERR("no modular square root");
2212   rc = gf_pywrap(yy);
2213 end:
2214   if (xx) MP_DROP(xx);
2215   return (rc);
2216 }
2217
2218 static PyObject *grmeth_quadsolve(PyObject *me, PyObject *arg)
2219 {
2220   PyObject *rc = 0;
2221   mp *xx = 0, *yy;
2222
2223   if (!PyArg_ParseTuple(arg, "O&:quadsolve", convgf, &xx)) goto end;
2224   if ((yy = gfreduce_quadsolve(GFREDUCE_PY(me), MP_NEW, xx)) == 0)
2225     VALERR("no solution found");
2226   rc = gf_pywrap(yy);
2227 end:
2228   if (xx) MP_DROP(xx);
2229   return (rc);
2230 }
2231
2232 static PyObject *grmeth_reduce(PyObject *me, PyObject *arg)
2233 {
2234   PyObject *z = 0;
2235   mp *yy = 0;
2236
2237   if (!PyArg_ParseTuple(arg, "O&:reduce", convgf, &yy)) goto end;
2238   z = gf_pywrap(gfreduce_do(GFREDUCE_PY(me), MP_NEW, yy));
2239 end:
2240   return (z);
2241 }
2242
2243 static void gfreduce_pydealloc(PyObject *me)
2244 {
2245   gfreduce_destroy(GFREDUCE_PY(me));
2246   FREEOBJ(me);
2247 }
2248
2249 static PyObject *gfreduce_pynew(PyTypeObject *ty,
2250                                  PyObject *arg, PyObject *kw)
2251 {
2252   gfreduce_pyobj *mr = 0;
2253   gfreduce r;
2254   char *kwlist[] = { "m", 0 };
2255   mp *xx = 0;
2256
2257   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convgf, &xx))
2258     goto end;
2259   if (MP_ZEROP(xx)) ZDIVERR("modulus is zero!");
2260   gfreduce_create(&r, xx);
2261   mr = (gfreduce_pyobj *)ty->tp_alloc(ty, 0);
2262   mr->mr = r;
2263 end:
2264   if (xx) MP_DROP(xx);
2265   return ((PyObject *)mr);
2266 }
2267
2268 static PyObject *grget_m(PyObject *me, void *hunoz)
2269   { return (gf_pywrap(MP_COPY(GFREDUCE_PY(me)->p))); }
2270
2271 static PyGetSetDef gfreduce_pygetset[] = {
2272 #define GETSETNAME(op, name) gr##op##_##name
2273   GET   (m,             "R.m -> reduction polynomial")
2274 #undef GETSETNAME
2275   { 0 }
2276 };
2277
2278 static PyMethodDef gfreduce_pymethods[] = {
2279 #define METHNAME(name) grmeth_##name
2280   METH  (reduce,        "R.reduce(X) -> X mod B.m")
2281   METH  (trace,        "R.trace(X) -> Tr(X) = x + x^2 + ... + x^{2^{m - 1}}")
2282   METH  (halftrace,   "R.halftrace(X) -> x + x^{2^2} + ... + x^{2^{m - 1}}")
2283   METH  (sqrt,          "R.sqrt(X) -> Y where Y^2 = X mod R")
2284   METH  (quadsolve,     "R.quadsolve(X) -> Y where Y^2 + Y = X mod R")
2285   METH  (exp,           "R.exp(X, N) -> X^N mod B.m")
2286 #undef METHNAME
2287   { 0 }
2288 };
2289
2290 static PyTypeObject *gfreduce_pytype, gfreduce_pytype_skel = {
2291   PyObject_HEAD_INIT(0) 0,              /* Header */
2292   "GFReduce",                           /* @tp_name@ */
2293   sizeof(gfreduce_pyobj),               /* @tp_basicsize@ */
2294   0,                                    /* @tp_itemsize@ */
2295
2296   gfreduce_pydealloc,                   /* @tp_dealloc@ */
2297   0,                                    /* @tp_print@ */
2298   0,                                    /* @tp_getattr@ */
2299   0,                                    /* @tp_setattr@ */
2300   0,                                    /* @tp_compare@ */
2301   0,                                    /* @tp_repr@ */
2302   0,                                    /* @tp_as_number@ */
2303   0,                                    /* @tp_as_sequence@ */
2304   0,                                    /* @tp_as_mapping@ */
2305   0,                                    /* @tp_hash@ */
2306   0,                                    /* @tp_call@ */
2307   0,                                    /* @tp_str@ */
2308   0,                                    /* @tp_getattro@ */
2309   0,                                    /* @tp_setattro@ */
2310   0,                                    /* @tp_as_buffer@ */
2311   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
2312     Py_TPFLAGS_BASETYPE,
2313
2314   /* @tp_doc@ */
2315 "A reduction context for reduction modulo sparse irreducible polynomials.",
2316
2317   0,                                    /* @tp_traverse@ */
2318   0,                                    /* @tp_clear@ */
2319   0,                                    /* @tp_richcompare@ */
2320   0,                                    /* @tp_weaklistoffset@ */
2321   0,                                    /* @tp_iter@ */
2322   0,                                    /* @tp_iternext@ */
2323   gfreduce_pymethods,                   /* @tp_methods@ */
2324   0,                                    /* @tp_members@ */
2325   gfreduce_pygetset,                    /* @tp_getset@ */
2326   0,                                    /* @tp_base@ */
2327   0,                                    /* @tp_dict@ */
2328   0,                                    /* @tp_descr_get@ */
2329   0,                                    /* @tp_descr_set@ */
2330   0,                                    /* @tp_dictoffset@ */
2331   0,                                    /* @tp_init@ */
2332   PyType_GenericAlloc,                  /* @tp_alloc@ */
2333   gfreduce_pynew,                       /* @tp_new@ */
2334   0,                                    /* @tp_free@ */
2335   0                                     /* @tp_is_gc@ */
2336 };
2337
2338 /*----- Normal/poly transformation ----------------------------------------*/
2339
2340 typedef struct gfn_pyobj {
2341   PyObject_HEAD
2342   mp *p;
2343   gfn ntop, pton;
2344 } gfn_pyobj;
2345
2346 static PyTypeObject *gfn_pytype, gfn_pytype_skel;
2347
2348 #define GFN_P(o) (((gfn_pyobj *)(o))->p)
2349 #define GFN_PTON(o) (&((gfn_pyobj *)(o))->pton)
2350 #define GFN_NTOP(o) (&((gfn_pyobj *)(o))->ntop)
2351
2352 static PyObject *gfn_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
2353 {
2354   mp *p = 0, *beta = 0;
2355   gfn_pyobj *gg = 0;
2356   char *kwlist[] = { "p", "beta", 0 };
2357
2358   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&O&:new", kwlist,
2359                                    convgf, &p, convgf, &beta))
2360     goto end;
2361   gg = PyObject_New(gfn_pyobj, ty);
2362   if (gfn_create(p, beta, &gg->ntop, &gg->pton)) {
2363     FREEOBJ(gg);
2364     gg = 0;
2365     VALERR("can't invert transformation matrix");
2366   }
2367   gg->p = MP_COPY(p);
2368 end:
2369   mp_drop(p);
2370   mp_drop(beta);
2371   return ((PyObject *)gg);
2372 }
2373
2374 static PyObject *gfnget_p(PyObject *me, void *hunoz)
2375   { return (gf_pywrap(MP_COPY(GFN_P(me)))); }
2376
2377 static PyObject *gfnget_beta(PyObject *me, void *hunoz)
2378 {
2379   gfn *n = GFN_NTOP(me);
2380   mp *x = n->r[n->n - 1];
2381   return (gf_pywrap(MP_COPY(x)));
2382 }
2383
2384 #define XFORMOP(name, NAME)                                             \
2385   static PyObject *gfnmeth_##name(PyObject *me, PyObject *arg)          \
2386   {                                                                     \
2387     mp *xx = 0;                                                         \
2388     mp *z = 0;                                                          \
2389                                                                         \
2390     if (!PyArg_ParseTuple(arg, "O&:" #name, convgf, &xx)) goto end;     \
2391     z = gfn_transform(GFN_##NAME(me), MP_NEW, xx);                      \
2392   end:                                                                  \
2393     mp_drop(xx);                                                        \
2394     if (!z) return (0);                                                 \
2395     return (mp_pywrap(z));                                              \
2396   }
2397 XFORMOP(pton, PTON)
2398 XFORMOP(ntop, NTOP)
2399 #undef XFORMOP
2400
2401 static void gfn_pydealloc(PyObject *me)
2402 {
2403   gfn_destroy(GFN_PTON(me));
2404   gfn_destroy(GFN_NTOP(me));
2405   FREEOBJ(me);
2406 }
2407
2408 static PyGetSetDef gfn_pygetset[] = {
2409 #define GETSETNAME(op, name) gfn##op##_##name
2410   GET   (p,             "X.p -> polynomial basis, as polynomial")
2411   GET   (beta,          "X.beta -> normal basis element, in poly form")
2412 #undef GETSETNAME
2413   { 0 }
2414 };
2415
2416 static PyMethodDef gfn_pymethods[] = {
2417 #define METHNAME(name) gfnmeth_##name
2418   METH  (pton,          "X.pton(A) -> normal-basis representation of A")
2419   METH  (ntop,          "X.ntop(A) -> polynomial-basis representation of A")
2420 #undef METHNAME
2421   { 0 }
2422 };
2423
2424 static PyTypeObject gfn_pytype_skel = {
2425   PyObject_HEAD_INIT(0) 0,              /* Header */
2426   "GFN",                                /* @tp_name@ */
2427   sizeof(gfn_pyobj),                    /* @tp_basicsize@ */
2428   0,                                    /* @tp_itemsize@ */
2429
2430   gfn_pydealloc,                        /* @tp_dealloc@ */
2431   0,                                    /* @tp_print@ */
2432   0,                                    /* @tp_getattr@ */
2433   0,                                    /* @tp_setattr@ */
2434   0,                                    /* @tp_compare@ */
2435   0,                                    /* @tp_repr@ */
2436   0,                                    /* @tp_as_number@ */
2437   0,                                    /* @tp_as_sequence@ */
2438   0,                                    /* @tp_as_mapping@ */
2439   0,                                    /* @tp_hash@ */
2440   0,                                    /* @tp_call@ */
2441   0,                                    /* @tp_str@ */
2442   0,                                    /* @tp_getattro@ */
2443   0,                                    /* @tp_setattro@ */
2444   0,                                    /* @tp_as_buffer@ */
2445   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
2446     Py_TPFLAGS_BASETYPE,
2447
2448   /* @tp_doc@ */
2449 "An object for transforming elements of binary fields between polynomial\n\
2450 and normal basis representations.",
2451
2452   0,                                    /* @tp_traverse@ */
2453   0,                                    /* @tp_clear@ */
2454   0,                                    /* @tp_richcompare@ */
2455   0,                                    /* @tp_weaklistoffset@ */
2456   0,                                    /* @tp_iter@ */
2457   0,                                    /* @tp_iternext@ */
2458   gfn_pymethods,                        /* @tp_methods@ */
2459   0,                                    /* @tp_members@ */
2460   gfn_pygetset,                         /* @tp_getset@ */
2461   0,                                    /* @tp_base@ */
2462   0,                                    /* @tp_dict@ */
2463   0,                                    /* @tp_descr_get@ */
2464   0,                                    /* @tp_descr_set@ */
2465   0,                                    /* @tp_dictoffset@ */
2466   0,                                    /* @tp_init@ */
2467   PyType_GenericAlloc,                  /* @tp_alloc@ */
2468   gfn_pynew,                            /* @tp_new@ */
2469   0,                                    /* @tp_free@ */
2470   0                                     /* @tp_is_gc@ */
2471 };
2472
2473 /*----- Glue --------------------------------------------------------------*/
2474
2475 static PyMethodDef methods[] = {
2476 #define METHNAME(func) meth_##func
2477   KWMETH(_MP_fromstring,        "\
2478 fromstring(STR, radix = 0) -> (X, REST)\n\
2479 \n\
2480 Parse STR as a large integer, according to radix.  If radix is zero,\n\
2481 read a prefix from STR to decide radix: allow `0' for octal, `0x' for hex\n\
2482 or `R_' for other radix R.")
2483   KWMETH(_GF_fromstring,        "\
2484 fromstring(STR, radix = 0) -> (X, REST)\n\
2485 \n\
2486 Parse STR as a binary polynomial, according to radix.  If radix is zero,\n\
2487 read a prefix from STR to decide radix: allow `0' for octal, `0x' for hex\n\
2488 or `R_' for other radix R.")
2489   METH  (_MP_factorial,         "\
2490 factorial(I) -> I!: compute factorial")
2491   METH  (_MP_fibonacci,         "\
2492 fibonacci(I) -> F(I): compute Fibonacci number")
2493   METH  (_MP_loadl,             "\
2494 loadl(STR) -> X: read little-endian bytes")
2495   METH  (_MP_loadb,             "\
2496 loadb(STR) -> X: read big-endian bytes")
2497   METH  (_MP_loadl2c,           "\
2498 loadl2c(STR) -> X: read little-endian bytes, two's complement")
2499   METH  (_MP_loadb2c,           "\
2500 loadb2c(STR) -> X: read big-endian bytes, two's complement")
2501   METH  (_MP_frombuf,           "\
2502 frombuf(STR) -> (X, REST): read buffer format")
2503   METH  (_GF_loadl,             "\
2504 loadl(STR) -> X: read little-endian bytes")
2505   METH  (_GF_loadb,             "\
2506 loadb(STR) -> X: read big-endian bytes")
2507   METH  (_GF_frombuf,           "\
2508 frombuf(STR) -> (X, REST): read buffer format")
2509 #undef METHNAME
2510   { 0 }
2511 };
2512
2513 void mp_pyinit(void)
2514 {
2515   INITTYPE(mp, root);
2516   INITTYPE(gf, root);
2517   INITTYPE(mpmul, root);
2518   INITTYPE(mpmont, root);
2519   INITTYPE(mpbarrett, root);
2520   INITTYPE(mpreduce, root);
2521   INITTYPE(mpcrt, root);
2522   INITTYPE(gfreduce, root);
2523   INITTYPE(gfn, root);
2524   addmethods(methods);
2525 }
2526
2527 void mp_pyinsert(PyObject *mod)
2528 {
2529   INSERT("MP", mp_pytype);
2530   INSERT("MPMul", mpmul_pytype);
2531   INSERT("MPMont", mpmont_pytype);
2532   INSERT("MPBarrett", mpbarrett_pytype);
2533   INSERT("MPReduce", mpreduce_pytype);
2534   INSERT("MPCRT", mpcrt_pytype);
2535   INSERT("GF", gf_pytype);
2536   INSERT("GFReduce", gfreduce_pytype);
2537   INSERT("GFN", gfn_pytype);
2538 }
2539
2540 /*----- That's all, folks -------------------------------------------------*/