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

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