Find a (partial) list of matrix unity roots-CodePudding

I am looking for a way to build a list of matrix square roots of 1 and -1. Only interested in dimensions 2 to 4, and I don't even need a complete set, just a practical usable one (meaning I could stop the loop at some point). Here's my (very) inefficient solution

def racines(n, dim, field='real'):
    if field == 'real':
        unites = [1, 0, -1]
    elif field == 'complex':
        unites = [1, 0, -1, 1j, -1j]
    if dim == 2:
        liste = [np.array([[a,b], [c,d]]) for a in unites
                 for b in unites for c in unites 
                 for d in unites]
    elif dim == 3:
        liste = [np.array([[a, b, c], [d, e, f], [g, h, i]])
                 for a in unites for b in unites for c in unites 
                 for d in unites for e in unites for f in unites
                 for g in unites for h in unites for i in unites]
    elif dim == 4:
        liste = [np.array([[a,b,c,d], [e,f,g,h], [i,j,k,l], [m,o,p,q]])
                 for a in unites for b in unites for c in unites for d in unites
                 for e in unites for f in unites for g in unites for h in unites
                 for i in unites for j in unites for k in unites for l in unites
                 for m in unites for o in unites for p in unites for q in unites]
    liste = [m for m in liste if np.linalg.det(m) != 0]
    liste = [m for m in liste if np.array_equal(np.matmul(m, m), np.identity(dim, dtype='int')*n)]
return(liste)

Works well until 4d where of course it becomes too much to handle... How to make it faster ?

CodePudding user response：

You can do it recursively. For exemple if I want the lists of n numbers between 0 and 4 you could use

def get_lists(size, max_it, accumulator = [], lists = []):
    if size ==0:
        lists.append(accumulator)
    else:
        for i in range(5):
            if len(lists) == max_it:
                break
            l = accumulator.copy()
            l.append(i)
            get_lists(size-1,max_it, l, lists)

            
def lists(size,max_it=-1):
    l = []  
    get_lists(size,max_it, [], l)
    return l        
print(lists(3,10))

with max_it the amount of different list you want (set to -1 if you want them all)

CodePudding user response：

For those who would need it, here's what did the trick. Still not pretty, but it works a lot faster in 4d

def racines(n, dim, field='real'):
    unites = [1, 0, -1]
    if field == 'complex':
        unites  = [1j, -1j]
    elif field != 'real':
        print('Unknown field !')
        return([])
    if dim == 2:
        liste = [np.array([[a,b], [c,d]]) for a in unites
                 for b in unites for c in unites 
                 for d in unites]
    elif dim == 3:
        liste = [np.array([[a, b, c], [d, e, f], [g, h, i]])
                 for a in unites for b in unites for c in unites 
                 for d in unites for e in unites for f in unites
                 for g in unites for h in unites for i in unites]
    elif dim == 4:
        roots_2d = racines(n, 2, field=field)
        neg_roots_2d = racines(-n, 2, field=field)
        other = [np.identity(2, dtype='int'), np.identity(2, dtype='int')*-1, np.zeros([2,2], dtype='int')]
        unites = roots_2d   neg_roots_2d   other
        liste = [np.hstack([np.vstack([a,b]), np.vstack([c,d])])
                 for a in unites for b in unites for c in unites for d in unites]
    else:
        print('No can do !')
        return([])
    liste = [m for m in liste if np.linalg.det(m) != 0]
    liste = [m for m in liste if np.array_equal(np.matmul(m, m), np.identity(dim, dtype='int')*n)]
    return(liste)
# res = racines(1, 4, field='real') # 2.43s
# res = racines(1, 4, field='complex') # 2min 18s