GitVita
New repositoryImport repositoryNew issue
Your repositories Your dashboard Sign in

microsoft/qdk

Public

mirrored fromhttps://github.com/microsoft/qdkAvailable

Watch0Fork0Star0
CodeCommitsIssuesPull requestsActionsInsightsSecurity
v1.9.0

Branches

Tags

  • No tags available.
0Branches0Tags
Go to file
Add file
Code

Clone

HTTPS

Download ZIP

qdk/compiler/qsc_eval/src

compiler/qsc_eval/src/state.rs

612lines · modecode

RawDownload
Latest commit unavailable.
1// Copyright (c) Microsoft Corporation.
2// Licensed under the MIT License.
3
4#[cfg(test)]
5mod tests;
6
7use num_bigint::BigUint;
8use num_complex::{Complex, Complex64};
9use std::{f64::consts::FRAC_1_SQRT_2, fmt::Write};
10
11#[must_use]
12pub fn format_state_id(id: &BigUint, qubit_count: usize) -> String {
13 format!("|{}⟩", fmt_basis_state_label(id, qubit_count))
14}
15
16#[must_use]
17pub fn get_phase(c: &Complex<f64>) -> f64 {
18 f64::atan2(c.im, c.re)
19}
20
21#[must_use]
22pub fn fmt_complex(c: &Complex<f64>) -> String {
23 // Format -0 as 0
24 // Also using Unicode Minus Sign instead of ASCII Hyphen-Minus
25 // and Unicode Mathematical Italic Small I instead of ASCII i.
26 format!(
27 "{}{:.4}{}{:.4}𝑖",
28 if c.re <= -0.00005 { "−" } else { "" },
29 c.re.abs(),
30 if c.im <= -0.00005 { "−" } else { "+" },
31 c.im.abs()
32 )
33}
34
35#[must_use]
36pub fn fmt_basis_state_label(id: &BigUint, qubit_count: usize) -> String {
37 // This will generate a bit string that shows the qubits in the order
38 // of allocation, left to right.
39 format!("{:0>qubit_count$}", id.to_str_radix(2))
40}
41
42#[must_use]
43fn is_significant(x: f64) -> bool {
44 x.abs() > 1e-9
45}
46
47#[must_use]
48fn is_fractional_part_significant(x: f64) -> bool {
49 is_significant(x - x.round())
50}
51
52/// Represents a non-zero rational number in the form ``numerator/denominator``
53/// Sign of the number is separated for easier composition and rendering
54#[derive(Copy, Clone, Debug)]
55struct RationalNumber {
56 sign: i64, // 1 if the number is positive, -1 if negative
57 numerator: i64, // Positive numerator
58 denominator: i64, // Positive denominator
59}
60
61impl RationalNumber {
62 fn new(numerator: i64, denominator: i64) -> Self {
63 Self {
64 sign: numerator.signum() * denominator.signum(),
65 numerator: numerator.abs(),
66 denominator: denominator.abs(),
67 }
68 }
69
70 fn abs(&self) -> Self {
71 Self {
72 sign: self.sign.abs(),
73 numerator: self.numerator,
74 denominator: self.denominator,
75 }
76 }
77
78 const DENOMINATORS: [i64; 36] = [
79 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 12, 14, 15, 16, 18, 20, 21, 24, 25, 27, 28, 30, 32, 35, 36,
80 40, 42, 45, 48, 49, 50, 54, 56, 60, 63, 64,
81 ];
82
83 /// Try to recognize a float number as a "nice" rational.
84 fn recognize(x: f64) -> Option<Self> {
85 for denominator in Self::DENOMINATORS {
86 #[allow(clippy::cast_precision_loss)] // We only use fixes set of denominators
87 let numerator: f64 = x * (denominator as f64);
88 if numerator.abs() <= 100.0 && !is_fractional_part_significant(numerator) {
89 #[allow(clippy::cast_possible_truncation)] // We only allow small numerators
90 let rounded_numerator: i64 = numerator.round() as i64;
91 return Some(Self::new(rounded_numerator, denominator));
92 }
93 }
94 None
95 }
96}
97
98/// Represents a non-zero algebraic number in the form ``fraction·√root``,
99/// Sign of the number is separated for easier composition and rendering
100#[derive(Debug)]
101struct AlgebraicNumber {
102 sign: i64, // 1 if the number is positive, -1 if negative
103 fraction: RationalNumber, // Positive rational number numerator/denominator
104 root: i64, // Component under square root
105}
106
107impl AlgebraicNumber {
108 fn new(fraction: &RationalNumber, root: i64) -> Self {
109 Self {
110 sign: fraction.sign,
111 fraction: fraction.abs(),
112 root,
113 }
114 }
115
116 const ROOTS: [i64; 6] = [1, 2, 3, 5, 6, 7];
117
118 /// Try to recognize a float number as a "nice" algebraic.
119 fn recognize(x: f64) -> Option<Self> {
120 for root in Self::ROOTS {
121 #[allow(clippy::cast_precision_loss)] // We only use fixes set of roots
122 let divided_by_root: f64 = x / (root as f64).sqrt();
123 if let Some(fraction) = RationalNumber::recognize(divided_by_root) {
124 return Some(Self::new(&fraction, root));
125 }
126 }
127 None
128 }
129}
130
131/// Represents a non-zero decimal number as an ``f64`` floating point value.
132/// Sign of the number is separated for easier composition and rendering
133#[derive(Debug)]
134struct DecimalNumber {
135 sign: i64, // 1 if the number is positive, -1 if negative
136 value: f64, // Positive floating point value
137}
138
139impl DecimalNumber {
140 fn new(value: f64) -> Self {
141 Self {
142 sign: if value >= 0.0 { 1 } else { -1 },
143 value: value.abs(),
144 }
145 }
146
147 // Tries to recognize a decimal number and always succeeds.
148 fn recognize(x: f64) -> Self {
149 Self::new(x)
150 }
151}
152
153/// Represents a real number, which can be algebraic, decimal or zero.
154#[derive(Debug)]
155enum RealNumber {
156 Algebraic(AlgebraicNumber),
157 Decimal(DecimalNumber),
158 Zero,
159}
160
161impl RealNumber {
162 fn sign(&self) -> i64 {
163 match self {
164 Self::Algebraic(algebraic) => algebraic.sign,
165 Self::Decimal(decimal) => decimal.sign,
166 Self::Zero => 0,
167 }
168 }
169
170 fn negate(&self) -> Self {
171 match self {
172 Self::Algebraic(algebraic) => Self::Algebraic(AlgebraicNumber {
173 sign: -algebraic.sign,
174 fraction: algebraic.fraction,
175 root: algebraic.root,
176 }),
177 Self::Decimal(decimal) => Self::Decimal(DecimalNumber {
178 sign: -decimal.sign,
179 value: decimal.value,
180 }),
181 Self::Zero => Self::Zero,
182 }
183 }
184
185 fn abs(&self) -> Self {
186 match self {
187 Self::Algebraic(algebraic) => Self::Algebraic(AlgebraicNumber {
188 sign: 1,
189 fraction: algebraic.fraction,
190 root: algebraic.root,
191 }),
192 Self::Decimal(decimal) => Self::Decimal(DecimalNumber {
193 sign: 1,
194 value: decimal.value,
195 }),
196 Self::Zero => Self::Zero,
197 }
198 }
199
200 /// Try to recognize a real number as zero, algebraic, or decimal if all else fails.
201 fn recognize(x: f64) -> Self {
202 if !is_significant(x) {
203 Self::Zero
204 } else if let Some(algebraic_number) = AlgebraicNumber::recognize(x) {
205 Self::Algebraic(algebraic_number)
206 } else {
207 Self::Decimal(DecimalNumber::recognize(x))
208 }
209 }
210}
211
212/// Represents a non-zero complex numbers in the polar form: ``magnitude·𝒆^(𝝅·𝒊·phase_multiplier)``
213/// Sign of the number is separated for easier composition and rendering
214#[derive(Debug)]
215struct PolarForm {
216 sign: i64, // For this form the sign is always 1
217 magnitude: AlgebraicNumber, // magnitude of the number
218 phase_multiplier: RationalNumber, // to be multiplied by 𝝅·𝒊 to get phase
219}
220
221impl PolarForm {
222 fn new(magnitude: &RationalNumber, pi_num: i64, pi_den: i64) -> Self {
223 Self {
224 sign: 1,
225 magnitude: AlgebraicNumber::new(magnitude, 1),
226 phase_multiplier: RationalNumber::new(pi_num, pi_den),
227 }
228 }
229
230 const PI_FRACTIONS: [(i64, i64); 16] = [
231 (1, 3),
232 (2, 3),
233 (1, 4),
234 (3, 4),
235 (1, 8),
236 (3, 8),
237 (5, 8),
238 (7, 8),
239 (1, 16),
240 (3, 16),
241 (5, 16),
242 (7, 16),
243 (9, 16),
244 (11, 16),
245 (13, 16),
246 (15, 16),
247 ];
248
249 /// Try to recognize a complex number and represent it in the polar form.
250 fn recognize(re: f64, im: f64) -> Option<Self> {
251 if !is_significant(re.abs()) && !is_significant(im.abs()) {
252 // 0 is better represented in Cartesian form, not polar.
253 return None;
254 }
255 for (pi_num, pi_den) in Self::PI_FRACTIONS {
256 #[allow(clippy::cast_precision_loss)] // We only use fixes set of fractions
257 let angle: f64 = std::f64::consts::PI * (pi_num as f64) / (pi_den as f64);
258 let sin: f64 = angle.sin();
259 let cos: f64 = angle.cos();
260 if is_significant(re / cos - im / sin) {
261 continue;
262 }
263 // We recognized the angle. Now try to recognize magnitude.
264 // It's OK to take abs value as we are only interested in magnitude
265 if let Some(magnitude) = RationalNumber::recognize((re / cos).abs()) {
266 return Some(Self::new(
267 &magnitude,
268 if im >= 0.0 { pi_num } else { pi_num - pi_den },
269 pi_den,
270 ));
271 }
272 }
273 None
274 }
275}
276
277/// Represents a non-zero complex number in the Cartesian form: ``real_part+𝒊·imaginary_part``
278/// Sign of the number is separated for easier composition and rendering
279#[derive(Debug)]
280struct CartesianForm {
281 sign: i64, // 1 the common sign is a "+", -1 is "-", 0 means the number is 0
282 real_part: RealNumber, // Real part
283 imaginary_part: RealNumber, // Imaginary part
284}
285
286impl CartesianForm {
287 fn new(real_part: RealNumber, imaginary_part: RealNumber) -> Self {
288 let sign = real_part.sign();
289 match sign {
290 0 => Self {
291 sign: imaginary_part.sign(),
292 real_part: RealNumber::Zero,
293 imaginary_part: imaginary_part.abs(),
294 },
295 1.. => Self {
296 sign,
297 real_part,
298 imaginary_part,
299 },
300 _ => Self {
301 sign,
302 real_part: real_part.abs(),
303 imaginary_part: imaginary_part.negate(),
304 },
305 }
306 }
307
308 /// Try to recognize a complex number and represent it in the Cartesian form.
309 fn recognize(re: f64, im: f64) -> Self {
310 Self::new(RealNumber::recognize(re), RealNumber::recognize(im))
311 }
312}
313
314/// Represents a non-zero complex number which can be in either polar or Cartesian form.
315#[derive(Debug)]
316enum ComplexNumber {
317 Polar(PolarForm),
318 Cartesian(CartesianForm),
319}
320
321impl ComplexNumber {
322 fn recognize(re: f64, im: f64) -> Self {
323 if let Some(exponent) = PolarForm::recognize(re, im) {
324 Self::Polar(exponent)
325 } else {
326 Self::Cartesian(CartesianForm::recognize(re, im))
327 }
328 }
329}
330
331/// Represents one term of a quantum state, which corresponds to one basis vector.
332struct Term {
333 basis_vector: BigUint,
334 coordinate: ComplexNumber,
335}
336
337fn get_terms_for_state(state: &Vec<(BigUint, Complex64)>) -> Vec<Term> {
338 let mut result: Vec<Term> = Vec::with_capacity(state.len());
339 for (basis_vector, coefficient) in state {
340 result.push(Term {
341 basis_vector: basis_vector.clone(),
342 coordinate: ComplexNumber::recognize(coefficient.re, coefficient.im),
343 });
344 }
345 result
346}
347
348/// Get the state represented as a formula in the LaTeX format if possible.
349/// `None` is returned if the resulting formula is not nice, i.e.
350/// if the formula consists of more than 16 terms or if more than two coefficients are not recognized.
351#[must_use]
352pub fn get_state_latex(state: &Vec<(BigUint, Complex64)>, qubit_count: usize) -> Option<String> {
353 if state.len() > 16 {
354 return None;
355 }
356
357 let terms: Vec<Term> = get_terms_for_state(state);
358 let mut bad_term_count: i64 = 0;
359 for term in &terms {
360 if let ComplexNumber::Cartesian(cartesian) = &term.coordinate {
361 if let RealNumber::Decimal(_) = &cartesian.real_part {
362 bad_term_count += 1;
363 }
364 if let RealNumber::Decimal(_) = &cartesian.imaginary_part {
365 bad_term_count += 1;
366 };
367 if bad_term_count > 2 {
368 return None;
369 }
370 }
371 }
372
373 let mut latex: String = String::with_capacity(200);
374 latex.push_str("$|\\psi\\rangle = ");
375 let mut is_first: bool = true;
376 for term in terms {
377 write_latex_for_term(&mut latex, &term, !is_first);
378 let basis_label = fmt_basis_state_label(&term.basis_vector, qubit_count);
379 write!(latex, "|{basis_label}\\rangle").expect("Expected to write basis label.");
380 is_first = false;
381 }
382 latex.push('$');
383 latex.shrink_to_fit();
384
385 Some(latex)
386}
387
388fn is_close_enough(val: &Complex64, target: &Complex64) -> bool {
389 (val.re - target.re).abs() < 1e-9 && (val.im - target.im).abs() < 1e-9
390}
391
392// Quick and dirty matching for the most common matrix elements we care about rendering
393// LaTeX for, e.g. 1/sqrt(2), -i/sqrt(2), etc.
394// Anything not in this list gets a standard rendering.
395fn get_latex_for_simple_term(val: &Complex64) -> Option<String> {
396 if is_close_enough(val, &Complex64::new(FRAC_1_SQRT_2, 0.0)) {
397 return Some("\\frac{1}{\\sqrt{2}}".to_string());
398 }
399 if is_close_enough(val, &Complex64::new(0.0, FRAC_1_SQRT_2)) {
400 return Some("\\frac{i}{\\sqrt{2}}".to_string());
401 }
402 if is_close_enough(val, &Complex64::new(-FRAC_1_SQRT_2, 0.0)) {
403 return Some("-\\frac{1}{\\sqrt{2}}".to_string());
404 }
405 if is_close_enough(val, &Complex64::new(0.0, -FRAC_1_SQRT_2)) {
406 return Some("-\\frac{i}{\\sqrt{2}}".to_string());
407 }
408 if is_close_enough(val, &Complex64::new(0.0, 0.5)) {
409 return Some("\\frac{i}{2}".to_string());
410 }
411 if is_close_enough(val, &Complex64::new(0.0, -0.5)) {
412 return Some("-\\frac{i}{2}".to_string());
413 }
414 if is_close_enough(val, &Complex64::new(0.5, 0.5)) {
415 return Some("\\frac{1}{2} + \\frac{i}{2}".to_string());
416 }
417 if is_close_enough(val, &Complex64::new(-0.5, -0.5)) {
418 return Some("-\\frac{1}{2} - \\frac{i}{2}".to_string());
419 }
420 if is_close_enough(val, &Complex64::new(-0.5, 0.5)) {
421 return Some("-\\frac{1}{2} + \\frac{i}{2}".to_string());
422 }
423 if is_close_enough(val, &Complex64::new(0.5, -0.5)) {
424 return Some("\\frac{1}{2} - \\frac{i}{2}".to_string());
425 }
426 None
427}
428
429#[must_use]
430pub fn get_matrix_latex(matrix: &Vec<Vec<Complex64>>) -> String {
431 let mut latex: String = String::with_capacity(500);
432 latex.push_str("$ \\begin{bmatrix} ");
433 for row in matrix {
434 let mut is_first: bool = true;
435 for element in row {
436 if !is_first {
437 latex.push_str(" & ");
438 }
439 is_first = false;
440
441 if let Some(simple_latex) = get_latex_for_simple_term(element) {
442 latex.push_str(&simple_latex);
443 continue;
444 }
445
446 let cpl = ComplexNumber::recognize(element.re, element.im);
447 write_latex_for_complex_number(&mut latex, &cpl);
448 }
449 latex.push_str(" \\\\ ");
450 }
451 latex.push_str("\\end{bmatrix} $");
452 latex.shrink_to_fit();
453 latex
454}
455
456/// Write latex for a standalone complex number
457/// '-', 0 and 1 are always rendered, but '+' is not.
458fn write_latex_for_complex_number(latex: &mut String, number: &ComplexNumber) {
459 match number {
460 ComplexNumber::Cartesian(cartesian_form) => {
461 write_latex_for_cartesian_form(latex, cartesian_form, false, true);
462 }
463 ComplexNumber::Polar(polar_form) => {
464 write_latex_for_polar_form(latex, polar_form, false);
465 }
466 }
467}
468
469/// Write latex for one term of quantum state.
470/// Latex is rendered for coefficient only (not for basis vector).
471/// + is rendered only if ``render_plus`` is true.
472fn write_latex_for_term(latex: &mut String, term: &Term, render_plus: bool) {
473 match &term.coordinate {
474 ComplexNumber::Cartesian(cartesian_form) => {
475 write_latex_for_cartesian_form(latex, cartesian_form, render_plus, false);
476 }
477 ComplexNumber::Polar(polar_form) => {
478 write_latex_for_polar_form(latex, polar_form, render_plus);
479 }
480 }
481}
482
483/// Write latex for polar form of a complex number.
484/// The sign is always + and is rendered only if ``render_plus`` is true.
485fn write_latex_for_polar_form(latex: &mut String, polar_form: &PolarForm, render_plus: bool) {
486 if polar_form.sign < 0 {
487 latex.push('-');
488 } else if render_plus {
489 latex.push('+');
490 }
491 write_latex_for_algebraic_number(latex, &polar_form.magnitude, false);
492 latex.push_str(" e^{");
493 if polar_form.phase_multiplier.sign < 0 {
494 latex.push('-');
495 }
496 let pi_num: i64 = polar_form.phase_multiplier.numerator;
497 if pi_num != 1 {
498 write!(latex, "{pi_num}").expect("Expected to write phase numerator.");
499 }
500 latex.push_str(" i \\pi");
501 let pi_den = polar_form.phase_multiplier.denominator;
502 if pi_den != 1 {
503 write!(latex, " / {pi_den}").expect("expected to write phase denominator.");
504 }
505 latex.push('}');
506}
507
508/// Write latex for cartesian form of a complex number.
509/// Common + is rendered only if ``render_plus`` is true.
510/// Brackets are used if both real and imaginary parts are present.
511/// If only one part is present, its sign is used as common.
512/// If both components are present, real part sign is used as common.
513/// 1 is rendered if ``render_one`` is true
514/// + is rendered if ``render_plus`` is true.
515fn write_latex_for_cartesian_form(
516 latex: &mut String,
517 cartesian_form: &CartesianForm,
518 render_plus: bool,
519 render_one: bool,
520) {
521 if cartesian_form.sign < 0 {
522 latex.push('-');
523 } else if render_plus {
524 latex.push('+');
525 }
526 if let RealNumber::Zero = cartesian_form.real_part {
527 if let RealNumber::Zero = cartesian_form.imaginary_part {
528 latex.push('0');
529 } else {
530 // Only imaginary part present
531 write_latex_for_real_number(latex, &cartesian_form.imaginary_part, false);
532 latex.push('i');
533 }
534 } else if let RealNumber::Zero = cartesian_form.imaginary_part {
535 // Only real part present
536 write_latex_for_real_number(latex, &cartesian_form.real_part, render_one);
537 } else {
538 // Both real and imaginary parts present
539 latex.push_str("\\left( ");
540 write_latex_for_real_number(latex, &cartesian_form.real_part, true);
541 latex.push(if cartesian_form.imaginary_part.sign() < 0 {
542 '-'
543 } else {
544 '+'
545 });
546 write_latex_for_real_number(latex, &cartesian_form.imaginary_part, false);
547 latex.push_str("i \\right)");
548 }
549}
550
551/// Write latex for real number. Note that the sign is not rendered.
552/// 1 is only rendered if ``render_one`` is true.
553fn write_latex_for_real_number(latex: &mut String, number: &RealNumber, render_one: bool) {
554 match number {
555 RealNumber::Algebraic(algebraic_number) => {
556 write_latex_for_algebraic_number(latex, algebraic_number, render_one);
557 }
558 RealNumber::Decimal(decimal_number) => {
559 write_latex_for_decimal_number(latex, decimal_number, render_one);
560 }
561 RealNumber::Zero => {
562 latex.push('0');
563 }
564 }
565}
566
567/// Write latex for decimal number. Note that the sign is not rendered.
568/// 1 is only rendered if ``render_one`` is true.
569fn write_latex_for_decimal_number(latex: &mut String, number: &DecimalNumber, render_one: bool) {
570 if render_one || is_significant(number.value - 1.0) {
571 // Using round() instead of neater string formatting({:.4})
572 // because we do not want trailing zeros (we need 0.5 and not 0.5000)
573 write!(latex, "{}", (number.value * 10000.0).round() / 10000.0)
574 .expect("Expected to write decimal value.");
575 }
576}
577
578/// Write latex for algebraic number. Note that the sign is not rendered.
579/// 1 is only rendered if ``render_one`` is true.
580fn write_latex_for_algebraic_number(
581 latex: &mut String,
582 number: &AlgebraicNumber,
583 render_one: bool,
584) {
585 let actually_needs_1: bool = render_one || number.fraction.denominator != 1;
586 if number.fraction.denominator != 1 {
587 latex.push_str("\\frac{");
588 }
589 if number.root == 1 {
590 if number.fraction.numerator == 1 {
591 if actually_needs_1 {
592 latex.push('1');
593 }
594 } else {
595 write!(latex, "{}", number.fraction.numerator).expect("Expected to write numerator.");
596 }
597 } else if number.fraction.numerator == 1 {
598 write!(latex, "\\sqrt{{{}}}", number.root)
599 .expect("Expected to write square root expression.");
600 } else {
601 write!(
602 latex,
603 "{} \\sqrt{{{}}}",
604 number.fraction.numerator, number.root
605 )
606 .expect("expected to write numerator and square root expression.");
607 }
608 if number.fraction.denominator != 1 {
609 write!(latex, "}}{{{}}}", number.fraction.denominator)
610 .expect("Expected to write denominator.");
611 }
612}
613