GMPY: divm() memory leak revisited

[mensanator]

Since the following discussion took place (unresolved),

<http://groups.google.c om/group/comp.lang.pytho n/browse_frm/thread/c3bd08ef3e4c478 a/2b54deb522c9b9d 9?lnk=st&q=divm ()+memory+leak+ group:comp.lang .python+author: mensanator&rnum =1&hl=en#2b54de b522c9b9d9>

I've kept it in the back of my mind as I've been
learning to use the base gmp library in c. Now I believe
I know what the problem is.

First, divm() is not a gmp function. It is a derived
function created in gmpy. It is derived from the
gmp invert() function (which I know from my testing
does not leak memory). So it's a gmpy specific bug.

Second, I learned from the gmp c library that temporary
mpz objects must be freed to prevent memory leak. Aha!
The gmpy source code is probably not freeing some temporary
mpz it created.

Third, we have the smoking gun:

divm(a,b,m): returns x such that b*x==a modulo m,
or else raises a ZeroDivisionErr or
exception if no such value x exists
(a, b and m must be mpz objects,
or else get coerced to mpz)

Of course, to "coerce" means to create temporary variables
to pass to the gmp library. It would appear that these
temporary variables are not being freed. Now if I'm right,
then I can prove this by eliminating the need to coerce
the operands by passing mpz's to the divm() function.

#
# if the parameters are already mpz's...
#
z = gmpy.mpz(812875 70543)
x = gmpy.mpz(858993 4592)
y = gmpy.mpz(348678 4401)
tot = 0
while True:
n = input('How many more divm: ')
if n<=0: break
print '%d more...' % n,
#
# ...then they won't need to be coerced
#
for i in xrange(n): gmpy.divm(z,x,y )
tot += n
print '...total %d' % tot
With coercing, I get

C:\Python23\use r\the_full_mont y>python gmpytest.py
How many more divm: 10000000
10000000 more...Fatal Python error: mp_allocate failure
abnormal program termination
peak Commit Charge (K): 792556

Without needing to coerce, the test ran to completion with
flat memory usage.

Unfortunately, c is still somewhat greek to me, but even so,
the problem appears obvious.

static PyObject *
Pygmpy_divm(PyO bject *self, PyObject *args)
{
PympzObject *num, *den, *mod, *res;
if(!PyArg_Parse Tuple(args, "O&O&O&",
Pympz_convert_a rg, &num,
Pympz_convert_a rg, &den,
Pympz_convert_a rg, &mod))
{
return last_try("divm" , 3, 3, args);
}
if(!(res = Pympz_new()))
return NULL;
if(mpz_invert(r es->z, den->z, mod->z)) { /* inverse exists */
mpz_mul(res->z, res->z, num->z);
mpz_mod(res->z, res->z, mod->z);
if(options.ZM_c b && mpz_sgn(res->z)==0) {
PyObject* result;
if(options.debu g)
fprintf(stderr, "calling %p from %s for %p %p %p %p\n",
options.ZM_cb, "divm", res, num, den, mod);
result = PyObject_CallFu nction(options. ZM_cb, "sOOOO",
"divm",
res, num, den, mod);
if(result != Py_None) {
Py_DECREF((PyOb ject*)res);
return result;
}
}
return (PyObject*)res;
} else {
PyObject* result = 0;
if(options.ZD_c b) {
result = PyObject_CallFu nction(options. ZD_cb,
"sOOO", "divm", num, den, mod);
} else {
PyErr_SetString (PyExc_ZeroDivi sionError, "not invertible");
}
Py_DECREF((PyOb ject*)res);
return result;
}
}

Note that 4 PympzObjects get created but only res gets passed to
Py_DECREF (which seems to be the method by which it's freed, not
100% sure about this). But I notice that other functions that
coerce variables call Py_DECREF on each of the coerced variables:

static PyObject *
Pygmpy_gcd(PyOb ject *self, PyObject *args)
{
PympzObject *a, *b, *c;
TWO_ARG_CONVERT ED("gcd", Pympz_convert_a rg,&a,&b);
assert(Pympz_Ch eck((PyObject*) a));
assert(Pympz_Ch eck((PyObject*) b));
if(!(c = Pympz_new())) {
Py_DECREF((PyOb ject*)a); Py_DECREF((PyOb ject*)b);
return NULL;
}
mpz_gcd(c->z, a->z, b->z);
Py_DECREF((PyOb ject*)a); Py_DECREF((PyOb ject*)b);
return (PyObject*)c;
}

Here the PympzObject c is not freed because it is returned
as the result of the function, but the coerced variables a and b
are.

So the fix may be to simply add

Py_DECREF((PyOb ject*)num);
Py_DECREF((PyOb ject*)den);
Py_DECREF((PyOb ject*)mod);

to divm(). (Assuming it's that simple, I could have overlooked
something.)

Unfortunately, I don't have any means of testing this theory.

I did see the reference to

<http://www.vrplumber.c om/programming/mstoolkit>

which I will be trying eventually, but in case I can't get that
to work I wanted to have this posted in case someone else wants
to take a crack at it.

[mensanator]

me********@aol. com wrote:

Since the following discussion took place (unresolved),

<http://groups.google.c om/group/comp.lang.pytho n/browse_frm/thread/c3bd08ef3e4c478 a/2b54deb522c9b9d 9?lnk=st&q=divm ()+memory+leak+ group:comp.lang .python+author: mensanator&rnum =1&hl=en#2b54de b522c9b9d9>

I've kept it in the back of my mind as I've been
learning to use the base gmp library in c. Now I believe
I know what the problem is.

First, divm() is not a gmp function. It is a derived
function created in gmpy. It is derived from the
gmp invert() function (which I know from my testing
does not leak memory). So it's a gmpy specific bug.

Second, I learned from the gmp c library that temporary
mpz objects must be freed to prevent memory leak. Aha!
The gmpy source code is probably not freeing some temporary
mpz it created.

Third, we have the smoking gun:

divm(a,b,m): returns x such that b*x==a modulo m,
or else raises a ZeroDivisionErr or
exception if no such value x exists
(a, b and m must be mpz objects,
or else get coerced to mpz)

Of course, to "coerce" means to create temporary variables
to pass to the gmp library. It would appear that these
temporary variables are not being freed. Now if I'm right,
then I can prove this by eliminating the need to coerce
the operands by passing mpz's to the divm() function.

#
# if the parameters are already mpz's...
#
z = gmpy.mpz(812875 70543)
x = gmpy.mpz(858993 4592)
y = gmpy.mpz(348678 4401)
tot = 0
while True:
n = input('How many more divm: ')
if n<=0: break
print '%d more...' % n,
#
# ...then they won't need to be coerced
#
for i in xrange(n): gmpy.divm(z,x,y )
tot += n
print '...total %d' % tot
With coercing, I get

C:\Python23\use r\the_full_mont y>python gmpytest.py
How many more divm: 10000000
10000000 more...Fatal Python error: mp_allocate failure
abnormal program termination
peak Commit Charge (K): 792556

Without needing to coerce, the test ran to completion with
flat memory usage.

Unfortunately, c is still somewhat greek to me, but even so,
the problem appears obvious.

static PyObject *
Pygmpy_divm(PyO bject *self, PyObject *args)
{
PympzObject *num, *den, *mod, *res;
if(!PyArg_Parse Tuple(args, "O&O&O&",
Pympz_convert_a rg, &num,
Pympz_convert_a rg, &den,
Pympz_convert_a rg, &mod))
{
return last_try("divm" , 3, 3, args);
}
if(!(res = Pympz_new()))
return NULL;
if(mpz_invert(r es->z, den->z, mod->z)) { /* inverse exists */
mpz_mul(res->z, res->z, num->z);
mpz_mod(res->z, res->z, mod->z);
if(options.ZM_c b && mpz_sgn(res->z)==0) {
PyObject* result;
if(options.debu g)
fprintf(stderr, "calling %p from %s for %p %p %p %p\n",
options.ZM_cb, "divm", res, num, den, mod);
result = PyObject_CallFu nction(options. ZM_cb, "sOOOO",
"divm",
res, num, den, mod);
if(result != Py_None) {
Py_DECREF((PyOb ject*)res);
return result;
}
}
return (PyObject*)res;
} else {
PyObject* result = 0;
if(options.ZD_c b) {
result = PyObject_CallFu nction(options. ZD_cb,
"sOOO", "divm", num, den, mod);
} else {
PyErr_SetString (PyExc_ZeroDivi sionError, "not invertible");
}
Py_DECREF((PyOb ject*)res);
return result;
}
}

Note that 4 PympzObjects get created but only res gets passed to
Py_DECREF (which seems to be the method by which it's freed, not
100% sure about this). But I notice that other functions that
coerce variables call Py_DECREF on each of the coerced variables:

static PyObject *
Pygmpy_gcd(PyOb ject *self, PyObject *args)
{
PympzObject *a, *b, *c;
TWO_ARG_CONVERT ED("gcd", Pympz_convert_a rg,&a,&b);
assert(Pympz_Ch eck((PyObject*) a));
assert(Pympz_Ch eck((PyObject*) b));
if(!(c = Pympz_new())) {
Py_DECREF((PyOb ject*)a); Py_DECREF((PyOb ject*)b);
return NULL;
}
mpz_gcd(c->z, a->z, b->z);
Py_DECREF((PyOb ject*)a); Py_DECREF((PyOb ject*)b);
return (PyObject*)c;
}

Here the PympzObject c is not freed because it is returned
as the result of the function, but the coerced variables a and b
are.

So the fix may be to simply add

Py_DECREF((PyOb ject*)num);
Py_DECREF((PyOb ject*)den);
Py_DECREF((PyOb ject*)mod);

to divm(). (Assuming it's that simple, I could have overlooked
something.)

Unfortunately, I don't have any means of testing this theory.

I did see the reference to

<http://www.vrplumber.c om/programming/mstoolkit>

which I will be trying eventually, but in case I can't get that
to work I wanted to have this posted in case someone else wants
to take a crack at it.

I completely forgot to mention the fourth thing I discovered.

The linear congruence algorithm used in divm() is wrong. To wit:

divm(6,12,14)
Traceback (most recent call last):
File "<pyshell#4 >", line 1, in -toplevel-
divm(6,12,14)
ZeroDivisionErr or: not invertible

Sure, (12,14) is not invertible, but that is not a requirement for
solving a linear congruence. All that's required is that gcd(12,14)
divides 6, which it obviously does. The divm() function cannot
handle the case when b and m are not coprime. A properly
written linear congruence algorithm will work around the fact
that (12,14) is not invertible. To wit:
linear_congruen ce(12,14,6)

mpz(4)

The correct algorithm can't solve a problem that's not solvable,
but it shouldn't let non-invertability cause an exception. Here's
the correct algorithm which divm() ought to incorprorate.

from gmpy import *
def linear_congruen ce(x,y,z):
#
# xa == z (mod y)
#
g = gcd(x,y)
d = divmod(z,g)
if d[1]==0:
#
# gcd(x,y) divides z, solution exists
#
if g==1:
#
# x,y coprime, modular inverse exists
#
a = invert(x,y)*z % y
else:
#
# x,y not coprime, no modular inverse
# ...but wait, if we get here g divides z and also
# by definition, x & y, so divide x,y,z by g to create
# a new congruence with x,y now coprime and invert() valid
#
x = x/g
y = y/g
z = z/g
a = invert(x,y)*z % y
else:
#
# g doesn't divide z, no solution
#
a = -1
return a

[Alex Martelli]

me********@aol. com <me********@aol .com> wrote:
...

Unfortunately, I don't have any means of testing this theory.

Yep -- I reproduced the memory leak you mentioned, and easily fixed it
(exactly as you suggest) in the current CVS version of gmpy (meant to be
"1.01 release candidate"). I need to fix some other pending bugs, then
remind myself of how the (expl.del) one makes a release on Sourceforge,
then it shd be fine (except for Windows -- I have no Windows development
system around to do anything... Mac and Linux only).

One of the bugs I must still get around to examining is the algorithmic
one you mention in your next post, btw.

Thanks for your diagnostic and debugging help. BTW, should you wish to
mail me privately, I'm "aleaxit" (my favourite userid, with or w/o the
trailing "it" depending on service -- see www.aleax.com to understand
why that is the case), and the best way to reach me these days is
through "gmail.com" .
Alex

[mensanator]

Alex Martelli wrote:

me********@aol. com <me********@aol .com> wrote:
...

Unfortunately, I don't have any means of testing this theory.
Yep -- I reproduced the memory leak you mentioned, and easily fixed it
(exactly as you suggest) in the current CVS version of gmpy (meant to be
"1.01 release candidate"). I need to fix some other pending bugs, then
remind myself of how the (expl.del) one makes a release on Sourceforge,
then it shd be fine

That's great!
(except for Windows -- I have no Windows development
system around to do anything... Mac and Linux only).
Well, I hope someone makes a Windows binary once you publish ver 1.01.
I downloaded all that MSVC command line stuff (over 800 MB!) but am not
looking forward to trying to get it to work.

One of the bugs I must still get around to examining is the algorithmic
one you mention in your next post, btw.

Thanks for your diagnostic and debugging help. BTW, should you wish to
mail me privately, I'm "aleaxit" (my favourite userid, with or w/o the
trailing "it" depending on service -- see www.aleax.com to understand
why that is the case), and the best way to reach me these days is
through "gmail.com" .
Alex