def chop(x, tol=1.e-14):
  r = np.real(x)
  if abs(r)<tol:
    r = 0

  i = np.imag(x)
  if abs(i)<tol:
    i = 0

  if i == 0:
    n = r
  elif r == 0:
    n = i*1j
  else:
    n = r + i*1j
  return(n)

import numpy as np
a = 1e-16 + 1e-16*1j
b = 1e-16 + 1j
c = 1 + 1e-16*1j
chop(a), chop(b), chop(c)

n = 4
w4 = np.exp(2*np.pi*1j/n)
w4 = chop(w4)
w4

F4 = np.zeros([4,4], dtype=complex)
for i in range(4):
    for j in range(4):
        F4[i,j] = w4**(i*j)
F4

I2 = np.array([[1,0],[0,1]]); I2

D2_4 = np.zeros([2,2], dtype=complex)
for i in range(2):
  D2_4[i,i] = w4**(i)
D2_4

M1 = np.concatenate([I2,D2_4],axis=1)
M2 = np.concatenate([I2,-D2_4],axis=1)
M3 = np.concatenate([M1,M2],axis=0); M3

n = 2
w2 = np.exp(2*np.pi*1j/n)
w2 = chop(w2)
w2

F2 = np.zeros([2,2], dtype=complex)
for i in range(2):
    for j in range(2):
        F2[i,j] = w2**(i*j)
F2

Z2 = np.zeros([2,2], dtype=complex)
Z2

K1 = np.concatenate([F2,Z2],axis=1)
K2 = np.concatenate([Z2,F2],axis=1)
K3 = np.concatenate([K1,K2],axis=0); K3

P4 = np.array([[1,0,0,0],
               [0,0,1,0],
               [0,1,0,0],
               [0,0,0,1]])
P4

M3.dot(K3).dot(P4)

F4

np.all(F4 == M3.dot(K3).dot(P4))