1use num::Complex;
2use num::Zero;
3
4pub type Complex64 = Complex<f64>;
5
6pub type Qubit = Vec<Complex64>;
7
8pub type Gate = Vec<Vec<Complex64>>;
9
10pub type BinaryChars = Vec<char>;
11
12pub struct State {
13 number_of_qubits: u32,
14 pub index: usize,
15 pub amp: Complex64,
16 pub prob: f64,
17}
18
19impl State {
20 pub fn to_binary_chars(&self, qb: &[u32]) -> BinaryChars {
21 let v = to_binary_chars(self.index, self.number_of_qubits as usize);
22
23 let mut bin = vec![];
24 for i in qb {
25 bin.push(v[*i as usize]);
26 }
27
28 bin
29 }
30}
31
32impl std::fmt::Display for State {
33 fn fmt(&self, dest: &mut std::fmt::Formatter) -> std::fmt::Result {
34 let bits: String = format!("{:>0n$b}", self.index, n = self.number_of_qubits as usize);
35 write!(
36 dest,
37 "[{}]({:>+.4} {:>+.4}): {:>.4}",
38 bits, self.amp.re, self.amp.im, self.prob,
39 )
40 }
41}
42
43pub struct Q {
44 qb: Qubit,
45}
46
47impl Q {
48 pub fn new() -> Q {
49 Q { qb: vec![] }
50 }
51
52 pub fn add(&mut self, qb: Qubit) -> u32 {
53 if self.qb.is_empty() {
54 self.qb = qb;
55 return 0;
56 }
57
58 self.tensor(qb);
59 self.number_of_qubits() - 1
60 }
61
62 pub fn zero(&mut self) -> u32 {
63 self.add(vec![Complex::new(1.0, 0.0), Complex::new(0.0, 0.0)])
64 }
65
66 pub fn zeros(&mut self, n: u32) -> Vec<u32> {
67 let mut list = vec![];
68
69 for _ in 0..n {
70 list.push(self.zero());
71 }
72
73 list
74 }
75
76 pub fn zero_log2(&mut self, n: u32) -> Vec<u32> {
77 let log2n = ((n as f64).log2() as u32) + 1;
78 self.zeros(log2n)
79 }
80
81 fn tensor(&mut self, qb: Qubit) {
82 let mut v: Qubit = vec![];
83
84 for x in &self.qb {
85 for y in &qb {
86 v.push(x * y);
87 }
88 }
89
90 self.qb = v
91 }
92
93 pub fn number_of_qubits(&self) -> u32 {
94 (self.qb.len() as f64).log2() as u32
95 }
96
97 pub fn x(&mut self, qb: &[u32]) {
98 self.apply_with(x(), qb)
99 }
100
101 pub fn h(&mut self, qb: &[u32]) {
102 self.apply_with(h(), qb)
103 }
104
105 fn apply_with(&mut self, g: Gate, qb: &[u32]) {
106 let list: Vec<Gate> = gate_list(self.number_of_qubits(), g, qb);
107 let g: Gate = tensor_with(&list);
108
109 self.apply(g)
110 }
111
112 pub fn apply(&mut self, g: Gate) {
113 let mut v: Qubit = vec![];
114
115 for item in &g {
116 let mut e = Complex::new(0.0, 0.0);
117
118 for (j, _) in item.iter().enumerate() {
119 e += item[j] * self.qb[j];
120 }
121
122 v.push(e);
123 }
124
125 self.qb = v
126 }
127
128 pub fn cmodexp2(&mut self, a: u32, n: u32, r0: &[u32], r1: &[u32]) {
129 let nob = self.number_of_qubits();
130 for (i, c) in r0.iter().enumerate() {
131 self.apply(cmodexp2(nob, a, i as u32, n, *c, r1))
132 }
133 }
134
135 pub fn iqft(&mut self, qb: &[u32]) {
136 let len = qb.len();
137
138 for i in (0..len).rev() {
139 let mut k = (len - i) as i32;
140
141 for j in ((i + 1)..len).rev() {
142 let theta = -2.0 * std::f64::consts::PI / (2.0_f64.powf(k as f64));
143 self.cr(theta, qb[j], qb[i]);
144 k -= 1;
145 }
146
147 self.h(&[qb[i]]);
148 }
149 }
150
151 pub fn cr(&mut self, theta: f64, control: u32, target: u32) {
152 let n = self.number_of_qubits();
153 let g: Gate = cr(theta, n, control, target);
154 self.apply(g)
155 }
156
157 pub fn state(&self) -> Vec<State> {
158 let mut list = vec![];
159 let nob = self.number_of_qubits();
160
161 for i in 0..self.qb.len() {
162 let r = round(self.qb[i]);
163 if r.is_zero() {
164 continue;
165 }
166
167 list.push(State {
168 number_of_qubits: nob,
169 index: i,
170 amp: r,
171 prob: r.norm().powf(2.0),
172 });
173 }
174
175 list
176 }
177}
178
179impl Default for Q {
180 fn default() -> Self {
181 Self::new()
182 }
183}
184
185fn gate_list(nob: u32, g: Gate, qb: &[u32]) -> Vec<Gate> {
186 let mut list: Vec<Gate> = vec![];
187
188 for i in 0..nob {
189 let mut found = false;
190
191 for j in qb {
192 if i == *j {
193 found = true;
194 break;
195 }
196 }
197
198 if found {
199 list.push(g.to_vec());
200 continue;
201 }
202
203 list.push(id(1));
204 }
205
206 list
207}
208
209fn tensor_with(list: &[Gate]) -> Gate {
210 let mut g: Gate = list[0].to_vec();
211
212 for i in list.iter().skip(1) {
213 g = tensor(g, i.to_vec());
214 }
215
216 g
217}
218
219fn tensor(m: Gate, n: Gate) -> Gate {
220 let mut g: Gate = vec![];
221
222 for (i, _) in m.iter().enumerate() {
223 for (k, _) in n.iter().enumerate() {
224 let mut v = vec![];
225
226 for (j, _) in m[i].iter().enumerate() {
227 for (l, _) in n[k].iter().enumerate() {
228 v.push(m[i][j] * n[k][l]);
229 }
230 }
231
232 g.push(v);
233 }
234 }
235
236 g
237}
238
239pub fn x() -> Gate {
240 vec![
241 vec![Complex::new(0.0, 0.0), Complex::new(1.0, 0.0)],
242 vec![Complex::new(1.0, 0.0), Complex::new(0.0, 0.0)],
243 ]
244}
245
246pub fn h() -> Gate {
247 let e = Complex::new(1.0 / std::f64::consts::SQRT_2, 0.0);
248 vec![vec![e, e], vec![e, -1.0 * e]]
249}
250
251pub fn id(nob: u32) -> Gate {
252 let s = 1 << nob;
253 let mut g = vec![vec![Complex::new(0.0, 0.0); s]; s];
254
255 for (i, row) in g.iter_mut().enumerate() {
256 row[i] = Complex::new(1.0, 0.0);
257 }
258
259 g
260}
261
262pub fn cr(theta: f64, nob: u32, control: u32, target: u32) -> Gate {
263 let mut g: Gate = id(nob);
264 let e = Complex::new(0.0, theta).exp();
265 let mask = (1 << (nob - 1 - control)) as usize;
266
267 for (i, v) in g.iter_mut().enumerate() {
268 if (i & mask) == mask && (i & (1 << (nob - 1 - target))) != 0 {
269 v[i] = e * v[i];
270 }
271 }
272
273 g
274}
275
276pub fn cmodexp2(nob: u32, a: u32, j: u32, n: u32, control: u32, target: &[u32]) -> Gate {
277 let r1len = target.len() as u32;
278 let a2jmodn = super::number::modexp2(a, j, n);
279
280 let mut index: Vec<usize> = (0..(1 << nob)).collect();
281 for (i, idx) in index.iter_mut().enumerate() {
282 if (i >> (nob - 1 - control)) & 1 == 0 {
283 continue;
284 }
285
286 let mask = (1 << r1len) - 1;
287 let k = (i & mask) as u32;
288 if k > (n - 1) {
289 continue;
290 }
291
292 let a2jkmodn = (a2jmodn * k) % n;
294 *idx = (i >> r1len << r1len) | a2jkmodn as usize;
295 }
296
297 let id: Gate = id(nob);
298 let mut g: Gate = vec![vec![]; id.len()];
299 for (i, ii) in index.iter().enumerate() {
300 g[*ii] = id[i].to_vec();
301 }
302
303 g
304}
305
306fn round(c: Complex64) -> Complex64 {
307 let mut round = c;
308 if c.re.abs() < 1e-13 {
309 round.re = 0.0;
310 }
311
312 if c.im.abs() < 1e-13 {
313 round.im = 0.0;
314 }
315
316 round
317}
318
319fn to_binary_chars(i: usize, nob: usize) -> BinaryChars {
320 format!("{:>0n$b}", i, n = nob).chars().collect()
321}
322
323#[test]
324fn test_to_binary_chars() {
325 assert_eq!(to_binary_chars(3, 5), vec!['0', '0', '0', '1', '1']);
326 assert_eq!(to_binary_chars(7, 5), vec!['0', '0', '1', '1', '1']);
327 assert_eq!(to_binary_chars(15, 5), vec!['0', '1', '1', '1', '1']);
328 assert_eq!(to_binary_chars(31, 5), vec!['1', '1', '1', '1', '1']);
329}
330
331#[test]
332fn test_is_eigen_vector() {
333 let n = 15;
334 let a = 7;
335 let t = 3;
336
337 let mut qsim = Q::new();
338 let r0 = qsim.zeros(t);
339 let r1 = qsim.zero_log2(n);
340
341 qsim.x(&[r1[r1.len() - 1]]);
342 qsim.h(&r0);
343 qsim.cmodexp2(a, n, &r0, &r1);
344 qsim.iqft(&r0);
345
346 let mut us = std::collections::HashMap::new();
347 for state in qsim.state().iter() {
348 let m1 = state.to_binary_chars(&r1);
349 let ui: String = m1.iter().collect();
350
351 let v = match us.get(&ui) {
352 Some(vv) => state.amp + vv,
353 None => state.amp,
354 };
355
356 us.insert(ui, v);
357 }
358
359 let expected = vec![
360 ("0001", Complex::new(1.0, 0.0)),
361 ("0100", Complex::new(0.0, 0.0)),
362 ("0111", Complex::new(0.0, 0.0)),
363 ("1101", Complex::new(0.0, 0.0)),
364 ];
365
366 for (i, e) in expected {
367 let v = match us.get(i) {
368 Some(v) => *v,
369 None => panic!("{} not found", i),
370 };
371
372 assert!((v - e).re.abs() < 1e-13);
373 assert!((v - e).im.abs() < 1e-13);
374 }
375}