HepLib
Functions.cpp
Go to the documentation of this file.
1 
6 #include "SD.h"
7 #include <cmath>
8 
9 namespace HepLib::SD {
10 
11  namespace {
12  /*-----------------------------------------------------*/
13  // Functions used in GiNaC
14  /*-----------------------------------------------------*/
15  static void print_VEO(const ex & ex1_in, const ex & ex2_in, const print_context & c) {
16  static ex log10 = log(numeric("10"));
17  if(is_zero(ex1_in) && is_zero(ex2_in)) {
18  c.s << "0(0)";
19  } else if(is_zero(ex1_in)) {
20  numeric ex2 = ex_to<numeric>(Rationalize(ex2_in));
21  int n2 = floor(ex_to<numeric>(log(ex2)/log10).to_double());
22  double d2 = round(ex_to<numeric>(ex2/power(10,n2)).to_double());
23  if(d2==10) {d2=1; n2++;}
24  c.s << "0(" << d2 << ")";
25  if(n2!=0) c.s << "E" << n2;
26  } else if(is_zero(ex2_in)) {
27  numeric ex1 = ex_to<numeric>(Rationalize(ex1_in));
28  bool neg = bool(ex1<0);
29  if(neg) ex1 = 0-ex1;
30  int n1 = floor(ex_to<numeric>(log(ex1)/log10).to_double());
31  int n12 = VEO_Digits;
32  ostringstream oss;
33  oss << NN(ex1/power(10,n1));
34  string sd1 = oss.str();
35  if(sd1.length()<2+n12) {
36  oss.clear();
37  oss.str("");
38  oss << NN(ex1/power(10,n1)+power(10,-4-n12));
39  sd1 = oss.str();
40  if(sd1.length()<2+n12) throw Error("VEO: sd1.length()<2+n12");
41  } else {
42  int ci = stoi(sd1.substr(2+n12,1));
43  if(ci>4) {
44  oss.clear();
45  oss.str("");
46  oss << NN(ex1/power(10,n1)+power(10,-n12));
47  sd1 = oss.str();
48  }
49  }
50  sd1 = sd1.substr(0, 2+n12);
51 
52  if(neg) c.s << "(-";
53  c.s << sd1 << "(0)";
54  if(n1!=0) c.s << "E" << n1;
55  if(neg) c.s << ")";
56  } else {
57  numeric ex1 = ex_to<numeric>(Rationalize(ex1_in));
58  numeric ex2 = ex_to<numeric>(Rationalize(ex2_in));
59  bool neg = bool(ex1<0);
60  if(neg) ex1 = 0-ex1;
61 
62  int n1 = floor(ex_to<numeric>(log(ex1)/log10).to_double());
63  int n2 = floor(ex_to<numeric>(log(ex2)/log10).to_double());
64  int d2 = round(ex_to<numeric>(ex2/power(10,n2)).to_double());
65  if(d2==10) {d2=1; n2++;}
66 
67  if(n2>n1) {
68  if(n2==n1+1 && (ex1>5*power(10,n1))) c.s << "1";
69  else c.s << "0";
70  c.s << "(" << d2 << ")";
71  if(n2!=0) c.s << "E" << n2;
72  } else {
73  int n12 = n1-n2;
74  if(n12>VEO_Digits) {
75  n12 = VEO_Digits;
76  d2 = 0;
77  }
78  ostringstream oss;
79  oss << NN(ex1/power(10,n1));
80  string sd1 = oss.str();
81  if(sd1.length()<3+n12) {
82  oss.clear();
83  oss.str("");
84  oss << NN(ex1/power(10,n1)+power(10,-5-n12));
85  sd1 = oss.str();
86  if(sd1.length()<3+n12) throw Error("VEO2: sd1.length()<3+n12");
87  } else {
88  int ci = stoi(sd1.substr(2+n12,1));
89  if(ci>4) {
90  oss.clear();
91  oss.str("");
92  oss << NN(ex1/power(10,n1)+power(10,-n12));
93  sd1 = oss.str();
94  }
95  }
96  sd1 = sd1.substr(0, 2+n12);
97  if(sd1[sd1.length()-1]=='.') sd1 = sd1.substr(0, sd1.length()-1);
98  if(neg) c.s << "(-";
99  c.s << sd1 << "(" << d2 << ")";
100  if(n1!=0) c.s << "E" << n1;
101  if(neg) c.s << ")";
102  }
103  }
104  return;
105  }
106 
107  static void print_VEO2(const ex & ex1_in, const ex & ex2_in, const print_context & c) {
108  static ex log10 = log(numeric("10"));
109  if(is_zero(ex1_in) && is_zero(ex2_in)) {
110  c.s << "00(00)";
111  } else if(is_zero(ex1_in)) {
112  numeric ex2 = ex_to<numeric>(Rationalize(ex2_in));
113  int n2 = floor(ex_to<numeric>(log(ex2)/log10).to_double())-1;
114  double d2 = round(ex_to<numeric>(ex2/power(10,n2)).to_double());
115  if(d2==100) {d2=10; n2++;}
116  c.s << "00(" << d2 << ")";
117  if(n2!=0) c.s << "E" << n2;
118  } else if(is_zero(ex2_in)) {
119  numeric ex1 = ex_to<numeric>(Rationalize(ex1_in));
120  bool neg = bool(ex1<0);
121  if(neg) ex1 = 0-ex1;
122  int n1 = floor(ex_to<numeric>(log(ex1)/log10).to_double());
123  int n12 = VEO_Digits;
124  ostringstream oss;
125  oss << NN(ex1/power(10,n1));
126  string sd1 = oss.str();
127  if(sd1.length()<2+n12) {
128  oss.clear();
129  oss.str("");
130  oss << NN(ex1/power(10,n1)+power(10,-4-n12));
131  sd1 = oss.str();
132  if(sd1.length()<2+n12) throw Error("VEO2: sd1.length()<2+n12");
133  } else {
134  int ci = stoi(sd1.substr(2+n12,1));
135  if(ci>4) {
136  oss.clear();
137  oss.str("");
138  oss << NN(ex1/power(10,n1)+power(10,-n12));
139  sd1 = oss.str();
140  }
141  }
142  sd1 = sd1.substr(0, 2+n12);
143 
144  if(neg) c.s << "(-";
145  c.s << sd1 << "(0)";
146  if(n1!=0) c.s << "E" << n1;
147  if(neg) c.s << ")";
148  } else {
149  numeric ex1 = ex_to<numeric>(Rationalize(ex1_in));
150  numeric ex2 = ex_to<numeric>(Rationalize(ex2_in));
151  bool neg = bool(ex1<0);
152  if(neg) ex1 = 0-ex1;
153 
154  int n1 = floor(ex_to<numeric>(log(ex1)/log10).to_double());
155  int n2 = floor(ex_to<numeric>(log(ex2)/log10).to_double());
156  int d2 = round(ex_to<numeric>(ex2/power(10,n2-1)).to_double());
157  if(d2==100) {d2=10; n2++;}
158 
159  if(n2>=n1+2) {
160  if(n2==n1+2 && (ex1>5*power(10,n1-1))) c.s << "01";
161  else c.s << "00";
162  c.s << "(" << d2 << ")";
163  if(n2!=0) c.s << "E" << n2;
164  } else if(n2>n1) {
165  int d1 = round(ex_to<numeric>(ex1/power(10,n2-1)).to_double());
166  if(d1<10) c.s << "0";
167  c.s << d1;
168  c.s << "(" << d2 << ")";
169  if(n2!=0) c.s << "E" << n2;
170  } else {
171  int n12 = n1-n2;
172  if(n12>VEO_Digits) {
173  n12 = VEO_Digits;
174  d2 = 0;
175  }
176  ostringstream oss;
177  oss << NN(ex1/power(10,n1));
178  string sd1 = oss.str();
179  if(sd1.length()<4+n12) {
180  oss.clear();
181  oss.str("");
182  oss << NN(ex1/power(10,n1)+power(10,-6-n12));
183  sd1 = oss.str();
184  if(sd1.length()<4+n12) throw Error("VEO2: sd1.length()<4+n12");
185  } else {
186  int ci = stoi(sd1.substr(3+n12,1));
187  if(ci>4) {
188  oss.clear();
189  oss.str("");
190  oss << NN(ex1/power(10,n1)+power(10,-1-n12));
191  sd1 = oss.str();
192  }
193  }
194  sd1 = sd1.substr(0, 3+n12);
195  if(sd1[sd1.length()-1]=='.') sd1 = sd1.substr(0, sd1.length()-1);
196  if(neg) c.s << "(-";
197  if(sd1.length()==3 && sd1[1]=='.') {
198  c.s << "0." << sd1[0] << sd1[2] << "(" << d2 << ")";
199  if(n1!=-1) c.s << "E" << n1+1;
200  } else {
201  c.s << sd1 << "(" << d2 << ")";
202  if(n1!=0) c.s << "E" << n1;
203  }
204  if(neg) c.s << ")";
205  }
206  }
207  return;
208  }
209 
210  }
211 
212  namespace {
213  static ex conjVE(const ex & x, const ex & y) { return VE(x,y).hold(); }
214  static ex zp1D(const ex & x, unsigned diff_param) {return 0;}
215  static ex zd1D(const ex & x, const symbol & s) {return 0;}
216  static ex zp2D(const ex & x, const ex & y, unsigned diff_param) {return 0;}
217  static ex zd2D(const ex & x, const ex & y, const symbol & s) {return 0;}
218  static ex dCT(const ex & x, const symbol & s) {return x.diff(s);}
219  static ex pCT(const ex & x, unsigned diff_param) {return 1;}
220  }
221 
222  REGISTER_FUNCTION(CV, do_not_evalf_params()) // for use symbol
223  REGISTER_FUNCTION(fabs, dummy())
224  REGISTER_FUNCTION(PL, do_not_evalf_params().expl_derivative_func(zd1D).derivative_func(zp1D))
225  REGISTER_FUNCTION(FTX, do_not_evalf_params().expl_derivative_func(zd2D).derivative_func(zp2D))
226  REGISTER_FUNCTION(CT, do_not_evalf_params().expl_derivative_func(dCT).derivative_func(pCT))
227  REGISTER_FUNCTION(WRA, do_not_evalf_params().expl_derivative_func(zd1D).derivative_func(zp1D))
228  REGISTER_FUNCTION(VE, conjugate_func(conjVE))
229  REGISTER_FUNCTION(VEO, print_func<print_dflt>(print_VEO))
230  REGISTER_FUNCTION(VEO2, print_func<print_dflt>(print_VEO2))
231 }
SD.h
SecDec header file.
HepLib::SD
namespace for Numerical integration with Sector Decomposition method
Definition: AsyMB.cpp:10
HepLib::SD::REGISTER_FUNCTION
REGISTER_FUNCTION(CV, do_not_evalf_params()) REGISTER_FUNCTION(PL
HepLib::SD::do_not_evalf_params
do_not_evalf_params().expl_derivative_func(zd1D).derivative_func(zp1D)) REGISTER_FUNCTION(FTX
HepLib::NN
ex NN(ex expr, int digits)
the nuerical evaluation
Definition: BASIC.cpp:1278
HepLib::Rationalize
ex Rationalize(const ex &e, int dn)
Definition: BASIC.cpp:2283
Generated by doxygen 1.9.1