1 /*
2 * R : A Computer Language for Statistical Data Analysis
3 * Copyright (C) 2000-2016 The R Core Team.
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or
8 * (at your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful,
11 * but WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 * GNU General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, a copy is available at
17 * https://www.R-project.org/Licenses/
18 */
19
20 /* for mantelhaen.test */
21
22 #include <R.h>
23 #include <Rmath.h>
24
25 static void
int_d2x2xk(int K,double * m,double * n,double * t,double * d)26 int_d2x2xk(int K, double *m, double *n, double *t, double *d)
27 {
28 int i, j, l, w, y, z;
29 double u, **c;
30
31 c = (double **) R_alloc(K + 1, sizeof(double *));
32 l = y = z = 0;
33 c[0] = (double *) R_alloc(1, sizeof(double));
34 c[0][0] = 1;
35 for(i = 0; i < K; i++) {
36 y = imax2(0, (int)(*t - *n));
37 z = imin2((int)*m, (int)*t);
38 c[i + 1] = (double *) R_alloc(l + z - y + 1, sizeof(double));
39 for(j = 0; j <= l + z - y; j++) c[i + 1][j] = 0;
40 for(j = 0; j <= z - y; j++) {
41 u = dhyper(j + y, *m, *n, *t, FALSE);
42 for(w = 0; w <= l; w++) c[i + 1][w + j] += c[i][w] * u;
43 }
44 l = l + z - y;
45 m++; n++; t++;
46 }
47
48 u = 0;
49 for(j = 0; j <= l; j++) u += c[K][j];
50 for(j = 0; j <= l; j++) d[j] = c[K][j] / u;
51 }
52
53 #include <Rinternals.h>
54
d2x2xk(SEXP sK,SEXP m,SEXP n,SEXP t,SEXP srn)55 SEXP d2x2xk(SEXP sK, SEXP m, SEXP n, SEXP t, SEXP srn)
56 {
57 int K = asInteger(sK), rn = asInteger(srn);
58 m = PROTECT(coerceVector(m, REALSXP));
59 n = PROTECT(coerceVector(n, REALSXP));
60 t = PROTECT(coerceVector(t, REALSXP));
61 SEXP ans = PROTECT(allocVector(REALSXP, rn));
62 int_d2x2xk(K, REAL(m), REAL(n), REAL(t), REAL(ans));
63 UNPROTECT(4);
64 return ans;
65 }
66