Update ControlledModExp2 with bitwise operators

itsubaki · Jan 8, 2025 · 53ac15c · 53ac15c
1 parent 462a421
commit 53ac15c
Show file tree

Hide file tree

Showing 3 changed files with 27 additions and 30 deletions.
diff --git a/quantum/gate/gate.go b/quantum/gate/gate.go
@@ -1,10 +1,8 @@
 package gate
 
 import (
-	"fmt"
 	"math"
 	"math/cmplx"
-	"strconv"
 
 	"github.com/itsubaki/q/math/matrix"
 	"github.com/itsubaki/q/math/number"
@@ -329,36 +327,34 @@ func QFT(n int) matrix.Matrix {
 // len(t) must be larger than log2(N).
 func ControlledModExp2(n, a, j, N, c int, t []int) matrix.Matrix {
 	m := I(n)
-	d, _ := m.Dimension()
-
-	r0len, r1len := n-len(t), len(t)
+	r1len := len(t)
 	a2jmodN := number.ModExp2(a, j, N)
 
-	idx := make([]int64, d)
+	d, _ := m.Dimension()
+	idx := make([]int, d)
 	for i := range d {
-		bits := []rune(fmt.Sprintf("%0*b", n, i))
-
-		if bits[c] == '0' {
-			idx[i] = int64(i)
+		if (i>>(n-1-c))&1 == 0 {
+			// control bit is 0, then do nothing
+			idx[i] = i
 			continue
 		}
 
-		k := number.Must(strconv.ParseInt(string(bits[r0len:]), 2, 0))
-		if k > int64(N-1) {
-			idx[i] = int64(i)
+		// r1len bits of i
+		mask := (1 << r1len) - 1
+		k := i & mask
+		if k > N-1 {
+			idx[i] = i
 			continue
 		}
 
-		a2jkmodN := (int64(a2jmodN) * k) % int64(N)
-		a2jkmodNs := []rune(fmt.Sprintf("%0*b", r1len, a2jkmodN))
-		newbits := append(bits[:r0len], []rune(a2jkmodNs)...)
-
-		idx[i] = number.Must(strconv.ParseInt(string(newbits), 2, 0))
+		// r0len bits of i + a2jkmodN bits
+		a2jkmodN := a2jmodN * k % N
+		idx[i] = (i >> r1len << r1len) | a2jkmodN
 	}
 
 	g := make(matrix.Matrix, d)
-	for i, ii := range idx {
-		g[ii] = m[i]
+	for i, j := range idx {
+		g[j] = m[i]
 	}
 
 	return g

diff --git a/quantum/gate/gate_test.go b/quantum/gate/gate_test.go
@@ -113,13 +113,13 @@ func ExampleControlledModExp2() {
 			continue
 		}
 
-		for ii, e := range r {
+		for l, e := range r {
 			if e == complex(0, 0) {
 				continue
 			}
 
 			// decimal number representation of a^2^j * k mod N
-			a2jkmodNs := fmt.Sprintf("%0*s", n, strconv.FormatInt(int64(ii), 2)[1:])
+			a2jkmodNs := fmt.Sprintf("%0*s", n, strconv.FormatInt(int64(l), 2)[1:])
 			a2jkmodN := number.Must(strconv.ParseInt(a2jkmodNs, 2, 64))
 			got := (int64(number.ModExp2(a, j, N)) * k) % int64(N)
 

diff --git a/quantum/qubit/qubit.go b/quantum/qubit/qubit.go
@@ -156,12 +156,12 @@ func (q *Qubit) Probability() []float64 {
 // Measure returns a measured qubit.
 func (q *Qubit) Measure(index int) *Qubit {
 	n := q.NumberOfBit()
+	mask := 1 << (n - 1 - index)
 
-	zidx, oidx := make([]int, 0), make([]int, 0)
-	zprop := make([]float64, 0)
+	zidx, zprop := make([]int, 0), make([]float64, 0)
+	oidx := make([]int, 0)
 	for i, p := range q.Probability() {
-		bits := []rune(fmt.Sprintf("%0*b", n, i))
-		if bits[index] == '0' {
+		if i&mask == 0 {
 			zidx, zprop = append(zidx, i), append(zprop, p)
 			continue
 		}
@@ -235,10 +235,9 @@ func (q *Qubit) State(index ...[]int) []State {
 			continue
 		}
 
-		bit := fmt.Sprintf("%0*b", n, i)
 		var bin []string
 		for _, idx := range index {
-			bin = append(bin, take(bit, idx))
+			bin = append(bin, take(n, i, idx))
 		}
 
 		state = append(state, NewState(amp, bin...))
@@ -261,10 +260,12 @@ func round(a complex128, eps ...float64) complex128 {
 	return a
 }
 
-func take(binary string, index []int) string {
+func take(n, i int, index []int) string {
+	s := fmt.Sprintf("%0*b", n, i)
+
 	var sb strings.Builder
 	for _, i := range index {
-		sb.WriteString(binary[i : i+1])
+		sb.WriteString(s[i : i+1])
 	}
 
 	return sb.String()