[maximus-decimus]

For fun, I ported your 3rd solution to J

  f =: 3 : 0
  'k a' =: y
  n =: #a
  A =: (n, n) $ a
  ^ 0j1 * k * 0.5 ^~ (*: A) + (*: |: A)
  )

I'm curious how the speed compares to Numpy, but I don't have a python environment installed. It ran in 6.6 seconds on my computer for n=10,000, but it used a ridiculous 10 gigabytes of memory. When you consider a 10,000x10,000 matrix has 100M elements, that you need two of them since you add two together and that a short has 4 bytes, you technically should only need 800 MB for the whole thing.