diff --git a/README.md b/README.md
index 1895fa7..76e4a18 100644
--- a/README.md
+++ b/README.md
@@ -27,9 +27,11 @@ TFinder: Similar to Tester, but only print the circuit's matrix.
 AllInit: Similar to Tester, but run for each possible initial state (2^N of them).
 
 ## Commands
+i, j, k are indicies of qubit. (0-indexed)
 
 - Not i
 - Hadamard i
+    + You can also use "H" instead of "Hadamard"
 - CNot i j
 - Swap i j
 - CCNot i j k
@@ -40,10 +42,10 @@ AllInit: Similar to Tester, but run for each possible initial state (2^N of them
 - TDag i
 - SqrtNot i
 - SqrtNotDag i
-- SqrtSwap i j (Not implemented yet)
+- SqrtSwap i j
+    + Not implemented yet
 - Rot i deg
     + Rotate qubit counterclockwise by degree, not rad 
 - Measure n
-    + Measures the joint qubit state n times using standard basis
-    + Only works when using Tester
-    + Will end command read
\ No newline at end of file
+    + Measures the joint qubit state n times using the standard basis
+    + Only works when using Tester
\ No newline at end of file
diff --git a/build.gradle.kts b/build.gradle.kts
index 602f8ca..4948533 100644
--- a/build.gradle.kts
+++ b/build.gradle.kts
@@ -6,7 +6,7 @@ plugins {
     id("com.github.johnrengelman.shadow") version "6.0.0"
 }
 group = "com.dedztbh"
-version = "1.0"
+version = "1.1"
 
 val ejmlVersion = "0.39"
 
@@ -17,6 +17,7 @@ repositories {
 dependencies {
     implementation("org.ejml:ejml-core:${ejmlVersion}")
     implementation("org.ejml:ejml-zdense:${ejmlVersion}")
+    testImplementation("org.ejml:ejml-cdense:${ejmlVersion}")
 }
 
 tasks.withType<KotlinCompile> {
diff --git a/src/main/kotlin/main.kt b/src/main/kotlin/main.kt
index d6af037..4996f68 100644
--- a/src/main/kotlin/main.kt
+++ b/src/main/kotlin/main.kt
@@ -1,6 +1,4 @@
 import operator.Operator
-import util.read
-import util.reader
 import java.io.File
 
 /**
diff --git a/src/main/kotlin/matrix/const.kt b/src/main/kotlin/matrix/const.kt
new file mode 100644
index 0000000..175243d
--- /dev/null
+++ b/src/main/kotlin/matrix/const.kt
@@ -0,0 +1,61 @@
+package matrix
+
+/**
+ * Created by DEDZTBH on 2020/09/25.
+ * Project KuantumCircuitSim
+ */
+
+/** 2^(-1/2) */
+const val HALF_AMPL = 0.70710678118654757273731092936941422522068023681640625
+
+/** Don't change these constant matrices! */
+val NOT = CMatrix(
+    arrayOf(
+        doubleArrayOf(0.0, 0.0, 1.0, 0.0),
+        doubleArrayOf(1.0, 0.0, 0.0, 0.0)
+    )
+)
+val H = CMatrix(
+    arrayOf(
+        doubleArrayOf(HALF_AMPL, 0.0, HALF_AMPL, 0.0),
+        doubleArrayOf(HALF_AMPL, 0.0, -HALF_AMPL, 0.0)
+    )
+)
+val I1 = COps.identity(1)
+val I2 = COps.identity(2)
+
+val KET0 = CMatrix(arrayOf(doubleArrayOf(1.0, 0.0), doubleArrayOf(0.0, 0.0)))
+val KET1 = CMatrix(arrayOf(doubleArrayOf(0.0, 0.0), doubleArrayOf(1.0, 0.0)))
+
+val KETBRA0 = COps.diag(1.0, 0.0, 0.0, 0.0)
+val KETBRA1 = COps.diag(0.0, 0.0, 1.0, 0.0)
+val SQRT_NOT = CMatrix(
+    arrayOf(
+        doubleArrayOf(0.5, 0.5, 0.5, -0.5),
+        doubleArrayOf(0.5, -0.5, 0.5, 0.5),
+    )
+)
+val SQRT_NOT_DAG = CMatrix(2, 2).also {
+    COps.transposeConjugate(SQRT_NOT, it)
+}
+val Z = CMatrix(
+    arrayOf(
+        doubleArrayOf(1.0, 0.0, 0.0, 0.0),
+        doubleArrayOf(0.0, 0.0, -1.0, 0.0),
+    )
+)
+val S = CMatrix(
+    arrayOf(
+        doubleArrayOf(1.0, 0.0, 0.0, 0.0),
+        doubleArrayOf(0.0, 0.0, 0.0, 1.0),
+    )
+)
+val T = CMatrix(
+    arrayOf(
+        doubleArrayOf(1.0, 0.0, 0.0, 0.0),
+        doubleArrayOf(0.0, 0.0, HALF_AMPL, HALF_AMPL),
+    )
+)
+val TDag = CMatrix(2, 2).also {
+    COps.transposeConjugate(T, it)
+}
diff --git a/src/main/kotlin/util/matrix.kt b/src/main/kotlin/matrix/op.kt
similarity index 53%
rename from src/main/kotlin/util/matrix.kt
rename to src/main/kotlin/matrix/op.kt
index 8f2c325..3758fea 100644
--- a/src/main/kotlin/util/matrix.kt
+++ b/src/main/kotlin/matrix/op.kt
@@ -1,11 +1,8 @@
-package util
+package matrix
 
 import org.ejml.UtilEjml
-import org.ejml.data.Complex_F64
 import org.ejml.data.MatrixSparse
 import org.ejml.data.MatrixType
-import org.ejml.data.ZMatrixRMaj
-import org.ejml.dense.row.CommonOps_ZDRM
 import org.ejml.ops.MatrixIO
 import java.io.PrintStream
 import java.text.DecimalFormat
@@ -17,16 +14,16 @@ import java.text.DecimalFormat
  */
 
 /**
- * Kronecker product of two ZDRM matrices
+ * Kronecker product of two complex matrices
  */
-infix fun ZMatrixRMaj.kron(B: ZMatrixRMaj) = let { A ->
+infix fun CMatrix.kron(B: CMatrix) = let { A ->
     val numColsC = A.numCols * B.numCols
     val numRowsC = A.numRows * B.numRows
 
-    val C = ZMatrixRMaj(numRowsC, numColsC)
+    val C = CMatrix(numRowsC, numColsC)
 
-    val acomplex = Complex_F64()
-    val bcomplex = Complex_F64()
+    val acomplex = CNumber()
+    val bcomplex = CNumber()
     for (i in 0 until A.numRows) {
         for (j in 0 until A.numCols) {
             A.get(i, j, acomplex)
@@ -46,77 +43,16 @@ infix fun ZMatrixRMaj.kron(B: ZMatrixRMaj) = let { A ->
 /**
  * Convenient pure ops
  */
-operator fun ZMatrixRMaj.plus(B: ZMatrixRMaj) =
-    ZMatrixRMaj(numRows, numCols).also { CommonOps_ZDRM.add(this, B, it) }
+operator fun CMatrix.plus(B: CMatrix) =
+    CMatrix(numRows, numCols).also { COps.add(this, B, it) }
 
-operator fun ZMatrixRMaj.times(B: ZMatrixRMaj) =
-    ZMatrixRMaj(numRows, B.numCols).also { CommonOps_ZDRM.mult(this, B, it) }
+operator fun CMatrix.times(B: CMatrix) =
+    CMatrix(numRows, B.numCols).also { COps.mult(this, B, it) }
 
 
-/**
- * Useful constants
- */
-typealias Matrix = ZMatrixRMaj
-typealias Ops = CommonOps_ZDRM
-
-/** 2^(-1/2) */
-const val HALF_AMPL = 0.70710678118654757273731092936941422522068023681640625
-
-/** Don't change these constant matrices! */
-val NOT = Matrix(
-    arrayOf(
-        doubleArrayOf(0.0, 0.0, 1.0, 0.0),
-        doubleArrayOf(1.0, 0.0, 0.0, 0.0)
-    )
-)
-val H = Matrix(
-    arrayOf(
-        doubleArrayOf(HALF_AMPL, 0.0, HALF_AMPL, 0.0),
-        doubleArrayOf(HALF_AMPL, 0.0, -HALF_AMPL, 0.0)
-    )
-)
-val I1 = Ops.identity(1)
-val I2 = Ops.identity(2)
-
-val KET0 = Matrix(arrayOf(doubleArrayOf(1.0, 0.0), doubleArrayOf(0.0, 0.0)))
-val KET1 = Matrix(arrayOf(doubleArrayOf(0.0, 0.0), doubleArrayOf(1.0, 0.0)))
-
-val KETBRA0 = Ops.diag(1.0, 0.0, 0.0, 0.0)
-val KETBRA1 = Ops.diag(0.0, 0.0, 1.0, 0.0)
-val SQRT_NOT = Matrix(
-    arrayOf(
-        doubleArrayOf(0.5, 0.5, 0.5, -0.5),
-        doubleArrayOf(0.5, -0.5, 0.5, 0.5),
-    )
-)
-val SQRT_NOT_DAG = ZMatrixRMaj(2, 2).also {
-    Ops.transposeConjugate(SQRT_NOT, it)
-}
-val Z = Matrix(
-    arrayOf(
-        doubleArrayOf(1.0, 0.0, 0.0, 0.0),
-        doubleArrayOf(0.0, 0.0, -1.0, 0.0),
-    )
-)
-val S = Matrix(
-    arrayOf(
-        doubleArrayOf(1.0, 0.0, 0.0, 0.0),
-        doubleArrayOf(0.0, 0.0, 0.0, 1.0),
-    )
-)
-val T = Matrix(
-    arrayOf(
-        doubleArrayOf(1.0, 0.0, 0.0, 0.0),
-        doubleArrayOf(0.0, 0.0, HALF_AMPL, HALF_AMPL),
-    )
-)
-val TDag = ZMatrixRMaj(2, 2).also {
-    Ops.transposeConjugate(T, it)
-}
-
 /** Printing util */
 object MyMatrixIO {
-    fun getMatrixType(mat: Matrix): String {
+    fun getMatrixType(mat: CMatrix): String {
         return if (mat.type == MatrixType.UNSPECIFIED) {
             mat.javaClass.simpleName
         } else {
@@ -124,7 +60,7 @@ object MyMatrixIO {
         }
     }
 
-    fun printTypeSize(out: PrintStream, mat: Matrix) {
+    fun printTypeSize(out: PrintStream, mat: CMatrix) {
         if (mat is MatrixSparse) {
             val m = mat as MatrixSparse
             out.println(
@@ -145,16 +81,16 @@ object MyMatrixIO {
     }
 }
 
-fun ZMatrixRMaj.printFancy2(
+fun CMatrix.printFancy2(
     out: PrintStream = System.out,
     length: Int = MatrixIO.DEFAULT_LENGTH,
     allssr: List<String>
-) = this.let { mat ->
+) = let { mat ->
     MyMatrixIO.printTypeSize(out, mat)
     val format = DecimalFormat("#")
     val builder = StringBuilder(length)
     val cols = mat.numCols
-    val c = Complex_F64()
+    val c = CNumber()
     var i = 0
     for (y in 0 until mat.numRows) {
         for (x in 0 until cols) {
diff --git a/src/main/kotlin/matrix/type.kt b/src/main/kotlin/matrix/type.kt
new file mode 100644
index 0000000..4a6df61
--- /dev/null
+++ b/src/main/kotlin/matrix/type.kt
@@ -0,0 +1,14 @@
+package matrix
+
+import org.ejml.data.Complex_F64
+import org.ejml.data.ZMatrixRMaj
+import org.ejml.dense.row.CommonOps_ZDRM
+
+/**
+ * Created by DEDZTBH on 2020/09/25.
+ * Project KuantumCircuitSim
+ */
+
+typealias CMatrix = ZMatrixRMaj
+typealias COps = CommonOps_ZDRM
+typealias CNumber = Complex_F64
\ No newline at end of file
diff --git a/src/main/kotlin/operator/AllInit.kt b/src/main/kotlin/operator/AllInit.kt
index fc671f4..d1cee21 100644
--- a/src/main/kotlin/operator/AllInit.kt
+++ b/src/main/kotlin/operator/AllInit.kt
@@ -1,6 +1,6 @@
 package operator
 
-import util.*
+import matrix.*
 
 
 /**
diff --git a/src/main/kotlin/operator/TFinder.kt b/src/main/kotlin/operator/TFinder.kt
index 934d44d..3cd5ea7 100644
--- a/src/main/kotlin/operator/TFinder.kt
+++ b/src/main/kotlin/operator/TFinder.kt
@@ -1,6 +1,8 @@
 package operator
 
-import util.*
+import matrix.*
+import readDouble
+import readInt
 import kotlin.math.PI
 import kotlin.math.cos
 import kotlin.math.pow
@@ -25,7 +27,7 @@ open class TFinder(val N: Int) : Operator {
     val allssr = alls.map { arr -> "|${arr.reversed().joinToString("")}>" }
 
     val jointStateSize = alls.size
-    val IN2 = Ops.identity(jointStateSize)
+    val IN2 = COps.identity(jointStateSize)
 
     var opMatrix = IN2
 
@@ -35,19 +37,19 @@ open class TFinder(val N: Int) : Operator {
         }
     }
 
-    val matrix0CtrlCache = Array(N) { HashMap<Matrix, Matrix>() }
-    fun get0CtrlMatrix(i: Int, mat: Matrix, cache: Boolean = true): Matrix {
+    val matrix0CtrlCache = Array(N) { HashMap<CMatrix, CMatrix>() }
+    fun get0CtrlMatrix(i: Int, mat: CMatrix, cache: Boolean = true): CMatrix {
         if (cache) matrix0CtrlCache[i][mat]?.let { return it }
         val res = IKronTable[i] kron mat kron IKronTable[N - i - 1]
         if (cache) matrix0CtrlCache[i][mat] = res
         return res
     }
 
-    val matrix1CtrlCache = Array(N) { Array(N) { HashMap<Matrix, Matrix>() } }
+    val matrix1CtrlCache = Array(N) { Array(N) { HashMap<CMatrix, CMatrix>() } }
 
     /** Control matrix generation based on this great article:
      * http://www.sakkaris.com/tutorials/quantum_control_gates.html */
-    fun get1CtrlMatrix(i: Int, j: Int, mat: Matrix, cache: Boolean = true): Matrix {
+    fun get1CtrlMatrix(i: Int, j: Int, mat: CMatrix, cache: Boolean = true): CMatrix {
         if (cache) matrix1CtrlCache[i][j][mat]?.let { return it }
         val res = get0CtrlMatrix(i, KETBRA0) + when {
             i < j -> IKronTable[i] kron KETBRA1 kron
@@ -62,7 +64,7 @@ open class TFinder(val N: Int) : Operator {
         return res
     }
 
-    fun getCCNotMatrix(i: Int, j: Int, k: Int): Matrix {
+    fun getCCNotMatrix(i: Int, j: Int, k: Int): CMatrix {
         /** using Sleator-Weinfurter construction
          * matrix is actually reverse order as graph */
         val cnotij = get1CtrlMatrix(i, j, NOT)
@@ -94,7 +96,7 @@ open class TFinder(val N: Int) : Operator {
         val i = readInt()
         val newOp = when (cmd) {
             "Not" -> get0CtrlMatrix(i, NOT)
-            "Hadamard" -> get0CtrlMatrix(i, H)
+            "Hadamard", "H" -> get0CtrlMatrix(i, H)
             "CNot" -> get1CtrlMatrix(i, readInt(), NOT)
             "Swap" -> {
                 /** https://algassert.com/post/1717
@@ -126,7 +128,7 @@ open class TFinder(val N: Int) : Operator {
                 val rad = readDouble() * PI / 180
                 val sine = sin(rad)
                 val cosine = cos(rad)
-                val rotMat = Matrix(
+                val rotMat = CMatrix(
                     arrayOf(
                         doubleArrayOf(cosine, 0.0, -sine, 0.0),
                         doubleArrayOf(sine, 0.0, cosine, 0.0)
diff --git a/src/main/kotlin/operator/Tester.kt b/src/main/kotlin/operator/Tester.kt
index 0d9b6b9..1063c7b 100644
--- a/src/main/kotlin/operator/Tester.kt
+++ b/src/main/kotlin/operator/Tester.kt
@@ -1,7 +1,11 @@
 package operator
 
-import org.ejml.data.Complex_F64
-import util.*
+import matrix.CMatrix
+import matrix.CNumber
+import matrix.printFancy2
+import matrix.times
+import readInt
+import upperBound
 import kotlin.random.Random
 
 /**
@@ -11,7 +15,7 @@ import kotlin.random.Random
 class Tester(N: Int) : TFinder(N) {
     /** 2^N by 1 column vector */
     val jointState =
-        Matrix(jointStateSize, 1).apply { set(0, 0, 1.0, 0.0) }
+        CMatrix(jointStateSize, 1).apply { set(0, 0, 1.0, 0.0) }
 
     override fun runCmd(cmd: String): Int {
         if (cmd == "Measure") {
@@ -19,7 +23,7 @@ class Tester(N: Int) : TFinder(N) {
             val results = opMatrix * jointState
             val probs = mutableListOf(0.0)
             val labels = mutableListOf<Int>()
-            val a = Complex_F64()
+            val a = CNumber()
             for (j in 0 until results.numRows) {
                 results.get(0, j, a)
                 if (a.magnitude2 > 0) {
@@ -34,7 +38,7 @@ class Tester(N: Int) : TFinder(N) {
                 println(allssr[labels[index]])
             }
             println()
-            return 1
+            return 0
         }
         return super.runCmd(cmd)
     }
diff --git a/src/main/kotlin/util.kt b/src/main/kotlin/util.kt
new file mode 100644
index 0000000..56f4578
--- /dev/null
+++ b/src/main/kotlin/util.kt
@@ -0,0 +1,37 @@
+import java.io.BufferedReader
+import java.util.*
+
+/**
+ * Created by DEDZTBH on 2020/09/25.
+ * Project KuantumCircuitSim
+ */
+
+/** IO */
+var _tokenizer: StringTokenizer = StringTokenizer("")
+fun read(): String {
+    while (!_tokenizer.hasMoreTokens())
+        _tokenizer = StringTokenizer(reader.readLine() ?: return "", " ")
+    return _tokenizer.nextToken()
+}
+
+fun readInt() = read().toInt()
+fun readDouble() = read().toDouble()
+fun readFloat() = read().toFloat()
+
+lateinit var reader: BufferedReader
+
+
+/** Algorithm */
+fun <T : Comparable<T>> List<T>.upperBound(key: T): Int {
+    var low = 0
+    var high = size - 1
+    while (low < high) {
+        val mid = low + (high - low + 1) / 2
+        if (this[mid] <= key) {
+            low = mid
+        } else {
+            high = mid - 1
+        }
+    }
+    return low
+}
\ No newline at end of file
diff --git a/src/main/kotlin/util/algo.kt b/src/main/kotlin/util/algo.kt
deleted file mode 100644
index c54ff06..0000000
--- a/src/main/kotlin/util/algo.kt
+++ /dev/null
@@ -1,20 +0,0 @@
-package util
-
-/**
- * Created by DEDZTBH on 2020/09/22.
- * Project KuantumCircuitSim
- */
-
-fun <T : Comparable<T>> List<T>.upperBound(key: T): Int {
-    var low = 0
-    var high = size - 1
-    while (low < high) {
-        val mid = low + (high - low + 1) / 2
-        if (this[mid] <= key) {
-            low = mid
-        } else {
-            high = mid - 1
-        }
-    }
-    return low
-}
\ No newline at end of file
diff --git a/src/main/kotlin/util/io.kt b/src/main/kotlin/util/io.kt
deleted file mode 100644
index c1d8cfd..0000000
--- a/src/main/kotlin/util/io.kt
+++ /dev/null
@@ -1,20 +0,0 @@
-package util
-
-import java.io.BufferedReader
-import java.util.*
-
-/**
- * Created by DEDZTBH on 2020/09/22.
- * Project KuantumCircuitSim
- */
-var _tokenizer: StringTokenizer = StringTokenizer("")
-fun read(): String {
-    while (!_tokenizer.hasMoreTokens())
-        _tokenizer = StringTokenizer(reader.readLine() ?: return "", " ")
-    return _tokenizer.nextToken()
-}
-
-fun readInt() = read().toInt()
-fun readDouble() = read().toDouble()
-
-lateinit var reader: BufferedReader
\ No newline at end of file
diff --git a/src/test/kotlin/performance.kt b/src/test/kotlin/performance.kt
new file mode 100644
index 0000000..379cd5d
--- /dev/null
+++ b/src/test/kotlin/performance.kt
@@ -0,0 +1,86 @@
+import org.ejml.data.CMatrixRMaj
+import org.ejml.data.Complex_F32
+import org.ejml.data.Complex_F64
+import org.ejml.data.ZMatrixRMaj
+import org.ejml.dense.row.CommonOps_CDRM
+import org.ejml.dense.row.CommonOps_ZDRM
+import kotlin.random.Random
+import kotlin.system.measureTimeMillis
+
+/**
+ * Created by DEDZTBH on 2020/09/25.
+ * Project KuantumCircuitSim
+ */
+
+infix fun ZMatrixRMaj.kron(B: ZMatrixRMaj) = let { A ->
+    val numColsC = A.numCols * B.numCols
+    val numRowsC = A.numRows * B.numRows
+
+    val C = ZMatrixRMaj(numRowsC, numColsC)
+
+    val acomplex = Complex_F64()
+    val bcomplex = Complex_F64()
+    for (i in 0 until A.numRows) {
+        for (j in 0 until A.numCols) {
+            A.get(i, j, acomplex)
+            for (rowB in 0 until B.numRows) {
+                for (colB in 0 until B.numCols) {
+                    B.get(rowB, colB, bcomplex)
+                    val aB = acomplex.times(bcomplex)
+                    C.set(i * B.numRows + rowB, j * B.numCols + colB, aB.real, aB.imaginary)
+                }
+            }
+        }
+    }
+
+    C
+}
+
+infix fun CMatrixRMaj.kron(B: CMatrixRMaj) = let { A ->
+    val numColsC = A.numCols * B.numCols
+    val numRowsC = A.numRows * B.numRows
+
+    val C = CMatrixRMaj(numRowsC, numColsC)
+
+    val acomplex = Complex_F32()
+    val bcomplex = Complex_F32()
+    for (i in 0 until A.numRows) {
+        for (j in 0 until A.numCols) {
+            A.get(i, j, acomplex)
+            for (rowB in 0 until B.numRows) {
+                for (colB in 0 until B.numCols) {
+                    B.get(rowB, colB, bcomplex)
+                    val aB = acomplex.times(bcomplex)
+                    C.set(i * B.numRows + rowB, j * B.numCols + colB, aB.real, aB.imaginary)
+                }
+            }
+        }
+    }
+
+    C
+}
+
+fun main() {
+    val N = 100
+    val t = 5
+
+    val ad = ZMatrixRMaj(Array(N) { DoubleArray(N * 2) { Random.nextDouble() } })
+    val bd = ZMatrixRMaj(Array(N) { DoubleArray(N * 2) { Random.nextDouble() } })
+
+    val af = CMatrixRMaj(Array(N) { FloatArray(N * 2) { Random.nextFloat() } })
+    val bf = CMatrixRMaj(Array(N) { FloatArray(N * 2) { Random.nextFloat() } })
+
+    LongArray(t) {
+        measureTimeMillis {
+            val c = ZMatrixRMaj(ad.numRows, bd.numCols).also { CommonOps_ZDRM.mult(ad, bd, it) }
+            val d = ad kron bd
+        }
+    }.average().let(::println)
+
+    LongArray(t) {
+        measureTimeMillis {
+            val c = CMatrixRMaj(af.numRows, bf.numCols).also { CommonOps_CDRM.mult(af, bf, it) }
+            val d = af kron bf
+        }
+    }.average().let(::println)
+}
\ No newline at end of file