4 {
5
6
7 static float srtopi = 2.0/sqrt(2.0*
M_PI);
8 static float upl = 100.0;
9
10 float prob = 0.0;
11 if(ndof <= 0) return prob;
12 if(chisq < 0.0) return prob;
13 if(ndof <= 60) {
14
15 if(chisq > upl) return prob;
16 float sum =
exp(-0.5*chisq);
17 float term = sum;
18
19 int m = ndof/2;
20 if(2*m == ndof) {
21 if(m == 1) return sum;
22 for(int i = 2; i <= m; i++) {
23 term = 0.5*term*chisq/(i-1);
24 sum += term;
25 }
26 return sum;
27
28 } else {
29
30 float srty = sqrt(chisq);
31 float temp = srty/M_SQRT2;
32 prob = erfc(temp);
33 if(ndof == 1) return prob;
34 if(ndof == 3) return (srtopi*srty*sum+prob);
35 m--;
36 for(int i = 1; i <= m; i++) {
37 term = term*chisq/(2*i+1);
38 sum += term;
39 }
40 return (srtopi*srty*sum + prob);
41
42 }
43
44 } else {
45
46 float srty = sqrt(chisq) - sqrt(ndof - 0.5);
47 if(srty < 12.0) prob=0.5*erfc(srty);
48 return prob;
49
50 }
51
52
53
54}
EvtComplex exp(const EvtComplex &c)