← Back to branch summary

~showard314/ubuntu/karmic/r-base/remove_start_comments

« back to all changes in this revision

Viewing changes to src/library/stats/src/swilk.c

  • Committer: Bazaar Package Importer
  • Author(s): Dirk Eddelbuettel
  • Date: 2009-01-19 12:40:24 UTC
  • mfrom: (5.1.4 sid)
  • Revision ID: james.westby@ubuntu.com-20090119124024-abxsf4e0y7713w9m
Tags: 2.8.1-2
http://bugs.debian.org/512324
debian/control: Add another Build-Depends: exclusion for the 
'kfreebsd-i386 kfreebsd-amd64 hurd-i386' architecture to openjdk-6-jdk.
Thanks to Petr Salinger for the heads-up.               (Closes: 512324)

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/library/stats/src/swilk.c
27
27
    const static float zero = 0.f;
28
28
    const static float one = 1.f;
29
29
    const static float two = 2.f;
30
 
    const static float three = 3.f;
31
30
 
32
 
    const static float z90 = 1.2816f;
33
 
    const static float z95 = 1.6449f;
34
 
    const static float z99 = 2.3263f;
35
 
    const static float zm = 1.7509f;
36
 
    const static float zss = .56268f;
37
 
    const static float bf1 = .8378f;
38
 
    const static double xx90 = .556;
39
 
    const static double xx95 = .622;
40
 
    const static float sqrth = .70711f;/* sqrt(1/2) = .7071068 */
41
31
    const static float small = 1e-19f;
42
 
    const static float pi6 = 1.909859f;
43
 
    const static float stqr = 1.047198f;
44
32
 
45
33
    /* polynomial coefficients */
46
34
    const static float g[2] = { -2.273f,.459f };
68
56
    float a1, a2, an, bf, ld, m, s, sa, xi, sx, xx, y, w1;
69
57
    float fac, asa, an25, ssa, z90f, sax, zfm, z95f, zsd, z99f, rsn, ssx, xsx;
70
58
 
71
 
    /* Parameter adjustments */
72
 
    --a;
73
 
 
74
59
    *pw = 1.;
75
60
    if (*w >= 0.) {
76
61
        *w = 1.;
77
62
    }
 
63
    if (*n < 3) {       *ifault = 1; return;
 
64
    }
 
65
 
78
66
    an = (float) (*n);
79
67
    nn2 = *n / 2;
80
 
    if (*n2 < nn2) {
81
 
        *ifault = 3; return;
82
 
    }
83
 
    if (*n < 3) {
84
 
        *ifault = 1; return;
85
 
    }
86
 
 
87
 
/*      If INIT is false, calculate coefficients a[] for the test */
 
68
    if (*n2 < nn2) {    *ifault = 3; return;
 
69
    }
 
70
    if (*n1 < 3) {      *ifault = 1; return;
 
71
    }
 
72
    ncens = *n - *n1;
 
73
    if (ncens < 0 || (ncens > 0 && *n < 20)) {  *ifault = 4; return;
 
74
    }
 
75
    if (ncens > 0) {
 
76
        delta = (float) ncens / an;
 
77
        if (delta > .8f) {      *ifault = 5; return;
 
78
        }
 
79
    } /* just for -Wall:*/ else { delta = 0.f; }
 
80
 
 
81
    --a; /* so we can keep using 1-based indices */
 
82
 
 
83
/*      If INIT is false (always when called from R),
 
84
 *      calculate coefficients a[] for the test statistic W */
88
85
    if (! (*init)) {
89
86
        if (*n == 3) {
 
87
            const static float sqrth = .70710678f;/* = sqrt(1/2), was .70711f */
90
88
            a[1] = sqrth;
91
89
        } else {
92
90
            an25 = an + .25;
119
117
        }
120
118
        *init = (1);
121
119
    }
122
 
    if (*n1 < 3) {
123
 
        *ifault = 1;    return;
124
 
    }
125
 
    ncens = *n - *n1;
126
 
    if (ncens < 0 || (ncens > 0 && *n < 20)) {
127
 
        *ifault = 4;    return;
128
 
    }
129
 
    delta = (float) ncens / an;
130
 
    if (delta > .8f) {
131
 
        *ifault = 5;    return;
132
 
    }
133
120
 
134
 
/*      If W input as negative, calculate significance level of -W */
 
121
/*      If W is input as negative, calculate significance level of -W */
135
122
 
136
123
    if (*w < zero) {
137
124
        w1 = 1. + *w;
148
135
 
149
136
/*      Check for correct sort order on range - scaled X */
150
137
 
151
 
    /* *ifault = 7; <-- a no-op, since it is set 0, below, in ANY CASE! */
 
138
    /* *ifault = 7; <-- a no-op, since it is changed below, in ANY CASE! */
152
139
    *ifault = 0;
153
140
    xx = x[0] / range;
154
141
    sx = xx;
158
145
        xi = x[i] / range;
159
146
        if (xx - xi > small) {
160
147
            /* Fortran had:      print *, "ANYTHING"
161
 
             * but do NOT; it *does* happen with sorted x (on Intel GNU/linux):
 
148
             * but do NOT; it *does* happen with sorted x (on Intel GNU/linux 32bit):
162
149
             *  shapiro.test(c(-1.7, -1,-1,-.73,-.61,-.5,-.24, .45,.62,.81,1))
163
150
             */
164
151
            *ifault = 7;
202
189
/*      Calculate significance level for W */
203
190
 
204
191
    if (*n == 3) {/* exact P value : */
 
192
        const static double pi6 = 1.90985931710274;/* = 6/pi, was  1.909859f */
 
193
        const static double stqr= 1.04719755119660;/* = asin(sqrt(3/4)), was 1.047198f */
205
194
        *pw = pi6 * (asin(sqrt(*w)) - stqr);
 
195
        if(*pw < 0.) *pw = 0.;
206
196
        return;
207
197
    }
208
198
    y = log(w1);
210
200
    if (*n <= 11) {
211
201
        gamma = poly(g, 2, an);
212
202
        if (y >= gamma) {
213
 
            *pw = small;/* FIXME: rather use an even smaller value, or NA ? */
 
203
            *pw = 1e-99;/* an "obvious" value, was 'small' which was 1e-19f */
214
204
            return;
215
205
        }
216
206
        y = -log(gamma - y);
222
212
    }
223
213
    /*DBG printf("c(w1=%g, w=%g, y=%g, m=%g, s=%g)\n",w1,*w,y,m,s); */
224
214
 
225
 
    if (ncens > 0) {/* <==>  n > n1 */
 
215
    if (ncens > 0) {/* <==>  n > n1 --- not happening currently when called from R */
226
216
 
227
217
/*      Censoring by proportion NCENS/N.
228
218
        Calculate mean and sd of normal equivalent deviate of W. */
229
219
 
 
220
        const static float three = 3.f;
 
221
 
 
222
        const static float z90 = 1.2816f;
 
223
        const static float z95 = 1.6449f;
 
224
        const static float z99 = 2.3263f;
 
225
        const static float zm = 1.7509f;
 
226
        const static float zss = .56268f;
 
227
        const static float bf1 = .8378f;
 
228
 
 
229
        const static double xx90 = .556;
 
230
        const static double xx95 = .622;
 
231
 
230
232
        ld = -log(delta);
231
233
        bf = one + xx * bf1;
232
234
        r__1 = pow(xx90, (double) xx);