From 2e7e9ecdc17823281cf2d54daf4fefa1d9aec183 Mon Sep 17 00:00:00 2001 From: John Lapeyre <1969884+jlapeyre@users.noreply.github.com> Date: Sat, 1 Nov 2025 14:09:11 -0400 Subject: [PATCH] Add field up-to-phase algorithm WIP --- examples/interface.rs | 4 +- src/config.rs | 9 ++- src/gridsynth.rs | 112 ++++++++++++++++++++++++++++++++------ src/main.rs | 9 +++ src/math.rs | 10 +++- src/odgp.rs | 6 +- src/ring/d_root_two.rs | 25 +++++++++ src/ring/mod.rs | 4 +- src/ring/z_root_two.rs | 26 +++++++++ src/to_upright.rs | 7 +-- tests/integration_test.rs | 3 +- 11 files changed, 183 insertions(+), 32 deletions(-) diff --git a/examples/interface.rs b/examples/interface.rs index 901f870..b41e81f 100644 --- a/examples/interface.rs +++ b/examples/interface.rs @@ -17,7 +17,9 @@ fn main() { let theta = pi / 8.0; let epsilon = 1e-10; let seed = 1234; - let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed); + // This flag is currently ignored + let up_to_phase = false; + let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed, up_to_phase); let gates = gridsynth_gates(&mut gridsynth_config); println!("{}", gates); } diff --git a/src/config.rs b/src/config.rs index 72e9fc9..4ae70c1 100644 --- a/src/config.rs +++ b/src/config.rs @@ -20,6 +20,7 @@ pub struct GridSynthConfig { pub verbose: bool, pub measure_time: bool, pub diophantine_data: DiophantineData, + pub up_to_phase: bool, } pub fn parse_decimal_with_exponent(input: &str) -> Option<(IBig, IBig)> { @@ -69,7 +70,12 @@ pub fn parse_decimal_with_exponent(input: &str) -> Option<(IBig, IBig)> { } /// Creates the default config to easily call the code from other rust packages. -pub fn config_from_theta_epsilon(theta: f64, epsilon: f64, seed: u64) -> GridSynthConfig { +pub fn config_from_theta_epsilon( + theta: f64, + epsilon: f64, + seed: u64, + up_to_phase: bool, +) -> GridSynthConfig { let (theta_num, theta_den) = parse_decimal_with_exponent(&theta.to_string()).unwrap(); let theta = ib_to_bf_prec(theta_num) / ib_to_bf_prec(theta_den); let (epsilon_num, epsilon_den) = parse_decimal_with_exponent(&epsilon.to_string()).unwrap(); @@ -97,5 +103,6 @@ pub fn config_from_theta_epsilon(theta: f64, epsilon: f64, seed: u64) -> GridSyn verbose, measure_time: time, diophantine_data, + up_to_phase, } } diff --git a/src/gridsynth.rs b/src/gridsynth.rs index 1e764f0..147a2de 100644 --- a/src/gridsynth.rs +++ b/src/gridsynth.rs @@ -5,8 +5,9 @@ use crate::common::{cos_fbig, fb_with_prec, get_prec_bits, ib_to_bf_prec, sin_fb use crate::config::GridSynthConfig; use crate::diophantine::diophantine_dyadic; use crate::math::solve_quadratic; +use crate::math::sqrt_fbig; use crate::region::Ellipse; -use crate::ring::{DOmega, DRootTwo}; +use crate::ring::{d_root_two, z_root_two, DOmega, DRootTwo}; use crate::synthesis_of_clifford_t::decompose_domega_unitary; use crate::tdgp::solve_tdgp; use crate::tdgp::Region; @@ -39,10 +40,39 @@ fn matrix_multiply_2x2( result } +// PhaseMode::Exact synthesize gate including exact phase +// PhaseMode::Shifted synthesize gate with a fixed phase factor of `exp(i pi/8)` + +// If we don't care about phase, then it is enough to check both `U` and `exp(i pi/8) U`. +// +// To synthesize up to a phase, we run both `PhaseMode::Exact` and +// `PhaseMode::Shifted` and keep the one with lower T count. (We don't compute the best +// exact solution and then the best with the phase factor. Rather we interleave candidates +// from each to avoid doing more work than necessary. +// +// The following comments assume we are checking `exp(i pi/8) U`. +// The pair 2 ± √2 enter in some places as scale factors. +// +// omega = exp(-i pi/4) +// delta = 1 + omega +// |delta|^2 = 2 + 2cos(pi/4) = 2 + √2 +// From Lemma 9.6 in R + S, we must scale the epsilon region by +// |delta| = √(2 + √2) +// +// We scale the UnitDisk by the root-2 conjugate of |delta|: +// |delta^●| = √(2 - √2) +// See Algorithm 9.8, page 20 of R+S. +#[derive(Debug)] +pub enum PhaseMode { + Exact, // no scaling + Shifted, // do scaling +} + #[derive(Debug)] pub struct EpsilonRegion { _theta: FBig, _epsilon: FBig, + phase: PhaseMode, d: FBig, z_x: FBig, z_y: FBig, @@ -50,30 +80,56 @@ pub struct EpsilonRegion { } impl EpsilonRegion { - pub fn new(theta: FBig, epsilon: FBig) -> Self { + pub fn new(theta: FBig, epsilon: FBig, phase: PhaseMode) -> Self { let ctx: Context = Context::::new(get_prec_bits()); let one = fb_with_prec(FBig::try_from(1.0).unwrap()); let two = fb_with_prec(FBig::try_from(2.0).unwrap()); let four = fb_with_prec(FBig::try_from(4.0).unwrap()); let epsilon_squared = fb_with_prec(&epsilon * &epsilon); let half_eps_sq = fb_with_prec(&epsilon_squared / &four); - let d = fb_with_prec(ctx.sqrt((one - half_eps_sq).repr()).value()); + let delta_squared_real: FBig = fb_with_prec(z_root_two::DELTA_SQUARED.to_real()); + // The PhaseMode::Shifted arm includes a factor of |delta| = √(2 + √2) + let d = match phase { + PhaseMode::Exact => fb_with_prec(ctx.sqrt((one - half_eps_sq).repr()).value()), + PhaseMode::Shifted => fb_with_prec( + ctx.sqrt(((one - half_eps_sq) * sqrt_fbig(&delta_squared_real.clone())).repr()) + .value(), + ), + }; 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); - let epsilon_neg4: FBig = fb_with_prec(epsilon.clone().powi(IBig::from(-4))); - let epsilon_neg2: FBig = fb_with_prec(epsilon.clone().powi(IBig::from(-2))); let d1: Matrix2> = Matrix2::new(z_x.clone(), neg_z_y.clone(), z_y.clone(), z_x.clone()); + + // The PhaseMode::Shifted arms of the following two expressions divide + // the unscaled matrix by a factor |delta|^2 = 2 + √2. + let npow4 = -4; + let npow2 = -2; + let epsilon_neg4: FBig = match phase { + PhaseMode::Exact => fb_with_prec(epsilon.clone().powi(IBig::from(npow4))), + PhaseMode::Shifted => { + fb_with_prec(epsilon.clone().powi(IBig::from(npow4))) / &delta_squared_real + } + }; + + let epsilon_neg2: FBig = match phase { + PhaseMode::Exact => fb_with_prec(epsilon.clone().powi(IBig::from(npow2))), + PhaseMode::Shifted => { + fb_with_prec(epsilon.clone().powi(IBig::from(npow2))) / &delta_squared_real + } + }; + + let zero: FBig = ib_to_bf_prec(IBig::ZERO); let d2: Matrix2> = Matrix2::new( 64 * epsilon_neg4, zero.clone(), zero.clone(), 4 * epsilon_neg2, ); + let d3: Matrix2> = Matrix2::new(z_x.clone(), z_y.clone(), neg_z_y, z_x.clone()); let px = fb_with_prec(&d * &z_x); @@ -85,6 +141,7 @@ impl EpsilonRegion { Self { _theta: theta, _epsilon: epsilon, + phase, d, z_x, z_y, @@ -97,18 +154,29 @@ impl Region for EpsilonRegion { fn ellipse(&self) -> Ellipse { self.ellipse.clone() } + + // Return true if `u` is inside shaded region in figure in eq (14) in R + S + // The radius is 1 in the figure. + // For "up to phase" it is scaled by |δ|^2 = 2 + √2. fn inside(&self, u: &DOmega) -> bool { let cos_term1 = fb_with_prec(&self.z_x * u.real()); let cos_term2 = fb_with_prec(&self.z_y * u.imag()); let cos_similarity = fb_with_prec(&cos_term1 + &cos_term2); - DRootTwo::from_domega(u.conj() * u) <= DRootTwo::from_int(IBig::ONE) - && cos_similarity >= self.d + let scale = match self.phase { + PhaseMode::Exact => d_root_two::ONE, + PhaseMode::Shifted => d_root_two::DELTA_SQUARED, + }; + DRootTwo::from_domega(u.conj() * u) <= scale && cos_similarity >= self.d } fn intersect(&self, u0: &DOmega, v: &DOmega) -> Option<(FBig, FBig)> { let a = v.conj() * v; let b = 2 * (v.conj() * u0); - let c = u0.conj() * u0 - DOmega::from_int(IBig::ONE); + let scale = match self.phase { + PhaseMode::Exact => DOmega::from_int(IBig::ONE), + PhaseMode::Shifted => DOmega::from_zroottwo(&z_root_two::DELTA_SQUARED), + }; + let c = u0.conj() * u0 - scale; let vz_term1 = fb_with_prec(&self.z_x * v.real()); let vz_term2 = fb_with_prec(&self.z_y * v.imag()); let vz = fb_with_prec(&vz_term1 + &vz_term2); @@ -136,11 +204,12 @@ impl Region for EpsilonRegion { #[derive(Debug)] pub struct UnitDisk { + phase: PhaseMode, ellipse: Ellipse, } impl UnitDisk { - pub fn new() -> Self { + pub fn new(phase: PhaseMode) -> Self { let ellipse = Ellipse::from( ib_to_bf_prec(IBig::ONE), ib_to_bf_prec(IBig::ZERO), @@ -149,7 +218,7 @@ impl UnitDisk { ib_to_bf_prec(IBig::ZERO), ib_to_bf_prec(IBig::ZERO), ); - Self { ellipse } + Self { phase, ellipse } } pub fn ellipse(&self) -> &Ellipse { @@ -157,24 +226,35 @@ impl UnitDisk { } } +// We might not want this impl Default for UnitDisk { fn default() -> Self { - Self::new() + Self::new(PhaseMode::Exact) } } +// For PhaseMode::Shifted, we scale the unit disk by the root-2 conjugate +// of |delta|^2. (This is LAMBDA_M) impl Region for UnitDisk { fn ellipse(&self) -> Ellipse { self.ellipse.clone() } fn inside(&self, u: &DOmega) -> bool { - DRootTwo::from_domega(u.conj() * u) <= DRootTwo::from_int(IBig::ONE) + let scale = match self.phase { + PhaseMode::Exact => d_root_two::ONE, + PhaseMode::Shifted => d_root_two::DELTA_SQUARED_M, + }; + DRootTwo::from_domega(u.conj() * u) <= scale } fn intersect(&self, u0: &DOmega, v: &DOmega) -> Option<(FBig, FBig)> { let a = v.conj() * v; let b = 2 * (v.conj() * u0); - let c = u0.conj() * u0 - IBig::ONE; + let shift_ = match self.phase { + PhaseMode::Exact => DOmega::from_zroottwo(&z_root_two::ONE), + PhaseMode::Shifted => DOmega::from_zroottwo(&z_root_two::DELTA_SQUARED_M), + }; + let c = u0.conj() * u0 - shift_; solve_quadratic(a.real(), b.real(), c.real()) } } @@ -256,8 +336,8 @@ fn setup_regions_and_transform( crate::region::Rectangle, ), ) { - let epsilon_region = EpsilonRegion::new(theta, epsilon); - let unit_disk = UnitDisk::new(); + let epsilon_region = EpsilonRegion::new(theta, epsilon, PhaseMode::Exact); + let unit_disk = UnitDisk::new(PhaseMode::Exact); let start_upright = if measure_time { Some(Instant::now()) diff --git a/src/main.rs b/src/main.rs index b0bb0ea..9374b81 100644 --- a/src/main.rs +++ b/src/main.rs @@ -93,6 +93,13 @@ fn build_command() -> Command { .short('g') .action(clap::ArgAction::SetTrue), ) + // We use `phase` rather than, say, `up_to_phase` to agree with original gridsynth. + .arg( + Arg::new("phase") + .long("phase") + .short('p') + .action(clap::ArgAction::SetTrue), + ) } fn parse_arguments(matches: &clap::ArgMatches) -> GridSynthConfig { @@ -127,6 +134,7 @@ fn parse_arguments(matches: &clap::ArgMatches) -> GridSynthConfig { .unwrap(); let verbose = matches.get_flag("verbose"); let measure_time = matches.get_flag("time"); + let up_to_phase = matches.get_flag("phase"); let seed = matches.get_one::("seed").unwrap().parse().unwrap(); let rng: StdRng = SeedableRng::seed_from_u64(seed); @@ -142,5 +150,6 @@ fn parse_arguments(matches: &clap::ArgMatches) -> GridSynthConfig { verbose, measure_time, diophantine_data, + up_to_phase, } } diff --git a/src/math.rs b/src/math.rs index cb3955a..68703da 100644 --- a/src/math.rs +++ b/src/math.rs @@ -1,7 +1,7 @@ // Copyright (c) 2024-2025 Shun Yamamoto and Nobuyuki Yoshioka, and IBM // Licensed under the MIT License. See LICENSE file in the project root for full license information. -use crate::common::{get_prec_bits, ib_to_bf_prec}; +use crate::common::{fb_with_prec, get_prec_bits, ib_to_bf_prec}; use dashu_float::round::mode; use dashu_float::Context; use dashu_float::{round::mode::HalfEven, FBig}; @@ -16,6 +16,14 @@ pub fn sqrt2() -> FBig { a } +// This may be wasteful because of the allocation +pub fn sqrt_fbig(x: &FBig) -> FBig { + let ctx: Context = Context::::new(get_prec_bits()); + let x = fb_with_prec(x.clone()); + let sx: FBig = ctx.sqrt(x.repr()).value(); + sx +} + pub fn ntz(n: &IBig) -> i64 { match n.trailing_zeros() { Some(k) => k as i64, diff --git a/src/odgp.rs b/src/odgp.rs index 30a8b38..3cdd75e 100644 --- a/src/odgp.rs +++ b/src/odgp.rs @@ -4,16 +4,12 @@ use crate::common::ib_to_bf_prec; use crate::math::{floorlog, pow_sqrt2, sqrt2}; use crate::region::Interval; +use crate::ring::z_root_two::LAMBDA; use crate::ring::{DRootTwo, ZRootTwo}; use dashu_float::round::mode::HalfEven; use dashu_float::FBig; use dashu_int::IBig; -const LAMBDA: ZRootTwo = ZRootTwo { - a: IBig::ONE, - b: IBig::ONE, -}; - pub fn solve_odgp(i: Interval, j: Interval) -> impl Iterator { // Can't return two different iterator types. So we can't do this check. // I checked with dbg! to confirm that omitting this check is ok: diff --git a/src/ring/d_root_two.rs b/src/ring/d_root_two.rs index d5d5162..7279d93 100644 --- a/src/ring/d_root_two.rs +++ b/src/ring/d_root_two.rs @@ -1,6 +1,7 @@ // Copyright (c) 2024-2025 Shun Yamamoto and Nobuyuki Yoshioka, and IBM // Licensed under the MIT License. See LICENSE file in the project root for full license information. +use dashu_base::Sign; use dashu_float::round::mode::HalfEven; use dashu_float::FBig; use dashu_int::IBig; @@ -17,6 +18,30 @@ pub struct DRootTwo { pub(crate) k: i64, } +pub const ONE: DRootTwo = DRootTwo { + alpha: ZRootTwo { + a: IBig::ONE, + b: IBig::ZERO, + }, + k: 0, +}; + +pub const DELTA_SQUARED: DRootTwo = DRootTwo { + alpha: ZRootTwo { + a: IBig::from_parts_const(Sign::Positive, 2), + b: IBig::ONE, + }, + k: 0, +}; + +pub const DELTA_SQUARED_M: DRootTwo = DRootTwo { + alpha: ZRootTwo { + a: IBig::from_parts_const(Sign::Positive, 2), + b: IBig::NEG_ONE, + }, + k: 0, +}; + impl DRootTwo { pub fn new(alpha: ZRootTwo, k: i64) -> Self { Self { alpha, k } diff --git a/src/ring/mod.rs b/src/ring/mod.rs index 6d45bce..8f8dfca 100644 --- a/src/ring/mod.rs +++ b/src/ring/mod.rs @@ -2,9 +2,9 @@ // Licensed under the MIT License. See LICENSE file in the project root for full license information. mod d_omega; -mod d_root_two; +pub mod d_root_two; mod z_omega; -mod z_root_two; +pub mod z_root_two; pub use d_omega::DOmega; pub use d_root_two::DRootTwo; diff --git a/src/ring/z_root_two.rs b/src/ring/z_root_two.rs index 24567a1..7b64481 100644 --- a/src/ring/z_root_two.rs +++ b/src/ring/z_root_two.rs @@ -4,6 +4,7 @@ use crate::common::ib_to_bf_prec; use crate::math::{floorsqrt, rounddiv, sign, sqrt2}; use crate::ring::ZOmega; +use dashu_base::Sign; use dashu_float::round::mode::HalfEven; use dashu_float::FBig; use dashu_int::IBig; @@ -18,6 +19,31 @@ pub struct ZRootTwo { pub(crate) b: IBig, } +pub const ONE: ZRootTwo = ZRootTwo { + a: IBig::ONE, + b: IBig::ZERO, +}; + +pub const DELTA_SQUARED: ZRootTwo = ZRootTwo { + a: IBig::from_parts_const(Sign::Positive, 2), + b: IBig::ONE, +}; + +// Absolute value squared of root two conjugate of δ. +pub const DELTA_SQUARED_M: ZRootTwo = ZRootTwo { + a: IBig::from_parts_const(Sign::Positive, 2), + b: IBig::NEG_ONE, +}; + +// See Definition 3.5 on pg 3 of R+S +// Careful! There is a different, unrelated defintion of lambda in +// the discussion in Definition 9.1 on page 19 of R+S. It is in fact +// the fixed phase factor used in the up-to-phase algorithm. +pub const LAMBDA: ZRootTwo = ZRootTwo { + a: IBig::ONE, + b: IBig::ONE, +}; + impl ZRootTwo { pub fn new(a: IBig, b: IBig) -> Self { Self { a, b } diff --git a/src/to_upright.rs b/src/to_upright.rs index 6b2043c..c05c17a 100644 --- a/src/to_upright.rs +++ b/src/to_upright.rs @@ -8,14 +8,11 @@ use crate::common::ib_to_bf_prec; use crate::grid_op::{EllipsePair, GridOp}; use crate::math::{floorsqrt, log}; use crate::region::{Ellipse, Rectangle}; -use crate::ring::{ZOmega, ZRootTwo}; +use crate::ring::z_root_two::LAMBDA; +use crate::ring::ZOmega; use crate::tdgp::Region; use dashu_int::IBig; use num_traits::Pow; -const LAMBDA: ZRootTwo = ZRootTwo { - a: IBig::ONE, - b: IBig::ONE, -}; fn reduction( ellipse_pair: EllipsePair, diff --git a/tests/integration_test.rs b/tests/integration_test.rs index 4ead8b8..e8ace5d 100644 --- a/tests/integration_test.rs +++ b/tests/integration_test.rs @@ -10,7 +10,8 @@ fn simple_test() { let theta = pi / 8.0; let epsilon = 1e-10; let seed = 1234; - let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed); + let up_to_phase = false; + let mut gridsynth_config = config_from_theta_epsilon(theta, epsilon, seed, up_to_phase); let gates = gridsynth_gates(&mut gridsynth_config); let expected_gates = "HTHTSHTSHTHTSHTHTSHTHTHTSHTSHTHTHTHTHTHTSHTSHTHTSHTSHTSHTSHTHTSHTSHTSHTHTHTHTHTHTSHTSHTHTSHTSHTSHTHTHTSHTSHTSHTSHTSHTSHTSHTHTHTHTHTSHTSHTSHTSHTSHTSHTHTHTHTHTSHTHTSHTHTHTSHTSHTSHTHTSHTSHTHTSHTHTSHTSHTHTSHTHTHTSHTSHTSHTSHTHTHTHTSHTHTHTSHTHTSHTHTHTSHTHTSHTHTSHTXSSWWW"; assert_eq!(gates, expected_gates);