Skip to content

Instantly share code, notes, and snippets.

@veprbl

veprbl/blackhat_0.9.9_fix.patch

Last active August 29, 2015 14:14
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save veprbl/b93e5487dfb8a8e6df2a to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save veprbl/b93e5487dfb8a8e6df2a to your computer and use it in GitHub Desktop.
Download ZIP
Raw
blackhat_0.9.9_fix.patch
diff --git a/Makefile.in b/Makefile.in
index 7f18f5c..6d9b237 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -876,10 +876,6 @@ slow:
install-data-local:
mkdir -p $(pkgdatadir)
-@HAS_SVNVERSION_TRUE@ if ! test `svnversion $(srcdir)` = "exported" ; then\
-@HAS_SVNVERSION_TRUE@ svnversion $(srcdir) > $(pkgdatadir)/svnversion ;\
-@HAS_SVNVERSION_TRUE@ svn info $(srcdir) | grep URL | grep -o "privatesvn/.*" | cut -d / -f2- > $(pkgdatadir)/svnurl ;\
-@HAS_SVNVERSION_TRUE@ fi
if test -f $(srcdir)/share/svnversion_dist ; then cp $(srcdir)/share/svnversion_dist $(pkgdatadir)/svnversion_dist; fi
if test -f $(srcdir)/share/svnurl_dist ; then cp $(srcdir)/share/svnurl_dist $(pkgdatadir)/svnurl_dist ; fi
diff --git a/my_programs/Makefile.in b/my_programs/Makefile.in
index 2ace3bd..ad9f89a 100644
--- a/my_programs/Makefile.in
+++ b/my_programs/Makefile.in
@@ -425,7 +425,7 @@ uninstall-am:
mostlyclean-libtool pdf pdf-am ps ps-am uninstall uninstall-am
-include $(DEPDIR)/*.Po
+#include $(DEPDIR)/*.Po
%.o: $(abs_srcdir)/%.cpp
test -f $(DEPDIR)/$*.Po || mkdir -p .deps ; echo "#dummy" > $(DEPDIR)/$*.Po
diff --git a/src/BH_typedefs.h b/src/BH_typedefs.h
index f5abd95..80f1a07 100644
--- a/src/BH_typedefs.h
+++ b/src/BH_typedefs.h
@@ -5,7 +5,94 @@
#define BH_TYPEDEFS_H_
#include <complex>
+
+#if defined(_LIBCPP_VERSION)
+#include <cfloat>
+
+#ifdef _QD_DD_REAL_H
+#error dd_real.h must not be included before this
+#endif
+
+#ifdef _QD_QD_REAL_H
+#error qd_real.h must not be included before this
+#endif
+
+#include "qd/qd_config.h"
+#ifdef QD_INLINE
+#define QD_INLINE_BACKUP
+#undef QD_INLINE
+#endif
+
+class dd_real;
+extern int to_int(const dd_real &a);
+
+#define sloppy_add(a, b) \
+ not_a_method1() {throw "Not implemented";}; \
+ static dd_real sloppy_add(a, b); \
+ operator int() {return to_int(*this);};
+#include "qd/dd_real.h"
+#undef sloppy_add
+
+class qd_real;
+extern int to_int(const qd_real &a);
+
+#define accurate_mul(a, b) \
+ not_a_method2() {throw "Not implemented";}; \
+ static qd_real accurate_mul(a, b); \
+ operator int() {return to_int(*this);};
#include "qd/qd_real.h"
+#undef accurate_mul
+
+#ifdef QD_INLINE_BACKUP
+#undef QD_INLINE_BACKUP
+#define QD_INLINE 1
+#include "qd/dd_inline.h"
+#include "qd/qd_inline.h"
+#endif
+
+inline bool operator==(const qd_real &a, int b) {
+ return (a[0] == b && a[1] == 0.0 && a[2] == 0.0 && a[3] == 0.0);
+}
+inline bool operator==(const dd_real &a, int b) {
+ return (a.x[0] == b && a.x[1] == 0.0);
+}
+
+_LIBCPP_BEGIN_NAMESPACE_STD
+
+ inline _LIBCPP_INLINE_VISIBILITY qd_real copysign(qd_real __x, qd_real __y) _NOEXCEPT
+ {return (__x.is_positive() != __y.is_positive()) ? __x : -__x;}
+ inline _LIBCPP_INLINE_VISIBILITY dd_real copysign(dd_real __x, dd_real __y) _NOEXCEPT
+ {return (__x.is_positive() != __y.is_positive()) ? __x : -__x;}
+ inline _LIBCPP_INLINE_VISIBILITY qd_real fmax(qd_real __x, qd_real __y) _NOEXCEPT
+ {return max(__x, __y);}
+ inline _LIBCPP_INLINE_VISIBILITY dd_real fmax(dd_real __x, dd_real __y) _NOEXCEPT
+ {return max(__x, __y);}
+ inline _LIBCPP_INLINE_VISIBILITY qd_real hypot(qd_real __x, qd_real __y) _NOEXCEPT
+ {return sqrt(__x*__x + __y*__y);}
+ inline _LIBCPP_INLINE_VISIBILITY dd_real hypot(dd_real __x, dd_real __y) _NOEXCEPT
+ {return sqrt(__x*__x + __y*__y);}
+ inline _LIBCPP_INLINE_VISIBILITY qd_real logb(qd_real __x) _NOEXCEPT
+ {return std::floor(to_double(log(__x) / std::log(FLT_RADIX)));}
+ inline _LIBCPP_INLINE_VISIBILITY dd_real logb(dd_real __x) _NOEXCEPT
+ {return std::floor(to_double(log(__x) / std::log(FLT_RADIX)));}
+ // assume FLT_RADIX == 2
+ static_assert(FLT_RADIX == 2, "FLT_RADIX == 2");
+ inline _LIBCPP_INLINE_VISIBILITY qd_real scalbn(qd_real __x, int __y) _NOEXCEPT
+ {return ldexp(__x, __y);}
+ inline _LIBCPP_INLINE_VISIBILITY dd_real scalbn(dd_real __x, int __y) _NOEXCEPT
+ {return ldexp(__x, __y);}
+
+_LIBCPP_END_NAMESPACE_STD
+
+using std::__1::copysign;
+using std::__1::hypot;
+using std::__1::fmax;
+using std::__1::logb;
+using std::__1::scalbn;
+#else
+#include "qd/dd_real.h"
+#include "qd/qd_real.h"
+#endif // if defined(_LIBCPP_VERSION)
namespace BH {
diff --git a/src/BerendsGiele.cpp b/src/BerendsGiele.cpp
index d54b84d..edd7191 100644
--- a/src/BerendsGiele.cpp
+++ b/src/BerendsGiele.cpp
@@ -13,6 +13,7 @@
*/
+#include "BH_typedefs.h"
#include "mom_conf.h"
#include "BerendsGiele.h"
#include "BerendsGiele_impl.h"
diff --git a/src/IR_checked_amplitudes.cpp b/src/IR_checked_amplitudes.cpp
index 906a1f6..beccc77 100644
--- a/src/IR_checked_amplitudes.cpp
+++ b/src/IR_checked_amplitudes.cpp
@@ -1,5 +1,6 @@
//#include "OneLoopHelAmpl.h"
//#include "BH_A0.h"
+#include "BH_typedefs.h"
#include "qd_suppl.h"
#include <typeinfo>
#include "IR_checked.h"
diff --git a/src/IR_checked_amplitudes_ep.cpp b/src/IR_checked_amplitudes_ep.cpp
index 867e77f..e50e7cc 100644
--- a/src/IR_checked_amplitudes_ep.cpp
+++ b/src/IR_checked_amplitudes_ep.cpp
@@ -6,6 +6,7 @@
*/
//#include "OneLoopHelAmpl.h"
//#include "BH_A0.h"
+#include "BH_typedefs.h"
#include "qd_suppl.h"
#include <typeinfo>
#include "IR_checked.h"
diff --git a/src/IR_checked_cut_part.cpp b/src/IR_checked_cut_part.cpp
index 1127da9..ddc4d83 100644
--- a/src/IR_checked_cut_part.cpp
+++ b/src/IR_checked_cut_part.cpp
@@ -97,7 +97,7 @@ template <class T> bool IR_check_Cut_Part(complex<T> Atree,complex<T> single_pol
}
int int_part=int(to_double(imag_ratio));
- return ((0<imag_ratio-int_part)&&(imag_ratio-int_part<T(2)*tolerance)&&(inverse_tree>tolerance/T(1000)));
+ return ((T(0)<imag_ratio-T(int_part))&&(imag_ratio-T(int_part)<T(2)*tolerance)&&(inverse_tree>tolerance/T(1000)));
}
diff --git a/src/Interface/BH_interface.cpp b/src/Interface/BH_interface.cpp
index 7476c17..301bc62 100644
--- a/src/Interface/BH_interface.cpp
+++ b/src/Interface/BH_interface.cpp
@@ -5,6 +5,7 @@
* Author: daniel
*/
+#include "BH_typedefs.h"
#include "BH_interface.h"
#include <iostream>
#include <string>
diff --git a/src/Series.cpp b/src/Series.cpp
index f14638a..ac80a78 100644
--- a/src/Series.cpp
+++ b/src/Series.cpp
@@ -399,12 +399,12 @@ template std::ostream& operator<<(std::ostream&, const Series<TYPE>&);
template <> TYPE Series<TYPE>::zero=TYPE(0);\
template <> TYPE Series<TYPE>::infinity=std::numeric_limits<TYPE>::quiet_NaN();
-INSTANTIATE(C)
-INSTANTIATE(CHP)
-INSTANTIATE(CVHP)
INSTANTIATE_STATIC(C)
INSTANTIATE_STATIC(CHP)
INSTANTIATE_STATIC(CVHP)
+INSTANTIATE(C)
+INSTANTIATE(CHP)
+INSTANTIATE(CVHP)
#if BH_USE_GMP
INSTANTIATE(std::complex<RGMP>)
diff --git a/src/Tree.h b/src/Tree.h
index 31bc827..6de76c3 100644
--- a/src/Tree.h
+++ b/src/Tree.h
@@ -49,7 +49,7 @@ template<class T> std::complex<T>
const std::vector<int>& arg /* indices of arguments */,
const std::vector<particle_ID>& leg,
int start, int end /* indices into the vectors */,
- int offshellMass = -1,
+ int offshellMass = defaultMass,
const std::vector<int>& massValue = empty);
template<class T> std::complex<T>
@@ -63,7 +63,7 @@ template<class T> std::complex<T>
const std::vector<int>& vectorK /* momenta */ = empty,
const std::vector<int>& polarization = empty,
const std::vector<int>& coupleTo /* quark flavor */ = empty,
- int offshellMass = -1,
+ int offshellMass = defaultMass,
const std::vector<int>& massValue = empty);
diff --git a/src/cached_OLHA.cpp b/src/cached_OLHA.cpp
index 38b1d97..219d625 100644
--- a/src/cached_OLHA.cpp
+++ b/src/cached_OLHA.cpp
@@ -271,14 +271,14 @@ template <class OLHA> Cached_OLHA_factory_impl<OLHA>::~Cached_OLHA_factory_impl(
template <class OLHA > void Cached_OLHA_factory_impl<OLHA>::print_state(){
_MESSAGE("=-=-=-=-=-=-=-=-=-=-= Cached_OLHA_factory =-=-=-=-=-=-=-=-=-=-= " );
- for ( map<const pro_cs,Cached_OLHA*>::iterator it=d_amplitudes.begin();it!=d_amplitudes.end();it++){
+ for ( map<pro_cs,Cached_OLHA*>::iterator it=d_amplitudes.begin();it!=d_amplitudes.end();it++){
(*it).second->print_stat();
};
_MESSAGE("=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= " );
}
template <class OLHA > void Cached_OLHA_factory_impl<OLHA>::fill_process_list(std::vector<std::pair<process,color_structure> >& li){
- for ( map<const pro_cs,Cached_OLHA*>::iterator it=d_amplitudes.begin();it!=d_amplitudes.end();it++){
+ for ( map<pro_cs,Cached_OLHA*>::iterator it=d_amplitudes.begin();it!=d_amplitudes.end();it++){
li.push_back( std::make_pair((*it).second->get_process(),(*it).second->get_color_structure() ) );
};
}
@@ -331,7 +331,7 @@ template <class OLHA> Cached_OLHA_user* Cached_OLHA_factory_impl<OLHA>::new_OLHA
Cached_OLHA* new_COLHA = new Cached_OLHA(new_OLHAb);
d_amplitudes.insert(pair<pro_cs,Cached_OLHA*>(pc,new_COLHA));
if(conjQ) return new Cached_OLHA_user_conjugate(new_COLHA,new_COLHA->add(ind),sign);
- else return Cached_OLHA_user_normal(new_COLHA,new_COLHA->add(ind));
+ else return new Cached_OLHA_user_normal(new_COLHA,new_COLHA->add(ind));
} else {
if(conjQ) return new Cached_OLHA_user_conjugate((*it).second,(*it).second->add(ind),sign);
else return new Cached_OLHA_user_normal((*it).second,(*it).second->add(ind));
@@ -417,7 +417,7 @@ template <> Cached_OLHA_user* Cached_OLHA_factory_impl<OneLoopHelAmpl>::new_OLHA
}
template <class OLHA> void Cached_OLHA_factory_impl<OLHA>::refresh(momentum_configuration<R>& mc, int mu_index){
- for ( map<const pro_cs,Cached_OLHA*>::iterator it=d_amplitudes.begin();it!=d_amplitudes.end();it++){
+ for ( map<pro_cs,Cached_OLHA*>::iterator it=d_amplitudes.begin();it!=d_amplitudes.end();it++){
(*it).second->refresh(mc,mu_index);
};
}
diff --git a/src/cached_THA.cpp b/src/cached_THA.cpp
index 377bb63..a2d3e66 100644
--- a/src/cached_THA.cpp
+++ b/src/cached_THA.cpp
@@ -323,7 +323,7 @@ Cached_THA_user* Cached_THA_factory::new_THA(process pro,const std::vector<int>&
/* end eval param caching */
- map<process,Cached_THA*>::iterator it = d_amplitudes.find(pro);
+ map<const process,Cached_THA*>::iterator it = d_amplitudes.find(pro);
if (it == d_amplitudes.end()){
// new part
diff --git a/src/cut_Darren.h b/src/cut_Darren.h
index 0b42769..d057ec1 100644
--- a/src/cut_Darren.h
+++ b/src/cut_Darren.h
@@ -344,9 +344,9 @@ public:
virtual CHP get_sub_terms(momentum_configuration<dd_real>& mc, const std::vector<int>& ind, coeffparam<RHP,TriangleSpecs::CPOINTS>& tp){return CHP(0,0);}
virtual CVHP get_sub_terms(momentum_configuration<qd_real>& mc, const std::vector<int>& ind, coeffparam<RVHP,TriangleSpecs::CPOINTS>& tp){return CVHP(0,0);}
- virtual C get_sub_terms(const eval_param<R>& ep, coeffparam<R,TriangleSpecs::CPOINTS>& tp) =0;
- virtual CHP get_sub_terms(const eval_param<RHP>& ep, coeffparam<RHP,TriangleSpecs::CPOINTS>& tp) =0;
- virtual CVHP get_sub_terms(const eval_param<RVHP>& ep, coeffparam<RVHP,TriangleSpecs::CPOINTS>& tp) =0;
+ virtual C get_sub_terms(const eval_param<R>& ep, coeffparam<R,TriangleSpecs::CPOINTS>& tp) {throw "Not implemented";};
+ virtual CHP get_sub_terms(const eval_param<RHP>& ep, coeffparam<RHP,TriangleSpecs::CPOINTS>& tp) {throw "Not implemented";};
+ virtual CVHP get_sub_terms(const eval_param<RVHP>& ep, coeffparam<RVHP,TriangleSpecs::CPOINTS>& tp) {throw "Not implemented";};
#ifdef BH_USE_GMP
virtual CGMP eval(momentum_configuration<RGMP>&,const std::vector<int>&);
virtual CGMP eval(const eval_param<RGMP>& ep);
diff --git a/src/cut_eval/C_2q1g2l_eval.cpp b/src/cut_eval/C_2q1g2l_eval.cpp
index b60dab0..7082bce 100644
--- a/src/cut_eval/C_2q1g2l_eval.cpp
+++ b/src/cut_eval/C_2q1g2l_eval.cpp
@@ -6,6 +6,7 @@
*/
+#include "BH_typedefs.h"
#include "C_2q1g2l_eval.h"
#include "tree_amp.h"
diff --git a/src/cut_eval/C_2q3g2l_eval.cpp b/src/cut_eval/C_2q3g2l_eval.cpp
index 1bd3400..79d285d 100644
--- a/src/cut_eval/C_2q3g2l_eval.cpp
+++ b/src/cut_eval/C_2q3g2l_eval.cpp
@@ -6,6 +6,7 @@
*/
+#include "BH_typedefs.h"
#include "C_2q3g2l_eval.h"
#include "tree_amp.h"
diff --git a/src/eval_param.cpp b/src/eval_param.cpp
index 1fbffdb..969abea 100644
--- a/src/eval_param.cpp
+++ b/src/eval_param.cpp
@@ -358,15 +358,6 @@ template std::vector<Cmom<RHP> >* extend_momenta<R,RHP>(const eval_param<R>& ep)
template std::vector<Cmom<RVHP> >* extend_momenta<R,RVHP>(const eval_param<R>& ep);
template std::vector<Cmom<RVHP> >* extend_momenta<RHP,RVHP>(const eval_param<RHP>& ep);
-template const Cmom<R> eval_param<R>::_ep_quark_ref;
-template const Cmom<RHP> eval_param<RHP>::_ep_quark_ref;
-template const Cmom<RVHP> eval_param<RVHP>::_ep_quark_ref;
-
-
-template void eval_param<R>::update(momentum_configuration<R>& mc, const std::vector<int>& ind);
-template void eval_param<RHP>::update(momentum_configuration<RHP>& mc, const std::vector<int>& ind);
-template void eval_param<RVHP>::update(momentum_configuration<RVHP>& mc, const std::vector<int>& ind);
-
template std::ostream& operator<<(std::ostream& s,const eval_param<R>& ep);
diff --git a/src/eval_param.h b/src/eval_param.h
index 67088b4..cab221f 100644
--- a/src/eval_param.h
+++ b/src/eval_param.h
@@ -8,6 +8,7 @@
#ifndef EVAL_PARAM_H_
#define EVAL_PARAM_H_
+#include "BH_typedefs.h"
#include "spinor.h"
#include "mom_conf.h"
#include "process.h"
diff --git a/src/extend.cpp b/src/extend.cpp
index b8e7c50..4f8da4c 100644
--- a/src/extend.cpp
+++ b/src/extend.cpp
@@ -21,6 +21,7 @@
*/
+#include "BH_typedefs.h"
#include "mom_conf.h"
#include "BH_debug.h"
@@ -62,7 +63,7 @@ inline bool IsZero(const RVHP& value)
{return(value < zeroThresholdVHP);}
template<class T> inline T Sign(const complex<T>& value)
-{if (IsZero(imag(value)) and real(value) < 0) return T(-1); else return T(1);}
+{if (IsZero(imag(value)) and real(value) < T(0)) return T(-1); else return T(1);}
/* Transform the vector "p" from the frame in which "transform" is at rest
(if it is positive-mass) or has zero energy component (if it is
@@ -71,7 +72,7 @@ template<class Te,class T0> static momentum<complex<Te> >
Boost(const momentum<complex<T0> >& p,
const momentum<complex<Te> >& transform)
{complex<Te> massSq = transform*transform;
- if (IsZero(imag(massSq)) and real(massSq) < 0)
+ if (IsZero(imag(massSq)) and real(massSq) < Te(0))
{complex<Te> mass = sqrt(-transform*transform);
momentum<complex<Te> > base(Te(0),transform.X(),transform.Y(),transform.Z());
return Boost(p,momentum<complex<Te> >(base*base/mass,Te(0),Te(0),Te(0))
@@ -363,7 +364,7 @@ template <class T_low,class T_high> momentum_configuration<T_high> extend_real(c
momenta.push_back(Cmom<T_high>(lambdat<T_high>(mc.Lt(ind[i])),lambda<T_high>(mc.L(ind[i]))));
}
else{
- if(real(mc.p(ind[i]).E())<=0){
+ if(real(mc.p(ind[i]).E())<=T_low(0)){
momenta.push_back(Cmom<T_high>(lambdat<T_high>(mc.Lt(ind[i])),
lambda<T_high>(-conj(Z4*lambdat<T_high>(mc.Lt(ind[i]))),conj(Z3*lambdat<T_high>(mc.Lt(ind[i]))))));
}
@@ -389,7 +390,7 @@ template <class T_low,class T_high> momentum_configuration<T_high> extend_real(c
knm1c=Cmom<T_high>(lambdat<T_high>(mc.Lt(ind[n-2])),lambda<T_high>(mc.L(ind[n-2])));
}
else{
- if(real(mc.p(ind[n-2]).E())<=0){
+ if(real(mc.p(ind[n-2]).E())<=T_low(0)){
knm1c=Cmom<T_high>(lambdat<T_high>(mc.Lt(ind[n-2])),
lambda<T_high>(-conj(Z4*lambdat<T_high>(mc.Lt(ind[n-2]))),conj(Z3*lambdat<T_high>(mc.Lt(ind[n-2])))));
}
@@ -414,7 +415,7 @@ template <class T_low,class T_high> momentum_configuration<T_high> extend_real(c
knm1f=Cmom<T_high>(sqrtt*lambdat<T_high>(mc.Lt(ind[n-2])),sqrtt*lambda<T_high>(mc.L(ind[n-2])));
}
else{
- if(real(mc.p(ind[n-2]).E())<=0){
+ if(real(mc.p(ind[n-2]).E())<=T_low(0)){
knm1f=Cmom<T_high>(sqrtt*lambdat<T_high>(mc.Lt(ind[n-2])),
sqrtt*lambda<T_high>(-conj(Z4*lambdat<T_high>(mc.Lt(ind[n-2]))),conj(Z3*lambdat<T_high>(mc.Lt(ind[n-2])))));
}
diff --git a/src/fact_momenta.cpp b/src/fact_momenta.cpp
index 49eaab9..abf02a7 100644
--- a/src/fact_momenta.cpp
+++ b/src/fact_momenta.cpp
@@ -148,7 +148,7 @@ template <class T> momentum_configuration<T>collkinematicsCS(int n,
vector<Cmom<T> > tempv(2);
tempv = decay(Ktot,0);
- if (tempv.at(0).E() == T(0))
+ if (tempv.at(0).E() == complex<T>(0,0))
return (collkinematicsCS(n,a,b,z,invmass));
// too many iterations, try again
@@ -239,7 +239,7 @@ template <class T> momentum_configuration<T>collkinematics(int n,
vector<Cmom<T> > tempv(2);
tempv = decay(Ktot,0);
- if (tempv.at(0).E() == T(0))
+ if (tempv.at(0).E() == complex<T>(0,0))
return (collkinematics(n,a,b,z,invmass));
// too many iterations, try again
@@ -343,7 +343,7 @@ template <class T> momentum_configuration<T> softkinematics(int n,
vector<Cmom<T> > tempv(2);
tempv = decay(Ktot,0);
- if (tempv.at(0).E() == T(0))
+ if (tempv.at(0).E() == complex<T>(0,0))
return (softkinematics(n,a,b,c,invmass));
// too many iterations, try again
@@ -466,7 +466,7 @@ template <class T> momentum_configuration<T> soft4kinematics(int n,
vector<Cmom<T> > tempv(2);
tempv = decay(Ktot,0);
- if (tempv.at(0).E() == T(0))
+ if (tempv.at(0).E() == complex<T>(0,0))
return (soft4kinematics(n,a,b,c,d,invmass));
temparray.insert(temparray.begin()+n1,tempv.at(0));
@@ -528,7 +528,7 @@ template <class T> vector<Cmom<T> > decay(Cmom<T> Ktot, int trial){
T alpha = auxy*auxy + auxz*auxz;
T beta = auxy*real(Ktot.Y()) + auxz*real(Ktot.Z());
T square = real(Ktot.E()*Ktot.E()-Ktot.X()*Ktot.X()-Ktot.Y()*Ktot.Y()-Ktot.Z()*Ktot.Z());
- T gamma = beta - 1/T(2)*(square);
+ T gamma = beta - T(0.5)*(square);
T delta = real(Ktot.E()*Ktot.E()) - real(Ktot.X()*Ktot.X());
T discrim = -alpha*delta + gamma*gamma;
if (discrim < 0.)
@@ -541,7 +541,7 @@ template <class T> vector<Cmom<T> > decay(Cmom<T> Ktot, int trial){
_DEBUG("Real decay achieved!\n");
T auxx = (real(Ktot.X())*gamma + real(Ktot.E())*sqrt(discrim))/delta;
T auxE = (auxx*real(Ktot.X()) + auxy*real(Ktot.Y()) + auxz*real(Ktot.Z()) - \
- 1/T(2)*(square))/real(Ktot.E());
+ T(0.5)*(square))/real(Ktot.E());
// 7/14/2008 -- check for dangerous configuration (instability in FHZ)
if (fabs(auxx-auxE) < 0.01)
{
@@ -570,7 +570,7 @@ template <class T> Cmom<T> randmom(T mass, short signE){
}
template <class T> T randdouble(void){
- return static_cast<T>(static_cast<T>(rand())/(RAND_MAX+T(1)));
+ return static_cast<T>(static_cast<T>(rand())/(T(RAND_MAX+1)));
}
// explicit instantiations
diff --git a/src/from_file.cpp b/src/from_file.cpp
index 9c90aaf..a747e9d 100644
--- a/src/from_file.cpp
+++ b/src/from_file.cpp
@@ -314,8 +314,8 @@ std::string GetParentDataDirectory(){
BH_DEBUG_MESSAGE3("Using ",BH_SOURCE_PATH+datafiles," as the parent data path.");
return BH_SOURCE_PATH+datafiles;
} else {
- BH_DEBUG_MESSAGE3("Using ",BH_INSTALL_PATH+datafiles," as the parent data path.");
- return BH_INSTALL_PATH+datafiles;
+ BH_DEBUG_MESSAGE3("Using ",BH_INSTALL_PATH+std::string("/share/blackhat")+datafiles," as the parent data path.");
+ return BH_INSTALL_PATH+std::string("/share/blackhat")+datafiles;
}
}
}
diff --git a/src/genkey.cc b/src/genkey.cc
index fb8b365..e866153 100644
--- a/src/genkey.cc
+++ b/src/genkey.cc
@@ -2,6 +2,7 @@
/* David A. Kosower, May 22, 2008 */
+#include "BH_typedefs.h"
#include "mom_conf.h"
#include "key.h"
#include "BH_error.h"
diff --git a/src/integrals.cpp b/src/integrals.cpp
index 2839f7a..19df52b 100644
--- a/src/integrals.cpp
+++ b/src/integrals.cpp
@@ -894,10 +894,10 @@ template <class T> complex<T> I4w4m(int order, momentum_configuration<T>& k,
cout << "imaginary rho" << endl;
#endif
- complex<T> arg1a = (1-lambda1+lambda2+rho)/T(2);
- complex<T> arg2a = (1-lambda1+lambda2-rho)/T(2);
- complex<T> arg1b = -(1-lambda1-lambda2-rho)/lambda1/T(2);
- complex<T> arg2b = -(1-lambda1-lambda2+rho)/lambda1/T(2);
+ complex<T> arg1a = (T(1)-lambda1+lambda2+rho)/T(2);
+ complex<T> arg2a = (T(1)-lambda1+lambda2-rho)/T(2);
+ complex<T> arg1b = -(T(1)-lambda1-lambda2-rho)/lambda1/T(2);
+ complex<T> arg2b = -(T(1)-lambda1-lambda2+rho)/lambda1/T(2);
complex<T> x1 = s*(arg1a-T(1))/m3sq, x2 = s*(arg2a-T(1))/m3sq;
complex<T> ln1 = ln(x1),
@@ -983,7 +983,7 @@ template <class T> complex<T> I4w4m(int order, momentum_configuration<T>& k,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
@@ -1338,7 +1338,7 @@ template <class T> complex<T> I4w4m(int order,
// But we must treat real & complex case for rho separately
// because of branch cut handling
- if (rhoSq >= 0)
+ if (rhoSq >= T(0))
{//T rho = sqrt(rhoSq);
complex<T> rho;
@@ -1417,7 +1417,7 @@ template <class T> complex<T> I4w4m(int order,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
@@ -1739,10 +1739,10 @@ template <class T> complex<T> I4w4m(int order,
cout << "imaginary rho" << endl;
#endif
- complex<T> arg1a = (1-lambda1+lambda2+rho)/T(2);
- complex<T> arg2a = (1-lambda1+lambda2-rho)/T(2);
- complex<T> arg1b = -(1-lambda1-lambda2-rho)/lambda1/T(2);
- complex<T> arg2b = -(1-lambda1-lambda2+rho)/lambda1/T(2);
+ complex<T> arg1a = (T(1)-lambda1+lambda2+rho)/T(2);
+ complex<T> arg2a = (T(1)-lambda1+lambda2-rho)/T(2);
+ complex<T> arg1b = -(T(1)-lambda1-lambda2-rho)/lambda1/T(2);
+ complex<T> arg2b = -(T(1)-lambda1-lambda2+rho)/lambda1/T(2);
complex<T> x1 = s*(arg1a-T(1))/m3sq, x2 = s*(arg2a-T(1))/m3sq;
complex<T> ln1 = ln(x1),
@@ -1828,7 +1828,7 @@ template <class T> complex<T> I4w4m(int order,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
@@ -2070,7 +2070,7 @@ template <class T> complex<T> I3w3m(int order, momentum_configuration<T>& k,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
diff --git a/src/integrals_ep.cpp b/src/integrals_ep.cpp
index f1dbd75..1973515 100644
--- a/src/integrals_ep.cpp
+++ b/src/integrals_ep.cpp
@@ -95,12 +95,12 @@ template <> inline RGMP epsilon<RGMP>(){ return exp10(-RGMP(RGMP::get_current_nb
// \theta(r); it is a macro rather than a function because we want to avoid
// evaluation of v unless r >= 0 (the evaluation might be invalid)
-#define Theta(r,v) ( (r) >= 0 ? (v) : 0 )
+#define Theta(r,v) ( (r) >= T(0) ? (v) : T(0) )
// \theta(r1)-\theta(r2)
#define Theta2(r1,r2,v) \
- ( ((r1) >= 0 and (r2) < 0) ? (v) \
- : ((r1) < 0 and (r2) >= 0) ? (-v) : 0 )
+ ( ((r1) >= T(0) and (r2) < T(0)) ? (v) \
+ : ((r1) < T(0) and (r2) >= T(0)) ? (-v) : T(0) )
// 3-mass triangle Gram det
template <class T> T Delta(T s1, T s2,T s3)
@@ -197,8 +197,8 @@ template <class T> complex<T> CLi2(T arg,T branchSide)
template <class T> inline T eta(T a, T b, T ab)
// args are Im a, Im b, and Im (a b) compared to eta in Denner, Nierste, & Scharf
// return value doesn't include 2 Pi I factor
-{return( a < 0 and b < 0 and ab >= 0 ? 1
- : (a >= 0 and b >= 0 and ab < 0) ? -1 : 0 );}
+{return( a < T(0) and b < T(0) and ab >= T(0) ? T(1)
+ : (a >= T(0) and b >= T(0) and ab < T(0)) ? T(-1) : T(0) );}
// 11/20/08
@@ -369,7 +369,7 @@ template <class T> complex<T> I4w4m(int order,
// But we must treat real & complex case for rho separately
// because of branch cut handling
- if (rhoSq >= 0)
+ if (rhoSq >= T(0))
{//T rho = sqrt(rhoSq);
complex<T> rho;
@@ -448,7 +448,7 @@ template <class T> complex<T> I4w4m(int order,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
@@ -481,11 +481,11 @@ template <class T> complex<T> I4w4m(int order,
T crit1, // coefficient of I eps in x_i
crit2; // coefficient of I eps' from m^2->m^2+I eps', in x_i;
- crit1 = 1/T(2)* m2sq/(m3sq* t* sqrt(-Delta3));
+ crit1 = T(0.5)* m2sq/(m3sq* t* sqrt(-Delta3));
- crit2 = 1/T(2)*(m1sq*square(m3sq) - m2sq*m3sq*m4sq + m2sq*m3sq*t - square(m3sq)*t - m2sq*m4sq*t +
+ crit2 = T(0.5)*(m1sq*square(m3sq) - m2sq*m3sq*m4sq + m2sq*m3sq*t - square(m3sq)*t - m2sq*m4sq*t +
m3sq*m4sq*t - m3sq*square(t) + s*square(t))/(square(m3sq)*square(t))
- +1/T(2)*(square(m1sq*m3sq)*m3sq - T(2)*m1sq*m2sq*square(m3sq)*m4sq + square(m2sq*m4sq)*m3sq +
+ +T(0.5)*(square(m1sq*m3sq)*m3sq - T(2)*m1sq*m2sq*square(m3sq)*m4sq + square(m2sq*m4sq)*m3sq +
m1sq*m2sq*square(m3sq)*t - m1sq*square(m3sq)*m3sq*t - m1sq*m2sq*m3sq*m4sq*t -
square(m2sq)*m3sq*m4sq*t + m1sq*square(m3sq)*m4sq*t + m2sq*square(m3sq)*m4sq*t +
square(m2sq*m4sq)*t - m2sq*m3sq*square(m4sq)*t - m1sq*square(m3sq)*s*t -
@@ -499,9 +499,9 @@ template <class T> complex<T> I4w4m(int order,
#endif
tag += " C";
- tag += Sign(crit1) > 0 ? "+" : "-";
- tag += Sign(crit2) > 0 ? "+" : "-";
- tag += Sign(crit1+crit2)/Sign(crit1) > 0 ? "v" : "x";
+ tag += Sign(crit1) > T(0) ? "+" : "-";
+ tag += Sign(crit2) > T(0) ? "+" : "-";
+ tag += Sign(crit1+crit2)/Sign(crit1) > T(0) ? "v" : "x";
tag += " ";
if (m1sq > 0) count += 1;
@@ -577,19 +577,19 @@ template <class T> complex<T> I4w4m(int order,
cout << "signs: " << Sign(Re(s1)) << " " << Sign(Re(s2)) << " " << Sign(Re(s3)) << endl;
#endif
tag += " ";
- tag += Re(s1) > 0 ? "+" : "-";
- tag += Re(s2) > 0 ? "+" : "-";
- tag += Re(s3) > 0 ? "+" : "-";
+ tag += Re(s1) > T(0) ? "+" : "-";
+ tag += Re(s2) > T(0) ? "+" : "-";
+ tag += Re(s3) > T(0) ? "+" : "-";
tag += " ";
T crit1, // coefficient of I eps in x_i
crit2; // coefficient of I eps' from m^2->m^2+I eps', in x_i;
- crit1 = 1/T(2)* m2sq/(m3sq* t* Re(sqrt(-Delta3)));
+ crit1 = T(0.5)* m2sq/(m3sq* t* Re(sqrt(-Delta3)));
- crit2 = 1/T(2)*(m1sq*square(m3sq) - m2sq*m3sq*m4sq + m2sq*m3sq*t - square(m3sq)*t - m2sq*m4sq*t +
+ crit2 = T(0.5)*(m1sq*square(m3sq) - m2sq*m3sq*m4sq + m2sq*m3sq*t - square(m3sq)*t - m2sq*m4sq*t +
m3sq*m4sq*t - m3sq*square(t) + s*square(t))/(square(m3sq)*square(t))
- +1/T(2)*(square(m1sq*m3sq)*m3sq - T(2)*m1sq*m2sq*square(m3sq)*m4sq + square(m2sq*m4sq)*m3sq +
+ +T(0.5)*(square(m1sq*m3sq)*m3sq - T(2)*m1sq*m2sq*square(m3sq)*m4sq + square(m2sq*m4sq)*m3sq +
m1sq*m2sq*square(m3sq)*t - m1sq*square(m3sq)*m3sq*t - m1sq*m2sq*m3sq*m4sq*t -
square(m2sq)*m3sq*m4sq*t + m1sq*square(m3sq)*m4sq*t + m2sq*square(m3sq)*m4sq*t +
square(m2sq*m4sq)*t - m2sq*m3sq*square(m4sq)*t - m1sq*square(m3sq)*s*t -
@@ -603,9 +603,9 @@ template <class T> complex<T> I4w4m(int order,
#endif
tag += " C";
- tag += Sign(crit1) > 0 ? "+" : "-";
- tag += Sign(crit2) > 0 ? "+" : "-";
- tag += Sign(crit1+crit2)/Sign(crit1) > 0 ? "v" : "x";
+ tag += Sign(crit1) > T(0) ? "+" : "-";
+ tag += Sign(crit2) > T(0) ? "+" : "-";
+ tag += Sign(crit1+crit2)/Sign(crit1) > T(0) ? "v" : "x";
tag += " ";
if (m1sq > 0) count += 1;
@@ -621,7 +621,7 @@ template <class T> complex<T> I4w4m(int order,
for (int j = 0; j < 3; j += 1)
{complex<T> arg = (ImSqrtDelta+delta[j])/(-ImSqrtDelta+delta[j]);
// One term has the wrong sign of imaginary part; but why?
- tag += Im( ((Li2(arg)-Li2(T(1)/arg))/ImSqrtDelta) ) > 0 ? "+" : "-";
+ tag += Im( ((Li2(arg)-Li2(T(1)/arg))/ImSqrtDelta) ) > T(0) ? "+" : "-";
#if DumpTerms
cout << j << ": " << arg << "; "
<< ((Li2(arg))/ImSqrtDelta)
@@ -674,11 +674,11 @@ template <class T> complex<T> I4w4m(int order,
T crit1, // coefficient of I eps in x_i
crit2; // coefficient of I eps' from m^2->m^2+I eps', in x_i;
- crit1 = 1/T(2)* m2sq/(m3sq* t* Re(rho));
+ crit1 = T(0.5)* m2sq/(m3sq* t* Re(rho));
- crit2 = 1/T(2)*(m1sq*square(m3sq) - m2sq*m3sq*m4sq + m2sq*m3sq*t - square(m3sq)*t - m2sq*m4sq*t +
+ crit2 = T(0.5)*(m1sq*square(m3sq) - m2sq*m3sq*m4sq + m2sq*m3sq*t - square(m3sq)*t - m2sq*m4sq*t +
m3sq*m4sq*t - m3sq*square(t) + s*square(t))/(square(m3sq)*square(t))
- +1/T(2)*(square(m1sq*m3sq)*m3sq - T(2)*m1sq*m2sq*square(m3sq)*m4sq + square(m2sq*m4sq)*m3sq +
+ +T(0.5)*(square(m1sq*m3sq)*m3sq - T(2)*m1sq*m2sq*square(m3sq)*m4sq + square(m2sq*m4sq)*m3sq +
m1sq*m2sq*square(m3sq)*t - m1sq*square(m3sq)*m3sq*t - m1sq*m2sq*m3sq*m4sq*t -
square(m2sq)*m3sq*m4sq*t + m1sq*square(m3sq)*m4sq*t + m2sq*square(m3sq)*m4sq*t +
square(m2sq*m4sq)*t - m2sq*m3sq*square(m4sq)*t - m1sq*square(m3sq)*s*t -
@@ -692,9 +692,9 @@ template <class T> complex<T> I4w4m(int order,
#endif
tag += " C";
- tag += Sign(crit1) > 0 ? "+" : "-";
- tag += Sign(crit2) > 0 ? "+" : "-";
- tag += Sign(crit1+crit2)/Sign(crit1) > 0 ? "v" : "x";
+ tag += Sign(crit1) > T(0) ? "+" : "-";
+ tag += Sign(crit2) > T(0) ? "+" : "-";
+ tag += Sign(crit1+crit2)/Sign(crit1) > T(0) ? "v" : "x";
tag += " ";
#if DumpArguments
@@ -770,10 +770,10 @@ template <class T> complex<T> I4w4m(int order,
cout << "imaginary rho" << endl;
#endif
- complex<T> arg1a = (1-lambda1+lambda2+rho)/T(2);
- complex<T> arg2a = (1-lambda1+lambda2-rho)/T(2);
- complex<T> arg1b = -(1-lambda1-lambda2-rho)/lambda1/T(2);
- complex<T> arg2b = -(1-lambda1-lambda2+rho)/lambda1/T(2);
+ complex<T> arg1a = (T(1)-lambda1+lambda2+rho)/T(2);
+ complex<T> arg2a = (T(1)-lambda1+lambda2-rho)/T(2);
+ complex<T> arg1b = -(T(1)-lambda1-lambda2-rho)/lambda1/T(2);
+ complex<T> arg2b = -(T(1)-lambda1-lambda2+rho)/lambda1/T(2);
complex<T> x1 = s*(arg1a-T(1))/m3sq, x2 = s*(arg2a-T(1))/m3sq;
complex<T> ln1 = ln(x1),
@@ -859,7 +859,7 @@ template <class T> complex<T> I4w4m(int order,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
@@ -1050,7 +1050,7 @@ template <class T> complex<T> I3w1m(int order,
case -1:
return(T(1)/s*CLnM(s,musq));
case 0:
- return(-1/(T(2)*s)*square(CLnM(s,musq)));
+ return(-T(1)/(T(2)*s)*square(CLnM(s,musq)));
default:
if (order < -2) return complex<T>(0,0);
#if _USE_GCC
@@ -1095,7 +1095,7 @@ template <class T> complex<T> I3w3m(int order,
cout << "Delta3: " << Delta3 << endl;
cout << "signs: " << Sign(s1) << " " << Sign(s2) << " " << Sign(s3) << endl;
#endif
- if (Delta3 >= 0)
+ if (Delta3 >= T(0))
{T sqrtDelta = sqrt(Delta3);
/* In this case, the arguments of the dilogs are generic complex
numbers, and we don't need to do anything special */
diff --git a/src/mom_conf.cpp b/src/mom_conf.cpp
index 4f46cdd..13948d2 100644
--- a/src/mom_conf.cpp
+++ b/src/mom_conf.cpp
@@ -2,6 +2,7 @@
\brief Implementation of the classe mom_conf
*/
+#include "BH_typedefs.h"
#include "mom_conf.h"
#include "BH_error.h"
#include <cstdarg>
diff --git a/src/mom_conf.h b/src/mom_conf.h
index 8f30a83..6454ca4 100644
--- a/src/mom_conf.h
+++ b/src/mom_conf.h
@@ -33,7 +33,10 @@
#if _USE_GCC
#include <ext/hash_map>
-using namespace __gnu_cxx;
+using __gnu_cxx::hash;
+using __gnu_cxx::hash_map;
+using __gnu_cxx::vector;
+using __gnu_cxx::pair;
namespace __gnu_cxx {
template<> struct hash<std::string>{
diff --git a/src/polylog_HP.hpp b/src/polylog_HP.hpp
index 5275ba2..02f0d83 100644
--- a/src/polylog_HP.hpp
+++ b/src/polylog_HP.hpp
@@ -99,11 +99,11 @@ const RHP Bernoulli2_HP[]= {
x = (x-RHP("1.0"))/x;}
else if (x > 0.5)
// PolyLog[2,x] -> -PolyLog[2,1-x] + Pi^2/6 - Log[x] Log[1-x]
- {added = PiSquaredOver6_HP - log(x)*log(1-x);
+ {added = PiSquaredOver6_HP - log(x)*log(RHP(1)-x);
factor = RHP("-1.0");
x = RHP("1.0")-x;}
- else if (x > 0) {}
- else if (x >= -1)
+ else if (x > RHP(0)) {}
+ else if (x >= RHP(-1))
// PolyLog[2, y] -> -Log[1 - y]^2/2 - PolyLog[2, y/(-1 + y)]
{added = RHP("-0.5")*sq(log(RHP("1.0")-x));
factor = RHP("-1.0");
@@ -206,7 +206,7 @@ RHP("-2.8382249570693706959264156336481764738284680928012882128228531714e78")
assert (sizeof(Bernoulli2_HP)/sizeof(Bernoulli2_HP[0]) >= CLi2TermLimit_HP);
- if (Re(z*conj(z)) > 1)
+ if (Re(z*conj(z)) > RHP(1))
{added = -PiSquaredOver6_HP-sq(log(-z))/RHP(2);
factor = RHP("-1");
z = RHP("1.")/z;}
@@ -326,7 +326,7 @@ RHP("-4.8412600798208880508789196709963412761130549942324620385115856258e129")
RHP factor = RHP("1.0");
RHP result;
- if (x < 0) {x = -x; factor = RHP("-1.0");}
+ if (x < RHP(0)) {x = -x; factor = RHP("-1.0");}
while (x > TwoPi) x -= TwoPi;
// shift to range [0,2 Pi/3] (larger upper limit prevents too many recursions)
diff --git a/src/polylog_VHP.hpp b/src/polylog_VHP.hpp
index da768fb..9f41382 100644
--- a/src/polylog_VHP.hpp
+++ b/src/polylog_VHP.hpp
@@ -146,11 +146,11 @@ RVHP ReLi2(RVHP x)
x = (x-RVHP("1.0"))/x;}
else if (x > 0.5)
// PolyLog[2,x] -> -PolyLog[2,1-x] + Pi^2/6 - Log[x] Log[1-x]
- {added = PiSquaredOver6_VHP - log(x)*log(1-x);
+ {added = PiSquaredOver6_VHP - log(x)*log(RVHP(1)-x);
factor = RVHP("-1.0");
x = RVHP("1.0")-x;}
- else if (x > 0) {}
- else if (x >= -1)
+ else if (x > RVHP(0)) {}
+ else if (x >= RVHP(-1))
// PolyLog[2, y] -> -Log[1 - y]^2/2 - PolyLog[2, y/(-1 + y)]
{added = RVHP("-0.5")*sq(log(RVHP("1.0")-x));
factor = RVHP("-1.0");
@@ -251,7 +251,7 @@ RVHP("-2.8382249570693706959264156336481764738284680928012882128228531714e78")
assert (sizeof(Bernoulli2_VHP)/sizeof(Bernoulli2_VHP[0]) >= CLi2TermLimit_VHP);
- if (Re(z*conj(z)) > 1)
+ if (Re(z*conj(z)) > RVHP(1))
{added = -PiSquaredOver6_VHP-sq(log(-z))/RVHP(2);
factor = RVHP("-1");
z = RVHP("1.")/z;}
@@ -371,7 +371,7 @@ RVHP("-4.8412600798208880508789196709963412761130549942324620385115856258e129")
RVHP factor = RVHP("1.0");
RVHP result;
- if (x < 0) {x = -x; factor = RVHP("-1.0");}
+ if (x < RVHP(0)) {x = -x; factor = RVHP("-1.0");}
while (x > TwoPi) x -= TwoPi;
// shift to range [0,2 Pi/3] (larger upper limit prevents too many recursions)
diff --git a/src/qd_suppl.cpp b/src/qd_suppl.cpp
index 83e71f8..9619f19 100644
--- a/src/qd_suppl.cpp
+++ b/src/qd_suppl.cpp
@@ -1,3 +1,4 @@
+#include "BH_typedefs.h"
#include "qd_suppl.h"
#include "BH_typedefs.h"
namespace BH {
diff --git a/src/ratext/pentagon_ratext.hpp b/src/ratext/pentagon_ratext.hpp
index 23ecbaf..c77330a 100644
--- a/src/ratext/pentagon_ratext.hpp
+++ b/src/ratext/pentagon_ratext.hpp
@@ -389,7 +389,7 @@ template <class BASE, class RatPentSpecs> template <class T> void pentagon_Rat<B
complex<T> m0sol=gamma*(pow(tsol_at_m0sol,2)*v1K5-T(0.5)*a1*tsol_at_m0sol+alp1*alp2*v2K5)/v2K5;
// Always put the branch cut in the same place for the sqrt(mass) computation
- if(abs(imag(m0sol))<DeltaZero<T>()){
+ if(fabs(imag(m0sol)) < DeltaZero<T>()){
m0sol=complex<T>(real(m0sol),0);
}
diff --git a/src/ratext/ratext_part.hpp b/src/ratext/ratext_part.hpp
index b8c2a67..e369354 100644
--- a/src/ratext/ratext_part.hpp
+++ b/src/ratext/ratext_part.hpp
@@ -51,7 +51,7 @@ template <class PENT,class BOX,class TRI,class BUB> template <class T> complex<T
#if _CONSERVATIVE_ERROR_ESTIMATE==1
T treeres=abs(get_tree(mc,ind)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<T>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=T(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>T(1)){treeres=T(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(T(10)));
}
@@ -284,7 +284,7 @@ template <class PENT,class BOX,class TRI,class BUB> complex<RHP> ratext_part<PEN
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RHP treeres=abs(get_tree(mc,ind)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RHP(1)){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RHP(10.)));
}
@@ -401,7 +401,7 @@ template <class PENT,class BOX,class TRI,class BUB> template <class T> complex<T
#if _CONSERVATIVE_ERROR_ESTIMATE==1
T treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<T>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=T(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>T(1)){treeres=T(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(T(10)));
}
@@ -635,7 +635,7 @@ template <class PENT,class BOX,class TRI,class BUB> complex<RHP> ratext_part<PEN
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RHP treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RHP(1)){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RHP(10)));
}
diff --git a/src/rational_eval/R_2q4g_eval.cpp b/src/rational_eval/R_2q4g_eval.cpp
index d6beb11..febc24c 100644
--- a/src/rational_eval/R_2q4g_eval.cpp
+++ b/src/rational_eval/R_2q4g_eval.cpp
@@ -5,6 +5,7 @@
* Author: daniel
*/
+#include "BH_typedefs.h"
#include "R_2q4g_eval.h"
#include "eval_param.h"
diff --git a/src/rational_eval/R_4g_eval.cpp b/src/rational_eval/R_4g_eval.cpp
index 21db739..31616a3 100644
--- a/src/rational_eval/R_4g_eval.cpp
+++ b/src/rational_eval/R_4g_eval.cpp
@@ -1,4 +1,5 @@
+#include "BH_typedefs.h"
#include "R_4g_eval.h"
#include "eval_param.h"
diff --git a/src/rational_eval/R_4q_eval.cpp b/src/rational_eval/R_4q_eval.cpp
index 0e9eef2..2941521 100644
--- a/src/rational_eval/R_4q_eval.cpp
+++ b/src/rational_eval/R_4q_eval.cpp
@@ -1,4 +1,5 @@
+#include "BH_typedefs.h"
#include "R_4g_eval.h"
#include "eval_param.h"
diff --git a/src/rational_eval/R_galln_eval.cpp b/src/rational_eval/R_galln_eval.cpp
index 394cabf..2655ff7 100644
--- a/src/rational_eval/R_galln_eval.cpp
+++ b/src/rational_eval/R_galln_eval.cpp
@@ -4,6 +4,7 @@
// all-p and all-m amplitudes, all-but-one p and all-but-one m amplitudes
// from Bern, Dixon, Kosower, PRD 72:125003 (2005) - "the last of the Mohicans"
+#include "BH_typedefs.h"
#include "R_alln_eval.h"
using namespace std;
diff --git a/src/splib.cpp b/src/splib.cpp
index f4159c8..2b3eeaa 100644
--- a/src/splib.cpp
+++ b/src/splib.cpp
@@ -2,6 +2,7 @@
\brief Implemetation for SpLib
This file implements the four-momenta and the spinor products. The conventions are the same as in the Mathematica package S\@M (aka SP\@M) so thatn a direct numerical comparaison is straight forward.
*/
+#include "BH_typedefs.h"
#include <iostream>
#include <iomanip>
#include <complex>
@@ -152,7 +153,7 @@ template <class T> Cmom<T> Cmom<T>::operator/=(const complex<T>& c){
}
template <class T> Cmom<T> Cmom<T>::operator*=(const complex<T>& c){
- if (c==T(0.)) {
+ if (c==complex<T>(0., 0.)) {
p=momentum<complex<T> >(complex<T>(0.,0.),complex<T>(0.,0.),complex<T>(0.,0.),complex<T>(0.,0.));
La=lambda<T>(0,0);
Lat=lambdat<T>(0,0);
diff --git a/src/splitamp_loop.cpp b/src/splitamp_loop.cpp
index b192448..ebe1d89 100644
--- a/src/splitamp_loop.cpp
+++ b/src/splitamp_loop.cpp
@@ -16,6 +16,7 @@
// primitive versus color-ordered amplitudes
// 1/6/09 - bug fixed
+#include "BH_typedefs.h"
#include "factorization.h"
#include "particles.h"
#include <cmath>
diff --git a/src/splitamp_tree.cpp b/src/splitamp_tree.cpp
index 3d37fdd..26d3ba8 100644
--- a/src/splitamp_tree.cpp
+++ b/src/splitamp_tree.cpp
@@ -8,6 +8,7 @@
// 6/23/08 - update to conform with new particle ID
// 12/21/08 -- added gluinos
+#include "BH_typedefs.h"
#include "factorization.h"
#include "BH_utilities.h"
diff --git a/src/standard_cut_part.hpp b/src/standard_cut_part.hpp
index aef8a08..1e5c401 100644
--- a/src/standard_cut_part.hpp
+++ b/src/standard_cut_part.hpp
@@ -59,7 +59,7 @@ template <class BOX,class TRI,class BUB> template <class T> SeriesC<T> standard_
std::complex<T> Ctreeres=get_tree(mc,ind); // Get the absolute power of the tree
T treeres=abs(Ctreeres); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<T>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=T(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>T(1)){treeres=T(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(T(10)));
}
@@ -353,7 +353,7 @@ template <class BOX,class TRI,class BUB> SeriesC<RHP> standard_cut_part<BOX,TRI,
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RHP treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RHP(1)){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RHP(10.)));
}
@@ -465,7 +465,7 @@ template <class BOX,class TRI,class BUB> SeriesC<RVHP> standard_cut_part<BOX,TRI
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RVHP treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RVHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RVHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RVHP(1)){treeres=RVHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RVHP(10)));
}
@@ -893,7 +893,7 @@ template <class BOX,class TRI,class BUB> SeriesC<RHP> standard_cut_part<BOX,TRI,
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RHP treeres=abs(get_tree(mc,ind)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RHP(1)){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RHP(10)));
}
@@ -1071,7 +1071,7 @@ template <class BOX,class TRI,class BUB> SeriesC<RVHP> standard_cut_part<BOX,TRI
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RVHP treeres=abs(get_tree(mc,ind)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RVHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RVHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RVHP(1)){treeres=RVHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RVHP(10)));
}
diff --git a/src/standard_cut_part_ep.hpp b/src/standard_cut_part_ep.hpp
index 27b818d..2194ac2 100644
--- a/src/standard_cut_part_ep.hpp
+++ b/src/standard_cut_part_ep.hpp
@@ -64,7 +64,7 @@ template <class BOX,class TRI,class BUB> template <class T> SeriesC<T> standard_
#if _CONSERVATIVE_ERROR_ESTIMATE==1
T treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<T>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=T(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>T(1)){treeres=T(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(T(10)));
}
@@ -372,7 +372,7 @@ template <class BOX,class TRI,class BUB> SeriesC<RHP> standard_cut_part<BOX,TRI,
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RHP treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RHP(1)){treeres=RHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RHP(10)));
}
@@ -488,7 +488,7 @@ template <class BOX,class TRI,class BUB> SeriesC<RVHP> standard_cut_part<BOX,TRI
#if _CONSERVATIVE_ERROR_ESTIMATE==1
RVHP treeres=abs(get_tree(ep)); // Get the absolute power of the tree
if(treeres > BH::DeltaZero<RVHP>()){ // If the tree is zero then there is no error_shift
- if(treeres>1){treeres=RVHP(1)/treeres;} // we turn all exponents into positive powers
+ if(treeres>RVHP(1)){treeres=RVHP(1)/treeres;} // we turn all exponents into positive powers
acc_shift=to_double(log(treeres)/log(RVHP(10)));
}
diff --git a/src/tree1.cc b/src/tree1.cc
index 2c25d6b..c8cb8e4 100644
--- a/src/tree1.cc
+++ b/src/tree1.cc
@@ -16,6 +16,7 @@
#define BRENDSGIELE_H_
#define BERENDSGIELE_IMPL_H_
//#include "BlackHat.h"
+#include "BH_typedefs.h"
#include <vector>
#include <map>
#include <ctime>
@@ -1485,6 +1486,9 @@ template <class T> momentum<complex<T> > GenerateMomentum(const T& dummy) {
return momentum<complex<T> >(sqrt(-bare*bare),bare.X(),bare.Y(),bare.Z());
}
+template
+momentum<complex<double> > GenerateMomentum<double>(const double& dummy);
+
#if 0
// Helicity ignored here
ParticleID FlavoredQuarkID(int flavor) {
@@ -2044,8 +2045,8 @@ template<class T> complex<T>
#endif
const vector<ParticleID>& leg,
int start, int end /* indices into the vectors */,
- int offshellMass = defaultMass,
- const vector<int>& massValue = empty)
+ int offshellMass,
+ const vector<int>& massValue)
#if 0
const vector<int>& leptons = empty
/* of leptons coupling via photon to the current,
diff --git a/src/tree2.cc b/src/tree2.cc
index 76de4aa..ba3bb7c 100644
--- a/src/tree2.cc
+++ b/src/tree2.cc
@@ -7,6 +7,7 @@
*/
+#include "BH_typedefs.h"
#include <vector>
#include <map>
#include "spinor.h"
@@ -748,10 +749,10 @@ template<class T> complex<T>
const vector<int>& arg /* indices of arguments */,
const vector<particle_ID>& leg /* helicities and particle ids */,
// Vector bosons
- const vector<int>& vectorK /* momenta */ = empty,
- const vector<int>& polarization = empty,
- const vector<int>& coupleTo /* quark flavor */ = empty,
- const vector<int>& massValue = empty)
+ const vector<int>& vectorK /* momenta */,
+ const vector<int>& polarization,
+ const vector<int>& coupleTo /* quark flavor */,
+ const vector<int>& massValue)
{ // until the right routine is available:
vector<int> helicity = Helicities(leg);
vector<int> idcode = ParticleCode(leg);
diff --git a/src/tree3.cc b/src/tree3.cc
index cfb4605..b641903 100644
--- a/src/tree3.cc
+++ b/src/tree3.cc
@@ -10,6 +10,7 @@
#define BRENDSGIELE_H_
#define BERENDSGIELE_IMPL_H_ // Don't include standard files
+#include "BH_typedefs.h"
#include <vector>
#include <map>
#include "spinor.h"
@@ -881,11 +882,11 @@ template<class T> complex<T>
of helicity L or R -- for helicity +/- the direction is
deduced from the helicity */,
// Vector bosons
- const vector<int>& vectorK /* momenta */ = empty,
- const vector<int>& polarization = empty,
- const vector<int>& coupleTo /* quark flavor */ = empty,
- int offshellMass = defaultMass,
- const vector<int>& massValue = empty)
+ const vector<int>& vectorK /* momenta */,
+ const vector<int>& polarization,
+ const vector<int>& coupleTo /* quark flavor */,
+ int offshellMass,
+ const vector<int>& massValue)
{ // Get reference momentum
size_t ref /* index of reference momentum */;
const string refKey(RefTag);
diff --git a/src/treebase.cc b/src/treebase.cc
index 7a49823..8ea374b 100644
--- a/src/treebase.cc
+++ b/src/treebase.cc
@@ -8,6 +8,7 @@ Entry point for tree routines.
#define BRENDSGIELE_H_
#define BERENDSGIELE_IMPL_H_
+#include "BH_typedefs.h"
#include <vector>
#include "spinor.h"
#include "mom_conf.h"
@@ -31,10 +32,10 @@ template<class T> complex<T> A(momentum_configuration<T>& k,
const vector<particle_ID>& leg,
int algorithm,
// Electroweak vectors
- const vector<int>& vectorK /* momenta */ = empty,
- const vector<int>& polarization = empty,
- const vector<int>& coupleTo /* quark flavor */ = empty,
- const vector<int>& massValue = empty)
+ const vector<int>& vectorK /* momenta */,
+ const vector<int>& polarization,
+ const vector<int>& coupleTo /* quark flavor */,
+ const vector<int>& massValue)
{vector<int> arg(leg.size()); // "leg" determines the number of arguments
for (int i = 0; i < arg.size(); i += 1) arg[i] = arg0[i];
// cout << "A: " << leg.size() << ": ";
@@ -59,7 +60,7 @@ template<class T> complex<T> A(momentum_configuration<T>& k,
// const vector<int>& vectorK /* momenta */ = empty,
// const vector<int>& polarization = empty,
// const vector<int>& coupleTo /* quark flavor */ = empty,
- const vector<int>& massValue = empty)
+ const vector<int>& massValue)
{vector<particle_ID> leg;
vector<int> vectorK(0), polarization(0), coupleTo(0);
// Assume leptons are paired sequentially, all couple to 1st quark flavor
diff --git a/src/util.cc b/src/util.cc
index 8af6b7c..6710c7b 100644
--- a/src/util.cc
+++ b/src/util.cc
@@ -9,6 +9,7 @@
*/
+#include "BH_typedefs.h"
#include <vector>
#include "spinor.h"
#include "mom_conf.h"
@@ -197,7 +198,7 @@ map<int,int> MassIndexCount(const vector<int>& massValue, int start, int end)
template<class T> int
MomentumSum(momentum_configuration<T>& k,
const vector<int>& v, int start, int end,
- const vector<int>& extraK = empty)
+ const vector<int>& extraK)
{/* For massless complex momenta, we must be careful to preserve the
separate lambda and lambda-tilde spinors, which will in general not
be identical to those produced by the "insert" below, even if k[sum]
@@ -236,7 +237,7 @@ template<class T> int
template<class T> /* static */ int FlatSum(momentum_configuration<T>& k,
int ref /* index of reference momentum */,
const vector<int>& v, int start, int end,
- const vector<int>& extraK = empty)
+ const vector<int>& extraK)
{// cout << "In FS " << start << ", " << end << endl;
/* For massless complex momenta, we must be careful to preserve the
separate lambda and lambda-tilde spinors, which will in general not
@@ -318,7 +319,7 @@ template<class T> int
NegativeFlatSum(momentum_configuration<T>& k,
int ref /* index of reference momentum */,
const vector<int>& v, int s1, int e1,
- int ev, const vector<int>& extraK = empty)
+ int ev, const vector<int>& extraK)
// With added vector momentum index "ev"
{// cout << "In NFS " << s1 << ", " << e1 << "; " << s2 << ", " << e2 << endl;
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment