Add ABC function

cmd/counting/main.go

@@ -53,9 +53,15 @@ func main() {
 	if err != nil {
 		panic(err)
 	}
-	defer f.Close()
+	defer func() {
+		if err := f.Close(); err != nil {
+			fmt.Println(err)
+		}
+	}()
 
-	pprof.StartCPUProfile(f)
+	if err := pprof.StartCPUProfile(f); err != nil {
+		panic(err)
+	}
 	defer pprof.StopCPUProfile()
 
 	// flags

quantum/gate/gate.go

@@ -27,6 +27,21 @@ func U(theta, phi, lambda float64) *matrix.Matrix {
 	)
 }
 
+// SU returns a special unitary gate.
+func SU(theta, phi, lambda float64) *matrix.Matrix {
+	phase := cmplx.Exp(complex(0, -1*(phi+lambda)/2))
+	return U(theta, phi, lambda).Mul(phase)
+}
+
+// ABC returns alpha, A, B, C such that ABC = I, exp(i * alpha) * AXBXC = U(theta, phi, lambda).
+func ABC(theta, phi, lambda float64) (float64, *matrix.Matrix, *matrix.Matrix, *matrix.Matrix) {
+	alpha := (phi + lambda) / 2
+	A := matrix.MatMul(RZ(phi), RY(theta/2))
+	B := matrix.MatMul(RY(-theta/2), RZ(-(lambda+phi)/2))
+	C := RZ((lambda - phi) / 2)
+	return alpha, A, B, C
+}
+
 // I returns an identity gate.
 func I(n ...int) *matrix.Matrix {
 	return matrix.TensorProductN(matrix.New(

quantum/gate/gate_test.go

@@ -3,6 +3,7 @@ package gate_test
 import (
 	"fmt"
 	"math"
+	"math/cmplx"
 	"testing"
 
 	"github.com/itsubaki/q/math/matrix"
@@ -105,6 +106,63 @@ func ExampleTensorProduct() {
 	// [(0+0i) (1+0i) (0+0i) (0+0i)]
 }
 
+func ExampleSU() {
+	u := gate.U(math.Pi/4, math.Pi/2, math.Pi/8)
+	su := gate.SU(math.Pi/4, math.Pi/2, math.Pi/8)
+
+	usud := matrix.MatMul(u, su.Dagger())
+	fmt.Println(usud.Mul(1 / usud.Data[0]).Equal(gate.I()))
+
+	det := func(m *matrix.Matrix) complex128 {
+		return m.Data[0]*m.Data[3] - m.Data[1]*m.Data[2]
+	}
+
+	fmt.Printf("%.4f\n", det(u))
+	fmt.Printf("%.4f\n", det(su))
+
+	// Output:
+	// true
+	// (-0.3827+0.9239i)
+	// (1.0000+0.0000i)
+}
+
+func ExampleABC() {
+	theta, phi, lambda := math.Pi/2, math.Pi/4, math.Pi/8
+
+	alpha, a, b, c := gate.ABC(theta, phi, lambda)
+	fmt.Println(matrix.MatMul(a, b, c).Equal(gate.I()))
+
+	axbxc := matrix.MatMul(a, gate.X(), b, gate.X(), c)
+	fmt.Println(axbxc.Equal(gate.SU(theta, phi, lambda)))
+
+	axbxcp := axbxc.Mul(cmplx.Exp(complex(0, alpha)))
+	fmt.Println(axbxcp.Equal(gate.U(theta, phi, lambda)))
+
+	// Output:
+	// true
+	// true
+	// true
+}
+
+func FuzzABC(f *testing.F) {
+	f.Add(0.0, 0.0, 0.0)
+	f.Add(math.Pi, math.Pi, math.Pi)
+	f.Add(math.Pi/2, math.Pi/4, math.Pi/8)
+
+	f.Fuzz(func(t *testing.T, theta, phi, lambda float64) {
+		alpha, a, b, c := gate.ABC(theta, phi, lambda)
+		if !matrix.MatMul(a, b, c).Equal(gate.I()) {
+			t.Fatalf("ABC != I (theta=%v phi=%v lambda=%v)", theta, phi, lambda)
+		}
+
+		phase := cmplx.Exp(complex(0, alpha))
+		axbxc := matrix.MatMul(a, gate.X(), b, gate.X(), c).Mul(phase)
+		if !axbxc.Equal(gate.U(theta, phi, lambda)) {
+			t.Fatalf("AXBXC != U (theta=%v phi=%v lambda=%v)", theta, phi, lambda)
+		}
+	})
+}
+
 func TestU(t *testing.T) {
 	cases := []struct {
 		in, want *matrix.Matrix