4 static void f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
5 double *expect, double *percnt, double *emin, double
6 *prt, double *pre, double *fact, int *ico, int
7 *iro, int *kyy, int *idif, int *irn, int *key,
8 int *ldkey, int *ipoin, double *stp, int *ldstp,
9 int *ifrq, double *dlp, double *dsp, double *tm,
10 int *key2, int *iwk, double *rwk);
11 static void f3xact(int *nrow, int *irow, int *ncol, int *icol,
12 double *dlp, int *mm, double *fact, int *ico, int
13 *iro, int *it, int *lb, int *nr, int *nt, int
14 *nu, int *itc, int *ist, double *stv, double *alen,
16 static void f4xact(int *nrow, int *irow, int *ncol, int *icol,
17 double *dsp, double *fact, int *icstk, int *ncstk,
18 int *lstk, int *mstk, int *nstk, int *nrstk, int
19 *irstk, double *ystk, const double *tol);
20 static void f5xact(double *pastp, const double *tol, int *kval, int *key,
21 int *ldkey, int *ipoin, double *stp, int *ldstp,
22 int *ifrq, int *npoin, int *nr, int *nl, int
23 *ifreq, int *itop, int *ipsh);
24 static void f6xact(int *nrow, int *irow, int *iflag, int *kyy,
25 int *key, int *ldkey, int *last, int *ipn);
26 static void f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
28 static void f8xact(int *irow, int *is, int *i1, int *izero, int *myNew);
29 static double f9xact(int *n, int *mm, int *ir, double *fact);
30 static void f10act(int *nrow, int *irow, int *ncol, int *icol,
31 double *val, int *xmin, double *fact, int *nd,
33 static void f11act(int *irow, int *i1, int *i2, int *myNew);
34 static void prterr(int icode, char *mes);
35 static int iwork(int iwkmax, int *iwkpt, int number, int itype);
36 // void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
37 // double *expect, double *percnt, double *emin, double *prt,
38 // double *pre, /* myNew in C : */ int *workspace);
39 static void isort(int *n, int *ix);
40 static double gammds(double *y, double *p, int *ifault);
41 static double alogam(double *x, int *ifault);
46 /* The only public function : */
47 void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
48 double *expect, double *percnt, double *emin, double *prt,
49 double *pre, /* myNew in C : */ int *workspace) {
52 ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
53 THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
54 VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
55 -----------------------------------------------------------------------
57 Purpose: Computes Fisher's exact test probabilities and a hybrid
58 approximation to Fisher exact test probabilities for a
59 contingency table using the network algorithm.
60 Usage: CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
63 NROW - The number of rows in the table. (Input)
64 NCOL - The number of columns in the table. (Input)
65 TABLE - NROW by NCOL matrix containing the contingency
67 LDTABL - Leading dimension of TABLE exactly as specified
68 in the dimension statement in the calling
70 EXPECT - Expected value used in the hybrid algorithm for
71 deciding when to use asymptotic theory
72 probabilities. (Input)
73 If EXPECT <= 0.0 then asymptotic theory probabilities
74 are not used and Fisher exact test probabilities are
75 computed. Otherwise, if PERCNT or more of the cells in
76 the remaining table have estimated expected values of
77 EXPECT or more, with no remaining cell having expected
78 value less than EMIN, then asymptotic chi-squared
79 probabilities are used. See the algorithm section of the
80 manual document for details.
81 Use EXPECT = 5.0 to obtain the 'Cochran' condition.
82 PERCNT - Percentage of remaining cells that must have
83 estimated expected values greater than EXPECT
84 before asymptotic probabilities can be used. (Input)
85 See argument EXPECT for details.
86 Use PERCNT = 80.0 to obtain the 'Cochran' condition.
87 EMIN - Minimum cell estimated expected value allowed for
88 asymptotic chi-squared probabilities to be used. (Input)
89 See argument EXPECT for details.
90 Use EMIN = 1.0 to obtain the 'Cochran' condition.
91 PRT - Probability of the observed table for fixed
92 marginal totals. (Output)
93 PRE - Table p-value. (Output)
94 PRE is the probability of a more extreme table,
95 where `extreme' is in a probabilistic sense.
96 If EXPECT < 0 then the Fisher exact probability
97 is returned. Otherwise, an approximation to the
98 Fisher exact probability is computed based upon
99 asymptotic chi-squared probabilities for ``large''
100 table expected values. The user defines ``large''
101 through the arguments EXPECT, PERCNT, and EMIN.
104 1. For many problems one megabyte or more of workspace can be
105 required. If the environment supports it, the user should begin
106 by increasing the workspace used to 200,000 units.
107 2. In FEXACT, LDSTP = 30*LDKEY. The proportion of table space used
108 by STP may be changed by changing the line MULT = 30 below to
110 3. FEXACT may be converted to single precision by setting IREAL = 3,
111 and converting all DOUBLE PRECISION specifications (except the
112 specifications for RWRK, IWRK, and DWRK) to REAL. This will
113 require changing the names and specifications of the intrinsic
114 functions ALOG, AMAX1, AMIN1, EXP, and REAL. In addition, the
115 machine specific constants will need to be changed, and the name
116 DWRK will need to be changed to RWRK in the call to F2XACT.
117 4. Machine specific constants are specified and documented in F2XACT.
118 A missing value code is specified in both FEXACT and F2XACT.
119 5. Although not a restriction, is is not generally practical to call
120 this routine with large tables which are not sparse and in
121 which the 'hybrid' algorithm has little effect. For example,
122 although it is feasible to compute exact probabilities for the
127 computing exact probabilities for a similar table which has been
128 enlarged by the addition of an extra row (or column) may not be
130 -----------------------------------------------------------------------
133 /* CONSTANT Parameters : */
135 /* To increase the length of the table of paste path lengths relative
136 to the length of the hash table, increase MULT.
139 /* AMISS is a missing value indicator which is returned when the
140 probability is not defined.
142 const double amiss = -12345.;
144 Set IREAL = 4 for DOUBLE PRECISION
145 Set IREAL = 3 for SINGLE PRECISION
150 /* System generated locals */
152 /* Local variables */
153 int nco, nro, ntot, numb, iiwk, irwk;
154 int i, j, k, kk, ldkey, ldstp, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
155 int i3a, i3b, i3c, i9a, iwkmax, iwkpt;
157 /* Workspace Allocation (freed at end) */
159 iwkmax = 2 * (int) (*workspace / 2);
160 // equiv = (double *) R_alloc(iwkmax / 2, sizeof(double));
161 equiv = (double *) calloc(iwkmax / 2, sizeof(double));
163 /* The check could never happen with Calloc!
164 equiv = Calloc(iwkmax / 2, double);
166 prterr(0, "Can not allocate specified workspace");
170 #define iwrk ((int *)equiv)
171 #define rwrk ((float *)equiv)
173 /* Parameter adjustments */
174 table -= *ldtabl + 1;
180 prterr(1, "NROW must be less than or equal to LDTABL.");
183 for (i = 1; i <= *nrow; ++i) {
184 for (j = 1; j <= *ncol; ++j) {
185 if (table[i + j * *ldtabl] < 0.)
186 prterr(2, "All elements of TABLE must be positive.");
187 ntot = (int) (ntot + table[i + j * *ldtabl]);
191 prterr(3, "All elements of TABLE are zero.\n"
192 "PRT and PRE are set to missing values.");
198 nco = max(*nrow, *ncol);
199 nro = *nrow + *ncol - nco;/* = min(*nrow, *ncol) */
200 k = *nrow + *ncol + 1;
204 i1 = iwork(iwkmax, &iwkpt, ikh, i_real);
205 i2 = iwork(iwkmax, &iwkpt, nco, i_int);
206 i3 = iwork(iwkmax, &iwkpt, nco, i_int);
207 i3a = iwork(iwkmax, &iwkpt, nco, i_int);
208 i3b = iwork(iwkmax, &iwkpt, nro, i_int);
209 i3c = iwork(iwkmax, &iwkpt, nro, i_int);
210 ikh = max(k * 5 + (kk << 1), nco * 7 + 800);
211 iiwk= iwork(iwkmax, &iwkpt, ikh, i_int);
212 ikh = max(nco + 401, k);
213 irwk= iwork(iwkmax, &iwkpt, ikh, i_real);
216 What follows below splits the remaining amount iwkmax - iwkpt of
217 (int) workspace into hash tables as follows.
219 INT 2 * ldkey i4 i5 i11
220 REAL 2 * ldkey i8 i9 i10
223 Hence, we need ldkey times
224 3 * 2 + 3 * 2 * s + 2 * mult * s + 6 * mult
225 chunks of integer memory, where s = sizeof(REAL) / sizeof(INT).
226 If doubles are used and are twice as long as ints, this gives
228 so that the value of ldkey can be obtained by dividing available
229 (int) workspace by this number.
231 In fact, because iwork() can actually s * n + s - 1 int chunks
232 when allocating a REAL, we use ldkey = available / numb - 1.
235 Can we always assume that sizeof(double) / sizeof(int) is 2?
238 if (i_real == 4) { /* Double precision reals */
239 numb = 18 + 10 * mult;
240 } else { /* Single precision reals */
241 numb = (mult << 3) + 12;
243 ldkey = (iwkmax - iwkpt) / numb - 1;
244 ldstp = mult * ldkey;
245 ikh = ldkey << 1; i4 = iwork(iwkmax, &iwkpt, ikh, i_int);
246 ikh = ldkey << 1; i5 = iwork(iwkmax, &iwkpt, ikh, i_int);
247 ikh = ldstp << 1; i6 = iwork(iwkmax, &iwkpt, ikh, i_real);
248 ikh = ldstp * 6; i7 = iwork(iwkmax, &iwkpt, ikh, i_int);
249 ikh = ldkey << 1; i8 = iwork(iwkmax, &iwkpt, ikh, i_real);
250 ikh = ldkey << 1; i9 = iwork(iwkmax, &iwkpt, ikh, i_real);
251 ikh = ldkey << 1; i9a = iwork(iwkmax, &iwkpt, ikh, i_real);
252 ikh = ldkey << 1; i10 = iwork(iwkmax, &iwkpt, ikh, i_int);
254 /* To convert to double precision, change RWRK to DWRK in the next CALL.
296 -----------------------------------------------------------------------
298 Purpose: Computes Fisher's exact test for a contingency table,
299 routine with workspace variables specified.
300 Usage: F2XACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
301 EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF,
302 IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ,
303 DLP, DSP, TM, KEY2, IWK, RWK)
304 -----------------------------------------------------------------------
307 f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
308 double *expect, double *percnt, double *emin, double *prt,
309 double *pre, double *fact, int *ico, int *iro, int *kyy,
310 int *idif, int *irn, int *key, int *ldkey, int *ipoin,
311 double *stp, int *ldstp, int *ifrq, double *dlp, double *dsp,
312 double *tm, int *key2, int *iwk, double *rwk)
314 /* IMAX is the largest representable int on the machine. */
315 const int imax = SINT_MAX;
316 // const int imax = 2147483647; //xx: I DONÂ¥T like this, and
317 // thanks to the hint from Jason Turner I don't do it anymore. (R.D-U).
319 /* AMISS is a missing value indicator which is returned when the
320 probability is not defined. */
321 const double amiss = -12345.;
323 /* TOL is chosen as the square root of the smallest relative spacing. */
325 const static double tol = 3.45254e-7;
327 static double tol = 3.45254e-7;
329 /* EMX is a large positive value used in comparing expected values. */
330 const static double emx = 1e30;
332 /* Local variables {{any really need to be static ???}} */
333 static int kval, kmax, jkey, last, ipsh, itmp, itop, jstp, ntot,
334 jstp2, jstp3, jstp4, i, ii, j, k, n, iflag, ncell, ifreq, chisq,
335 ikkey, ikstp, ikstp2, k1, kb, kd, ks,
336 i31, i32, i33, i34, i35, i36, i37, i38, i39,
337 i41, i42, i43, i44, i45, i46, i47, i48, i310, i311,
338 nco, nrb, ipn, ipo, itp, nro, nro2;
339 static double dspt, dd, df,ddf, drn,dro, emn, obs, obs2, obs3,
348 /* Parameter adjustments */
349 table -= *ldtabl + 1;
367 /* Check table dimensions */
369 prterr(1, "NROW must be less than or equal to LDTABL.");
371 prterr(4, "NCOL must be at least 2");
373 /* Initialize KEY array */
374 for (i = 1; i <= *ldkey << 1; ++i) {
378 /* Initialize parameters */
391 /* nco := max(nrow, ncol) : */
396 /* Initialize pointers for workspace */
410 k = *nrow + *ncol + 1;
420 /* Compute row marginals and total */
422 for (i = 1; i <= *nrow; ++i) {
424 for (j = 1; j <= *ncol; ++j) {
425 if (table[i + j * *ldtabl] < -1e-4)
426 prterr(2, "All elements of TABLE must be positive.");
427 iro[i] += (int) table[i + j * *ldtabl];
433 prterr(3, "All elements of TABLE are zero.\n"
434 "PRT and PRE are set to missing values.");
439 /* Column marginals */
440 for (i = 1; i <= *ncol; ++i) {
442 for (j = 1; j <= *nrow; ++j)
443 ico[i] += (int) table[j + i * *ldtabl];
447 isort(nrow, &iro[1]);
448 isort(ncol, &ico[1]);
450 /* Determine row and column marginals.
451 Define max(nrow,ncol) =: nco >= nro := min(nrow,ncol)
454 Swap marginals if necessary to ico[1:nco] & iro[1:nro]
459 for (i = 1; i <= nco; ++i) {
469 /* Get multiplers for stack */
471 for (i = 2; i <= nro; ++i) {
472 /* Hash table multipliers */
473 if (iro[i - 1] + 1 <= imax / kyy[i - 1]) {
474 kyy[i] = kyy[i - 1] * (iro[i - 1] + 1);
479 /* Maximum product */
480 if (iro[nro - 1] + 1 <= imax / kyy[nro - 1]) {
481 kmax = (iro[nro] + 1) * kyy[nro - 1];
484 prterr(5, "The hash table key cannot be computed because "
486 "is larger than the largest representable int.\n"
487 "The algorithm cannot proceed.\n"
488 "Reduce the workspace size, or use `exact = FALSE'.");
492 /* Compute log factorials */
495 if(ntot >= 2) fact[2] = log(2.);
496 /* MM: old code assuming log() to be SLOW */
497 for (i = 3; i <= ntot; i += 2) {
498 fact[i] = fact[i - 1] + log((double) i);
501 fact[j] = fact[i] + fact[2] + fact[j / 2] - fact[j / 2 - 1];
503 /* Compute obs := observed path length */
506 for (j = 1; j <= nco; ++j) {
508 for (i = 1; i <= nro; ++i) {
510 dd += fact[(int) table[j + i * *ldtabl]];
511 ntot += (int) table[j + i * *ldtabl];
513 dd += fact[(int) table[i + j * *ldtabl]];
514 ntot += (int) table[i + j * *ldtabl];
517 obs += fact[ico[j]] - dd;
519 /* Denominator of observed table: DRO */
520 dro = f9xact(&nro, &ntot, &iro[1], fact);
521 *prt = exp(obs - dro);
522 /* Initialize pointers */
527 jstp2 = *ldstp * 3 + 1;
528 jstp3 = (*ldstp << 2) + 1;
529 jstp4 = *ldstp * 5 + 1;
532 ikstp2 = *ldstp << 1;
537 ifrq[ikstp2 + 1] = -1;
545 /* IDIF is the difference in going to the daughter */
546 for (i = 1; i <= nro; ++i)
549 /* Generate the first daughter */
552 ntot = min(n, iro[kd]);
558 while (n > 0 && kd != 1);
567 for (i = kb + 1; i <= nco; ++i)
572 /* Arc to daughter length=ICO(KB) */
573 for (i = 1; i <= nro; ++i)
574 irn[i] = iro[i] - idif[i];
579 if (irn[1] > irn[2]) {
584 } else if (nro == 3) {
588 if (irn[2] > irn[3]) {
601 } else if (ii > irn[2]) {
604 } else if (irn[2] > irn[3]) {
610 for (j = 2; j <= nro; ++j) {
614 while (ii < irn[i]) {
623 /* Adjust start for zero */
624 for (i = 1; i <= nro; ++i) {
635 /* Some table values */
636 ddf = f9xact(&nro, &n, &idif[1], fact);
637 drn = f9xact(&nro2, &ntot, &irn[nrb], fact) - dro + ddf;
640 kval = irn[1] + irn[2] * kyy[2];
641 for (i = 3; i <= nro; ++i) {
642 kval += irn[i] * kyy[i];
644 /* Get hash table entry */
645 i = kval % (*ldkey << 1) + 1;
646 /* Search for unused location */
647 for (itp = i; itp <= *ldkey << 1; ++itp) {
659 for (itp = 1; itp <= i - 1; ++itp) {
671 prterr(6, "LDKEY is too small.\n"
672 "It is not possible to give the value of LDKEY required,\n"
673 "but you could try doubling LDKEY (and possibly LDSTP).");
675 prterr(6, "LDKEY is too small for this problem.\n"
676 "Try increasing the size of the workspace.");
682 ipn = ipoin[ipo + ikkey];
683 pastp = stp[ipn + ikstp];
684 ifreq = ifrq[ipn + ikstp];
685 /* Compute shortest and longest path */
687 obs2 = obs - fact[ico[kb + 1]] - fact[ico[kb + 2]] - ddf;
688 for (i = 3; i <= k1; ++i) {
689 obs2 -= fact[ico[kb + i]];
692 dspt = obs - obs2 - ddf;
693 /* Compute longest path */
695 f3xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dlp[itp],
696 &ntot, fact, &iwk[i31], &iwk[i32], &iwk[i33],
697 &iwk[i34], &iwk[i35], &iwk[i36], &iwk[i37],
698 &iwk[i38], &iwk[i39], &rwk[i310], &rwk[i311], &tol);
699 dlp[itp] = min(0., dlp[itp]);
700 /* Compute shortest path */
702 f4xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dsp[itp], fact,
703 &iwk[i47], &iwk[i41], &iwk[i42], &iwk[i43],
704 &iwk[i44], &iwk[i45], &iwk[i46], &rwk[i48], &tol);
705 dsp[itp] = min(0., dsp[itp] - dspt);
706 /* Use chi-squared approximation? */
707 if ((irn[nrb] * ico[kb + 1]) > ntot * emn) {
709 for (i = 0; i < nro2; ++i)
710 for (j = 1; j <= k1; ++j)
711 if (irn[nrb + i] * ico[kb + j] >= ntot * *expect)
714 if (ncell * 100 >= k1 * nro2 * *percnt) {
716 for (i = 0; i < nro2; ++i)
717 tmp += (fact[irn[nrb + i]] -
718 fact[irn[nrb + i] - 1]);
720 for (j = 1; j <= k1; ++j)
721 tmp += (nro2 - 1) * (fact[ico[kb + j]] -
722 fact[ico[kb + j] - 1]);
723 df = (double) ((nro2 - 1) * (k1 - 1));
724 tmp += df * 1.83787706640934548356065947281;
725 tmp -= (nro2 * k1 - 1) * (fact[ntot] - fact[ntot - 1]);
726 tm[itp] = (obs - dro) * -2. - tmp;
728 /* tm(itp) set to a flag value */
735 obs3 = obs2 - dlp[itp];
737 if (tm[itp] == -9876.) {
744 obs2 = obs - drn - dro;
749 /* Process node with new PASTP */
752 *pre += (double) ifreq * exp(pastp + drn);
753 } else if (pastp < obs2) {
755 df = (double) ((nro2 - 1) * (k1 - 1));
757 pv = pgamma(fmax2(0., tmp + (pastp + drn) * 2.) / 2.,
758 df / 2., /*scale = */ 1.,
759 /*lower_tail = */FALSE, /*log_p = */ FALSE);
761 d1 = max(0., tmp + (pastp + drn) * 2.) / 2.;
763 pv = 1. - gammds(&d1, &d2, &ifault);
765 *pre += (double) ifreq * exp(pastp + drn) * pv;
767 /* Put daughter on queue */
769 f5xact(&d1, &tol, &kval, &key[jkey], ldkey, &ipoin[jkey],
770 &stp[jstp], ldstp, &ifrq[jstp], &ifrq[jstp2],
771 &ifrq[jstp3], &ifrq[jstp4], &ifreq, &itop, &ipsh);
775 /* Get next PASTP on chain */
776 ipn = ifrq[ipn + ikstp2];
778 pastp = stp[ipn + ikstp];
779 ifreq = ifrq[ipn + ikstp];
782 /* Generate a new daughter node */
783 f7xact(&kmax, &iro[1], &idif[1], &kd, &ks, &iflag);
789 /* Go get a new mother from stage K */
792 f6xact(&nro, &iro[1], &iflag, &kyy[1], &key[ikkey + 1], ldkey,
794 /* Update pointers */
797 /* else iflag == 3 : no additional nodes to process */
803 jkey = *ldkey - jkey + 2;
804 jstp = *ldstp - jstp + 2;
805 jstp2 = (*ldstp << 1) + jstp;
806 for (i = 1; i <= *ldkey << 1; ++i)
813 -----------------------------------------------------------------------
815 Purpose: Computes the shortest path length for a given table.
816 Usage: F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, IRO,
817 IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, TOL)
819 NROW - The number of rows in the table. (Input)
820 IROW - Vector of length NROW containing the row sums
821 for the table. (Input)
822 NCOL - The number of columns in the table. (Input)
823 ICOL - Vector of length K containing the column sums
824 for the table. (Input)
825 DLP - The longest path for the table. (Output)
826 MM - The total count in the table. (Output)
827 FACT - Vector containing the logarithms of factorials. (Input)
828 ICO - Work vector of length MAX(NROW,NCOL).
829 IRO - Work vector of length MAX(NROW,NCOL).
830 IT - Work vector of length MAX(NROW,NCOL).
831 LB - Work vector of length MAX(NROW,NCOL).
832 NR - Work vector of length MAX(NROW,NCOL).
833 NT - Work vector of length MAX(NROW,NCOL).
834 NU - Work vector of length MAX(NROW,NCOL).
835 ITC - Work vector of length 400.
836 IST - Work vector of length 400.
837 STV - Work vector of length 400.
838 ALEN - Work vector of length MAX(NROW,NCOL).
839 TOL - Tolerance. (Input)
840 -----------------------------------------------------------------------
844 f3xact(int *nrow, int *irow, int *ncol, int *icol, double *dlp,
845 int *mm, double *fact, int *ico, int *iro, int *it,
846 int *lb, int *nr, int *nt, int *nu, int *itc, int *ist,
847 double *stv, double *alen, const double *tol)
849 /* Initialized data */
850 static int ldst = 200;
854 /* Local variables */
858 static int n11, n12, ii, nn, ks, ic1, ic2, nc1, nn1;
861 static int nct, ipn, irl, key, lev, itp, nro;
863 static int nrt, kyy, nc1s;
865 /* Parameter adjustments */
880 for (i = 0; i <= *ncol; ++i) {
883 for (i = 1; i <= 400; ++i) {
889 *dlp -= fact[icol[1]];
890 for (i = 2; i <= *ncol; ++i) {
891 *dlp -= fact[icol[i]];
899 *dlp = *dlp - fact[irow[1]] - fact[irow[2]];
900 for (i = 3; i <= *nrow; ++i) {
901 *dlp -= fact[irow[i]];
907 if (*nrow * *ncol == 4) {
908 n11 = (irow[1] + 1) * (icol[1] + 1) / (*mm + 2);
910 *dlp = *dlp - fact[n11] - fact[n12] - fact[icol[1] - n11]
911 - fact[icol[2] - n12];
914 /* Test for optimal table */
917 if (irow[*nrow] <= irow[1] + *ncol) {
918 f10act(nrow, &irow[1], ncol, &icol[1], &val, &xmin, fact,
919 &lb[1], &nu[1], &nr[1]);
922 if (icol[*ncol] <= icol[1] + *nrow) {
923 f10act(ncol, &icol[1], nrow, &irow[1], &val, &xmin, fact,
924 &lb[1], &nu[1], &nr[1]);
932 /* Setup for dynamic programming */
935 if (*nrow >= *ncol) {
938 for (i = 1; i <= *nrow; ++i) {
943 for (i = 2; i <= *ncol; ++i) {
945 nt[i] = nt[i - 1] - ico[i];
952 for (i = 2; i <= *nrow; ++i) {
954 nt[i] = nt[i - 1] - ico[i];
956 for (i = 1; i <= *ncol; ++i)
959 /* Initialize pointers */
967 LnewNode: /* Setup to generate new node */
973 lb[1] = (int) ((double) ((nrt + 1) * (nct + 1)) /
974 (double) (nn + nr1 * nc1s + 1) - *tol) - 1;
975 nu[1] = (int) ((double) ((nrt + nc1s) * (nct + nr1)) /
976 (double) (nn + nr1 + nc1s)) - lb[1] + 1;
979 LoopNode: /* Generate a node */
991 alen[lev] = alen[lev - 1] + fact[lb[lev]];
998 lb[lev] = (int) ((double) ((nrt + 1) * (nct + 1)) /
999 (double) (nn1 + nr1 * nc1 + 1) - *tol);
1000 nu[lev] = (int) ((double) ((nrt + nc1) * (nct + nr1)) /
1001 (double) (nn1 + nr1 + nc1) - lb[lev] + 1);
1002 nr[lev] = nrt - lb[lev];
1005 alen[nco] = alen[lev] + fact[nr[lev]];
1008 v = val + alen[nco];
1010 /* Only 1 row left */
1011 v = v + fact[ico[1] - lb[1]] + fact[ico[2] - lb[2]];
1012 for (i = 3; i <= nco; ++i) {
1013 v += fact[ico[i] - lb[i]];
1018 } else if (nro == 3 && nco == 2) {
1019 /* 3 rows and 2 columns */
1020 nn1 = nn - iro[irl] + 2;
1021 ic1 = ico[1] - lb[1];
1022 ic2 = ico[2] - lb[2];
1023 n11 = (iro[irl + 1] + 1) * (ic1 + 1) / nn1;
1024 n12 = iro[irl + 1] - n11;
1025 v = v + fact[n11] + fact[n12] + fact[ic1 - n11]
1031 /* Column marginals are new node */
1032 for (i = 1; i <= nco; ++i) {
1033 it[i] = ico[i] - lb[i];
1035 /* Sort column marginals */
1037 if (it[1] > it[2]) {
1042 } else if (nco == 3) {
1046 if (it[2] > it[3]) {
1059 } else if (ii > it[2]) {
1062 } else if (it[2] > it[3]) {
1068 isort(&nco, &it[1]);
1070 /* Compute hash value */
1071 key = it[1] * kyy + it[2];
1072 for (i = 3; i <= nco; ++i) {
1073 key = it[i] + key * kyy;
1076 //PROBLEM "Bug in FEXACT: gave negative key" RECOVER(NULL_ENTRY);
1077 printf("Bug in FEXACT: gave negative key \n"); //xx:another one of my ugly kluges (R.D-U)
1080 ipn = key % ldst + 1;
1082 /* Find empty position */
1083 for (itp = ipn, ii = ks + ipn; itp <= ldst; ++itp, ++ii) {
1086 } else if (ist[ii] == key) {
1091 for (itp = 1, ii = ks + 1; itp <= ipn - 1; ++itp, ++ii) {
1094 } else if (ist[ii] == key) {
1099 prterr(30, "Stack length exceeded in f3xact.\n"
1100 "This problem should not occur.");
1102 L180: /* Push onto stack */
1110 L190: /* Marginals already on stack */
1111 stv[ii] = min(v, stv[ii]);
1116 L200: /* Pop item from stack */
1119 itp = itc[nitc + k] + k;
1124 /* Compute marginals */
1125 for (i = nco; i >= 2; --i) {
1130 /* Set up nt array */
1131 nt[1] = nn - ico[1];
1132 for (i = 2; i <= nco; ++i)
1133 nt[i] = nt[i - 1] - ico[i];
1135 /* Test for optimality (L90) */
1137 if (iro[nro] <= iro[irl] + nco) {
1138 f10act(&nro, &iro[irl], &nco, &ico[1], &val, &xmin, fact,
1139 &lb[1], &nu[1], &nr[1]);
1141 if (!xmin && ico[nco] <= ico[1] + nro)
1142 f10act(&nco, &ico[1], &nro, &iro[irl], &val, &xmin, fact,
1143 &lb[1], &nu[1], &nr[1]);
1151 } else if (nro > 2 && nst > 0) {
1152 /* Go to next level */
1167 -----------------------------------------------------------------------
1169 Purpose: Computes the longest path length for a given table.
1170 Usage: CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK,
1171 NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK,
1174 NROW - The number of rows in the table. (Input)
1175 IROW - Vector of length NROW containing the row sums for the
1177 NCOL - The number of columns in the table. (Input)
1178 ICOL - Vector of length K containing the column sums for the
1180 DSP - The shortest path for the table. (Output)
1181 FACT - Vector containing the logarithms of factorials. (Input)
1182 ICSTK - NCOL by NROW+NCOL+1 work array.
1183 NCSTK - Work vector of length NROW+NCOL+1.
1184 LSTK - Work vector of length NROW+NCOL+1.
1185 MSTK - Work vector of length NROW+NCOL+1.
1186 NSTK - Work vector of length NROW+NCOL+1.
1187 NRSTK - Work vector of length NROW+NCOL+1.
1188 IRSTK - NROW by MAX(NROW,NCOL) work array.
1189 YSTK - Work vector of length NROW+NCOL+1.
1190 TOL - Tolerance. (Input)
1191 -----------------------------------------------------------------------
1195 f4xact(int *nrow, int *irow, int *ncol, int *icol, double *dsp,
1196 double *fact, int *icstk, int *ncstk, int *lstk, int *mstk,
1197 int *nstk, int *nrstk, int *irstk, double *ystk, const double *tol)
1199 /* System generated locals */
1202 /* Local variables */
1203 int i, j, k, l, m, n, mn, ic1, ir1, ict, irt, istk, nco, nro;
1206 /* Parameter adjustments */
1219 /* Take care of the easy cases first */
1221 for (i = 1; i <= *ncol; ++i) {
1222 *dsp -= fact[icol[i]];
1227 for (i = 1; i <= *nrow; ++i) {
1228 *dsp -= fact[irow[i]];
1232 if (*nrow * *ncol == 4) {
1233 if (irow[2] <= icol[2]) {
1234 *dsp = *dsp - fact[irow[2]] - fact[icol[1]]
1235 - fact[icol[2] - irow[2]];
1237 *dsp = *dsp - fact[icol[2]] - fact[irow[1]]
1238 - fact[irow[2] - icol[2]];
1242 /* initialization before loop */
1243 for (i = 1; i <= *nrow; ++i) {
1244 irstk[i + *nrow] = irow[*nrow - i + 1];
1246 for (j = 1; j <= *ncol; ++j) {
1247 icstk[j + *ncol] = icol[*ncol - j + 1];
1262 ir1 = irstk[istk * *nrow + 1];
1263 ic1 = icstk[istk * *ncol + 1];
1270 } else if (ir1 < ic1) {
1291 irt = irstk[i + istk * *nrow];
1292 ict = icstk[j + istk * *ncol];
1301 f11act(&irstk[istk * *nrow + 1], &i, &nro,
1302 &irstk[(istk + 1) * *nrow + 1]);
1303 f11act(&icstk[istk * *ncol + 1], &j, &nco,
1304 &icstk[(istk + 1) * *ncol + 1]);
1305 } else if (irt > ict) {
1307 f11act(&icstk[istk * *ncol + 1], &j, &nco,
1308 &icstk[(istk + 1) * *ncol + 1]);
1310 f8xact(&irstk[istk * *nrow + 1], &ikh, &i,
1311 &nro, &irstk[(istk + 1) * *nrow + 1]);
1314 f11act(&irstk[istk * *nrow + 1], &i, &nro,
1315 &irstk[(istk + 1) * *nrow + 1]);
1317 f8xact(&icstk[istk * *ncol + 1], &ikh, &j,
1318 &nco, &icstk[(istk + 1) * *ncol + 1]);
1322 for (k = 1; k <= nco; ++k) {
1323 y += fact[icstk[k + (istk + 1) * *ncol]];
1328 for (k = 1; k <= nro; ++k) {
1329 y += fact[irstk[k + (istk + 1) * *nrow]];
1342 } while(1);/* end do */
1347 if (*dsp - amx <= *tol) {
1357 if (*dsp - amx <= *tol) {
1366 if (l > mstk[istk]) goto L100;
1373 if (irstk[l + istk * *nrow] <
1374 irstk[l - 1 + istk * *nrow]) goto L60;
1377 if (icstk[l + istk * *ncol] <
1378 icstk[l - 1 + istk * *ncol]) goto L60;
1384 -----------------------------------------------------------------------
1386 Purpose: Put node on stack in network algorithm.
1387 Usage: CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP,
1388 LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP,
1391 PASTP - The past path length. (Input)
1392 TOL - Tolerance for equivalence of past path lengths. (Input)
1393 KVAL - Key value. (Input)
1394 KEY - Vector of length LDKEY containing the key values. (in/out)
1395 LDKEY - Length of vector KEY. (Input)
1396 IPOIN - Vector of length LDKEY pointing to the
1397 linked list of past path lengths. (in/out)
1398 STP - Vector of length LSDTP containing the
1399 linked lists of past path lengths. (in/out)
1400 LDSTP - Length of vector STP. (Input)
1401 IFRQ - Vector of length LDSTP containing the past path
1402 frequencies. (in/out)
1403 NPOIN - Vector of length LDSTP containing the pointers to
1404 the next past path length. (in/out)
1405 NR - Vector of length LDSTP containing the right object
1406 pointers in the tree of past path lengths. (in/out)
1407 NL - Vector of length LDSTP containing the left object
1408 pointers in the tree of past path lengths. (in/out)
1409 IFREQ - Frequency of the current path length. (Input)
1410 ITOP - Pointer to the top of STP. (Input)
1411 IPSH - Option parameter. (Input)
1412 If IPSH is true, the past path length is found in the
1413 table KEY. Otherwise the location of the past path
1414 length is assumed known and to have been found in
1415 a previous call. ==>>>>> USING "static" variables
1416 -----------------------------------------------------------------------
1420 f5xact(double *pastp, const double *tol, int *kval, int *key, int *ldkey,
1421 int *ipoin, double *stp, int *ldstp, int *ifrq, int *npoin,
1422 int *nr, int *nl, int *ifreq, int *itop, int *ipsh)
1424 /* Local variables */
1425 static int itmp, ird, ipn, itp;
1426 double test1, test2;
1428 /* Parameter adjustments */
1439 /* Convert KVAL to int in range 1, ..., LDKEY. */
1440 ird = *kval % *ldkey + 1;
1441 /* Search for an unused location */
1442 for (itp = ird; itp <= *ldkey; ++itp) {
1443 if (key[itp] == *kval) {
1450 for (itp = 1; itp <= ird - 1; ++itp) {
1451 if (key[itp] == *kval) {
1458 /* Return if KEY array is full */
1460 prterr(6, "LDKEY is too small for this problem.\n"
1461 "It is not possible to estimate the value of LDKEY "
1463 "but twice the current value may be sufficient.");
1465 prterr(6, "LDKEY is too small for this problem.\n"
1466 "Try increasing the size of the workspace.");
1473 /* Return if STP array full */
1474 if (*itop > *ldstp) {
1476 prterr(7, "LDSTP is too small for this problem.\n"
1477 "It is not possible to estimate the value of LDSTP "
1479 "but twice the current value may be sufficient.");
1481 prterr(7, "LDSTP is too small for this problem.\n"
1482 "Try increasing the size of the workspace.");
1484 /* Update STP, etc. */
1488 stp[*itop] = *pastp;
1489 ifrq[*itop] = *ifreq;
1493 /* Find location, if any, of pastp */
1496 test1 = *pastp - *tol;
1497 test2 = *pastp + *tol;
1500 if (stp[ipn] < test1) {
1505 } else if (stp[ipn] > test2) {
1511 ifrq[ipn] += *ifreq;
1514 /* Return if STP array full */
1516 if (*itop > *ldstp) {
1518 prterr(7, "LDSTP is too small for this problem.\n"
1519 "It is not possible to estimate the value of LDSTP "
1521 "but twice the current value may be sufficient.");
1523 prterr(7, "LDSTP is too small for this problem.\n"
1524 "Try increasing the size of the workspace.");
1527 /* Find location to add value */
1532 if (stp[ipn] < test1) {
1540 } else if (stp[ipn] > test2) {
1549 /* Update STP, etc. */
1550 npoin[*itop] = npoin[itmp];
1551 npoin[itmp] = *itop;
1552 stp[*itop] = *pastp;
1553 ifrq[*itop] = *ifreq;
1559 -----------------------------------------------------------------------
1561 Purpose: Pop a node off the stack.
1562 Usage: CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, IPN)
1564 NROW - The number of rows in the table. (Input)
1565 IROW - Vector of length nrow containing the row sums on
1567 IFLAG - Set to 3 if there are no additional nodes to process.
1569 KYY - Constant mutlipliers used in forming the hash
1571 KEY - Vector of length LDKEY containing the hash table
1573 LDKEY - Length of vector KEY. (Input)
1574 LAST - Index of the last key popped off the stack. (In/out)
1575 IPN - Pointer to the linked list of past path lengths. (Output)
1576 -----------------------------------------------------------------------
1579 f6xact(int *nrow, int *irow, int *iflag, int *kyy, int *key, int
1580 *ldkey, int *last, int *ipn)
1584 /* Parameter adjustments */
1592 if (*last <= *ldkey) {
1593 if (key[*last] < 0) {
1596 /* Get KVAL from the stack */
1599 for (j = *nrow; j >= 2; --j) {
1600 irow[j] = kval / kyy[j];
1601 kval -= irow[j] * kyy[j];
1613 -----------------------------------------------------------------------
1615 Purpose: Generate the new nodes for given marinal totals.
1616 Usage: CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG)
1618 NROW - The number of rows in the table. (Input)
1619 IMAX - The row marginal totals. (Input)
1620 IDIF - The column counts for the new column. (in/out)
1621 K - Indicator for the row to decrement. (in/out)
1622 KS - Indicator for the row to increment. (in/out)
1623 IFLAG - Status indicator. (Output)
1624 If IFLAG is zero, a new table was generated. For
1625 IFLAG = 1, no additional tables could be generated.
1626 -----------------------------------------------------------------------
1630 f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
1636 /* Parameter adjustments */
1642 /* Find node which can be incremented, ks */
1646 } while (idif[*ks] == imax[*ks]);
1648 /* Find node to decrement (>ks) */
1649 if (idif[*k] > 0 && *k > *ks) {
1653 } while(imax[*k] == 0);
1657 /* Find node to increment (>=ks) */
1658 while (idif[m] >= imax[m]) {
1664 if (idif[m] == imax[m]) {
1671 /* Check for finish */
1672 for (k1 = *k + 1; k1 <= *nrow; ++k1) {
1681 /* Reallocate counts */
1683 for (i = 1; i <= *k; ++i) {
1691 m = min(mm, imax[*k]);
1694 } while (mm > 0 && *k != 1);
1696 /* Check that all counts reallocated */
1713 } while (idif[*ks] >= imax[*ks]);
1718 -----------------------------------------------------------------------
1720 Purpose: Routine for reducing a vector when there is a zero
1722 Usage: CALL F8XACT (IROW, IS, I1, IZERO, NEW)
1724 IROW - Vector containing the row counts. (Input)
1725 IS - Indicator. (Input)
1726 I1 - Indicator. (Input)
1727 IZERO - Position of the zero. (Input)
1728 NEW - Vector of new row counts. (Output)
1729 -----------------------------------------------------------------------
1733 f8xact(int *irow, int *is, int *i1, int *izero, int *myNew)
1737 /* Parameter adjustments */
1742 for (i = 1; i < *i1; ++i)
1745 for (i = *i1; i <= *izero - 1; ++i) {
1746 if (*is >= irow[i + 1])
1748 myNew[i] = irow[i + 1];
1762 -----------------------------------------------------------------------
1764 Purpose: Computes the log of a multinomial coefficient.
1765 Usage: F9XACT(N, MM, IR, FACT)
1767 N - Length of IR. (Input)
1768 MM - Number for factorial in numerator. (Input)
1769 IR - Vector of length N containing the numbers for
1770 the denominator of the factorial. (Input)
1771 FACT - Table of log factorials. (Input)
1772 F9XACT - The log of the multinomal coefficient. (Output)
1773 -----------------------------------------------------------------------
1777 f9xact(int *n, int *mm, int *ir, double *fact)
1783 for (k = 0; k < *n; k++)
1789 -----------------------------------------------------------------------
1791 Purpose: Computes the shortest path length for special tables.
1792 Usage: F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, NE, M)
1794 NROW - The number of rows in the table. (Input)
1795 IROW - Vector of length NROW containing the row totals. (Input)
1796 NCOL - The number of columns in the table. (Input)
1797 ICO - Vector of length NCOL containing the column totals.(Input)
1798 VAL - The shortest path. (Output)
1799 XMIN - Set to true if shortest path obtained. (Output)
1800 FACT - Vector containing the logarithms of factorials. (Input)
1801 ND - Workspace vector of length NROW. (Input)
1802 NE - Workspace vector of length NCOL. (Input)
1803 M - Workspace vector of length NCOL. (Input)
1805 Chapter: STAT/LIBRARY Categorical and Discrete Data Analysis
1806 -----------------------------------------------------------------------
1810 f10act(int *nrow, int *irow, int *ncol, int *icol, double *val,
1811 int *xmin, double *fact, int *nd, int *ne, int *m)
1813 /* Local variables */
1814 int i, is, ix, nrw1;
1816 /* Parameter adjustments */
1824 for (i = 1; i <= *nrow - 1; ++i)
1827 is = icol[1] / *nrow;
1828 ix = icol[1] - *nrow * is;
1834 for (i = 2; i <= *ncol; ++i) {
1835 ix = icol[i] / *nrow;
1838 ix = icol[i] - *nrow * ix;
1844 for (i = *nrow - 2; i >= 1; --i)
1849 for (i = *nrow; i >= 2; --i) {
1850 ix = ix + is + nd[nrw1 - i] - irow[i];
1855 for (i = 1; i <= *ncol; ++i) {
1858 *val = *val + is * fact[ix + 1] + (*nrow - is) * fact[ix];
1866 -----------------------------------------------------------------------
1868 Purpose: Routine for revising row totals.
1869 Usage: CALL F11ACT (IROW, I1, I2, NEW)
1871 IROW - Vector containing the row totals. (Input)
1872 I1 - Indicator. (Input)
1873 I2 - Indicator. (Input)
1874 NEW - Vector containing the row totals. (Output)
1875 -----------------------------------------------------------------------
1878 f11act(int *irow, int *i1, int *i2, int *myNew)
1882 /* Parameter adjustments */
1886 for (i = 1; i <= (*i1 - 1); ++i) myNew[i] = irow[i];
1887 for (i = *i1; i <= *i2; ++i) myNew[i] = irow[i + 1];
1893 -----------------------------------------------------------------------
1895 Purpose: Print an error message and stop.
1896 Usage: prterr(icode, mes)
1898 icode - Integer code for the error message. (Input)
1899 mes - Character string containing the error message. (Input)
1900 -----------------------------------------------------------------------
1903 prterr(int icode, char *mes)
1905 // PROBLEM "FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY);
1906 // printf("FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY));
1907 printf("FEXACT error %d.\n", icode); //xx:another one of my ugly kluges
1912 -----------------------------------------------------------------------
1914 Purpose: Routine for allocating workspace.
1915 Usage: iwork (iwkmax, iwkpt, number, itype)
1917 iwkmax - Maximum (int) amount of workspace. (Input)
1918 iwkpt - Amount of (int) workspace currently allocated. (in/out)
1919 number - Number of elements of workspace desired. (Input)
1920 itype - Workspace type. (Input)
1925 iwork(): Index in rwrk, dwrk, or iwrk of the beginning of
1926 the first free element in the workspace array. (Output)
1927 -----------------------------------------------------------------------
1930 iwork(int iwkmax, int *iwkpt, int number, int itype)
1935 if (itype == 2 || itype == 3)
1940 *iwkpt += (number << 1);
1943 if (*iwkpt > iwkmax)
1944 prterr(40, "Out of workspace.");
1951 void isort(int *n, int *ix)
1954 -----------------------------------------------------------------------
1956 Purpose: Shell sort for an int vector.
1957 Usage: CALL ISORT (N, IX)
1959 N - Lenth of vector IX. (Input)
1960 IX - Vector to be sorted. (in/out)
1961 -----------------------------------------------------------------------
1963 static int ikey, i, j, m, il[10], kl, it, iu[10], ku;
1965 /* Parameter adjustments */
1981 /* Find element in first half */
1985 if (ix[ikey] > ix[i]) {
1989 /* Find element in second half */
1992 if (ix[j] > ix[ikey]) {
2005 /* Save upper and lower subscripts of the array yet to be sorted */
2007 if (j - kl < ku - j) {
2021 prterr(20, "This should never occur.");
2023 /* Use another segment */
2034 double gammds(double *y, double *p, int *ifault)
2037 -----------------------------------------------------------------------
2039 Purpose: Cumulative distribution for the gamma distribution.
2040 Usage: PGAMMA (Q, ALPHA,IFAULT)
2042 Q - Value at which the distribution is desired. (Input)
2043 ALPHA - Parameter in the gamma distribution. (Input)
2044 IFAULT - Error indicator. (Output)
2047 1 An argument is misspecified.
2048 2 A numerical error has occurred.
2049 PGAMMA - The cdf for the gamma distribution with parameter alpha
2050 evaluated at Q. (Output)
2051 -----------------------------------------------------------------------
2053 Algorithm AS 147 APPL. Statist. (1980) VOL. 29, P. 113
2055 Computes the incomplete gamma integral for positive parameters Y, P
2056 using and infinite series.
2059 static double a, c, f, g;
2062 /* Checks for the admissibility of arguments and value of F */
2065 if (*y <= 0. || *p <= 0.) {
2071 ALOGAM is natural log of gamma function no need to test ifail as
2072 an error is impossible
2076 f = exp(*p * log(*y) - alogam(&a, &ifail) - *y);
2098 -----------------------------------------------------------------------
2100 Purpose: Value of the log-gamma function.
2101 Usage: ALOGAM (X, IFAULT)
2103 X - Value at which the log-gamma function is to be evaluated.
2105 IFAULT - Error indicator. (Output)
2109 ALGAMA - The value of the log-gamma function at XX. (Output)
2110 -----------------------------------------------------------------------
2112 Algorithm ACM 291, Comm. ACM. (1966) Vol. 9, P. 684
2114 Evaluates natural logarithm of gamma(x) for X greater than zero.
2117 double alogam(double *x, int *ifault)
2119 /* Initialized data */
2121 static double a1 = .918938533204673;
2122 static double a2 = 5.95238095238e-4;
2123 static double a3 = 7.93650793651e-4;
2124 static double a4 = .002777777777778;
2125 static double a5 = .083333333333333;
2127 /* Local variables */
2128 static double f, y, z;
2152 return(f + (y - .5) * log(y) - y + a1 +
2153 (((-a2 * z + a3) * z - a4) * z + a5) / y);
2157 #endif /* not USING_R */