1

I am trying to find a way for modular matrix inversion. I found the code here:

def generalizedEuclidianAlgorithm(a, b):
    if b > a:
        #print a, b
        return generalizedEuclidianAlgorithm(b,a);
    elif b == 0:
        return (1, 0);
    else:
        #print a,b
        (x, y) = generalizedEuclidianAlgorithm(b, a % b);
        return (y, x - (a / b) * y)

def inversemodp(a, p):
    a = a % p
    if (a == 0):
        print "a is 0 mod p"
        return 0
    (x,y) = generalizedEuclidianAlgorithm(p, a % p);
    return y % p

def identitymatrix(n):
    return [[long(x == y) for x in range(0, n)] for y in range(0, n)]

def inversematrix(matrix, q):
    n = len(matrix)
    A = np.matrix([[ matrix[j, i] for i in range(0,n)] for j in range(0, n)], dtype = long)
    Ainv = np.matrix(identitymatrix(n), dtype = long)
    for i in range(0, n):
        factor = inversemodp(A[i,i], q)
        A[i] = A[i] * factor % q
        Ainv[i] = Ainv[i] * factor % q
        for j in range(0, n):
            if (i != j):
                factor = A[j, i]
                A[j] = (A[j] - factor * A[i]) % q
                Ainv[j] = (Ainv[j] - factor * Ainv[i]) % q
                # print A, Ainv
                # print i, j, factor
    return Ainv

When I test with small prime q and matrix elements, the result is correct. However, when I test with large prime q and matrix consisting of large elements (e.g, 1024 bits), it printed out the error:

A = np.matrix([[ matrix[j, i] for i in range(0,n)] for j in range(0, n)], dtype = long)

File "C:\Python27\lib\site-packages\numpy\matrixlib\defmatrix.py", line 257, in __new__

arr = N.array(data, dtype=dtype, copy=copy)

Is it because of dtype of np.matrix? If so, why "long" datatype cannot support this case?

Could you please show me how to solve the error. Thanks in advance

Community
  • 1
  • 1
santa
  • 193
  • 10
  • `dtype=long` is turned into a fixed-width integer dtype (e.g. `int64` on my machine). You probably want `dtype=object` instead. – Mark Dickinson Jul 11 '14 at 17:32
  • Is there a reason for doing it by hand? Use the numpy library inside the [scipy](http://www.scipy.org/) package. A lot of what you're trying to do has already been implemented and optimized there. – Jesuisme Jul 11 '14 at 16:35
  • 1
    @Jesuisme I don't think numpy support modular matrix inversion – santa Jul 12 '14 at 03:19
  • @MarkDickinson: thank you. When I change to dtype=object, it results errors like 'cannot support *: long and matrix' at: -A[i] = A[i] * factor % q, -Ainv[i] = Ainv[i] * factor % q, -A[j] = (A[j] - factor * A[i]) % q, -Ainv[j] = (Ainv[j] - factor * Ainv[i]) % q – santa Jul 12 '14 at 03:34
  • @MarkDickinson: sorry for the messy comment because I don't know how to write a code in the comment – santa Jul 12 '14 at 03:37
  • 1
    Interesting. It looks like this is a bug in NumPy. I can reproduce your problem with NumPy 1.6.2, but not with NumPy 1.8.1. See https://github.com/numpy/numpy/issues/2469. A workaround is to use `numpy.array` instead of `numpy.matrix`, which means you need to use `numpy.dot` to do a matrix multiplication. There are many reasons to avoid the `matrix` type in NumPy; looks like you've just found another one. :-) – Mark Dickinson Jul 12 '14 at 08:02

1 Answers1

2

I have just figured out that np.array causes many problems. One solution is just to treat the matrix as 2-dimensional list. The following is my code

def generalizedEuclidianAlgorithm(a, b):
    if b > a:
        return generalizedEuclidianAlgorithm(b,a);
    elif b == 0:
        return (1, 0);
    else:
        (x, y) = generalizedEuclidianAlgorithm(b, a % b);
        return (y, x - (a / b) * y)

def inversemodp(a, p):
    a = a % p
    if (a == 0):
        print "a is 0 mod p"
        return 0
    (x,y) = generalizedEuclidianAlgorithm(p, a % p);
    return y % p

def identitymatrix(n):
    return [[long(x == y) for x in range(0, n)] for y in range(0, n)

def multiply_vector_scalar (vector, scalar, q):
    kq = []
    for i in range (0, len(vector)):
        kq.append (vector[i] * scalar %q)
    return kq

def minus_vector_scalar1(vector1, scalar, vector2, q):
    kq = []
    for i in range (0, len(vector1)):
        kq.append ((vector1[i] - scalar * vector2[i]) %q)
    return kq

def inversematrix1(matrix, q):
    n = len(matrix)

    A =[]
    for j in range (0, n):
        temp = []
        for i in range (0, n):
            temp.append (matrix[j][i])
        A.append(temp)

    Ainv = identitymatrix(n)

    for i in range(0, n):
        factor = inversemodp(A[i][i], q)
        A[i] = multiply_vector_scalar(A[i],factor,q)
        Ainv[i] = multiply_vector_scalar(Ainv[i],factor,q)
        for j in range(0, n):
            if (i != j):
                factor = A[j][i]
                A[j] = minus_vector_scalar1(A[j],factor,A[i],q)
                Ainv[j] = minus_vector_scalar1(Ainv[j],factor,Ainv[i],q)
    return Ainv
santa
  • 193
  • 10