diff --git a/examples/interface.rs b/examples/interface.rs index 1fb3926..7fef5e2 100644 --- a/examples/interface.rs +++ b/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); } diff --git a/src/config.rs b/src/config.rs index 4a20628..29f891c 100644 --- a/src/config.rs +++ b/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, } 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, } } diff --git a/src/gridsynth.rs b/src/gridsynth.rs index cd6c0fb..9723b4d 100644 --- a/src/gridsynth.rs +++ b/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) -> Matrix2>> { + 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 = fb_with_prec(cos_fbig(&neg_theta_half)); + let z_y: FBig = fb_with_prec(sin_fbig(&neg_theta_half)); + let neg_z_y: FBig = -fb_with_prec(z_y.clone()); + let zero: FBig = 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, epsilon: &FBig) -> 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 = 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, @@ -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 } } diff --git a/src/main.rs b/src/main.rs index b0bb0ea..73bea01 100644 --- a/src/main.rs +++ b/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::("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, } } diff --git a/src/unitary.rs b/src/unitary.rs index 5b3c1ca..248fafc 100644 --- a/src/unitary.rs +++ b/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>> { + 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(), diff --git a/tests/integration_test.rs b/tests/integration_test.rs index 30cce44..9c87533 100644 --- a/tests/integration_test.rs +++ b/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)); + } +}