fisher2.c

   1 #include "fisher2.h"
   2
   3
   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,
  15                   const double *tol);
  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,
  27                    int *iflag);
  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,
  32                   int *ne, int *m);
  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);
  42
  43
  44 /* The only public function : */
  45 void fexact(int *nrow, int *ncol, double *table, int *ldtabl,
  46        double *expect, double *percnt, double *emin, double *prt,
  47        double *pre, /* myNew in C : */ int *workspace) {
  48
  49 /*
  50   ALGORITHM 643, COLLECTED ALGORITHMS FROM ACM.
  51   THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE,
  52   VOL. 19, NO. 4, DECEMBER, 1993, PP. 484-488.
  53   -----------------------------------------------------------------------
  54   Name:       FEXACT
  55   Purpose:    Computes Fisher's exact test probabilities and a hybrid
  56               approximation to Fisher exact test probabilities for a
  57               contingency table using the network algorithm.
  58   Usage:      CALL FEXACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
  59                            EMIN, PRT, PRE)
  60   Arguments:
  61     NROW    - The number of rows in the table.                  (Input)
  62     NCOL    - The number of columns in the table.               (Input)
  63     TABLE   - NROW by NCOL matrix containing the contingency
  64               table.                                            (Input)
  65     LDTABL  - Leading dimension of TABLE exactly as specified
  66               in the dimension statement in the calling
  67               program.                                          (Input)
  68     EXPECT  - Expected value used in the hybrid algorithm for
  69               deciding when to use asymptotic theory
  70               probabilities.                                    (Input)
  71               If EXPECT <= 0.0 then asymptotic theory probabilities
  72               are not used and Fisher exact test probabilities are
  73               computed.  Otherwise, if PERCNT or more of the cells in
  74               the remaining table have estimated expected values of
  75               EXPECT or more, with no remaining cell having expected
  76               value less than EMIN, then asymptotic chi-squared
  77               probabilities are used.  See the algorithm section of the
  78               manual document for details.
  79               Use EXPECT = 5.0 to obtain the 'Cochran' condition.
  80     PERCNT  - Percentage of remaining cells that must have
  81               estimated expected values greater than EXPECT
  82               before asymptotic probabilities can be used.      (Input)
  83               See argument EXPECT for details.
  84               Use PERCNT = 80.0 to obtain the 'Cochran' condition.
  85     EMIN    - Minimum cell estimated expected value allowed for
  86               asymptotic chi-squared probabilities to be used.  (Input)
  87               See argument EXPECT for details.
  88               Use EMIN = 1.0 to obtain the 'Cochran' condition.
  89     PRT     - Probability of the observed table for fixed
  90               marginal totals.                                  (Output)
  91     PRE     - Table p-value.                                    (Output)
  92               PRE is the probability of a more extreme table,
  93               where `extreme' is in a probabilistic sense.
  94               If EXPECT < 0 then the Fisher exact probability
  95               is returned.  Otherwise, an approximation to the
  96               Fisher exact probability is computed based upon
  97               asymptotic chi-squared probabilities for ``large''
  98               table expected values.  The user defines ``large''
  99               through the arguments EXPECT, PERCNT, and EMIN.
 100
 101   Remarks:
 102   1. For many problems one megabyte or more of workspace can be
 103      required.  If the environment supports it, the user should begin
 104      by increasing the workspace used to 200,000 units.
 105   2. In FEXACT, LDSTP = 30*LDKEY.  The proportion of table space used
 106      by STP may be changed by changing the line MULT = 30 below to
 107      another value.
 108   3. FEXACT may be converted to single precision by setting IREAL = 3,
 109      and converting all DOUBLE PRECISION specifications (except the
 110      specifications for RWRK, IWRK, and DWRK) to REAL.  This will
 111      require changing the names and specifications of the intrinsic
 112      functions ALOG, AMAX1, AMIN1, EXP, and REAL.  In addition, the
 113      machine specific constants will need to be changed, and the name
 114      DWRK will need to be changed to RWRK in the call to F2XACT.
 115   4. Machine specific constants are specified and documented in F2XACT.
 116      A missing value code is specified in both FEXACT and F2XACT.
 117   5. Although not a restriction, is is not generally practical to call
 118      this routine with large tables which are not sparse and in
 119      which the 'hybrid' algorithm has little effect.  For example,
 120      although it is feasible to compute exact probabilities for the
 121      table
 122             1 8 5 4 4 2 2
 123             5 3 3 4 3 1 0
 124            10 1 4 0 0 0 0,
 125      computing exact probabilities for a similar table which has been
 126      enlarged by the addition of an extra row (or column) may not be
 127      feasible.
 128   -----------------------------------------------------------------------
 129   */
 130
 131     /* CONSTANT Parameters : */
 132
 133     /* To increase the length of the table of paste path lengths relative
 134        to the length of the hash table, increase MULT.
 135     */
 136     const int mult = 30;
 137     /* AMISS is a missing value indicator which is returned when the
 138        probability is not defined.
 139     */
 140     const double amiss = -12345.;
 141     /*
 142       Set IREAL = 4 for DOUBLE PRECISION
 143       Set IREAL = 3 for SINGLE PRECISION
 144     */
 145 #define i_real 4
 146 #define i_int  2
 147
 148     /* System generated locals */
 149     int ikh;
 150     /* Local variables */
 151     int nco, nro, ntot, numb, iiwk, irwk;
 152     int i, j, k, kk, ldkey, ldstp, i1, i2, i3, i4, i5, i6, i7, i8, i9, i10;
 153     int i3a, i3b, i3c, i9a, iwkmax, iwkpt;
 154
 155     /* Workspace Allocation (freed at end) */
 156     double *equiv;
 157     iwkmax = 2 * (int) (*workspace / 2);
 158 //    equiv = (double *) R_alloc(iwkmax / 2, sizeof(double));
 159     equiv = (double *) calloc(iwkmax / 2, sizeof(double));
 160
 161     /* The check could never happen with Calloc!
 162     equiv = Calloc(iwkmax / 2, double);
 163     if (!equiv) {
 164         prterr(0, "Can not allocate specified workspace");
 165     } */
 166
 167 #define dwrk (equiv)
 168 #define iwrk ((int *)equiv)
 169 #define rwrk ((float *)equiv)
 170
 171     /* Parameter adjustments */
 172     table -= *ldtabl + 1;
 173
 174     /* Function Body */
 175     iwkpt = 0;
 176
 177     if (*nrow > *ldtabl)
 178         prterr(1, "NROW must be less than or equal to LDTABL.");
 179
 180     ntot = 0;
 181     for (i = 1; i <= *nrow; ++i) {
 182         for (j = 1; j <= *ncol; ++j) {
 183             if (table[i + j * *ldtabl] < 0.)
 184                 prterr(2, "All elements of TABLE must be positive.");
 185             ntot = (int) (ntot + table[i + j * *ldtabl]);
 186         }
 187     }
 188     if (ntot == 0) {
 189         prterr(3, "All elements of TABLE are zero.\n"
 190                "PRT and PRE are set to missing values.");
 191         *prt = amiss;
 192         *pre = amiss;
 193         goto L_End;
 194     }
 195
 196     nco = max(*nrow, *ncol);
 197     nro = *nrow + *ncol - nco;/* = min(*nrow, *ncol) */
 198     k = *nrow + *ncol + 1;
 199     kk = k * nco;
 200
 201     ikh = ntot + 1;
 202     i1  = iwork(iwkmax, &iwkpt, ikh, i_real);
 203     i2  = iwork(iwkmax, &iwkpt, nco, i_int);
 204     i3  = iwork(iwkmax, &iwkpt, nco, i_int);
 205     i3a = iwork(iwkmax, &iwkpt, nco, i_int);
 206     i3b = iwork(iwkmax, &iwkpt, nro, i_int);
 207     i3c = iwork(iwkmax, &iwkpt, nro, i_int);
 208     ikh = max(k * 5 + (kk << 1), nco * 7 + 800);
 209     iiwk= iwork(iwkmax, &iwkpt, ikh, i_int);
 210     ikh = max(nco + 401, k);
 211     irwk= iwork(iwkmax, &iwkpt, ikh, i_real);
 212
 213     /* NOTE:
 214        What follows below splits the remaining amount iwkmax - iwkpt of
 215        (int) workspace into hash tables as follows.
 216            type  size       index
 217            INT   2 * ldkey  i4 i5 i11
 218            REAL  2 * ldkey  i8 i9 i10
 219            REAL  2 * ldstp  i6
 220            INT   6 * ldstp  i7
 221        Hence, we need ldkey times
 222            3 * 2 + 3 * 2 * s + 2 * mult * s + 6 * mult
 223        chunks of integer memory, where s = sizeof(REAL) / sizeof(INT).
 224        If doubles are used and are twice as long as ints, this gives
 225            18 + 10 * mult
 226        so that the value of ldkey can be obtained by dividing available
 227        (int) workspace by this number.
 228
 229        In fact, because iwork() can actually s * n + s - 1 int chunks
 230        when allocating a REAL, we use ldkey = available / numb - 1.
 231
 232        FIXME:
 233        Can we always assume that sizeof(double) / sizeof(int) is 2?
 234        */
 235
 236     if (i_real == 4) {          /* Double precision reals */
 237         numb = 18 + 10 * mult;
 238     } else {                    /* Single precision reals */
 239         numb = (mult << 3) + 12;
 240     }
 241     ldkey = (iwkmax - iwkpt) / numb - 1;
 242     ldstp = mult * ldkey;
 243     ikh = ldkey << 1;   i4  = iwork(iwkmax, &iwkpt, ikh, i_int);
 244     ikh = ldkey << 1;   i5  = iwork(iwkmax, &iwkpt, ikh, i_int);
 245     ikh = ldstp << 1;   i6  = iwork(iwkmax, &iwkpt, ikh, i_real);
 246     ikh = ldstp * 6;    i7  = iwork(iwkmax, &iwkpt, ikh, i_int);
 247     ikh = ldkey << 1;   i8  = iwork(iwkmax, &iwkpt, ikh, i_real);
 248     ikh = ldkey << 1;   i9  = iwork(iwkmax, &iwkpt, ikh, i_real);
 249     ikh = ldkey << 1;   i9a = iwork(iwkmax, &iwkpt, ikh, i_real);
 250     ikh = ldkey << 1;   i10 = iwork(iwkmax, &iwkpt, ikh, i_int);
 251
 252     /* To convert to double precision, change RWRK to DWRK in the next CALL.
 253      */
 254     f2xact(nrow,
 255            ncol,
 256            &table[*ldtabl + 1],
 257            ldtabl,
 258            expect,
 259            percnt,
 260            emin,
 261            prt,
 262            pre,
 263            dwrk + i1,
 264            iwrk + i2,
 265            iwrk + i3,
 266            iwrk + i3a,
 267            iwrk + i3b,
 268            iwrk + i3c,
 269            iwrk + i4,
 270            &ldkey,
 271            iwrk + i5,
 272            dwrk + i6,
 273            &ldstp,
 274            iwrk + i7,
 275            dwrk + i8,
 276            dwrk + i9,
 277            dwrk + i9a,
 278            iwrk + i10,
 279            iwrk + iiwk,
 280            dwrk + irwk);
 281
 282 L_End:
 283     /* Free(equiv); */
 284     free(equiv);
 285     return;
 286 }
 287
 288 #undef rwrk
 289 #undef iwrk
 290 #undef dwrk
 291
 292
 293 /*
 294   -----------------------------------------------------------------------
 295   Name:         F2XACT
 296   Purpose:      Computes Fisher's exact test for a contingency table,
 297                 routine with workspace variables specified.
 298   Usage:        F2XACT (NROW, NCOL, TABLE, LDTABL, EXPECT, PERCNT,
 299                         EMIN, PRT, PRE, FACT, ICO, IRO, KYY, IDIF,
 300                         IRN, KEY, LDKEY, IPOIN, STP, LDSTP, IFRQ,
 301                         DLP, DSP, TM, KEY2, IWK, RWK)
 302   -----------------------------------------------------------------------
 303   */
 304 void
 305 f2xact(int *nrow, int *ncol, double *table, int *ldtabl,
 306        double *expect, double *percnt, double *emin, double *prt,
 307        double *pre, double *fact, int *ico, int *iro, int *kyy,
 308        int *idif, int *irn, int *key, int *ldkey, int *ipoin,
 309        double *stp, int *ldstp, int *ifrq, double *dlp, double *dsp,
 310        double *tm, int *key2, int *iwk, double *rwk)
 311 {
 312     /* IMAX is the largest representable int on the machine. */
 313     const int imax = SINT_MAX;
 314 //      const int imax = 2147483647; //xx: I DON¥T like this, and
 315 // thanks to the hint from Jason Turner I don't do it anymore. (R.D-U).
 316
 317     /* AMISS is a missing value indicator which is returned when the
 318        probability is not defined. */
 319     const double amiss = -12345.;
 320
 321     /* TOL is chosen as the square root of the smallest relative spacing. */
 322 #ifndef Macintosh
 323     const  static double tol = 3.45254e-7;
 324 #else
 325     static double tol = 3.45254e-7;
 326 #endif
 327     /* EMX is a large positive value used in comparing expected values. */
 328     const static double emx = 1e30;
 329
 330     /* Local variables {{any really need to be static ???}} */
 331     static int kval, kmax, jkey, last, ipsh, itmp, itop, jstp, ntot,
 332         jstp2, jstp3, jstp4, i, ii, j, k, n, iflag, ncell, ifreq, chisq,
 333         ikkey, ikstp, ikstp2, k1, kb, kd, ks,
 334         i31, i32, i33, i34, i35, i36, i37, i38, i39,
 335         i41, i42, i43, i44, i45, i46, i47, i48, i310, i311,
 336         nco, nrb, ipn, ipo, itp, nro, nro2;
 337     static double dspt, dd, df,ddf, drn,dro, emn, obs, obs2, obs3,
 338         pastp, pv, tmp;
 339     double d1;
 340 #ifndef USING_R
 341     double d2;
 342     static int ifault;
 343 #endif
 344     int nr_gt_nc=0;
 345
 346     /* Parameter adjustments */
 347     table -= *ldtabl + 1;
 348     --ico;
 349     --iro;
 350     --kyy;
 351     --idif;
 352     --irn;
 353     --key;
 354     --ipoin;
 355     --stp;
 356     --ifrq;
 357     --dlp;
 358     --dsp;
 359     --tm;
 360     --key2;
 361     --iwk;
 362     --rwk;
 363
 364
 365     /* Check table dimensions */
 366     if (*nrow > *ldtabl)
 367         prterr(1, "NROW must be less than or equal to LDTABL.");
 368     if (*ncol <= 1)
 369         prterr(4, "NCOL must be at least 2");
 370
 371     /* Initialize KEY array */
 372     for (i = 1; i <= *ldkey << 1; ++i) {
 373         key[i] = -9999;
 374         key2[i] = -9999;
 375     }
 376     /* Initialize parameters */
 377     *pre = 0.;
 378     itop = 0;
 379     if (*expect > 0.)
 380         emn = *emin;
 381     else
 382         emn = emx;
 383  if (*nrow > *ncol){
 384     nr_gt_nc =  1;
 385 }
 386 else{
 387          nr_gt_nc =  0;
 388 }
 389     /* nco := max(nrow, ncol) : */
 390     if(nr_gt_nc)
 391         nco = *nrow;
 392     else
 393         nco = *ncol;
 394     /* Initialize pointers for workspace */
 395     /* f3xact */
 396     i31 = 1;
 397     i32 = i31 + nco;
 398     i33 = i32 + nco;
 399     i34 = i33 + nco;
 400     i35 = i34 + nco;
 401     i36 = i35 + nco;
 402     i37 = i36 + nco;
 403     i38 = i37 + nco;
 404     i39 = i38 + 400;
 405     i310 = 1;
 406     i311 = 401;
 407     /* f4xact */
 408     k = *nrow + *ncol + 1;
 409     i41 = 1;
 410     i42 = i41 + k;
 411     i43 = i42 + k;
 412     i44 = i43 + k;
 413     i45 = i44 + k;
 414     i46 = i45 + k;
 415     i47 = i46 + k * nco;
 416     i48 = 1;
 417
 418     /* Compute row marginals and total */
 419     ntot = 0;
 420     for (i = 1; i <= *nrow; ++i) {
 421         iro[i] = 0;
 422         for (j = 1; j <= *ncol; ++j) {
 423             if (table[i + j * *ldtabl] < -1e-4)
 424                 prterr(2, "All elements of TABLE must be positive.");
 425             iro[i] += (int) table[i + j * *ldtabl];
 426         }
 427         ntot += iro[i];
 428     }
 429
 430     if (ntot == 0) {
 431         prterr(3, "All elements of TABLE are zero.\n"
 432                "PRT and PRE are set to missing values.");
 433         *pre = *prt = amiss;
 434         return;
 435     }
 436
 437     /* Column marginals */
 438     for (i = 1; i <= *ncol; ++i) {
 439         ico[i] = 0;
 440         for (j = 1; j <= *nrow; ++j)
 441             ico[i] += (int) table[j + i * *ldtabl];
 442     }
 443
 444     /* sort marginals */
 445     isort(nrow, &iro[1]);
 446     isort(ncol, &ico[1]);
 447
 448     /*  Determine row and column marginals.
 449         Define max(nrow,ncol) =: nco >= nro := min(nrow,ncol)
 450         nco is defined above
 451
 452         Swap marginals if necessary to  ico[1:nco] & iro[1:nro]
 453      */
 454     if (nr_gt_nc) {
 455         nro = *ncol;
 456         /* Swap marginals */
 457         for (i = 1; i <= nco; ++i) {
 458             itmp = iro[i];
 459             if (i <= nro)
 460                 iro[i] = ico[i];
 461             ico[i] = itmp;
 462         }
 463     } else
 464         nro = *nrow;
 465
 466
 467     /* Get multiplers for stack */
 468     kyy[1] = 1;
 469     for (i = 2; i <= nro; ++i) {
 470         /* Hash table multipliers */
 471         if (iro[i - 1] + 1 <= imax / kyy[i - 1]) {
 472             kyy[i] = kyy[i - 1] * (iro[i - 1] + 1);
 473             j /= kyy[i - 1];
 474         } else
 475             goto L_ERR_5;
 476     }
 477     /* Maximum product */
 478     if (iro[nro - 1] + 1 <= imax / kyy[nro - 1]) {
 479         kmax = (iro[nro] + 1) * kyy[nro - 1];
 480     } else {
 481     L_ERR_5:
 482         prterr(5, "The hash table key cannot be computed because "
 483                "the largest key\n"
 484                "is larger than the largest representable int.\n"
 485                "The algorithm cannot proceed.\n"
 486                "Reduce the workspace size, or use `exact = FALSE'.");
 487         return;
 488     }
 489
 490     /* Compute log factorials */
 491     fact[0] = 0.;
 492     fact[1] = 0.;
 493     if(ntot >= 2) fact[2] = log(2.);
 494 /* MM: old code assuming log() to be SLOW */
 495     for (i = 3; i <= ntot; i += 2) {
 496         fact[i] = fact[i - 1] + log((double) i);
 497         j = i + 1;
 498         if (j <= ntot)
 499             fact[j] = fact[i] + fact[2] + fact[j / 2] - fact[j / 2 - 1];
 500     }
 501     /* Compute obs := observed path length */
 502     obs = tol;
 503     ntot = 0;
 504     for (j = 1; j <= nco; ++j) {
 505         dd = 0.;
 506         for (i = 1; i <= nro; ++i) {
 507             if (nr_gt_nc) {
 508                 dd += fact[(int) table[j + i * *ldtabl]];
 509                 ntot +=    (int) table[j + i * *ldtabl];
 510             } else {
 511                 dd += fact[(int) table[i + j * *ldtabl]];
 512                 ntot +=    (int) table[i + j * *ldtabl];
 513             }
 514         }
 515         obs += fact[ico[j]] - dd;
 516     }
 517     /* Denominator of observed table: DRO */
 518     dro = f9xact(&nro, &ntot, &iro[1], fact);
 519     *prt = exp(obs - dro);
 520     /* Initialize pointers */
 521     k = nco;
 522     last = *ldkey + 1;
 523     jkey = *ldkey + 1;
 524     jstp = *ldstp + 1;
 525     jstp2 = *ldstp * 3 + 1;
 526     jstp3 = (*ldstp << 2) + 1;
 527     jstp4 = *ldstp * 5 + 1;
 528     ikkey = 0;
 529     ikstp = 0;
 530     ikstp2 = *ldstp << 1;
 531     ipo = 1;
 532     ipoin[1] = 1;
 533     stp[1] = 0.;
 534     ifrq[1] = 1;
 535     ifrq[ikstp2 + 1] = -1;
 536
 537 Outer_Loop:
 538     kb = nco - k + 1;
 539     ks = 0;
 540     n = ico[kb];
 541     kd = nro + 1;
 542     kmax = nro;
 543     /* IDIF is the difference in going to the daughter */
 544     for (i = 1; i <= nro; ++i)
 545         idif[i] = 0;
 546
 547     /* Generate the first daughter */
 548     do {
 549         --kd;
 550         ntot = min(n, iro[kd]);
 551         idif[kd] = ntot;
 552         if (idif[kmax] == 0)
 553             --kmax;
 554         n -= ntot;
 555     }
 556     while (n > 0 && kd != 1);
 557
 558     if (n != 0) {
 559         goto L310;
 560     }
 561
 562     k1 = k - 1;
 563     n = ico[kb];
 564     ntot = 0;
 565     for (i = kb + 1; i <= nco; ++i)
 566         ntot += ico[i];
 567
 568
 569 L150:
 570     /* Arc to daughter length=ICO(KB) */
 571     for (i = 1; i <= nro; ++i)
 572         irn[i] = iro[i] - idif[i];
 573
 574     /* Sort irn */
 575     if (k1 > 1) {
 576         if (nro == 2) {
 577             if (irn[1] > irn[2]) {
 578                 ii = irn[1];
 579                 irn[1] = irn[2];
 580                 irn[2] = ii;
 581             }
 582         } else if (nro == 3) {
 583             ii = irn[1];
 584             if (ii > irn[3]) {
 585                 if (ii > irn[2]) {
 586                     if (irn[2] > irn[3]) {
 587                         irn[1] = irn[3];
 588                         irn[3] = ii;
 589                     } else {
 590                         irn[1] = irn[2];
 591                         irn[2] = irn[3];
 592                         irn[3] = ii;
 593                     }
 594                 } else {
 595                     irn[1] = irn[3];
 596                     irn[3] = irn[2];
 597                     irn[2] = ii;
 598                 }
 599             } else if (ii > irn[2]) {
 600                 irn[1] = irn[2];
 601                 irn[2] = ii;
 602             } else if (irn[2] > irn[3]) {
 603                 ii = irn[2];
 604                 irn[2] = irn[3];
 605                 irn[3] = ii;
 606             }
 607         } else {
 608             for (j = 2; j <= nro; ++j) {
 609                 i = j - 1;
 610                 ii = irn[j];
 611
 612                 while (ii < irn[i]) {
 613                     irn[i + 1] = irn[i];
 614                     --i;
 615                     if (i == 0)
 616                         break;
 617                 }
 618                 irn[i + 1] = ii;
 619             }
 620         }
 621         /* Adjust start for zero */
 622         for (i = 1; i <= nro; ++i) {
 623             if (irn[i] != 0)
 624                 break;
 625         }
 626
 627         nrb = i;
 628         nro2 = nro - i + 1;
 629     } else {
 630         nrb = 1;
 631         nro2 = nro;
 632     }
 633     /* Some table values */
 634     ddf = f9xact(&nro, &n, &idif[1], fact);
 635     drn = f9xact(&nro2, &ntot, &irn[nrb], fact) - dro + ddf;
 636     /* Get hash value */
 637     if (k1 > 1) {
 638         kval = irn[1] + irn[2] * kyy[2];
 639         for (i = 3; i <= nro; ++i) {
 640             kval += irn[i] * kyy[i];
 641         }
 642         /* Get hash table entry */
 643         i = kval % (*ldkey << 1) + 1;
 644         /* Search for unused location */
 645         for (itp = i; itp <= *ldkey << 1; ++itp) {
 646             ii = key2[itp];
 647             if (ii == kval) {
 648                 goto L240;
 649             } else if (ii < 0) {
 650                 key2[itp] = kval;
 651                 dlp[itp] = 1.;
 652                 dsp[itp] = 1.;
 653                 goto L240;
 654             }
 655         }
 656
 657         for (itp = 1; itp <= i - 1; ++itp) {
 658             ii = key2[itp];
 659             if (ii == kval) {
 660                 goto L240;
 661             } else if (ii < 0) {
 662                 key2[itp] = kval;
 663                 dlp[itp] = 1.;
 664                 goto L240;
 665             }
 666         }
 667
 668         /* KH
 669            prterr(6, "LDKEY is too small.\n"
 670            "It is not possible to give the value of LDKEY required,\n"
 671            "but you could try doubling LDKEY (and possibly LDSTP).");
 672            */
 673         prterr(6, "LDKEY is too small for this problem.\n"
 674                "Try increasing the size of the workspace.");
 675     }
 676
 677 L240:
 678     ipsh = (1);
 679     /* Recover pastp */
 680     ipn = ipoin[ipo + ikkey];
 681     pastp = stp[ipn + ikstp];
 682     ifreq = ifrq[ipn + ikstp];
 683     /* Compute shortest and longest path */
 684     if (k1 > 1) {
 685         obs2 = obs - fact[ico[kb + 1]] - fact[ico[kb + 2]] - ddf;
 686         for (i = 3; i <= k1; ++i) {
 687             obs2 -= fact[ico[kb + i]];
 688         }
 689         if (dlp[itp] > 0.) {
 690             dspt = obs - obs2 - ddf;
 691             /* Compute longest path */
 692             dlp[itp] = 0.;
 693             f3xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dlp[itp],
 694                    &ntot, fact, &iwk[i31], &iwk[i32], &iwk[i33],
 695                    &iwk[i34], &iwk[i35], &iwk[i36], &iwk[i37],
 696                    &iwk[i38], &iwk[i39], &rwk[i310], &rwk[i311], &tol);
 697             dlp[itp] = min(0., dlp[itp]);
 698             /* Compute shortest path */
 699             dsp[itp] = dspt;
 700             f4xact(&nro2, &irn[nrb], &k1, &ico[kb + 1], &dsp[itp], fact,
 701                    &iwk[i47], &iwk[i41], &iwk[i42], &iwk[i43],
 702                    &iwk[i44], &iwk[i45], &iwk[i46], &rwk[i48], &tol);
 703             dsp[itp] = min(0., dsp[itp] - dspt);
 704             /* Use chi-squared approximation? */
 705             if ((irn[nrb] * ico[kb + 1]) > ntot * emn) {
 706                 ncell = 0.;
 707                 for (i = 0; i < nro2; ++i)
 708                     for (j = 1; j <= k1; ++j)
 709                         if (irn[nrb + i] * ico[kb + j] >= ntot * *expect)
 710                             ncell++;
 711
 712                 if (ncell * 100 >= k1 * nro2 * *percnt) {
 713                     tmp = 0.;
 714                     for (i = 0; i < nro2; ++i)
 715                         tmp += (fact[irn[nrb + i]] -
 716                                 fact[irn[nrb + i] - 1]);
 717                     tmp *= k1 - 1;
 718                     for (j = 1; j <= k1; ++j)
 719                         tmp += (nro2 - 1) * (fact[ico[kb + j]] -
 720                                              fact[ico[kb + j] - 1]);
 721                     df = (double) ((nro2 - 1) * (k1 - 1));
 722                     tmp += df * 1.83787706640934548356065947281;
 723                     tmp -= (nro2 * k1 - 1) * (fact[ntot] - fact[ntot - 1]);
 724                     tm[itp] = (obs - dro) * -2. - tmp;
 725                 } else {
 726                     /* tm(itp) set to a flag value */
 727                     tm[itp] = -9876.;
 728                 }
 729             } else {
 730                 tm[itp] = -9876.;
 731             }
 732         }
 733         obs3 = obs2 - dlp[itp];
 734         obs2 -= dsp[itp];
 735         if (tm[itp] == -9876.) {
 736             chisq = (0);
 737         } else {
 738             chisq = (1);
 739             tmp = tm[itp];
 740         }
 741     } else {
 742         obs2 = obs - drn - dro;
 743         obs3 = obs2;
 744     }
 745
 746 L300:
 747     /* Process node with new PASTP */
 748     if (pastp <= obs3) {
 749         /* Update pre */
 750         *pre += (double) ifreq * exp(pastp + drn);
 751     } else if (pastp < obs2) {
 752         if (chisq) {
 753             df = (double) ((nro2 - 1) * (k1 - 1));
 754 #ifdef USING_R
 755             pv = pgamma(fmax2(0., tmp + (pastp + drn) * 2.) / 2.,
 756                         df / 2., /*scale = */ 1.,
 757                         /*lower_tail = */FALSE, /*log_p = */ FALSE);
 758 #else
 759             d1 = max(0., tmp + (pastp + drn) * 2.) / 2.;
 760             d2 = df / 2.;
 761             pv = 1. - gammds(&d1, &d2, &ifault);
 762 #endif
 763             *pre += (double) ifreq * exp(pastp + drn) * pv;
 764         } else {
 765             /* Put daughter on queue */
 766             d1 = pastp + ddf;
 767             f5xact(&d1, &tol, &kval, &key[jkey], ldkey, &ipoin[jkey],
 768                    &stp[jstp], ldstp, &ifrq[jstp], &ifrq[jstp2],
 769                    &ifrq[jstp3], &ifrq[jstp4], &ifreq, &itop, &ipsh);
 770             ipsh = (0);
 771         }
 772     }
 773     /* Get next PASTP on chain */
 774     ipn = ifrq[ipn + ikstp2];
 775     if (ipn > 0) {
 776         pastp = stp[ipn + ikstp];
 777         ifreq = ifrq[ipn + ikstp];
 778         goto L300;
 779     }
 780     /* Generate a new daughter node */
 781     f7xact(&kmax, &iro[1], &idif[1], &kd, &ks, &iflag);
 782     if (iflag != 1) {
 783         goto L150;
 784     }
 785
 786 L310:
 787     /* Go get a new mother from stage K */
 788     do {
 789         iflag = 1;
 790         f6xact(&nro, &iro[1], &iflag, &kyy[1], &key[ikkey + 1], ldkey,
 791                &last, &ipo);
 792         /* Update pointers */
 793         if (iflag != 3)
 794             goto Outer_Loop;
 795         /* else  iflag == 3 : no additional nodes to process */
 796         --k;
 797         itop = 0;
 798         ikkey = jkey - 1;
 799         ikstp = jstp - 1;
 800         ikstp2 = jstp2 - 1;
 801         jkey = *ldkey - jkey + 2;
 802         jstp = *ldstp - jstp + 2;
 803         jstp2 = (*ldstp << 1) + jstp;
 804         for (i = 1; i <= *ldkey << 1; ++i)
 805             key2[i] = -9999;
 806
 807     } while (k >= 2);
 808 }
 809
 810 /*
 811   -----------------------------------------------------------------------
 812   Name:       F3XACT
 813   Purpose:    Computes the shortest path length for a given table.
 814   Usage:      F3XACT (NROW, IROW, NCOL, ICOL, DLP, MM, FACT, ICO, IRO,
 815                       IT, LB, NR, NT, NU, ITC, IST, STV, ALEN, TOL)
 816   Arguments:
 817     NROW    - The number of rows in the table.                  (Input)
 818     IROW    - Vector of length NROW containing the row sums
 819               for the table.                                    (Input)
 820     NCOL    - The number of columns in the table.               (Input)
 821     ICOL    - Vector of length K containing the column sums
 822               for the table.                                    (Input)
 823     DLP     - The longest path for the table.                   (Output)
 824     MM      - The total count in the table.                     (Output)
 825     FACT    - Vector containing the logarithms of factorials.   (Input)
 826     ICO     - Work vector of length MAX(NROW,NCOL).
 827     IRO     - Work vector of length MAX(NROW,NCOL).
 828     IT      - Work vector of length MAX(NROW,NCOL).
 829     LB      - Work vector of length MAX(NROW,NCOL).
 830     NR      - Work vector of length MAX(NROW,NCOL).
 831     NT      - Work vector of length MAX(NROW,NCOL).
 832     NU      - Work vector of length MAX(NROW,NCOL).
 833     ITC     - Work vector of length 400.
 834     IST     - Work vector of length 400.
 835     STV     - Work vector of length 400.
 836     ALEN    - Work vector of length MAX(NROW,NCOL).
 837     TOL     - Tolerance.                                        (Input)
 838   -----------------------------------------------------------------------
 839   */
 840
 841 void
 842 f3xact(int *nrow, int *irow, int *ncol, int *icol, double *dlp,
 843        int *mm, double *fact, int *ico, int *iro, int *it,
 844        int *lb, int *nr, int *nt, int *nu, int *itc, int *ist,
 845        double *stv, double *alen, const double *tol)
 846 {
 847     /* Initialized data */
 848     static int ldst = 200;
 849     static int nst = 0;
 850     static int nitc = 0;
 851
 852     /* Local variables */
 853     static int xmin;
 854     static int i, k;
 855     static double v;
 856     static int n11, n12, ii, nn, ks, ic1, ic2, nc1, nn1;
 857     static int nr1, nco;
 858     static double val;
 859     static int nct, ipn, irl, key, lev, itp, nro;
 860     static double vmn;
 861     static int nrt, kyy, nc1s;
 862
 863     /* Parameter adjustments */
 864     --stv;
 865     --ist;
 866     --itc;
 867     --nu;
 868     --nt;
 869     --nr;
 870     --lb;
 871     --it;
 872     --iro;
 873     --ico;
 874     --icol;
 875     --irow;
 876
 877     /* Function Body */
 878     for (i = 0; i <= *ncol; ++i) {
 879         alen[i] = 0.;
 880     }
 881     for (i = 1; i <= 400; ++i) {
 882         ist[i] = -1;
 883     }
 884     /* nrow is 1 */
 885     if (*nrow <= 1) {
 886         if (*nrow > 0) {
 887             *dlp -= fact[icol[1]];
 888             for (i = 2; i <= *ncol; ++i) {
 889                 *dlp -= fact[icol[i]];
 890             }
 891         }
 892         return;
 893     }
 894     /* ncol is 1 */
 895     if (*ncol <= 1) {
 896         if (*ncol > 0) {
 897             *dlp = *dlp - fact[irow[1]] - fact[irow[2]];
 898             for (i = 3; i <= *nrow; ++i) {
 899                 *dlp -= fact[irow[i]];
 900             }
 901         }
 902         return;
 903     }
 904     /* 2 by 2 table */
 905     if (*nrow * *ncol == 4) {
 906         n11 = (irow[1] + 1) * (icol[1] + 1) / (*mm + 2);
 907         n12 = irow[1] - n11;
 908         *dlp = *dlp - fact[n11] - fact[n12] - fact[icol[1] - n11]
 909             - fact[icol[2] - n12];
 910         return;
 911     }
 912     /* Test for optimal table */
 913     val = 0.;
 914     xmin = (0);
 915     if (irow[*nrow] <= irow[1] + *ncol) {
 916         f10act(nrow, &irow[1], ncol, &icol[1], &val, &xmin, fact,
 917                &lb[1], &nu[1], &nr[1]);
 918     }
 919     if (! xmin) {
 920         if (icol[*ncol] <= icol[1] + *nrow) {
 921             f10act(ncol, &icol[1], nrow, &irow[1], &val, &xmin, fact,
 922                    &lb[1], &nu[1], &nr[1]);
 923         }
 924     }
 925
 926     if (xmin) {
 927         *dlp -= val;
 928         return;
 929     }
 930     /* Setup for dynamic programming */
 931     nn = *mm;
 932     /* Minimize ncol */
 933     if (*nrow >= *ncol) {
 934         nro = *nrow;
 935         nco = *ncol;
 936         for (i = 1; i <= *nrow; ++i) {
 937             iro[i] = irow[i];
 938         }
 939         ico[1] = icol[1];
 940         nt[1] = nn - ico[1];
 941         for (i = 2; i <= *ncol; ++i) {
 942             ico[i] = icol[i];
 943             nt[i] = nt[i - 1] - ico[i];
 944         }
 945     } else {
 946         nro = *ncol;
 947         nco = *nrow;
 948         ico[1] = irow[1];
 949         nt[1] = nn - ico[1];
 950         for (i = 2; i <= *nrow; ++i) {
 951             ico[i] = irow[i];
 952             nt[i] = nt[i - 1] - ico[i];
 953         }
 954         for (i = 1; i <= *ncol; ++i)
 955             iro[i] = icol[i];
 956     }
 957     /* Initialize pointers */
 958     vmn = 1e10;
 959     nc1s = nco - 1;
 960     irl = 1;
 961     ks = 0;
 962     k = ldst;
 963     kyy = ico[nco] + 1;
 964
 965 LnewNode: /* Setup to generate new node */
 966
 967     lev = 1;
 968     nr1 = nro - 1;
 969     nrt = iro[irl];
 970     nct = ico[1];
 971     lb[1] = (int) ((double) ((nrt + 1) * (nct + 1)) /
 972                     (double) (nn + nr1 * nc1s + 1) - *tol) - 1;
 973     nu[1] = (int) ((double) ((nrt + nc1s) * (nct + nr1)) /
 974                     (double) (nn + nr1 + nc1s)) - lb[1] + 1;
 975     nr[1] = nrt - lb[1];
 976
 977 LoopNode: /* Generate a node */
 978     --nu[lev];
 979     if (nu[lev] == 0) {
 980         if (lev == 1)
 981             goto L200;
 982
 983         --lev;
 984         goto LoopNode;
 985     }
 986     ++lb[lev];
 987     --nr[lev];
 988 L120:
 989     alen[lev] = alen[lev - 1] + fact[lb[lev]];
 990     if (lev < nc1s) {
 991         nn1 = nt[lev];
 992         nrt = nr[lev];
 993         ++lev;
 994         nc1 = nco - lev;
 995         nct = ico[lev];
 996         lb[lev] = (int) ((double) ((nrt + 1) * (nct + 1)) /
 997                           (double) (nn1 + nr1 * nc1 + 1) - *tol);
 998         nu[lev] = (int) ((double) ((nrt + nc1) * (nct + nr1)) /
 999                           (double) (nn1 + nr1 + nc1) - lb[lev] + 1);
1000         nr[lev] = nrt - lb[lev];
1001         goto L120;
1002     }
1003     alen[nco] = alen[lev] + fact[nr[lev]];
1004     lb[nco] = nr[lev];
1005
1006     v = val + alen[nco];
1007     if (nro == 2) {
1008         /* Only 1 row left */
1009         v = v + fact[ico[1] - lb[1]] + fact[ico[2] - lb[2]];
1010         for (i = 3; i <= nco; ++i) {
1011             v += fact[ico[i] - lb[i]];
1012         }
1013         if (v < vmn) {
1014             vmn = v;
1015         }
1016     } else if (nro == 3 && nco == 2) {
1017         /* 3 rows and 2 columns */
1018         nn1 = nn - iro[irl] + 2;
1019         ic1 = ico[1] - lb[1];
1020         ic2 = ico[2] - lb[2];
1021         n11 = (iro[irl + 1] + 1) * (ic1 + 1) / nn1;
1022         n12 = iro[irl + 1] - n11;
1023         v = v + fact[n11] + fact[n12] + fact[ic1 - n11]
1024             + fact[ic2 - n12];
1025         if (v < vmn) {
1026             vmn = v;
1027         }
1028     } else {
1029         /* Column marginals are new node */
1030         for (i = 1; i <= nco; ++i) {
1031             it[i] = ico[i] - lb[i];
1032         }
1033         /* Sort column marginals */
1034         if (nco == 2) {
1035             if (it[1] > it[2]) {
1036                 ii = it[1];
1037                 it[1] = it[2];
1038                 it[2] = ii;
1039             }
1040         } else if (nco == 3) {
1041             ii = it[1];
1042             if (ii > it[3]) {
1043                 if (ii > it[2]) {
1044                     if (it[2] > it[3]) {
1045                         it[1] = it[3];
1046                         it[3] = ii;
1047                     } else {
1048                         it[1] = it[2];
1049                         it[2] = it[3];
1050                         it[3] = ii;
1051                     }
1052                 } else {
1053                     it[1] = it[3];
1054                     it[3] = it[2];
1055                     it[2] = ii;
1056                 }
1057             } else if (ii > it[2]) {
1058                 it[1] = it[2];
1059                 it[2] = ii;
1060             } else if (it[2] > it[3]) {
1061                 ii = it[2];
1062                 it[2] = it[3];
1063                 it[3] = ii;
1064             }
1065         } else {
1066             isort(&nco, &it[1]);
1067         }
1068         /* Compute hash value */
1069         key = it[1] * kyy + it[2];
1070         for (i = 3; i <= nco; ++i) {
1071             key = it[i] + key * kyy;
1072         }
1073         if(key < 0)
1074                 //PROBLEM "Bug in FEXACT: gave negative key" RECOVER(NULL_ENTRY);
1075                 printf("Bug in FEXACT: gave negative key \n"); //xx:another one of my ugly kluges (R.D-U)
1076
1077         /* Table index */
1078         ipn = key % ldst + 1;
1079
1080         /* Find empty position */
1081         for (itp = ipn, ii = ks + ipn; itp <= ldst; ++itp, ++ii) {
1082             if (ist[ii] < 0) {
1083                 goto L180;
1084             } else if (ist[ii] == key) {
1085                 goto L190;
1086             }
1087         }
1088
1089         for (itp = 1, ii = ks + 1; itp <= ipn - 1; ++itp, ++ii) {
1090             if (ist[ii] < 0) {
1091                 goto L180;
1092             } else if (ist[ii] == key) {
1093                 goto L190;
1094             }
1095         }
1096
1097         prterr(30, "Stack length exceeded in f3xact.\n"
1098                "This problem should not occur.");
1099
1100 L180: /* Push onto stack */
1101         ist[ii] = key;
1102         stv[ii] = v;
1103         ++nst;
1104         ii = nst + ks;
1105         itc[ii] = itp;
1106         goto LoopNode;
1107
1108 L190: /* Marginals already on stack */
1109         stv[ii] = min(v, stv[ii]);
1110     }
1111     goto LoopNode;
1112
1113
1114 L200: /* Pop item from stack */
1115     if (nitc > 0) {
1116         /* Stack index */
1117         itp = itc[nitc + k] + k;
1118         --nitc;
1119         val = stv[itp];
1120         key = ist[itp];
1121         ist[itp] = -1;
1122         /* Compute marginals */
1123         for (i = nco; i >= 2; --i) {
1124             ico[i] = key % kyy;
1125             key /= kyy;
1126         }
1127         ico[1] = key;
1128         /* Set up nt array */
1129         nt[1] = nn - ico[1];
1130         for (i = 2; i <= nco; ++i)
1131             nt[i] = nt[i - 1] - ico[i];
1132
1133         /* Test for optimality (L90) */
1134         xmin = (0);
1135         if (iro[nro] <= iro[irl] + nco) {
1136             f10act(&nro, &iro[irl], &nco, &ico[1], &val, &xmin, fact,
1137                    &lb[1], &nu[1], &nr[1]);
1138         }
1139         if (!xmin && ico[nco] <= ico[1] + nro)
1140             f10act(&nco, &ico[1], &nro, &iro[irl], &val, &xmin, fact,
1141                    &lb[1], &nu[1], &nr[1]);
1142         if (xmin) {
1143             if (vmn > val)
1144                 vmn = val;
1145             goto L200;
1146         }
1147         else goto LnewNode;
1148
1149     } else if (nro > 2 && nst > 0) {
1150         /* Go to next level */
1151         nitc = nst;
1152         nst = 0;
1153         k = ks;
1154         ks = ldst - ks;
1155         nn -= iro[irl];
1156         ++irl;
1157         --nro;
1158         goto L200;
1159     }
1160
1161     *dlp -= vmn;
1162 }
1163
1164 /*
1165   -----------------------------------------------------------------------
1166   Name:       F4XACT
1167   Purpose:    Computes the longest path length for a given table.
1168   Usage:      CALL F4XACT (NROW, IROW, NCOL, ICOL, DSP, FACT, ICSTK,
1169                           NCSTK, LSTK, MSTK, NSTK, NRSTK, IRSTK, YSTK,
1170                           TOL)
1171   Arguments:
1172      NROW   - The number of rows in the table.  (Input)
1173      IROW   - Vector of length NROW containing the row sums for the
1174               table.  (Input)
1175      NCOL   - The number of columns in the table.  (Input)
1176      ICOL   - Vector of length K containing the column sums for the
1177               table.  (Input)
1178      DSP    - The shortest path for the table.  (Output)
1179      FACT   - Vector containing the logarithms of factorials.  (Input)
1180      ICSTK  - NCOL by NROW+NCOL+1 work array.
1181      NCSTK  - Work vector of length NROW+NCOL+1.
1182      LSTK   - Work vector of length NROW+NCOL+1.
1183      MSTK   - Work vector of length NROW+NCOL+1.
1184      NSTK   - Work vector of length NROW+NCOL+1.
1185      NRSTK  - Work vector of length NROW+NCOL+1.
1186      IRSTK  - NROW by MAX(NROW,NCOL) work array.
1187      YSTK   - Work vector of length NROW+NCOL+1.
1188      TOL    - Tolerance.  (Input)
1189   -----------------------------------------------------------------------
1190   */
1191
1192 void
1193 f4xact(int *nrow, int *irow, int *ncol, int *icol, double *dsp,
1194        double *fact, int *icstk, int *ncstk, int *lstk, int *mstk,
1195        int *nstk, int *nrstk, int *irstk, double *ystk, const double *tol)
1196 {
1197     /* System generated locals */
1198     int ikh;
1199
1200     /* Local variables */
1201     int i, j, k, l, m, n, mn, ic1, ir1, ict, irt, istk, nco, nro;
1202     double y, amx;
1203
1204     /* Parameter adjustments */
1205     irstk -= *nrow + 1;
1206     --irow;
1207     icstk -= *ncol + 1;
1208     --icol;
1209     --ncstk;
1210     --lstk;
1211     --mstk;
1212     --nstk;
1213     --nrstk;
1214     --ystk;
1215
1216     /* Function Body */
1217     /* Take care of the easy cases first */
1218     if (*nrow == 1) {
1219         for (i = 1; i <= *ncol; ++i) {
1220             *dsp -= fact[icol[i]];
1221         }
1222         return;
1223     }
1224     if (*ncol == 1) {
1225         for (i = 1; i <= *nrow; ++i) {
1226             *dsp -= fact[irow[i]];
1227         }
1228         return;
1229     }
1230     if (*nrow * *ncol == 4) {
1231         if (irow[2] <= icol[2]) {
1232             *dsp = *dsp - fact[irow[2]] - fact[icol[1]]
1233                 - fact[icol[2] - irow[2]];
1234         } else {
1235             *dsp = *dsp - fact[icol[2]] - fact[irow[1]]
1236                 - fact[irow[2] - icol[2]];
1237         }
1238         return;
1239     }
1240     /* initialization before loop */
1241     for (i = 1; i <= *nrow; ++i) {
1242         irstk[i + *nrow] = irow[*nrow - i + 1];
1243     }
1244     for (j = 1; j <= *ncol; ++j) {
1245         icstk[j + *ncol] = icol[*ncol - j + 1];
1246     }
1247
1248     nro = *nrow;
1249     nco = *ncol;
1250     nrstk[1] = nro;
1251     ncstk[1] = nco;
1252     ystk[1] = 0.;
1253     y = 0.;
1254     istk = 1;
1255     l = 1;
1256     amx = 0.;
1257
1258     /* First LOOP */
1259     do {
1260         ir1 = irstk[istk * *nrow + 1];
1261         ic1 = icstk[istk * *ncol + 1];
1262         if (ir1 > ic1) {
1263             if (nro >= nco) {
1264                 m = nco - 1;    n = 2;
1265             } else {
1266                 m = nro;        n = 1;
1267             }
1268         } else if (ir1 < ic1) {
1269             if (nro <= nco) {
1270                 m = nro - 1;    n = 1;
1271             } else {
1272                 m = nco;        n = 2;
1273             }
1274         } else {
1275             if (nro <= nco) {
1276                 m = nro - 1;    n = 1;
1277             } else {
1278                 m = nco - 1;    n = 2;
1279             }
1280         }
1281
1282     L60:
1283         if (n == 1) {
1284             i = l; j = 1;
1285         } else {
1286             i = 1; j = l;
1287         }
1288
1289         irt = irstk[i + istk * *nrow];
1290         ict = icstk[j + istk * *ncol];
1291         mn = irt;
1292         if (mn > ict) {
1293             mn = ict;
1294         }
1295         y += fact[mn];
1296         if (irt == ict) {
1297             --nro;
1298             --nco;
1299             f11act(&irstk[istk * *nrow + 1], &i, &nro,
1300                    &irstk[(istk + 1) * *nrow + 1]);
1301             f11act(&icstk[istk * *ncol + 1], &j, &nco,
1302                    &icstk[(istk + 1) * *ncol + 1]);
1303         } else if (irt > ict) {
1304             --nco;
1305             f11act(&icstk[istk * *ncol + 1], &j, &nco,
1306                    &icstk[(istk + 1) * *ncol + 1]);
1307             ikh = irt - ict;
1308             f8xact(&irstk[istk * *nrow + 1], &ikh, &i,
1309                    &nro, &irstk[(istk + 1) * *nrow + 1]);
1310         } else {
1311             --nro;
1312             f11act(&irstk[istk * *nrow + 1], &i, &nro,
1313                    &irstk[(istk + 1) * *nrow + 1]);
1314             ikh = ict - irt;
1315             f8xact(&icstk[istk * *ncol + 1], &ikh, &j,
1316                    &nco, &icstk[(istk + 1) * *ncol + 1]);
1317         }
1318
1319         if (nro == 1) {
1320             for (k = 1; k <= nco; ++k) {
1321                 y += fact[icstk[k + (istk + 1) * *ncol]];
1322             }
1323             break;
1324         }
1325         if (nco == 1) {
1326             for (k = 1; k <= nro; ++k) {
1327                 y += fact[irstk[k + (istk + 1) * *nrow]];
1328             }
1329             break;
1330         }
1331
1332         lstk[istk] = l;
1333         mstk[istk] = m;
1334         nstk[istk] = n;
1335         ++istk;
1336         nrstk[istk] = nro;
1337         ncstk[istk] = nco;
1338         ystk[istk] = y;
1339         l = 1;
1340     } while(1);/* end do */
1341
1342 /* L90:*/
1343     if (y > amx) {
1344         amx = y;
1345         if (*dsp - amx <= *tol) {
1346             *dsp = 0.;
1347             return;
1348         }
1349     }
1350
1351 L100:
1352     --istk;
1353     if (istk == 0) {
1354         *dsp -= amx;
1355         if (*dsp - amx <= *tol) {
1356             *dsp = 0.;
1357         }
1358         return;
1359     }
1360     l = lstk[istk] + 1;
1361
1362 /* L110: */
1363     for(;; ++l) {
1364         if (l > mstk[istk])     goto L100;
1365
1366         n = nstk[istk];
1367         nro = nrstk[istk];
1368         nco = ncstk[istk];
1369         y = ystk[istk];
1370         if (n == 1) {
1371             if (irstk[l     + istk * *nrow] <
1372                 irstk[l - 1 + istk * *nrow])    goto L60;
1373         }
1374         else if (n == 2) {
1375             if (icstk[l     + istk * *ncol] <
1376                 icstk[l - 1 + istk * *ncol])    goto L60;
1377         }
1378     }
1379 }
1380
1381 /*
1382   -----------------------------------------------------------------------
1383   Name:       F5XACT
1384   Purpose:    Put node on stack in network algorithm.
1385   Usage:      CALL F5XACT (PASTP, TOL, KVAL, KEY, LDKEY, IPOIN, STP,
1386                           LDSTP, IFRQ, NPOIN, NR, NL, IFREQ, ITOP,
1387                           IPSH)
1388   Arguments:
1389      PASTP  - The past path length.                             (Input)
1390      TOL    - Tolerance for equivalence of past path lengths.   (Input)
1391      KVAL   - Key value.                                        (Input)
1392      KEY    - Vector of length LDKEY containing the key values. (in/out)
1393      LDKEY  - Length of vector KEY.                             (Input)
1394      IPOIN  - Vector of length LDKEY pointing to the
1395               linked list of past path lengths.                 (in/out)
1396      STP    - Vector of length LSDTP containing the
1397               linked lists of past path lengths.                (in/out)
1398      LDSTP  - Length of vector STP.                             (Input)
1399      IFRQ   - Vector of length LDSTP containing the past path
1400               frequencies.                                      (in/out)
1401      NPOIN  - Vector of length LDSTP containing the pointers to
1402               the next past path length.                        (in/out)
1403      NR     - Vector of length LDSTP containing the right object
1404               pointers in the tree of past path lengths.        (in/out)
1405      NL     - Vector of length LDSTP containing the left object
1406               pointers in the tree of past path lengths.        (in/out)
1407      IFREQ  - Frequency of the current path length.             (Input)
1408      ITOP   - Pointer to the top of STP.                        (Input)
1409      IPSH   - Option parameter.                                 (Input)
1410               If IPSH is true, the past path length is found in the
1411               table KEY.  Otherwise the location of the past path
1412               length is assumed known and to have been found in
1413               a previous call. ==>>>>> USING "static" variables
1414   -----------------------------------------------------------------------
1415   */
1416
1417 void
1418 f5xact(double *pastp, const double *tol, int *kval, int *key, int *ldkey,
1419        int *ipoin, double *stp, int *ldstp, int *ifrq, int *npoin,
1420        int *nr, int *nl, int *ifreq, int *itop, int *ipsh)
1421 {
1422     /* Local variables */
1423     static int itmp, ird, ipn, itp;
1424     double test1, test2;
1425
1426     /* Parameter adjustments */
1427     --nl;
1428     --nr;
1429     --npoin;
1430     --ifrq;
1431     --stp;
1432     --ipoin;
1433     --key;
1434
1435     /* Function Body */
1436     if (*ipsh) {
1437         /* Convert KVAL to int in range 1, ..., LDKEY. */
1438         ird = *kval % *ldkey + 1;
1439         /* Search for an unused location */
1440         for (itp = ird; itp <= *ldkey; ++itp) {
1441             if (key[itp] == *kval) {
1442                 goto L40;
1443             }
1444             if (key[itp] < 0) {
1445                 goto L30;
1446             }
1447         }
1448         for (itp = 1; itp <= ird - 1; ++itp) {
1449             if (key[itp] == *kval) {
1450                 goto L40;
1451             }
1452             if (key[itp] < 0) {
1453                 goto L30;
1454             }
1455         }
1456         /* Return if KEY array is full */
1457         /* KH
1458           prterr(6, "LDKEY is too small for this problem.\n"
1459           "It is not possible to estimate the value of LDKEY "
1460           "required,\n"
1461           "but twice the current value may be sufficient.");
1462           */
1463         prterr(6, "LDKEY is too small for this problem.\n"
1464                "Try increasing the size of the workspace.");
1465
1466         /* Update KEY */
1467 L30:
1468         key[itp] = *kval;
1469         ++(*itop);
1470         ipoin[itp] = *itop;
1471         /* Return if STP array full */
1472         if (*itop > *ldstp) {
1473             /* KH
1474                prterr(7, "LDSTP is too small for this problem.\n"
1475                "It is not possible to estimate the value of LDSTP "
1476                "required,\n"
1477                "but twice the current value may be sufficient.");
1478                */
1479             prterr(7, "LDSTP is too small for this problem.\n"
1480                    "Try increasing the size of the workspace.");
1481         }
1482         /* Update STP, etc. */
1483         npoin[*itop] = -1;
1484         nr[*itop] = -1;
1485         nl[*itop] = -1;
1486         stp[*itop] = *pastp;
1487         ifrq[*itop] = *ifreq;
1488         return;
1489     }
1490
1491     /* Find location, if any, of pastp */
1492 L40:
1493     ipn = ipoin[itp];
1494     test1 = *pastp - *tol;
1495     test2 = *pastp + *tol;
1496
1497 L50:
1498     if (stp[ipn] < test1) {
1499         ipn = nl[ipn];
1500         if (ipn > 0) {
1501             goto L50;
1502         }
1503     } else if (stp[ipn] > test2) {
1504         ipn = nr[ipn];
1505         if (ipn > 0) {
1506             goto L50;
1507         }
1508     } else {
1509         ifrq[ipn] += *ifreq;
1510         return;
1511     }
1512     /* Return if STP array full */
1513     ++(*itop);
1514     if (*itop > *ldstp) {
1515         /*
1516           prterr(7, "LDSTP is too small for this problem.\n"
1517           "It is not possible to estimate the value of LDSTP "
1518           "required,\n"
1519           "but twice the current value may be sufficient.");
1520           */
1521         prterr(7, "LDSTP is too small for this problem.\n"
1522                "Try increasing the size of the workspace.");
1523         return;
1524     }
1525     /* Find location to add value */
1526     ipn = ipoin[itp];
1527     itmp = ipn;
1528
1529 L60:
1530     if (stp[ipn] < test1) {
1531         itmp = ipn;
1532         ipn = nl[ipn];
1533         if (ipn > 0) {
1534             goto L60;
1535         } else {
1536             nl[itmp] = *itop;
1537         }
1538     } else if (stp[ipn] > test2) {
1539         itmp = ipn;
1540         ipn = nr[ipn];
1541         if (ipn > 0) {
1542             goto L60;
1543         } else {
1544             nr[itmp] = *itop;
1545         }
1546     }
1547     /* Update STP, etc. */
1548     npoin[*itop] = npoin[itmp];
1549     npoin[itmp] = *itop;
1550     stp[*itop] = *pastp;
1551     ifrq[*itop] = *ifreq;
1552     nl[*itop] = -1;
1553     nr[*itop] = -1;
1554 }
1555
1556 /*
1557   -----------------------------------------------------------------------
1558   Name:       F6XACT
1559   Purpose:    Pop a node off the stack.
1560   Usage:      CALL F6XACT (NROW, IROW, IFLAG, KYY, KEY, LDKEY, LAST, IPN)
1561   Arguments:
1562     NROW    - The number of rows in the table.                  (Input)
1563     IROW    - Vector of length nrow containing the row sums on
1564               output.                                           (Output)
1565     IFLAG   - Set to 3 if there are no additional nodes to process.
1566                                                                 (Output)
1567     KYY     - Constant mutlipliers used in forming the hash
1568               table key.                                        (Input)
1569     KEY     - Vector of length LDKEY containing the hash table
1570               keys.                                             (In/out)
1571     LDKEY   - Length of vector KEY.                             (Input)
1572     LAST    - Index of the last key popped off the stack.       (In/out)
1573     IPN     - Pointer to the linked list of past path lengths.  (Output)
1574   -----------------------------------------------------------------------
1575   */
1576 void
1577 f6xact(int *nrow, int *irow, int *iflag, int *kyy, int *key, int
1578        *ldkey, int *last, int *ipn)
1579 {
1580     int kval, j;
1581
1582     /* Parameter adjustments */
1583     --key;
1584     --kyy;
1585     --irow;
1586
1587     /* Function Body */
1588 L10:
1589     ++(*last);
1590     if (*last <= *ldkey) {
1591         if (key[*last] < 0) {
1592             goto L10;
1593         }
1594         /* Get KVAL from the stack */
1595         kval = key[*last];
1596         key[*last] = -9999;
1597         for (j = *nrow; j >= 2; --j) {
1598             irow[j] = kval / kyy[j];
1599             kval -= irow[j] * kyy[j];
1600         }
1601         irow[1] = kval;
1602         *ipn = *last;
1603     } else {
1604         *last = 0;
1605         *iflag = 3;
1606     }
1607     return;
1608 }
1609
1610 /*
1611   -----------------------------------------------------------------------
1612   Name:       F7XACT
1613   Purpose:    Generate the new nodes for given marinal totals.
1614   Usage:      CALL F7XACT (NROW, IMAX, IDIF, K, KS, IFLAG)
1615   Arguments:
1616     NROW    - The number of rows in the table.                  (Input)
1617     IMAX    - The row marginal totals.                          (Input)
1618     IDIF    - The column counts for the new column.             (in/out)
1619     K       - Indicator for the row to decrement.               (in/out)
1620     KS      - Indicator for the row to increment.               (in/out)
1621     IFLAG   - Status indicator.                                 (Output)
1622               If IFLAG is zero, a new table was generated.  For
1623               IFLAG = 1, no additional tables could be generated.
1624   -----------------------------------------------------------------------
1625   */
1626
1627 void
1628 f7xact(int *nrow, int *imax, int *idif, int *k, int *ks,
1629        int *iflag)
1630
1631 {
1632     int i, m, k1, mm;
1633
1634     /* Parameter adjustments */
1635     --idif;
1636     --imax;
1637
1638     /* Function Body */
1639     *iflag = 0;
1640     /* Find node which can be incremented, ks */
1641     if (*ks == 0)
1642         do {
1643             ++(*ks);
1644         } while (idif[*ks] == imax[*ks]);
1645
1646     /* Find node to decrement (>ks) */
1647     if (idif[*k] > 0 && *k > *ks) {
1648         --idif[*k];
1649         do {
1650             --(*k);
1651         } while(imax[*k] == 0);
1652
1653         m = *k;
1654
1655         /* Find node to increment (>=ks) */
1656         while (idif[m] >= imax[m]) {
1657             --m;
1658         }
1659         ++idif[m];
1660         /* Change ks */
1661         if (m == *ks) {
1662             if (idif[m] == imax[m]) {
1663                 *ks = *k;
1664             }
1665         }
1666     }
1667     else {
1668  Loop:
1669         /* Check for finish */
1670         for (k1 = *k + 1; k1 <= *nrow; ++k1) {
1671             if (idif[k1] > 0) {
1672                 goto L70;
1673             }
1674         }
1675         *iflag = 1;
1676         return;
1677
1678  L70:
1679         /* Reallocate counts */
1680         mm = 1;
1681         for (i = 1; i <= *k; ++i) {
1682             mm += idif[i];
1683             idif[i] = 0;
1684         }
1685         *k = k1;
1686
1687         do {
1688             --(*k);
1689             m = min(mm, imax[*k]);
1690             idif[*k] = m;
1691             mm -= m;
1692         } while (mm > 0 && *k != 1);
1693
1694         /* Check that all counts reallocated */
1695         if (mm > 0) {
1696             if (k1 != *nrow) {
1697                 *k = k1;
1698                 goto Loop;
1699             }
1700             *iflag = 1;
1701             return;
1702         }
1703         /* Get ks */
1704         --idif[k1];
1705         *ks = 0;
1706         do {
1707             ++(*ks);
1708             if (*ks > *k) {
1709                 return;
1710             }
1711         } while (idif[*ks] >= imax[*ks]);
1712     }
1713 }
1714
1715 /*
1716   -----------------------------------------------------------------------
1717   Name:       F8XACT
1718   Purpose:    Routine for reducing a vector when there is a zero
1719               element.
1720   Usage:      CALL F8XACT (IROW, IS, I1, IZERO, NEW)
1721   Arguments:
1722      IROW   - Vector containing the row counts.                 (Input)
1723      IS     - Indicator.                                        (Input)
1724      I1     - Indicator.                                        (Input)
1725      IZERO  - Position of the zero.                             (Input)
1726      NEW    - Vector of new row counts.                         (Output)
1727   -----------------------------------------------------------------------
1728   */
1729
1730 void
1731 f8xact(int *irow, int *is, int *i1, int *izero, int *myNew)
1732 {
1733     int i;
1734
1735     /* Parameter adjustments */
1736     --myNew;
1737     --irow;
1738
1739     /* Function Body */
1740     for (i = 1; i < *i1; ++i)
1741         myNew[i] = irow[i];
1742
1743     for (i = *i1; i <= *izero - 1; ++i) {
1744         if (*is >= irow[i + 1])
1745             break;
1746         myNew[i] = irow[i + 1];
1747     }
1748
1749     myNew[i] = *is;
1750
1751     for(;;) {
1752         ++i;
1753         if (i > *izero)
1754             return;
1755         myNew[i] = irow[i];
1756     }
1757 }
1758
1759 /*
1760   -----------------------------------------------------------------------
1761   Name:       F9XACT
1762   Purpose:    Computes the log of a multinomial coefficient.
1763   Usage:      F9XACT(N, MM, IR, FACT)
1764   Arguments:
1765      N      - Length of IR.                                     (Input)
1766      MM     - Number for factorial in numerator.                (Input)
1767      IR     - Vector of length N containing the numbers for
1768               the denominator of the factorial.                 (Input)
1769      FACT   - Table of log factorials.                          (Input)
1770      F9XACT - The log of the multinomal coefficient.            (Output)
1771   -----------------------------------------------------------------------
1772   */
1773
1774 double
1775 f9xact(int *n, int *mm, int *ir, double *fact)
1776 {
1777     double d;
1778     int k;
1779
1780     d = fact[*mm];
1781     for (k = 0; k < *n; k++)
1782         d -= fact[ir[k]];
1783     return d;
1784 }
1785
1786 /*
1787   -----------------------------------------------------------------------
1788   Name:     F10ACT
1789   Purpose:  Computes the shortest path length for special tables.
1790   Usage:    F10ACT (NROW, IROW, NCOL, ICOL, VAL, XMIN, FACT, ND, NE, M)
1791   Arguments:
1792      NROW   - The number of rows in the table.                  (Input)
1793      IROW   - Vector of length NROW containing the row totals.  (Input)
1794      NCOL   - The number of columns in the table.               (Input)
1795      ICO    - Vector of length NCOL containing the column totals.(Input)
1796      VAL    - The shortest path.                                (Output)
1797      XMIN   - Set to true if shortest path obtained.            (Output)
1798      FACT   - Vector containing the logarithms of factorials.   (Input)
1799      ND     - Workspace vector of length NROW.                  (Input)
1800      NE     - Workspace vector of length NCOL.                  (Input)
1801      M      - Workspace vector of length NCOL.                  (Input)
1802
1803   Chapter:    STAT/LIBRARY Categorical and Discrete Data Analysis
1804   -----------------------------------------------------------------------
1805   */
1806
1807 void
1808 f10act(int *nrow, int *irow, int *ncol, int *icol, double *val,
1809        int *xmin, double *fact, int *nd, int *ne, int *m)
1810 {
1811     /* Local variables */
1812     int i, is, ix, nrw1;
1813
1814     /* Parameter adjustments */
1815     --m;
1816     --ne;
1817     --nd;
1818     --icol;
1819     --irow;
1820
1821     /* Function Body */
1822     for (i = 1; i <= *nrow - 1; ++i)
1823         nd[i] = 0;
1824
1825     is = icol[1] / *nrow;
1826     ix = icol[1] - *nrow * is;
1827     ne[1] = is;
1828     m[1] = ix;
1829     if (ix != 0)
1830         ++nd[ix];
1831
1832     for (i = 2; i <= *ncol; ++i) {
1833         ix = icol[i] / *nrow;
1834         ne[i] = ix;
1835         is += ix;
1836         ix = icol[i] - *nrow * ix;
1837         m[i] = ix;
1838         if (ix != 0)
1839             ++nd[ix];
1840     }
1841
1842     for (i = *nrow - 2; i >= 1; --i)
1843         nd[i] += nd[i + 1];
1844
1845     ix = 0;
1846     nrw1 = *nrow + 1;
1847     for (i = *nrow; i >= 2; --i) {
1848         ix = ix + is + nd[nrw1 - i] - irow[i];
1849         if (ix < 0)
1850             return;
1851     }
1852
1853     for (i = 1; i <= *ncol; ++i) {
1854         ix = ne[i];
1855         is = m[i];
1856         *val = *val + is * fact[ix + 1] + (*nrow - is) * fact[ix];
1857     }
1858     *xmin = (1);
1859
1860     return;
1861 }
1862
1863 /*
1864   -----------------------------------------------------------------------
1865   Name:       F11ACT
1866   Purpose:    Routine for revising row totals.
1867   Usage:      CALL F11ACT (IROW, I1, I2, NEW)
1868   Arguments:
1869      IROW   - Vector containing the row totals. (Input)
1870      I1     - Indicator.                        (Input)
1871      I2     - Indicator.                        (Input)
1872      NEW    - Vector containing the row totals. (Output)
1873   -----------------------------------------------------------------------
1874   */
1875 void
1876 f11act(int *irow, int *i1, int *i2, int *myNew)
1877 {
1878     int i;
1879
1880     /* Parameter adjustments */
1881     --myNew;
1882     --irow;
1883
1884     for (i = 1; i <= (*i1 - 1); ++i)    myNew[i] = irow[i];
1885     for (i = *i1; i <= *i2; ++i)        myNew[i] = irow[i + 1];
1886
1887     return;
1888 }
1889
1890 /*
1891   -----------------------------------------------------------------------
1892   Name:       prterr
1893   Purpose:    Print an error message and stop.
1894   Usage:      prterr(icode, mes)
1895   Arguments:
1896      icode  - Integer code for the error message.               (Input)
1897      mes    - Character string containing the error message.    (Input)
1898   -----------------------------------------------------------------------
1899   */
1900 void
1901 prterr(int icode, char *mes)
1902 {
1903 //    PROBLEM "FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY);
1904 //   printf("FEXACT error %d.\n%s", icode, mes RECOVER(NULL_ENTRY));
1905    printf("FEXACT error %d.\n", icode); //xx:another one of my ugly kluges
1906    return;
1907 }
1908
1909 /*
1910   -----------------------------------------------------------------------
1911   Name:       iwork
1912   Purpose:    Routine for allocating workspace.
1913   Usage:      iwork (iwkmax, iwkpt, number, itype)
1914   Arguments:
1915      iwkmax - Maximum (int) amount of workspace.                (Input)
1916      iwkpt  - Amount of (int) workspace currently allocated.    (in/out)
1917      number - Number of elements of workspace desired.          (Input)
1918      itype  - Workspace type.                                   (Input)
1919               ITYPE  TYPE
1920                 2    integer
1921                 3    float
1922                 4    double
1923      iwork(): Index in rwrk, dwrk, or iwrk of the beginning of
1924               the first free element in the workspace array.    (Output)
1925   -----------------------------------------------------------------------
1926   */
1927 int
1928 iwork(int iwkmax, int *iwkpt, int number, int itype)
1929 {
1930     int i;
1931
1932     i = *iwkpt;
1933     if (itype == 2 || itype == 3)
1934         *iwkpt += number;
1935     else { /* double */
1936         if (i % 2 != 0)
1937             ++i;
1938         *iwkpt += (number << 1);
1939         i /= 2;
1940     }
1941     if (*iwkpt > iwkmax)
1942         prterr(40, "Out of workspace.");
1943
1944     return i;
1945 }
1946
1947 #ifndef USING_R
1948
1949 void isort(int *n, int *ix)
1950 {
1951 /*
1952   -----------------------------------------------------------------------
1953   Name:       ISORT
1954   Purpose:    Shell sort for an int vector.
1955   Usage:      CALL ISORT (N, IX)
1956   Arguments:
1957      N      - Lenth of vector IX.       (Input)
1958      IX     - Vector to be sorted.      (in/out)
1959   -----------------------------------------------------------------------
1960   */
1961     static int ikey, i, j, m, il[10], kl, it, iu[10], ku;
1962
1963     /* Parameter adjustments */
1964     --ix;
1965
1966     /* Function Body */
1967     m = 1;
1968     i = 1;
1969     j = *n;
1970
1971 L10:
1972     if (i >= j) {
1973         goto L40;
1974     }
1975     kl = i;
1976     ku = j;
1977     ikey = i;
1978     ++j;
1979     /* Find element in first half */
1980 L20:
1981     ++i;
1982     if (i < j) {
1983         if (ix[ikey] > ix[i]) {
1984             goto L20;
1985         }
1986     }
1987     /* Find element in second half */
1988 L30:
1989     --j;
1990     if (ix[j] > ix[ikey]) {
1991         goto L30;
1992     }
1993     /* Interchange */
1994     if (i < j) {
1995         it = ix[i];
1996         ix[i] = ix[j];
1997         ix[j] = it;
1998         goto L20;
1999     }
2000     it = ix[ikey];
2001     ix[ikey] = ix[j];
2002     ix[j] = it;
2003     /* Save upper and lower subscripts of the array yet to be sorted */
2004     if (m < 11) {
2005         if (j - kl < ku - j) {
2006             il[m - 1] = j + 1;
2007             iu[m - 1] = ku;
2008             i = kl;
2009             --j;
2010         } else {
2011             il[m - 1] = kl;
2012             iu[m - 1] = j - 1;
2013             i = j + 1;
2014             j = ku;
2015         }
2016         ++m;
2017         goto L10;
2018     } else {
2019         prterr(20, "This should never occur.");
2020     }
2021     /* Use another segment */
2022 L40:
2023     --m;
2024     if (m == 0) {
2025         return;
2026     }
2027     i = il[m - 1];
2028     j = iu[m - 1];
2029     goto L10;
2030 }
2031
2032 double gammds(double *y, double *p, int *ifault)
2033 {
2034 /*
2035   -----------------------------------------------------------------------
2036   Name:       GAMMDS
2037   Purpose:    Cumulative distribution for the gamma distribution.
2038   Usage:      PGAMMA (Q, ALPHA,IFAULT)
2039   Arguments:
2040      Q      - Value at which the distribution is desired.  (Input)
2041      ALPHA  - Parameter in the gamma distribution.  (Input)
2042      IFAULT - Error indicator.  (Output)
2043                IFAULT  DEFINITION
2044                  0     No error
2045                  1     An argument is misspecified.
2046                  2     A numerical error has occurred.
2047      PGAMMA - The cdf for the gamma distribution with parameter alpha
2048               evaluated at Q.  (Output)
2049   -----------------------------------------------------------------------
2050
2051   Algorithm AS 147 APPL. Statist. (1980) VOL. 29, P. 113
2052
2053   Computes the incomplete gamma integral for positive parameters Y, P
2054   using and infinite series.
2055   */
2056
2057     static double a, c, f, g;
2058     static int ifail;
2059
2060     /* Checks for the admissibility of arguments and value of F */
2061     *ifault = 1;
2062     g = 0.;
2063     if (*y <= 0. || *p <= 0.) {
2064         return g;
2065     }
2066     *ifault = 2;
2067
2068     /*
2069       ALOGAM is natural log of gamma function no need to test ifail as
2070       an error is impossible
2071       */
2072
2073     a = *p + 1.;
2074     f = exp(*p * log(*y) - alogam(&a, &ifail) - *y);
2075     if (f == 0.) {
2076         return g;
2077     }
2078     *ifault = 0;
2079
2080     /* Series begins */
2081     c = 1.;
2082     g = 1.;
2083     a = *p;
2084 L10:
2085     a += 1.;
2086     c = c * *y / a;
2087     g += c;
2088     if (c / g > 1e-6) {
2089         goto L10;
2090     }
2091     g *= f;
2092     return g;
2093 }
2094
2095 /*
2096   -----------------------------------------------------------------------
2097   Name:       ALOGAM
2098   Purpose:    Value of the log-gamma function.
2099   Usage:      ALOGAM (X, IFAULT)
2100   Arguments:
2101      X      - Value at which the log-gamma function is to be evaluated.
2102               (Input)
2103      IFAULT  - Error indicator.  (Output)
2104                IFAULT  DEFINITION
2105                  0     No error
2106                  1     X < 0
2107      ALGAMA - The value of the log-gamma function at XX.  (Output)
2108   -----------------------------------------------------------------------
2109
2110   Algorithm ACM 291, Comm. ACM. (1966) Vol. 9, P. 684
2111
2112   Evaluates natural logarithm of gamma(x) for X greater than zero.
2113   */
2114
2115 double alogam(double *x, int *ifault)
2116 {
2117     /* Initialized data */
2118         //printf("alogam x = %f\t%d\n",*x,*ifault);
2119     static double a1 = .918938533204673;
2120     static double a2 = 5.95238095238e-4;
2121     static double a3 = 7.93650793651e-4;
2122     static double a4 = .002777777777778;
2123     static double a5 = .083333333333333;
2124
2125     /* Local variables */
2126     static double f, y, z;
2127
2128     *ifault = 1;
2129     if (*x < 0.) {
2130         return(0.);
2131     }
2132     *ifault = 0;
2133     y = *x;
2134     f = 0.;
2135     if (y >= 7.) {
2136         goto L30;
2137     }
2138     f = y;
2139 L10:
2140     y += 1.;
2141     if (y >= 7.) {
2142         goto L20;
2143     }
2144     f *= y;
2145     goto L10;
2146 L20:
2147     f = -log(f);
2148 L30:
2149     z = 1. / (y * y);
2150
2151         //printf("returning %f\n",(f + (y - .5) * log(y) - y + a1 + (((-a2 * z + a3) * z - a4) * z + a5) / y));
2152     return(f + (y - .5) * log(y) - y + a1 +
2153            (((-a2 * z + a3) * z - a4) * z + a5) / y);
2154 }
2155
2156
2157 #endif /* not USING_R */
2158