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

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