Library makefilakefile modified --> LINKFLAG without OPTI

csl/include/libmakefile_data.h

@@ -49,10 +49,17 @@ void print_libmakefile_data(std::ostream &out, bool quad, bool clang) {
                 "\t\t-freciprocal-math\n\n";
 #endif
     }
-    if (!quad)
-      out << "CFLAGS = -std=c++17 -O2 $(MATH_OPTI) -Wall -Wextra -Wpedantic -Wno-deprecated-declarations -fPIC\n";
-    else
-      out << "CFLAGS = -std=c++17 -O2 $(MATH_OPTI) -Wall -Wextra -Wpedantic -Wno-deprecated-declarations -fPIC -DQUAD=1 -DQUADSIZE=16\n";
+    out << "DEFAULTFLAGS = -std=c++17 -Wall -Wextra -Wpedantic -Wno-deprecated-declarations -fPIC\n";
+    out << "OPTIFLAGS = -O2 $(MATH_OPTI)\n";
+    if (!quad) {
+      out << "CXXFLAGS  = $(DEFAULTFLAGS) $(OPTIFLAGS)\n";
+      out << "LINKFLAGS = $(DEFAULTFLAGS)\n";
+    }
+    else {
+      out << "QUADFLAGS = -DQUAD=1 -DQUADSIZE=16\n";
+      out << "CXXFLAGS  = $(DEFAULTFLAGS) $(OPTIFLAGS) $(QUADFLAGS)\n";
+      out << "LINKFLAGS = $(DEFAULTFLAGS) $(QUADFLAGS)\n";
+    }
     out << "\n";
     out << "SRCDIR = src\n";
     out << "INCDIR = include\n";
@@ -76,16 +83,16 @@ void print_libmakefile_data(std::ostream &out, bool quad, bool clang) {
     out << "lib: $(NAMELIB).so $(NAMELIB).a\n";
     out << "\n";
     out << "$(OBJDIR)/%.o: $(SRCDIR)/%.cpp\n";
-    out << "	$(CXX) $(CFLAGS) -c $< -o $@ $(INCPATH)\n";
+    out << "	$(CXX) $(CXXFLAGS) -c $< -o $@ $(INCPATH)\n";
     out << "\n";
     out << "$(SOBJDIR)/%.o: $(SCRIPTDIR)/%.cpp\n";
-    out << "	$(CXX) $(CFLAGS) -c $< -o $@ $(INCPATH)\n";
+    out << "	$(CXX) $(CXXFLAGS) -c $< -o $@ $(INCPATH)\n";
     out << "\n";
     out << "%.x: $(SOBJDIR)/%.o $(OBJ)\n";
-    out << "	$(CXX) $(CFLAGS) -o $(BINDIR)/$@ $< $(OBJ) $(INCPATH) $(LIBPATH) $(LIBS)\n";
+    out << "	$(CXX) $(LINKFLAGS) -o $(BINDIR)/$@ $< $(OBJ) $(INCPATH) $(LIBPATH) $(LIBS)\n";
     out << "\n";
     out << "$(NAMELIB).so: $(OBJ)\n";
-    out << "	$(CXX) $(CFLAGS) -shared -o $(LIBDIR)/$@ $(OBJ) $(LIBPATH) $(LIBS)\n";
+    out << "	$(CXX) $(CXXFLAGS) -shared -o $(LIBDIR)/$@ $(OBJ) $(LIBPATH) $(LIBS)\n";
     out << "$(NAMELIB).a: $(OBJ)\n";
     out << "	ar rcs $(LIBDIR)/$@ $(OBJ)\n";
     out << "\n";

csl/src/operations.cpp

@@ -2776,6 +2776,16 @@ bool Pow::mergeTerms()
             argument[1] = CSL_1;
         }
     }
+    else if (argument[0]->getType() == csl::Type::IntFraction
+            and argument[1]->isInteger()
+            and std::fabs(argument[1]->evaluateScalar()) < 4) {
+        const long long num = round(
+                csl::pow_s(argument[0]->getNum(), argument[1])->evaluateScalar());
+        const long long denom = round(
+                csl::pow_s(argument[0]->getDenom(), argument[1])->evaluateScalar());
+        argument[0] = csl::intfraction_s(num, denom);
+        argument[1] = CSL_1;
+    }
     else if (argument[0]->getPrimaryType() == csl::PrimaryType::Numerical and
             argument[1]->getType() == csl::Type::IntFraction and
             std::abs(argument[1]->getNum()) == 1 and

marty/include/interactionTerm.h

@@ -408,7 +408,9 @@ class InteractionTerm {
     std::vector<permutation> permutations;
 };
 
+int matchBOnA(csl::Expr const& A, csl::Expr &B);
 bool hardComparison(csl::Expr const&, csl::Expr const&);
+bool hardOrdering(csl::Expr const&, csl::Expr const&);
 
 }
 

marty/include/model.h

@@ -241,7 +241,7 @@ public:
             std::pair<size_t, size_t> const &treeCoupling,
             std::pair<size_t, size_t> const &loopCoupling,
             Insertion                 const &mediator,
-            FeynOptions               const &feynOptions = {}
+            FeynOptions                      feynOptions = {}
             );
 
     WilsonSet computeWilsonPenguins_4Fermions(

marty/include/wilson.h

@@ -78,6 +78,7 @@ class WilsonOperator {
     static bool hardComparison(csl::Expr const& A,
                                csl::Expr const& B);
 
+    bool operator< (WilsonOperator const& other) const;
     bool operator==(WilsonOperator const& other) const;
 
     private:
@@ -129,6 +130,7 @@ struct WilsonSet: public std::vector<Wilson> {
     }
 
     void merge();
+    void sort();
 
     Kinematics  kinematics {};
     FeynOptions options {};

marty/src/dimensionalRegularization.cpp

@@ -79,6 +79,7 @@ namespace mty::simpli {
                 }
             }
         });
+        csl::Replace(res, csl::DMinko, csl::int_s(4));
     }
 
 }

marty/src/insertion.cpp

@@ -158,7 +158,7 @@ Insertion Mediator(Insertion const& init)
 {
     Insertion other(init);
     other.setMediator(true);
-    return other;
+    return OffShell(other);
 }
 
 Insertion Left(Insertion const &init)

marty/src/interactionTerm.cpp

@@ -806,9 +806,7 @@ ostream& operator<<(ostream              & out,
     return out;
 }
 
-static bool hardComparison_impl(
-        csl::Expr const &A,
-        csl::Expr       &B)
+int matchBOnA(csl::Expr const& A, csl::Expr &B)
 {
     std::vector<csl::Expr> tensorsInA;
     std::vector<csl::Expr> tensorsInB;
@@ -823,7 +821,7 @@ static bool hardComparison_impl(
             tensorsInB.push_back(el);
     });
     if (tensorsInA.size() != tensorsInB.size()) {
-        return false;
+        return tensorsInA.size() < tensorsInB.size();
     }
     std::sort(tensorsInA.begin(), tensorsInA.end());
     std::sort(tensorsInB.begin(), tensorsInB.end());
@@ -831,7 +829,7 @@ static bool hardComparison_impl(
     for (size_t i = tensorsInA.size(); i --> 0 ;)
         if (tensorsInA[i]->getParent_info() 
                 != tensorsInB[i]->getParent_info()) {
-            return false;
+            return tensorsInA[i]->getName() < tensorsInB[i]->getName();
         }
         else {
             csl::IndexStructure Astruct = tensorsInA[i]->getIndexStructure();
@@ -890,8 +888,17 @@ static bool hardComparison_impl(
         ++index;
     }
     csl::DeepRefresh(B);
-    const auto res = A->compareWithDummy(B.get());
-    return res;
+    return -1;
+}
+
+static bool hardComparison_impl(
+        csl::Expr const &A,
+        csl::Expr       &B)
+{
+    const int match = matchBOnA(A, B);
+    if (match != -1)
+        return false;
+    return A->compareWithDummy(B.get());
 }
 
 bool hardComparison(
@@ -904,4 +911,24 @@ bool hardComparison(
     return hardComparison_impl(A, B_renameIndices);
 }
 
+static bool hardOrdering_impl(
+        csl::Expr const &A,
+        csl::Expr       &B)
+{
+    const int match = matchBOnA(A, B);
+    if (match != -1)
+        return match;
+    return A < B;
+}
+
+bool hardOrdering(
+        csl::Expr const& A,
+        csl::Expr const& B
+        )
+{
+    auto B_renameIndices = csl::DeepCopy(B);
+    csl::RenameIndices(B_renameIndices);
+    return hardOrdering_impl(A, B_renameIndices);
+}
+
 } // End of namespace mty

marty/src/model.cpp

@@ -351,6 +351,7 @@ csl::Expr Model::computeSquaredAmplitude(
                     csl::GetComplexConjugate(ampl[i].coef.getCoefficient()), 
                     ampl[j].coef.getCoefficient()
                     });
+            csl::Replace(prod, csl::DMinko, 4);
             res.push_back(prod);
             if (i != j) {
                 res.push_back(csl::GetHermitianConjugate(prod, &dirac4));
@@ -439,6 +440,7 @@ WilsonSet Model::getWilsonCoefficients(
         w.coef.setCoefficient(a);
     }
     WilsonSet res(wilsons.begin(), wilsons.end());
+    res.sort();
     res.options    = ampl.getOptions();
     res.kinematics = ampl.getKinematics();
 
@@ -521,39 +523,42 @@ WilsonSet Model::computeSingleWilsonPenguin_4Fermions(
         std::pair<size_t, size_t> const &treeCoupling,
         std::pair<size_t, size_t> const &loopCoupling,
         Insertion                 const &mediator,
-        FeynOptions               const &feynOptions
+        FeynOptions                      feynOptions
         )
 {
     csl::ScopedProperty verbose(&mty::option::verboseAmplitude, false);
     auto const &insertions = kinematics.getInsertions();
     std::vector<Insertion> treeInsertions = {
-        OffShell(insertions[treeCoupling.first]), 
-        OffShell(insertions[treeCoupling.second]),
-        OffShell(Mediator(mediator))
+        insertions[treeCoupling.first], 
+        insertions[treeCoupling.second],
+        Mediator(mediator)
     };
     std::vector<Insertion> loopInsertions = {
-        OffShell(insertions[loopCoupling.first]), 
-        OffShell(insertions[loopCoupling.second]),
+        insertions[loopCoupling.first], 
+        insertions[loopCoupling.second],
         mediator.isIncoming() ? 
-            OffShell(Mediator(Outgoing(mediator)))
-            : OffShell(Mediator(Incoming(mediator)))
+            Mediator(Outgoing(mediator))
+            : Mediator(Incoming(mediator))
     };
     auto treeAmplitude = computePartialAmplitude(
             TreeLevel, treeInsertions, feynOptions);
     if (treeAmplitude.empty()) {
         return {};
     }
+    const bool massless = mediator.getField()->getMass() == CSL_0;
     csl::ScopedProperty prop(
             &mty::option::keepOnlyFirstMassInLoop, 
-            mediator.getField()->getMass() != CSL_0 
-                && mty::option::useMassiveSimplifications
+            !massless && mty::option::useMassiveSimplifications
             );
     auto loopAmplitude = computePartialAmplitude(
             OneLoop, loopInsertions, feynOptions);
     if (loopAmplitude.empty()) {
         return {};
     }
-
+    if (massless) {
+        auto wil = getWilsonCoefficients(loopAmplitude);
+        Display(wil);
+    }
     Amplitude connexion = connectAmplitudes(
             treeAmplitude, loopAmplitude, feynOptions);
     feynOptions.applyFilters(connexion.getDiagrams(), true);
@@ -579,7 +584,6 @@ WilsonSet Model::computeWilsonPenguins_4Fermions(
     }};
     std::vector<Amplitude> amplitudes;
     amplitudes.reserve(bosons.size());
-    feynOptions.setTopology(Topology::Triangle | Topology::Mass);
     WilsonSet res;
     for (const mty::Particle &mediator : bosons) {
         for (const auto &[first, second, other1, other2] : pairs) {
@@ -639,6 +643,7 @@ WilsonSet Model::computeWilsonCoefficients_4Fermions(
         addWilson(wil, res, false);
 
     res.merge();
+    res.sort();
 
     return res;
 }

marty/src/wilson.cpp

@@ -150,6 +150,12 @@ bool WilsonOperator::hardComparison(csl::Expr const& A,
     return mty::hardComparison(A, B);
 }
 
+bool WilsonOperator::operator<(WilsonOperator const& other) const
+{
+    return hardOrdering(DeepRefreshed(op), 
+                        DeepRefreshed(other.op));
+}
+
 bool WilsonOperator::operator==(WilsonOperator const& other) const
 {
     return hardComparison(DeepRefreshed(op), 
@@ -175,6 +181,13 @@ void WilsonSet::merge()
     erase(last, end());
 }
 
+void WilsonSet::sort()
+{
+    std::sort(begin(), end(), [&](Wilson const &A, Wilson const &B) {
+        return A.op < B.op;
+    });
+}
+
 ///////////////////////////////////////////////////
 /*************************************************/
 // Other functions                               //