GMPY: divm() still broken (was: module conflict? gmpy and operator)

[mensanator at aol.com]

Drat, gmpy is still broken.

Prior to ver 1.01, divm(4,8,20) would produce a
'not invertible' error even though this is a valid
linear congruence (because gcd(8,20)=4 divides 4).
The issue was that the modular inverse function invert()
requires that 8 & 20 be co-prime (have gcd=1). Linear
congruence doesn't require that, only that the gcd
divides 4. But if the gcd divides 4, then all three
parameters are divisible by 4, so we can convert the
non-invertible (4,8,20) into the invertible (1,2,5).

That feature was added to the 1.01 version of divm()
in addition to fixing the memory leak.

And although it works, it has an "interesting" side
effect: completely corrupting the gmpy module.


IDLE 1.1a3
>>> from gmpy import *
>>> x = mpz('8')
>>> y = mpz('20')
>>> z = mpz('4')
>>> x,y,z
(mpz(8), mpz(20), mpz(4))
>>> divm(z,x,y)
mpz(3)

Correct answer...

>>> x,y,z
(mpz(2), mpz(5), mpz(1))

...but variables have changed (that's not supposed to
happen). 2,5,1 results from dividing x,y,z by gcd(8,20)
in order to obtain a valid modular inverse.

Invoking divm(z,x,y) again still produces the correct
answer.

>>> divm(z,x,y)
mpz(3)

But if x,y,z are actually 2,5,1, this should fail with
a 'not invertible' error. So even though it is being
called with x,y,z, the original values are still being
used somehow.

Calling divm() with literals...

>>> divm(4,8,20)
mpz(3)

...also produces the correct answer, but the literals
had to be coerced to mpzs.

Earlier, we saw that x,y,z somehow have two different
values simultaneously. With the literals coerced to mpzs,
we have corrupted gmpy so that it cannot do any further
coercion on 8:

>>> mpz(8)
mpz(2)

8 is now stuck with value 2 when coerced to mpz.

So even though divm(4,8,20) should work, it can't any
more since failure to coerce 4 & 8 means divm() is
actually trying to evaluate (1,2,20).

>>> divm(4,8,20)

Traceback (most recent call last):
  File "<pyshell#21>", line 1, in -toplevel-
    divm(4,8,20)
ZeroDivisionError: not invertible


My speculation is that the three mpz_divexact functions
shown below must be doing some kind of bizarre pointer
magic that permits x,y,z to have two different values
simultaneously while also corrupting the coercion of
literals. It's legal in the gmp program to have a result
overwrite an operand, but my guess is that shouldn't
be done.

static char doc_divm[]="\
divm(a,b,m): returns x such that b*x==a modulo m, or else raises\n\
a ZeroDivisionError exception if no such value x exists\n\
(a, b and m must be mpz objects, or else get coerced to mpz)\n\
";
static PyObject *
Pygmpy_divm(PyObject *self, PyObject *args)
{
    PympzObject *num, *den, *mod, *res;
    int ok;

    if(!PyArg_ParseTuple(args, "O&O&O&",
        Pympz_convert_arg, &num,
        Pympz_convert_arg, &den,
        Pympz_convert_arg, &mod))
    {
        return last_try("divm", 3, 3, args);
    }
    if(!(res = Pympz_new()))
        return NULL;
    if(mpz_invert(res->z, den->z, mod->z)) { /* inverse exists */
      ok = 1;
    } else {
      /* last-ditch attempt: do num, den AND mod have a gcd>1 ? */
      PympzObject *gcd;
      if(!(gcd = Pympz_new()))
          return NULL;
      mpz_gcd(gcd->z, num->z, den->z);
      mpz_gcd(gcd->z, gcd->z, mod->z);
      mpz_divexact(num->z, num->z, gcd->z);
      mpz_divexact(den->z, den->z, gcd->z);
      mpz_divexact(mod->z, mod->z, gcd->z);
      Py_DECREF((PyObject*)gcd);
      ok = mpz_invert(res->z, den->z, mod->z);
    }

    if (ok) {
        mpz_mul(res->z, res->z, num->z);
        mpz_mod(res->z, res->z, mod->z);
        if(options.ZM_cb && mpz_sgn(res->z)==0) {
            PyObject* result;
            if(options.debug)
                fprintf(stderr, "calling %p from %s for %p %p %p %p\n",
                    options.ZM_cb, "divm", res, num, den, mod);
            result = PyObject_CallFunction(options.ZM_cb, "sOOOO",
"divm",
                                           res, num, den, mod);
            if(result != Py_None) {
                Py_DECREF((PyObject*)res);
                return result;
            }
        }
        Py_DECREF((PyObject*)num);
        Py_DECREF((PyObject*)den);
        Py_DECREF((PyObject*)mod);
        return (PyObject*)res;
    } else {
        PyObject* result = 0;
        if(options.ZD_cb) {
            result = PyObject_CallFunction(options.ZD_cb,
                "sOOO", "divm", num, den, mod);
        } else {
            PyErr_SetString(PyExc_ZeroDivisionError, "not invertible");
        }
        Py_DECREF((PyObject*)num);
        Py_DECREF((PyObject*)den);
        Py_DECREF((PyObject*)mod);
        Py_DECREF((PyObject*)res);
        return result;
    }
}


And despite using gmpy 1.01 for six months, I've never
encountered this problem before because in the algorithm
where I use it, my operands are guaranteed to be invertible
by construction, so I have never provoked the corrupting
influence of divm().