Optionally check the solution at the end of the run

examples/interface.rs

@@ -19,6 +19,6 @@ fn main() {
     let seed = 1234;
     let verbose = false;
     let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed, verbose);
-    let gates = gridsynth_gates(&mut gridsynth_config);
-    println!("{}", gates);
+    let res = gridsynth_gates(&mut gridsynth_config);
+    println!("{}", res.gates);
 }

src/config.rs

@@ -20,6 +20,28 @@ pub struct GridSynthConfig {
     pub verbose: bool,
     pub measure_time: bool,
     pub diophantine_data: DiophantineData,
+    pub check_solution: bool,
+}
+
+impl GridSynthConfig {
+    /// Turns on or off checking solutions at the end of the run
+    pub fn with_check_solution(self, check_solution: bool) -> Self {
+        Self {
+            check_solution,
+            ..self
+        }
+    }
+}
+
+/// The result of running the gridsynth algorithm
+pub struct GridSynthResult {
+    /// List of gates.
+    pub gates: String,
+
+    // /// The global phase factor.
+    // pub global_phase: bool,
+    /// If correctness is performed, stores the result.
+    pub is_correct: Option<bool>,
 }
 
 pub fn parse_decimal_with_exponent(input: &str) -> Option<(IBig, IBig)> {
@@ -108,5 +130,6 @@ pub fn config_from_theta_epsilon(
         verbose,
         measure_time: time,
         diophantine_data,
+        check_solution: false,
     }
 }

src/gridsynth.rs

@@ -2,7 +2,7 @@
 // Licensed under the MIT License. See LICENSE file in the project root for full license information.
 
 use crate::common::{cos_fbig, fb_with_prec, get_prec_bits, ib_to_bf_prec, sin_fbig};
-use crate::config::GridSynthConfig;
+use crate::config::{GridSynthConfig, GridSynthResult};
 use crate::diophantine::diophantine_dyadic;
 use crate::math::solve_quadratic;
 use crate::region::Ellipse;
@@ -12,6 +12,7 @@ use crate::tdgp::solve_tdgp;
 use crate::tdgp::Region;
 use crate::to_upright::to_upright_set_pair;
 use crate::unitary::DOmegaUnitary;
+use dashu_base::SquareRoot;
 use dashu_float::round::mode::{self, HalfEven};
 use dashu_float::{Context, FBig};
 use dashu_int::IBig;
@@ -20,6 +21,7 @@ use dashu_int::IBig;
 use log::debug;
 
 use nalgebra::{Matrix2, Vector2};
+use num::Complex;
 use std::cmp::Ordering;
 use std::time::{Duration, Instant};
 
@@ -42,6 +44,46 @@ fn matrix_multiply_2x2(
     result
 }
 
+fn rotation_mat(theta: &FBig<HalfEven>) -> Matrix2<Complex<FBig<HalfEven>>> {
+    let two = fb_with_prec(FBig::try_from(2.0).unwrap());
+    let theta_half = fb_with_prec(theta / &two);
+    let neg_theta_half = -fb_with_prec(theta_half);
+    let z_x: FBig<HalfEven> = fb_with_prec(cos_fbig(&neg_theta_half));
+    let z_y: FBig<HalfEven> = fb_with_prec(sin_fbig(&neg_theta_half));
+    let neg_z_y: FBig<HalfEven> = -fb_with_prec(z_y.clone());
+    let zero: FBig<HalfEven> = ib_to_bf_prec(IBig::ZERO);
+
+    Matrix2::new(
+        Complex::new(z_x.clone(), z_y.clone()),
+        Complex::new(zero.clone(), zero.clone()),
+        Complex::new(zero.clone(), zero.clone()),
+        Complex::new(z_x.clone(), neg_z_y.clone()),
+    )
+}
+
+/// Checks correctness of the synthesized circuit.
+fn check_solution(gates: &str, theta: &FBig<HalfEven>, epsilon: &FBig<HalfEven>) -> bool {
+    let expected = rotation_mat(theta);
+    let synthesized = DOmegaUnitary::from_gates(gates).to_complex_matrix();
+
+    // x = e^{-i theta / 2}
+    let x = expected[(0, 0)].clone();
+    let u = synthesized[(0, 0)].clone();
+
+    // This computes the eigenvalues of A^* A, where A = expected - synthesized.
+    // The operator norm is the square root of this.
+    let eig: FBig<HalfEven> = 2 - 2 * (&x.re * &u.re + &x.im * &u.im);
+
+    // Due to small numerical imprecisions when computing cos(\theta / 2), it may happen that we
+    // get a slightly negative value.
+    let eig = eig.max(FBig::from(0));
+
+    // Compute the norm.
+    let norm = fb_with_prec(eig.sqrt());
+
+    norm < *epsilon
+}
+
 #[derive(Debug)]
 pub struct EpsilonRegion {
     _theta: FBig<HalfEven>,
@@ -371,7 +413,7 @@ fn gridsynth(config: &mut GridSynthConfig) -> DOmegaUnitary {
     search_for_solution(&epsilon_region, &unit_disk, &transformed, config)
 }
 
-pub fn gridsynth_gates(config: &mut GridSynthConfig) -> String {
+pub fn gridsynth_gates(config: &mut GridSynthConfig) -> GridSynthResult {
     let start_total = if config.measure_time {
         Some(Instant::now())
     } else {
@@ -403,5 +445,13 @@ pub fn gridsynth_gates(config: &mut GridSynthConfig) -> String {
             );
         }
     }
-    gates
+
+    // Peform validation check, if required.
+    let is_correct = if config.check_solution {
+        Some(check_solution(&gates, &config.theta, &config.epsilon))
+    } else {
+        None
+    };
+
+    GridSynthResult { gates, is_correct }
 }

src/main.rs

@@ -47,14 +47,14 @@ fn main() {
         None
     };
 
-    let gates = gridsynth_gates(&mut args);
+    let res = gridsynth_gates(&mut args);
 
     if let Some(start_time) = start {
         let elapsed = start_time.elapsed();
         info!("Elapsed time: {:.3?}", elapsed);
     }
 
-    println!("{}", gates);
+    println!("{}", res.gates);
 }
 
 fn build_command() -> Command {
@@ -93,6 +93,11 @@ fn build_command() -> Command {
                 .short('g')
                 .action(clap::ArgAction::SetTrue),
         )
+        .arg(
+            Arg::new("check")
+                .long("check")
+                .action(clap::ArgAction::SetTrue),
+        )
 }
 
 fn parse_arguments(matches: &clap::ArgMatches) -> GridSynthConfig {
@@ -127,6 +132,7 @@ fn parse_arguments(matches: &clap::ArgMatches) -> GridSynthConfig {
         .unwrap();
     let verbose = matches.get_flag("verbose");
     let measure_time = matches.get_flag("time");
+    let check_solution = matches.get_flag("check");
 
     let seed = matches.get_one::<String>("seed").unwrap().parse().unwrap();
     let rng: StdRng = SeedableRng::seed_from_u64(seed);
@@ -142,5 +148,6 @@ fn parse_arguments(matches: &clap::ArgMatches) -> GridSynthConfig {
         verbose,
         measure_time,
         diophantine_data,
+        check_solution,
     }
 }

src/unitary.rs

@@ -2,7 +2,10 @@
 // Licensed under the MIT License. See LICENSE file in the project root for full license information.
 
 use crate::ring::DOmega;
+use dashu_float::{round::mode::HalfEven, FBig};
 use dashu_int::IBig;
+use nalgebra::Matrix2;
+use num::Complex;
 use once_cell::sync::OnceCell;
 use std::fmt::{Debug, Display, Formatter, Result};
 
@@ -55,6 +58,17 @@ impl DOmegaUnitary {
         })
     }
 
+    /// Returns 2x2 nalgebra matrix
+    pub fn to_complex_matrix(&self) -> Matrix2<Complex<FBig<HalfEven>>> {
+        let mat = self.to_matrix();
+        Matrix2::new(
+            Complex::new(mat[0][0].real().clone(), mat[0][0].imag().clone()),
+            Complex::new(mat[0][1].real().clone(), mat[0][1].imag().clone()),
+            Complex::new(mat[1][0].real().clone(), mat[1][0].imag().clone()),
+            Complex::new(mat[1][1].real().clone(), mat[1][1].imag().clone()),
+        )
+    }
+
     pub fn mul_by_t_from_left(&self) -> Self {
         Self::new(
             self.z.clone(),
@@ -161,6 +175,7 @@ impl DOmegaUnitary {
         let mut unitary = Self::identity();
         for g in gates.chars().rev() {
             unitary = match g {
+                'I' => unitary,
                 'H' => unitary.renew_denomexp(unitary.k() + 1).mul_by_h_from_left(),
                 'T' => unitary.mul_by_t_from_left(),
                 'S' => unitary.mul_by_s_from_left(),

tests/integration_test.rs

@@ -23,7 +23,7 @@ fn simple_test() {
     for (seed, expected_gates) in test_inputs {
         clear_caches();
         let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed, verbose);
-        let gates = gridsynth_gates(&mut gridsynth_config);
+        let gates = gridsynth_gates(&mut gridsynth_config).gates;
         assert_eq!(gates, expected_gates, "Test failed for seed: {}", seed);
     }
 }
@@ -42,7 +42,7 @@ fn pi_over_two_test() {
         for seed in seeds {
             clear_caches();
             let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed, verbose);
-            let gates = gridsynth_gates(&mut gridsynth_config);
+            let gates = gridsynth_gates(&mut gridsynth_config).gates;
             let expected_gates = "SWWWWWWW";
             assert_eq!(gates, expected_gates);
         }
@@ -59,8 +59,36 @@ fn slowness_from_allocation_test() {
     let seed = 1234;
     let verbose = false;
     let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed, verbose);
-    let gates = gridsynth_gates(&mut gridsynth_config);
+    let gates = gridsynth_gates(&mut gridsynth_config).gates;
     let expected_gates = "SHTSHTHTSHTSHTHTSHTHTHTSHTSHTSHTHTSHTSHTHTHTSHTSHTSHTHTHTSHTHTSHTSHTSHTSHTSHTSHTSHTHTHTHTSHTHTHTHTHTHTSHTSHTHTHTSHTSHTHTHTHTSHTSHTHTSHTSHTSHTHTSHTHTHTSHTSHTHTSHTSHTSHTHTSHTHTHTSHTHTHTSHTSHTSHTSHTHTHTHTHTSHTSHTSHTHTSHTSHTHTSHTHTSHTHTHTHTSHTSHTHTSHTSHTSHTHTSHTSHTSHTHTSHTHTHTSHTHTHTSHTSHTHTSHTSHTHTSHTSHTSHTHTHTSHTSHTSHTSHTHTHTHTSHTHTHTHTHTSHTSHTHTSHSWW";
 
     assert_eq!(gates, expected_gates);
 }
+
+#[test]
+#[serial]
+fn test_correct_decomposition() {
+    let epsilon = 1e-8;
+
+    let verbose = false;
+    let seed = 0u64;
+
+    let thetas = (0..32).map(|k| k as f64 * std::f64::consts::PI / 16.0);
+
+    for theta in thetas {
+        clear_caches();
+        let mut gridsynth_config =
+            config_from_theta_epsilon(theta, epsilon, seed, verbose).with_check_solution(true);
+
+        let res = gridsynth_gates(&mut gridsynth_config);
+
+        // not printed, unless cargo test is run with -- -no-capture
+        println!(
+            "theta = {theta}, gates = {}, correct = {:?}",
+            res.gates, res.is_correct
+        );
+
+        // Check that the check result exits and is valid.
+        assert!(res.is_correct.is_some_and(|v| v));
+    }
+}