Quantum phase, measurement, and information loss / exact Python
Gaussian rational matrix arithmetic, three-setting tomography, partial dephasing, operation order, coherent records, and packet refinement, amplitude damping, modulus-12 arithmetic interference, and three-qubit syndrome recovery. Python 3.10 or later; standard library only.
Read the worked derivation · Download Python · Source and results ZIP · Revision and source identities
Exact generated records: measurement-settings.csv · partial-dephasing.csv · results.json · measurement-comparison.svg · amplitude-damping.csv · arithmetic-phase.csv · syndrome-recovery.csv
Run python -B quantum_phase_measurement.py. Add --output quantum-results to reproduce the records.
1 """Exact finite quantum phase, measurement, and information-loss example.
2
3 Python 3.10+, standard library only. All calculations use Q(i); no complex
4 floating-point arithmetic, sampling, simulator, or Lean replay is involved.
5 Run with --output DIRECTORY to reproduce the CSV, JSON, and SVG records.
6 """
7 from __future__ import annotations
8
9 import argparse
10 import csv
11 from dataclasses import dataclass
12 from fractions import Fraction as Q
13 import json
14 from pathlib import Path
15
16
17 @dataclass(frozen=True)
18 class C:
19 """A Gaussian rational, stored as exact real and imaginary fractions."""
20 re: Q = Q(0)
21 im: Q = Q(0)
22
23 def __post_init__(self):
24 object.__setattr__(self, "re", Q(self.re))
25 object.__setattr__(self, "im", Q(self.im))
26
27 def __add__(self, other):
28 other = other if isinstance(other, C) else C(other)
29 return C(self.re + other.re, self.im + other.im)
30
31 __radd__ = __add__
32
33 def __neg__(self):
34 return C(-self.re, -self.im)
35
36 def __sub__(self, other):
37 return self + -other
38
39 def __mul__(self, other):
40 other = other if isinstance(other, C) else C(other)
41 return C(self.re*other.re - self.im*other.im,
42 self.re*other.im + self.im*other.re)
43
44 __rmul__ = __mul__
45
46 def __truediv__(self, real):
47 return C(self.re/Q(real), self.im/Q(real))
48
49 def conj(self):
50 return C(self.re, -self.im)
51
52 def abs2(self):
53 return self.re**2 + self.im**2
54
55 def __str__(self):
56 if not self.im:
57 return str(self.re)
58 return f"{self.re}{'+' if self.im > 0 else '-'}{abs(self.im)}i"
59
60
61 def require(condition, message):
62 if not condition:
63 raise AssertionError(message)
64
65
66 def matrix(rows):
67 return tuple(tuple(x if isinstance(x, C) else C(x) for x in row)
68 for row in rows)
69
70
71 def multiply(a, b):
72 return tuple(tuple(sum((a[i][k]*b[k][j] for k in range(len(b))), C())
73 for j in range(len(b[0]))) for i in range(len(a)))
74
75
76 def scale(t, a):
77 return tuple(tuple(t*x for x in row) for row in a)
78
79
80 def add(a, b):
81 return tuple(tuple(x+y for x, y in zip(ar, br)) for ar, br in zip(a, b))
82
83
84 def adjoint(a):
85 return tuple(tuple(a[j][i].conj() for j in range(len(a)))
86 for i in range(len(a[0])))
87
88
89 def trace(a):
90 return sum((a[i][i] for i in range(len(a))), C())
91
92
93 def purity(a):
94 result = trace(multiply(a, a))
95 require(result.im == 0, "Purity must be real")
96 return result.re
97
98
99 def outer(v):
100 return tuple(tuple(x*y.conj() for y in v) for x in v)
101
102
103 def state(a=Q(3, 5), b=Q(4, 5), phase=C(1)):
104 a, b = Q(a), Q(b)
105 if a*a+b*b != 1 or phase.abs2() != 1:
106 raise ValueError("Use normalized real amplitudes and a unit phase")
107 return outer((C(a), b*phase))
108
109
110 def dephase(rho, groups=None):
111 """Keep exactly the blocks belonging to a declared address partition."""
112 n = len(rho)
113 if groups is None:
114 groups = tuple((i,) for i in range(n))
115 addresses = [i for group in groups for i in group]
116 if sorted(addresses) != list(range(n)):
117 raise ValueError("Packet groups must partition all matrix addresses")
118 packet = {i: k for k, group in enumerate(groups) for i in group}
119 return tuple(tuple(rho[i][j] if packet[i] == packet[j] else C()
120 for j in range(n)) for i in range(n))
121
122
123 def partial_dephase(rho, t):
124 """D_t = t id + (1-t) D, with t the retained coherence fraction."""
125 t = Q(t)
126 if not 0 <= t <= 1:
127 raise ValueError("Retained coherence t must lie in [0, 1]")
128 return add(scale(t, rho), scale(1-t, dephase(rho)))
129
130
131 H0 = matrix(((1, 1), (1, -1)))
132 X_PLUS = scale(Q(1, 2), matrix(((1, 1), (1, 1))))
133 Y_PLUS = scale(Q(1, 2), matrix(((1, C(0, -1)), (C(0, 1), 1))))
134 Z_PLUS = matrix(((1, 0), (0, 0)))
135 EFFECTS = (X_PLUS, Y_PLUS, Z_PLUS)
136
137
138 def hadamard(rho):
139 # H = H0/sqrt(2); the density-matrix conjugation has only rational factors.
140 return scale(Q(1, 2), multiply(multiply(H0, rho), H0))
141
142
143 def probabilities(rho):
144 """Plus probabilities in the X, Y, Z settings, evaluated by trace."""
145 values = tuple(trace(multiply(effect, rho)) for effect in EFFECTS)
146 require(all(value.im == 0 for value in values), "Probabilities must be real")
147 return tuple(value.re for value in values)
148
149
150 def reconstruct(px, py, pz):
151 """Linear tomography from ideal ensemble probabilities, not one shot."""
152 c = C(Q(px)-Q(1, 2), Q(1, 2)-Q(py))
153 return matrix(((Q(pz), c), (c.conj(), 1-Q(pz))))
154
155
156 def coherent_record(rho):
157 """|j> -> |j>|j>, in system-first order 00, 01, 10, 11."""
158 return tuple(tuple(rho[i//3][j//3] if i in (0, 3) and j in (0, 3)
159 else C() for j in range(4)) for i in range(4))
160
161
162 def uncompute_record(joint):
163 # CNOT is its own inverse and swaps the addresses 10 and 11.
164 order = (0, 1, 3, 2)
165 return tuple(tuple(joint[order[i]][order[j]] for j in range(4))
166 for i in range(4))
167
168
169 def forget_record(joint):
170 return tuple(tuple(sum((joint[2*i+k][2*j+k] for k in range(2)), C())
171 for j in range(2)) for i in range(2))
172
173
174 def circle(u):
175 u = Q(u)
176 return C((1-u*u)/(1+u*u), 2*u/(1+u*u))
177
178
179 def amplitude_damp(rho, survival=Q(4, 5), jump=Q(3, 5)):
180 """Kraus channel with gamma=jump**2 and sqrt(1-gamma)=survival."""
181 survival, jump = Q(survival), Q(jump)
182 if min(survival, jump) < 0 or survival**2+jump**2 != 1:
183 raise ValueError('Use nonnegative survival/jump amplitudes of squared norm one')
184 k0 = matrix(((1, 0), (0, survival)))
185 k1 = matrix(((0, jump), (0, 0)))
186 return add(multiply(multiply(k0, rho), adjoint(k0)),
187 multiply(multiply(k1, rho), adjoint(k1)))
188
189
190 ARITHMETIC_C = (4, 0, 2, 0, -2, 0, -4, 0, -2, 0, 2, 0)
191 ARITHMETIC_P = matrix(tuple(tuple(Q(ARITHMETIC_C[(n-m) % 12], 12)
192 for m in range(12)) for n in range(12)))
193
194
195 def arithmetic_state(phase=C(1)):
196 """Unnormalized v_+ + phase*v_-; its squared norm is 48."""
197 if phase.abs2() != 1:
198 raise ValueError('Use a unit relative phase')
199 return tuple(C(abs(x)) if x >= 0 else abs(x)*phase for x in ARITHMETIC_C)
200
201
202 def arithmetic_probability(phase):
203 v = arithmetic_state(phase)
204 projected = multiply(ARITHMETIC_P, tuple((x,) for x in v))
205 return sum(row[0].abs2() for row in projected)/48
206
207
208 SYNDROME_MASKS = (0, 4, 2, 1)
209 SYNDROME_GROUPS = tuple((mask, 7 ^ mask) for mask in SYNDROME_MASKS)
210
211
212 def encode(rho):
213 return tuple(tuple(rho[i//7][j//7] if i in (0, 7) and j in (0, 7)
214 else C() for j in range(8)) for i in range(8))
215
216
217 def bit_flip(rho, mask):
218 if mask not in range(8):
219 raise ValueError('Use a three-bit error mask')
220 return tuple(tuple(rho[i ^ mask][j ^ mask] for j in range(8)) for i in range(8))
221
222
223 def recover_logical(rho):
224 """Syndrome measurement, matching bit flip, then logical decoding."""
225 return tuple(tuple(sum((rho[group[i]][group[j]] for group in SYNDROME_GROUPS), C())
226 for j in range(2)) for i in range(2))
227
228
229 def check():
230 counts = dict(states=0, partial_channels=0, records=0, controls=0,
231 damping=0, arithmetic=0, syndrome=0)
232 phases = [circle(Q(k, 4)) for k in range(-6, 7)] + [C(-1)]
233 for k in range(-6, 7):
234 amplitudes = circle(Q(k, 4))
235 a, b = amplitudes.re, amplitudes.im
236 for phase in phases:
237 rho = state(a, b, phase)
238 require(adjoint(rho) == rho and trace(rho) == C(1), "Density matrix")
239 require(multiply(rho, rho) == rho, "Pure-state idempotence")
240 require(reconstruct(*probabilities(rho)) == rho, "Tomography")
241 require(hadamard(hadamard(rho)) == rho, "Hadamard inverse")
242 require(hadamard(rho)[0][0].re == probabilities(rho)[0], "X setting")
243 require(dephase(dephase(rho)) == dephase(rho), "Dephasing projection")
244 counts['states'] += 1
245 for t in (Q(0), Q(1, 4), Q(1, 2), Q(3, 4), Q(1)):
246 channel = partial_dephase(rho, t)
247 require(channel[0][1] == t*rho[0][1], "Coherence transfer")
248 require(purity(channel) == a**4+b**4+2*t*t*a*a*b*b, "Purity law")
249 # Independent random-unitary implementation, Z rho Z.
250 z = matrix(((1, 0), (0, -1)))
251 random_phase = add(scale((1+t)/2, rho),
252 scale((1-t)/2, multiply(multiply(z, rho), z)))
253 require(channel == random_phase, "Random-unitary channel")
254 require(reconstruct(*probabilities(channel)) == channel, "Mixed tomography")
255 counts['partial_channels'] += 1
256 joint = coherent_record(rho)
257 require(purity(joint) == 1, "Coherent record stays pure")
258 require(forget_record(joint) == dephase(rho), "Partial trace")
259 require(forget_record(uncompute_record(joint)) == rho, "Uncomputation")
260 require(forget_record(uncompute_record(dephase(joint))) == dephase(rho),
261 "Classical record does not restore phase")
262 counts['records'] += 1
263 for survival, jump in ((Q(1), Q(0)), (Q(4, 5), Q(3, 5)), (Q(0), Q(1))):
264 damped = amplitude_damp(rho, survival, jump)
265 closed = matrix(((a*a+jump*jump*b*b, survival*rho[0][1]),
266 (survival*rho[1][0], survival*survival*b*b)))
267 require(damped == closed and trace(damped) == C(1), 'Damping transfer')
268 require(purity(damped) == 1-2*jump**2*survival**2*b**4, 'Damping purity')
269 counts['damping'] += 1
270
271 require(multiply(ARITHMETIC_P, ARITHMETIC_P) == ARITHMETIC_P,
272 'Arithmetic selector is a projector')
273 require(adjoint(ARITHMETIC_P) == ARITHMETIC_P and trace(ARITHMETIC_P) == C(4),
274 'Arithmetic orthogonal projector rank')
275 require(sum(x*abs(x) for x in ARITHMETIC_C) == 0, 'Signed/unsigned orthogonality')
276 counts['arithmetic'] += 1
277 for phase in phases:
278 require(arithmetic_probability(phase) == (1-phase.re)/2, 'Arithmetic phase response')
279 rho12 = scale(Q(1, 48), outer(arithmetic_state(phase)))
280 require(trace(multiply(ARITHMETIC_P, dephase(rho12))) == C(Q(1, 3)),
281 'Arithmetic complete dephasing')
282 counts['arithmetic'] += 1
283 logical = state(phase=phase)
284 encoded = encode(logical)
285 noisy = scale(0, encoded)
286 for mask, weight in zip(SYNDROME_MASKS, (Q(1, 2), Q(1, 4), Q(1, 8), Q(1, 8))):
287 error = bit_flip(encoded, mask)
288 require(recover_logical(error) == logical, 'Single-error recovery')
289 require(dephase(error, SYNDROME_GROUPS) == error, 'Syndrome preserves logical coherence')
290 noisy = add(noisy, scale(weight, error))
291 counts['syndrome'] += 1
292 require(recover_logical(noisy) == logical, 'Mixed-error recovery')
293 require(recover_logical(dephase(noisy)) == dephase(logical), 'Fine measurement loses logical phase')
294 counts['syndrome'] += 1
295 require(recover_logical(bit_flip(encode(state()), 6)) != state(), 'Two-flip failure')
296 phase_error = state(phase=C(-1))
297 require(recover_logical(encode(phase_error)) == phase_error != state(), 'Logical phase-error failure')
298 counts['syndrome'] += 2
299
300 real, imaginary = state(), state(phase=C(0, 1))
301 require(dephase(real) == dephase(imaginary), "Equal summaries")
302 require(probabilities(real)[0] != probabilities(imaginary)[0], "Recovery obstruction")
303 require(probabilities(real) == (Q(49, 50), Q(1, 2), Q(9, 25)), "Real phase")
304 require(probabilities(imaginary) == (Q(1, 2), Q(49, 50), Q(9, 25)), "Imaginary phase")
305 oblique = state(phase=C(Q(3, 5), Q(4, 5)))
306 require(probabilities(oblique) == (Q(197, 250), Q(221, 250), Q(9, 25)), "Oblique phase")
307 require(purity(partial_dephase(real, Q(1, 2))) == Q(409, 625), "Half coherence")
308 before = dephase(hadamard(real))
309 after = hadamard(dephase(real))
310 require(before[0][0] == C(Q(49, 50)) and after[0][0] == C(Q(1, 2)), "Order witness")
311 # Attempting to double all coherences is not a channel on all input states.
312 invalid = matrix(((Q(1, 2), 1), (1, Q(1, 2))))
313 minus = (C(1), C(-1))
314 witness = sum((minus[i].conj()*invalid[i][j]*minus[j]
315 for i in range(2) for j in range(2)), C())
316 require(witness == C(-1), "Inverse fails positivity")
317 bell = matrix(((Q(1, 2), 0, 0, Q(1, 2)), (0, 0, 0, 0),
318 (0, 0, 0, 0), (Q(1, 2), 0, 0, Q(1, 2))))
319 require(dephase(bell, ((0, 3), (1, 2))) == bell, "Within-packet coherence")
320 require(purity(dephase(bell)) == Q(1, 2), "Refined measurement")
321 counts['controls'] = 10
322 return counts
323
324
325 def phase_rows():
326 phases = (('1', C(1)), ('i', C(0, 1)), ('-1', C(-1)),
327 ('-i', C(0, -1)), ('(3+4i)/5', C(Q(3, 5), Q(4, 5))))
328 return [(label, *probabilities(state(phase=z)),
329 purity(state(phase=z)), purity(dephase(state(phase=z))))
330 for label, z in phases]
331
332
333 FILES = ('measurement-settings.csv', 'partial-dephasing.csv', 'results.json',
334 'measurement-comparison.svg', 'amplitude-damping.csv',
335 'arithmetic-phase.csv', 'syndrome-recovery.csv')
336
337
338 def save(folder, counts):
339 folder = Path(folder)
340 folder.mkdir(parents=True, exist_ok=True)
341 def csv_file(name, header, rows):
342 with (folder/name).open('w', encoding='utf-8', newline='') as stream:
343 writer = csv.writer(stream, lineterminator='\n')
344 writer.writerow(header)
345 writer.writerows(rows)
346 csv_file(FILES[0], ('phase', 'p_X_plus', 'p_Y_plus', 'p_Z_plus',
347 'input_purity', 'dephased_purity'), phase_rows())
348 rows = []
349 for j in range(9):
350 t = Q(j, 8)
351 rho = partial_dephase(state(), t)
352 rows.append((t, *probabilities(rho), purity(rho), Q(288, 625)*(1-t)**2))
353 csv_file(FILES[1], ('retained_coherence', 'p_X_plus', 'p_Y_plus',
354 'p_Z_plus', 'purity', 'squared_HS_error'), rows)
355 phases = (('1', C(1)), ('i', C(0, 1)), ('-1', C(-1)),
356 ('-i', C(0, -1)), ('(3+4i)/5', C(Q(3, 5), Q(4, 5))))
357 csv_file('amplitude-damping.csv', ('phase', 'p_X_plus', 'p_Y_plus', 'p_Z_plus', 'purity'),
358 [(label, *probabilities(amplitude_damp(state(phase=z))),
359 purity(amplitude_damp(state(phase=z)))) for label, z in phases])
360 csv_file('arithmetic-phase.csv', ('phase', 'coherent_selected_probability', 'dephased_selected_probability'),
361 [(label, arithmetic_probability(z), Q(1, 3)) for label, z in phases])
362 csv_file('syndrome-recovery.csv', ('bit_flip_mask', 'first_address', 'second_address',
363 'noise_probability', 'recovered_p0', 'recovered_upper_right'),
364 [(mask, f'{group[0]:03b}', f'{group[1]:03b}', weight,
365 recover_logical(bit_flip(encode(state()), mask))[0][0],
366 recover_logical(bit_flip(encode(state()), mask))[0][1])
367 for mask, group, weight in zip(SYNDROME_MASKS, SYNDROME_GROUPS,
368 (Q(1, 2), Q(1, 4), Q(1, 8), Q(1, 8)))])
369 result = dict(scope='Exact finite teaching example; no sampling or Lean replay.',
370 checks_by_group=counts, total_cases=sum(counts.values()),
371 phase_rows=[[str(x) for x in row] for row in phase_rows()],
372 order_p0=['49/50', '1/2'],
373 tomography_oblique=['197/250', '221/250', '9/25'],
374 coarse_bell_purity='1', fine_bell_purity='1/2',
375 damping_purity=str(purity(amplitude_damp(state()))),
376 arithmetic_selected=['0', '1', '1/3'],
377 syndrome_groups=SYNDROME_GROUPS)
378 (folder/FILES[2]).write_text(json.dumps(result, indent=2)+'\n', encoding='utf-8', newline='\n')
379 # Static scientific chart; exact labels carry the data, geometry is rounded.
380 svg = ['<svg xmlns="http://www.w3.org/2000/svg" viewBox="0 0 780 510" role="img" aria-labelledby="title desc">',
381 '<title id="title">One state family, three measurement settings</title>',
382 '<desc id="desc">Five phases have the same Z-plus probability 9/25. X and Y measurements distinguish them. Exact fractions label every bar.</desc>',
383 '<rect width="780" height="510" fill="#faf8f2"/>',
384 '<g font-family="sans-serif" fill="#173b32"><text x="30" y="34" font-size="23">Relative phase changes the observable probabilities</text>',
385 '<text x="30" y="62" font-size="15">State: (3/5, 4z/5). Each labelled bar shows an exact plus probability.</text>']
386 colors = ('#245f50', '#95632a', '#536caa')
387 for index, row in enumerate(phase_rows()):
388 y = 102+index*77
389 svg.append(f'<text x="30" y="{y+18}" font-size="16">z = {row[0]}</text>')
390 for j, p in enumerate(row[1:4]):
391 x = 170+j*200
392 svg.append(f'<rect x="{x}" y="{y}" width="{float(p)*140:.3f}" height="22" fill="{colors[j]}"/>')
393 svg.append(f'<text x="{x}" y="{y+43}" font-size="14">{"XYZ"[j]}+: {p}</text>')
394 svg.append('</g></svg>')
395 (folder/FILES[3]).write_text('\n'.join(svg)+'\n', encoding='utf-8', newline='\n')
396
397
398 def main():
399 parser = argparse.ArgumentParser(description=__doc__)
400 parser.add_argument('--output', type=Path)
401 args = parser.parse_args()
402 counts = check()
403 print('Coordinate probabilities: 9/25, 16/25')
404 print('X-plus, phase 1: 49/50; phase i: 1/2')
405 print('Y-plus, phase 1: 1/2; phase i: 49/50')
406 print('Purity: original 1; dephased 337/625; half coherence 409/625')
407 print('Order witness, outcome 0: 49/50 versus 1/2')
408 print('Coherent record: reversible; classical record: phase lost')
409 print('Bell purity: parity packets 1; individual addresses 1/2')
410 print('Damping: ground population 369/625; purity 316897/390625')
411 print('Arithmetic selector: unsigned 0; signed 1; dephased 1/3')
412 print('Syndrome recovery: four single-error branches preserve the logical state')
413 print(f'PASS: {sum(counts.values())} exact cases')
414 if args.output:
415 save(args.output, counts)
416 print('Saved: '+', '.join(FILES))
417
418
419 if __name__ == '__main__':
420 main()