8 #define SINT_MAX INT_MAX
10 #define max(a, b) ((a) < (b) ? (b) : (a))
11 #define min(a, b) ((a) > (b) ? (b) : (a))
13 static void f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
14 double *expect, double *percnt, double *emin, double
15 *prt, double *pre, double *fact, int *ico, int
16 *iro, int *kyy, int *idif, int *irn, int *key,
17 int *ldkey, int *ipoin, double *stp, int *ldstp,
18 int *ifrq, double *dlp, double *dsp, double *tm,
19 int *key2, int *iwk, double *rwk);
20 static void f3xact(int *nrow, int *irow, int *ncol, int *icol,
21 double *dlp, int *mm, double *fact, int *ico, int
22 *iro, int *it, int *lb, int *nr, int *nt, int
23 *nu, int *itc, int *ist, double *stv, double *alen,
25 static void f4xact(int *nrow, int *irow, int *ncol, int *icol,
26 double *dsp, double *fact, int *icstk, int *ncstk,
27 int *lstk, int *mstk, int *nstk, int *nrstk, int
28 *irstk, double *ystk, const double *tol);
29 static void f5xact(double *pastp, const double *tol, int *kval, int *key,
30 int *ldkey, int *ipoin, double *stp, int *ldstp,
31 int *ifrq, int *npoin, int *nr, int *nl, int
32 *ifreq, int *itop, int *ipsh);
33 static void f6xact(int *nrow, int *irow, int *iflag, int *kyy,
34 int *key, int *ldkey, int *last, int *ipn);
35 static void f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
37 static void f8xact(int *irow, int *is, int *i1, int *izero, int *new);
38 static double f9xact(int *n, int *mm, int *ir, double *fact);
39 static void f10act(int *nrow, int *irow, int *ncol, int *icol,
40 double *val, int *xmin, double *fact, int *nd,
42 static void f11act(int *irow, int *i1, int *i2, int *new);
43 static void prterr(int icode, char *mes);
44 static int iwork(int iwkmax, int *iwkpt, int number, int itype);
45 // void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
46 // double *expect, double *percnt, double *emin, double *prt,
47 // double *pre, /* new in C : */ int *workspace);
48 static void isort(int *n, int *ix);
49 static double gammds(double *y, double *p, int *ifault);
50 static double alogam(double *x, int *ifault);
55 /* The only public function : */
57 fexact(int *nrow, int *ncol, double *table, int *ldtabl,
58 double *expect, double *percnt, double *emin, double *prt,
59 double *pre, /* new in C : */ int *workspace)
63 ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
64 THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
65 VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
66 -----------------------------------------------------------------------
68 Purpose: Computes Fisher's exact test probabilities and a hybrid
69 approximation to Fisher exact test probabilities for a
70 contingency table using the network algorithm.
71 Usage: CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
74 NROW - The number of rows in the table. (Input)
75 NCOL - The number of columns in the table. (Input)
76 TABLE - NROW by NCOL matrix containing the contingency
78 LDTABL - Leading dimension of TABLE exactly as specified
79 in the dimension statement in the calling
81 EXPECT - Expected value used in the hybrid algorithm for
82 deciding when to use asymptotic theory
83 probabilities. (Input)
84 If EXPECT <= 0.0 then asymptotic theory probabilities
85 are not used and Fisher exact test probabilities are
86 computed. Otherwise, if PERCNT or more of the cells in
87 the remaining table have estimated expected values of
88 EXPECT or more, with no remaining cell having expected
89 value less than EMIN, then asymptotic chi-squared
90 probabilities are used. See the algorithm section of the
91 manual document for details.
92 Use EXPECT = 5.0 to obtain the 'Cochran' condition.
93 PERCNT - Percentage of remaining cells that must have
94 estimated expected values greater than EXPECT
95 before asymptotic probabilities can be used. (Input)
96 See argument EXPECT for details.
97 Use PERCNT = 80.0 to obtain the 'Cochran' condition.
98 EMIN - Minimum cell estimated expected value allowed for
99 asymptotic chi-squared probabilities to be used. (Input)
100 See argument EXPECT for details.
101 Use EMIN = 1.0 to obtain the 'Cochran' condition.
102 PRT - Probability of the observed table for fixed
103 marginal totals. (Output)
104 PRE - Table p-value. (Output)
105 PRE is the probability of a more extreme table,
106 where `extreme' is in a probabilistic sense.
107 If EXPECT < 0 then the Fisher exact probability
108 is returned. Otherwise, an approximation to the
109 Fisher exact probability is computed based upon
110 asymptotic chi-squared probabilities for ``large''
111 table expected values. The user defines ``large''
112 through the arguments EXPECT, PERCNT, and EMIN.
115 1. For many problems one megabyte or more of workspace can be
116 required. If the environment supports it, the user should begin
117 by increasing the workspace used to 200,000 units.
118 2. In FEXACT, LDSTP = 30*LDKEY. The proportion of table space used
119 by STP may be changed by changing the line MULT = 30 below to
121 3. FEXACT may be converted to single precision by setting IREAL = 3,
122 and converting all DOUBLE PRECISION specifications (except the
123 specifications for RWRK, IWRK, and DWRK) to REAL. This will
124 require changing the names and specifications of the intrinsic
125 functions ALOG, AMAX1, AMIN1, EXP, and REAL. In addition, the
126 machine specific constants will need to be changed, and the name
127 DWRK will need to be changed to RWRK in the call to F2XACT.
128 4. Machine specific constants are specified and documented in F2XACT.
129 A missing value code is specified in both FEXACT and F2XACT.
130 5. Although not a restriction, is is not generally practical to call
131 this routine with large tables which are not sparse and in
132 which the 'hybrid' algorithm has little effect. For example,
133 although it is feasible to compute exact probabilities for the
138 computing exact probabilities for a similar table which has been
139 enlarged by the addition of an extra row (or column) may not be
141 -----------------------------------------------------------------------
144 /* CONSTANT Parameters : */
146 /* To increase the length of the table of paste path lengths relative
147 to the length of the hash table, increase MULT.
150 /* AMISS is a missing value indicator which is returned when the
151 probability is not defined.
153 const double amiss = -12345.;
155 Set IREAL = 4 for DOUBLE PRECISION
156 Set IREAL = 3 for SINGLE PRECISION
161 /* System generated locals */
163 /* Local variables */
164 int nco, nro, ntot, numb, iiwk, irwk;
165 int i, j, k, kk, ldkey, ldstp, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
166 int i3a, i3b, i3c, i9a, iwkmax, iwkpt;
168 /* Workspace Allocation (freed at end) */
170 iwkmax = 2 * (int) (*workspace / 2);
171 // equiv = (double *) R_alloc(iwkmax / 2, sizeof(double));
172 equiv = (double *) calloc(iwkmax / 2, sizeof(double));
174 /* The check could never happen with Calloc!
175 equiv = Calloc(iwkmax / 2, double);
177 prterr(0, "Can not allocate specified workspace");
181 #define iwrk ((int *)equiv)
182 #define rwrk ((float *)equiv)
184 /* Parameter adjustments */
185 table -= *ldtabl + 1;
191 prterr(1, "NROW must be less than or equal to LDTABL.");
194 for (i = 1; i <= *nrow; ++i) {
195 for (j = 1; j <= *ncol; ++j) {
196 if (table[i + j * *ldtabl] < 0.)
197 prterr(2, "All elements of TABLE must be positive.");
198 ntot = (int) (ntot + table[i + j * *ldtabl]);
202 prterr(3, "All elements of TABLE are zero.\n"
203 "PRT and PRE are set to missing values.");
209 nco = max(*nrow, *ncol);
210 nro = *nrow + *ncol - nco;/* = min(*nrow, *ncol) */
211 k = *nrow + *ncol + 1;
215 i1 = iwork(iwkmax, &iwkpt, ikh, i_real);
216 i2 = iwork(iwkmax, &iwkpt, nco, i_int);
217 i3 = iwork(iwkmax, &iwkpt, nco, i_int);
218 i3a = iwork(iwkmax, &iwkpt, nco, i_int);
219 i3b = iwork(iwkmax, &iwkpt, nro, i_int);
220 i3c = iwork(iwkmax, &iwkpt, nro, i_int);
221 ikh = max(k * 5 + (kk << 1), nco * 7 + 800);
222 iiwk= iwork(iwkmax, &iwkpt, ikh, i_int);
223 ikh = max(nco + 401, k);
224 irwk= iwork(iwkmax, &iwkpt, ikh, i_real);
227 What follows below splits the remaining amount iwkmax - iwkpt of
228 (int) workspace into hash tables as follows.
230 INT 2 * ldkey i4 i5 i11
231 REAL 2 * ldkey i8 i9 i10
234 Hence, we need ldkey times
235 3 * 2 + 3 * 2 * s + 2 * mult * s + 6 * mult
236 chunks of integer memory, where s = sizeof(REAL) / sizeof(INT).
237 If doubles are used and are twice as long as ints, this gives
239 so that the value of ldkey can be obtained by dividing available
240 (int) workspace by this number.
242 In fact, because iwork() can actually s * n + s - 1 int chunks
243 when allocating a REAL, we use ldkey = available / numb - 1.
246 Can we always assume that sizeof(double) / sizeof(int) is 2?
249 if (i_real == 4) { /* Double precision reals */
250 numb = 18 + 10 * mult;
251 } else { /* Single precision reals */
252 numb = (mult << 3) + 12;
254 ldkey = (iwkmax - iwkpt) / numb - 1;
255 ldstp = mult * ldkey;
256 ikh = ldkey << 1; i4 = iwork(iwkmax, &iwkpt, ikh, i_int);
257 ikh = ldkey << 1; i5 = iwork(iwkmax, &iwkpt, ikh, i_int);
258 ikh = ldstp << 1; i6 = iwork(iwkmax, &iwkpt, ikh, i_real);
259 ikh = ldstp * 6; i7 = iwork(iwkmax, &iwkpt, ikh, i_int);
260 ikh = ldkey << 1; i8 = iwork(iwkmax, &iwkpt, ikh, i_real);
261 ikh = ldkey << 1; i9 = iwork(iwkmax, &iwkpt, ikh, i_real);
262 ikh = ldkey << 1; i9a = iwork(iwkmax, &iwkpt, ikh, i_real);
263 ikh = ldkey << 1; i10 = iwork(iwkmax, &iwkpt, ikh, i_int);
265 /* To convert to double precision, change RWRK to DWRK in the next CALL.
307 -----------------------------------------------------------------------
309 Purpose: Computes Fisher's exact test for a contingency table,
310 routine with workspace variables specified.
311 Usage: F2XACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
312 EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF,
313 IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ,
314 DLP, DSP, TM, KEY2, IWK, RWK)
315 -----------------------------------------------------------------------
318 f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
319 double *expect, double *percnt, double *emin, double *prt,
320 double *pre, double *fact, int *ico, int *iro, int *kyy,
321 int *idif, int *irn, int *key, int *ldkey, int *ipoin,
322 double *stp, int *ldstp, int *ifrq, double *dlp, double *dsp,
323 double *tm, int *key2, int *iwk, double *rwk)
325 /* IMAX is the largest representable int on the machine. */
326 const int imax = SINT_MAX;
327 // const int imax = 2147483647; //xx: I DON´T like this, and
328 // thanks to the hint from Jason Turner I don't do it anymore. (R.D-U).
330 /* AMISS is a missing value indicator which is returned when the
331 probability is not defined. */
332 const double amiss = -12345.;
334 /* TOL is chosen as the square root of the smallest relative spacing. */
336 const static double tol = 3.45254e-7;
338 static double tol = 3.45254e-7;
340 /* EMX is a large positive value used in comparing expected values. */
341 const static double emx = 1e30;
343 /* Local variables {{any really need to be static ???}} */
344 static int kval, kmax, jkey, last, ipsh, itmp, itop, jstp, ntot,
345 jstp2, jstp3, jstp4, i, ii, j, k, n, iflag, ncell, ifreq, chisq,
346 ikkey, ikstp, ikstp2, k1, kb, kd, ks,
347 i31, i32, i33, i34, i35, i36, i37, i38, i39,
348 i41, i42, i43, i44, i45, i46, i47, i48, i310, i311,
349 nco, nrb, ipn, ipo, itp, nro, nro2;
350 static double dspt, dd, df,ddf, drn,dro, emn, obs, obs2, obs3,
359 /* Parameter adjustments */
360 table -= *ldtabl + 1;
378 /* Check table dimensions */
380 prterr(1, "NROW must be less than or equal to LDTABL.");
382 prterr(4, "NCOL must be at least 2");
384 /* Initialize KEY array */
385 for (i = 1; i <= *ldkey << 1; ++i) {
389 /* Initialize parameters */
402 /* nco := max(nrow, ncol) : */
407 /* Initialize pointers for workspace */
421 k = *nrow + *ncol + 1;
431 /* Compute row marginals and total */
433 for (i = 1; i <= *nrow; ++i) {
435 for (j = 1; j <= *ncol; ++j) {
436 if (table[i + j * *ldtabl] < -1e-4)
437 prterr(2, "All elements of TABLE must be positive.");
438 iro[i] += (int) table[i + j * *ldtabl];
444 prterr(3, "All elements of TABLE are zero.\n"
445 "PRT and PRE are set to missing values.");
450 /* Column marginals */
451 for (i = 1; i <= *ncol; ++i) {
453 for (j = 1; j <= *nrow; ++j)
454 ico[i] += (int) table[j + i * *ldtabl];
458 isort(nrow, &iro[1]);
459 isort(ncol, &ico[1]);
461 /* Determine row and column marginals.
462 Define max(nrow,ncol) =: nco >= nro := min(nrow,ncol)
465 Swap marginals if necessary to ico[1:nco] & iro[1:nro]
470 for (i = 1; i <= nco; ++i) {
480 /* Get multiplers for stack */
482 for (i = 2; i <= nro; ++i) {
483 /* Hash table multipliers */
484 if (iro[i - 1] + 1 <= imax / kyy[i - 1]) {
485 kyy[i] = kyy[i - 1] * (iro[i - 1] + 1);
490 /* Maximum product */
491 if (iro[nro - 1] + 1 <= imax / kyy[nro - 1]) {
492 kmax = (iro[nro] + 1) * kyy[nro - 1];
495 prterr(5, "The hash table key cannot be computed because "
497 "is larger than the largest representable int.\n"
498 "The algorithm cannot proceed.\n"
499 "Reduce the workspace size, or use `exact = FALSE'.");
503 /* Compute log factorials */
506 if(ntot >= 2) fact[2] = log(2.);
507 /* MM: old code assuming log() to be SLOW */
508 for (i = 3; i <= ntot; i += 2) {
509 fact[i] = fact[i - 1] + log((double) i);
512 fact[j] = fact[i] + fact[2] + fact[j / 2] - fact[j / 2 - 1];
514 /* Compute obs := observed path length */
517 for (j = 1; j <= nco; ++j) {
519 for (i = 1; i <= nro; ++i) {
521 dd += fact[(int) table[j + i * *ldtabl]];
522 ntot += (int) table[j + i * *ldtabl];
524 dd += fact[(int) table[i + j * *ldtabl]];
525 ntot += (int) table[i + j * *ldtabl];
528 obs += fact[ico[j]] - dd;
530 /* Denominator of observed table: DRO */
531 dro = f9xact(&nro, &ntot, &iro[1], fact);
532 *prt = exp(obs - dro);
533 /* Initialize pointers */
538 jstp2 = *ldstp * 3 + 1;
539 jstp3 = (*ldstp << 2) + 1;
540 jstp4 = *ldstp * 5 + 1;
543 ikstp2 = *ldstp << 1;
548 ifrq[ikstp2 + 1] = -1;
556 /* IDIF is the difference in going to the daughter */
557 for (i = 1; i <= nro; ++i)
560 /* Generate the first daughter */
563 ntot = min(n, iro[kd]);
569 while (n > 0 && kd != 1);
578 for (i = kb + 1; i <= nco; ++i)
583 /* Arc to daughter length=ICO(KB) */
584 for (i = 1; i <= nro; ++i)
585 irn[i] = iro[i] - idif[i];
590 if (irn[1] > irn[2]) {
595 } else if (nro == 3) {
599 if (irn[2] > irn[3]) {
612 } else if (ii > irn[2]) {
615 } else if (irn[2] > irn[3]) {
621 for (j = 2; j <= nro; ++j) {
625 while (ii < irn[i]) {
634 /* Adjust start for zero */
635 for (i = 1; i <= nro; ++i) {
646 /* Some table values */
647 ddf = f9xact(&nro, &n, &idif[1], fact);
648 drn = f9xact(&nro2, &ntot, &irn[nrb], fact) - dro + ddf;
651 kval = irn[1] + irn[2] * kyy[2];
652 for (i = 3; i <= nro; ++i) {
653 kval += irn[i] * kyy[i];
655 /* Get hash table entry */
656 i = kval % (*ldkey << 1) + 1;
657 /* Search for unused location */
658 for (itp = i; itp <= *ldkey << 1; ++itp) {
670 for (itp = 1; itp <= i - 1; ++itp) {
682 prterr(6, "LDKEY is too small.\n"
683 "It is not possible to give the value of LDKEY required,\n"
684 "but you could try doubling LDKEY (and possibly LDSTP).");
686 prterr(6, "LDKEY is too small for this problem.\n"
687 "Try increasing the size of the workspace.");
693 ipn = ipoin[ipo + ikkey];
694 pastp = stp[ipn + ikstp];
695 ifreq = ifrq[ipn + ikstp];
696 /* Compute shortest and longest path */
698 obs2 = obs - fact[ico[kb + 1]] - fact[ico[kb + 2]] - ddf;
699 for (i = 3; i <= k1; ++i) {
700 obs2 -= fact[ico[kb + i]];
703 dspt = obs - obs2 - ddf;
704 /* Compute longest path */
706 f3xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dlp[itp],
707 &ntot, fact, &iwk[i31], &iwk[i32], &iwk[i33],
708 &iwk[i34], &iwk[i35], &iwk[i36], &iwk[i37],
709 &iwk[i38], &iwk[i39], &rwk[i310], &rwk[i311], &tol);
710 dlp[itp] = min(0., dlp[itp]);
711 /* Compute shortest path */
713 f4xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dsp[itp], fact,
714 &iwk[i47], &iwk[i41], &iwk[i42], &iwk[i43],
715 &iwk[i44], &iwk[i45], &iwk[i46], &rwk[i48], &tol);
716 dsp[itp] = min(0., dsp[itp] - dspt);
717 /* Use chi-squared approximation? */
718 if ((irn[nrb] * ico[kb + 1]) > ntot * emn) {
720 for (i = 0; i < nro2; ++i)
721 for (j = 1; j <= k1; ++j)
722 if (irn[nrb + i] * ico[kb + j] >= ntot * *expect)
725 if (ncell * 100 >= k1 * nro2 * *percnt) {
727 for (i = 0; i < nro2; ++i)
728 tmp += (fact[irn[nrb + i]] -
729 fact[irn[nrb + i] - 1]);
731 for (j = 1; j <= k1; ++j)
732 tmp += (nro2 - 1) * (fact[ico[kb + j]] -
733 fact[ico[kb + j] - 1]);
734 df = (double) ((nro2 - 1) * (k1 - 1));
735 tmp += df * 1.83787706640934548356065947281;
736 tmp -= (nro2 * k1 - 1) * (fact[ntot] - fact[ntot - 1]);
737 tm[itp] = (obs - dro) * -2. - tmp;
739 /* tm(itp) set to a flag value */
746 obs3 = obs2 - dlp[itp];
748 if (tm[itp] == -9876.) {
755 obs2 = obs - drn - dro;
760 /* Process node with new PASTP */
763 *pre += (double) ifreq * exp(pastp + drn);
764 } else if (pastp < obs2) {
766 df = (double) ((nro2 - 1) * (k1 - 1));
768 pv = pgamma(fmax2(0., tmp + (pastp + drn) * 2.) / 2.,
769 df / 2., /*scale = */ 1.,
770 /*lower_tail = */FALSE, /*log_p = */ FALSE);
772 d1 = max(0., tmp + (pastp + drn) * 2.) / 2.;
774 pv = 1. - gammds(&d1, &d2, &ifault);
776 *pre += (double) ifreq * exp(pastp + drn) * pv;
778 /* Put daughter on queue */
780 f5xact(&d1, &tol, &kval, &key[jkey], ldkey, &ipoin[jkey],
781 &stp[jstp], ldstp, &ifrq[jstp], &ifrq[jstp2],
782 &ifrq[jstp3], &ifrq[jstp4], &ifreq, &itop, &ipsh);
786 /* Get next PASTP on chain */
787 ipn = ifrq[ipn + ikstp2];
789 pastp = stp[ipn + ikstp];
790 ifreq = ifrq[ipn + ikstp];
793 /* Generate a new daughter node */
794 f7xact(&kmax, &iro[1], &idif[1], &kd, &ks, &iflag);
800 /* Go get a new mother from stage K */
803 f6xact(&nro, &iro[1], &iflag, &kyy[1], &key[ikkey + 1], ldkey,
805 /* Update pointers */
808 /* else iflag == 3 : no additional nodes to process */
814 jkey = *ldkey - jkey + 2;
815 jstp = *ldstp - jstp + 2;
816 jstp2 = (*ldstp << 1) + jstp;
817 for (i = 1; i <= *ldkey << 1; ++i)
824 -----------------------------------------------------------------------
826 Purpose: Computes the shortest path length for a given table.
827 Usage: F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, IRO,
828 IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, TOL)
830 NROW - The number of rows in the table. (Input)
831 IROW - Vector of length NROW containing the row sums
832 for the table. (Input)
833 NCOL - The number of columns in the table. (Input)
834 ICOL - Vector of length K containing the column sums
835 for the table. (Input)
836 DLP - The longest path for the table. (Output)
837 MM - The total count in the table. (Output)
838 FACT - Vector containing the logarithms of factorials. (Input)
839 ICO - Work vector of length MAX(NROW,NCOL).
840 IRO - Work vector of length MAX(NROW,NCOL).
841 IT - Work vector of length MAX(NROW,NCOL).
842 LB - Work vector of length MAX(NROW,NCOL).
843 NR - Work vector of length MAX(NROW,NCOL).
844 NT - Work vector of length MAX(NROW,NCOL).
845 NU - Work vector of length MAX(NROW,NCOL).
846 ITC - Work vector of length 400.
847 IST - Work vector of length 400.
848 STV - Work vector of length 400.
849 ALEN - Work vector of length MAX(NROW,NCOL).
850 TOL - Tolerance. (Input)
851 -----------------------------------------------------------------------
855 f3xact(int *nrow, int *irow, int *ncol, int *icol, double *dlp,
856 int *mm, double *fact, int *ico, int *iro, int *it,
857 int *lb, int *nr, int *nt, int *nu, int *itc, int *ist,
858 double *stv, double *alen, const double *tol)
860 /* Initialized data */
861 static int ldst = 200;
865 /* Local variables */
869 static int n11, n12, ii, nn, ks, ic1, ic2, nc1, nn1;
872 static int nct, ipn, irl, key, lev, itp, nro;
874 static int nrt, kyy, nc1s;
876 /* Parameter adjustments */
891 for (i = 0; i <= *ncol; ++i) {
894 for (i = 1; i <= 400; ++i) {
900 *dlp -= fact[icol[1]];
901 for (i = 2; i <= *ncol; ++i) {
902 *dlp -= fact[icol[i]];
910 *dlp = *dlp - fact[irow[1]] - fact[irow[2]];
911 for (i = 3; i <= *nrow; ++i) {
912 *dlp -= fact[irow[i]];
918 if (*nrow * *ncol == 4) {
919 n11 = (irow[1] + 1) * (icol[1] + 1) / (*mm + 2);
921 *dlp = *dlp - fact[n11] - fact[n12] - fact[icol[1] - n11]
922 - fact[icol[2] - n12];
925 /* Test for optimal table */
928 if (irow[*nrow] <= irow[1] + *ncol) {
929 f10act(nrow, &irow[1], ncol, &icol[1], &val, &xmin, fact,
930 &lb[1], &nu[1], &nr[1]);
933 if (icol[*ncol] <= icol[1] + *nrow) {
934 f10act(ncol, &icol[1], nrow, &irow[1], &val, &xmin, fact,
935 &lb[1], &nu[1], &nr[1]);
943 /* Setup for dynamic programming */
946 if (*nrow >= *ncol) {
949 for (i = 1; i <= *nrow; ++i) {
954 for (i = 2; i <= *ncol; ++i) {
956 nt[i] = nt[i - 1] - ico[i];
963 for (i = 2; i <= *nrow; ++i) {
965 nt[i] = nt[i - 1] - ico[i];
967 for (i = 1; i <= *ncol; ++i)
970 /* Initialize pointers */
978 LnewNode: /* Setup to generate new node */
984 lb[1] = (int) ((double) ((nrt + 1) * (nct + 1)) /
985 (double) (nn + nr1 * nc1s + 1) - *tol) - 1;
986 nu[1] = (int) ((double) ((nrt + nc1s) * (nct + nr1)) /
987 (double) (nn + nr1 + nc1s)) - lb[1] + 1;
990 LoopNode: /* Generate a node */
1002 alen[lev] = alen[lev - 1] + fact[lb[lev]];
1009 lb[lev] = (int) ((double) ((nrt + 1) * (nct + 1)) /
1010 (double) (nn1 + nr1 * nc1 + 1) - *tol);
1011 nu[lev] = (int) ((double) ((nrt + nc1) * (nct + nr1)) /
1012 (double) (nn1 + nr1 + nc1) - lb[lev] + 1);
1013 nr[lev] = nrt - lb[lev];
1016 alen[nco] = alen[lev] + fact[nr[lev]];
1019 v = val + alen[nco];
1021 /* Only 1 row left */
1022 v = v + fact[ico[1] - lb[1]] + fact[ico[2] - lb[2]];
1023 for (i = 3; i <= nco; ++i) {
1024 v += fact[ico[i] - lb[i]];
1029 } else if (nro == 3 && nco == 2) {
1030 /* 3 rows and 2 columns */
1031 nn1 = nn - iro[irl] + 2;
1032 ic1 = ico[1] - lb[1];
1033 ic2 = ico[2] - lb[2];
1034 n11 = (iro[irl + 1] + 1) * (ic1 + 1) / nn1;
1035 n12 = iro[irl + 1] - n11;
1036 v = v + fact[n11] + fact[n12] + fact[ic1 - n11]
1042 /* Column marginals are new node */
1043 for (i = 1; i <= nco; ++i) {
1044 it[i] = ico[i] - lb[i];
1046 /* Sort column marginals */
1048 if (it[1] > it[2]) {
1053 } else if (nco == 3) {
1057 if (it[2] > it[3]) {
1070 } else if (ii > it[2]) {
1073 } else if (it[2] > it[3]) {
1079 isort(&nco, &it[1]);
1081 /* Compute hash value */
1082 key = it[1] * kyy + it[2];
1083 for (i = 3; i <= nco; ++i) {
1084 key = it[i] + key * kyy;
1087 //PROBLEM "Bug in FEXACT: gave negative key" RECOVER(NULL_ENTRY);
1088 printf("Bug in FEXACT: gave negative key \n"); //xx:another one of my ugly kluges (R.D-U)
1091 ipn = key % ldst + 1;
1093 /* Find empty position */
1094 for (itp = ipn, ii = ks + ipn; itp <= ldst; ++itp, ++ii) {
1097 } else if (ist[ii] == key) {
1102 for (itp = 1, ii = ks + 1; itp <= ipn - 1; ++itp, ++ii) {
1105 } else if (ist[ii] == key) {
1110 prterr(30, "Stack length exceeded in f3xact.\n"
1111 "This problem should not occur.");
1113 L180: /* Push onto stack */
1121 L190: /* Marginals already on stack */
1122 stv[ii] = min(v, stv[ii]);
1127 L200: /* Pop item from stack */
1130 itp = itc[nitc + k] + k;
1135 /* Compute marginals */
1136 for (i = nco; i >= 2; --i) {
1141 /* Set up nt array */
1142 nt[1] = nn - ico[1];
1143 for (i = 2; i <= nco; ++i)
1144 nt[i] = nt[i - 1] - ico[i];
1146 /* Test for optimality (L90) */
1148 if (iro[nro] <= iro[irl] + nco) {
1149 f10act(&nro, &iro[irl], &nco, &ico[1], &val, &xmin, fact,
1150 &lb[1], &nu[1], &nr[1]);
1152 if (!xmin && ico[nco] <= ico[1] + nro)
1153 f10act(&nco, &ico[1], &nro, &iro[irl], &val, &xmin, fact,
1154 &lb[1], &nu[1], &nr[1]);
1162 } else if (nro > 2 && nst > 0) {
1163 /* Go to next level */
1178 -----------------------------------------------------------------------
1180 Purpose: Computes the longest path length for a given table.
1181 Usage: CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK,
1182 NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK,
1185 NROW - The number of rows in the table. (Input)
1186 IROW - Vector of length NROW containing the row sums for the
1188 NCOL - The number of columns in the table. (Input)
1189 ICOL - Vector of length K containing the column sums for the
1191 DSP - The shortest path for the table. (Output)
1192 FACT - Vector containing the logarithms of factorials. (Input)
1193 ICSTK - NCOL by NROW+NCOL+1 work array.
1194 NCSTK - Work vector of length NROW+NCOL+1.
1195 LSTK - Work vector of length NROW+NCOL+1.
1196 MSTK - Work vector of length NROW+NCOL+1.
1197 NSTK - Work vector of length NROW+NCOL+1.
1198 NRSTK - Work vector of length NROW+NCOL+1.
1199 IRSTK - NROW by MAX(NROW,NCOL) work array.
1200 YSTK - Work vector of length NROW+NCOL+1.
1201 TOL - Tolerance. (Input)
1202 -----------------------------------------------------------------------
1206 f4xact(int *nrow, int *irow, int *ncol, int *icol, double *dsp,
1207 double *fact, int *icstk, int *ncstk, int *lstk, int *mstk,
1208 int *nstk, int *nrstk, int *irstk, double *ystk, const double *tol)
1210 /* System generated locals */
1213 /* Local variables */
1214 int i, j, k, l, m, n, mn, ic1, ir1, ict, irt, istk, nco, nro;
1217 /* Parameter adjustments */
1230 /* Take care of the easy cases first */
1232 for (i = 1; i <= *ncol; ++i) {
1233 *dsp -= fact[icol[i]];
1238 for (i = 1; i <= *nrow; ++i) {
1239 *dsp -= fact[irow[i]];
1243 if (*nrow * *ncol == 4) {
1244 if (irow[2] <= icol[2]) {
1245 *dsp = *dsp - fact[irow[2]] - fact[icol[1]]
1246 - fact[icol[2] - irow[2]];
1248 *dsp = *dsp - fact[icol[2]] - fact[irow[1]]
1249 - fact[irow[2] - icol[2]];
1253 /* initialization before loop */
1254 for (i = 1; i <= *nrow; ++i) {
1255 irstk[i + *nrow] = irow[*nrow - i + 1];
1257 for (j = 1; j <= *ncol; ++j) {
1258 icstk[j + *ncol] = icol[*ncol - j + 1];
1273 ir1 = irstk[istk * *nrow + 1];
1274 ic1 = icstk[istk * *ncol + 1];
1281 } else if (ir1 < ic1) {
1302 irt = irstk[i + istk * *nrow];
1303 ict = icstk[j + istk * *ncol];
1312 f11act(&irstk[istk * *nrow + 1], &i, &nro,
1313 &irstk[(istk + 1) * *nrow + 1]);
1314 f11act(&icstk[istk * *ncol + 1], &j, &nco,
1315 &icstk[(istk + 1) * *ncol + 1]);
1316 } else if (irt > ict) {
1318 f11act(&icstk[istk * *ncol + 1], &j, &nco,
1319 &icstk[(istk + 1) * *ncol + 1]);
1321 f8xact(&irstk[istk * *nrow + 1], &ikh, &i,
1322 &nro, &irstk[(istk + 1) * *nrow + 1]);
1325 f11act(&irstk[istk * *nrow + 1], &i, &nro,
1326 &irstk[(istk + 1) * *nrow + 1]);
1328 f8xact(&icstk[istk * *ncol + 1], &ikh, &j,
1329 &nco, &icstk[(istk + 1) * *ncol + 1]);
1333 for (k = 1; k <= nco; ++k) {
1334 y += fact[icstk[k + (istk + 1) * *ncol]];
1339 for (k = 1; k <= nro; ++k) {
1340 y += fact[irstk[k + (istk + 1) * *nrow]];
1353 } while(1);/* end do */
1358 if (*dsp - amx <= *tol) {
1368 if (*dsp - amx <= *tol) {
1377 if (l > mstk[istk]) goto L100;
1384 if (irstk[l + istk * *nrow] <
1385 irstk[l - 1 + istk * *nrow]) goto L60;
1388 if (icstk[l + istk * *ncol] <
1389 icstk[l - 1 + istk * *ncol]) goto L60;
1395 -----------------------------------------------------------------------
1397 Purpose: Put node on stack in network algorithm.
1398 Usage: CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP,
1399 LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP,
1402 PASTP - The past path length. (Input)
1403 TOL - Tolerance for equivalence of past path lengths. (Input)
1404 KVAL - Key value. (Input)
1405 KEY - Vector of length LDKEY containing the key values. (in/out)
1406 LDKEY - Length of vector KEY. (Input)
1407 IPOIN - Vector of length LDKEY pointing to the
1408 linked list of past path lengths. (in/out)
1409 STP - Vector of length LSDTP containing the
1410 linked lists of past path lengths. (in/out)
1411 LDSTP - Length of vector STP. (Input)
1412 IFRQ - Vector of length LDSTP containing the past path
1413 frequencies. (in/out)
1414 NPOIN - Vector of length LDSTP containing the pointers to
1415 the next past path length. (in/out)
1416 NR - Vector of length LDSTP containing the right object
1417 pointers in the tree of past path lengths. (in/out)
1418 NL - Vector of length LDSTP containing the left object
1419 pointers in the tree of past path lengths. (in/out)
1420 IFREQ - Frequency of the current path length. (Input)
1421 ITOP - Pointer to the top of STP. (Input)
1422 IPSH - Option parameter. (Input)
1423 If IPSH is true, the past path length is found in the
1424 table KEY. Otherwise the location of the past path
1425 length is assumed known and to have been found in
1426 a previous call. ==>>>>> USING "static" variables
1427 -----------------------------------------------------------------------
1431 f5xact(double *pastp, const double *tol, int *kval, int *key, int *ldkey,
1432 int *ipoin, double *stp, int *ldstp, int *ifrq, int *npoin,
1433 int *nr, int *nl, int *ifreq, int *itop, int *ipsh)
1435 /* Local variables */
1436 static int itmp, ird, ipn, itp;
1437 double test1, test2;
1439 /* Parameter adjustments */
1450 /* Convert KVAL to int in range 1, ..., LDKEY. */
1451 ird = *kval % *ldkey + 1;
1452 /* Search for an unused location */
1453 for (itp = ird; itp <= *ldkey; ++itp) {
1454 if (key[itp] == *kval) {
1461 for (itp = 1; itp <= ird - 1; ++itp) {
1462 if (key[itp] == *kval) {
1469 /* Return if KEY array is full */
1471 prterr(6, "LDKEY is too small for this problem.\n"
1472 "It is not possible to estimate the value of LDKEY "
1474 "but twice the current value may be sufficient.");
1476 prterr(6, "LDKEY is too small for this problem.\n"
1477 "Try increasing the size of the workspace.");
1484 /* Return if STP array full */
1485 if (*itop > *ldstp) {
1487 prterr(7, "LDSTP is too small for this problem.\n"
1488 "It is not possible to estimate the value of LDSTP "
1490 "but twice the current value may be sufficient.");
1492 prterr(7, "LDSTP is too small for this problem.\n"
1493 "Try increasing the size of the workspace.");
1495 /* Update STP, etc. */
1499 stp[*itop] = *pastp;
1500 ifrq[*itop] = *ifreq;
1504 /* Find location, if any, of pastp */
1507 test1 = *pastp - *tol;
1508 test2 = *pastp + *tol;
1511 if (stp[ipn] < test1) {
1516 } else if (stp[ipn] > test2) {
1522 ifrq[ipn] += *ifreq;
1525 /* Return if STP array full */
1527 if (*itop > *ldstp) {
1529 prterr(7, "LDSTP is too small for this problem.\n"
1530 "It is not possible to estimate the value of LDSTP "
1532 "but twice the current value may be sufficient.");
1534 prterr(7, "LDSTP is too small for this problem.\n"
1535 "Try increasing the size of the workspace.");
1538 /* Find location to add value */
1543 if (stp[ipn] < test1) {
1551 } else if (stp[ipn] > test2) {
1560 /* Update STP, etc. */
1561 npoin[*itop] = npoin[itmp];
1562 npoin[itmp] = *itop;
1563 stp[*itop] = *pastp;
1564 ifrq[*itop] = *ifreq;
1570 -----------------------------------------------------------------------
1572 Purpose: Pop a node off the stack.
1573 Usage: CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, IPN)
1575 NROW - The number of rows in the table. (Input)
1576 IROW - Vector of length nrow containing the row sums on
1578 IFLAG - Set to 3 if there are no additional nodes to process.
1580 KYY - Constant mutlipliers used in forming the hash
1582 KEY - Vector of length LDKEY containing the hash table
1584 LDKEY - Length of vector KEY. (Input)
1585 LAST - Index of the last key popped off the stack. (In/out)
1586 IPN - Pointer to the linked list of past path lengths. (Output)
1587 -----------------------------------------------------------------------
1590 f6xact(int *nrow, int *irow, int *iflag, int *kyy, int *key, int
1591 *ldkey, int *last, int *ipn)
1595 /* Parameter adjustments */
1603 if (*last <= *ldkey) {
1604 if (key[*last] < 0) {
1607 /* Get KVAL from the stack */
1610 for (j = *nrow; j >= 2; --j) {
1611 irow[j] = kval / kyy[j];
1612 kval -= irow[j] * kyy[j];
1624 -----------------------------------------------------------------------
1626 Purpose: Generate the new nodes for given marinal totals.
1627 Usage: CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG)
1629 NROW - The number of rows in the table. (Input)
1630 IMAX - The row marginal totals. (Input)
1631 IDIF - The column counts for the new column. (in/out)
1632 K - Indicator for the row to decrement. (in/out)
1633 KS - Indicator for the row to increment. (in/out)
1634 IFLAG - Status indicator. (Output)
1635 If IFLAG is zero, a new table was generated. For
1636 IFLAG = 1, no additional tables could be generated.
1637 -----------------------------------------------------------------------
1641 f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
1647 /* Parameter adjustments */
1653 /* Find node which can be incremented, ks */
1657 } while (idif[*ks] == imax[*ks]);
1659 /* Find node to decrement (>ks) */
1660 if (idif[*k] > 0 && *k > *ks) {
1664 } while(imax[*k] == 0);
1668 /* Find node to increment (>=ks) */
1669 while (idif[m] >= imax[m]) {
1675 if (idif[m] == imax[m]) {
1682 /* Check for finish */
1683 for (k1 = *k + 1; k1 <= *nrow; ++k1) {
1692 /* Reallocate counts */
1694 for (i = 1; i <= *k; ++i) {
1702 m = min(mm, imax[*k]);
1705 } while (mm > 0 && *k != 1);
1707 /* Check that all counts reallocated */
1724 } while (idif[*ks] >= imax[*ks]);
1729 -----------------------------------------------------------------------
1731 Purpose: Routine for reducing a vector when there is a zero
1733 Usage: CALL F8XACT (IROW, IS, I1, IZERO, NEW)
1735 IROW - Vector containing the row counts. (Input)
1736 IS - Indicator. (Input)
1737 I1 - Indicator. (Input)
1738 IZERO - Position of the zero. (Input)
1739 NEW - Vector of new row counts. (Output)
1740 -----------------------------------------------------------------------
1744 f8xact(int *irow, int *is, int *i1, int *izero, int *new)
1748 /* Parameter adjustments */
1753 for (i = 1; i < *i1; ++i)
1756 for (i = *i1; i <= *izero - 1; ++i) {
1757 if (*is >= irow[i + 1])
1759 new[i] = irow[i + 1];
1773 -----------------------------------------------------------------------
1775 Purpose: Computes the log of a multinomial coefficient.
1776 Usage: F9XACT(N, MM, IR, FACT)
1778 N - Length of IR. (Input)
1779 MM - Number for factorial in numerator. (Input)
1780 IR - Vector of length N containing the numbers for
1781 the denominator of the factorial. (Input)
1782 FACT - Table of log factorials. (Input)
1783 F9XACT - The log of the multinomal coefficient. (Output)
1784 -----------------------------------------------------------------------
1788 f9xact(int *n, int *mm, int *ir, double *fact)
1794 for (k = 0; k < *n; k++)
1800 -----------------------------------------------------------------------
1802 Purpose: Computes the shortest path length for special tables.
1803 Usage: F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, NE, M)
1805 NROW - The number of rows in the table. (Input)
1806 IROW - Vector of length NROW containing the row totals. (Input)
1807 NCOL - The number of columns in the table. (Input)
1808 ICO - Vector of length NCOL containing the column totals.(Input)
1809 VAL - The shortest path. (Output)
1810 XMIN - Set to true if shortest path obtained. (Output)
1811 FACT - Vector containing the logarithms of factorials. (Input)
1812 ND - Workspace vector of length NROW. (Input)
1813 NE - Workspace vector of length NCOL. (Input)
1814 M - Workspace vector of length NCOL. (Input)
1816 Chapter: STAT/LIBRARY Categorical and Discrete Data Analysis
1817 -----------------------------------------------------------------------
1821 f10act(int *nrow, int *irow, int *ncol, int *icol, double *val,
1822 int *xmin, double *fact, int *nd, int *ne, int *m)
1824 /* Local variables */
1825 int i, is, ix, nrw1;
1827 /* Parameter adjustments */
1835 for (i = 1; i <= *nrow - 1; ++i)
1838 is = icol[1] / *nrow;
1839 ix = icol[1] - *nrow * is;
1845 for (i = 2; i <= *ncol; ++i) {
1846 ix = icol[i] / *nrow;
1849 ix = icol[i] - *nrow * ix;
1855 for (i = *nrow - 2; i >= 1; --i)
1860 for (i = *nrow; i >= 2; --i) {
1861 ix = ix + is + nd[nrw1 - i] - irow[i];
1866 for (i = 1; i <= *ncol; ++i) {
1869 *val = *val + is * fact[ix + 1] + (*nrow - is) * fact[ix];
1877 -----------------------------------------------------------------------
1879 Purpose: Routine for revising row totals.
1880 Usage: CALL F11ACT (IROW, I1, I2, NEW)
1882 IROW - Vector containing the row totals. (Input)
1883 I1 - Indicator. (Input)
1884 I2 - Indicator. (Input)
1885 NEW - Vector containing the row totals. (Output)
1886 -----------------------------------------------------------------------
1889 f11act(int *irow, int *i1, int *i2, int *new)
1893 /* Parameter adjustments */
1897 for (i = 1; i <= (*i1 - 1); ++i) new[i] = irow[i];
1898 for (i = *i1; i <= *i2; ++i) new[i] = irow[i + 1];
1904 -----------------------------------------------------------------------
1906 Purpose: Print an error message and stop.
1907 Usage: prterr(icode, mes)
1909 icode - Integer code for the error message. (Input)
1910 mes - Character string containing the error message. (Input)
1911 -----------------------------------------------------------------------
1914 prterr(int icode, char *mes)
1916 // PROBLEM "FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY);
1917 // printf("FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY));
1918 printf("FEXACT error %d.\n", icode); //xx:another one of my ugly kluges
1923 -----------------------------------------------------------------------
1925 Purpose: Routine for allocating workspace.
1926 Usage: iwork (iwkmax, iwkpt, number, itype)
1928 iwkmax - Maximum (int) amount of workspace. (Input)
1929 iwkpt - Amount of (int) workspace currently allocated. (in/out)
1930 number - Number of elements of workspace desired. (Input)
1931 itype - Workspace type. (Input)
1936 iwork(): Index in rwrk, dwrk, or iwrk of the beginning of
1937 the first free element in the workspace array. (Output)
1938 -----------------------------------------------------------------------
1941 iwork(int iwkmax, int *iwkpt, int number, int itype)
1946 if (itype == 2 || itype == 3)
1951 *iwkpt += (number << 1);
1954 if (*iwkpt > iwkmax)
1955 prterr(40, "Out of workspace.");
1962 void isort(int *n, int *ix)
1965 -----------------------------------------------------------------------
1967 Purpose: Shell sort for an int vector.
1968 Usage: CALL ISORT (N, IX)
1970 N - Lenth of vector IX. (Input)
1971 IX - Vector to be sorted. (in/out)
1972 -----------------------------------------------------------------------
1974 static int ikey, i, j, m, il[10], kl, it, iu[10], ku;
1976 /* Parameter adjustments */
1992 /* Find element in first half */
1996 if (ix[ikey] > ix[i]) {
2000 /* Find element in second half */
2003 if (ix[j] > ix[ikey]) {
2016 /* Save upper and lower subscripts of the array yet to be sorted */
2018 if (j - kl < ku - j) {
2032 prterr(20, "This should never occur.");
2034 /* Use another segment */
2045 double gammds(double *y, double *p, int *ifault)
2048 -----------------------------------------------------------------------
2050 Purpose: Cumulative distribution for the gamma distribution.
2051 Usage: PGAMMA (Q, ALPHA,IFAULT)
2053 Q - Value at which the distribution is desired. (Input)
2054 ALPHA - Parameter in the gamma distribution. (Input)
2055 IFAULT - Error indicator. (Output)
2058 1 An argument is misspecified.
2059 2 A numerical error has occurred.
2060 PGAMMA - The cdf for the gamma distribution with parameter alpha
2061 evaluated at Q. (Output)
2062 -----------------------------------------------------------------------
2064 Algorithm AS 147 APPL. Statist. (1980) VOL. 29, P. 113
2066 Computes the incomplete gamma integral for positive parameters Y, P
2067 using and infinite series.
2070 static double a, c, f, g;
2073 /* Checks for the admissibility of arguments and value of F */
2076 if (*y <= 0. || *p <= 0.) {
2082 ALOGAM is natural log of gamma function no need to test ifail as
2083 an error is impossible
2087 f = exp(*p * log(*y) - alogam(&a, &ifail) - *y);
2109 -----------------------------------------------------------------------
2111 Purpose: Value of the log-gamma function.
2112 Usage: ALOGAM (X, IFAULT)
2114 X - Value at which the log-gamma function is to be evaluated.
2116 IFAULT - Error indicator. (Output)
2120 ALGAMA - The value of the log-gamma function at XX. (Output)
2121 -----------------------------------------------------------------------
2123 Algorithm ACM 291, Comm. ACM. (1966) Vol. 9, P. 684
2125 Evaluates natural logarithm of gamma(x) for X greater than zero.
2128 double alogam(double *x, int *ifault)
2130 /* Initialized data */
2132 static double a1 = .918938533204673;
2133 static double a2 = 5.95238095238e-4;
2134 static double a3 = 7.93650793651e-4;
2135 static double a4 = .002777777777778;
2136 static double a5 = .083333333333333;
2138 /* Local variables */
2139 static double f, y, z;
2163 return(f + (y - .5) * log(y) - y + a1 +
2164 (((-a2 * z + a3) * z - a4) * z + a5) / y);
2167 #endif /* not USING_R */