

def dec2den(j,N,d):
    den = [0 for k in range(0,N)]
    jv = j
    for k in range(0,N):
        if jv >= d**(N-1-k):
            den[k] = jv//(d**(N-1-k))
            jv = jv - den[k]*d**(N-1-k)
    return den

def den2dec(local,d):
    # local = list with the local computational base state values 
    # d = individual qudit dimension
    N = len(local)
    j = 0
    for k in range(0,N):
        j += local[k]*d**(N-1-k)
    return j # value of the global computational basis index

def dec2den_opt(j,N,d): # optimized version
    den = [0 for k in range(0,N)]
    if j < d:
        den[N-1] = j; return den
    rem = [] # list for the remainder
    k = j
    while k >= d:
        rem.append(k%d)
        k = k//d
    rem.append(k%d)
    for j in range(0,len(rem)):
        den[N-1-j] = rem[j]
    return den

import numpy as np

def swap_qudits(d,N,t,c):
    # d = qudit dimension; N = No. of qudits; t = target qudit (0,...,N-1); c = control qudit (0,...,N-1)
    D = d**N; swap = np.zeros((D,D), dtype=int) # set all d^2N elements to 0
    for k in range (0,d**N):
        localr = dec2den_opt(k,N,d) # sequence of values for the local basis
        #locall = [0 for k in range(0,N)]
        #for s in range(0,N):
        #    locall[s] = localr[s]
        locall = localr.copy()
        #rt = localr[t]; 
        locall[c] = localr[t]
        #rc = localr[c]; 
        locall[t] = localr[c]
        globall = den2dec(locall,d)
        swap[globall,k] = 1
    return swap

d=2; N=2; t=0; c=1; print(swap_qudits(d,N,t,c)) # ok

d=3; N=2; t=0; c=1; print(swap_qudits(d,N,t,c)) # ok

d=2; N=3; t=0; c=2; print(swap_qudits(d,N,t,c)) # ok

