module confict? gmpy and operator

[mensanator]

## Holy Mother of Pearl!
##
## >>> for i in range(10):
## for j in range(10):
## print '%4d' % (gmpy.mpz(i)*gmpy.mpz(j)),
## print
##
##
## 0 0 0 0 0 0 0 0 0 0
## 0 1 2 3 1 5 6 7 2 9
## 0 2 4 6 2 10 12 14 4 18
## 0 3 6 9 3 15 18 21 6 27
## 0 1 2 3 1 5 6 7 2 9
## 0 5 10 15 5 25 30 35 10 45
## 0 6 12 18 6 30 36 42 12 54
## 0 7 14 21 7 35 42 49 14 63
## 0 2 4 6 2 10 12 14 4 18
## 0 9 18 27 9 45 54 63 18 81

No wonder I couldn't figure out why my program wasn't working,
gmpy mutiplication seems to have stopped working.

My guess is that when I did

import gmpy
import random
import operator

I somehow messed up gmpy's mutilpication. I brought in operator
because I needed to do

NN = sum(map(operator.__mul__,X,C))

which seemed to work ok. It was later when I got

8 * 3 = 6 (the 3 was an mpz)

that I got in trouble. So I got rid of operator and did

def product(a,b): return a*b

NN = sum(map(product,X,C))

which seems to work ok.

Am I correct that there is a conflict with

import gmpy
import operator

Would reversing the import order help or should I just avoid
mixing them?

 

[Tim Peters]

[[email protected]]

## Holy Mother of Pearl!
##
## >>> for i in range(10):
## for j in range(10):
## print '%4d' % (gmpy.mpz(i)*gmpy.mpz(j)),
## print
##
##
## 0 0 0 0 0 0 0 0 0 0
## 0 1 2 3 1 5 6 7 2 9
## 0 2 4 6 2 10 12 14 4 18
## 0 3 6 9 3 15 18 21 6 27
## 0 1 2 3 1 5 6 7 2 9
## 0 5 10 15 5 25 30 35 10 45
## 0 6 12 18 6 30 36 42 12 54
## 0 7 14 21 7 35 42 49 14 63
## 0 2 4 6 2 10 12 14 4 18
## 0 9 18 27 9 45 54 63 18 81

No wonder I couldn't figure out why my program wasn't working,
gmpy mutiplication seems to have stopped working.

My guess is that when I did

import gmpy
import random
import operator

I somehow messed up gmpy's mutilpication. I brought in operator
because I needed to do

NN = sum(map(operator.__mul__,X,C))

which seemed to work ok. It was later when I got

8 * 3 = 6 (the 3 was an mpz)

that I got in trouble. So I got rid of operator and did

def product(a,b): return a*b

NN = sum(map(product,X,C))

which seems to work ok.

Please try to give a complete program that illustrates the problem.
For example, this program shows no problem on my box (Windows Python
2.4.3):

"""
import gmpy
import random
import operator

NN = sum(map(operator.__mul__, [1, 2, 3], [4, 5, 6]))

for i in range(10):
for j in range(10):
print '%4d' % (gmpy.mpz(i)*gmpy.mpz(j)),
print
"""

The output was:

0 0 0 0 0 0 0 0 0 0
0 1 2 3 4 5 6 7 8 9
0 2 4 6 8 10 12 14 16 18
0 3 6 9 12 15 18 21 24 27
0 4 8 12 16 20 24 28 32 36
0 5 10 15 20 25 30 35 40 45
0 6 12 18 24 30 36 42 48 54
0 7 14 21 28 35 42 49 56 63
0 8 16 24 32 40 48 56 64 72
0 9 18 27 36 45 54 63 72 81

Am I correct that there is a conflict with

import gmpy
import operator

I don't see a problem of any kind, and there's no way to guess from
the above what you did that mattered. It's hard to believe that
merely importing any module (operator or otherwise) has any effect on
gmpy.

Would reversing the import order help or should I just avoid
mixing them?

Both are dubious.

 

[mensanator]

[[email protected]]

Please try to give a complete program that illustrates the problem.
For example, this program shows no problem on my box (Windows Python
2.4.3):

Sorry about that. I thought the problem was more obvious.

I don't see a problem of any kind, and there's no way to guess from
the above what you did that mattered. It's hard to believe that
merely importing any module (operator or otherwise) has any effect on
gmpy.

That's what I was asking about. I don't know enough about low level
stuff to know. I've been using gmpy for years with the only problem
being the memory leak that was fixed in the latest version. I've never
used operator before, so I just assumed that's where the problem was.

But I don't see why the program behaves the way it does either.

Here's the complete program:

import gmpy
import random
import operator

def product(a,b):
return a*b

def gcdlist(X):
mingcd = 999999
for i in xrange(1,len(X)):
g = gmpy.gcd(X[i-1],X)
if g<mingcd: mingcd = g
return mingcd

X = [8,16,20,24,40,72,84,92]
g = gcdlist(X)
s = sum(X)

print ' X:',X
print ' gcd:',g

N = 0
while N<s:
N = g * random.randint(s,3333)
print ' N:',N

if N>s:
C = [1 for i in X]
diff = N-s
done = False
i = -1
XL = -len(X)
while not done:
while diff>=X:
C += 1
diff -= X
if diff==0:
done = True
else:
i -= 1
if i<XL:
done = True
NN = sum(map(operator.__mul__,X,C))
print ' X:',X
print ' C:',C
print ' NN:',NN
print ' diff:',diff
print
if diff>0:
p = 0
q = 1
done = False
while not done:
gpq = gmpy.gcd(X[p],X[q])
if diff % gpq == 0:
done = True
else:
q += 1
if q==len(X):
p += 1
q = p + 1
print 'p: %d q: %d' % (p,q)
a = gmpy.divm(diff,X[p],X[q])
b = (X[p]*a - diff)/X[q]
print 'a: %d b: %d X[p]: %d X[q]: %d' %
(a,b,X[p],X[q])
if C[q]==b:
print 'non-zero adjustment'
print 'a: %d b: %d' % (a + X[q],b + X[p])
C[p] += a + X[q]
C[q] -= b + X[p]
else:
C[p] += a
C[q] -= b
NN = sum(map(operator.__mul__,X,C))
print ' X:',X
print ' C:',C
print ' NN:',NN
print

This is what the output should look like

X: [8, 16, 20, 24, 40, 72, 84, 92]
gcd: 4
N: 4492
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [1, 1, 1, 1, 1, 1, 2, 45]
NN: 4488
diff: 4

p: 0 q: 2
a: 3 b: 1 X[p]: 8 X[q]: 20
non-zero adjustment
a: 23 b: 9
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [24, 1, -8, 1, 1, 1, 2, 45]
NN: 4492

X is a list of integers selected so that the GCD of the entire list
is >2.
N is selected to be greater than sum(X) and divisible by GCD(X).
C is initialized to [1,1,1,1,1,1,1,1].
NN is the sum of each X element multiplied by its corresponding C
element. The values of C are then adjusted until the difference
between N and NN is 0.

If a difference of 0 is unobtainable (in the example diff=4 is smaller
than the smallest element of X), then the program solves a linear
congruence by first finding a pair of X whose GCD divides diff
(8 and 20). For the solution a=3, b=1, we add three to C[0] and
subtract 1 from C[2] making
C: [1, 1, 1, 1, 1, 1, 2, 45]
into
C: [4, 1, 0, 1, 1, 1, 2, 45]

But I don't want any 0s in the list, only positive and negative
integers. Thus, in this example, a 0 adjustment is made:
C: [24, 1, -8, 1, 1, 1, 2, 45]

Ok, now the problem. In the first three runs, diff was 0, so the
linear congruence was skipped. But in the fourth run, b=0
which indicates the error (8*3 evaluated to 6 resulting in b=0).


Note I tested gmpy after each run from the Idle prompt.
Once the failure occured, the problem continues to exist.
That's how I made that original demo, I wrote a simple test
AFTER gmpy had gotten funny.
X: [8, 16, 20, 24, 40, 72, 84, 92]
gcd: 4
N: 12192
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [1, 1, 2, 1, 2, 1, 1, 129]
NN: 12192
diff: 0
X: [8, 16, 20, 24, 40, 72, 84, 92]
gcd: 4
N: 7340
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [1, 1, 1, 1, 1, 1, 2, 76]
NN: 7340
diff: 0
X: [8, 16, 20, 24, 40, 72, 84, 92]
gcd: 4
N: 8072
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [2, 1, 1, 1, 1, 2, 1, 84]
NN: 8072
diff: 0
X: [8, 16, 20, 24, 40, 72, 84, 92]
gcd: 4
N: 6296
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [2, 1, 1, 1, 2, 1, 1, 65]
NN: 6292
diff: 4

p: 0 q: 2
a: 3 b: 0 X[p]: 8 X[q]: 20
X: [8, 16, 20, 24, 40, 72, 84, 92]
C: [mpz(5), 1, mpz(1), 1, 2, 1, 1, 65]
NN: 6286
6


Hey, maybe the gmpy divm() function isn't fixed after all:

IDLE 1.1a36

Can you verify this?

 

[mensanator]

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.1a3mpz(3)

Correct answer...
(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.
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...
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(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).

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().