>From 46c5054a5e4886ce3c2c181520f6e9faae43a75a Mon Sep 17 00:00:00 2001
From: Matthias Dellweg <dellweg@tp1.uni-duesseldorf.de>
Date: Thu, 12 Dec 2013 11:02:30 +0100
Subject: [PATCH] Implementation of Pauli matrices.

Signed-off-by: Matthias Dellweg <dellweg@tp1.uni-duesseldorf.de>
---
 ginac/CMakeLists.txt |   2 +
 ginac/Makefile.am    |   4 +-
 ginac/ginac.h        |   1 +
 ginac/pauli.cpp      | 369 +++++++++++++++++++++++++++++++++++++++++++++++++++
 ginac/pauli.h        | 139 +++++++++++++++++++
 5 files changed, 513 insertions(+), 2 deletions(-)
 create mode 100644 ginac/pauli.cpp
 create mode 100644 ginac/pauli.h

diff --git a/ginac/CMakeLists.txt b/ginac/CMakeLists.txt
index ed6ac11..0fa857e 100644
--- a/ginac/CMakeLists.txt
+++ b/ginac/CMakeLists.txt
@@ -36,6 +36,7 @@ set(ginaclib_sources
     parser/parse_context.cpp
     parser/parser_compat.cpp
     parser/parser.cpp
+    pauli.cpp
     polynomial/chinrem_gcd.cpp
     polynomial/collect_vargs.cpp
     polynomial/cra_garner.cpp
@@ -107,6 +108,7 @@ set(ginaclib_public_headers
     wildcard.h 
     parser/parser.h 
     parser/parse_context.h
+    pauli.h
 )
 
 set(ginaclib_private_headers
diff --git a/ginac/Makefile.am b/ginac/Makefile.am
index e149d0c..e130058 100644
--- a/ginac/Makefile.am
+++ b/ginac/Makefile.am
@@ -7,7 +7,7 @@ libginac_la_SOURCES = add.cpp archive.cpp basic.cpp clifford.cpp color.cpp \
   inifcns_trans.cpp inifcns_gamma.cpp inifcns_nstdsums.cpp \
   integral.cpp lst.cpp matrix.cpp mul.cpp ncmul.cpp normal.cpp numeric.cpp \
   operators.cpp power.cpp registrar.cpp relational.cpp remember.cpp \
-  pseries.cpp print.cpp symbol.cpp symmetry.cpp tensor.cpp \
+  pseries.cpp print.cpp symbol.cpp symmetry.cpp tensor.cpp pauli.cpp \
   utils.cpp wildcard.cpp \
   remember.h tostring.h utils.h crc32.h hash_seed.h compiler.h \
   parser/parse_binop_rhs.cpp \
@@ -63,7 +63,7 @@ ginacinclude_HEADERS = ginac.h add.h archive.h assertion.h basic.h class_info.h
   exprseq.h fail.h factor.h fderivative.h flags.h function.h hash_map.h idx.h indexed.h \
   inifcns.h integral.h lst.h matrix.h mul.h ncmul.h normal.h numeric.h operators.h \
   power.h print.h pseries.h ptr.h registrar.h relational.h structure.h \
-  symbol.h symmetry.h tensor.h version.h wildcard.h \
+  symbol.h symmetry.h tensor.h version.h wildcard.h pauli.h \
   parser/parser.h \
   parser/parse_context.h
 
diff --git a/ginac/ginac.h b/ginac/ginac.h
index 27d935c..31539d1 100644
--- a/ginac/ginac.h
+++ b/ginac/ginac.h
@@ -63,6 +63,7 @@
 #include "idx.h"
 #include "indexed.h"
 #include "tensor.h"
+#include "pauli.h"
 #include "color.h"
 #include "clifford.h"
 
diff --git a/ginac/pauli.cpp b/ginac/pauli.cpp
new file mode 100644
index 0000000..1d4736e
--- /dev/null
+++ b/ginac/pauli.cpp
@@ -0,0 +1,369 @@
+/** @file pauli.cpp
+ *
+ *  Implementation of GiNaC's pauli matrix objects. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#include "pauli.h"
+#include "idx.h"
+#include "ncmul.h"
+#include "symmetry.h"
+#include "operators.h"
+#include "numeric.h"
+#include "mul.h"
+#include "power.h" // for sqrt()
+#include "symbol.h"
+#include "archive.h"
+#include "utils.h"
+
+#include <iostream>
+#include <stdexcept>
+
+namespace GiNaC {
+
+GINAC_IMPLEMENT_REGISTERED_CLASS(pauli, indexed)
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(paulione, tensor,
+  print_func<print_dflt>(&paulione::do_print).
+  print_func<print_latex>(&paulione::do_print_latex))
+
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(paulisigma, tensor,
+  print_func<print_dflt>(&paulisigma::do_print).
+  print_func<print_latex>(&paulisigma::do_print))
+
+//////////
+// default constructors
+//////////
+
+pauli::pauli() : representation_label(0)
+{
+}
+
+DEFAULT_CTOR(paulione)
+DEFAULT_CTOR(paulisigma)
+
+//////////
+// other constructors
+//////////
+
+/** Construct object without any pauli index. This constructor is for
+ *  internal use only. Use the pauli_ONE() function instead.
+ *  @see pauli_ONE */
+pauli::pauli(const ex & b, unsigned char rl) : inherited(b), representation_label(rl)
+{
+}
+
+/** Construct object with one pauli index. This constructor is for internal
+ *  use only. Use the pauli_sigma() function instead.
+ *  @see pauli_sigma */
+pauli::pauli(const ex & b, const ex & i1, unsigned char rl) : inherited(b, i1), representation_label(rl)
+{
+}
+
+pauli::pauli(unsigned char rl, const exvector & v, bool discardable) : inherited(not_symmetric(), v, discardable), representation_label(rl)
+{
+}
+
+pauli::pauli(unsigned char rl, std::auto_ptr<exvector> vp) : inherited(not_symmetric(), vp), representation_label(rl)
+{
+}
+
+return_type_t pauli::return_type_tinfo() const
+{
+	return make_return_type_t<pauli>(representation_label);
+}
+
+//////////
+// archiving
+//////////
+
+void pauli::read_archive(const archive_node& n, lst& sym_lst)
+{
+	inherited::read_archive(n, sym_lst);
+	unsigned rl;
+	n.find_unsigned("label", rl);
+	representation_label = rl;
+}
+
+void pauli::archive(archive_node &n) const
+{
+	inherited::archive(n);
+	n.add_unsigned("label", representation_label);
+}
+
+GINAC_BIND_UNARCHIVER(pauli);
+GINAC_BIND_UNARCHIVER(paulione);
+GINAC_BIND_UNARCHIVER(paulisigma);
+
+//////////
+// functions overriding virtual functions from base classes
+//////////
+
+int pauli::compare_same_type(const basic & other) const
+{
+	GINAC_ASSERT(is_a<pauli>(other));
+	const pauli &o = static_cast<const pauli &>(other);
+
+	if (representation_label != o.representation_label) {
+		// different representation label
+		return representation_label < o.representation_label ? -1 : 1;
+	}
+
+	return inherited::compare_same_type(other);
+}
+
+bool pauli::match_same_type(const basic & other) const
+{
+	GINAC_ASSERT(is_a<pauli>(other));
+	const pauli &o = static_cast<const pauli &>(other);
+
+	return representation_label == o.representation_label;
+}
+
+DEFAULT_COMPARE(paulione)
+DEFAULT_COMPARE(paulisigma)
+
+DEFAULT_PRINT_LATEX(paulione, "ONE", "\\mathbb{1}")
+DEFAULT_PRINT(paulisigma, "\\sigma")
+
+/** Perform automatic simplification on noncommutative product of pauli
+ *  objects. This removes superfluous ONEs. */
+ex pauli::eval_ncmul(const exvector & v) const
+{
+	exvector s;
+	s.reserve(v.size());
+
+	// Remove superfluous ONEs
+	exvector::const_iterator it = v.begin(), itend = v.end();
+	while (it != itend) {
+		if (!is_a<paulione>(it->op(0)))
+			s.push_back(*it);
+		it++;
+	}
+
+	if (s.empty())
+		return pauli(paulione(), representation_label);
+	else
+		return hold_ncmul(s);
+}
+
+ex pauli::thiscontainer(const exvector & v) const
+{
+	return pauli(representation_label, v);
+}
+
+ex pauli::thiscontainer(std::auto_ptr<exvector> vp) const
+{
+	return pauli(representation_label, vp);
+}
+
+/** Contraction of pauli matrix with something else. */
+bool paulisigma::contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const
+{
+	GINAC_ASSERT(is_a<indexed>(*self));
+	GINAC_ASSERT(is_a<indexed>(*other));
+	GINAC_ASSERT(self->nops() == 2);
+	GINAC_ASSERT(is_a<paulisigma>(self->op(0)));
+	unsigned char rl = ex_to<pauli>(*self).get_representation_label();
+
+	if (is_exactly_a<paulisigma>(other->op(0))) {
+
+		// Contraction only makes sense if the represenation labels are equal
+		GINAC_ASSERT(is_a<pauli>(*other));
+		if (ex_to<pauli>(*other).get_representation_label() != rl)
+			return false;
+
+		// sigma.a sigma.a = 3 ONE
+		if (other - self == 1) {
+			*self = numeric(3);
+			*other = pauli_ONE(rl);
+			return true;
+
+		// sigma.a sigma.b sigma.a = -1 sigma.b
+		} else if (other - self == 2
+		        && is_a<pauli>(self[1])) {
+			*self = numeric(-1);
+			*other = _ex1;
+			return true;
+
+		// sigma.a S sigma.a = 2 Tr(S) - S
+		} else {
+			exvector::iterator it = self + 1;
+			while (it != other) {
+				if (!is_a<pauli>(*it)) {
+					return false;
+				}
+				it++;
+			}
+
+			it = self + 1;
+			ex S = _ex1;
+			while (it != other) {
+				S *= *it;
+				*it++ = _ex1;
+			}
+
+			*self = pauli_trace(S, rl) * pauli_ONE(rl) * 2 - S;
+			*other = _ex1;
+			return true;
+		}
+	}
+
+	return false;
+}
+
+//////////
+// global functions
+//////////
+
+ex pauli_ONE(unsigned char rl)
+{
+	static ex ONE = (new paulione)->setflag(status_flags::dynallocated);
+	return pauli(ONE, rl);
+}
+
+ex pauli_sigma(const ex & a, unsigned char rl)
+{
+	static ex t = (new paulisigma)->setflag(status_flags::dynallocated);
+
+	if (!is_a<idx>(a))
+		throw(std::invalid_argument("indices of pauli_sigma must be of type idx"));
+	if (!ex_to<idx>(a).get_dim().is_equal(3))
+		throw(std::invalid_argument("index dimension for pauli_sigma must be 3"));
+
+	return pauli(t, a, rl);
+}
+
+/** Check whether a given tinfo key (as returned by return_type_tinfo()
+ *  is that of a pauli object (with an arbitrary representation label). */
+static bool is_pauli_tinfo(const return_type_t& ti)
+{
+	return *(ti.tinfo) == typeid(pauli);
+}
+
+/** Extract representation label from tinfo key (as returned by
+ *  return_type_tinfo()). */
+static unsigned char get_representation_label(const return_type_t& ti)
+{
+	return (unsigned char)ti.rl;
+}
+
+ex pauli_trace(const ex & e, const std::set<unsigned char> & rls)
+{
+	if (is_a<pauli>(e)) {
+
+		unsigned char rl = ex_to<pauli>(e).get_representation_label();
+
+		// Are we taking the trace over this object's representation label?
+		if (rls.find(rl) == rls.end())
+			return e;
+
+		// Yes, all matrices are traceless, except for pauli_ONE
+		if (is_a<paulione>(e.op(0)))
+			return _ex2;
+		else
+			return _ex0;
+
+	} else if (is_exactly_a<mul>(e)) {
+
+		// Trace of product: pull out non-pauli factors
+		ex prod = _ex1;
+		for (size_t i=0; i<e.nops(); i++) {
+			const ex &o = e.op(i);
+			if (is_pauli_tinfo(o.return_type_tinfo()))
+				prod *= pauli_trace(o, rls);
+			else
+				prod *= o;
+		}
+		return prod;
+
+	} else if (is_exactly_a<ncmul>(e)) {
+
+		unsigned char rl = get_representation_label(e.return_type_tinfo());
+
+		// Are we taking the trace over this string's representation label?
+		if (rls.find(rl) == rls.end())
+			return e;
+
+		// Yes, expand product if necessary
+		ex e_expanded = e.expand();
+		if (!is_a<ncmul>(e_expanded))
+			return pauli_trace(e_expanded, rls);
+
+		size_t num = e.nops();
+
+		if (num == 2) {
+
+			// Tr sigma_a sigma_b = 2 delta_a_b
+			return delta_tensor(e.op(0).op(1), e.op(1).op(1)) * 2;
+
+		} else {
+
+			// Traces of 3 or more matrices are computed recursively:
+			// Tr sigma_a1 .. sigma_an =
+			//     2 delta_a(n-1)_an Tr sigma_a1 .. sigma_a(n-2)
+			//   + i epsilon_a(n-1)_an_k Tr sigma_a1 .. sigma_a(n-2) sigma_k
+			const ex &last_index = e.op(num - 1).op(1);
+			const ex &next_to_last_index = e.op(num - 2).op(1);
+			idx summation_index((new symbol)->setflag(status_flags::dynallocated), 3);
+
+			exvector v1;
+			v1.reserve(num - 2);
+			for (size_t i=0; i<num-2; i++)
+				v1.push_back(e.op(i));
+
+			exvector v2 = v1;
+			v2.push_back(pauli_sigma(summation_index, rl));
+
+			return delta_tensor(next_to_last_index, last_index) * pauli_trace(ncmul(v1), rl) * 2
+			       + epsilon_tensor(next_to_last_index, last_index, summation_index) * pauli_trace(ncmul(v2), rl) * I;
+		}
+
+	} else if (e.nops() > 0) {
+
+		// Trace maps to all other container classes (this includes sums)
+		pointer_to_map_function_1arg<const std::set<unsigned char> &> fcn(pauli_trace, rls);
+		return e.map(fcn);
+
+	} else
+		return _ex0;
+}
+
+ex pauli_trace(const ex & e, const lst & rll)
+{
+	// Convert list to set
+	std::set<unsigned char> rls;
+	for (lst::const_iterator i = rll.begin(); i != rll.end(); ++i) {
+		if (i->info(info_flags::nonnegint))
+			rls.insert(ex_to<numeric>(*i).to_int());
+	}
+
+	return pauli_trace(e, rls);
+}
+
+ex pauli_trace(const ex & e, unsigned char rl)
+{
+	// Convert label to set
+	std::set<unsigned char> rls;
+	rls.insert(rl);
+
+	return pauli_trace(e, rls);
+}
+
+} // namespace GiNaC
diff --git a/ginac/pauli.h b/ginac/pauli.h
new file mode 100644
index 0000000..f98589a
--- /dev/null
+++ b/ginac/pauli.h
@@ -0,0 +1,139 @@
+/** @file pauli.h
+ *
+ *  Interface to GiNaC's pauli matrix objects. */
+
+/*
+ *  GiNaC Copyright (C) 1999-2011 Johannes Gutenberg University Mainz, Germany
+ *
+ *  This program is free software; you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation; either version 2 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program; if not, write to the Free Software
+ *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
+ */
+
+#ifndef GINAC_PAULI_H
+#define GINAC_PAULI_H
+
+#include "indexed.h"
+#include "tensor.h"
+
+#include <set>
+
+namespace GiNaC {
+
+/** This class holds a sigma_matix or the unity element.
+ *  A representation label (an unsigned 8-bit integer) is used to distinguish
+ *  elements from different algebras (objects with different labels
+ *  commutate). These objects implement an abstract representation,
+ *  not a specific matrix representation. The indices used for pauli
+ *  objects should not have a variance. */
+class pauli : public indexed
+{
+	GINAC_DECLARE_REGISTERED_CLASS(pauli, indexed)
+	// other constructors
+public:
+	pauli(const ex & b, unsigned char rl = 0);
+	pauli(const ex & b, const ex & i1, unsigned char rl = 0);
+
+	// internal constructors
+	pauli(unsigned char rl, const exvector & v, bool discardable = false);
+	pauli(unsigned char rl, std::auto_ptr<exvector> vp);
+	void archive(archive_node& n) const;
+	void read_archive(const archive_node& n, lst& sym_lst);
+
+	// functions overriding virtual functions from base classes
+protected:
+	ex eval_ncmul(const exvector & v) const;
+	bool match_same_type(const basic & other) const;
+	ex thiscontainer(const exvector & v) const;
+	ex thiscontainer(std::auto_ptr<exvector> vp) const;
+	unsigned return_type() const { return return_types::noncommutative; }
+	return_type_t return_type_tinfo() const;
+
+	// non-virtual functions in this class
+public:
+	unsigned char get_representation_label() const {return representation_label;}
+
+	// member variables
+private:
+	unsigned char representation_label; /**< Representation label to distinguish independent pauli matrices coming from separated fermion lines */
+};
+GINAC_DECLARE_UNARCHIVER(pauli);
+
+
+/** This class represents the pauli unity element. */
+class paulione : public tensor
+{
+	GINAC_DECLARE_REGISTERED_CLASS(paulione, tensor)
+
+	// non-virtual functions in this class
+protected:
+	void do_print(const print_context & c, unsigned level) const;
+	void do_print_latex(const print_latex & c, unsigned level) const;
+};
+GINAC_DECLARE_UNARCHIVER(paulione);
+
+/** This class represents a pauli matrix. */
+class paulisigma : public tensor
+{
+	GINAC_DECLARE_REGISTERED_CLASS(paulisigma, tensor)
+
+	// functions overriding virtual functions from base classes
+public:
+	bool contract_with(exvector::iterator self, exvector::iterator other, exvector & v) const;
+
+	// non-virtual functions in this class
+protected:
+	void do_print(const print_context & c, unsigned level) const;
+	void do_print_latex(const print_latex & c, unsigned level) const;
+};
+GINAC_DECLARE_UNARCHIVER(paulisigma);
+
+
+// global functions
+
+/** Create the pauli unity element. This is an indexed object, although it
+ *  has no indices.
+ *
+ *  @param rl Representation label
+ *  @return newly constructed unity element */
+ex pauli_ONE(unsigned char rl = 0);
+
+/** Create a pauli matrix.
+ *
+ *  @param a Index
+ *  @param rl Representation label
+ *  @return newly constructed unity generator */
+ex pauli_sigma(const ex & a, unsigned char rl = 0);
+
+/** Calculate pauli traces over the specified set of representation labels.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rls Set of representation labels */
+ex pauli_trace(const ex & e, const std::set<unsigned char> & rls);
+
+/** Calculate pauli traces over the specified list of representation labels.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rll List of representation labels */
+ex pauli_trace(const ex & e, const lst & rll);
+
+/** Calculate the trace of an expression containing pauli objects with a
+ *  specified representation label.
+ *
+ *  @param e Expression to take the trace of
+ *  @param rl Representation label */
+ex pauli_trace(const ex & e, unsigned char rl = 0);
+
+} // namespace GiNaC
+
+#endif // ndef GINAC_PAULI_H
-- 
1.8.5.2

