ROOT  6.06/09
Reference Guide
testSpecFuncSiCi.cxx
Go to the documentation of this file.
1 #include "Math/SpecFunc.h"
2 #include <iostream>
3 #include <cmath>
4 
5 using std::cout;
6 using std::endl;
7 
8 // input for wolframalpha.com: Table[{N[x/10], N[CosIntegral[x/10], 18]}, {x, -100, 100}],
9 // delete imaginary part for x<0, remove entry for x=0
10 double civalues[200][2]=
11 {{-10.,-0.045456433004455373},
12 {-9.9,-0.036763956296836376},
13 {-9.8,-0.027519181109809388},
14 {-9.7,-0.017804097705837473},
15 {-9.6,-0.007707036058533522},
16 {-9.5,0.002678058835650657},
17 {-9.4,0.013252418669884902},
18 {-9.3,0.02391330446927572},
19 {-9.2,0.03455491341975988},
20 {-9.1,0.045069332542612144},
21 {-9.,0.055347531333133607},
22 {-8.9,0.06528038500431541},
23 {-8.8,0.074759719566173204},
24 {-8.7,0.083679369633444155},
25 {-8.6,0.091936239592257399},
26 {-8.5,0.099431358573421916},
27 {-8.4,0.106070919578643913},
28 {-8.3,0.11176729308811125},
29 {-8.2,0.116440005544566539},
30 {-8.1,0.120016673260596569},
31 {-8.,0.122433882532009557},
32 {-7.9,0.123638007059717843},
33 {-7.8,0.123585954183608557},
34 {-7.7,0.122245831911846316},
35 {-7.6,0.119597529284565836},
36 {-7.5,0.11563320323793427},
37 {-7.4,0.110357665828378243},
38 {-7.3,0.103788666432027628},
39 {-7.2,0.095957064345180793},
40 {-7.1,0.08690688807134784},
41 {-7.,0.076695278482184518},
42 {-6.9,0.065392313975951504},
43 {-6.8,0.053080716720199198},
44 {-6.7,0.039855440047043456},
45 {-6.6,0.025823138061263264},
46 {-6.5,0.011101519514930109},
47 {-6.4,-0.004181411011335063},
48 {-6.3,-0.019888220609842217},
49 {-6.2,-0.035873019273454991},
50 {-6.1,-0.051982528980021969},
51 {-6.,-0.068057243893247126},
52 {-5.9,-0.083932674118556494},
53 {-5.8,-0.099440664689378585},
54 {-5.7,-0.114410780761679062},
55 {-5.6,-0.128671749369807815},
56 {-5.5,-0.142052947551519255},
57 {-5.4,-0.154385926190724433},
58 {-5.3,-0.165505958558927253},
59 {-5.2,-0.175253602265659485},
60 {-5.1,-0.183476263159298894},
61 {-5.,-0.190029749656643879},
62 {-4.9,-0.194779806026237225},
63 {-4.8,-0.197603613309935236},
64 {-4.7,-0.198391246842472852},
65 {-4.6,-0.197047079722356195},
66 {-4.5,-0.193491122101738757},
67 {-4.4,-0.187660286800440685},
68 {-4.3,-0.179509572512633253},
69 {-4.2,-0.169013156767156739},
70 {-4.1,-0.15616539182812111},
71 {-4.,-0.140981697886930412},
72 {-3.9,-0.123499349207815143},
73 {-3.8,-0.103778150356897706},
74 {-3.7,-0.08190100128429845},
75 {-3.6,-0.057974351859800879},
76 {-3.5,-0.032128548512481116},
77 {-3.4,-0.004518077930741954},
78 {-3.3,0.024678284607958114},
79 {-3.2,0.055257411719942492},
80 {-3.1,0.086991831195536989},
81 {-3.,0.119629786008000328},
82 {-2.9,0.152895324159588311},
83 {-2.8,0.186488389643175768},
84 {-2.7,0.220084878632961618},
85 {-2.6,0.253336616062584192},
86 {-2.5,0.285871196365383495},
87 {-2.4,0.317291617436697984},
88 {-2.3,0.347175617540316224},
89 {-2.2,0.375074599049832154},
90 {-2.1,0.40051198784439639},
91 {-2.,0.422980828774864996},
92 {-1.9,0.441940349681598846},
93 {-1.8,0.456811129418336893},
94 {-1.7,0.46696836417695463},
95 {-1.6,0.471732516931877803},
96 {-1.5,0.470356317195399887},
97 {-1.4,0.462006585094677276},
98 {-1.3,0.445738567528534523},
99 {-1.2,0.420459182894240503},
100 {-1.1,0.384873377424650815},
101 {-1.,0.337403922900968135},
102 {-0.9,0.27606783046777286},
103 {-0.8,0.198278615952467177},
104 {-0.7,0.100514707008897833},
105 {-0.6,-0.022270706959279763},
106 {-0.5,-0.177784078806612901},
107 {-0.4,-0.378809346425244332},
108 {-0.3,-0.649172932971161745},
109 {-0.2,-1.042205595672781975},
110 {-0.1,-1.727868386657296639},
111 {0.1,-1.727868386657296639},
112 {0.2,-1.042205595672781975},
113 {0.3,-0.649172932971161745},
114 {0.4,-0.3788093464252443321},
115 {0.5,-0.1777840788066129013},
116 {0.6,-0.02227070695927976253},
117 {0.7,0.1005147070088978327},
118 {0.8,0.198278615952467177},
119 {0.9,0.2760678304677728602},
120 {1.,0.3374039229009681347},
121 {1.1,0.3848733774246508155},
122 {1.2,0.4204591828942405027},
123 {1.3,0.445738567528534523},
124 {1.4,0.4620065850946772763},
125 {1.5,0.4703563171953998867},
126 {1.6,0.4717325169318778034},
127 {1.7,0.4669683641769546303},
128 {1.8,0.4568111294183368931},
129 {1.9,0.4419403496815988459},
130 {2.,0.4229808287748649957},
131 {2.1,0.4005119878443963905},
132 {2.2,0.375074599049832154},
133 {2.3,0.3471756175403162244},
134 {2.4,0.3172916174366979837},
135 {2.5,0.2858711963653834954},
136 {2.6,0.2533366160625841922},
137 {2.7,0.2200848786329616183},
138 {2.8,0.1864883896431757677},
139 {2.9,0.1528953241595883109},
140 {3.,0.1196297860080003276},
141 {3.1,0.08699183119553698867},
142 {3.2,0.05525741171994249172},
143 {3.3,0.02467828460795811365},
144 {3.4,-0.004518077930741953583},
145 {3.5,-0.03212854851248111562},
146 {3.6,-0.05797435185980087899},
147 {3.7,-0.08190100128429845044},
148 {3.8,-0.103778150356897706},
149 {3.9,-0.1234993492078151427},
150 {4.,-0.1409816978869304116},
151 {4.1,-0.1561653918281211096},
152 {4.2,-0.1690131567671567387},
153 {4.3,-0.1795095725126332535},
154 {4.4,-0.1876602868004406848},
155 {4.5,-0.1934911221017387574},
156 {4.6,-0.1970470797223561954},
157 {4.7,-0.1983912468424728519},
158 {4.8,-0.1976036133099352364},
159 {4.9,-0.1947798060262372251},
160 {5.,-0.1900297496566438786},
161 {5.1,-0.1834762631592988937},
162 {5.2,-0.1752536022656594845},
163 {5.3,-0.1655059585589272527},
164 {5.4,-0.1543859261907244327},
165 {5.5,-0.1420529475515192553},
166 {5.6,-0.1286717493698078152},
167 {5.7,-0.1144107807616790621},
168 {5.8,-0.09944066468937858453},
169 {5.9,-0.08393267411855649376},
170 {6.,-0.0680572438932471262},
171 {6.1,-0.05198252898002196938},
172 {6.2,-0.03587301927345499141},
173 {6.3,-0.01988822060984221721},
174 {6.4,-0.004181411011335062901},
175 {6.5,0.01110151951493010868},
176 {6.6,0.0258231380612632636},
177 {6.7,0.03985544004704345647},
178 {6.8,0.05308071672019919828},
179 {6.9,0.06539231397595150445},
180 {7.,0.07669527848218451838},
181 {7.1,0.08690688807134783966},
182 {7.2,0.09595706434518079339},
183 {7.3,0.1037886664320276284},
184 {7.4,0.1103576658283782431},
185 {7.5,0.1156332032379342704},
186 {7.6,0.1195975292845658358},
187 {7.7,0.1222458319118463165},
188 {7.8,0.1235859541836085566},
189 {7.9,0.1236380070597178432},
190 {8.,0.1224338825320095573},
191 {8.1,0.1200166732605965689},
192 {8.2,0.1164400055445665387},
193 {8.3,0.1117672930881112503},
194 {8.4,0.1060709195786439133},
195 {8.5,0.09943135857342191604},
196 {8.6,0.09193623959225739866},
197 {8.7,0.08367936963344415541},
198 {8.8,0.07475971956617320393},
199 {8.9,0.06528038500431541014},
200 {9.,0.05534753133313360709},
201 {9.1,0.04506933254261214352},
202 {9.2,0.03455491341975987981},
203 {9.3,0.02391330446927571965},
204 {9.4,0.01325241866988490234},
205 {9.5,0.002678058835650657419},
206 {9.6,-0.007707036058533521531},
207 {9.7,-0.01780409770583747337},
208 {9.8,-0.02751918110980938821},
209 {9.9,-0.03676395629683637606},
210 {10.,-0.04545643300445537263}};
211 
212 // input for wolframalpha.com: Table[{N[x/10], N[SinIntegral[x/10], 18]}, {x, -100, 100}]
213 double sivalues[201][2]=
214 {{-10.,-1.658347594218874049},
215 {-9.9,-1.663384056595864698},
216 {-9.8,-1.667569616851386723},
217 {-9.7,-1.67084456972736341},
218 {-9.6,-1.673156980105444773},
219 {-9.5,-1.674463342281433091},
220 {-9.4,-1.674729172532594998},
221 {-9.3,-1.673929528316134798},
222 {-9.2,-1.672049447994015827},
223 {-9.1,-1.669084305598515127},
224 {-9.,-1.665040075829602495},
225 {-8.9,-1.659933505204107197},
226 {-8.8,-1.653792186051813822},
227 {-8.7,-1.646654530868780776},
228 {-8.6,-1.638569645386548639},
229 {-8.5,-1.629597099590385592},
230 {-8.4,-1.619806596812884436},
231 {-8.3,-1.609277541933421013},
232 {-8.2,-1.59809851062137326},
233 {-8.1,-1.586366622463643107},
234 {-8.,-1.574186821706942052},
235 {-7.9,-1.561671070214550172},
236 {-7.8,-1.548937458077995749},
237 {-7.7,-1.53610923812865961},
238 {-7.6,-1.523313791355258246},
239 {-7.5,-1.510681530943385878},
240 {-7.4,-1.498344753306055538},
241 {-7.3,-1.486436445063168036},
242 {-7.2,-1.475089055447246115},
243 {-7.1,-1.46443324405734104},
244 {-7.,-1.454596614248093591},
245 {-6.9,-1.445702442722501157},
246 {-6.8,-1.437868416091684784},
247 {-6.7,-1.43120538527026366},
248 {-6.6,-1.425816148589978626},
249 {-6.5,-1.421794274435881687},
250 {-6.4,-1.419222974038433272},
251 {-6.3,-1.418174034791727003},
252 {-6.2,-1.418706824114094167},
253 {-6.1,-1.420867373424620212},
254 {-6.,-1.424687551280506536},
255 {-5.9,-1.430184334109366196},
256 {-5.8,-1.437359182281820432},
257 {-5.7,-1.446197528508234603},
258 {-5.6,-1.456668384714831391},
259 {-5.5,-1.468724072665098669},
260 {-5.4,-1.482300082649290016},
261 {-5.3,-1.497315063575331138},
262 {-5.2,-1.513670946766480634},
263 {-5.1,-1.531253204712921712},
264 {-5.,-1.549931244944674137},
265 {-4.9,-1.569558938100652056},
266 {-4.8,-1.589975278172365569},
267 {-4.7,-1.611005171809781381},
268 {-4.6,-1.63246035250034989},
269 {-4.5,-1.654140414379243984},
270 {-4.4,-1.675833959408374161},
271 {-4.3,-1.697319850682468461},
272 {-4.2,-1.718368563690868599},
273 {-4.1,-1.738743626491768926},
274 {-4.,-1.758203138949053058},
275 {-3.9,-1.776501360447805439},
276 {-3.8,-1.793390354849570171},
277 {-3.7,-1.80862168087845377},
278 {-3.6,-1.821948115649503541},
279 {-3.5,-1.833125398665997048},
280 {-3.4,-1.841913983326143036},
281 {-3.3,-1.848080782795211419},
282 {-3.2,-1.851400897018440279},
283 {-3.1,-1.851659307674519836},
284 {-3.,-1.848652527999468256},
285 {-2.9,-1.842190194645858628},
286 {-2.8,-1.832096589081322327},
287 {-2.7,-1.818212076470866054},
288 {-2.6,-1.80039445052677016},
289 {-2.5,-1.778520173443826642},
290 {-2.4,-1.752485500761767386},
291 {-2.3,-1.722207481805503392},
292 {-2.2,-1.68762482724109852},
293 {-2.1,-1.64869863624441878},
294 {-2.,-1.605412976802694849},
295 {-1.9,-1.55777531374881852},
296 {-1.8,-1.505816780255578572},
297 {-1.7,-1.449592289683321133},
298 {-1.6,-1.389180485870438428},
299 {-1.5,-1.32468353117211968},
300 {-1.4,-1.256226732779217943},
301 {-1.3,-1.183958009076062919},
302 {-1.2,-1.10804719901371859},
303 {-1.1,-1.028685218673733522},
304 {-1.,-0.9460830703671830149},
305 {-0.9,-0.8604707107452929328},
306 {-0.8,-0.7720957854819965603},
307 {-0.7,-0.6812222391166113109},
308 {-0.6,-0.5881288096080800669},
309 {-0.5,-0.4931074180430666892},
310 {-0.4,-0.396461464751372883},
311 {-0.3,-0.2985040438070431614},
312 {-0.2,-0.1995560885262338214},
313 {-0.1,-0.09994446110827695016},
314 {0,0},
315 {0.1,0.09994446110827695016},
316 {0.2,0.1995560885262338214},
317 {0.3,0.2985040438070431614},
318 {0.4,0.396461464751372883},
319 {0.5,0.4931074180430666892},
320 {0.6,0.5881288096080800669},
321 {0.7,0.6812222391166113109},
322 {0.8,0.7720957854819965603},
323 {0.9,0.8604707107452929328},
324 {1.,0.9460830703671830149},
325 {1.1,1.028685218673733522},
326 {1.2,1.10804719901371859},
327 {1.3,1.183958009076062919},
328 {1.4,1.256226732779217943},
329 {1.5,1.32468353117211968},
330 {1.6,1.389180485870438428},
331 {1.7,1.449592289683321133},
332 {1.8,1.505816780255578572},
333 {1.9,1.55777531374881852},
334 {2.,1.605412976802694849},
335 {2.1,1.64869863624441878},
336 {2.2,1.68762482724109852},
337 {2.3,1.722207481805503392},
338 {2.4,1.752485500761767386},
339 {2.5,1.778520173443826642},
340 {2.6,1.80039445052677016},
341 {2.7,1.818212076470866054},
342 {2.8,1.832096589081322327},
343 {2.9,1.842190194645858628},
344 {3.,1.848652527999468256},
345 {3.1,1.851659307674519836},
346 {3.2,1.851400897018440279},
347 {3.3,1.848080782795211419},
348 {3.4,1.841913983326143036},
349 {3.5,1.833125398665997048},
350 {3.6,1.821948115649503541},
351 {3.7,1.80862168087845377},
352 {3.8,1.793390354849570171},
353 {3.9,1.776501360447805439},
354 {4.,1.758203138949053058},
355 {4.1,1.738743626491768926},
356 {4.2,1.718368563690868599},
357 {4.3,1.697319850682468461},
358 {4.4,1.675833959408374161},
359 {4.5,1.654140414379243984},
360 {4.6,1.63246035250034989},
361 {4.7,1.611005171809781381},
362 {4.8,1.589975278172365569},
363 {4.9,1.569558938100652056},
364 {5.,1.549931244944674137},
365 {5.1,1.531253204712921712},
366 {5.2,1.513670946766480634},
367 {5.3,1.497315063575331138},
368 {5.4,1.482300082649290016},
369 {5.5,1.468724072665098669},
370 {5.6,1.456668384714831391},
371 {5.7,1.446197528508234603},
372 {5.8,1.437359182281820432},
373 {5.9,1.430184334109366196},
374 {6.,1.424687551280506536},
375 {6.1,1.420867373424620212},
376 {6.2,1.418706824114094167},
377 {6.3,1.418174034791727003},
378 {6.4,1.419222974038433272},
379 {6.5,1.421794274435881687},
380 {6.6,1.425816148589978626},
381 {6.7,1.43120538527026366},
382 {6.8,1.437868416091684784},
383 {6.9,1.445702442722501157},
384 {7.,1.454596614248093591},
385 {7.1,1.46443324405734104},
386 {7.2,1.475089055447246115},
387 {7.3,1.486436445063168036},
388 {7.4,1.498344753306055538},
389 {7.5,1.510681530943385878},
390 {7.6,1.523313791355258246},
391 {7.7,1.53610923812865961},
392 {7.8,1.548937458077995749},
393 {7.9,1.561671070214550172},
394 {8.,1.574186821706942052},
395 {8.1,1.586366622463643107},
396 {8.2,1.59809851062137326},
397 {8.3,1.609277541933421013},
398 {8.4,1.619806596812884436},
399 {8.5,1.629597099590385592},
400 {8.6,1.638569645386548639},
401 {8.7,1.646654530868780776},
402 {8.8,1.653792186051813822},
403 {8.9,1.659933505204107197},
404 {9.,1.665040075829602495},
405 {9.1,1.669084305598515127},
406 {9.2,1.672049447994015827},
407 {9.3,1.673929528316134798},
408 {9.4,1.674729172532594998},
409 {9.5,1.674463342281433091},
410 {9.6,1.673156980105444773},
411 {9.7,1.67084456972736341},
412 {9.8,1.667569616851386723},
413 {9.9,1.663384056595864698},
414 {10.,1.658347594218874049}};
415 
416 int testSiCi() {
417  const double maxrelThreshold = 1E-13;
418  const double maxabsThreshold = 1.25E-15;
419 
420  int fail = 0;
421  double xmaxrel = -1, maxrel = 0, xmaxabs=-1, maxabs=0;
422  for (int i = 0; i < 201; ++i) {
423  double x = sivalues[i][0];
424  double SiMathematica = sivalues[i][1];
425  double SiRoot = ROOT::Math::sinint(x);
426  if (std::abs(SiRoot-SiMathematica) > maxabs) {
427  maxabs = std::fabs(SiRoot-SiMathematica);
428  xmaxabs = x;
429  }
430  if (std::abs(SiRoot/SiMathematica-1) > maxrel) {
431  maxrel = std::abs(SiRoot/SiMathematica-1);
432  xmaxrel = x;
433  }
434 
435 // std::cout << "x=" << x << ": Si(x) = " << SiMathematica
436 // << ", sinint(x) = " << SiRoot
437 // << ", rel. diff = " << SiRoot/SiMathematica-1
438 // << std::endl;
439  }
440  std::cout << "Si: Maximum relative deviation: " << maxrel
441  << " for x=" << xmaxrel
442  << (maxrel > maxabsThreshold ? " -> FAIL" : " -> pass") << std::endl;
443  std::cout << "Si: Maximum absolute deviation: " << maxabs
444  << " for x=" << xmaxabs
445  << (maxabs > maxabsThreshold ? " -> FAIL" : " -> pass") << std::endl;
446  if (maxrel > maxabsThreshold) fail += 1;
447  if (maxabs > maxabsThreshold) fail += 2;
448 
449  xmaxrel = -1; maxrel = 0; xmaxabs=-1; maxabs=0;
450  for (int i = 0; i < 200; ++i) {
451  double x = civalues[i][0];
452  double CiMathematica = civalues[i][1];
453  double CiRoot = ROOT::Math::cosint(x);
454  if (std::fabs(CiRoot-CiMathematica) > maxabs) {
455  maxabs = std::fabs(CiRoot-CiMathematica);
456  xmaxabs = x;
457  }
458  if (std::fabs(CiRoot/CiMathematica-1) > maxrel) {
459  maxrel = std::fabs(CiRoot/CiMathematica-1);
460  xmaxrel = x;
461  }
462 // std::cout << "x=" << x << ": Ci(x) = " << CiMathematica
463 // << ", cosint(x) = " << CiRoot
464 // << ", rel. diff = " << CiRoot/CiMathematica-1
465 // << std::endl;
466  }
467  std::cout << "Ci: Maximum relative deviation: " << maxrel
468  << " for x=" << xmaxrel
469  << (maxrel > maxrelThreshold ? " -> FAIL" : " -> pass") <<std::endl;
470  std::cout << "Ci: Maximum absolute deviation: " << maxabs
471  << " for x=" << xmaxabs
472  << (maxabs > maxabsThreshold ? " -> FAIL" : " -> pass")<< std::endl;
473 
474  if (maxrel > maxrelThreshold) fail += 4;
475  if (maxabs > maxabsThreshold) fail += 8;
476  return fail;
477 }
478 
480  int status = testSiCi();
481  std::cerr << "Test Special Function Sin integral and Cos integral :\t";
482  if (status) std::cout << "FAILED !!!" << std::endl;
483  else
484  std::cerr << "OK " << std::endl;
485  return status;
486 }
487 
488 int main() {
489  return testSpecFuncSiCi();
490 }
double cosint(double x)
Calculates the real part of the cosine integral (Ci).
double sinint(double x)
Calculates the sine integral.
TAlienJobStatus * status
Definition: TAlienJob.cxx:51
int main()
Double_t x[n]
Definition: legend1.C:17
int testSpecFuncSiCi()
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
double civalues[200][2]
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Double_t E()
Definition: TMath.h:54
int testSiCi()
double sivalues[201][2]