mv FC to HEP

HEP/DGamma.cpp

+/**
+ * @file
+ * @brief Functions for DGamma
+ */
+ 
+#include "HEP.h"
+
+namespace HepLib {
+
+    static ex expl_TR_diff2(const ex & arg, const symbol & s) {
+        auto ret = arg.diff(s);
+        if(!is_a<add>(ret)) ret = lst{ret};
+        ex res = 0;
+        for(auto item : ret) {
+            if(!is_a<mul>(item)) res += TR(item);
+            else {
+                ex c=1; ex v=1;
+                for(auto it : item) {
+                    if(!Index::has(it) && !DGamma::has(it)) c *= it;
+                    else v *= it;
+                }
+                res += c * TR(v);
+            }
+        }
+        return res;
+    }
+    
+    static ex expl_TR_diff(const ex & arg, const symbol & s) {
+        auto ret = arg.diff(s);
+        auto nd = numer_denom(ret.normal());
+        auto num = collect_common_factors(nd.op(0));
+        if(!is_a<mul>(num)) return TR(num)/nd.op(1);
+        ex c=1;
+        ex v=1;
+        for(auto it : num) {
+            if(!Index::has(it) && !DGamma::has(it)) c *= it;
+            else v *= it;
+        }
+        return c*TR(v)/nd.op(1);
+    }
+    
+    //-----------------------------------------------------------
+    // DGamma Class
+    //-----------------------------------------------------------
+    GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(DGamma, basic,
+        print_func<print_dflt>(&DGamma::print).
+        print_func<FormFormat>(&DGamma::form_print).
+        print_func<FCFormat>(&DGamma::fc_print)
+    )
+    
+    DEFAULT_CTOR(DGamma)
+    GINAC_BIND_UNARCHIVER(DGamma);
+    IMPLEMENT_HAS(DGamma)
+    IMPLEMENT_ALL(DGamma)
+        
+    return_type_t DGamma::return_type_tinfo() const {
+        return make_return_type_t<clifford>(rl);
+        //return make_return_type_t<DGamma>(rl);
+    }
+    
+    bool DGamma::match_same_type(const basic & other) const {
+        const DGamma &o = static_cast<const DGamma &>(other);
+        return rl == o.rl;
+    }
+    
+    unsigned DGamma::get_rl() {
+        return rl;
+    }
+    
+    int DGamma::compare_same_type(const basic &other) const {
+        const DGamma &o = static_cast<const DGamma &>(other);
+        if (rl != o.rl) return rl < o.rl ? -1 : 1;
+        return pi.compare(o.pi);
+    }
+    
+    DGamma::DGamma(int int_1567, unsigned _rl) : pi(int_1567), rl(_rl) { }
+    DGamma::DGamma(const Vector &p, unsigned _rl) : pi(p), rl(_rl) { }
+    DGamma::DGamma(const Index &i, unsigned _rl) : pi(i), rl(_rl) { }
+    DGamma::DGamma(const DGamma &g, unsigned _rl) : pi(g.pi), rl(_rl) { }
+    
+    size_t DGamma::nops() const { return 2; }
+    ex DGamma::op(size_t i) const {
+        if(i==0) return pi;
+        else if(i==1) return rl;
+        return 0;
+    }
+    
+    ex & DGamma::let_op(size_t i) {
+        static ex ex_rl = numeric(rl);
+        ensure_if_modifiable();
+        if(i==0) return pi;
+        else return ex_rl;
+    }
+    
+    ex DGamma::eval() const {
+        if(flags & status_flags::evaluated) return *this;
+        else if(is_zero(pi) || is_zero(pi-1) || is_zero(pi-5) || is_zero(pi-6) || is_zero(pi-7)) return this->hold();
+        else if(!is_a<Vector>(pi) && !is_a<Index>(pi)) return GAS(pi,rl);
+        else return this->hold();
+    }
+    
+    void DGamma::print(const print_dflt &c, unsigned level) const {
+        if(is_zero(pi-1)) {
+            c.s << "𝕚";
+            return;
+        }
+        c.s << "(" << "𝛾";
+        if(is_a<numeric>(pi)) c.s << pi;
+        else c.s << "." << pi;
+        c.s << ")";
+    }
+    
+    void DGamma::form_print(const FormFormat &c, unsigned level) const {
+        c << "g_(" << rl;
+        if(!is_zero(pi-1)) {
+            c << "," << pi;
+            if(is_a<numeric>(pi)) c << "_"; // 5_, 6_, 7_ in form
+        }
+        c << ")";
+    }
+    
+    void DGamma::fc_print(const FCFormat &c, unsigned level) const {
+        if(is_zero(pi-1)) {
+            c << "1";
+            return;
+        }
+        if(is_a<Vector>(pi)) c << "GSD";
+        else c << "GAD";
+        c << "(" << pi << ")";
+    }
+    
+    void DGamma::archive(archive_node & n) const {
+        inherited::archive(n);
+        n.add_ex("pi", pi);
+        n.add_unsigned("rl", rl);
+    }
+    
+    void DGamma::read_archive(const archive_node& n, lst& sym_lst) {
+        inherited::read_archive(n, sym_lst);
+        n.find_unsigned("rl", rl);
+        n.find_ex("pi", pi, sym_lst);
+    }
+    
+    ex DGamma::derivative(const symbol & s) const {
+        return 0;
+    }
+    
+    ex DGamma::conjugate() const {
+        if(is_a<Index>(pi) || is_a<Vector>(pi)) return *this;
+        else if(is_zero(pi-5)) return -1*DGamma(5, rl);
+        else if(is_zero(pi-6)) return DGamma(7, rl);
+        else if(is_zero(pi-7)) return DGamma(6, rl);
+        else if(is_zero(pi-1)) return DGamma(1, rl);
+        throw Error("invalid Dirac Gamma Found.");
+    }
+    
+    //-----------------------------------------------------------
+    // TR/GAS functions
+    //-----------------------------------------------------------
+    namespace {
+        void TR_form_print(const ex &arg, const print_context &c0) {
+            auto c = static_cast<const FormFormat &>(c0);
+            c << "(" << arg << ")";
+        }
+        void TR_fc_print(const ex &arg, const print_context &c0) {
+            auto c = static_cast<const FCFormat &>(c0);
+            c << "DiracTrace(" << arg << ")";
+        }
+        ex tr_conj(const ex & e) {
+            return TR(e.conjugate());
+        }
+        
+        void TTR_form_print(const ex &arg, const print_context &c0) {
+            auto c = static_cast<const FormFormat &>(c0);
+            if(!is_a<lst>(arg)) c << "TTR(" << arg << ")";
+            else {
+                bool first = true;
+                for(auto item : arg) {
+                    if(first) { first=false; c << "TTR(" << item; }
+                    else c << "," << item;
+                }
+                c << ")";
+            }
+        }
+        void TTR_fc_print(const ex &arg, const print_context &c0) {
+            auto c = static_cast<const FCFormat &>(c0);
+            if(!is_a<lst>(arg)) c << "SUNTrace(SUNT(" << arg << "))";
+            else {
+                bool first = true;
+                for(auto item : arg) {
+                    if(first) { first=false; c << "SUNTrace(SUNT(" << item; }
+                    else c << "," << item;
+                }
+                c << "))";
+            }
+        }
+        ex ttr_conj(const ex & e) {
+            lst argv;
+            if(!is_a<lst>(e)) return TTR(e);
+            else argv = ex_to<lst>(e);
+            lst as;
+            for(auto it=argv.rbegin(); it!=argv.rend(); ++it) as.append(*it);
+            return TTR(as);
+        }
+    }
+    
+    REGISTER_FUNCTION(TR, do_not_evalf_params().
+        conjugate_func(tr_conj).
+        print_func<FormFormat>(&TR_form_print).
+        print_func<FCFormat>(&TR_fc_print).
+        set_return_type(return_types::commutative).
+        expl_derivative_func(expl_TR_diff)
+    );
+    
+    REGISTER_FUNCTION(TTR, do_not_evalf_params().
+        conjugate_func(ttr_conj).
+        print_func<FormFormat>(&TTR_form_print).
+        print_func<FCFormat>(&TTR_fc_print)
+    );
+    
+    REGISTER_FUNCTION(HF, do_not_evalf_params());
+    
+    /**
+     * @brief function similar to GAD/GSD in FeynClac
+     * @param expr momentum/index or 1,5,6,7
+     * @param rl the represent line number
+     * @return expanded/translasted to Dirac Gamma objects
+     */
+    ex GAS(const ex &expr, unsigned rl) {
+        if(is_zero(expr-1)) return DGamma(1,rl);
+        else if(is_zero(expr-5)) return DGamma(5,rl);
+        else if(is_zero(expr-6)) return DGamma(6,rl);
+        else if(is_zero(expr-7)) return DGamma(7,rl);
+        
+        ex tmp = expand(expr);
+        lst ex_lst;
+        if(is_a<add>(tmp)) {
+            for(auto item : tmp) ex_lst.append(item);
+        } else ex_lst.append(tmp);
+        ex res = 0;
+        for(auto item : ex_lst) {
+            lst mul_lst;
+            if(is_a<mul>(item)) {
+                for(auto ii : item) mul_lst.append(ii);
+            } else mul_lst.append(item);
+            ex c=1, g=1;
+            for(auto ii : mul_lst) {
+                if(is_a<Vector>(ii)) {
+                    if(is_a<DGamma>(g)) throw Error("Something Wrong with GAS");
+                    g = DGamma(ex_to<Vector>(ii),rl);
+                } else if(is_a<Index>(ii)) {
+                    if(is_a<DGamma>(g)) throw Error("Something Wrong with GAS");
+                    g = DGamma(ex_to<Index>(ii),rl);
+                } else {
+                    c *= ii;
+                }
+            }
+            if(!is_a<DGamma>(g)) throw Error("Something Wrong with GAS");
+            res += c * g;
+        }
+        return res;
+    }
+    
+    
+
+
+}
+

HEP/Eps.cpp

+/**
+ * @file
+ * @brief Functions for Levi-Civita Tensor
+ */
+ 
+#include "HEP.h"
+
+namespace HepLib {
+
+    //-----------------------------------------------------------
+    // Eps Class
+    //-----------------------------------------------------------
+    GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Eps, basic,
+        print_func<print_dflt>(&Eps::print).
+        print_func<FormFormat>(&Eps::form_print).
+        print_func<FCFormat>(&Eps::fc_print)
+    )
+    
+    DEFAULT_CTOR(Eps)
+    GINAC_BIND_UNARCHIVER(Eps);
+    IMPLEMENT_HAS(Eps)
+    IMPLEMENT_ALL(Eps)
+
+    Eps::Eps(const Vector &x1, const Vector &x2, const Vector &x3, const Vector &x4) : pis{x1,x2,x3,x4} { }
+    Eps::Eps(const Vector &x1, const Vector &x2, const Vector &x3, const Index &x4) : pis{x1,x2,x3,x4} { }
+    Eps::Eps(const Vector &x1, const Vector &x2, const Index &x3, const Index &x4) : pis{x1,x2,x3,x4} { }
+    Eps::Eps(const Vector &x1, const Index &x2, const Index &x3, const Index &x4) : pis{x1,x2,x3,x4} { }
+    Eps::Eps(const Index &x1, const Index &x2, const Index &x3, const Index &x4) : pis{x1,x2,x3,x4} { }
+    Eps::Eps(vector<Vector> vs, vector<Index> is) {
+        int i=0;
+        for(auto vi : vs) pis[i++] = vi;
+        for(auto ii : is) pis[i++] = ii;
+    }
+
+    int Eps::compare_same_type(const basic &other) const {
+        const Eps &o = static_cast<const Eps &>(other);
+        for(int i=0; i<4; i++) {
+            auto c = pis[i].compare(o.pis[i]);
+            if(c!=0) return c;
+        }
+        return 0;
+    }
+    
+    ex Eps::eval() const {
+        if(flags & status_flags::evaluated) return *this;
+        bool ok = true;
+        int ii = 4;
+        for(int i=0; i<4; i++) {
+            if(!is_a<Vector>(pis[i]) && !is_a<Index>(pis[i])) ok = false;
+            else if(is_a<Vector>(pis[i]) && ii!=4) ok = false; // Vector after Index
+            else if(is_a<Index>(pis[i]) && ii==4) ii = i;
+            if(!ok) break;
+        }
+        if(!ok) return LC(pis[0],pis[1],pis[2],pis[3]);
+        
+        if(isSorted(ii,pis) && isSorted(4-ii,pis+ii)) return this->hold();
+        else {
+            ex pis2[4];
+            for(int i=0; i<4; i++) pis2[i] = pis[i];
+            int ac1 = ACSort(ii,pis2);
+            int ac2 = ACSort(4-ii,pis2+ii);
+            if(ac1 * ac2==0) return 0;
+            return ac1 * ac2 * LC(pis2[0],pis2[1],pis2[2],pis2[3]);
+        }
+    }
+    
+    void Eps::print(const print_dflt &c, unsigned level) const {
+        c.s << "𝜀" << "(";
+        for(int i=0; i<3; i++) c.s << pis[i] << ",";
+        c.s << pis[3] << ")";
+    }
+    
+    void Eps::form_print(const FormFormat &c, unsigned level) const {
+        c << "e_(";
+        for(int i=0; i<3; i++) c << pis[i] << ",";
+        c << pis[3] << ")";
+    }
+    
+    void Eps::fc_print(const FCFormat &c, unsigned level) const {
+        c << "LCD[";
+        bool first = true;
+        for(int i=3; i>=0; i--) {
+            if(is_a<Vector>(pis[i]) && first) {
+                c << "][";
+                first = false;
+            }
+            c << pis[i];
+            
+            if(i==0) c << "]";
+            else if(!first || !is_a<Vector>(pis[i-1])) c << ",";
+        }
+        if(first) c << "[]";
+    }
+    
+    size_t Eps::nops() const { return 4; }
+    ex Eps::op(size_t i) const {
+        return pis[i];
+    }
+    ex & Eps::let_op(size_t i) {
+        ensure_if_modifiable();
+        return pis[i];
+    }
+    
+    void Eps::archive(archive_node & n) const {
+        inherited::archive(n);
+        for(int i=0; i<4; i++) n.add_ex("pis"+to_string(i), pis[i]);
+    }
+    
+    void Eps::read_archive(const archive_node& n, lst& sym_lst) {
+        inherited::read_archive(n, sym_lst);
+        for(int i=0; i<4; i++) {
+            n.find_ex("pis"+to_string(i), pis[i], sym_lst);
+        }
+    }
+    
+    ex Eps::derivative(const symbol & s) const {
+        return 0;
+    }
+    
+    /**
+     * @brief function similar to LCD in FeynCalc
+     * @param pi1 vector/index in the 1st position
+     * @param pi2 vector/index in the 2nd position
+     * @param pi3 vector/index in the 3rd position
+     * @param pi4 vector/index in the 4th position
+     * @return expanded/translated to Eps objects
+     */
+    ex LC(ex pi1, ex pi2, ex pi3, ex pi4) {
+        bool isEps = true;
+        lst pis = lst {pi1, pi2, pi3, pi4};
+        for(auto pi : pis) {
+            if(!is_a<Vector>(pi) && !is_a<Index>(pi)) {
+                isEps = false;
+                break;
+            }
+        }
+        
+        if(isEps) {
+            vector<Vector> vs;
+            vector<Index> is;
+            ex sign = 1;
+            for(int i=0; i<4; i++) {
+                if(is_a<Vector>(pis.op(i))) {
+                    vs.push_back(ex_to<Vector>(pis.op(i)));
+                    sign *= pow(-1, is.size());
+                } else is.push_back(ex_to<Index>(pis.op(i)));
+            }
+            return sign*Eps(vs, is);
+        }
+        
+        for(int i=0; i<4; i++) {
+            auto pi = mma_collect_lst(pis.op(i), [](const ex & e)->bool{return Index::has(e) || Vector::has(e);});
+            for(auto item : pi) { // check 
+                if(!is_a<Vector>(item.op(1)) && !is_a<Index>(item.op(1))) {
+                    cout << "pi = " << pi << endl;
+                    throw Error("LC Error: there is no Index or Vector.");
+                }
+            }
+            pis.let_op(i) = pi;
+        }
+        
+        ex res = 0;
+        for(auto i0 : pis.op(0))
+        for(auto i1 : pis.op(1))
+        for(auto i2 : pis.op(2))
+        for(auto i3 : pis.op(3)) {
+            res += i0.op(0)*i1.op(0)*i2.op(0)*i3.op(0) * LC(i0.op(1), i1.op(1), i2.op(1), i3.op(1));
+        }
+        return res;
+    }
+
+}
+

HEP/HEP.h

+/**
+ * @file 
+ * @brief HEP header file
+ */
+ 
+#pragma once
+
+#include "IBP.h"
+
+/**
+ * @brief namespace for Index, Vector, Dirac Gamma, etc.
+ */
+namespace HepLib {
+
+    using namespace std;
+    using namespace GiNaC;
+    using namespace HepLib;
+    using namespace IBP;
+    
+    extern const Symbol NA;
+    extern const Symbol NF;
+    extern const Symbol gs;
+    extern const Symbol as;
+    extern const Symbol mu;
+    extern const Symbol nL;
+    extern const Symbol nH;
+    extern exmap SP_map;
+    
+    extern int form_trace_mode;
+    extern const int form_trace_all;
+    extern const int form_trace_each_all;
+    extern const int form_trace_each_each;
+    
+    extern int form_expand_mode;
+    extern const int form_expand_none;
+    extern const int form_expand_tr;
+    extern const int form_expand_ci;
+    extern const int form_expand_li;
+    extern const int form_expand_all;
+    
+    extern bool Apart_using_fermat;
+    extern bool form_using_su3;
+    extern bool form_using_dim4;
+    
+    class Index;
+    class Vector;
+    class Pair;
+        
+    /**
+     * @brief class for FormFormat Output
+     */
+    class FormFormat : public print_dflt {
+        GINAC_DECLARE_PRINT_CONTEXT(FormFormat, print_dflt)
+    public:
+        FormFormat(ostream &os, unsigned opt=0);
+        static void power_print(const power & p, const FormFormat & c, unsigned level=0);
+        
+        template<class T> const FormFormat & operator << (const T & v) const {
+            s << v;
+            return *this;
+        };
+        const FormFormat & operator << (const basic & v) const;
+        const FormFormat & operator << (const ex & v) const;
+        const FormFormat & operator << (const lst & v) const;
+        const FormFormat & operator<<(std::ostream& (*v)(std::ostream&)) const;
+        
+        /**
+         * @brief inner class for some static initializations
+         */
+        class _init {
+            public: _init();
+        };
+    private:
+        static _init FormFormat_init;
+    };
+        
+    /**
+     * @brief class for FCFormat Output
+     */
+    class FCFormat : public print_dflt {
+        GINAC_DECLARE_PRINT_CONTEXT(FCFormat, print_dflt)
+    public:
+        FCFormat(ostream &os, unsigned opt=0);
+        static void add_print(const add & a, const FCFormat & c, unsigned level=0);
+        static void mul_print(const mul & m, const FCFormat & c, unsigned level=0);
+        static void ncmul_print(const ncmul & p, const FCFormat & c, unsigned level=0);
+        
+        template<class T> const FCFormat & operator << (const T & v) const {
+            s << v;
+            return *this;
+        };
+        const FCFormat & operator << (const basic & v) const;
+        const FCFormat & operator << (const ex & v) const;
+        const FCFormat & operator << (const lst & v) const;
+        const FCFormat & operator<<(std::ostream& (*v)(std::ostream&)) const;
+        
+        const FCFormat & operator << (const matrix & v) const;
+        const FCFormat & operator << (const exvector & v) const;
+        const FCFormat & operator << (const exmap & v) const;
+        const FCFormat & operator << (const exset & v) const;
+        
+        /**
+         * @brief inner class for some static initializations
+         */
+        class _init {
+            public: _init();
+        };
+    private:
+        static _init FCFormat_init;
+    };
+    extern FCFormat fcout;
+    
+    /**
+     * @brief class for index object
+     */
+    class Index : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(Index, basic)
+    public:
+        enum Type {VD, CF, CA};
+        Index(const string &s, const Type type=Type::VD);
+        Pair operator() (const Index & i);
+        Pair operator() (const Vector & p);
+        Symbol name;
+        Type type;
+        void print(const print_context &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static bool hasc(const ex &e);
+        static bool hasv(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(Index);
+    
+    /**
+     * @brief class for vector object
+     */
+    class Vector : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(Vector, basic)
+    public:
+        Vector(const string &s);
+        Pair operator() (const Vector & p);
+        Pair operator() (const Index & mu);
+        Symbol name;
+        void print(const print_context &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(Vector);
+    
+    /**
+     * @brief class for SUNT object
+     */
+    class SUNT : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(SUNT, basic)
+    public:
+        SUNT(ex a, ex i, ex j);
+        ex aij[3]; // Index
+        size_t nops() const override;
+        ex op(size_t i) const override;
+        ex& let_op(size_t i) override;
+        void form_print(const FormFormat &c, unsigned level = 0) const;
+        void fc_print(const FCFormat &c, unsigned level = 0) const;
+        void print(const print_dflt &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+        ex conjugate() const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(SUNT);
+    
+    /**
+     * @brief class for SUNF object
+     */
+    class SUNF : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(SUNF, basic)
+    public:
+        SUNF(ex i, ex j, ex k);
+        ex ijk[3]; // Index
+        size_t nops() const override;
+        ex op(size_t i) const override;
+        ex& let_op(size_t i) override;
+        ex eval() const override;
+        void print(const print_dflt &c, unsigned level = 0) const;
+        void form_print(const FormFormat &c, unsigned level = 0) const;
+        void fc_print(const FCFormat &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(SUNF);
+    
+    /**
+     * @brief class for SUNF4 object
+     */
+    class SUNF4 : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(SUNF4, basic)
+    public:
+        SUNF4(ex i, ex j, ex k, ex l);
+        ex ijkl[4]; // Index
+        size_t nops() const override;
+        ex op(size_t i) const override;
+        ex& let_op(size_t i) override;
+        ex eval() const override;
+        void print(const print_dflt &c, unsigned level = 0) const;
+        void form_print(const FormFormat &c, unsigned level = 0) const;
+        void fc_print(const FCFormat &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(SUNF4);
+    
+    /**
+     * @brief class for Pair object
+     */
+    class Pair : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(Pair, basic)
+    public:
+        Pair(const Vector &p1, const Vector &p2);
+        Pair(const Index &i1, const Index &i2);
+        Pair(const Vector &p, const Index &i);
+        Pair(const Index &i, const Vector &p);
+        size_t nops() const override;
+        ex op(size_t i) const override;
+        ex& let_op(size_t i) override;
+        ex eval() const override;
+        void print(const print_dflt &c, unsigned level = 0) const;
+        void form_print(const FormFormat &c, unsigned level = 0) const;
+        void fc_print(const FCFormat &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+    private:
+        ex lr[2];
+    };
+    GINAC_DECLARE_UNARCHIVER(Pair);
+    ex SP(const ex &a, bool use_map=false);
+    ex SP(const ex &a, const ex &b, bool use_map=false);
+    ex sp(const ex & a, const ex & b);
+    ex sp(const ex & a);
+    ex& letSP(const ex &p1, const ex &p2);
+    ex& letSP(const ex &p);
+    void clearSP(const ex &p1, const ex &p2);
+    void clearSP(const ex &p);
+    void clearSP();
+    ex SP2sp(const ex & exin);
+    exmap sp_map();
+    
+    /**
+     * @brief class for Levi-Civita object
+     */
+    class Eps : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(Eps, basic)
+    public:
+        ex pis[4];
+        Eps(const Vector &p1, const Vector &p2, const Vector &p3, const Vector &p4);
+        Eps(const Vector &p1, const Vector &p2, const Vector &p3, const Index &i1);
+        Eps(const Vector &p1, const Vector &p2, const Index &i1, const Index &i2);
+        Eps(const Vector &p1, const Index &i1, const Index &i2, const Index &i3);
+        Eps(const Index &i1, const Index &i2, const Index &i3, const Index &i4);
+        Eps(vector<Vector> vs, vector<Index> is);
+        size_t nops() const override;
+        ex op(size_t i) const override;
+        ex & let_op(size_t i) override;
+        ex eval() const override;
+        void print(const print_dflt &c, unsigned level = 0) const;
+        void form_print(const FormFormat &c, unsigned level = 0) const;
+        void fc_print(const FCFormat &c, unsigned level = 0) const;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(Eps);
+    ex LC(ex pi1, ex pi2, ex pi3, ex pi4);
+    
+    /**
+     * @brief class for DGamma object
+     */
+    class DGamma : public basic {
+    GINAC_DECLARE_REGISTERED_CLASS(DGamma, basic)
+    public:
+        ex pi;
+        unsigned rl;
+        DGamma(const Vector &p, unsigned rl=0);
+        DGamma(const Index &i, unsigned rl=0);
+        DGamma(int int_1567, unsigned _rl=0);
+        DGamma(const DGamma &g, unsigned _rl);
+        void print(const print_dflt &c, unsigned level = 0) const;
+        void form_print(const FormFormat &c, unsigned level = 0) const;
+        void fc_print(const FCFormat &c, unsigned level = 0) const;
+        return_type_t return_type_tinfo() const override;
+        unsigned return_type() const override { return return_types::noncommutative; }
+        bool match_same_type(const basic & other) const override;
+        unsigned get_rl();
+        size_t nops() const override;
+        ex op(size_t i) const override;
+        ex& let_op(size_t i) override;
+        ex eval() const override;
+        void archive(archive_node & n) const override;
+        void read_archive(const archive_node& n, lst& sym_lst) override;
+        static bool has(const ex &e);
+        static lst all(const ex &e);
+        ex derivative(const symbol & s) const override;
+        ex conjugate() const override;
+    };
+    GINAC_DECLARE_UNARCHIVER(DGamma);
+    
+    //-----------------------------------------------------------
+    // TR/GAS functions
+    //-----------------------------------------------------------
+    DECLARE_FUNCTION_3P(Matrix)
+    DECLARE_FUNCTION_1P(TR)
+    DECLARE_FUNCTION_1P(TTR)
+    DECLARE_FUNCTION_1P(HF)
+    
+    inline ex GAS(const Vector &p, unsigned rl=0) { return DGamma(p,rl); }
+    inline ex GAS(const Index &i, unsigned rl=0) { return DGamma(i,rl); }
+    ex GAS(const ex &expr, unsigned rl=0);
+    
+    // Form, TIR, Apart
+    ex charge_conjugate(const ex &);
+    ex form(const ex &expr, int verb=0);
+    ex TIR(const ex &expr_in, const lst &loop_ps, const lst &ext_ps);
+    ex MatrixContract(const ex & expr_in);
+    ex Apart(const ex &expr_in, const lst &vars, exmap sgnmap={});
+    ex Apart(const ex &expr_in, const lst &loops, const lst & extmoms);
+    ex ApartIR2ex(const ex & expr_in);
+    ex ApartIR2F(const ex & expr_in);
+    ex F2ex(const ex & expr_in);
+    ex ApartIRC(const ex & expr_in);
+    void ApartIBP(int IBPmethod, exvector &io_vec, const lst & loops_exts=lst{}, const lst & cut_props=lst{}, std::function<lst(const Base &, const ex &)> uf=IBP::LoopUF);
+    
+    /**
+     * @brief ApartIR function with 1 argument
+     */
+    class ApartIR1_SERIAL { public: static unsigned serial; };
+    template<typename T1>
+    inline GiNaC::function ApartIR(const T1 & p1) {
+        return GiNaC::function(ApartIR1_SERIAL::serial, ex(p1));
+    }
+    
+    /**
+     * @brief ApartIR function with 2 arguments
+     */
+    class ApartIR2_SERIAL { public: static unsigned serial; };
+    template<typename T1, typename T2>
+    inline GiNaC::function ApartIR(const T1 & p1, const T2 & p2) {
+        return GiNaC::function(ApartIR2_SERIAL::serial, ex(p1), ex(p2));
+    }
+        
+        
+}
+

HEP/Pair.cpp

+/**
+ * @file
+ * @brief Pair class
+ */
+ 
+#include "HEP.h"
+
+namespace HepLib {
+
+    //-----------------------------------------------------------
+    // Pair Class
+    //-----------------------------------------------------------
+    GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(Pair, basic,
+        print_func<print_dflt>(&Pair::print).
+        print_func<FormFormat>(&Pair::form_print).
+        print_func<FCFormat>(&Pair::fc_print)
+    )
+    
+    DEFAULT_CTOR(Pair)
+    GINAC_BIND_UNARCHIVER(Pair);
+    IMPLEMENT_HAS(Pair)
+    IMPLEMENT_ALL(Pair)
+
+    /**
+     * @brief Pair constructor, ordered internally
+     * @param p1 Vector in p1.p2
+     * @param p2 Vector  in p1.p2
+     * @return p1.p2
+     */
+    Pair::Pair(const Vector &p1, const Vector &p2) { 
+        if(p1.compare(p2)>0) {
+            lr[0]=p1;lr[1]=p2; 
+        } else {
+            lr[0]=p2;lr[1]=p1; 
+        }
+    }
+    
+    /**
+     * @brief Pair constructor, ordered internally
+     * @param i1 Index in i1.i2
+     * @param i2 Index  in i1.i2
+     * @return i1.i2 i.e., g^{i1,i2}
+     */
+    Pair::Pair(const Index &i1, const Index &i2) { 
+        if(i1.compare(i2)>0) {
+            lr[0]=i1;lr[1]=i2; 
+        } else {
+            lr[0]=i2;lr[1]=i1; 
+        }
+    }
+    
+    /**
+     * @brief Pair constructor
+     * @param p Vector in p^i
+     * @param i Index  in p^i
+     * @return p^i
+     */
+    Pair::Pair(const Vector &p, const Index &i)  : lr{p, i} { }
+    
+    /**
+     * @brief Pair constructor
+     * @param i Index  in p^i
+     * @param p Vector in p^i
+     * @return p^i
+     */
+    Pair::Pair(const Index &i, const Vector &p) : lr{p, i} { }
+
+    int Pair::compare_same_type(const basic &other) const {
+        const Pair &o = static_cast<const Pair &>(other);
+        int c1 = lr[0].compare(o.lr[0]);
+        if(c1!=0) return c1;
+        int c2 = lr[1].compare(o.lr[1]);
+        return c2;
+    }
+    
+    /**
+     * @brief default print function
+     * @param c default print format
+     * @param level level in print function
+     */
+    void Pair::print(const print_dflt &c, unsigned level) const {
+        c.s << lr[0] << "." << lr[1];
+    }
+    
+    /**
+     * @brief FormFormat print function
+     * @param c FormFormat print format
+     * @param level level in print function
+     */
+    void Pair::form_print(const FormFormat &c, unsigned level) const {
+        if(is_a<Vector>(lr[0]) && is_a<Vector>(lr[1])) c << lr[0] << "." << lr[1];
+        else if(is_a<Vector>(lr[0]) && is_a<Index>(lr[1])) c << lr[0] << "(" << lr[1] << ")";
+        else if(is_a<Index>(lr[0]) && is_a<Index>(lr[1])) c << "d_(" << lr[0] << "," << lr[1] << ")";
+    }
+    
+    /**
+     * @brief FCFormat print function
+     * @param c FCFormat print format
+     * @param level level in print function
+     */
+    void Pair::fc_print(const FCFormat &c, unsigned level) const {
+        if(is_a<Vector>(lr[0]) && is_a<Vector>(lr[1])) c << "SPD(" << lr[0] << "," << lr[1] << ")";
+        else if(is_a<Vector>(lr[0]) && is_a<Index>(lr[1])) c << "FVD(" << lr[0] << "," << lr[1] << ")";
+        else if(is_a<Index>(lr[0]) && is_a<Index>(lr[1])) {
+            auto ii = ex_to<Index>(lr[0]);
+            if(ii.type == Index::Type::VD) c << "MTD(" << lr[0] << "," << lr[1] << ")";
+            else if(ii.type == Index::Type::CF) c << "Delta(" << lr[0] << "," << lr[1] << ")";
+            else if(ii.type == Index::Type::CA) c << "SUNDelta(" << lr[0] << "," << lr[1] << ")";
+            else throw Error("Pair::fc_print unexpected.");
+        }
+    }
+    
+    size_t Pair::nops() const { return 2; }
+    ex Pair::op(size_t i) const {
+        return lr[i];
+    }
+    ex & Pair::let_op(size_t i) {
+        ensure_if_modifiable();
+        return lr[i];
+    }
+    
+    ex Pair::eval() const {
+        if(flags & status_flags::evaluated) return *this;
+        else if(!is_a<Vector>(lr[0]) && !is_a<Index>(lr[0])) return SP(lr[0],lr[1]);
+        else if(!is_a<Vector>(lr[1]) && !is_a<Index>(lr[1])) return SP(lr[0],lr[1]);
+        else if((is_a<Vector>(lr[0]) && is_a<Vector>(lr[1]))) 
+            return Pair(ex_to<Vector>(lr[0]),ex_to<Vector>(lr[1])).hold();
+        else if((is_a<Index>(lr[0]) && is_a<Index>(lr[1]))) 
+            return Pair(ex_to<Index>(lr[0]),ex_to<Index>(lr[1])).hold();
+        else if((is_a<Index>(lr[0]) && is_a<Vector>(lr[1])))
+            return Pair(ex_to<Vector>(lr[1]),ex_to<Index>(lr[0])).hold();
+        else return this->hold();
+    }
+    
+    void Pair::archive(archive_node & n) const {
+        inherited::archive(n);
+        for(int i=0; i<2; i++) n.add_ex("lr"+to_string(i), lr[i]);
+    }
+    
+    void Pair::read_archive(const archive_node& n, lst& sym_lst) {
+        inherited::read_archive(n, sym_lst);
+        for(int i=0; i<2; i++) {
+            n.find_ex("lr"+to_string(i), lr[i], sym_lst);
+        }
+    }
+    
+    ex Pair::derivative(const symbol & s) const {
+        return 0;
+    }
+    
+    //-----------------------------------------------------------
+    // SP function - ScalarProduct
+    //-----------------------------------------------------------
+    ex SP(const ex & a, bool use_map) { return SP(a,a,use_map); }
+    
+    /**
+     * @brief Function similar to SPD/FVD in FeynCalc
+     * @param a the 1st vector/index
+     * @param b the 2nd vector/index
+     * @param use_map true for apply the SP_map in the function call
+     * @return expanded/translated Pair objects
+     */
+    ex SP(const ex & a, const ex & b, bool use_map) {
+        ex ret_sp;
+        if(is_a<Vector>(a) && is_a<Vector>(b)) ret_sp = Pair(ex_to<Vector>(a), ex_to<Vector>(b));
+        else if(is_a<Vector>(a) && is_a<Index>(b)) ret_sp = Pair(ex_to<Vector>(a), ex_to<Index>(b));
+        else if(is_a<Index>(a) && is_a<Vector>(b)) ret_sp = Pair(ex_to<Vector>(b), ex_to<Index>(a));
+        else if(is_a<Index>(a) && is_a<Index>(b)) ret_sp = Pair(ex_to<Index>(a), ex_to<Index>(b));
+        if(is_a<Pair>(ret_sp)) {
+            if(use_map) return ret_sp.subs(SP_map);
+            else return ret_sp;
+        }
+
+        lst alst, blst;
+        auto aex = a.expand();
+        auto bex = b.expand();
+        if(is_zero(aex) || is_zero(bex)) return 0;
+        if(is_a<add>(aex)) {
+            for(auto item : aex) alst.append(item);
+        } else alst.append(aex);
+        for(int i=0; i<alst.nops(); i++) {
+            if(!is_a<mul>(alst.op(i))) alst.let_op(i) = lst{alst.op(i)};
+            ex c=1;
+            ex v=1;
+            for(auto ii : alst.op(i)) {
+                if(is_a<Vector>(ii) || is_a<Index>(ii)) {
+                    if(!is_zero(v-1)) throw Error("Error Found in SP @1 : "+ex2str(alst));
+                    v = ii;
+                } else c *= ii;
+            }
+            if(is_zero(v-1)) throw Error("Error Found in SP @2 : "+ex2str(alst));
+            alst.let_op(i) = lst{c,v};
+        }
+        
+        if(is_a<add>(bex)) {
+            for(auto item : bex) blst.append(item);
+        } else blst.append(bex);
+        for(int i=0; i<blst.nops(); i++) {
+            if(!is_a<mul>(blst.op(i))) blst.let_op(i) = lst{blst.op(i)};
+            ex c=1;
+            ex v=1;
+            for(auto ii : blst.op(i)) {
+                if(is_a<Vector>(ii) || is_a<Index>(ii)) {
+                    if(!is_zero(v-1)) throw Error("Error Found in SP @3");
+                    v = ii;
+                } else c *= ii;
+            }
+            if(is_zero(v-1)) throw Error("Error Found in SP @4");
+            blst.let_op(i) = lst{c,v};
+        }
+        
+        ex res = 0;
+        for(auto ai : alst) {
+            for(auto bi : blst) res += ai.op(0) * bi.op(0) * SP(ai.op(1), bi.op(1), use_map);
+        }
+        return res;
+    }
+    
+    /**
+     * @brief translated the vector dot a.b to a*b, useful in SecDec
+     * @param a the 1st vector argument or a symbol
+     * @param b the 2nd vector or a symbol
+     * @return expression with vector name left only
+     */
+    ex sp(const ex & a, const ex & b) {
+        if(is_a<Vector>(a) && is_a<Vector>(b)) return ex_to<Vector>(a).name * ex_to<Vector>(b).name;
+        else throw Error("Error in sp(2).");
+    }
+    ex sp(const ex & a) {
+        if(is_a<Vector>(a)) return ex_to<Vector>(a).name * ex_to<Vector>(a).name;
+        else throw Error("Error in sp(1).");
+    }
+    
+    /**
+     * @brief return the reference of p1.p2
+     * @param p1 the 1st vector/index
+     * @param p2 the 2nd vector/index
+     * @return the refernce of p1.p2 in SP_map, can be used in assignment
+     */
+    ex & letSP(const ex &p1, const ex &p2) {
+        if(!(is_a<Vector>(p1) || is_a<Index>(p1)) || !(is_a<Vector>(p2) || is_a<Index>(p2)))
+            throw Error("Invalide arguments for letSP (2).");
+        return SP_map[SP(p1,p2)];
+    }
+    
+    /**
+     * @brief return the reference of p.p
+     * @param p the vector/index
+     * @return the refernce of p.p in SP_map, can be used in assignment
+     */
+    ex & letSP(const ex &p) {
+        if(!is_a<Vector>(p))
+            throw Error("Invalide arguments for letSP (1).");
+        return SP_map[SP(p)];
+    }
+    
+    /**
+     * @brief delete the assignment of p1.p2
+     * @param p1 the 1st vector/index
+     * @param p2 the 2nd vector/index
+     */
+    void clearSP(const ex &p1, const ex &p2) {
+        if(!(is_a<Vector>(p1) || is_a<Index>(p1)) || !(is_a<Vector>(p2) || is_a<Index>(p2)))
+            throw Error("Invalide arguments for resetSP (2).");
+        SP_map.erase(SP(p1,p2));
+    }
+    
+    /**
+     * @brief delete the assignment of p.p
+     * @param p vector/index
+     */
+    void clearSP(const ex &p) {
+        if(!is_a<Vector>(p))
+            throw Error("Invalide arguments for resetSP (1).");
+        SP_map.erase(SP(p));
+    }
+    
+    /**
+     * @brief delete all assignment in SP_map
+     */
+    void clearSP() {
+        SP_map.clear();
+    }
+    
+    /**
+     * @brief convert SP(a,b) to sp(a,b)
+     * @param exin the input expression
+     * @return expression with SP(a,b) replaced to sp(a,b)
+     */
+    ex SP2sp(const ex & exin) {
+        auto ret = exin.subs(SP_map);
+        lst vecs = Vector::all(ret);
+        exmap spmap;
+        for(auto vp1 : vecs) {
+            for(auto vp2 : vecs) {
+                spmap[SP(vp1,vp2)]=sp(vp1,vp2);
+            }
+        }
+        return ret.subs(spmap);
+    }
+    
+    /**
+     * @brief the SP_map with SP(a,b) replaced to sp(a,b)
+     * @return the SP_map with SP(a,b) replaced to sp(a,b)
+     */
+    exmap sp_map() {
+        exmap ret;
+        for(auto kv : SP_map) {
+            auto key = kv.first;
+            if(key.nops()!=2 || !is_a<Vector>(key.op(0)) || !is_a<Vector>(key.op(1))) continue;
+            auto p1 = ex_to<Vector>(key.op(0));
+            auto p2 = ex_to<Vector>(key.op(1));
+            ret[p1.name * p2.name] = kv.second.subs(SP_map);
+        }
+        return ret;
+    }
+
+}
+

HEP/TIR.cpp

+/**
+ * @file
+ * @brief Functions to perform Tensor Index Reduction
+ */
+ 
+#include "HEP.h"
+#include <cmath>
+
+namespace HepLib {
+
+    namespace {
+        
+        lst combs(const lst & items, int n) {
+            lst ret;
+            if(n<1 || items.nops()<1) return lst{ret};
+            if(n==1) {
+                for(auto it : items) ret.append(lst{it});
+                return ret;
+            }
+            if(items.nops()==1) {
+                lst comb;
+                for(int i=0; i<n; i++) comb.append(items.op(0));
+                ret.append(comb);
+                return ret;
+            }
+            auto its = items;
+            ex first = its.op(0);
+            its.remove_first();
+            //for(int i=0; i<=n; i++) {
+            for(int i=n; i>=0; i--) {
+                auto rets = combs(its, n-i);
+                for(auto it : rets) {
+                    lst item = ex_to<lst>(it);
+                    for(int j=0; j<i; j++) item.append(first);
+                    ret.append(item);
+                }
+            }
+            return ret;
+        }
+        
+    }
+
+    /**
+     * @brief Tensor Index Reduction
+     * @param expr_in expression 
+     * @param loop_ps lst contains loop vectors
+     * @param ext_ps lst constains external vectors
+     * @return TIR result
+     */
+    ex TIR(const ex &expr_in, const lst &loop_ps, const lst &ext_ps) {
+        for(auto pi : loop_ps) {
+            if(!is_a<Vector>(pi)) throw Error("TIR invalid 2nd argument");
+        }
+        for(auto pi : ext_ps) {
+            if(!is_a<Vector>(pi)) throw Error("TIR invalid 3rd argument");
+        }
+
+        if(expr_in.has(coVF(w))) throw Error("TIR error: expr_in has coVF already.");
+
+        auto expr = mma_collect(expr_in, [&](const ex & e)->bool{
+            if(!Index::hasv(e)) return false;
+            for(const_preorder_iterator i = e.preorder_begin(); i != e.preorder_end(); ++i) {
+                auto item = *i;
+                if(is_a<Pair>(item) && is_a<Index>(item.op(1)) && loop_ps.has(item.op(0))) return true;
+            }
+            return false;
+        }, false, true);
+
+        static exmap cache_map;
+        expr = MapFunction([ext_ps,loop_ps](const ex &e, MapFunction &self)->ex{
+            if(e.is_equal(coVF(1))) return 1;
+            else if(e.match(coVF(w))) {
+                ex map_key = lst{e,loop_ps,ext_ps};
+                if(cache_map.find(map_key)!=cache_map.end()) return cache_map[map_key];
+                lst vis, lps;
+                map<ex,int,ex_is_less> pc;
+                if(is_a<mul>(e.op(0))) {
+                    for(auto item : e.op(0)) {
+                        vis.append(item);
+                        lps.append(item.op(0));
+                        pc[item.op(0)]++;
+                    }
+                } else {
+                    vis.append(e.op(0));
+                    lps.append(e.op(0).op(0));
+                }
+                lps.sort();
+                lps.unique();
+                
+                auto visn = vis.nops();
+                if(lps.nops()<2) {
+                    ex eqL=1;
+                    lst is, bis;
+                    for(auto vi : vis) {
+                        eqL *= vi;
+                        is.append(vi.op(1));
+                    }
+
+                    //for(int er=((visn%2==0)?0:1); er<=visn; er+=2) {
+                    for(int er=visn; er>=((visn%2==0)?0:1); er-=2) {
+                        auto ep_comb = combs(ext_ps, er);
+                        for(auto item : ep_comb) {
+                            ex bi=1;
+                            for(int j=0; j<er; j++) bi *= SP(item.op(j), is.op(j));
+                            for(int j=er; j<visn; j=j+2) bi *= SP(is.op(j), is.op(j+1));
+                            bi = bi.symmetrize(is);
+                            bis.append(bi);
+                    }}
+                    
+                    int n = bis.nops();
+                    lst mat;
+                    for(auto bi : bis) {
+                        for(auto bj : bis) mat.append(bi*bj);
+                    }
+                    for(auto bj : bis) mat.append(eqL*bj);
+                    
+                    // we need to remap the dummy index, to avoid SP_map set the index object to 0
+                    if(true) {
+                        lst isu = is;
+                        isu.sort();
+                        isu.unique();
+                        exmap i2u, u2i;
+                        int c=0;
+                        for(auto item : isu) {
+                            auto ii = Index("TIR"+to_string(c++));
+                            i2u[item] = ii;
+                            u2i[ii] = item;
+                        }
+                        mat = ex_to<lst>(subs(mat,i2u));
+                        mat = ex_to<lst>(form(mat).subs(u2i));
+                    }
+
+                    lst rep_vs;
+                    ex tree = mat;
+                    for(const_preorder_iterator i = tree.preorder_begin(); i != tree.preorder_end(); ++i) {
+                        auto e = (*i);
+                        if(is_a<symbol>(e) || is_a<Pair>(e) || is_a<Eps>(e)) {
+                            rep_vs.append(e);
+                        }
+                    }
+                    rep_vs.sort();
+                    rep_vs.unique();
+                    sort_lst(rep_vs);
+                    
+                    exmap v2f;
+                    symtab st;
+                    int fvi = 0;
+                    for(auto vi : rep_vs) {
+                        auto name = "v" + to_string(fvi);
+                        v2f[vi] = Symbol(name);
+                        st[name] = vi;
+                        fvi++;
+                    }
+
+                    stringstream ss;
+                    for(int i=0; i<fvi; i++) ss << "&(J=v" << i << ");" << endl;
+                    Fermat fermat;
+                    fermat.Init();
+                    fermat.Execute(ss.str());
+                    ss.clear();
+                    ss.str("");
+                    
+                    ss << "Array m[" << n << "," << n+1 << "];" << endl;
+                    fermat.Execute(ss.str());
+                    ss.clear();
+                    ss.str("");
+                    
+                    ss << "[m]:=[(";
+                    for(auto item : mat) {
+                        ss << item.subs(v2f) << ",";
+                    }
+                    ss << ")];" << endl;
+                    ss << "Redrowech([m]);" << endl;
+                    auto tmp = ss.str();
+
+                    string_replace_all(tmp,",)]",")]");
+                    fermat.Execute(tmp);
+                    ss.clear();
+                    ss.str("");
+
+                    ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
+                    ss << "![m" << endl;
+                    auto ostr = fermat.Execute(ss.str());
+                    fermat.Exit();
+                    
+                    // make sure last char is 0
+                    if(ostr[ostr.length()-1]!='0') throw Error("TIR: last char is NOT 0.");
+                    ostr = ostr.substr(0, ostr.length()-1);
+                    string_trim(ostr);
+                    
+                    ostr.erase(0, ostr.find(":=")+2);
+                    string_replace_all(ostr, "[", "{");
+                    string_replace_all(ostr, "]", "}");
+                    Parser fp(st);
+                    auto mat2 = fp.Read(ostr);
+
+                    ex res = 0;
+                    for(int i=0; i<n; i++) {
+                        if(is_zero(mat2.op(i).op(i)) && !is_zero(mat2.op(i).op(n))) throw Error("Zero Determinant in TIR.");
+                        res += bis.op(i) * mat2.op(i).op(n);
+                    }
+                    res = res.subs(SP_map);
+                    cache_map[map_key] = res;
+                    return res;
+                } else {
+                    int cmin=10000, cmax=-1;
+                    ex pmin, pmax;
+                    for(auto lpi : lps) {
+                        auto cc = pc[lpi];
+                        if(cc<cmin) { cmin = cc; pmin = lpi; }
+                        if(cc>cmax) { cmax = cc; pmax = lpi; }
+                    }
+                    
+                    auto lp0 = pmin;
+                    auto ext_ps2 = ext_ps;
+                    for(auto lpi : lps)
+                        if(!is_zero(lpi-lp0)) ext_ps2.append(lpi);
+                    ex ret = TIR(e.op(0), lst{ lp0 }, ext_ps2);
+                    ret = TIR(ret, loop_ps, ext_ps);
+                    cache_map[map_key] = ret;
+                    return ret;
+                }
+            } else if (!e.has(coVF(w))) return e;
+            else return e.map(self);
+        })(expr);
+
+        return expr;
+        
+    }
+
+}
+