diff --git a/pkg/math/matrix/matrix.go b/pkg/math/matrix/matrix.go
index 8e4d50c..fd0daf7 100644
--- a/pkg/math/matrix/matrix.go
+++ b/pkg/math/matrix/matrix.go
@@ -95,7 +95,7 @@ func (m Matrix) Conjugate() Matrix {
 	return out
 }
 
-// Dagger returns conjugate transpose matrix
+// Dagger returns conjugate transpose matrix.
 func (m Matrix) Dagger() Matrix {
 	p, q := m.Dimension()
 
diff --git a/pkg/quantum/gate/gate.go b/pkg/quantum/gate/gate.go
index 172d6c6..35c14c1 100644
--- a/pkg/quantum/gate/gate.go
+++ b/pkg/quantum/gate/gate.go
@@ -27,6 +27,11 @@ func Empty(n ...int) []matrix.Matrix {
 	return make([]matrix.Matrix, n[0])
 }
 
+// Theta returns 2 * pi / 2**k
+func Theta(k int) float64 {
+	return 2 * math.Pi / math.Pow(2, float64(k))
+}
+
 func U(theta, phi, lambda float64) matrix.Matrix {
 	v := complex(theta/2, 0)
 	return matrix.Matrix{
@@ -35,39 +40,6 @@ func U(theta, phi, lambda float64) matrix.Matrix {
 	}
 }
 
-func RX(theta float64) matrix.Matrix {
-	v := complex(theta/2, 0)
-	return matrix.Matrix{
-		[]complex128{cmplx.Cos(v), -1i * cmplx.Sin(v)},
-		[]complex128{-1i * cmplx.Sin(v), cmplx.Cos(v)},
-	}
-}
-
-func RY(theta float64) matrix.Matrix {
-	v := complex(theta/2, 0)
-	return matrix.Matrix{
-		[]complex128{cmplx.Cos(v), -1 * cmplx.Sin(v)},
-		[]complex128{cmplx.Sin(v), cmplx.Cos(v)},
-	}
-}
-
-func RZ(theta float64) matrix.Matrix {
-	v := complex(0, theta/2)
-	return matrix.Matrix{
-		[]complex128{cmplx.Exp(-1 * v), 0},
-		[]complex128{0, cmplx.Exp(v)},
-	}
-}
-
-func R(k int) matrix.Matrix {
-	p := 2 * math.Pi / math.Pow(2, float64(k))
-	e := cmplx.Exp(complex(0, p))
-	return matrix.Matrix{
-		[]complex128{1, 0},
-		[]complex128{0, e},
-	}
-}
-
 func I(n ...int) matrix.Matrix {
 	return matrix.TensorProductN(matrix.Matrix{
 		[]complex128{1, 0},
@@ -119,6 +91,38 @@ func T(n ...int) matrix.Matrix {
 	}, n...)
 }
 
+func R(theta float64) matrix.Matrix {
+	e := cmplx.Exp(complex(0, theta))
+	return matrix.Matrix{
+		[]complex128{1, 0},
+		[]complex128{0, e},
+	}
+}
+
+func RX(theta float64) matrix.Matrix {
+	v := complex(theta/2, 0)
+	return matrix.Matrix{
+		[]complex128{cmplx.Cos(v), -1i * cmplx.Sin(v)},
+		[]complex128{-1i * cmplx.Sin(v), cmplx.Cos(v)},
+	}
+}
+
+func RY(theta float64) matrix.Matrix {
+	v := complex(theta/2, 0)
+	return matrix.Matrix{
+		[]complex128{cmplx.Cos(v), -1 * cmplx.Sin(v)},
+		[]complex128{cmplx.Sin(v), cmplx.Cos(v)},
+	}
+}
+
+func RZ(theta float64) matrix.Matrix {
+	v := complex(0, theta/2)
+	return matrix.Matrix{
+		[]complex128{cmplx.Exp(-1 * v), 0},
+		[]complex128{0, cmplx.Exp(v)},
+	}
+}
+
 func Controlled(u matrix.Matrix, n int, c []int, t int) matrix.Matrix {
 	g := I([]int{n}...)
 	f := fmt.Sprintf("%s%s%s", "%0", strconv.Itoa(n), "s")
@@ -288,13 +292,13 @@ func CS(n, c, t int) matrix.Matrix {
 	return ControlledS(n, []int{c}, t)
 }
 
-func ControlledR(k, n int, c []int, t int) matrix.Matrix {
+func ControlledR(theta float64, n int, c []int, t int) matrix.Matrix {
 	g := I([]int{n}...)
 	d, _ := g.Dimension()
 	f := fmt.Sprintf("%s%s%s", "%0", strconv.Itoa(n), "s")
 
-	p := 2 * math.Pi / math.Pow(2, float64(k))
-	e := cmplx.Exp(complex(0, p))
+	// exp(i * theta)
+	e := cmplx.Exp(complex(0, theta))
 
 	for i := 0; i < d; i++ {
 		bits := []rune(fmt.Sprintf(f, strconv.FormatInt(int64(i), 2)))
@@ -316,8 +320,8 @@ func ControlledR(k, n int, c []int, t int) matrix.Matrix {
 	return g
 }
 
-func CR(k, n, c, t int) matrix.Matrix {
-	return ControlledR(k, n, []int{c}, t)
+func CR(theta float64, n, c, t int) matrix.Matrix {
+	return ControlledR(theta, n, []int{c}, t)
 }
 
 // Swap returns unitary matrix of Swap operation.
@@ -355,7 +359,7 @@ func QFT(n int) matrix.Matrix {
 
 		k := 2
 		for j := i + 1; j < n; j++ {
-			g = g.Apply(CR(k, n, j, i))
+			g = g.Apply(CR(Theta(k), n, j, i))
 			k++
 		}
 	}
@@ -363,8 +367,8 @@ func QFT(n int) matrix.Matrix {
 	return g
 }
 
-// ControlledModExp2 returns unitary matrix of controlled modular exponentiation operation. |j>|k> -> |j>|a**(2**j) * k mod N>
-// len(t) must be larger than log2(N)
+// ControlledModExp2 returns unitary matrix of controlled modular exponentiation operation. |j>|k> -> |j>|a**(2**j) * k mod N>.
+// len(t) must be larger than log2(N).
 func ControlledModExp2(n, a, j, N, c int, t []int) matrix.Matrix {
 	m := I([]int{n}...)
 	d, _ := m.Dimension()
diff --git a/pkg/quantum/gate/gate_test.go b/pkg/quantum/gate/gate_test.go
index a68afea..7d3325d 100644
--- a/pkg/quantum/gate/gate_test.go
+++ b/pkg/quantum/gate/gate_test.go
@@ -208,7 +208,7 @@ func TestC(t *testing.T) {
 		{gate.C(gate.X(), 2, 0, 1), gate.CNOT(2, 0, 1)},
 		{gate.C(gate.Z(), 2, 0, 1), gate.CZ(2, 0, 1)},
 		{gate.C(gate.S(), 2, 0, 1), gate.CS(2, 0, 1)},
-		{gate.C(gate.R(100), 2, 0, 1), gate.CR(100, 2, 0, 1)},
+		{gate.C(gate.R(gate.Theta(4)), 2, 0, 1), gate.CR(gate.Theta(4), 2, 0, 1)},
 		{gate.C(gate.X(), 3, 0, 2), gate.CNOT(3, 0, 2)},
 		{gate.C(gate.X(), 3, 0, 1), gate.CNOT(3, 0, 1)},
 		{gate.C(gate.X(), 3, 1, 0), gate.CNOT(3, 1, 0)},
@@ -291,14 +291,14 @@ func TestIsUnitary(t *testing.T) {
 		{gate.Y(), true},
 		{gate.Z(), true},
 		{gate.U(1, 2, 3), true},
-		{gate.R(4), true},
+		{gate.R(gate.Theta(4)), true},
 		{gate.RX(1.23), true},
 		{gate.RY(1.23), true},
 		{gate.RZ(1.23), true},
 		{gate.ControlledS(2, []int{0}, 1), true},
-		{gate.ControlledR(10, 2, []int{0}, 1), true},
+		{gate.ControlledR(gate.Theta(4), 2, []int{0}, 1), true},
 		{gate.CS(2, 0, 1), true},
-		{gate.CR(10, 2, 0, 1), true},
+		{gate.CR(gate.Theta(4), 2, 0, 1), true},
 		{gate.QFT(2), true},
 		{gate.New(
 			[]complex128{1, 2},
diff --git a/q.go b/q.go
index 1eca225..bfc2dfa 100644
--- a/q.go
+++ b/q.go
@@ -115,12 +115,12 @@ func (q *Q) Reset(qb ...Qubit) {
 	}
 }
 
-func (q *Q) I(qb ...Qubit) *Q {
-	return q.Apply(gate.I(), qb...)
+func (q *Q) U(theta, phi, lambda float64, qb ...Qubit) *Q {
+	return q.Apply(gate.U(theta, phi, lambda), qb...)
 }
 
-func (q *Q) H(qb ...Qubit) *Q {
-	return q.Apply(gate.H(), qb...)
+func (q *Q) I(qb ...Qubit) *Q {
+	return q.Apply(gate.I(), qb...)
 }
 
 func (q *Q) X(qb ...Qubit) *Q {
@@ -135,6 +135,10 @@ func (q *Q) Z(qb ...Qubit) *Q {
 	return q.Apply(gate.Z(), qb...)
 }
 
+func (q *Q) H(qb ...Qubit) *Q {
+	return q.Apply(gate.H(), qb...)
+}
+
 func (q *Q) S(qb ...Qubit) *Q {
 	return q.Apply(gate.S(), qb...)
 }
@@ -143,8 +147,8 @@ func (q *Q) T(qb ...Qubit) *Q {
 	return q.Apply(gate.T(), qb...)
 }
 
-func (q *Q) U(theta, phi, lambda float64, qb ...Qubit) *Q {
-	return q.Apply(gate.U(theta, phi, lambda), qb...)
+func (q *Q) R(theta float64, qb ...Qubit) *Q {
+	return q.Apply(gate.R(theta), qb...)
 }
 
 func (q *Q) RX(theta float64, qb ...Qubit) *Q {
@@ -199,26 +203,27 @@ func (q *Q) C(m matrix.Matrix, control, target Qubit) *Q {
 	return q.Controlled(m, []Qubit{control}, target)
 }
 
-func (q *Q) ControlledR(k int, control []Qubit, target Qubit) *Q {
+func (q *Q) ControlledNot(control []Qubit, target Qubit) *Q {
 	n := q.NumberOfBit()
-	g := gate.ControlledR(k, n, Index(control...), target.Index())
+	g := gate.ControlledNot(n, Index(control...), target.Index())
 	q.qb.Apply(g)
 	return q
 }
 
-func (q *Q) CR(k int, control, target Qubit) *Q {
-	return q.ControlledR(k, []Qubit{control}, target)
+func (q *Q) CNOT(control, target Qubit) *Q {
+	return q.ControlledNot([]Qubit{control}, target)
 }
 
-func (q *Q) InverseCR(k int, control, target Qubit) *Q {
-	n := q.NumberOfBit()
-	g := gate.ControlledR(k, n, []int{control.Index()}, target.Index()).Dagger()
-	q.qb.Apply(g)
-	return q
+func (q *Q) CCNOT(control0, control1, target Qubit) *Q {
+	return q.ControlledNot([]Qubit{control0, control1}, target)
+}
+
+func (q *Q) CCCNOT(control0, control1, control2, target Qubit) *Q {
+	return q.ControlledNot([]Qubit{control0, control1, control2}, target)
 }
 
-func (q *Q) InvCR(k int, control, target Qubit) *Q {
-	return q.InverseCR(k, control, target)
+func (q *Q) Toffoli(control0, control1, target Qubit) *Q {
+	return q.CCNOT(control0, control1, target)
 }
 
 func (q *Q) ControlledZ(control []Qubit, target Qubit) *Q {
@@ -236,27 +241,30 @@ func (q *Q) CCZ(control0, control1, target Qubit) *Q {
 	return q.ControlledZ([]Qubit{control0, control1}, target)
 }
 
-func (q *Q) ControlledNot(control []Qubit, target Qubit) *Q {
+func (q *Q) ControlledR(theta float64, control []Qubit, target Qubit) *Q {
 	n := q.NumberOfBit()
-	g := gate.ControlledNot(n, Index(control...), target.Index())
+	g := gate.ControlledR(theta, n, Index(control...), target.Index())
 	q.qb.Apply(g)
 	return q
 }
 
-func (q *Q) CNOT(control, target Qubit) *Q {
-	return q.ControlledNot([]Qubit{control}, target)
+func (q *Q) CR(theta float64, control, target Qubit) *Q {
+	return q.ControlledR(theta, []Qubit{control}, target)
 }
 
-func (q *Q) CCNOT(control0, control1, target Qubit) *Q {
-	return q.ControlledNot([]Qubit{control0, control1}, target)
+func (q *Q) ControlledModExp2(a, j, N int, control Qubit, target []Qubit) *Q {
+	n := q.NumberOfBit()
+	g := gate.ControlledModExp2(n, a, j, N, control.Index(), Index(target...))
+	q.qb.Apply(g)
+	return q
 }
 
-func (q *Q) CCCNOT(control0, control1, control2, target Qubit) *Q {
-	return q.ControlledNot([]Qubit{control0, control1, control2}, target)
-}
+func (q *Q) CModExp2(a, N int, control []Qubit, target []Qubit) *Q {
+	for i := 0; i < len(control); i++ {
+		q.ControlledModExp2(a, i, N, control[i], target)
+	}
 
-func (q *Q) Toffoli(control0, control1, target Qubit) *Q {
-	return q.CCNOT(control0, control1, target)
+	return q
 }
 
 func (q *Q) Condition(condition bool, m matrix.Matrix, qb ...Qubit) *Q {
@@ -275,21 +283,6 @@ func (q *Q) ConditionZ(condition bool, qb ...Qubit) *Q {
 	return q.Condition(condition, gate.Z(), qb...)
 }
 
-func (q *Q) ControlledModExp2(a, j, N int, control Qubit, target []Qubit) *Q {
-	n := q.NumberOfBit()
-	g := gate.ControlledModExp2(n, a, j, N, control.Index(), Index(target...))
-	q.qb.Apply(g)
-	return q
-}
-
-func (q *Q) CModExp2(a, N int, control []Qubit, target []Qubit) *Q {
-	for i := 0; i < len(control); i++ {
-		q.ControlledModExp2(a, i, N, control[i], target)
-	}
-
-	return q
-}
-
 func (q *Q) Swap(qb ...Qubit) *Q {
 	n := q.NumberOfBit()
 	l := len(qb)
@@ -310,7 +303,7 @@ func (q *Q) QFT(qb ...Qubit) *Q {
 
 		k := 2
 		for j := i + 1; j < l; j++ {
-			q.CR(k, qb[j], qb[i])
+			q.CR(gate.Theta(k), qb[j], qb[i])
 			k++
 		}
 	}
@@ -323,7 +316,7 @@ func (q *Q) InverseQFT(qb ...Qubit) *Q {
 	for i := l - 1; i > -1; i-- {
 		k := l - i
 		for j := l - 1; j > i; j-- {
-			q.InvCR(k, qb[j], qb[i])
+			q.CR(-1*gate.Theta(k), qb[j], qb[i])
 			k--
 		}
 
diff --git a/q_test.go b/q_test.go
index 432a729..43e0977 100644
--- a/q_test.go
+++ b/q_test.go
@@ -221,6 +221,20 @@ func ExampleQ_Apply() {
 	// [11][  3]( 0.7071 0.0000i): 0.5000
 }
 
+func ExampleQ_U() {
+	qsim := q.New()
+
+	q0 := qsim.Zero()
+	qsim.U(math.Pi, 0, math.Pi, q0)
+
+	for _, s := range qsim.State() {
+		fmt.Println(s)
+	}
+
+	// Output:
+	// [1][  1]( 1.0000 0.0000i): 1.0000
+}
+
 func ExampleQ_I() {
 	qsim := q.New()
 
@@ -305,18 +319,18 @@ func ExampleQ_T() {
 	// [1][  1]( 0.7071 0.7071i): 1.0000
 }
 
-func ExampleQ_U() {
+func ExampleQ_R() {
 	qsim := q.New()
 
-	q0 := qsim.Zero()
-	qsim.U(math.Pi, 0, math.Pi, q0)
+	q0 := qsim.One()
+	qsim.R(2*math.Pi/4, q0)
 
 	for _, s := range qsim.State() {
 		fmt.Println(s)
 	}
 
 	// Output:
-	// [1][  1]( 1.0000 0.0000i): 1.0000
+	// [1][  1]( 0.0000 1.0000i): 1.0000
 }
 
 func ExampleQ_RX() {
@@ -874,10 +888,10 @@ func Example_qFT() {
 	q2 := qsim.Zero()
 
 	qsim.H(q0)
-	qsim.CR(2, q1, q0)
-	qsim.CR(3, q2, q0)
+	qsim.CR(gate.Theta(2), q1, q0)
+	qsim.CR(gate.Theta(3), q2, q0)
 	qsim.H(q1)
-	qsim.CR(2, q2, q1)
+	qsim.CR(gate.Theta(2), q2, q1)
 	qsim.H(q2)
 	qsim.Swap(q0, q2)
 
@@ -1014,16 +1028,17 @@ func Example_shorFactoring51() {
 	qsim.CNOT(q2, q6)
 	qsim.CNOT(q3, q7)
 
+	// inverse QFT
 	qsim.Swap(q0, q1, q2, q3)
 	qsim.H(q3)
-	qsim.CR(2, q3, q2)
+	qsim.CR(-1*gate.Theta(2), q3, q2)
 	qsim.H(q2)
-	qsim.CR(3, q3, q1)
-	qsim.CR(2, q2, q1)
+	qsim.CR(-1*gate.Theta(3), q3, q1)
+	qsim.CR(-1*gate.Theta(2), q2, q1)
 	qsim.H(q1)
-	qsim.CR(4, q3, q0)
-	qsim.CR(3, q2, q0)
-	qsim.CR(2, q1, q0)
+	qsim.CR(-1*gate.Theta(4), q3, q0)
+	qsim.CR(-1*gate.Theta(3), q2, q0)
+	qsim.CR(-1*gate.Theta(2), q1, q0)
 	qsim.H(q0)
 
 	m := qsim.Measure(q0, q1, q2, q3).BinaryString()
@@ -1070,11 +1085,12 @@ func Example_shorFactoring85() {
 	qsim.CNOT(q2, q6)
 	qsim.CNOT(q3, q7)
 
+	// inverse QFT
 	qsim.Swap(q0, q1, q2, q3)
 	qsim.H(q3)
-	qsim.CR(2, q3, q2).H(q2)
-	qsim.CR(3, q3, q1).CR(2, q2, q1).H(q1)
-	qsim.CR(4, q3, q0).CR(3, q2, q0).CR(2, q1, q0).H(q0)
+	qsim.CR(-1*gate.Theta(2), q3, q2).H(q2)
+	qsim.CR(-1*gate.Theta(3), q3, q1).CR(-1*gate.Theta(2), q2, q1).H(q1)
+	qsim.CR(-1*gate.Theta(4), q3, q0).CR(-1*gate.Theta(3), q2, q0).CR(-1*gate.Theta(2), q1, q0).H(q0)
 
 	m := qsim.Measure(q0, q1, q2, q3).BinaryString()
 	s, r, d, ok := number.FindOrder(a, N, fmt.Sprintf("0.%s", m))