Advanced Computing Platform for Theoretical Physics

Skip to content

GitLab

Commit 2e5a5f65 authored by Feng Feng's avatar Feng Feng
Browse files

cleanup

parent 8abf1d7e
Expand all Hide whitespace changes
Inline Side-by-side
BAK/CMakeLists.txt deleted 100644 → 0
View file @ 8abf1d7e
# usage
# cmake -DCMAKE_INSTALL_PREFIX=p0 -DINC_PATH="p1;p2" -DLIB_PATH="p3;p4" path
CMAKE_MINIMUM_REQUIRED(VERSION 3.10)
SET(CMAKE_C_COMPILER gcc)
SET(CMAKE_CXX_COMPILER g++)
SET(CMAKE_FORTRAN_COMPILER gfortran)
#-----------------------------------------------
# qgraf
#-----------------------------------------------
ADD_EXECUTABLE(qgraf "bin/qgraf-3.1.4.f")
#-----------------------------------------------
# INC_PATH & LIB_PATH
#-----------------------------------------------
SET(INC_PATH "${INC_PATH};${CMAKE_INSTALL_PREFIX}/include")
LIST(REMOVE_DUPLICATES INC_PATH)
SET(INC_FLAGS "")
FOREACH(dir ${INC_PATH})
IF(EXISTS ${dir})
INCLUDE_DIRECTORIES(${dir})
SET(INC_FLAGS "${INC_FLAGS} -I'${dir}'")
ENDIF()
ENDFOREACH()
STRING(FIND "${INC_FLAGS}" "${CMAKE_INSTALL_PREFIX}/include" INC_POS)
IF( "${INC_POS}" STREQUAL "-1" )
SET(INC_FLAGS "${INC_FLAGS} -I'${CMAKE_INSTALL_PREFIX}/include'")
ENDIF()
SET(LIB_PATH "${LIB_PATH};${CMAKE_INSTALL_PREFIX}/lib")
LIST(REMOVE_DUPLICATES LIB_PATH)
SET(LIB_FLAGS "-Wl,-rpath,.")
FOREACH(dir ${LIB_PATH})
IF(EXISTS ${dir})
LINK_DIRECTORIES(${dir})
SET(LIB_FLAGS "${LIB_FLAGS} -Wl,-rpath,'${dir}' -L'${dir}'")
ENDIF()
ENDFOREACH()
STRING(FIND "${LIB_FLAGS}" "${CMAKE_INSTALL_PREFIX}/lib" LIB_POS)
IF( "${LIB_POS}" STREQUAL "-1" )
SET(LIB_FLAGS "${LIB_FLAGS} -Wl,-rpath,'${CMAKE_INSTALL_PREFIX}/lib' -L'${CMAKE_INSTALL_PREFIX}/lib'")
ENDIF()
INCLUDE_DIRECTORIES(BEFORE include)
CONFIGURE_FILE(include/Init.cpp "${CMAKE_CURRENT_BINARY_DIR}/OTHER/Init.cpp" @ONLY)
CONFIGURE_FILE(include/HepLib.pc "${CMAKE_CURRENT_BINARY_DIR}/OTHER/HepLib.pc" @ONLY)
CONFIGURE_FILE(include/heplib++ "${CMAKE_CURRENT_BINARY_DIR}/OTHER/heplib++" @ONLY)
CONFIGURE_FILE(include/update "${CMAKE_CURRENT_BINARY_DIR}/OTHER/update" @ONLY)
#-----------------------------------------------
# External Libraries
#-----------------------------------------------
LINK_LIBRARIES(ginac cln Minuit2 mpfr qhullstatic quadmath gomp cubaq)
#-----------------------------------------------
PROJECT(HepLib VERSION 1.0 LANGUAGES C CXX Fortran)
AUX_SOURCE_DIRECTORY(BASIC BASIC_SRCS)
AUX_SOURCE_DIRECTORY(SD SD_SRCS)
AUX_SOURCE_DIRECTORY(HEP HEP_SRCS)
AUX_SOURCE_DIRECTORY(IBP IBP_SRCS)
AUX_SOURCE_DIRECTORY(EX EX_SRCS)
AUX_SOURCE_DIRECTORY("${CMAKE_CURRENT_BINARY_DIR}/OTHER" OTHER)
#-----------------------------------------------
# WSTP - Wolfram Symbolic Transfer Protocol
#-----------------------------------------------
IF(APPLE)
FIND_LIBRARY(HAS_WSTP WSTPi4 ${LINK_DIRECTORIES})
IF(HAS_WSTP)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lc++ -lWSTPi4 -framework Foundation")
ELSE()
STRING(REPLACE "BASIC/WSKernel.cpp" "" BASIC_SRCS "${BASIC_SRCS}")
ENDIF(HAS_WSTP)
ELSE()
FIND_LIBRARY(HAS_WSTP WSTP64i4 ${LINK_DIRECTORIES})
IF(HAS_WSTP)
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -lm -lpthread -lrt -lstdc++ -ldl -luuid -lWSTP64i4")
ELSE()
STRING(REPLACE "BASIC/WSKernel.cpp" "" BASIC_SRCS "${BASIC_SRCS}")
ENDIF(HAS_WSTP)
ENDIF(APPLE)
#-----------------------------------------------
# install library HepLib
#-----------------------------------------------
ADD_LIBRARY(HepLib SHARED ${BASIC_SRCS} ${SD_SRCS} ${HEP_SRCS} ${IBP_SRCS} ${EX_SRCS} ${OTHER})
SET_TARGET_PROPERTIES(HepLib PROPERTIES INSTALL_RPATH "${LIB_PATH}")
INSTALL(TARGETS HepLib DESTINATION lib)
#-----------------------------------------------
# install include/*.h & other *.h
#-----------------------------------------------
INSTALL(DIRECTORY ${CMAKE_SOURCE_DIR}/include
DESTINATION ${CMAKE_INSTALL_PREFIX}
FILES_MATCHING PATTERN "*.h")
INSTALL(FILES "${CMAKE_SOURCE_DIR}/EX/Wrap.h" DESTINATION include)
#-----------------------------------------------
# install bin/*
#-----------------------------------------------
FILE(GLOB cpp_files "bin/*.cpp")
FOREACH(exe_cpp ${cpp_files})
GET_FILENAME_COMPONENT(exe ${exe_cpp} NAME_WE)
ADD_EXECUTABLE(${exe} bin/${exe}.cpp)
SET_TARGET_PROPERTIES(${exe} PROPERTIES INSTALL_RPATH "${LIB_PATH}")
TARGET_LINK_LIBRARIES(${exe} HepLib)
INSTALL(TARGETS ${exe} DESTINATION bin)
ENDFOREACH()
INSTALL(TARGETS qgraf DESTINATION bin)
INSTALL(FILES "${CMAKE_CURRENT_BINARY_DIR}/OTHER/HepLib.pc" DESTINATION "lib/pkgconfig")
INSTALL(FILES "${CMAKE_CURRENT_BINARY_DIR}/OTHER/heplib++" DESTINATION "bin"
PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ)
INSTALL(FILES "${CMAKE_CURRENT_BINARY_DIR}/OTHER/update" DESTINATION "bin"
PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ)
# OpenMP
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp -std=gnu++14")
BAK/Laporta.cpp deleted 100644 → 0
View file @ 8abf1d7e
#include "IBP.h"
namespace HepLib::IBP {
lst Laporta::F2lst(ex f) {
if(!f.match(F(w))) {
cerr << ErrColor << "Not F(w) form: " << f << RESET << endl;
exit(1);
}
if(ccF[f].nops()>0) return ccF[f];
auto al = f.op(0);
lst sl, rl;
ex r=0, s=0, sid=0;
for(int i=0; i<al.nops(); i++) {
if(al.op(i) > ex(1)/2) {
sid += pow(2, i);
r += al.op(i);
rl.append(al.op(i));
} else {
s -= al.op(i);
sl.append(al.op(i));
}
}
lst ret = lst{sid, r+s, r, s, sl, rl};
if(Cuts.nops()>0) {
lst clst;
for(auto c : Cuts) {
int ci = ex_to<numeric>(c).to_int();
clst.append(al.op(ci));
}
ret.prepend(clst);
}
ccF[f] = ret;
ccFinv[ret] = f;
return ret;
}
lst Laporta::I2lst(ex ibp) {
if(ccI[ibp].nops()>0) return ccI[ibp];
lst fis;
exset fiset;
ibp.find(F(w), fiset);
for(auto &item : fiset) {
fis.append(F2lst(item));
}
for(int i=0; i<fis.nops(); i++) {
for(int j=i+1; j<fis.nops(); j++) {
if(ex_less(fis.op(i), fis.op(j))) {
auto tmp = fis.op(i);
fis.let_op(i) = fis.op(j);
fis.let_op(j) = tmp;
}
}
}
ccImax[ibp] = ex_to<lst>(fis.op(0));
fis.prepend(0);
fis.let_op(0) = fis.op(1);
fis.let_op(1) = fis.nops()-1;
ccI[ibp] = fis;
return fis;
}
ex Laporta::collectF(ex expr) {
if(!expr.subs(F(w)==0).is_zero()) {
cerr << ErrColor << "expr is NOT zero with F(w)==0: " << expr << RESET << endl;
exit(1);
}
exset fs;
expr.find(F(w), fs);
ex ret = 0;
for(auto item : fs) {
ret += normal(expr.coeff(item, 1)) * item;
}
return ret;
}
void Laporta::Prepare(lst loop, lst ext, lst prop, lst repl) {
if(Verbose > 0) cout << " IBP Prepare @ " << now() << endl;
lst sps;
for(auto it : loop) {
for(auto ii : loop) sps.append(it*ii);
for(auto ii : ext) sps.append(it*ii);
}
sps.sort();
sps.unique();
lst sp2s, s2sp, ss;
for(auto item : sps) {
symbol si;
ss.append(si);
sp2s.append(item==si);
s2sp.append(si==item);
}
lst eqns;
for(int i=0; i<prop.nops(); i++) {
auto eq = prop.op(i).expand();
eq = eq.subs(sp2s, subs_options::algebraic);
eq = eq.subs(repl, subs_options::algebraic);
eqns.append(eq == iWF(i));
}
auto s2p = lsolve(eqns, ss);
preIBPs.remove_all();
lst ns;
for(int i=0; i<prop.nops(); i++) ns.append(a(i));
for(auto l : loop) {
ex ibp = 0;
symbol sl;
for(int i=0; i<prop.nops(); i++) {
auto ns_tmp = ns;
ns_tmp.let_op(i) = ns.op(i) + 1;
auto dp = prop.op(i).subs(l==sl).diff(sl).subs(sl==l);
ibp -= a(i) * F(ns_tmp) * dp;
}
for(auto ii : ext) {
auto ibp_tmp = ibp * ii;
ibp_tmp = ibp_tmp.expand();
ibp_tmp = ibp_tmp.subs(sp2s, subs_options::algebraic);
ibp_tmp = ibp_tmp.subs(repl, subs_options::algebraic);
ibp_tmp = ibp_tmp.subs(s2p, subs_options::algebraic);
ex res = 0;
for(int i=0; i<prop.nops(); i++) {
auto ci = ibp_tmp.coeff(iWF(i), 1);
ci = MapFunction([i](const ex & e, MapFunction &self)->ex{
if(e.match(F(w))) {
auto tmp = e.op(0);
tmp.let_op(i) = tmp.op(i)-1;
return F(tmp);
} else if(!e.has(F(w))) return e;
else return e.map(self);
})(ci);
res += ci;
}
res += ibp_tmp.subs(lst{iWF(w)==0});
preIBPs.append(res);
}
for(auto ii : loop) {
auto ibp_tmp = ibp * ii;
ibp_tmp = ibp_tmp.expand();
ibp_tmp = ibp_tmp.subs(sp2s, subs_options::algebraic);
ibp_tmp = ibp_tmp.subs(repl, subs_options::algebraic);
ibp_tmp = ibp_tmp.subs(s2p, subs_options::algebraic);
ex res = 0;
for(int i=0; i<prop.nops(); i++) {
auto ci = ibp_tmp.coeff(iWF(i), 1);
ci = MapFunction([i](const ex &e, MapFunction &self)->ex {
if(e.match(F(w))) {
auto tmp = e.op(0);
tmp.let_op(i) = tmp.op(i)-1;
return F(tmp);
} else if(!e.has(F(w))) return e;
else return e.map(self);
})(ci);
res += ci;
}
res += ibp_tmp.subs(lst{iWF(w)==0});
if(ii==l) res += d*F(ns);
preIBPs.append(res);
}
}
}
void Laporta::Generate2(lst seed) {
lst ns, a2s;
for(int i=0; i<seed.nops(); i++) {
Symbol s("a"+to_string(i));
a2s.append(a(i)==s+Symbol("eps"));
ns.append(s);
}
lst seeds;
for(auto const & item : preIBPs) {
exset fs;
item.find(F(w), fs);
for(auto fi : fs) {
lst eqs;
for(int i=0; i<fi.op(0).nops(); i++) {
eqs.append(subs(fi.op(0).op(i)==seed.op(i),a2s));
}
auto sol = lsolve(eqs, ns);
seeds.append(subs(ns,sol));
auto ii = item.subs(a2s).subs(sol);
if(ii.is_zero()) continue;
IBPs.append(ii);
}
}
for(auto seed : seeds)
for(auto const & item : preIBPs) {
exset fs;
item.find(F(w), fs);
for(auto fi : fs) {
lst eqs;
for(int i=0; i<fi.op(0).nops(); i++) {
eqs.append(subs(fi.op(0).op(i)==seed.op(i),a2s));
}
auto sol = lsolve(eqs, ns);
auto ii = item.subs(a2s).subs(sol);
if(ii.is_zero()) continue;
IBPs.append(ii);
}
}
if(Cuts.nops()>1) {
int total = IBPs.nops();
for(int i=0; i<IBPs.nops(); i++) {
auto tmp = IBPs.op(i);
exset fs;
find(tmp, F(w), fs);
exmap repl;
for(auto fi : fs){
for(auto ic : Cuts) {
int j = ex_to<numeric>(ic).to_int();
if(fi.op(0).op(j)<=0) {
repl[fi]=0;
break;
}
}
}
IBPs.let_op(i) = tmp.subs(repl);
}
}
}
void Laporta::Generate(vector<lst> seeds) {
if(Verbose > 0) cout << " IBP Generate @ " << now() << endl;
IBPs.remove_all();
lst ns;
for(int i=0; i<seeds[0].nops(); i++) ns.append(a(i));
for(auto seed : seeds) {
for(auto const & item : preIBPs) {
auto ii = item.subs(ns, seed);
if(ii.is_zero()) continue;
IBPs.append(ii);
}}
if(Cuts.nops()>1) {
int total = IBPs.nops();
for(int i=0; i<IBPs.nops(); i++) {
auto tmp = IBPs.op(i);
exset fs;
find(tmp, F(w), fs);
exmap repl;
for(auto fi : fs){
for(auto ic : Cuts) {
int j = ex_to<numeric>(ic).to_int();
if(fi.op(0).op(j)<=0) {
repl[fi]=0;
break;
}
}
}
IBPs.let_op(i) = tmp.subs(repl);
}
}
if(Verbose > 0) cout << " - Number of Identities: " << IBPs.nops() << endl << flush;
}
void Laporta::Reduce() {
if(Verbose > 0) cout << now() << " - IBP Reduce ..." << endl << flush;
if(Verbose > 1) cout << " - Sorting ... @ " << now(false) << endl << flush;
vector<ex> ibps;
for(auto item : IBPs) {
if(item.is_zero()) continue;
I2lst(item);
//if(F2lst(F(lst{1,1,1,0,1,1,1})).op(2)+2<ccImax[item].op(2)) continue;
ibps.push_back(item);
}
exset fs;
find(IBPs, F(w), fs);
vector<ex> fvec;
for(auto item : fs) fvec.push_back(item);
// sort Fs
sort(fvec.begin(), fvec.end(), [this](auto &a, auto &b) {
return ex_less(F2lst(b),F2lst(a));
}); // larger first
// sort Identities
sort(ibps.begin(), ibps.end(), [this](auto &a, auto &b) {
lst ai = I2lst(a);
lst bi = I2lst(b);
return ex_less(bi,ai);
}); // larger first
map<ex, int, ex_is_less> f2idx;
for(int i=0; i<fvec.size(); i++) f2idx[fvec[i]] = i;
int cols = fvec.size();
int rows = ibps.size();
if(Verbose > 1) cout << " - Fermating @ " << now(false) << " :> " << flush;
Fermat fermat;
fermat.Init();
auto Variables = gather_symbols(ibps);
ostringstream ss;
for(auto j : Variables) ss << "&(J="<<j<<");" << endl;
fermat.Execute(ss.str());
ss.clear();
ss.str("");
ss << "Array mat[" << rows << "," << cols+1 << "] Sparse;" << endl;
for(int r=0; r<rows; r++) {
auto ibp = ibps[r];
exset fs;
find(ibp, F(w), fs);
for(auto f : fs) {
ss << "mat["<<(r+1)<<","<<(f2idx[f]+1)<<"] := " << ibp.coeff(f, 1) << ";" << endl;
fermat.Execute(ss.str());
ss.clear();
ss.str("");
}
}
ss << "Redrowech([mat]);" << endl;
fermat.Execute(ss.str());
ss.clear();
ss.str("");
ss << "&(U=1);" << endl; // ugly printing, the whitespace matters
ss << "![mat" << endl;
auto ostr = fermat.Execute(ss.str());
fermat.Exit();
if(Verbose > 1) cout << now(false) << endl;
// make sure last char is 0
if(ostr[ostr.length()-1]!='0') throw Error("Reduce: 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;
auto mat2 = fp.Read(ostr);
ex irows[rows];
for(int i=0; i<rows; i++) {
irows[i] = 0;
for(int j=0; j<fvec.size(); j++) {
if(is_zero(get_op(mat2,i,j))) continue;
irows[i] += get_op(mat2,i,j) * fvec[j];
}
}
if(Verbose > 1) cout << " - Finalizing ... @ " << now(false) << endl << flush;
ibps.clear();
for(int i=0; i<rows; i++) {
if(irows[i].is_zero()) continue;
ibps.push_back(irows[i]);
}
Rules.remove_all();
for(int i=0; i<ibps.size(); i++) {
auto ibp = ibps[ibps.size()-i-1];
auto ifmt = I2lst(ibp);
auto fx = ccFinv[ex_to<lst>(ifmt.op(0))];
auto c1 = ibp.coeff(fx, 1);
auto c0 = ibp.coeff(fx, 0);
auto res = -c0/c1;
if(!(c1-1).is_zero()) {
cout << "c1 = " << c1 << endl;
cout << "ibps.size()-i-1 : " << ibps.size()-i-1 << endl;
cout << ibp << endl;
cout << endl << endl;
}
Rules.append(fx==res);
}
}
}
BAK/WSKernel.cpp deleted 100644 → 0
View file @ 8abf1d7e
#include "WSKernel.h"
namespace HepLib {
const char * WSKernel::Error::what() const throw () {
return msg.c_str();
}
WSKernel::Error::Error(WSLINK lp) {
const char *err_msg = WSErrorMessage(lp);
msg = err_msg;
if(err_msg!=NULL) WSReleaseErrorMessage(lp, err_msg);
};
WSKernel::WSKernel(const string & open_str) {
int error;
ep = WSInitialize((WSParametersPointer)0);
if(ep == (WSENV)0) throw Error(lp);
string ostr = open_str;
if(ostr.length()<1) {
#ifdef __APPLE__
ostr = "-linkmode launch -linkname '/Applications/Mathematica.app/Contents/MacOS/WolframKernel -wstp'";
#else
ostr = "-linkmode launch -linkname 'math -wstp'";
#endif
}
lp = WSOpenString(ep, ostr.c_str(), &error);
//lp = WSOpenString(ep, "-linkcreate -linkprotocol IntraProcess", &error);
//lp = WSOpenArgcArgv(ep, argc, argv, &error);
if(lp == (WSLINK)0 || error != WSEOK) throw Error(lp);
}
WSKernel::~WSKernel() {
WSClose(lp);
WSDeinitialize(ep);
}
string WSKernel::Evaluate(const string & expr) {
WSPutFunction(lp, "EvaluatePacket", 1);
WSPutFunction(lp, "ToString", 1);
WSPutFunction(lp, "InputForm", 1);
WSPutFunction(lp, "ToExpression", 1);
WSPutString(lp, expr.c_str());
WSEndPacket(lp);
WSFlush(lp);
int pkt;
while( (pkt=WSNextPacket(lp)) && (pkt!= RETURNPKT) ) WSNewPacket(lp);
const char *mma_out;
if(!WSGetString(lp, &mma_out)) throw Error(lp);
char oo[1024];
strcpy(oo, mma_out);
string ret_string = oo;
if(!WSEndPacket(lp)) throw Error(lp);
WSReleaseString(lp, mma_out);
WSFlush(lp);
return ret_string;
}
}
BAK/WSKernel.h deleted 100644 → 0
View file @ 8abf1d7e
/**
* @file
* @brief WSKernel header file
*/
#pragma once
#include "wstp.h"
#include <iostream>
#include <cstring>
namespace HepLib {
using namespace std;
/**
* @brief interface to Wolfram Mathematica
*/
class WSKernel {
public:
/**
* @brief inner Error class for WSKernel
*/
class Error : public exception {
public:
string msg;
const char * what() const throw ();
Error(WSLINK lp);
};
WSKernel(const string & open_str="");
~WSKernel();
string Evaluate(const string & expr);
private:
WSENV ep;
WSLINK lp;
};
}
BAK/wstp.h deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
BASIC/BASIC.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
BASIC/BASIC.h deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
BASIC/Functions.cpp deleted 100644 → 0
View file @ 8abf1d7e
/**
* @file
* @brief Basic/Useful Functions
*/
#include "BASIC.h"
namespace HepLib {
REGISTER_FUNCTION(coCF, do_not_evalf_params())
REGISTER_FUNCTION(coVF, do_not_evalf_params())
REGISTER_FUNCTION(x, do_not_evalf_params())
REGISTER_FUNCTION(y, do_not_evalf_params())
REGISTER_FUNCTION(z, do_not_evalf_params())
// F function, up to 2 arguments
unsigned F1_SERIAL::serial = GiNaC::function::register_new(function_options("F",1).do_not_evalf_params().overloaded(2));
unsigned F2_SERIAL::serial = GiNaC::function::register_new(function_options("F",2).do_not_evalf_params().overloaded(2));
// WF function, up to 5 arguments
unsigned WF1_SERIAL::serial = GiNaC::function::register_new(function_options("WF",1).do_not_evalf_params().overloaded(5));
unsigned WF2_SERIAL::serial = GiNaC::function::register_new(function_options("WF",2).do_not_evalf_params().overloaded(5));
unsigned WF3_SERIAL::serial = GiNaC::function::register_new(function_options("WF",3).do_not_evalf_params().overloaded(5));
unsigned WF4_SERIAL::serial = GiNaC::function::register_new(function_options("WF",4).do_not_evalf_params().overloaded(5));
unsigned WF5_SERIAL::serial = GiNaC::function::register_new(function_options("WF",5).do_not_evalf_params().overloaded(5));
// iWF function, up to 5 arguments
unsigned iWF1_SERIAL::serial = GiNaC::function::register_new(function_options("iWF",1).do_not_evalf_params().overloaded(5));
unsigned iWF2_SERIAL::serial = GiNaC::function::register_new(function_options("iWF",2).do_not_evalf_params().overloaded(5));
unsigned iWF3_SERIAL::serial = GiNaC::function::register_new(function_options("iWF",3).do_not_evalf_params().overloaded(5));
unsigned iWF4_SERIAL::serial = GiNaC::function::register_new(function_options("iWF",4).do_not_evalf_params().overloaded(5));
unsigned iWF5_SERIAL::serial = GiNaC::function::register_new(function_options("iWF",5).do_not_evalf_params().overloaded(5));
/*-----------------------------------------------------*/
// MapFunction Class
/*-----------------------------------------------------*/
MapFunction::MapFunction(std::function<ex(const ex &, MapFunction &)> func) : Function(func) { }
ex MapFunction::operator()(const ex &e) {
return Function(e, *this);
}
ex MapFunction::subs(const ex & expr, const ex & pat, std::function<ex(const ex &)>f) {
MapFunction map([pat,f](const ex & e, MapFunction &self)->ex{
if(!e.has(pat)) return e;
else if(e.match(pat)) return f(e);
else return e.map(self);
});
return map(expr);
}
/*-----------------------------------------------------*/
// Parser Class
/*-----------------------------------------------------*/
// copy from GiNaC parser, note the alignas(2)
namespace {
alignas(2) static ex sqrt_reader(const exvector& ev) {
return GiNaC::sqrt(ev[0]);
}
alignas(2) static ex pow_reader(const exvector& ev) {
return GiNaC::pow(ev[0], ev[1]);
}
alignas(2) static ex power_reader(const exvector& ev) {
return GiNaC::power(ev[0], ev[1]);
}
alignas(2) static ex lst_reader(const exvector& ev) {
return GiNaC::lst(ev.begin(), ev.end());
}
class functions_hack : public GiNaC::function {
public:
static const std::vector<function_options>& get_registered_functions() {
return function::registered_functions();
}
};
static reader_func encode_serial_as_reader_func(unsigned serial) {
uintptr_t u = (uintptr_t)serial;
u = (u << 1) | (uintptr_t)1;
reader_func ptr = (reader_func)((void *)u);
return ptr;
}
const prototype_table& ginac_reader() {
using std::make_pair;
static bool initialized = false;
static prototype_table reader;
if (!initialized) {
reader[make_pair("sqrt", 1)] = sqrt_reader;
reader[make_pair("pow", 2)] = pow_reader;
reader[make_pair("power", 2)] = power_reader;
reader[make_pair("lst", 0)] = lst_reader;
enum { log, exp, sin, cos, tan, asin, acos, atan, sinh, cosh, tanh, asinh, acosh, atanh, atan2,
Li2, Li3, zetaderiv, Li, S, H, lgamma, tgamma, beta, factorial, binomial, Order,
NFUNCTIONS
};
auto it = functions_hack::get_registered_functions().begin();
unsigned serial = 0;
for ( ; serial<NFUNCTIONS; ++it, ++serial ) {
prototype proto = make_pair(it->get_name(), it->get_nparams());
reader[proto] = encode_serial_as_reader_func(serial);
}
initialized = true;
}
return reader;
}
}
// note iSymbol will override Symbol if they share same name
Parser::Parser() : STable(Symbol::Table), FTable(ginac_reader()) {
for(auto item : iSymbol::Table) STable[item.first] = item.second;
}
Parser::Parser(symtab st) : STable(Symbol::Table), FTable(ginac_reader()) {
for(auto item : iSymbol::Table) STable[item.first] = item.second;
for(auto item : st) STable[item.first] = item.second;
}
ex Parser::Read(string instr, bool s2S) {
unsigned serial = 0;
for (auto & it : functions_hack::get_registered_functions()) {
prototype proto = make_pair(it.get_name(), it.get_nparams());
FTable[proto] = encode_serial_as_reader_func(serial);
++serial;
}
parser ginac_parser(STable, false, FTable);
ex ret = ginac_parser(instr);
if(!s2S) return ret;
// check & replace symbol to Symbol object
auto st = ginac_parser.get_syms();
bool redo = false;
exmap repl;
for(auto kv : st) {
if(is_exactly_a<symbol>(kv.second) && STable.find(kv.first)==STable.end()) {
string ss = kv.first;
STable[ss] = Symbol(ss);
redo = true;
repl[kv.second] = Symbol(ss);
}
}
if(redo) return ret.subs(repl);
else return ret;
}
ex Parser::ReadFile(string filename, bool s2S) {
ifstream ifs(filename);
string ostr((istreambuf_iterator<char>(ifs)), (istreambuf_iterator<char>()));
ifs.close();
return Read(ostr,s2S);
}
/*-----------------------------------------------------*/
// string Functions
/*-----------------------------------------------------*/
void string_replace_all(string &str, const string &from, const string &to) {
size_t start_pos = 0;
while((start_pos = str.find(from, start_pos)) != string::npos) {
str.replace(start_pos, from.length(), to);
start_pos += to.length();
}
}
void string_trim(string &ostr) {
const char* WhiteSpace = " \t\v\r\n";
if(!ostr.empty()) {
ostr.erase(0, ostr.find_first_not_of(WhiteSpace));
ostr.erase(ostr.find_last_not_of(WhiteSpace)+1);
}
}
void Combinations(int n, int m, std::function<void(const int*)> f) {
if(m<1 || m>n) return;
std::string bitmask(m,1);
bitmask.resize(n,0);
do {
int is[m]; int j=0;
for (int i=0; i<n; ++i) {
if(bitmask[i]) { is[j]=i; j++; }
}
f(is);
} while (std::prev_permutation(bitmask.begin(), bitmask.end()));
}
void CombinationsR(int n, int m, std::function<void(const int*)> f) {
if(m<1) return;
Combinations(n+m-1, n-1, [n,m,f](const int* isr) {
int is[m];
int mi=m;
for(int ni=0; ni<n; ni++) {
int c=0;
if(ni==n-1) c=(n+m-2)-isr[n-2];
else if(ni==0) c=isr[ni];
else c=isr[ni]-isr[ni-1]-1;
for(int j=0; j<c; j++) is[--mi] = n-1-ni;
}
f(is);
});
}
void Permutations(int n, std::function<void(const int*)> f) {
if(n<1) return;
int pis[n];
for(int i=0; i<n; i++) pis[i]=i;
do { f(pis); } while(std::next_permutation(pis,pis+n));
}
void Permutations(int n, int m, std::function<void(const int*)> f) {
if(m<1 || m>n) return;
Combinations(n, m, [m,f](const int *ns1) {
Permutations(m, [m,f,ns1](const int *ns2) {
int ns[m];
for(int i=0; i<m; i++) ns[i] = ns1[ns2[i]];
f(ns);
});
});
}
namespace {
struct Generator {
public:
int *a;
Generator(int s, int v) : cSlots(s) , cValues(v) {
a = new int[s];
for(int i = 0; i<cSlots-1; i++) a[i]=1-1;
a[cSlots-1]=0-1;
nextInd = cSlots;
}
~Generator() { delete a; }
bool doNext() {
for (;;) {
if (a[nextInd-1]==cValues-1) {
nextInd--;
if(nextInd==0) return false;
} else {
a[nextInd-1]++;
while (nextInd<cSlots) {
nextInd++;
a[nextInd-1]=1-1;
}
return true;
}
}
}
private:
int cSlots;
int cValues;
int nextInd;
};
}
void PermutationsR(int n, int m, std::function<void(const int*)> f) {
Generator g(m,n);
while (g.doNext()) f(g.a);
}
bool isSorted(const lst & exs) {
for(int i=0; i<exs.nops()-1; i++) {
if(exs.op(i).is_equal(exs.op(i+1))) continue;
if(!ex_less(exs.op(i),exs.op(i+1))) return false;
}
return true;
}
bool isSorted(int n, const ex exs[]) {
for(int i=0; i<n-1; i++) {
if(exs[i].is_equal(exs[i+1])) continue;
if(!ex_less(exs[i],exs[i+1])) return false;
}
return true;
}
int ACSort(lst & exs) {
int ac = 0;
int n = exs.nops();
for(int i=0; i<n-1; i++)
for(int j=n-1; j>i; j--)
if(ex_less(exs.op(j),exs.op(j-1))) {
auto tmp = exs.op(j-1);
exs.let_op(j-1) = exs.op(j);
exs.let_op(j) = tmp;
ac++;
}
for(int i=0; i<n-1; i++) if(exs.op(i).is_equal(exs.op(i+1))) return 0;
return (ac%2==1) ? -1 : 1;
}
int ACSort(int n, ex exs[]) {
int ac = 0;
for(int i=0; i<n-1; i++)
for(int j=n-1; j>i; j--)
if(ex_less(exs[j],exs[j-1])) {
auto tmp = exs[j-1];
exs[j-1] = exs[j];
exs[j] = tmp;
ac++;
}
for(int i=0; i<n-1; i++) if(exs[i].is_equal(exs[i+1])) return 0;
return (ac%2==1) ? -1 : 1;
}
}
BASIC/Process.cpp deleted 100644 → 0
View file @ 8abf1d7e
/**
* @file
* @brief Class/Function for Process
*/
#include "BASIC.h"
#include <fstream>
inline bool file_exists(const char* fn) {
return (access(fn,F_OK)!=-1);
}
namespace HepLib {
//-----------------------------------------------------------
// Fermat Class
//-----------------------------------------------------------
#define ENTER endl<<endl<<endl
Fermat::~Fermat() { Exit(); }
void Fermat::Exit() {
if(getpid() != pid) return; // happens @ fork child process
if(inited) {
ostringstream script;
script << "&q;" << endl << "&x;" << ENTER;
string istr = script.str();
write(P2C[1], istr.c_str(), istr.length());
int st;
waitpid(fpid, &st, WUNTRACED);
inited = false;
exited = true;
}
}
void Fermat::Init(string fer_path) {
if(inited) return;
inited = true;
pid = getpid();
if (pipe(P2C)==-1 || pipe(C2P)==-1) {
throw Error("pipe failed in Fermat::Init.");
}
fpid = fork();
if (fpid == 0) { // child process
setpgid(0,0);
close(P2C[1]);
close(C2P[0]);
dup2(C2P[1], 1);
close(C2P[1]);
dup2(P2C[0], 0);
close(P2C[0]);
execlp(fer_path.c_str(), fer_path.c_str(), NULL);
exit(0);
}
// parent process
close(P2C[0]); // P2C[1] for write
close(C2P[1]); // C2P[0] for read
ostringstream script;
script << "&(_d=90000);" << endl; // width of the display on the window
script << "&d" << endl << "0;" << endl; // off floating point representation
script << "&(_t=0);" << endl; // off a certain fast probabalistic algorithm
script << "&(t=0);" << endl; // off timing
script << "&(_s=0);" << endl;
script << "&(_o=1000);" << endl; // http://home.bway.net/lewis/fer64mono.html
script << "&(M=' ');" << endl; // prompt
script << "!('" << Sentinel << "');" << ENTER;
string istr = script.str();
write(P2C[1], istr.c_str(), istr.length());
string ostr;
int n = 1024;
char buffer[n+1]; // make sure the last one is '\0'
int nio;
while(true) {
for(int i=0; i<n+1; i++) buffer[i] = '\0';
nio = read(C2P[0], buffer, n);
if(nio>0) ostr += buffer;
else throw Error(ostr);
auto cpos = ostr.find(Sentinel);
if(cpos!=string::npos) {
const char* WhiteSpace = " \t\v\r\n";
auto lpos = ostr.find_last_not_of(WhiteSpace);
if(ostr[lpos]!='0') read(C2P[0], buffer, n); // last 0, due to Sentinel
ostr.erase(cpos);
break;
}
}
string_replace_all(ostr, "*** entry > 30 or < 5 means turn off mono multiply.", "");
string_replace_all(ostr, "*** Fermat Warning. Early exit from mod_multivar_Chinese.", "");
if(ostr.find("***")!=string::npos) {
cout << "Fermat script: " << endl << istr << endl << endl;
throw Error(ostr.c_str());
}
}
// out string still contains the last number
string Fermat::Execute(string expr) {
if(exited) throw Error("Fermat has already exited.");
if(getpid() != pid) throw Error("Fermat: can not Execute on child process.");
ostringstream script;
script << expr << endl;
script << "!('" << Sentinel << "')" << ENTER;
string istr = script.str();
write(P2C[1], istr.c_str(), istr.length());
string ostr;
int n = buffer_size;
char buffer[n+1]; // make sure the last one is '\0'
int nio;
while(true) {
for(int i=0; i<n+1; i++) buffer[i] = '\0';
nio = read(C2P[0], buffer, n);
if(nio>0) ostr += buffer;
else throw Error(ostr);
auto cpos = ostr.find(Sentinel);
if(cpos!=string::npos) {
const char* WhiteSpace = " \t\v\r\n";
auto lpos = ostr.find_last_not_of(WhiteSpace);
if(ostr[lpos]!='0') read(C2P[0], buffer, n); // last 0, due to Sentinel
ostr.erase(cpos);
break;
}
}
string_trim(ostr);
if(ostr.find("***")!=string::npos) throw Error(ostr.c_str());
return ostr;
}
//-----------------------------------------------------------
// Form Class
//-----------------------------------------------------------
Form::~Form() { Exit(); }
void Form::Exit() {
if(getpid() != pid) return; // happens @ fork child process
if(inited && !exited) {
string exit_cmd = "\n.end\n" + Prompt +"\n";
write(io[0][1], exit_cmd.c_str(), exit_cmd.length());
char buffer[8];
read(io[1][0], buffer, 8);
int st;
waitpid(fpid, &st, WUNTRACED);
}
inited = false;
exited = true;
}
void Form::Init(string form_path) {
if(inited) return;
inited = true;
pid = getpid();
if (pipe(io[0])==-1 || pipe(io[1])==-1 || pipe(stdo)==-1) {
inited = false;
exited = true;
throw Error("pipe failed in Form::Init.");
}
fpid = fork();
if (fpid == 0) {
setpgid(0,0);
close(io[0][1]);
close(io[1][0]);
close(stdo[0]);
dup2(stdo[1], 1);
auto cpid = getpid(); // current process id
ostringstream oss;
oss << "init-" << cpid << ".frm";
std::ofstream ofs;
ofs.open(oss.str().c_str(), ios::out);
if (!ofs) {
inited = false;
exited = true;
throw Error("failed to open init.frm file!");
}
ofs << "Off Statistics;" << endl;
ofs << "#ifndef `PIPES_'" << endl;
ofs << " #message \"No pipes found\";" << endl;
ofs << " .end;" << endl;
ofs << "#endif" << endl;
ofs << "#if (`PIPES_' <= 0)" << endl;
ofs << " #message \"No pipes found\";" << endl;
ofs << " .end;" << endl;
ofs << "#endif" << endl;
ofs << "#procedure put(mexp)" << endl;
ofs << " #write \""<<Sentinel<<"\\n\"" << endl;
ofs << " #toexternal \"%E\", `mexp'" << endl;
ofs << " #toexternal \""<<Sentinel<<"\"" << endl;
ofs << " #toexternal \"\\n\"" << endl;
ofs << "#endprocedure" << endl;
ofs << "#setexternal `PIPE1_';" << endl;
ofs << "#toexternal \"OK\"" << endl;
ofs << "Local [o]=0;" << endl;
ofs << ".sort" << endl;
ofs << "Format Mathematica;" << endl;
ofs << "#prompt \"" << Prompt << "\"" << endl;
ofs << "#fromexternal-" << endl;
ofs << ".end" << endl;
ofs.close();
execlp(form_path.c_str(), form_path.c_str(),
"-pipe",
(to_string(io[0][0])+","+to_string(io[1][1])).c_str(),
"-M", ("init-"+to_string(cpid)).c_str(),
NULL);
exit(0);
}
close(io[0][0]);
close(io[1][1]);
close(stdo[1]);
char buffer[1024];
read(io[1][0], buffer, sizeof(buffer));
char* p = strstr(buffer, "\n");
if(p==NULL){
inited = false;
exited = true;
cout << "the return is: <|" << buffer << "|>" << endl;
throw Error("Init Failed: Expect a Line break!");
}
sprintf(p, ",%d\n\n", fpid);
write(io[0][1], buffer, strlen(buffer));
read(io[1][0], buffer, sizeof(buffer));
p = strstr(buffer, "OK");
if(p==NULL || p!=buffer) {
inited = false;
exited = true;
throw Error("Init Failed: Expect OK!");
}
ostringstream oss;
oss << "init-" << fpid << ".frm";
if(file_exists(oss.str().c_str())) remove(oss.str().c_str());
if(true) { // read the terminal output
string estr;
int n = 1024;
char buffer[n+1]; // make sure the last one is '\0'
int nio;
while(true) {
for(int i=0; i<n+1; i++) buffer[i] = '\0';
nio = read(stdo[0], buffer, n);
if(nio>0) estr += buffer;
auto cpos = estr.find(Sentinel);
if(cpos!=string::npos) break;
cpos = estr.find("-->");
if(cpos!=string::npos) {
inited = false;
exited = true;
throw Error(estr);
}
}
}
}
string Form::Execute(string script, const string & out_var) {
if(exited) throw Error("Form has already exited.");
if(getpid() != pid) throw Error("Form: can not Execute on child process.");
string istr = script;
istr += "\n.sort\n#call put(";
istr += out_var;
istr += ")\n.sort\n";
istr += Prompt + "\n"; // prompt
write(io[0][1], istr.c_str(), istr.length());
string ostr;
int n = buffer_size;
char buffer[n+1]; // make sure the last one is '\0'
int nio;
string estr;
while(true) {
for(int i=0; i<n+1; i++) buffer[i] = '\0';
nio = read(stdo[0], buffer, n);
if(nio>0) estr += buffer;
auto cpos = estr.find(Sentinel);
if(cpos!=string::npos) break;
cpos = estr.find("-->");
if(cpos!=string::npos) {
inited = false;
exited = true;
throw Error(estr);
}
}
while(true) {
for(int i=0; i<n+1; i++) buffer[i] = '\0';
nio = read(io[1][0], buffer, n);
if(nio>0) ostr += buffer;
else {
string estr;
while(true) {
for(int i=0; i<n+1; i++) buffer[i] = '\0';
nio = read(stdo[0], buffer, n);
if(nio<=0) break;
estr += buffer;
}
inited = false;
exited = true;
throw Error(estr.c_str());
}
auto cpos = ostr.find(Sentinel);
if(cpos!=string::npos) {
ostr.replace(cpos, Sentinel.length(), "");
break;
}
}
return ostr;
}
}
CMakeLists.txt deleted 100644 → 0
View file @ 8abf1d7e
# usage
# cmake -DCMAKE_INSTALL_PREFIX=p0 -DINC_PATH="p1;p2" -DLIB_PATH="p3;p4" path
CMAKE_MINIMUM_REQUIRED(VERSION 3.10)
SET(CMAKE_C_COMPILER gcc)
SET(CMAKE_CXX_COMPILER g++)
SET(CMAKE_FORTRAN_COMPILER gfortran)
#-----------------------------------------------
# qgraf
#-----------------------------------------------
ADD_EXECUTABLE(qgraf "Qgraf/qgraf-3.1.4.f")
#-----------------------------------------------
# INC_PATH & LIB_PATH
#-----------------------------------------------
SET(INC_PATH "${INC_PATH};${CMAKE_INSTALL_PREFIX}/include")
LIST(REMOVE_DUPLICATES INC_PATH)
SET(INC_FLAGS "")
FOREACH(dir ${INC_PATH})
IF(EXISTS ${dir})
INCLUDE_DIRECTORIES(${dir})
SET(INC_FLAGS "${INC_FLAGS} -I'${dir}'")
ENDIF()
ENDFOREACH()
STRING(FIND "${INC_FLAGS}" "${CMAKE_INSTALL_PREFIX}/include" INC_POS)
IF( "${INC_POS}" STREQUAL "-1" )
SET(INC_FLAGS "${INC_FLAGS} -I'${CMAKE_INSTALL_PREFIX}/include'")
ENDIF()
SET(LIB_PATH "${LIB_PATH};${CMAKE_INSTALL_PREFIX}/lib")
LIST(REMOVE_DUPLICATES LIB_PATH)
SET(LIB_FLAGS "-Wl,-rpath,.")
FOREACH(dir ${LIB_PATH})
IF(EXISTS ${dir})
LINK_DIRECTORIES(${dir})
SET(LIB_FLAGS "${LIB_FLAGS} -Wl,-rpath,'${dir}' -L'${dir}'")
ENDIF()
ENDFOREACH()
STRING(FIND "${LIB_FLAGS}" "${CMAKE_INSTALL_PREFIX}/lib" LIB_POS)
IF( "${LIB_POS}" STREQUAL "-1" )
SET(LIB_FLAGS "${LIB_FLAGS} -Wl,-rpath,'${CMAKE_INSTALL_PREFIX}/lib' -L'${CMAKE_INSTALL_PREFIX}/lib'")
ENDIF()
#-----------------------------------------------
# HepLib - Include & Configure
#-----------------------------------------------
INCLUDE_DIRECTORIES(BEFORE MB)
INCLUDE_DIRECTORIES(BEFORE DE)
INCLUDE_DIRECTORIES(BEFORE IBP)
INCLUDE_DIRECTORIES(BEFORE SD)
INCLUDE_DIRECTORIES(BEFORE Hep)
INCLUDE_DIRECTORIES(BEFORE Qgraf)
INCLUDE_DIRECTORIES(BEFORE EX)
INCLUDE_DIRECTORIES(BEFORE BASIC)
INCLUDE_DIRECTORIES(BEFORE INC)
CONFIGURE_FILE(INC/Init.cpp "${CMAKE_CURRENT_BINARY_DIR}/OINC/Init.cpp" @ONLY)
CONFIGURE_FILE(INC/HepLib.pc "${CMAKE_CURRENT_BINARY_DIR}/OINC/HepLib.pc" @ONLY)
CONFIGURE_FILE(INC/heplib++ "${CMAKE_CURRENT_BINARY_DIR}/OINC/heplib++" @ONLY)
CONFIGURE_FILE(INC/update "${CMAKE_CURRENT_BINARY_DIR}/OINC/update" @ONLY)
#-----------------------------------------------
# External Libraries
#-----------------------------------------------
LINK_LIBRARIES(ginac cln Minuit2 mpfr qhullstatic quadmath gomp cubaq)
#-----------------------------------------------
PROJECT(HepLib VERSION 1.0 LANGUAGES C CXX Fortran)
AUX_SOURCE_DIRECTORY(BASIC BASIC_SRCS)
AUX_SOURCE_DIRECTORY(SD SD_SRCS)
AUX_SOURCE_DIRECTORY(Hep Hep_SRCS)
AUX_SOURCE_DIRECTORY(Qgraf Qgraf_SRCS)
AUX_SOURCE_DIRECTORY(IBP IBP_SRCS)
AUX_SOURCE_DIRECTORY(EX EX_SRCS)
AUX_SOURCE_DIRECTORY("${CMAKE_CURRENT_BINARY_DIR}/OINC" OINC_SRCS)
#-----------------------------------------------
# install library HepLib
#-----------------------------------------------
ADD_LIBRARY(HepLib SHARED ${BASIC_SRCS} ${Hep_SRCS} ${Qgraf_SRCS} ${IBP_SRCS} ${SD_SRCS} ${EX_SRCS} ${OINC_SRCS})
SET_TARGET_PROPERTIES(HepLib PROPERTIES INSTALL_RPATH "${LIB_PATH}")
INSTALL(TARGETS HepLib DESTINATION lib)
#-----------------------------------------------
# install header
#-----------------------------------------------
INSTALL(FILES "${CMAKE_SOURCE_DIR}/INC/HepLib.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/BASIC/BASIC.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/MB/MB.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/IBP/IBP.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/SD/SD.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/SD/NFunctions.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/SD/mpreal.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/Hep/HEP.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/Qgraf/QGRAF.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/EX/HepLibW.h" DESTINATION include)
INSTALL(FILES "${CMAKE_SOURCE_DIR}/EX/HepLib.i" DESTINATION include)
#-----------------------------------------------
# install bin/*
#-----------------------------------------------
FILE(GLOB cpp_files "bin/*.cpp")
FOREACH(exe_cpp ${cpp_files})
GET_FILENAME_COMPONENT(exe ${exe_cpp} NAME_WE)
ADD_EXECUTABLE(${exe} bin/${exe}.cpp)
SET_TARGET_PROPERTIES(${exe} PROPERTIES INSTALL_RPATH "${LIB_PATH}")
TARGET_LINK_LIBRARIES(${exe} HepLib)
INSTALL(TARGETS ${exe} DESTINATION bin)
ENDFOREACH()
INSTALL(TARGETS qgraf DESTINATION bin)
INSTALL(FILES "${CMAKE_CURRENT_BINARY_DIR}/OINC/HepLib.pc" DESTINATION "lib/pkgconfig")
INSTALL(FILES "${CMAKE_CURRENT_BINARY_DIR}/OINC/heplib++" DESTINATION "bin"
PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ)
INSTALL(FILES "${CMAKE_CURRENT_BINARY_DIR}/OINC/update" DESTINATION "bin"
PERMISSIONS OWNER_EXECUTE OWNER_WRITE OWNER_READ GROUP_EXECUTE GROUP_READ)
#-----------------------------------------------
# OpenMP support
#-----------------------------------------------
SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")
SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fopenmp")
DE/det.cpp deleted 100644 → 0
View file @ 8abf1d7e
#include "functions.cpp"
int main(int argc, char *argv[]) {
const char *x_name = "x";
const char *ep_name = "ep";
const char *m_path = "mat.m";
const char *f_path = "fmat.m";
const char *n_path = "nmat.m";
const char *s_path = "smat.m";
const char *t_path = "";
numeric ep_ex;
bool is_ep = false;
for (int opt; (opt = getopt(argc, argv, "x:e:m:f:n:s:t:N:")) != -1;) {
switch (opt) {
case 'x': x_name = optarg; break;
case 'e': ep_name = optarg; break;
case 'm': m_path = optarg; break;
case 'f': f_path = optarg; break;
case 'n': n_path = optarg; break;
case 's': s_path = optarg; break;
case 't': t_path = optarg; break;
case 'N': ep_ex = numeric(optarg); is_ep = true; break;
default: return 1;
}
}
argc -= optind;
argv += optind;
if((argc == 1) && !strcmp(argv[0], "fuchsify")) {
symbol x(x_name), ep(ep_name);
parser reader(symtab({{ep_name, ep}, {x_name, x}}));
matrix mat = load_matrix(m_path, reader);
if(is_ep) matrix_map_inplace(mat, [&](auto &&e) { return e.subs(ep == ep_ex); });
cout << "---------------------------------------" << endl;
cout << mat.rows() << " x " << mat.cols() << " matrix loaded from " << m_path << endl;
auto mt = fuchsify(mat, x);
save_matrix(f_path, mt.first);
if(strlen(t_path)>0) save_matrix(t_path, mt.second);
cout << "---------------------------------------" << endl;
cout << "f-matrix saved to " << f_path << endl;
if(strlen(t_path)>0) cout << "t-matrix saved to " << t_path << endl;
} else if((argc == 1) && !strcmp(argv[0], "normalize")) {
symbol x(x_name), ep(ep_name);
parser reader(symtab({{ep_name, ep}, {x_name, x}}));
matrix mat = load_matrix(f_path, reader);
if(is_ep) matrix_map_inplace(mat, [&](auto &&e) { return e.subs(ep == ep_ex); });
cout << "---------------------------------------" << endl;
cout << mat.rows() << " x " << mat.cols() << " matrix loaded from " << f_path << endl;
auto mt = normalize(mat, x);
save_matrix(n_path, mt.first);
if(strlen(t_path)>0) save_matrix(t_path, mt.second);
cout << "---------------------------------------" << endl;
cout << "n-matrix saved to " << n_path << endl;
if(strlen(t_path)>0) cout << "t-matrix saved to " << t_path << endl;
} else if ((argc == 1 || argc == 2) && !strcmp(argv[0], "shearing")) {
symbol x(x_name), ep(ep_name);
parser reader(symtab({{ep_name, ep}, {x_name, x}}));
matrix mat = load_matrix(f_path, reader);
if(is_ep) matrix_map_inplace(mat, [&](auto &&e) { return e.subs(ep == ep_ex); });
int epN = -19790923;
if(argc == 2) epN = stoi(argv[1]);
cout << "---------------------------------------" << endl;
cout << mat.rows() << " x " << mat.cols() << " matrix loaded from " << f_path << endl;
auto mt = shearing(mat, x, ep, epN);
save_matrix(n_path, mt.first);
if(strlen(t_path)>0) save_matrix(t_path, mt.second);
cout << "---------------------------------------" << endl;
cout << "n-matrix saved to " << n_path << endl;
if(strlen(t_path)>0) cout << "t-matrix saved to " << t_path << endl;
} else if ((argc > 1) && !strcmp(argv[0], "dess")) {
Digits = 100;
symbol x(x_name), ep(ep_name);
parser reader(symtab({{ep_name, ep}, {x_name, x}}));
int xN = 10, epN = -19790923;
ex x0 = x;
if(argc>1) xN = stoi(argv[1]);
if(argc>2) epN = stoi(argv[2]);
if(argc>3) x0 = numeric(argv[3]);
matrix mat = load_matrix(n_path, reader);
cout << "---------------------------------------" << endl;
cout << mat.rows() << " x " << mat.cols() << " matrix loaded from " << f_path << endl;
auto smat = dess(mat, x, ep, x0, xN, epN);
cout << "---------------------------------------" << endl;
save_matrix(s_path, smat);
cout << "s-matrix saved to " << s_path << endl;
} else if ((argc == 2) && !strcmp(argv[0], "show")) {
symbol x(x_name), ep(ep_name);
parser reader(symtab({{ep_name, ep}, {x_name, x}}));
m_path = argv[1];
matrix mat = load_matrix(m_path, reader);
if(is_ep) matrix_map_inplace(mat, [&](auto &&e) { return e.subs(ep == ep_ex); });
cout << "---------------------------------------" << endl;
cout << mat.rows() << " x " << mat.cols() << " matrix loaded from " << m_path << endl;
cout << "---------------------------------------" << endl;
int pr = prank(mat, x);
matrix a0 = a0_matrix(mat, x, pr);
a0 = normal(a0);
cout << "matrix information:" << endl;
cout << " Poincare Rank: " << pr << endl;
auto ev_map = eigenvalues(a0);
vector<vector<ex>> ev_groups;
for(auto &kv : ev_map) {
auto ev = kv.first;
bool in_g = false;
for(auto &gi : ev_groups) {
if(in_g) break;
for(auto &ev_i : gi) {
auto diff_ev = normal(ev_i-ev);
if(is_a<numeric>(diff_ev) && ex_to<numeric>(diff_ev).is_integer()) {
in_g = true;
gi.push_back(ev);
break;
}
}
}
if(!in_g) {
vector<ex> gi_new;
gi_new.push_back(ev);
ev_groups.push_back(gi_new);
}
}
cout << " eigen values of A0 summary:" << endl;
for(auto &gi : ev_groups) {
sort(gi.begin(),gi.end(),[&](const auto &a, const auto &b){
return normal(a-b).info(info_flags::positive);
});
ostringstream ostr;
cout << " ";
for(auto &ev_i : gi) {
ostr << ev_i << "[" << ev_map[ev_i] << "], ";
}
string str = ostr.str();
cout << str.substr(0,str.size()-3) << endl;
}
cout << "---------------------------------------" << endl;
cout << " jordan block of A0:" << endl;
ostringstream ostr;
int pn = 0;
for(auto kv : jordan(a0).second) {
pn++;
ostr << kv.first << "[" << kv.second << "], ";
if(pn>=4) {
string str = ostr.str();
ostr.str("");
ostr.clear();
cout << " " << str.substr(0,str.size()-3) << endl;
pn = 0;
}
}
} else if(false) {
cout << "---------------------------------------" << endl;
cout << "missing or invalid arguments!" << endl;
} else {
cout << "det fuchsify / det shearing / det dess xN epN x0" << endl;
}
cout << "\r" << flush;
cout << "---------------------------------------" << endl;
return 0;
}
DE/fuchsia.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
DE/functions.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
EX/HepLib.i deleted 100644 → 0
View file @ 8abf1d7e
/**
* @file
* @brief Wrap Class for SWIG
*/
%module(directors="1") HepLib
%feature("director") MapFunction;
%feature("allowexcept");
%{
#define SWIG_FILE_WITH_INIT
#include "HepLibW.h"
%}
//--------------------------------------------------------------------
%include std_string.i
%include exception.i
%include std_vector.i
%exception {
try {
$action
} catch(const std::exception& e) {
SWIG_exception(SWIG_RuntimeError, e.what());
} catch (const std::string& e) {
SWIG_exception(SWIG_RuntimeError, e.c_str());
} catch(...) {
SWIG_exception(SWIG_RuntimeError, "HepLib unkonwn exception.");
}
}
namespace std {
%template(expr_vector) vector<expr>;
%template(int_vector) vector<int>;
}
class MapFunction {
public:
MapFunction();
virtual ~MapFunction();
virtual expr map(const expr &e);
expr operator() (const expr &e);
};
class expr {
public:
expr(int i);
expr(const std::string &s);
expr operator+(const expr &e);
expr operator-(const expr &e);
expr operator*(const expr &e);
expr operator/(const expr &e);
expr operator-();
expr operator==(const expr &e);
expr operator+(const int i);
expr operator-(const int i);
expr operator*(const int i);
expr operator/(const int i);
unsigned int nops();
expr op(unsigned int i);
void let_op(unsigned int i, expr e);
expr expand();
expr normal();
expr factor();
expr series(const expr &s, int o);
expr subs(const std::vector<expr> &ev);
expr subs(const expr &e);
bool match(const expr &e);
bool isSymbol();
bool isVector();
bool isIndex();
bool isPair();
bool isDGamma();
bool info(std::string sflags);
expr map(MapFunction &mf);
std::string str();
std::string __str__();
%pythoncode %{
def __radd__(self, other):
return expr(other) + self
def __rsub__(self, other):
return expr(other) - self
def __rmul__(self, other):
return expr(other) * self
def __rdiv__(self, other):
return expr(other) / self
%}
};
extern expr expand(const expr &e);
extern expr normal(const expr &e);
extern expr factor(const expr &e);
extern expr series(const expr &e, const expr &s, int o);
extern expr subs(const expr &e, const std::vector<expr> &ev);
extern expr subs(const expr &e1, const expr &e2);
extern expr pow(const expr &e1, const expr &e2);
extern expr pow(const expr &e, const int n);
extern expr Index(const std::string &s);
extern expr Vector(const std::string &s);
extern expr Symbol(const std::string &s);
extern expr SP(const expr &e1, const expr &e2);
extern expr GAS(const expr &e);
extern expr TR(const expr &e);
extern expr form(const expr &e);
extern void letSP(const expr &e1, const expr &e2, const expr &e12);
extern expr call(const std::string func, const std::vector<expr> &ev);
extern expr call(const std::string func, const expr &e);
extern expr w(const int wi);
extern expr x(const int i);
extern expr y(const int i);
extern expr z(const int i);
// Integral
%warnfilter(509) Integral;
class Integral {
public:
int epN;
int epsN;
int verb;
Integral();
void Functions(const std::vector<expr> &ev);
void Propagators(const std::vector<expr> &ev, const std::vector<expr> &loops={}, const std::vector<expr> &tloops={});
void Replacements(const std::vector<expr> &ev, int lt=0);
void Exponents(const std::vector<expr> &ev);
void Exponents(const std::vector<int> &ev);
void Evaluate();
std::string str();
std::string __str__();
};
EX/HepLibW.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
EX/HepLibW.h deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
EX/exfactor.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
HEP/Apart.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
HEP/Basic.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
HEP/DGamma.cpp deleted 100644 → 0
View file @ 8abf1d7e
This diff is collapsed.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Provided by Institute of Theoretical Physics, Chinese Academy of Sciences. Please contact support@itp.ac.cn for more information.