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