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

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