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