19 #include "local_proto.h" 23 static double Cdhc_alnfac(
int j);
24 static double Cdhc_correc(
int i,
int n);
28 void Cdhc_nscor1(
double s[],
int n,
int n2,
double work[],
int *ifault)
30 double ani, c, c1, d, scor;
47 d = c1 - log((
double)n);
50 for (i = 0; i < n2; ++i) {
51 ani = (double)n - i - 1;
53 for (scor = 0.0, j = 0; j <
NSTEP; ++j)
54 scor += work[0 * NSTEP + j] *
55 exp(work[1 * NSTEP + j] + work[2 * NSTEP + j] * i
56 + work[3 * NSTEP + j] * ani + c);
58 d += log((
double)(i + 1.0) / ani);
67 double xstart = -9.0, pi2 = -0.918938533, xx;
73 for (i = 0; i <
NSTEP; ++i) {
74 work[0 * NSTEP + i] = xx;
75 work[1 * NSTEP + i] = pi2 - xx * xx * 0.5;
78 xx = xstart +
H * (i + 1.0);
88 static double Cdhc_alnfac(
int j)
90 static double r[7] = { 0.0, 0.0, 0.69314718056, 1.79175946923,
91 3.17805383035, 4.78749174278, 6.57925121101
103 return (w - 0.5) * log(w) - w + 0.918938522305 +
104 (((4.0 - 3.0 * z) * z - 14.0) * z + 420.0) / (5040.0 * w);
113 static double eps[4] = { 0.419885, 0.450536, 0.456936, 0.468488 };
114 static double dl1[4] = { 0.112063, 0.121770, 0.239299, 0.215159 };
115 static double dl2[4] = { 0.080122, 0.111348, -0.211867, -0.115049 };
116 static double gam[4] = { 0.474798, 0.469051, 0.208597, 0.259784 };
117 static double lam[4] = { 0.282765, 0.304856, 0.407708, 0.414093 };
118 static double bb = -0.283833, d = -0.106136, b1 = 0.5641896;
139 k = (n2 < 3) ? n2 : 3;
140 for (i = 0; i < k; ++i) {
141 e1 = (1.0 + i - eps[i]) / (n + gam[i]);
142 e2 = pow(e1, lam[i]);
143 s[i] = e1 + e2 * (dl1[i] + e2 * dl2[i]) / n - Cdhc_correc(1 + i, n);
148 for (i = 3; i < n2; ++i) {
149 l1 = lam[3] + bb / (1.0 + i + d);
150 e1 = (1.0 + i - eps[3]) / (n + gam[3]);
152 s[i] = e1 + e2 * (dl1[3] + e2 * dl2[3]) / n - Cdhc_correc(1 + i, n);
157 for (i = 0; i < n2; ++i)
168 static double Cdhc_correc(
int i,
int n)
170 static double c1[7] = { 9.5, 28.7, 1.9, 0.0, -7.0, -6.2, -1.6 };
171 static double c2[7] = { -6.195e3, -9.569e3, -6.728e3, -17.614e3,
172 -8.278e3, -3.570e3, 1.075e3
174 static double c3[7] = { 9.338e4, 1.7516e5, 4.1040e5, 2.157e6,
175 2.376e6, 2.065e6, 2.065e6
177 static double mic = 1.0e-6, c14 = 1.9e-5;
185 else if (i != 4 && n > 20)
187 else if (i == 4 && n > 40)
191 an = 1.0 / (double)(n * n);
192 return (c1[i - 1] + an * (c2[i - 1] + an * c3[i - 1])) * mic;
void Cdhc_nscor1(double s[], int n, int n2, double work[], int *ifault)
double Cdhc_alnorm(double x, int upper)
void Cdhc_nscor2(double s[], int n, int n2, int *ifault)