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

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