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 static PyObject *meth__MP_factorial(PyObject *me, PyObject *arg)
 931 {
 932   unsigned long i;
 933   mp *x;
 934   if (!PyArg_ParseTuple(arg, "OO&:factorial", &me, convulong, &i))
 935     return (0);
 936   x = mp_factorial(i);
 937   return mp_pywrap(x);
 938 }
 939
 940 #define LOADOP(pre, py, name)                                           \
 941   static PyObject *meth__##py##_##name(PyObject *me, PyObject *arg)     \
 942   {                                                                     \
 943     char *p;                                                            \
 944     int len;                                                            \
 945     if (!PyArg_ParseTuple(arg, "Os#:" #name, &me, &p, &len)) return (0); \
 946     return (pre##_pywrap(mp_##name(MP_NEW, p, len)));                   \
 947   }
 948 LOADOP(mp, MP, loadl)
 949 LOADOP(mp, MP, loadb)
 950 LOADOP(mp, MP, loadl2c)
 951 LOADOP(mp, MP, loadb2c)
 952 LOADOP(gf, GF, loadl)
 953 LOADOP(gf, GF, loadb)
 954 #undef LOADOP
 955
 956 /*----- Products of small integers ----------------------------------------*/
 957
 958 typedef struct mpmul_pyobj {
 959   PyObject_HEAD
 960   int livep;
 961   mpmul mm;
 962 } mpmul_pyobj;
 963
 964 #define MPMUL_LIVEP(o) (((mpmul_pyobj *)(o))->livep)
 965 #define MPMUL_PY(o) (&((mpmul_pyobj *)(o))->mm)
 966
 967 static void mpmul_pydealloc(PyObject *me)
 968 {
 969   if (MPMUL_LIVEP(me))
 970     mp_drop(mpmul_done(MPMUL_PY(me)));
 971   FREEOBJ(me);
 972 }
 973
 974 static PyObject *mmmeth_factor(PyObject *me, PyObject *arg)
 975 {
 976   PyObject *q, *i;
 977   mp *x;
 978
 979   if (!MPMUL_LIVEP(me)) VALERR("MPMul object invalid");
 980   if (PyTuple_Size(arg) != 1)
 981     i = PyObject_GetIter(arg);
 982   else {
 983     if ((q = PyTuple_GetItem(arg, 0)) == 0) goto end;
 984     if ((i = PyObject_GetIter(q)) == 0) {
 985       PyErr_Clear(); /* that's ok */
 986       i = PyObject_GetIter(arg);
 987     }
 988   }
 989   if (!i) goto end;
 990   while ((q = PyIter_Next(i)) != 0) {
 991     x = getmp(q); Py_DECREF(q); if (!x) {
 992       Py_DECREF(i);
 993       goto end;
 994     }
 995     mpmul_add(MPMUL_PY(me), x);
 996     MP_DROP(x);
 997   }
 998   Py_DECREF(i);
 999   RETURN_ME;
1000 end:
1001   return (0);
1002 }
1003
1004 static PyObject *mmmeth_done(PyObject *me, PyObject *arg)
1005 {
1006   mp *x;
1007
1008   if (!PyArg_ParseTuple(arg, ":done")) goto end;
1009   if (!MPMUL_LIVEP(me)) VALERR("MPMul object invalid");
1010   x = mpmul_done(MPMUL_PY(me));
1011   MPMUL_LIVEP(me) = 0;
1012   return (mp_pywrap(x));
1013 end:
1014   return (0);
1015 }
1016
1017 static PyObject *mpmul_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1018 {
1019   mpmul_pyobj *mm;
1020
1021   if (kw) TYERR("keyword arguments not allowed here");
1022   mm = (mpmul_pyobj *)ty->tp_alloc(ty, 0);
1023   mpmul_init(&mm->mm);
1024   mm->livep = 1;
1025   if (mmmeth_factor((PyObject *)mm, arg) == 0) {
1026     Py_DECREF(mm);
1027     goto end;
1028   }
1029   return ((PyObject *)mm);
1030 end:
1031   return (0);
1032 }
1033
1034 static PyObject *mmget_livep(PyObject *me, void *hunoz)
1035   { return (getbool(MPMUL_LIVEP(me))); }
1036
1037 static PyGetSetDef mpmul_pygetset[] = {
1038 #define GETSETNAME(op, name) mm##op##_##name
1039   GET   (livep,                 "MM.livep -> flag: object still valid?")
1040 #undef GETSETNAME
1041   { 0 }
1042 };
1043
1044 static PyMethodDef mpmul_pymethods[] = {
1045 #define METHNAME(name) mmmeth_##name
1046   METH  (factor,                "MM.factor(ITERABLE) or MM.factor(I, ...)")
1047   METH  (done,                  "MM.done() -> PRODUCT")
1048 #undef METHNAME
1049   { 0 }
1050 };
1051
1052 static PyTypeObject *mpmul_pytype, mpmul_pytype_skel = {
1053   PyObject_HEAD_INIT(0) 0,              /* Header */
1054   "catacomb.MPMul",                     /* @tp_name@ */
1055   sizeof(mpmul_pyobj),                  /* @tp_basicsize@ */
1056   0,                                    /* @tp_itemsize@ */
1057
1058   mpmul_pydealloc,                      /* @tp_dealloc@ */
1059   0,                                    /* @tp_print@ */
1060   0,                                    /* @tp_getattr@ */
1061   0,                                    /* @tp_setattr@ */
1062   0,                                    /* @tp_compare@ */
1063   0,                                    /* @tp_repr@ */
1064   0,                                    /* @tp_as_number@ */
1065   0,                                    /* @tp_as_sequence@ */
1066   0,                                    /* @tp_as_mapping@ */
1067   0,                                    /* @tp_hash@ */
1068   0,                                    /* @tp_call@ */
1069   0,                                    /* @tp_str@ */
1070   0,                                    /* @tp_getattro@ */
1071   0,                                    /* @tp_setattro@ */
1072   0,                                    /* @tp_as_buffer@ */
1073   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1074     Py_TPFLAGS_BASETYPE,
1075
1076   /* @tp_doc@ */
1077 "An object for multiplying many small integers.",
1078
1079   0,                                    /* @tp_traverse@ */
1080   0,                                    /* @tp_clear@ */
1081   0,                                    /* @tp_richcompare@ */
1082   0,                                    /* @tp_weaklistoffset@ */
1083   0,                                    /* @tp_iter@ */
1084   0,                                    /* @tp_iternext@ */
1085   mpmul_pymethods,                      /* @tp_methods@ */
1086   0,                                    /* @tp_members@ */
1087   mpmul_pygetset,                       /* @tp_getset@ */
1088   0,                                    /* @tp_base@ */
1089   0,                                    /* @tp_dict@ */
1090   0,                                    /* @tp_descr_get@ */
1091   0,                                    /* @tp_descr_set@ */
1092   0,                                    /* @tp_dictoffset@ */
1093   0,                                    /* @tp_init@ */
1094   PyType_GenericAlloc,                  /* @tp_alloc@ */
1095   mpmul_pynew,                          /* @tp_new@ */
1096   0,                                    /* @tp_free@ */
1097   0                                     /* @tp_is_gc@ */
1098 };
1099
1100 /*----- Montgomery reduction ----------------------------------------------*/
1101
1102 typedef struct mpmont_pyobj {
1103   PyObject_HEAD
1104   mpmont mm;
1105 } mpmont_pyobj;
1106
1107 #define MPMONT_PY(o) (&((mpmont_pyobj *)(o))->mm)
1108
1109 static PyObject *mmmeth_int(PyObject *me, PyObject *arg)
1110 {
1111   PyObject *z = 0;
1112   mp *yy = 0;
1113   mpmont *mm = MPMONT_PY(me);
1114
1115   if (!PyArg_ParseTuple(arg, "O&:in", convmp, &yy))
1116     goto end;
1117   mp_div(0, &yy, yy, mm->m);
1118   z = mp_pywrap(mpmont_mul(mm, MP_NEW, yy, mm->r2));
1119 end:
1120   if (yy) MP_DROP(yy);
1121   return (z);
1122 }
1123
1124 static PyObject *mmmeth_mul(PyObject *me, PyObject *arg)
1125 {
1126   PyObject *rc = 0;
1127   mp *yy = 0, *zz = 0;
1128
1129   if (!PyArg_ParseTuple(arg, "O&O&:mul", convmp, &yy, convmp, &zz))
1130     goto end;
1131   rc = mp_pywrap(mpmont_mul(MPMONT_PY(me), MP_NEW, yy, zz));
1132 end:
1133   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1134   return (rc);
1135 }
1136
1137 static PyObject *mmmeth_exp(PyObject *me, PyObject *arg)
1138 {
1139   PyObject *rc = 0;
1140   mp *yy = 0, *zz = 0;
1141
1142   if (!PyArg_ParseTuple(arg, "O&O&:exp", convmp, &yy, convmp, &zz))
1143     goto end;
1144   if (MP_NEGP(zz)) {
1145     if ((yy = mp_modinv_checked(yy, yy, MPMONT_PY(me)->m)) == 0) goto end;
1146     zz = mp_neg(zz, zz);
1147   }
1148   rc = mp_pywrap(mpmont_exp(MPMONT_PY(me), MP_NEW, yy, zz));
1149 end:
1150   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1151   return (rc);
1152 }
1153
1154 static PyObject *mmmeth_expr(PyObject *me, PyObject *arg)
1155 {
1156   PyObject *rc = 0;
1157   mp *yy = 0, *zz = 0;
1158
1159   if (!PyArg_ParseTuple(arg, "O&O&:expr", convmp, &yy, convmp, &zz))
1160     goto end;
1161   if (MP_NEGP(zz)) {
1162     yy = mpmont_reduce(MPMONT_PY(me), yy, yy);
1163     if ((yy = mp_modinv_checked(yy, yy, MPMONT_PY(me)->m)) == 0) goto end;
1164     yy = mpmont_mul(MPMONT_PY(me), yy, yy, MPMONT_PY(me)->r2);
1165     zz = mp_neg(zz, zz);
1166   }
1167   rc = mp_pywrap(mpmont_expr(MPMONT_PY(me), MP_NEW, yy, zz));
1168 end:
1169   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1170   return (rc);
1171 }
1172
1173 static PyObject *mm_mexpr_id(PyObject *me)
1174   { return mp_pywrap(MP_COPY(MPMONT_PY(me)->r)); }
1175
1176 static int mm_mexpr_fill(void *p, PyObject *me, PyObject *x, PyObject *y)
1177 {
1178   mp *xx = 0, *yy = 0;
1179   mp_expfactor *f = p;
1180   mpmont *mm = MPMONT_PY(me);
1181
1182   if ((xx = getmp(x)) == 0 || (yy = getmp(y)) == 0)
1183     goto fail;
1184   if (MP_NEGP(yy)) {
1185     xx = mpmont_reduce(mm, xx, xx);
1186     if ((xx = mp_modinv_checked(xx, xx, yy)) == 0)
1187       goto fail;
1188     xx = mpmont_mul(mm, xx, xx, mm->r2);
1189     yy = mp_neg(yy, yy);
1190   }
1191   f->base = xx;
1192   f->exp = yy;
1193   return (0);
1194
1195 fail:
1196   mp_drop(xx); mp_drop(yy);
1197   return (-1);
1198 }
1199
1200 static PyObject *mm_mexpr(PyObject *me, void *v, int n)
1201   { return mp_pywrap(mpmont_mexpr(MPMONT_PY(me), MP_NEW, v, n)); }
1202
1203 static void mp_mexp_drop(void *p)
1204 {
1205   mp_expfactor *f = p;
1206   mp_drop(f->base);
1207   mp_drop(f->exp);
1208 }
1209
1210 static PyObject *mmmeth_mexpr(PyObject *me, PyObject *arg)
1211 {
1212   return mexp_common(me, arg, sizeof(mp_expfactor),
1213                      mm_mexpr_id, mm_mexpr_fill, mm_mexpr, mp_mexp_drop);
1214 }
1215
1216 static PyObject *mp_mexp_id(PyObject *me)
1217   { return mp_pywrap(MP_ONE); }
1218
1219 static int mp_mexp_fill(void *p, PyObject *me, PyObject *x, PyObject *y)
1220 {
1221   mp *xx = 0, *yy = 0;
1222   mp_expfactor *f = p;
1223
1224   if ((xx = getmp(x)) == 0 || (yy = getmp(y)) == 0)
1225     goto fail;
1226   if (MP_NEGP(yy)) {
1227     if ((xx = mp_modinv_checked(xx, xx, yy)) == 0)
1228       goto fail;
1229     yy = mp_neg(yy, yy);
1230   }
1231   f->base = xx;
1232   f->exp = yy;
1233   return (0);
1234
1235 fail:
1236   mp_drop(xx); mp_drop(yy);
1237   return (-1);
1238 }
1239
1240 static PyObject *mm_mexp(PyObject *me, void *v, int n)
1241   { return mp_pywrap(mpmont_mexp(MPMONT_PY(me), MP_NEW, v, n)); }
1242
1243 static PyObject *mmmeth_mexp(PyObject *me, PyObject *arg)
1244 {
1245   return mexp_common(me, arg, sizeof(mp_expfactor),
1246                      mp_mexp_id, mp_mexp_fill, mm_mexp, mp_mexp_drop);
1247 }
1248
1249 #define mmmeth_ext mmmeth_reduce
1250 static PyObject *mmmeth_reduce(PyObject *me, PyObject *arg)
1251 {
1252   PyObject *z = 0;
1253   mp *yy = 0;
1254
1255   if (!PyArg_ParseTuple(arg, "O&", convmp, &yy)) goto end;
1256   z = mp_pywrap(mpmont_reduce(MPMONT_PY(me), MP_NEW, yy));
1257 end:
1258   return (z);
1259 }
1260
1261 static void mpmont_pydealloc(PyObject *me)
1262 {
1263   mpmont_destroy(MPMONT_PY(me));
1264   FREEOBJ(me);
1265 }
1266
1267 static PyObject *mpmont_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1268 {
1269   mpmont_pyobj *mm = 0;
1270   char *kwlist[] = { "m", 0 };
1271   mp *xx = 0;
1272
1273   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convmp, &xx))
1274     goto end;
1275   if (!MP_POSP(xx) || !MP_ODDP(xx)) VALERR("m must be positive and odd");
1276   mm = (mpmont_pyobj *)ty->tp_alloc(ty, 0);
1277   mpmont_create(&mm->mm, xx);
1278 end:
1279   if (xx) MP_DROP(xx);
1280   return ((PyObject *)mm);
1281 }
1282
1283 static PyObject *mmget_m(PyObject *me, void *hunoz)
1284   { return (mp_pywrap(MP_COPY(MPMONT_PY(me)->m))); }
1285
1286 static PyObject *mmget_r(PyObject *me, void *hunoz)
1287   { return (mp_pywrap(MP_COPY(MPMONT_PY(me)->r))); }
1288
1289 static PyObject *mmget_r2(PyObject *me, void *hunoz)
1290   { return (mp_pywrap(MP_COPY(MPMONT_PY(me)->r2))); }
1291
1292 static PyGetSetDef mpmont_pygetset[] = {
1293 #define GETSETNAME(op, name) mm##op##_##name
1294   GET   (m,             "M.m -> modulus for reduction")
1295   GET   (r,             "M.r -> multiplicative identity")
1296   GET   (r2,            "M.r2 -> M.r^2, Montgomerization factor")
1297 #undef GETSETNAME
1298   { 0 }
1299 };
1300
1301 static PyMethodDef mpmont_pymethods[] = {
1302 #define METHNAME(name) mmmeth_##name
1303   METH  (int,           "M.out(X) -> XR")
1304   METH  (mul,           "M.mul(XR, YR) -> ZR where Z = X Y")
1305   METH  (expr,          "M.expr(XR, N) -> ZR where Z = X^N mod M.m")
1306   METH  (mexpr,         "\
1307 B.mexp([(XR0, N0), (XR1, N1), ...]) = ZR where Z = X0^N0 X1^N1 mod B.m\n\
1308 \t(the list may be flattened if this more convenient.)")
1309   METH  (reduce,        "M.reduce(XR) -> X")
1310   METH  (ext,           "M.ext(XR) -> X")
1311   METH  (exp,           "M.exp(X, N) -> X^N mod M.m")
1312   METH  (mexp,          "\
1313 B.mexp([(X0, N0), (X1, N1), ...]) = X0^N0 X1^N1 mod B.m\n\
1314 \t(the list may be flattened if this more convenient.)")
1315 #undef METHNAME
1316   { 0 }
1317 };
1318
1319 static PyTypeObject *mpmont_pytype, mpmont_pytype_skel = {
1320   PyObject_HEAD_INIT(0) 0,              /* Header */
1321   "catacomb.MPMont",                    /* @tp_name@ */
1322   sizeof(mpmont_pyobj),                 /* @tp_basicsize@ */
1323   0,                                    /* @tp_itemsize@ */
1324
1325   mpmont_pydealloc,                     /* @tp_dealloc@ */
1326   0,                                    /* @tp_print@ */
1327   0,                                    /* @tp_getattr@ */
1328   0,                                    /* @tp_setattr@ */
1329   0,                                    /* @tp_compare@ */
1330   0,                                    /* @tp_repr@ */
1331   0,                                    /* @tp_as_number@ */
1332   0,                                    /* @tp_as_sequence@ */
1333   0,                                    /* @tp_as_mapping@ */
1334   0,                                    /* @tp_hash@ */
1335   0,                                    /* @tp_call@ */
1336   0,                                    /* @tp_str@ */
1337   0,                                    /* @tp_getattro@ */
1338   0,                                    /* @tp_setattro@ */
1339   0,                                    /* @tp_as_buffer@ */
1340   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1341     Py_TPFLAGS_BASETYPE,
1342
1343   /* @tp_doc@ */
1344 "A Montgomery reduction context.",
1345
1346   0,                                    /* @tp_traverse@ */
1347   0,                                    /* @tp_clear@ */
1348   0,                                    /* @tp_richcompare@ */
1349   0,                                    /* @tp_weaklistoffset@ */
1350   0,                                    /* @tp_iter@ */
1351   0,                                    /* @tp_iternext@ */
1352   mpmont_pymethods,                     /* @tp_methods@ */
1353   0,                                    /* @tp_members@ */
1354   mpmont_pygetset,                      /* @tp_getset@ */
1355   0,                                    /* @tp_base@ */
1356   0,                                    /* @tp_dict@ */
1357   0,                                    /* @tp_descr_get@ */
1358   0,                                    /* @tp_descr_set@ */
1359   0,                                    /* @tp_dictoffset@ */
1360   0,                                    /* @tp_init@ */
1361   PyType_GenericAlloc,                  /* @tp_alloc@ */
1362   mpmont_pynew,                         /* @tp_new@ */
1363   0,                                    /* @tp_free@ */
1364   0                                     /* @tp_is_gc@ */
1365 };
1366
1367 /*----- Barrett reduction -------------------------------------------------*/
1368
1369 typedef struct mpbarrett_pyobj {
1370   PyObject_HEAD
1371   mpbarrett mb;
1372 } mpbarrett_pyobj;
1373
1374 #define MPBARRETT_PY(o) (&((mpbarrett_pyobj *)(o))->mb)
1375
1376 static PyObject *mbmeth_exp(PyObject *me, PyObject *arg)
1377 {
1378   PyObject *rc = 0;
1379   mp *yy = 0, *zz = 0;
1380
1381   if (!PyArg_ParseTuple(arg, "O&O&:exp", convmp, &yy, convmp, &zz))
1382     goto end;
1383   if (MP_NEGP(zz)) {
1384     if ((yy = mp_modinv_checked(yy, yy, MPBARRETT_PY(me)->m)) == 0) goto end;
1385     zz = mp_neg(zz, zz);
1386   }
1387   rc = mp_pywrap(mpbarrett_exp(MPBARRETT_PY(me), MP_NEW, yy, zz));
1388 end:
1389   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1390   return (rc);
1391 }
1392
1393 static PyObject *mb_mexp(PyObject *me, void *v, int n)
1394   { return mp_pywrap(mpbarrett_mexp(MPBARRETT_PY(me), MP_NEW, v, n)); }
1395
1396 static PyObject *mbmeth_mexp(PyObject *me, PyObject *arg)
1397 {
1398   return mexp_common(me, arg, sizeof(mp_expfactor),
1399                      mp_mexp_id, mp_mexp_fill, mb_mexp, mp_mexp_drop);
1400 }
1401
1402 static PyObject *mbmeth_reduce(PyObject *me, PyObject *arg)
1403 {
1404   PyObject *z = 0;
1405   mp *yy = 0;
1406
1407   if (!PyArg_ParseTuple(arg, "O&:reduce", convmp, &yy))
1408     goto end;
1409   z = mp_pywrap(mpbarrett_reduce(MPBARRETT_PY(me), MP_NEW, yy));
1410 end:
1411   return (z);
1412 }
1413
1414 static void mpbarrett_pydealloc(PyObject *me)
1415 {
1416   mpbarrett_destroy(MPBARRETT_PY(me));
1417   FREEOBJ(me);
1418 }
1419
1420 static PyObject *mpbarrett_pynew(PyTypeObject *ty,
1421                                  PyObject *arg, PyObject *kw)
1422 {
1423   mpbarrett_pyobj *mb = 0;
1424   char *kwlist[] = { "m", 0 };
1425   mp *xx = 0;
1426
1427   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convmp, &xx))
1428     goto end;
1429   if (!MP_POSP(xx)) VALERR("m must be positive");
1430   mb = (mpbarrett_pyobj *)ty->tp_alloc(ty, 0);
1431   mpbarrett_create(&mb->mb, xx);
1432 end:
1433   if (xx) MP_DROP(xx);
1434   return ((PyObject *)mb);
1435 }
1436
1437 static PyObject *mbget_m(PyObject *me, void *hunoz)
1438   { return (mp_pywrap(MP_COPY(MPBARRETT_PY(me)->m))); }
1439
1440 static PyGetSetDef mpbarrett_pygetset[] = {
1441 #define GETSETNAME(op, name) mb##op##_##name
1442   GET   (m,             "B.m -> modulus for reduction")
1443 #undef GETSETNAME
1444   { 0 }
1445 };
1446
1447 static PyMethodDef mpbarrett_pymethods[] = {
1448 #define METHNAME(name) mbmeth_##name
1449   METH  (reduce,        "B.reduce(X) -> X mod B.m")
1450   METH  (exp,           "B.exp(X, N) -> X^N mod B.m")
1451   METH  (mexp,          "\
1452 B.mexp([(X0, N0), (X1, N1), ...]) = X0^N0 X1^N1 mod B.m\n\
1453 \t(the list may be flattened if this more convenient.)")
1454 #undef METHNAME
1455   { 0 }
1456 };
1457
1458 static PyTypeObject *mpbarrett_pytype, mpbarrett_pytype_skel = {
1459   PyObject_HEAD_INIT(0) 0,              /* Header */
1460   "catacomb.MPBarrett",                 /* @tp_name@ */
1461   sizeof(mpbarrett_pyobj),              /* @tp_basicsize@ */
1462   0,                                    /* @tp_itemsize@ */
1463
1464   mpbarrett_pydealloc,                  /* @tp_dealloc@ */
1465   0,                                    /* @tp_print@ */
1466   0,                                    /* @tp_getattr@ */
1467   0,                                    /* @tp_setattr@ */
1468   0,                                    /* @tp_compare@ */
1469   0,                                    /* @tp_repr@ */
1470   0,                                    /* @tp_as_number@ */
1471   0,                                    /* @tp_as_sequence@ */
1472   0,                                    /* @tp_as_mapping@ */
1473   0,                                    /* @tp_hash@ */
1474   0,                                    /* @tp_call@ */
1475   0,                                    /* @tp_str@ */
1476   0,                                    /* @tp_getattro@ */
1477   0,                                    /* @tp_setattro@ */
1478   0,                                    /* @tp_as_buffer@ */
1479   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1480     Py_TPFLAGS_BASETYPE,
1481
1482   /* @tp_doc@ */
1483 "A Barrett reduction context.",
1484
1485   0,                                    /* @tp_traverse@ */
1486   0,                                    /* @tp_clear@ */
1487   0,                                    /* @tp_richcompare@ */
1488   0,                                    /* @tp_weaklistoffset@ */
1489   0,                                    /* @tp_iter@ */
1490   0,                                    /* @tp_iternext@ */
1491   mpbarrett_pymethods,                  /* @tp_methods@ */
1492   0,                                    /* @tp_members@ */
1493   mpbarrett_pygetset,                   /* @tp_getset@ */
1494   0,                                    /* @tp_base@ */
1495   0,                                    /* @tp_dict@ */
1496   0,                                    /* @tp_descr_get@ */
1497   0,                                    /* @tp_descr_set@ */
1498   0,                                    /* @tp_dictoffset@ */
1499   0,                                    /* @tp_init@ */
1500   PyType_GenericAlloc,                  /* @tp_alloc@ */
1501   mpbarrett_pynew,                      /* @tp_new@ */
1502   0,                                    /* @tp_free@ */
1503   0                                     /* @tp_is_gc@ */
1504 };
1505
1506 /*----- Nice prime reduction ----------------------------------------------*/
1507
1508 typedef struct mpreduce_pyobj {
1509   PyObject_HEAD
1510   mpreduce mr;
1511 } mpreduce_pyobj;
1512
1513 #define MPREDUCE_PY(o) (&((mpreduce_pyobj *)(o))->mr)
1514
1515 static PyObject *mrmeth_exp(PyObject *me, PyObject *arg)
1516 {
1517   PyObject *rc = 0;
1518   mp *yy = 0, *zz = 0;
1519
1520   if (!PyArg_ParseTuple(arg, "O&O&:exp", convmp, &yy, convmp, &zz))
1521     goto end;
1522   if (MP_NEGP(zz)) {
1523     if ((yy = mp_modinv_checked(yy, yy, MPREDUCE_PY(me)->p)) == 0) goto end;
1524     zz = mp_neg(zz, zz);
1525   }
1526   rc = mp_pywrap(mpreduce_exp(MPREDUCE_PY(me), MP_NEW, yy, zz));
1527 end:
1528   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
1529   return (rc);
1530 }
1531
1532 static PyObject *mrmeth_reduce(PyObject *me, PyObject *arg)
1533 {
1534   PyObject *z = 0;
1535   mp *yy = 0;
1536
1537   if (!PyArg_ParseTuple(arg, "O&:reduce", convmp, &yy)) goto end;
1538   z = mp_pywrap(mpreduce_do(MPREDUCE_PY(me), MP_NEW, yy));
1539 end:
1540   return (z);
1541 }
1542
1543 static void mpreduce_pydealloc(PyObject *me)
1544 {
1545   mpreduce_destroy(MPREDUCE_PY(me));
1546   FREEOBJ(me);
1547 }
1548
1549 static PyObject *mpreduce_pynew(PyTypeObject *ty,
1550                                  PyObject *arg, PyObject *kw)
1551 {
1552   mpreduce_pyobj *mr = 0;
1553   mpreduce r;
1554   char *kwlist[] = { "m", 0 };
1555   mp *xx = 0;
1556
1557   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convmp, &xx))
1558     goto end;
1559   if (!MP_POSP(xx)) VALERR("m must be positive");
1560   if (mpreduce_create(&r, xx)) VALERR("bad modulus (must be 2^k - ...)");
1561   mr = (mpreduce_pyobj *)ty->tp_alloc(ty, 0);
1562   mr->mr = r;
1563 end:
1564   if (xx) MP_DROP(xx);
1565   return ((PyObject *)mr);
1566 }
1567
1568 static PyObject *mrget_m(PyObject *me, void *hunoz)
1569   { return (mp_pywrap(MP_COPY(MPREDUCE_PY(me)->p))); }
1570
1571 static PyGetSetDef mpreduce_pygetset[] = {
1572 #define GETSETNAME(op, name) mr##op##_##name
1573   GET   (m,             "R.m -> modulus for reduction")
1574 #undef GETSETNAME
1575   { 0 }
1576 };
1577
1578 static PyMethodDef mpreduce_pymethods[] = {
1579 #define METHNAME(name) mrmeth_##name
1580   METH  (reduce,        "R.reduce(X) -> X mod B.m")
1581   METH  (exp,           "R.exp(X, N) -> X^N mod B.m")
1582 #undef METHNAME
1583   { 0 }
1584 };
1585
1586 static PyTypeObject *mpreduce_pytype, mpreduce_pytype_skel = {
1587   PyObject_HEAD_INIT(0) 0,              /* Header */
1588   "catacomb.MPReduce",                  /* @tp_name@ */
1589   sizeof(mpreduce_pyobj),               /* @tp_basicsize@ */
1590   0,                                    /* @tp_itemsize@ */
1591
1592   mpreduce_pydealloc,                   /* @tp_dealloc@ */
1593   0,                                    /* @tp_print@ */
1594   0,                                    /* @tp_getattr@ */
1595   0,                                    /* @tp_setattr@ */
1596   0,                                    /* @tp_compare@ */
1597   0,                                    /* @tp_repr@ */
1598   0,                                    /* @tp_as_number@ */
1599   0,                                    /* @tp_as_sequence@ */
1600   0,                                    /* @tp_as_mapping@ */
1601   0,                                    /* @tp_hash@ */
1602   0,                                    /* @tp_call@ */
1603   0,                                    /* @tp_str@ */
1604   0,                                    /* @tp_getattro@ */
1605   0,                                    /* @tp_setattro@ */
1606   0,                                    /* @tp_as_buffer@ */
1607   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1608     Py_TPFLAGS_BASETYPE,
1609
1610   /* @tp_doc@ */
1611 "A reduction context for reduction modulo primes of special form.",
1612
1613   0,                                    /* @tp_traverse@ */
1614   0,                                    /* @tp_clear@ */
1615   0,                                    /* @tp_richcompare@ */
1616   0,                                    /* @tp_weaklistoffset@ */
1617   0,                                    /* @tp_iter@ */
1618   0,                                    /* @tp_iternext@ */
1619   mpreduce_pymethods,                   /* @tp_methods@ */
1620   0,                                    /* @tp_members@ */
1621   mpreduce_pygetset,                    /* @tp_getset@ */
1622   0,                                    /* @tp_base@ */
1623   0,                                    /* @tp_dict@ */
1624   0,                                    /* @tp_descr_get@ */
1625   0,                                    /* @tp_descr_set@ */
1626   0,                                    /* @tp_dictoffset@ */
1627   0,                                    /* @tp_init@ */
1628   PyType_GenericAlloc,                  /* @tp_alloc@ */
1629   mpreduce_pynew,                       /* @tp_new@ */
1630   0,                                    /* @tp_free@ */
1631   0                                     /* @tp_is_gc@ */
1632 };
1633
1634 /*----- Chinese Remainder Theorem solution --------------------------------*/
1635
1636 typedef struct mpcrt_pyobj {
1637   PyObject_HEAD
1638   mpcrt c;
1639 } mpcrt_pyobj;
1640
1641 #define MPCRT_PY(o) (&((mpcrt_pyobj *)(o))->c)
1642
1643 static PyObject *mcmeth_solve(PyObject *me, PyObject *arg)
1644 {
1645   mpcrt *c = MPCRT_PY(me);
1646   PyObject *q = 0, *x, *z = 0;
1647   mp *xx;
1648   mp **v = 0;
1649   int i = 0, n = c->k;
1650
1651   Py_INCREF(me);
1652   if (PyTuple_Size(arg) == n)
1653     q = arg;
1654   else if (!PyArg_ParseTuple(arg, "O:solve", &q))
1655     goto end;
1656   Py_INCREF(q);
1657   if (!PySequence_Check(q)) TYERR("want a sequence of residues");
1658   if (PySequence_Size(q) != n) VALERR("residue count mismatch");
1659   v = xmalloc(n * sizeof(*v));
1660   for (i = 0; i < n; i++) {
1661     if ((x = PySequence_GetItem(q, i)) == 0) goto end;
1662     xx = getmp(x); Py_DECREF(x); if (!xx) goto end;
1663     v[i] = xx;
1664   }
1665   z = mp_pywrap(mpcrt_solve(c, MP_NEW, v));
1666 end:
1667   if (v) {
1668     n = i;
1669     for (i = 0; i < n; i++)
1670       MP_DROP(v[i]);
1671     xfree(v);
1672   }
1673   Py_DECREF(me);
1674   Py_XDECREF(q);
1675   return (z);
1676 }
1677
1678 static void mpcrt_pydealloc(PyObject *me)
1679 {
1680   mpcrt *c = MPCRT_PY(me);
1681   mpcrt_destroy(c);
1682   xfree(c->v);
1683 }
1684
1685 static PyObject *mpcrt_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1686 {
1687   mpcrt_mod *v = 0;
1688   int n, i = 0;
1689   char *kwlist[] = { "mv", 0 };
1690   PyObject *q = 0, *x;
1691   mp *xx;
1692   mpcrt_pyobj *c = 0;
1693
1694   if (PyTuple_Size(arg) > 1)
1695     q = arg;
1696   else if (!PyArg_ParseTupleAndKeywords(arg, kw, "O:new", kwlist, &q))
1697     goto end;
1698   Py_INCREF(q);
1699   if (!PySequence_Check(q)) TYERR("want a sequence of moduli");
1700   n = PySequence_Size(q);
1701   if (PyErr_Occurred()) goto end;
1702   if (!n) VALERR("want at least one modulus");
1703   v = xmalloc(n * sizeof(*v));
1704   for (i = 0; i < n; i++) {
1705     if ((x = PySequence_GetItem(q, i)) == 0) goto end;
1706     xx = getmp(x); Py_DECREF(x); if (!xx) goto end;
1707     v[i].m = xx; v[i].n = 0; v[i].ni = 0; v[i].nni = 0;
1708   }
1709   c = (mpcrt_pyobj *)ty->tp_alloc(ty, 0);
1710   mpcrt_create(&c->c, v, n, 0);
1711   Py_DECREF(q);
1712   return ((PyObject *)c);
1713
1714 end:
1715   if (v) {
1716     n = i;
1717     for (i = 0; i < n; i++)
1718       MP_DROP(v[i].m);
1719     xfree(v);
1720   }
1721   Py_XDECREF(q);
1722   return (0);
1723 }
1724
1725 static PyObject *mcget_product(PyObject *me, void *hunoz)
1726   { return (mp_pywrap(MP_COPY(MPCRT_PY(me)->mb.m))); }
1727
1728 static PyObject *mcget_moduli(PyObject *me, void *hunoz)
1729 {
1730   int i;
1731   PyObject *q;
1732   mpcrt *c = MPCRT_PY(me);
1733
1734   if ((q = PyList_New(c->k)) == 0) return (0);
1735   for (i = 0; i < c->k; i++)
1736     PyList_SetItem(q, i, mp_pywrap(c->v[i].m));
1737   return (q);
1738 }
1739
1740 static PyGetSetDef mpcrt_pygetset[] = {
1741 #define GETSETNAME(op, name) mc##op##_##name
1742   GET   (product,       "C.product -> product of moduli")
1743   GET   (moduli,        "C.moduli -> list of individual moduli")
1744 #undef GETSETNAME
1745   { 0 }
1746 };
1747
1748 static PyMethodDef mpcrt_pymethods[] = {
1749 #define METHNAME(name) mcmeth_##name
1750   METH  (solve,         "C.solve([R0, R1]) -> X mod C.product")
1751 #undef METHNAME
1752   { 0 }
1753 };
1754
1755 static PyTypeObject *mpcrt_pytype, mpcrt_pytype_skel = {
1756   PyObject_HEAD_INIT(0) 0,              /* Header */
1757   "catacomb.MPCRT",                     /* @tp_name@ */
1758   sizeof(mpcrt_pyobj),                  /* @tp_basicsize@ */
1759   0,                                    /* @tp_itemsize@ */
1760
1761   mpcrt_pydealloc,                      /* @tp_dealloc@ */
1762   0,                                    /* @tp_print@ */
1763   0,                                    /* @tp_getattr@ */
1764   0,                                    /* @tp_setattr@ */
1765   0,                                    /* @tp_compare@ */
1766   0,                                    /* @tp_repr@ */
1767   0,                                    /* @tp_as_number@ */
1768   0,                                    /* @tp_as_sequence@ */
1769   0,                                    /* @tp_as_mapping@ */
1770   0,                                    /* @tp_hash@ */
1771   0,                                    /* @tp_call@ */
1772   0,                                    /* @tp_str@ */
1773   0,                                    /* @tp_getattro@ */
1774   0,                                    /* @tp_setattro@ */
1775   0,                                    /* @tp_as_buffer@ */
1776   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
1777     Py_TPFLAGS_BASETYPE,
1778
1779   /* @tp_doc@ */
1780 "A context for the solution of Chinese Remainder Theorem problems.",
1781
1782   0,                                    /* @tp_traverse@ */
1783   0,                                    /* @tp_clear@ */
1784   0,                                    /* @tp_richcompare@ */
1785   0,                                    /* @tp_weaklistoffset@ */
1786   0,                                    /* @tp_iter@ */
1787   0,                                    /* @tp_iternext@ */
1788   mpcrt_pymethods,                      /* @tp_methods@ */
1789   0,                                    /* @tp_members@ */
1790   mpcrt_pygetset,                       /* @tp_getset@ */
1791   0,                                    /* @tp_base@ */
1792   0,                                    /* @tp_dict@ */
1793   0,                                    /* @tp_descr_get@ */
1794   0,                                    /* @tp_descr_set@ */
1795   0,                                    /* @tp_dictoffset@ */
1796   0,                                    /* @tp_init@ */
1797   PyType_GenericAlloc,                  /* @tp_alloc@ */
1798   mpcrt_pynew,                          /* @tp_new@ */
1799   0,                                    /* @tp_free@ */
1800   0                                     /* @tp_is_gc@ */
1801 };
1802
1803 /*----- Binary polynomials ------------------------------------------------*/
1804
1805 static PyObject *gf_pyrepr(PyObject *o)
1806   { return mp_topystring(MP_X(o), 16, "GF(", "0x", "L)"); }
1807
1808 static PyObject *gf_pyrichcompare(PyObject *x, PyObject *y, int op)
1809 {
1810   mp *xx, *yy;
1811   int xl, yl;
1812   int b;
1813
1814   if (gfbinop(x, y, &xx, &yy)) RETURN_NOTIMPL;
1815   switch (op) {
1816     case Py_EQ: b = MP_EQ(xx, yy); break;
1817     case Py_NE: b = !MP_EQ(xx, yy); break;
1818     default:
1819       xl = mp_bits(xx);
1820       yl = mp_bits(yy);
1821       switch (op) {
1822         case Py_LT: b = xl < yl; break;
1823         case Py_LE: b = xl <= yl; break;
1824         case Py_GT: b = xl > yl; break;
1825         case Py_GE: b = xl >= yl; break;
1826         default: abort();
1827       }
1828       break;
1829   }
1830   MP_DROP(xx); MP_DROP(yy);
1831   return (getbool(b));
1832 }
1833
1834 static PyObject *gf_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
1835 {
1836   PyObject *x;
1837   mp *z;
1838   mp_pyobj *zz = 0;
1839   int radix = 0;
1840   char *kwlist[] = { "x", "radix", 0 };
1841
1842   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O|i:gf", kwlist, &x, &radix))
1843     goto end;
1844   if (GF_PYCHECK(x)) RETURN_OBJ(x);
1845   if (!good_radix_p(radix, 1)) VALERR("radix out of range");
1846   if ((z = mp_frompyobject(x, radix)) == 0) {
1847     PyErr_Format(PyExc_TypeError, "can't convert %.100s to gf",
1848                  x->ob_type->tp_name);
1849     goto end;
1850   }
1851   if (MP_NEGP(z)) {
1852     MP_DROP(z);
1853     VALERR("gf cannot be negative");
1854   }
1855   zz = (mp_pyobj *)ty->tp_alloc(ty, 0);
1856   zz->x = z;
1857 end:
1858   return ((PyObject *)zz);
1859 }
1860
1861 static long gf_pyhash(PyObject *me)
1862 {
1863   long i = mp_tolong(MP_X(me));
1864   i ^= 0xc7ecd67c; /* random perturbance */
1865   if (i == -1)
1866     i = -2;
1867   return (i);
1868 }
1869
1870 static PyObject *gf_pyexp(PyObject *x, PyObject *y, PyObject *z)
1871 {
1872   mp *xx = 0, *yy = 0, *zz = 0;
1873   mp *r = 0;
1874   PyObject *rc = 0;
1875
1876   if ((xx = tomp(x)) == 0 || (yy = tomp(y)) == 0 ||
1877       (z && z != Py_None && (zz = tomp(z)) == 0)) {
1878     mp_drop(xx); mp_drop(yy); mp_drop(zz);
1879     RETURN_NOTIMPL;
1880   }
1881   if (!z || z == Py_None) {
1882     if (MP_NEGP(yy)) VALERR("negative exponent");
1883     r = gf_exp(MP_NEW, xx, yy);
1884   } else {
1885     gfreduce gr;
1886     if (MP_ZEROP(zz)) ZDIVERR("zero modulus");
1887     if (MP_NEGP(yy)) {
1888       if ((xx = gf_modinv_checked(xx, xx, zz)) == 0) goto end;
1889       yy = mp_neg(yy, yy);
1890     }
1891     gfreduce_create(&gr, zz);
1892     r = gfreduce_exp(&gr, MP_NEW, xx, yy);
1893     gfreduce_destroy(&gr);
1894   }
1895   rc = gf_pywrap(r);
1896 end:
1897   mp_drop(xx); mp_drop(yy); mp_drop(zz);
1898   return (rc);
1899 }
1900
1901 static PyObject *gfmeth_sqr(PyObject *me, PyObject *arg)
1902 {
1903   if (!PyArg_ParseTuple(arg, ":sqr")) return (0);
1904   return (gf_pywrap(gf_sqr(MP_NEW, MP_X(me))));
1905 }
1906
1907 static PyObject *gfmeth_gcd(PyObject *me, PyObject *arg)
1908 {
1909   PyObject *z = 0;
1910   mp *yy = 0, *zz = MP_NEW;
1911
1912   if (!PyArg_ParseTuple(arg, "O&:gcd", convgf, &yy)) goto end;
1913   gf_gcd(&zz, 0, 0, MP_X(me), yy);
1914   z = gf_pywrap(zz);
1915 end:
1916   if (yy) MP_DROP(yy);
1917   return (z);
1918 }
1919
1920 static PyObject *gfmeth_gcdx(PyObject *me, PyObject *arg)
1921 {
1922   PyObject *z = 0;
1923   mp *yy = 0, *zz = MP_NEW, *uu = MP_NEW, *vv = MP_NEW;
1924
1925   if (!PyArg_ParseTuple(arg, "O&:gcdx", convgf, &yy))
1926     goto end;
1927   gf_gcd(&zz, &uu, &vv, MP_X(me), yy);
1928   z = Py_BuildValue("(NNN)",
1929                     gf_pywrap(zz), gf_pywrap(uu), gf_pywrap(vv));
1930 end:
1931   if (yy) MP_DROP(yy);
1932   return (z);
1933 }
1934
1935 static PyObject *gfmeth_modinv(PyObject *me, PyObject *arg)
1936 {
1937   PyObject *z = 0;
1938   mp *yy = 0, *zz = MP_NEW;
1939
1940   if (!PyArg_ParseTuple(arg, "O&:modinv", convgf, &yy) ||
1941       (zz = gf_modinv_checked(MP_NEW, yy, MP_X(me))) == 0)
1942     goto end;
1943   z = gf_pywrap(zz);
1944 end:
1945   if (yy) MP_DROP(yy);
1946   return (z);
1947 }
1948
1949 static PyObject *gfmeth_irreduciblep(PyObject *me, PyObject *arg)
1950 {
1951   if (!PyArg_ParseTuple(arg, ":irreduciblep")) return (0);
1952   return getbool(gf_irreduciblep(MP_X(me)));
1953 }
1954
1955 static PyObject *gfget_degree(PyObject *me, void *hunoz)
1956   { return (PyInt_FromLong(mp_bits(MP_X(me)) - 1)); }
1957
1958 static PyGetSetDef gf_pygetset[] = {
1959 #define GETSETNAME(op, name) gf##op##_##name
1960   GET   (degree,        "X.degree -> polynomial degree of X")
1961 #undef GETSETNAME
1962 #define GETSETNAME(op, name) mp##op##_##name
1963   GET   (nbits,         "X.nbits -> bit length of X")
1964   GET   (noctets,       "X.noctets -> octet length of X")
1965 #undef GETSETNAME
1966   { 0 }
1967 };
1968
1969 static PyMethodDef gf_pymethods[] = {
1970 #define METHNAME(func) gfmeth_##func
1971   METH  (setbit,        "X.setbit(N) -> X with bit N set")
1972   METH  (clearbit,      "X.clearbit(N) -> X with bit N clear")
1973   METH  (testbit,       "X.testbit(N) -> true/false if bit N set/clear in X")
1974   METH  (sqr,           "X.sqr() -> X^2")
1975   METH  (gcd,           "X.gcd(Y) -> gcd(X, Y)")
1976   METH  (gcdx,
1977          "X.gcdx(Y) -> (gcd(X, Y), U, V) with X U + Y V = gcd(X, Y)")
1978   METH  (modinv,        "X.modinv(Y) -> multiplicative inverse of Y mod X")
1979   METH  (irreduciblep,  "X.irreduciblep() -> true/false")
1980 #undef METHNAME
1981 #define METHNAME(func) mpmeth_##func
1982   KWMETH(tostring,      "X.tostring(radix = 10) -> STR")
1983   KWMETH(storel,        "X.storel(len = -1) -> little-endian bytes")
1984   KWMETH(storeb,        "X.storeb(len = -1) -> big-endian bytes")
1985   KWMETH(storel2c,
1986          "X.storel2c(len = -1) -> little-endian bytes, two's complement")
1987   KWMETH(storeb2c,
1988          "X.storeb2c(len = -1) -> big-endian bytes, two's complement")
1989   METH  (tobuf,         "X.tobuf() -> buffer format")
1990 #undef METHNAME
1991   { 0 }
1992 };
1993
1994 static PyNumberMethods gf_pynumber = {
1995   gf_pyadd,                             /* @nb_add@ */
1996   gf_pysub,                             /* @nb_subtract@ */
1997   gf_pymul,                             /* @nb_multiply@ */
1998   0,                                    /* @nb_divide@ */
1999   gf_pymod,                             /* @nb_remainder@ */
2000   gf_pydivmod,                          /* @nb_divmod@ */
2001   gf_pyexp,                             /* @nb_power@ */
2002   mp_pyid,                              /* @nb_negative@ */
2003   mp_pyid,                              /* @nb_positive@ */
2004   mp_pyid,                              /* @nb_absolute@ */
2005   mp_pynonzerop,                        /* @nb_nonzero@ */
2006   0 /* doesn't make any sense */,       /* @nb_invert@ */
2007   gf_pylsl,                             /* @nb_lshift@ */
2008   gf_pylsr,                             /* @nb_rshift@ */
2009   gf_pyand,                             /* @nb_and@ */
2010   gf_pyxor,                             /* @nb_xor@ */
2011   gf_pyor,                              /* @nb_or@ */
2012   gf_pycoerce,                          /* @nb_coerce@ */
2013   mp_pyint,                             /* @nb_int@ */
2014   mp_pylong,                            /* @nb_long@ */
2015   0 /* doesn't make any sense */,       /* @nb_float@ */
2016   mp_pyoct,                             /* @nb_oct@ */
2017   mp_pyhex,                             /* @nb_hex@ */
2018
2019   0,                                    /* @nb_inplace_add@ */
2020   0,                                    /* @nb_inplace_subtract@ */
2021   0,                                    /* @nb_inplace_multiply@ */
2022   0,                                    /* @nb_inplace_divide@ */
2023   0,                                    /* @nb_inplace_remainder@ */
2024   0,                                    /* @nb_inplace_power@ */
2025   0,                                    /* @nb_inplace_lshift@ */
2026   0,                                    /* @nb_inplace_rshift@ */
2027   0,                                    /* @nb_inplace_and@ */
2028   0,                                    /* @nb_inplace_xor@ */
2029   0,                                    /* @nb_inplace_or@ */
2030
2031   gf_pydiv,                             /* @nb_floor_divide@ */
2032   0,                                    /* @nb_true_divide@ */
2033   0,                                    /* @nb_inplace_floor_divide@ */
2034   0,                                    /* @nb_inplace_true_divide@ */
2035 };
2036
2037 static PyTypeObject gf_pytype_skel = {
2038   PyObject_HEAD_INIT(0) 0,              /* Header */
2039   "catacomb.GF",                        /* @tp_name@ */
2040   sizeof(mp_pyobj),                     /* @tp_basicsize@ */
2041   0,                                    /* @tp_itemsize@ */
2042
2043   mp_pydealloc,                         /* @tp_dealloc@ */
2044   0,                                    /* @tp_print@ */
2045   0,                                    /* @tp_getattr@ */
2046   0,                                    /* @tp_setattr@ */
2047   0,                                    /* @tp_compare@ */
2048   gf_pyrepr,                            /* @tp_repr@ */
2049   &gf_pynumber,                         /* @tp_as_number@ */
2050   0,                                    /* @tp_as_sequence@ */
2051   0,                                    /* @tp_as_mapping@ */
2052   gf_pyhash,                            /* @tp_hash@ */
2053   0,                                    /* @tp_call@ */
2054   mp_pyhex,                             /* @tp_str@ */
2055   0,                                    /* @tp_getattro@ */
2056   0,                                    /* @tp_setattro@ */
2057   0,                                    /* @tp_as_buffer@ */
2058   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
2059     Py_TPFLAGS_CHECKTYPES |
2060     Py_TPFLAGS_BASETYPE,
2061
2062   /* @tp_doc@ */
2063 "Binary polynomials.  Support almost all the standard arithmetic\n\
2064 operations.\n\
2065 \n\
2066 Constructor gf(X, radix = R) attempts to convert X to a `gf'.  If\n\
2067 X is a string, it's read in radix-R form, or we look for a prefix\n\
2068 if R = 0.  Other acceptable things are ints and longs.\n\
2069 \n\
2070 The name is hopelessly wrong from a technical point of view, but\n\
2071 but it's much easier to type than `p2' or `c2' or whatever.\n\
2072 \n\
2073 Notes:\n\
2074 \n\
2075   * Use `//' for division.  GFs don't have `/' division.",
2076
2077   0,                                    /* @tp_traverse@ */
2078   0,                                    /* @tp_clear@ */
2079   gf_pyrichcompare,                     /* @tp_richcompare@ */
2080   0,                                    /* @tp_weaklistoffset@ */
2081   0,                                    /* @tp_iter@ */
2082   0,                                    /* @tp_iternext@ */
2083   gf_pymethods,                         /* @tp_methods@ */
2084   0,                                    /* @tp_members@ */
2085   gf_pygetset,                          /* @tp_getset@ */
2086   0,                                    /* @tp_base@ */
2087   0,                                    /* @tp_dict@ */
2088   0,                                    /* @tp_descr_get@ */
2089   0,                                    /* @tp_descr_set@ */
2090   0,                                    /* @tp_dictoffset@ */
2091   0,                                    /* @tp_init@ */
2092   PyType_GenericAlloc,                  /* @tp_alloc@ */
2093   gf_pynew,                             /* @tp_new@ */
2094   0,                                    /* @tp_free@ */
2095   0                                     /* @tp_is_gc@ */
2096 };
2097
2098 static PyObject *meth__GF_fromstring(PyObject *me,
2099                                     PyObject *arg, PyObject *kw)
2100 {
2101   int r = 0;
2102   char *p;
2103   int len;
2104   PyObject *z = 0;
2105   mp *zz;
2106   mptext_stringctx sc;
2107   char *kwlist[] = { "class", "x", "radix", 0 };
2108
2109   if (!PyArg_ParseTupleAndKeywords(arg, kw, "Os#|i:fromstring",
2110                                    kwlist, &me, &p, &len, &r))
2111     goto end;
2112   if (!good_radix_p(r, 1)) VALERR("bad radix");
2113   sc.buf = p; sc.lim = p + len;
2114   if ((zz = mp_read(MP_NEW, r, &mptext_stringops, &sc)) == 0 ||
2115       MP_NEGP(zz)) {
2116     if (zz) MP_DROP(zz);
2117     SYNERR("bad binary polynomial");
2118   }
2119   z = Py_BuildValue("(Ns#)", gf_pywrap(zz), sc.buf, (int)(sc.lim - sc.buf));
2120 end:
2121   return (z);
2122 }
2123
2124 /*----- Sparse poly reduction ---------------------------------------------*/
2125
2126 typedef struct gfreduce_pyobj {
2127   PyObject_HEAD
2128   gfreduce mr;
2129 } gfreduce_pyobj;
2130
2131 #define GFREDUCE_PY(o) (&((gfreduce_pyobj *)(o))->mr)
2132
2133 static PyObject *grmeth_exp(PyObject *me, PyObject *arg)
2134 {
2135   PyObject *rc = 0;
2136   mp *yy = 0, *zz = 0;
2137
2138   if (!PyArg_ParseTuple(arg, "O&O&:exp", convgf, &yy, convgf, &zz))
2139     goto end;
2140   if (MP_NEGP(zz)) {
2141     if ((yy = gf_modinv_checked(yy, yy, GFREDUCE_PY(me)->p)) == 0) goto end;
2142     zz = mp_neg(zz, zz);
2143   }
2144   rc = gf_pywrap(gfreduce_exp(GFREDUCE_PY(me), MP_NEW, yy, zz));
2145 end:
2146   if (yy) MP_DROP(yy); if (zz) MP_DROP(zz);
2147   return (rc);
2148 }
2149
2150 static PyObject *grmeth_trace(PyObject *me, PyObject *arg)
2151 {
2152   PyObject *rc = 0;
2153   mp *xx = 0;
2154
2155   if (!PyArg_ParseTuple(arg, "O&:trace", convgf, &xx)) goto end;
2156   rc = PyInt_FromLong(gfreduce_trace(GFREDUCE_PY(me), xx));
2157 end:
2158   if (xx) MP_DROP(xx);
2159   return (rc);
2160 }
2161
2162 static PyObject *grmeth_halftrace(PyObject *me, PyObject *arg)
2163 {
2164   PyObject *rc = 0;
2165   mp *xx = 0;
2166
2167   if (!PyArg_ParseTuple(arg, "O&:halftrace", convgf, &xx)) goto end;
2168   rc = gf_pywrap(gfreduce_halftrace(GFREDUCE_PY(me), MP_NEW, xx));
2169 end:
2170   if (xx) MP_DROP(xx);
2171   return (rc);
2172 }
2173
2174 static PyObject *grmeth_sqrt(PyObject *me, PyObject *arg)
2175 {
2176   PyObject *rc = 0;
2177   mp *xx = 0, *yy;
2178
2179   if (!PyArg_ParseTuple(arg, "O&:sqrt", convgf, &xx)) goto end;
2180   if ((yy = gfreduce_sqrt(GFREDUCE_PY(me), MP_NEW, xx)) == 0)
2181     VALERR("no modular square root");
2182   rc = gf_pywrap(yy);
2183 end:
2184   if (xx) MP_DROP(xx);
2185   return (rc);
2186 }
2187
2188 static PyObject *grmeth_quadsolve(PyObject *me, PyObject *arg)
2189 {
2190   PyObject *rc = 0;
2191   mp *xx = 0, *yy;
2192
2193   if (!PyArg_ParseTuple(arg, "O&:quadsolve", convgf, &xx)) goto end;
2194   if ((yy = gfreduce_quadsolve(GFREDUCE_PY(me), MP_NEW, xx)) == 0)
2195     VALERR("no solution found");
2196   rc = gf_pywrap(yy);
2197 end:
2198   if (xx) MP_DROP(xx);
2199   return (rc);
2200 }
2201
2202 static PyObject *grmeth_reduce(PyObject *me, PyObject *arg)
2203 {
2204   PyObject *z = 0;
2205   mp *yy = 0;
2206
2207   if (!PyArg_ParseTuple(arg, "O&:reduce", convgf, &yy)) goto end;
2208   z = gf_pywrap(gfreduce_do(GFREDUCE_PY(me), MP_NEW, yy));
2209 end:
2210   return (z);
2211 }
2212
2213 static void gfreduce_pydealloc(PyObject *me)
2214 {
2215   gfreduce_destroy(GFREDUCE_PY(me));
2216   FREEOBJ(me);
2217 }
2218
2219 static PyObject *gfreduce_pynew(PyTypeObject *ty,
2220                                  PyObject *arg, PyObject *kw)
2221 {
2222   gfreduce_pyobj *mr = 0;
2223   gfreduce r;
2224   char *kwlist[] = { "m", 0 };
2225   mp *xx = 0;
2226
2227   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&:new", kwlist, convgf, &xx))
2228     goto end;
2229   if (MP_ZEROP(xx)) ZDIVERR("modulus is zero!");
2230   gfreduce_create(&r, xx);
2231   mr = (gfreduce_pyobj *)ty->tp_alloc(ty, 0);
2232   mr->mr = r;
2233 end:
2234   if (xx) MP_DROP(xx);
2235   return ((PyObject *)mr);
2236 }
2237
2238 static PyObject *grget_m(PyObject *me, void *hunoz)
2239   { return (gf_pywrap(MP_COPY(GFREDUCE_PY(me)->p))); }
2240
2241 static PyGetSetDef gfreduce_pygetset[] = {
2242 #define GETSETNAME(op, name) gr##op##_##name
2243   GET   (m,             "R.m -> reduction polynomial")
2244 #undef GETSETNAME
2245   { 0 }
2246 };
2247
2248 static PyMethodDef gfreduce_pymethods[] = {
2249 #define METHNAME(name) grmeth_##name
2250   METH  (reduce,        "R.reduce(X) -> X mod B.m")
2251   METH  (trace,        "R.trace(X) -> Tr(X) = x + x^2 + ... + x^{2^{m - 1}}")
2252   METH  (halftrace,   "R.halftrace(X) -> x + x^{2^2} + ... + x^{2^{m - 1}}")
2253   METH  (sqrt,          "R.sqrt(X) -> Y where Y^2 = X mod R")
2254   METH  (quadsolve,     "R.quadsolve(X) -> Y where Y^2 + Y = X mod R")
2255   METH  (exp,           "R.exp(X, N) -> X^N mod B.m")
2256 #undef METHNAME
2257   { 0 }
2258 };
2259
2260 static PyTypeObject *gfreduce_pytype, gfreduce_pytype_skel = {
2261   PyObject_HEAD_INIT(0) 0,              /* Header */
2262   "catacomb.GFReduce",                  /* @tp_name@ */
2263   sizeof(gfreduce_pyobj),               /* @tp_basicsize@ */
2264   0,                                    /* @tp_itemsize@ */
2265
2266   gfreduce_pydealloc,                   /* @tp_dealloc@ */
2267   0,                                    /* @tp_print@ */
2268   0,                                    /* @tp_getattr@ */
2269   0,                                    /* @tp_setattr@ */
2270   0,                                    /* @tp_compare@ */
2271   0,                                    /* @tp_repr@ */
2272   0,                                    /* @tp_as_number@ */
2273   0,                                    /* @tp_as_sequence@ */
2274   0,                                    /* @tp_as_mapping@ */
2275   0,                                    /* @tp_hash@ */
2276   0,                                    /* @tp_call@ */
2277   0,                                    /* @tp_str@ */
2278   0,                                    /* @tp_getattro@ */
2279   0,                                    /* @tp_setattro@ */
2280   0,                                    /* @tp_as_buffer@ */
2281   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
2282     Py_TPFLAGS_BASETYPE,
2283
2284   /* @tp_doc@ */
2285 "A reduction context for reduction modulo sparse irreducible polynomials.",
2286
2287   0,                                    /* @tp_traverse@ */
2288   0,                                    /* @tp_clear@ */
2289   0,                                    /* @tp_richcompare@ */
2290   0,                                    /* @tp_weaklistoffset@ */
2291   0,                                    /* @tp_iter@ */
2292   0,                                    /* @tp_iternext@ */
2293   gfreduce_pymethods,                   /* @tp_methods@ */
2294   0,                                    /* @tp_members@ */
2295   gfreduce_pygetset,                    /* @tp_getset@ */
2296   0,                                    /* @tp_base@ */
2297   0,                                    /* @tp_dict@ */
2298   0,                                    /* @tp_descr_get@ */
2299   0,                                    /* @tp_descr_set@ */
2300   0,                                    /* @tp_dictoffset@ */
2301   0,                                    /* @tp_init@ */
2302   PyType_GenericAlloc,                  /* @tp_alloc@ */
2303   gfreduce_pynew,                       /* @tp_new@ */
2304   0,                                    /* @tp_free@ */
2305   0                                     /* @tp_is_gc@ */
2306 };
2307
2308 /*----- Normal/poly transformation ----------------------------------------*/
2309
2310 typedef struct gfn_pyobj {
2311   PyObject_HEAD
2312   mp *p;
2313   gfn ntop, pton;
2314 } gfn_pyobj;
2315
2316 static PyTypeObject *gfn_pytype, gfn_pytype_skel;
2317
2318 #define GFN_P(o) (((gfn_pyobj *)(o))->p)
2319 #define GFN_PTON(o) (&((gfn_pyobj *)(o))->pton)
2320 #define GFN_NTOP(o) (&((gfn_pyobj *)(o))->ntop)
2321
2322 static PyObject *gfn_pynew(PyTypeObject *ty, PyObject *arg, PyObject *kw)
2323 {
2324   mp *p = 0, *beta = 0;
2325   gfn_pyobj *gg = 0;
2326   char *kwlist[] = { "p", "beta", 0 };
2327
2328   if (!PyArg_ParseTupleAndKeywords(arg, kw, "O&O&:new", kwlist,
2329                                    convgf, &p, convgf, &beta))
2330     goto end;
2331   gg = PyObject_New(gfn_pyobj, ty);
2332   if (gfn_create(p, beta, &gg->ntop, &gg->pton)) {
2333     FREEOBJ(gg);
2334     gg = 0;
2335     VALERR("can't invert transformation matrix");
2336   }
2337   gg->p = MP_COPY(p);
2338 end:
2339   mp_drop(p);
2340   mp_drop(beta);
2341   return ((PyObject *)gg);
2342 }
2343
2344 static PyObject *gfnget_p(PyObject *me, void *hunoz)
2345   { return (gf_pywrap(MP_COPY(GFN_P(me)))); }
2346
2347 static PyObject *gfnget_beta(PyObject *me, void *hunoz)
2348 {
2349   gfn *n = GFN_NTOP(me);
2350   mp *x = n->r[n->n - 1];
2351   return (gf_pywrap(MP_COPY(x)));
2352 }
2353
2354 #define XFORMOP(name, NAME)                                             \
2355   static PyObject *gfnmeth_##name(PyObject *me, PyObject *arg)          \
2356   {                                                                     \
2357     mp *xx = 0;                                                         \
2358     mp *z = 0;                                                          \
2359                                                                         \
2360     if (!PyArg_ParseTuple(arg, "O&:" #name, convgf, &xx)) goto end;     \
2361     z = gfn_transform(GFN_##NAME(me), MP_NEW, xx);                      \
2362   end:                                                                  \
2363     mp_drop(xx);                                                        \
2364     if (!z) return (0);                                                 \
2365     return (mp_pywrap(z));                                              \
2366   }
2367 XFORMOP(pton, PTON)
2368 XFORMOP(ntop, NTOP)
2369 #undef XFORMOP
2370
2371 static void gfn_pydealloc(PyObject *me)
2372 {
2373   gfn_destroy(GFN_PTON(me));
2374   gfn_destroy(GFN_NTOP(me));
2375   FREEOBJ(me);
2376 }
2377
2378 static PyGetSetDef gfn_pygetset[] = {
2379 #define GETSETNAME(op, name) gfn##op##_##name
2380   GET   (p,             "X.p -> polynomial basis, as polynomial")
2381   GET   (beta,          "X.beta -> normal basis element, in poly form")
2382 #undef GETSETNAME
2383   { 0 }
2384 };
2385
2386 static PyMethodDef gfn_pymethods[] = {
2387 #define METHNAME(name) gfnmeth_##name
2388   METH  (pton,          "X.pton(A) -> normal-basis representation of A")
2389   METH  (ntop,          "X.ntop(A) -> polynomial-basis representation of A")
2390 #undef METHNAME
2391   { 0 }
2392 };
2393
2394 static PyTypeObject gfn_pytype_skel = {
2395   PyObject_HEAD_INIT(0) 0,              /* Header */
2396   "catacomb.GFN",                       /* @tp_name@ */
2397   sizeof(gfn_pyobj),                    /* @tp_basicsize@ */
2398   0,                                    /* @tp_itemsize@ */
2399
2400   gfn_pydealloc,                        /* @tp_dealloc@ */
2401   0,                                    /* @tp_print@ */
2402   0,                                    /* @tp_getattr@ */
2403   0,                                    /* @tp_setattr@ */
2404   0,                                    /* @tp_compare@ */
2405   0,                                    /* @tp_repr@ */
2406   0,                                    /* @tp_as_number@ */
2407   0,                                    /* @tp_as_sequence@ */
2408   0,                                    /* @tp_as_mapping@ */
2409   0,                                    /* @tp_hash@ */
2410   0,                                    /* @tp_call@ */
2411   0,                                    /* @tp_str@ */
2412   0,                                    /* @tp_getattro@ */
2413   0,                                    /* @tp_setattro@ */
2414   0,                                    /* @tp_as_buffer@ */
2415   Py_TPFLAGS_DEFAULT |                  /* @tp_flags@ */
2416     Py_TPFLAGS_BASETYPE,
2417
2418   /* @tp_doc@ */
2419 "An object for transforming elements of binary fields between polynomial\n\
2420 and normal basis representations.",
2421
2422   0,                                    /* @tp_traverse@ */
2423   0,                                    /* @tp_clear@ */
2424   0,                                    /* @tp_richcompare@ */
2425   0,                                    /* @tp_weaklistoffset@ */
2426   0,                                    /* @tp_iter@ */
2427   0,                                    /* @tp_iternext@ */
2428   gfn_pymethods,                        /* @tp_methods@ */
2429   0,                                    /* @tp_members@ */
2430   gfn_pygetset,                         /* @tp_getset@ */
2431   0,                                    /* @tp_base@ */
2432   0,                                    /* @tp_dict@ */
2433   0,                                    /* @tp_descr_get@ */
2434   0,                                    /* @tp_descr_set@ */
2435   0,                                    /* @tp_dictoffset@ */
2436   0,                                    /* @tp_init@ */
2437   PyType_GenericAlloc,                  /* @tp_alloc@ */
2438   gfn_pynew,                            /* @tp_new@ */
2439   0,                                    /* @tp_free@ */
2440   0                                     /* @tp_is_gc@ */
2441 };
2442
2443 /*----- Glue --------------------------------------------------------------*/
2444
2445 static PyMethodDef methods[] = {
2446 #define METHNAME(func) meth_##func
2447   KWMETH(_MP_fromstring,        "\
2448 fromstring(STR, radix = 0) -> (X, REST)\n\
2449 \n\
2450 Parse STR as a large integer, according to radix.  If radix is zero,\n\
2451 read a prefix from STR to decide radix: allow `0' for octal, `0x' for hex\n\
2452 or `R_' for other radix R.")
2453   KWMETH(_GF_fromstring,        "\
2454 fromstring(STR, radix = 0) -> (X, REST)\n\
2455 \n\
2456 Parse STR as a binary polynomial, according to radix.  If radix is zero,\n\
2457 read a prefix from STR to decide radix: allow `0' for octal, `0x' for hex\n\
2458 or `R_' for other radix R.")
2459   METH  (_MP_factorial,         "\
2460 factorial(I) -> I!: compute factorial")
2461   METH  (_MP_loadl,             "\
2462 loadl(STR) -> X: read little-endian bytes")
2463   METH  (_MP_loadb,             "\
2464 loadb(STR) -> X: read big-endian bytes")
2465   METH  (_MP_loadl2c,           "\
2466 loadl2c(STR) -> X: read little-endian bytes, two's complement")
2467   METH  (_MP_loadb2c,           "\
2468 loadb2c(STR) -> X: read big-endian bytes, two's complement")
2469   METH  (_MP_frombuf,           "\
2470 frombuf(STR) -> (X, REST): read buffer format")
2471   METH  (_GF_loadl,             "\
2472 loadl(STR) -> X: read little-endian bytes")
2473   METH  (_GF_loadb,             "\
2474 loadb(STR) -> X: read big-endian bytes")
2475   METH  (_GF_frombuf,           "\
2476 frombuf(STR) -> (X, REST): read buffer format")
2477 #undef METHNAME
2478   { 0 }
2479 };
2480
2481 void mp_pyinit(void)
2482 {
2483   INITTYPE(mp, root);
2484   INITTYPE(gf, root);
2485   INITTYPE(mpmul, root);
2486   INITTYPE(mpmont, root);
2487   INITTYPE(mpbarrett, root);
2488   INITTYPE(mpreduce, root);
2489   INITTYPE(mpcrt, root);
2490   INITTYPE(gfreduce, root);
2491   INITTYPE(gfn, root);
2492   addmethods(methods);
2493 }
2494
2495 void mp_pyinsert(PyObject *mod)
2496 {
2497   INSERT("MP", mp_pytype);
2498   INSERT("MPMul", mpmul_pytype);
2499   INSERT("MPMont", mpmont_pytype);
2500   INSERT("MPBarrett", mpbarrett_pytype);
2501   INSERT("MPReduce", mpreduce_pytype);
2502   INSERT("MPCRT", mpcrt_pytype);
2503   INSERT("GF", gf_pytype);
2504   INSERT("GFReduce", gfreduce_pytype);
2505   INSERT("GFN", gfn_pytype);
2506 }
2507
2508 /*----- That's all, folks -------------------------------------------------*/