HepLib
Loading...
Searching...
No Matches
Functions.cpp
Go to the documentation of this file.
1
6#include "SD.h"
7#include <cmath>
8
9namespace 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::VEO_Digits
int VEO_Digits
Definition Init.cpp:180
HepLib::SD::do_not_evalf_params
do_not_evalf_params().expl_derivative_func(zd2D).derivative_func(zp2D)) REGISTER_FUNCTION(CT
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:2293
HepLib::REGISTER_FUNCTION
REGISTER_FUNCTION(GMat, do_not_evalf_params().print_func< FCFormat >(&GMat_fc_print).conjugate_func(mat_conj).set_return_type(return_types::commutative)) bool IsZero(const ex &e)
Definition Basic.cpp:1002
Generated by doxygen 1.9.8