fisher2.c

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