1"""Pauli operators and states"""
2
3from sympy import I, Mul, Add, Pow, exp, Integer
4from sympy.physics.quantum import Operator, Ket, Bra
5from sympy.physics.quantum import ComplexSpace
6from sympy.matrices import Matrix
7from sympy.functions.special.tensor_functions import KroneckerDelta
8
9__all__ = [
10    'SigmaX', 'SigmaY', 'SigmaZ', 'SigmaMinus', 'SigmaPlus', 'SigmaZKet',
11    'SigmaZBra', 'qsimplify_pauli'
12]
13
14
15class SigmaOpBase(Operator):
16    """Pauli sigma operator, base class"""
17
18    @property
19    def name(self):
20        return self.args[0]
21
22    @property
23    def use_name(self):
24        return bool(self.args[0]) is not False
25
26    @classmethod
27    def default_args(self):
28        return (False,)
29
30    def __new__(cls, *args, **hints):
31        return Operator.__new__(cls, *args, **hints)
32
33    def _eval_commutator_BosonOp(self, other, **hints):
34        return Integer(0)
35
36
37class SigmaX(SigmaOpBase):
38    """Pauli sigma x operator
39
40    Parameters
41    ==========
42
43    name : str
44        An optional string that labels the operator. Pauli operators with
45        different names commute.
46
47    Examples
48    ========
49
50    >>> from sympy.physics.quantum import represent
51    >>> from sympy.physics.quantum.pauli import SigmaX
52    >>> sx = SigmaX()
53    >>> sx
54    SigmaX()
55    >>> represent(sx)
56    Matrix([
57    [0, 1],
58    [1, 0]])
59    """
60
61    def __new__(cls, *args, **hints):
62        return SigmaOpBase.__new__(cls, *args, **hints)
63
64    def _eval_commutator_SigmaY(self, other, **hints):
65        if self.name != other.name:
66            return Integer(0)
67        else:
68            return 2 * I * SigmaZ(self.name)
69
70    def _eval_commutator_SigmaZ(self, other, **hints):
71        if self.name != other.name:
72            return Integer(0)
73        else:
74            return - 2 * I * SigmaY(self.name)
75
76    def _eval_commutator_BosonOp(self, other, **hints):
77        return Integer(0)
78
79    def _eval_anticommutator_SigmaY(self, other, **hints):
80        return Integer(0)
81
82    def _eval_anticommutator_SigmaZ(self, other, **hints):
83        return Integer(0)
84
85    def _eval_adjoint(self):
86        return self
87
88    def _print_contents_latex(self, printer, *args):
89        if self.use_name:
90            return r'{\sigma_x^{(%s)}}' % str(self.name)
91        else:
92            return r'{\sigma_x}'
93
94    def _print_contents(self, printer, *args):
95        return 'SigmaX()'
96
97    def _eval_power(self, e):
98        if e.is_Integer and e.is_positive:
99            return SigmaX(self.name).__pow__(int(e) % 2)
100
101    def _represent_default_basis(self, **options):
102        format = options.get('format', 'sympy')
103        if format == 'sympy':
104            return Matrix([[0, 1], [1, 0]])
105        else:
106            raise NotImplementedError('Representation in format ' +
107                                      format + ' not implemented.')
108
109
110class SigmaY(SigmaOpBase):
111    """Pauli sigma y operator
112
113    Parameters
114    ==========
115
116    name : str
117        An optional string that labels the operator. Pauli operators with
118        different names commute.
119
120    Examples
121    ========
122
123    >>> from sympy.physics.quantum import represent
124    >>> from sympy.physics.quantum.pauli import SigmaY
125    >>> sy = SigmaY()
126    >>> sy
127    SigmaY()
128    >>> represent(sy)
129    Matrix([
130    [0, -I],
131    [I,  0]])
132    """
133
134    def __new__(cls, *args, **hints):
135        return SigmaOpBase.__new__(cls, *args)
136
137    def _eval_commutator_SigmaZ(self, other, **hints):
138        if self.name != other.name:
139            return Integer(0)
140        else:
141            return 2 * I * SigmaX(self.name)
142
143    def _eval_commutator_SigmaX(self, other, **hints):
144        if self.name != other.name:
145            return Integer(0)
146        else:
147            return - 2 * I * SigmaZ(self.name)
148
149    def _eval_anticommutator_SigmaX(self, other, **hints):
150        return Integer(0)
151
152    def _eval_anticommutator_SigmaZ(self, other, **hints):
153        return Integer(0)
154
155    def _eval_adjoint(self):
156        return self
157
158    def _print_contents_latex(self, printer, *args):
159        if self.use_name:
160            return r'{\sigma_y^{(%s)}}' % str(self.name)
161        else:
162            return r'{\sigma_y}'
163
164    def _print_contents(self, printer, *args):
165        return 'SigmaY()'
166
167    def _eval_power(self, e):
168        if e.is_Integer and e.is_positive:
169            return SigmaY(self.name).__pow__(int(e) % 2)
170
171    def _represent_default_basis(self, **options):
172        format = options.get('format', 'sympy')
173        if format == 'sympy':
174            return Matrix([[0, -I], [I, 0]])
175        else:
176            raise NotImplementedError('Representation in format ' +
177                                      format + ' not implemented.')
178
179
180class SigmaZ(SigmaOpBase):
181    """Pauli sigma z operator
182
183    Parameters
184    ==========
185
186    name : str
187        An optional string that labels the operator. Pauli operators with
188        different names commute.
189
190    Examples
191    ========
192
193    >>> from sympy.physics.quantum import represent
194    >>> from sympy.physics.quantum.pauli import SigmaZ
195    >>> sz = SigmaZ()
196    >>> sz ** 3
197    SigmaZ()
198    >>> represent(sz)
199    Matrix([
200    [1,  0],
201    [0, -1]])
202    """
203
204    def __new__(cls, *args, **hints):
205        return SigmaOpBase.__new__(cls, *args)
206
207    def _eval_commutator_SigmaX(self, other, **hints):
208        if self.name != other.name:
209            return Integer(0)
210        else:
211            return 2 * I * SigmaY(self.name)
212
213    def _eval_commutator_SigmaY(self, other, **hints):
214        if self.name != other.name:
215            return Integer(0)
216        else:
217            return - 2 * I * SigmaX(self.name)
218
219    def _eval_anticommutator_SigmaX(self, other, **hints):
220        return Integer(0)
221
222    def _eval_anticommutator_SigmaY(self, other, **hints):
223        return Integer(0)
224
225    def _eval_adjoint(self):
226        return self
227
228    def _print_contents_latex(self, printer, *args):
229        if self.use_name:
230            return r'{\sigma_z^{(%s)}}' % str(self.name)
231        else:
232            return r'{\sigma_z}'
233
234    def _print_contents(self, printer, *args):
235        return 'SigmaZ()'
236
237    def _eval_power(self, e):
238        if e.is_Integer and e.is_positive:
239            return SigmaZ(self.name).__pow__(int(e) % 2)
240
241    def _represent_default_basis(self, **options):
242        format = options.get('format', 'sympy')
243        if format == 'sympy':
244            return Matrix([[1, 0], [0, -1]])
245        else:
246            raise NotImplementedError('Representation in format ' +
247                                      format + ' not implemented.')
248
249
250class SigmaMinus(SigmaOpBase):
251    """Pauli sigma minus operator
252
253    Parameters
254    ==========
255
256    name : str
257        An optional string that labels the operator. Pauli operators with
258        different names commute.
259
260    Examples
261    ========
262
263    >>> from sympy.physics.quantum import represent, Dagger
264    >>> from sympy.physics.quantum.pauli import SigmaMinus
265    >>> sm = SigmaMinus()
266    >>> sm
267    SigmaMinus()
268    >>> Dagger(sm)
269    SigmaPlus()
270    >>> represent(sm)
271    Matrix([
272    [0, 0],
273    [1, 0]])
274    """
275
276    def __new__(cls, *args, **hints):
277        return SigmaOpBase.__new__(cls, *args)
278
279    def _eval_commutator_SigmaX(self, other, **hints):
280        if self.name != other.name:
281            return Integer(0)
282        else:
283            return -SigmaZ(self.name)
284
285    def _eval_commutator_SigmaY(self, other, **hints):
286        if self.name != other.name:
287            return Integer(0)
288        else:
289            return I * SigmaZ(self.name)
290
291    def _eval_commutator_SigmaZ(self, other, **hints):
292        return 2 * self
293
294    def _eval_commutator_SigmaMinus(self, other, **hints):
295        return SigmaZ(self.name)
296
297    def _eval_anticommutator_SigmaZ(self, other, **hints):
298        return Integer(0)
299
300    def _eval_anticommutator_SigmaX(self, other, **hints):
301        return Integer(1)
302
303    def _eval_anticommutator_SigmaY(self, other, **hints):
304        return - I * Integer(1)
305
306    def _eval_anticommutator_SigmaPlus(self, other, **hints):
307        return Integer(1)
308
309    def _eval_adjoint(self):
310        return SigmaPlus(self.name)
311
312    def _eval_power(self, e):
313        if e.is_Integer and e.is_positive:
314            return Integer(0)
315
316    def _print_contents_latex(self, printer, *args):
317        if self.use_name:
318            return r'{\sigma_-^{(%s)}}' % str(self.name)
319        else:
320            return r'{\sigma_-}'
321
322    def _print_contents(self, printer, *args):
323        return 'SigmaMinus()'
324
325    def _represent_default_basis(self, **options):
326        format = options.get('format', 'sympy')
327        if format == 'sympy':
328            return Matrix([[0, 0], [1, 0]])
329        else:
330            raise NotImplementedError('Representation in format ' +
331                                      format + ' not implemented.')
332
333
334class SigmaPlus(SigmaOpBase):
335    """Pauli sigma plus operator
336
337    Parameters
338    ==========
339
340    name : str
341        An optional string that labels the operator. Pauli operators with
342        different names commute.
343
344    Examples
345    ========
346
347    >>> from sympy.physics.quantum import represent, Dagger
348    >>> from sympy.physics.quantum.pauli import SigmaPlus
349    >>> sp = SigmaPlus()
350    >>> sp
351    SigmaPlus()
352    >>> Dagger(sp)
353    SigmaMinus()
354    >>> represent(sp)
355    Matrix([
356    [0, 1],
357    [0, 0]])
358    """
359
360    def __new__(cls, *args, **hints):
361        return SigmaOpBase.__new__(cls, *args)
362
363    def _eval_commutator_SigmaX(self, other, **hints):
364        if self.name != other.name:
365            return Integer(0)
366        else:
367            return SigmaZ(self.name)
368
369    def _eval_commutator_SigmaY(self, other, **hints):
370        if self.name != other.name:
371            return Integer(0)
372        else:
373            return I * SigmaZ(self.name)
374
375    def _eval_commutator_SigmaZ(self, other, **hints):
376        if self.name != other.name:
377            return Integer(0)
378        else:
379            return -2 * self
380
381    def _eval_commutator_SigmaMinus(self, other, **hints):
382        return SigmaZ(self.name)
383
384    def _eval_anticommutator_SigmaZ(self, other, **hints):
385        return Integer(0)
386
387    def _eval_anticommutator_SigmaX(self, other, **hints):
388        return Integer(1)
389
390    def _eval_anticommutator_SigmaY(self, other, **hints):
391        return I * Integer(1)
392
393    def _eval_anticommutator_SigmaMinus(self, other, **hints):
394        return Integer(1)
395
396    def _eval_adjoint(self):
397        return SigmaMinus(self.name)
398
399    def _eval_mul(self, other):
400        return self * other
401
402    def _eval_power(self, e):
403        if e.is_Integer and e.is_positive:
404            return Integer(0)
405
406    def _print_contents_latex(self, printer, *args):
407        if self.use_name:
408            return r'{\sigma_+^{(%s)}}' % str(self.name)
409        else:
410            return r'{\sigma_+}'
411
412    def _print_contents(self, printer, *args):
413        return 'SigmaPlus()'
414
415    def _represent_default_basis(self, **options):
416        format = options.get('format', 'sympy')
417        if format == 'sympy':
418            return Matrix([[0, 1], [0, 0]])
419        else:
420            raise NotImplementedError('Representation in format ' +
421                                      format + ' not implemented.')
422
423
424class SigmaZKet(Ket):
425    """Ket for a two-level system quantum system.
426
427    Parameters
428    ==========
429
430    n : Number
431        The state number (0 or 1).
432
433    """
434
435    def __new__(cls, n):
436        if n not in [0, 1]:
437            raise ValueError("n must be 0 or 1")
438        return Ket.__new__(cls, n)
439
440    @property
441    def n(self):
442        return self.label[0]
443
444    @classmethod
445    def dual_class(self):
446        return SigmaZBra
447
448    @classmethod
449    def _eval_hilbert_space(cls, label):
450        return ComplexSpace(2)
451
452    def _eval_innerproduct_SigmaZBra(self, bra, **hints):
453        return KroneckerDelta(self.n, bra.n)
454
455    def _apply_operator_SigmaZ(self, op, **options):
456        if self.n == 0:
457            return self
458        else:
459            return Integer(-1) * self
460
461    def _apply_operator_SigmaX(self, op, **options):
462        return SigmaZKet(1) if self.n == 0 else SigmaZKet(0)
463
464    def _apply_operator_SigmaY(self, op, **options):
465        return I * SigmaZKet(1) if self.n == 0 else (-I) * SigmaZKet(0)
466
467    def _apply_operator_SigmaMinus(self, op, **options):
468        if self.n == 0:
469            return SigmaZKet(1)
470        else:
471            return Integer(0)
472
473    def _apply_operator_SigmaPlus(self, op, **options):
474        if self.n == 0:
475            return Integer(0)
476        else:
477            return SigmaZKet(0)
478
479    def _represent_default_basis(self, **options):
480        format = options.get('format', 'sympy')
481        if format == 'sympy':
482            return Matrix([[1], [0]]) if self.n == 0 else Matrix([[0], [1]])
483        else:
484            raise NotImplementedError('Representation in format ' +
485                                      format + ' not implemented.')
486
487
488class SigmaZBra(Bra):
489    """Bra for a two-level quantum system.
490
491    Parameters
492    ==========
493
494    n : Number
495        The state number (0 or 1).
496
497    """
498
499    def __new__(cls, n):
500        if n not in [0, 1]:
501            raise ValueError("n must be 0 or 1")
502        return Bra.__new__(cls, n)
503
504    @property
505    def n(self):
506        return self.label[0]
507
508    @classmethod
509    def dual_class(self):
510        return SigmaZKet
511
512
513def _qsimplify_pauli_product(a, b):
514    """
515    Internal helper function for simplifying products of Pauli operators.
516    """
517    if not (isinstance(a, SigmaOpBase) and isinstance(b, SigmaOpBase)):
518        return Mul(a, b)
519
520    if a.name != b.name:
521        # Pauli matrices with different labels commute; sort by name
522        if a.name < b.name:
523            return Mul(a, b)
524        else:
525            return Mul(b, a)
526
527    elif isinstance(a, SigmaX):
528
529        if isinstance(b, SigmaX):
530            return Integer(1)
531
532        if isinstance(b, SigmaY):
533            return I * SigmaZ(a.name)
534
535        if isinstance(b, SigmaZ):
536            return - I * SigmaY(a.name)
537
538        if isinstance(b, SigmaMinus):
539            return (Integer(1)/2 + SigmaZ(a.name)/2)
540
541        if isinstance(b, SigmaPlus):
542            return (Integer(1)/2 - SigmaZ(a.name)/2)
543
544    elif isinstance(a, SigmaY):
545
546        if isinstance(b, SigmaX):
547            return - I * SigmaZ(a.name)
548
549        if isinstance(b, SigmaY):
550            return Integer(1)
551
552        if isinstance(b, SigmaZ):
553            return I * SigmaX(a.name)
554
555        if isinstance(b, SigmaMinus):
556            return -I * (Integer(1) + SigmaZ(a.name))/2
557
558        if isinstance(b, SigmaPlus):
559            return I * (Integer(1) - SigmaZ(a.name))/2
560
561    elif isinstance(a, SigmaZ):
562
563        if isinstance(b, SigmaX):
564            return I * SigmaY(a.name)
565
566        if isinstance(b, SigmaY):
567            return - I * SigmaX(a.name)
568
569        if isinstance(b, SigmaZ):
570            return Integer(1)
571
572        if isinstance(b, SigmaMinus):
573            return - SigmaMinus(a.name)
574
575        if isinstance(b, SigmaPlus):
576            return SigmaPlus(a.name)
577
578    elif isinstance(a, SigmaMinus):
579
580        if isinstance(b, SigmaX):
581            return (Integer(1) - SigmaZ(a.name))/2
582
583        if isinstance(b, SigmaY):
584            return - I * (Integer(1) - SigmaZ(a.name))/2
585
586        if isinstance(b, SigmaZ):
587            # (SigmaX(a.name) - I * SigmaY(a.name))/2
588            return SigmaMinus(b.name)
589
590        if isinstance(b, SigmaMinus):
591            return Integer(0)
592
593        if isinstance(b, SigmaPlus):
594            return Integer(1)/2 - SigmaZ(a.name)/2
595
596    elif isinstance(a, SigmaPlus):
597
598        if isinstance(b, SigmaX):
599            return (Integer(1) + SigmaZ(a.name))/2
600
601        if isinstance(b, SigmaY):
602            return I * (Integer(1) + SigmaZ(a.name))/2
603
604        if isinstance(b, SigmaZ):
605            #-(SigmaX(a.name) + I * SigmaY(a.name))/2
606            return -SigmaPlus(a.name)
607
608        if isinstance(b, SigmaMinus):
609            return (Integer(1) + SigmaZ(a.name))/2
610
611        if isinstance(b, SigmaPlus):
612            return Integer(0)
613
614    else:
615        return a * b
616
617
618def qsimplify_pauli(e):
619    """
620    Simplify an expression that includes products of pauli operators.
621
622    Parameters
623    ==========
624
625    e : expression
626        An expression that contains products of Pauli operators that is
627        to be simplified.
628
629    Examples
630    ========
631
632    >>> from sympy.physics.quantum.pauli import SigmaX, SigmaY
633    >>> from sympy.physics.quantum.pauli import qsimplify_pauli
634    >>> sx, sy = SigmaX(), SigmaY()
635    >>> sx * sy
636    SigmaX()*SigmaY()
637    >>> qsimplify_pauli(sx * sy)
638    I*SigmaZ()
639    """
640    if isinstance(e, Operator):
641        return e
642
643    if isinstance(e, (Add, Pow, exp)):
644        t = type(e)
645        return t(*(qsimplify_pauli(arg) for arg in e.args))
646
647    if isinstance(e, Mul):
648
649        c, nc = e.args_cnc()
650
651        nc_s = []
652        while nc:
653            curr = nc.pop(0)
654
655            while (len(nc) and
656                   isinstance(curr, SigmaOpBase) and
657                   isinstance(nc[0], SigmaOpBase) and
658                   curr.name == nc[0].name):
659
660                x = nc.pop(0)
661                y = _qsimplify_pauli_product(curr, x)
662                c1, nc1 = y.args_cnc()
663                curr = Mul(*nc1)
664                c = c + c1
665
666            nc_s.append(curr)
667
668        return Mul(*c) * Mul(*nc_s)
669
670    return e
671