src/language/stats/crosstabs.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000, 2006 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 /* FIXME:
  18
  19    - Pearson's R (but not Spearman!) is off a little.
  20    - T values for Spearman's R and Pearson's R are wrong.
  21    - How to calculate significance of symmetric and directional measures?
  22    - Asymmetric ASEs and T values for lambda are wrong.
  23    - ASE of Goodman and Kruskal's tau is not calculated.
  24    - ASE of symmetric somers' d is wrong.
  25    - Approx. T of uncertainty coefficient is wrong.
  26
  27 */
  28
  29 #include <config.h>
  30
  31 #include <ctype.h>
  32 #include <gsl/gsl_cdf.h>
  33 #include <stdlib.h>
  34 #include <stdio.h>
  35
  36 #include <data/case.h>
  37 #include <data/casegrouper.h>
  38 #include <data/casereader.h>
  39 #include <data/data-out.h>
  40 #include <data/dictionary.h>
  41 #include <data/format.h>
  42 #include <data/procedure.h>
  43 #include <data/value-labels.h>
  44 #include <data/variable.h>
  45 #include <language/command.h>
  46 #include <language/dictionary/split-file.h>
  47 #include <language/lexer/lexer.h>
  48 #include <language/lexer/variable-parser.h>
  49 #include <libpspp/array.h>
  50 #include <libpspp/assertion.h>
  51 #include <libpspp/compiler.h>
  52 #include <libpspp/hash.h>
  53 #include <libpspp/message.h>
  54 #include <libpspp/misc.h>
  55 #include <libpspp/pool.h>
  56 #include <libpspp/str.h>
  57 #include <output/output.h>
  58 #include <output/table.h>
  59
  60 #include "minmax.h"
  61 #include "xalloc.h"
  62 #include "xmalloca.h"
  63
  64 #include "gettext.h"
  65 #define _(msgid) gettext (msgid)
  66 #define N_(msgid) msgid
  67
  68 /* (headers) */
  69
  70 /* (specification)
  71    crosstabs (crs_):
  72      *^tables=custom;
  73      +variables=custom;
  74      missing=miss:!table/include/report;
  75      +write[wr_]=none,cells,all;
  76      +format=fmt:!labels/nolabels/novallabs,
  77              val:!avalue/dvalue,
  78              indx:!noindex/index,
  79              tabl:!tables/notables,
  80              box:!box/nobox,
  81              pivot:!pivot/nopivot;
  82      +cells[cl_]=count,none,expected,row,column,total,residual,sresidual,
  83                  asresidual,all;
  84      +statistics[st_]=chisq,phi,cc,lambda,uc,none,btau,ctau,risk,gamma,d,
  85                       kappa,eta,corr,all.
  86 */
  87 /* (declarations) */
  88 /* (functions) */
  89
  90 /* Number of chi-square statistics. */
  91 #define N_CHISQ 5
  92
  93 /* Number of symmetric statistics. */
  94 #define N_SYMMETRIC 9
  95
  96 /* Number of directional statistics. */
  97 #define N_DIRECTIONAL 13
  98
  99 /* A single table entry for general mode. */
 100 struct table_entry
 101   {
 102     int table;          /* Flattened table number. */
 103     union
 104       {
 105         double freq;    /* Frequency count. */
 106         double *data;   /* Crosstabulation table for integer mode. */
 107       }
 108     u;
 109     union value values[1];      /* Values. */
 110   };
 111
 112 /* A crosstabulation. */
 113 struct crosstab
 114   {
 115     int nvar;                   /* Number of variables. */
 116     double missing;             /* Missing cases count. */
 117     int ofs;                    /* Integer mode: Offset into sorted_tab[]. */
 118     const struct variable *vars[2];     /* At least two variables; sorted by
 119                                    larger indices first. */
 120   };
 121
 122 /* Integer mode variable info. */
 123 struct var_range
 124   {
 125     int min;                    /* Minimum value. */
 126     int max;                    /* Maximum value + 1. */
 127     int count;                  /* max - min. */
 128   };
 129
 130 static inline struct var_range *
 131 get_var_range (const struct variable *v)
 132 {
 133   return var_get_aux (v);
 134 }
 135
 136 /* Indexes into crosstab.v. */
 137 enum
 138   {
 139     ROW_VAR = 0,
 140     COL_VAR = 1
 141   };
 142
 143 /* General mode crosstabulation table. */
 144 static struct hsh_table *gen_tab;       /* Hash table. */
 145 static int n_sorted_tab;                /* Number of entries in sorted_tab. */
 146 static struct table_entry **sorted_tab; /* Sorted table. */
 147
 148 /* Variables specifies on VARIABLES. */
 149 static const struct variable **variables;
 150 static size_t variables_cnt;
 151
 152 /* TABLES. */
 153 static struct crosstab **xtab;
 154 static int nxtab;
 155
 156 /* Integer or general mode? */
 157 enum
 158   {
 159     INTEGER,
 160     GENERAL
 161   };
 162 static int mode;
 163
 164 /* CELLS. */
 165 static int num_cells;           /* Number of cells requested. */
 166 static int cells[8];            /* Cells requested. */
 167
 168 /* WRITE. */
 169 static int write_style;         /* One of WR_* that specifies the WRITE style. */
 170
 171 /* Command parsing info. */
 172 static struct cmd_crosstabs cmd;
 173
 174 /* Pools. */
 175 static struct pool *pl_tc;      /* For table cells. */
 176 static struct pool *pl_col;     /* For column data. */
 177
 178 static int internal_cmd_crosstabs (struct lexer *lexer, struct dataset *ds);
 179 static void precalc (struct casereader *, const struct dataset *);
 180 static void calc_general (struct ccase *, const struct dataset *);
 181 static void calc_integer (struct ccase *, const struct dataset *);
 182 static void postcalc (void);
 183 static void submit (struct tab_table *);
 184
 185 static void format_short (char *s, const struct fmt_spec *fp,
 186                          const union value *v);
 187
 188 /* Parse and execute CROSSTABS, then clean up. */
 189 int
 190 cmd_crosstabs (struct lexer *lexer, struct dataset *ds)
 191 {
 192   int result = internal_cmd_crosstabs (lexer, ds);
 193   int i;
 194
 195   free (variables);
 196   pool_destroy (pl_tc);
 197   pool_destroy (pl_col);
 198
 199   for (i = 0; i < nxtab; i++)
 200     free (xtab[i]);
 201   free (xtab);
 202
 203   return result;
 204 }
 205
 206 /* Parses and executes the CROSSTABS procedure. */
 207 static int
 208 internal_cmd_crosstabs (struct lexer *lexer, struct dataset *ds)
 209 {
 210   struct casegrouper *grouper;
 211   struct casereader *input, *group;
 212   bool ok;
 213   int i;
 214
 215   variables = NULL;
 216   variables_cnt = 0;
 217   xtab = NULL;
 218   nxtab = 0;
 219   pl_tc = pool_create ();
 220   pl_col = pool_create ();
 221
 222   if (!parse_crosstabs (lexer, ds, &cmd, NULL))
 223     return CMD_FAILURE;
 224
 225   mode = variables ? INTEGER : GENERAL;
 226
 227   /* CELLS. */
 228   if (!cmd.sbc_cells)
 229     {
 230       cmd.a_cells[CRS_CL_COUNT] = 1;
 231     }
 232   else
 233     {
 234       int count = 0;
 235
 236       for (i = 0; i < CRS_CL_count; i++)
 237         if (cmd.a_cells[i])
 238           count++;
 239       if (count == 0)
 240         {
 241           cmd.a_cells[CRS_CL_COUNT] = 1;
 242           cmd.a_cells[CRS_CL_ROW] = 1;
 243           cmd.a_cells[CRS_CL_COLUMN] = 1;
 244           cmd.a_cells[CRS_CL_TOTAL] = 1;
 245         }
 246       if (cmd.a_cells[CRS_CL_ALL])
 247         {
 248           for (i = 0; i < CRS_CL_count; i++)
 249             cmd.a_cells[i] = 1;
 250           cmd.a_cells[CRS_CL_ALL] = 0;
 251         }
 252       cmd.a_cells[CRS_CL_NONE] = 0;
 253     }
 254   for (num_cells = i = 0; i < CRS_CL_count; i++)
 255     if (cmd.a_cells[i])
 256       cells[num_cells++] = i;
 257
 258   /* STATISTICS. */
 259   if (cmd.sbc_statistics)
 260     {
 261       int i;
 262       int count = 0;
 263
 264       for (i = 0; i < CRS_ST_count; i++)
 265         if (cmd.a_statistics[i])
 266           count++;
 267       if (count == 0)
 268         cmd.a_statistics[CRS_ST_CHISQ] = 1;
 269       if (cmd.a_statistics[CRS_ST_ALL])
 270         for (i = 0; i < CRS_ST_count; i++)
 271           cmd.a_statistics[i] = 1;
 272     }
 273
 274   /* MISSING. */
 275   if (cmd.miss == CRS_REPORT && mode == GENERAL)
 276     {
 277       msg (SE, _("Missing mode REPORT not allowed in general mode.  "
 278                  "Assuming MISSING=TABLE."));
 279       cmd.miss = CRS_TABLE;
 280     }
 281
 282   /* WRITE. */
 283   if (cmd.a_write[CRS_WR_ALL] && cmd.a_write[CRS_WR_CELLS])
 284     cmd.a_write[CRS_WR_ALL] = 0;
 285   if (cmd.a_write[CRS_WR_ALL] && mode == GENERAL)
 286     {
 287       msg (SE, _("Write mode ALL not allowed in general mode.  "
 288                  "Assuming WRITE=CELLS."));
 289       cmd.a_write[CRS_WR_CELLS] = 1;
 290     }
 291   if (cmd.sbc_write
 292       && (cmd.a_write[CRS_WR_NONE]
 293           + cmd.a_write[CRS_WR_ALL]
 294           + cmd.a_write[CRS_WR_CELLS] == 0))
 295     cmd.a_write[CRS_WR_CELLS] = 1;
 296   if (cmd.a_write[CRS_WR_CELLS])
 297     write_style = CRS_WR_CELLS;
 298   else if (cmd.a_write[CRS_WR_ALL])
 299     write_style = CRS_WR_ALL;
 300   else
 301     write_style = CRS_WR_NONE;
 302
 303   input = casereader_create_filter_weight (proc_open (ds), dataset_dict (ds),
 304                                            NULL, NULL);
 305   grouper = casegrouper_create_splits (input, dataset_dict (ds));
 306   while (casegrouper_get_next_group (grouper, &group))
 307     {
 308       struct ccase c;
 309
 310       precalc (group, ds);
 311
 312       for (; casereader_read (group, &c); case_destroy (&c))
 313         {
 314           if (mode == GENERAL)
 315             calc_general (&c, ds);
 316           else
 317             calc_integer (&c, ds);
 318         }
 319       casereader_destroy (group);
 320
 321       postcalc ();
 322     }
 323   ok = casegrouper_destroy (grouper);
 324   ok = proc_commit (ds) && ok;
 325
 326   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 327 }
 328
 329 /* Parses the TABLES subcommand. */
 330 static int
 331 crs_custom_tables (struct lexer *lexer, struct dataset *ds, struct cmd_crosstabs *cmd UNUSED, void *aux UNUSED)
 332 {
 333   struct const_var_set *var_set;
 334   int n_by;
 335   const struct variable ***by = NULL;
 336   size_t *by_nvar = NULL;
 337   size_t nx = 1;
 338   int success = 0;
 339
 340   /* Ensure that this is a TABLES subcommand. */
 341   if (!lex_match_id (lexer, "TABLES")
 342       && (lex_token (lexer) != T_ID ||
 343           dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) == NULL)
 344       && lex_token (lexer) != T_ALL)
 345     return 2;
 346   lex_match (lexer, '=');
 347
 348   if (variables != NULL)
 349     var_set = const_var_set_create_from_array (variables, variables_cnt);
 350   else
 351     var_set = const_var_set_create_from_dict (dataset_dict (ds));
 352   assert (var_set != NULL);
 353
 354   for (n_by = 0; ;)
 355     {
 356       by = xnrealloc (by, n_by + 1, sizeof *by);
 357       by_nvar = xnrealloc (by_nvar, n_by + 1, sizeof *by_nvar);
 358       if (!parse_const_var_set_vars (lexer, var_set, &by[n_by], &by_nvar[n_by],
 359                                PV_NO_DUPLICATE | PV_NO_SCRATCH))
 360         goto done;
 361       if (xalloc_oversized (nx, by_nvar[n_by]))
 362         {
 363           msg (SE, _("Too many cross-tabulation variables or dimensions."));
 364           goto done;
 365         }
 366       nx *= by_nvar[n_by];
 367       n_by++;
 368
 369       if (!lex_match (lexer, T_BY))
 370         {
 371           if (n_by < 2)
 372             {
 373               lex_error (lexer, _("expecting BY"));
 374               goto done;
 375             }
 376           else
 377             break;
 378         }
 379     }
 380
 381   {
 382     int *by_iter = xcalloc (n_by, sizeof *by_iter);
 383     int i;
 384
 385     xtab = xnrealloc (xtab, nxtab + nx, sizeof *xtab);
 386     for (i = 0; i < nx; i++)
 387       {
 388         struct crosstab *x;
 389
 390         x = xmalloc (sizeof *x + sizeof (struct variable *) * (n_by - 2));
 391         x->nvar = n_by;
 392         x->missing = 0.;
 393
 394         {
 395           int i;
 396
 397           for (i = 0; i < n_by; i++)
 398             x->vars[i] = by[i][by_iter[i]];
 399         }
 400
 401         {
 402           int i;
 403
 404           for (i = n_by - 1; i >= 0; i--)
 405             {
 406               if (++by_iter[i] < by_nvar[i])
 407                 break;
 408               by_iter[i] = 0;
 409             }
 410         }
 411
 412         xtab[nxtab++] = x;
 413       }
 414     free (by_iter);
 415   }
 416   success = 1;
 417
 418  done:
 419   /* All return paths lead here. */
 420   {
 421     int i;
 422
 423     for (i = 0; i < n_by; i++)
 424       free (by[i]);
 425     free (by);
 426     free (by_nvar);
 427   }
 428
 429   const_var_set_destroy (var_set);
 430
 431   return success;
 432 }
 433
 434 /* Parses the VARIABLES subcommand. */
 435 static int
 436 crs_custom_variables (struct lexer *lexer, struct dataset *ds, struct cmd_crosstabs *cmd UNUSED, void *aux UNUSED)
 437 {
 438   if (nxtab)
 439     {
 440       msg (SE, _("VARIABLES must be specified before TABLES."));
 441       return 0;
 442     }
 443
 444   lex_match (lexer, '=');
 445
 446   for (;;)
 447     {
 448       size_t orig_nv = variables_cnt;
 449       size_t i;
 450
 451       long min, max;
 452
 453       if (!parse_variables_const (lexer, dataset_dict (ds),
 454                             &variables, &variables_cnt,
 455                             (PV_APPEND | PV_NUMERIC
 456                              | PV_NO_DUPLICATE | PV_NO_SCRATCH)))
 457         return 0;
 458
 459       if (lex_token (lexer) != '(')
 460         {
 461           lex_error (lexer, "expecting `('");
 462           goto lossage;
 463         }
 464       lex_get (lexer);
 465
 466       if (!lex_force_int (lexer))
 467         goto lossage;
 468       min = lex_integer (lexer);
 469       lex_get (lexer);
 470
 471       lex_match (lexer, ',');
 472
 473       if (!lex_force_int (lexer))
 474         goto lossage;
 475       max = lex_integer (lexer);
 476       if (max < min)
 477         {
 478           msg (SE, _("Maximum value (%ld) less than minimum value (%ld)."),
 479                max, min);
 480           goto lossage;
 481         }
 482       lex_get (lexer);
 483
 484       if (lex_token (lexer) != ')')
 485         {
 486           lex_error (lexer, "expecting `)'");
 487           goto lossage;
 488         }
 489       lex_get (lexer);
 490
 491       for (i = orig_nv; i < variables_cnt; i++)
 492         {
 493           struct var_range *vr = xmalloc (sizeof *vr);
 494           vr->min = min;
 495           vr->max = max + 1.;
 496           vr->count = max - min + 1;
 497           var_attach_aux (variables[i], vr, var_dtor_free);
 498         }
 499
 500       if (lex_token (lexer) == '/')
 501         break;
 502     }
 503
 504   return 1;
 505
 506  lossage:
 507   free (variables);
 508   variables = NULL;
 509   return 0;
 510 }
 511 \f
 512 /* Data file processing. */
 513
 514 static int compare_table_entry (const void *, const void *, const void *);
 515 static unsigned hash_table_entry (const void *, const void *);
 516
 517 /* Set up the crosstabulation tables for processing. */
 518 static void
 519 precalc (struct casereader *input, const struct dataset *ds)
 520 {
 521   struct ccase c;
 522
 523   if (casereader_peek (input, 0, &c))
 524     {
 525       output_split_file_values (ds, &c);
 526       case_destroy (&c);
 527     }
 528
 529   if (mode == GENERAL)
 530     {
 531       gen_tab = hsh_create (512, compare_table_entry, hash_table_entry,
 532                             NULL, NULL);
 533     }
 534   else
 535     {
 536       int i;
 537
 538       sorted_tab = NULL;
 539       n_sorted_tab = 0;
 540
 541       for (i = 0; i < nxtab; i++)
 542         {
 543           struct crosstab *x = xtab[i];
 544           int count = 1;
 545           int *v;
 546           int j;
 547
 548           x->ofs = n_sorted_tab;
 549
 550           for (j = 2; j < x->nvar; j++)
 551             count *= get_var_range (x->vars[j - 2])->count;
 552
 553           sorted_tab = xnrealloc (sorted_tab,
 554                                   n_sorted_tab + count, sizeof *sorted_tab);
 555           v = xmalloca (sizeof *v * x->nvar);
 556           for (j = 2; j < x->nvar; j++)
 557             v[j] = get_var_range (x->vars[j])->min;
 558           for (j = 0; j < count; j++)
 559             {
 560               struct table_entry *te;
 561               int k;
 562
 563               te = sorted_tab[n_sorted_tab++]
 564                 = xmalloc (sizeof *te + sizeof (union value) * (x->nvar - 1));
 565               te->table = i;
 566
 567               {
 568                 int row_cnt = get_var_range (x->vars[0])->count;
 569                 int col_cnt = get_var_range (x->vars[1])->count;
 570                 const int mat_size = row_cnt * col_cnt;
 571                 int m;
 572
 573                 te->u.data = xnmalloc (mat_size, sizeof *te->u.data);
 574                 for (m = 0; m < mat_size; m++)
 575                   te->u.data[m] = 0.;
 576               }
 577
 578               for (k = 2; k < x->nvar; k++)
 579                 te->values[k].f = v[k];
 580               for (k = 2; k < x->nvar; k++)
 581                 {
 582                   struct var_range *vr = get_var_range (x->vars[k]);
 583                   if (++v[k] >= vr->max)
 584                     v[k] = vr->min;
 585                   else
 586                     break;
 587                 }
 588             }
 589           freea (v);
 590         }
 591
 592       sorted_tab = xnrealloc (sorted_tab,
 593                               n_sorted_tab + 1, sizeof *sorted_tab);
 594       sorted_tab[n_sorted_tab] = NULL;
 595     }
 596
 597 }
 598
 599 /* Form crosstabulations for general mode. */
 600 static void
 601 calc_general (struct ccase *c, const struct dataset *ds)
 602 {
 603   /* Missing values to exclude. */
 604   enum mv_class exclude = (cmd.miss == CRS_TABLE ? MV_ANY
 605                            : cmd.miss == CRS_INCLUDE ? MV_SYSTEM
 606                            : MV_NEVER);
 607
 608   /* Case weight. */
 609   double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 610
 611   /* Flattened current table index. */
 612   int t;
 613
 614   for (t = 0; t < nxtab; t++)
 615     {
 616       struct crosstab *x = xtab[t];
 617       const size_t entry_size = (sizeof (struct table_entry)
 618                                  + sizeof (union value) * (x->nvar - 1));
 619       struct table_entry *te = xmalloca (entry_size);
 620
 621       /* Construct table entry for the current record and table. */
 622       te->table = t;
 623       {
 624         int j;
 625
 626         assert (x != NULL);
 627         for (j = 0; j < x->nvar; j++)
 628           {
 629             const union value *v = case_data (c, x->vars[j]);
 630             if (var_is_value_missing (x->vars[j], v, exclude))
 631               {
 632                 x->missing += weight;
 633                 goto next_crosstab;
 634               }
 635
 636             if (var_is_numeric (x->vars[j]))
 637               te->values[j].f = case_num (c, x->vars[j]);
 638             else
 639               {
 640                 memcpy (te->values[j].s, case_str (c, x->vars[j]),
 641                         var_get_width (x->vars[j]));
 642
 643                 /* Necessary in order to simplify comparisons. */
 644                 memset (&te->values[j].s[var_get_width (x->vars[j])], 0,
 645                         sizeof (union value) - var_get_width (x->vars[j]));
 646               }
 647           }
 648       }
 649
 650       /* Add record to hash table. */
 651       {
 652         struct table_entry **tepp
 653           = (struct table_entry **) hsh_probe (gen_tab, te);
 654         if (*tepp == NULL)
 655           {
 656             struct table_entry *tep = pool_alloc (pl_tc, entry_size);
 657
 658             te->u.freq = weight;
 659             memcpy (tep, te, entry_size);
 660
 661             *tepp = tep;
 662           }
 663         else
 664           (*tepp)->u.freq += weight;
 665       }
 666
 667     next_crosstab:
 668       freea (te);
 669     }
 670 }
 671
 672 static void
 673 calc_integer (struct ccase *c, const struct dataset *ds)
 674 {
 675   bool bad_warn = true;
 676
 677   /* Case weight. */
 678   double weight = dict_get_case_weight (dataset_dict (ds), c, &bad_warn);
 679
 680   /* Flattened current table index. */
 681   int t;
 682
 683   for (t = 0; t < nxtab; t++)
 684     {
 685       struct crosstab *x = xtab[t];
 686       int i, fact, ofs;
 687
 688       fact = i = 1;
 689       ofs = x->ofs;
 690       for (i = 0; i < x->nvar; i++)
 691         {
 692           const struct variable *const v = x->vars[i];
 693           struct var_range *vr = get_var_range (v);
 694           double value = case_num (c, v);
 695
 696           /* Note that the first test also rules out SYSMIS. */
 697           if ((value < vr->min || value >= vr->max)
 698               || (cmd.miss == CRS_TABLE
 699                   && var_is_num_missing (v, value, MV_USER)))
 700             {
 701               x->missing += weight;
 702               goto next_crosstab;
 703             }
 704
 705           if (i > 1)
 706             {
 707               ofs += fact * ((int) value - vr->min);
 708               fact *= vr->count;
 709             }
 710         }
 711
 712       {
 713         const struct variable *row_var = x->vars[ROW_VAR];
 714         const int row = case_num (c, row_var) - get_var_range (row_var)->min;
 715
 716         const struct variable *col_var = x->vars[COL_VAR];
 717         const int col = case_num (c, col_var) - get_var_range (col_var)->min;
 718
 719         const int col_dim = get_var_range (col_var)->count;
 720
 721         sorted_tab[ofs]->u.data[col + row * col_dim] += weight;
 722       }
 723
 724     next_crosstab: ;
 725     }
 726 }
 727
 728 /* Compare the table_entry's at A and B and return a strcmp()-type
 729    result. */
 730 static int
 731 compare_table_entry (const void *a_, const void *b_, const void *aux UNUSED)
 732 {
 733   const struct table_entry *a = a_;
 734   const struct table_entry *b = b_;
 735
 736   if (a->table > b->table)
 737     return 1;
 738   else if (a->table < b->table)
 739     return -1;
 740
 741   {
 742     const struct crosstab *x = xtab[a->table];
 743     int i;
 744
 745     for (i = x->nvar - 1; i >= 0; i--)
 746       if (var_is_numeric (x->vars[i]))
 747         {
 748           const double diffnum = a->values[i].f - b->values[i].f;
 749           if (diffnum < 0)
 750             return -1;
 751           else if (diffnum > 0)
 752             return 1;
 753         }
 754       else
 755         {
 756           const int diffstr = strncmp (a->values[i].s, b->values[i].s,
 757                                        var_get_width (x->vars[i]));
 758           if (diffstr)
 759             return diffstr;
 760         }
 761   }
 762
 763   return 0;
 764 }
 765
 766 /* Calculate a hash value from table_entry A. */
 767 static unsigned
 768 hash_table_entry (const void *a_, const void *aux UNUSED)
 769 {
 770   const struct table_entry *a = a_;
 771   unsigned long hash;
 772   int i;
 773
 774   hash = a->table;
 775   for (i = 0; i < xtab[a->table]->nvar; i++)
 776     hash ^= hsh_hash_bytes (&a->values[i], sizeof a->values[i]);
 777
 778   return hash;
 779 }
 780 \f
 781 /* Post-data reading calculations. */
 782
 783 static struct table_entry **find_pivot_extent (struct table_entry **,
 784                                                int *cnt, int pivot);
 785 static void enum_var_values (struct table_entry **entries, int entry_cnt,
 786                              int var_idx,
 787                              union value **values, int *value_cnt);
 788 static void output_pivot_table (struct table_entry **, struct table_entry **,
 789                                 double **, double **, double **,
 790                                 int *, int *, int *);
 791 static void make_summary_table (void);
 792
 793 static void
 794 postcalc (void)
 795 {
 796   if (mode == GENERAL)
 797     {
 798       n_sorted_tab = hsh_count (gen_tab);
 799       sorted_tab = (struct table_entry **) hsh_sort (gen_tab);
 800     }
 801
 802   make_summary_table ();
 803
 804   /* Identify all the individual crosstabulation tables, and deal with
 805      them. */
 806   {
 807     struct table_entry **pb = sorted_tab, **pe; /* Pivot begin, pivot end. */
 808     int pc = n_sorted_tab;                      /* Pivot count. */
 809
 810     double *mat = NULL, *row_tot = NULL, *col_tot = NULL;
 811     int maxrows = 0, maxcols = 0, maxcells = 0;
 812
 813     for (;;)
 814       {
 815         pe = find_pivot_extent (pb, &pc, cmd.pivot == CRS_PIVOT);
 816         if (pe == NULL)
 817           break;
 818
 819         output_pivot_table (pb, pe, &mat, &row_tot, &col_tot,
 820                             &maxrows, &maxcols, &maxcells);
 821
 822         pb = pe;
 823       }
 824     free (mat);
 825     free (row_tot);
 826     free (col_tot);
 827   }
 828
 829   hsh_destroy (gen_tab);
 830   if (mode == INTEGER)
 831     {
 832       int i;
 833       for (i = 0; i < n_sorted_tab; i++)
 834         {
 835           free (sorted_tab[i]->u.data);
 836           free (sorted_tab[i]);
 837         }
 838       free (sorted_tab);
 839     }
 840 }
 841
 842 static void insert_summary (struct tab_table *, int tab_index, double valid);
 843
 844 /* Output a table summarizing the cases processed. */
 845 static void
 846 make_summary_table (void)
 847 {
 848   struct tab_table *summary;
 849
 850   struct table_entry **pb = sorted_tab, **pe;
 851   int pc = n_sorted_tab;
 852   int cur_tab = 0;
 853
 854   summary = tab_create (7, 3 + nxtab, 1);
 855   tab_title (summary, _("Summary."));
 856   tab_headers (summary, 1, 0, 3, 0);
 857   tab_joint_text (summary, 1, 0, 6, 0, TAB_CENTER, _("Cases"));
 858   tab_joint_text (summary, 1, 1, 2, 1, TAB_CENTER, _("Valid"));
 859   tab_joint_text (summary, 3, 1, 4, 1, TAB_CENTER, _("Missing"));
 860   tab_joint_text (summary, 5, 1, 6, 1, TAB_CENTER, _("Total"));
 861   tab_hline (summary, TAL_1, 1, 6, 1);
 862   tab_hline (summary, TAL_1, 1, 6, 2);
 863   tab_vline (summary, TAL_1, 3, 1, 1);
 864   tab_vline (summary, TAL_1, 5, 1, 1);
 865   {
 866     int i;
 867
 868     for (i = 0; i < 3; i++)
 869       {
 870         tab_text (summary, 1 + i * 2, 2, TAB_RIGHT, _("N"));
 871         tab_text (summary, 2 + i * 2, 2, TAB_RIGHT, _("Percent"));
 872       }
 873   }
 874   tab_offset (summary, 0, 3);
 875
 876   for (;;)
 877     {
 878       double valid;
 879
 880       pe = find_pivot_extent (pb, &pc, cmd.pivot == CRS_PIVOT);
 881       if (pe == NULL)
 882         break;
 883
 884       while (cur_tab < (*pb)->table)
 885         insert_summary (summary, cur_tab++, 0.);
 886
 887       if (mode == GENERAL)
 888         for (valid = 0.; pb < pe; pb++)
 889           valid += (*pb)->u.freq;
 890       else
 891         {
 892           const struct crosstab *const x = xtab[(*pb)->table];
 893           const int n_cols = get_var_range (x->vars[COL_VAR])->count;
 894           const int n_rows = get_var_range (x->vars[ROW_VAR])->count;
 895           const int count = n_cols * n_rows;
 896
 897           for (valid = 0.; pb < pe; pb++)
 898             {
 899               const double *data = (*pb)->u.data;
 900               int i;
 901
 902               for (i = 0; i < count; i++)
 903                 valid += *data++;
 904             }
 905         }
 906       insert_summary (summary, cur_tab++, valid);
 907
 908       pb = pe;
 909     }
 910
 911   while (cur_tab < nxtab)
 912     insert_summary (summary, cur_tab++, 0.);
 913
 914   submit (summary);
 915 }
 916
 917 /* Inserts a line into T describing the crosstabulation at index
 918    TAB_INDEX, which has VALID valid observations. */
 919 static void
 920 insert_summary (struct tab_table *t, int tab_index, double valid)
 921 {
 922   struct crosstab *x = xtab[tab_index];
 923
 924   tab_hline (t, TAL_1, 0, 6, 0);
 925
 926   /* Crosstabulation name. */
 927   {
 928     char *buf = xmalloca (128 * x->nvar);
 929     char *cp = buf;
 930     int i;
 931
 932     for (i = 0; i < x->nvar; i++)
 933       {
 934         if (i > 0)
 935           cp = stpcpy (cp, " * ");
 936
 937         cp = stpcpy (cp, var_to_string (x->vars[i]));
 938       }
 939     tab_text (t, 0, 0, TAB_LEFT, buf);
 940
 941     freea (buf);
 942   }
 943
 944   /* Counts and percentages. */
 945   {
 946     double n[3];
 947     int i;
 948
 949     n[0] = valid;
 950     n[1] = x->missing;
 951     n[2] = n[0] + n[1];
 952
 953
 954     for (i = 0; i < 3; i++)
 955       {
 956         tab_float (t, i * 2 + 1, 0, TAB_RIGHT, n[i], 8, 0);
 957         tab_text (t, i * 2 + 2, 0, TAB_RIGHT | TAT_PRINTF, "%.1f%%",
 958                   n[i] / n[2] * 100.);
 959       }
 960   }
 961
 962   tab_next_row (t);
 963 }
 964 \f
 965 /* Output. */
 966
 967 /* Tables. */
 968 static struct tab_table *table; /* Crosstabulation table. */
 969 static struct tab_table *chisq; /* Chi-square table. */
 970 static struct tab_table *sym;           /* Symmetric measures table. */
 971 static struct tab_table *risk;          /* Risk estimate table. */
 972 static struct tab_table *direct;        /* Directional measures table. */
 973
 974 /* Statistics. */
 975 static int chisq_fisher;        /* Did any rows include Fisher's exact test? */
 976
 977 /* Column values, number of columns. */
 978 static union value *cols;
 979 static int n_cols;
 980
 981 /* Row values, number of rows. */
 982 static union value *rows;
 983 static int n_rows;
 984
 985 /* Number of statistically interesting columns/rows (columns/rows with
 986    data in them). */
 987 static int ns_cols, ns_rows;
 988
 989 /* Crosstabulation. */
 990 static const struct crosstab *x;
 991
 992 /* Number of variables from the crosstabulation to consider.  This is
 993    either x->nvar, if pivoting is on, or 2, if pivoting is off. */
 994 static int nvar;
 995
 996 /* Matrix contents. */
 997 static double *mat;             /* Matrix proper. */
 998 static double *row_tot;         /* Row totals. */
 999 static double *col_tot;         /* Column totals. */
1000 static double W;                /* Grand total. */
1001
1002 static void display_dimensions (struct tab_table *, int first_difference,
1003                                 struct table_entry *);
1004 static void display_crosstabulation (void);
1005 static void display_chisq (void);
1006 static void display_symmetric (void);
1007 static void display_risk (void);
1008 static void display_directional (void);
1009 static void crosstabs_dim (struct tab_table *, struct outp_driver *);
1010 static void table_value_missing (struct tab_table *table, int c, int r,
1011                                  unsigned char opt, const union value *v,
1012                                  const struct variable *var);
1013 static void delete_missing (void);
1014
1015 /* Output pivot table beginning at PB and continuing until PE,
1016    exclusive.  For efficiency, *MATP is a pointer to a matrix that can
1017    hold *MAXROWS entries. */
1018 static void
1019 output_pivot_table (struct table_entry **pb, struct table_entry **pe,
1020                     double **matp, double **row_totp, double **col_totp,
1021                     int *maxrows, int *maxcols, int *maxcells)
1022 {
1023   /* Subtable. */
1024   struct table_entry **tb = pb, **te;   /* Table begin, table end. */
1025   int tc = pe - pb;             /* Table count. */
1026
1027   /* Table entry for header comparison. */
1028   struct table_entry *cmp = NULL;
1029
1030   x = xtab[(*pb)->table];
1031   enum_var_values (pb, pe - pb, COL_VAR, &cols, &n_cols);
1032
1033   nvar = cmd.pivot == CRS_PIVOT ? x->nvar : 2;
1034
1035   /* Crosstabulation table initialization. */
1036   if (num_cells)
1037     {
1038       table = tab_create (nvar + n_cols,
1039                           (pe - pb) / n_cols * 3 / 2 * num_cells + 10, 1);
1040       tab_headers (table, nvar - 1, 0, 2, 0);
1041
1042       /* First header line. */
1043       tab_joint_text (table, nvar - 1, 0, (nvar - 1) + (n_cols - 1), 0,
1044                       TAB_CENTER | TAT_TITLE, var_get_name (x->vars[COL_VAR]));
1045
1046       tab_hline (table, TAL_1, nvar - 1, nvar + n_cols - 2, 1);
1047
1048       /* Second header line. */
1049       {
1050         int i;
1051
1052         for (i = 2; i < nvar; i++)
1053           tab_joint_text (table, nvar - i - 1, 0, nvar - i - 1, 1,
1054                           TAB_RIGHT | TAT_TITLE, var_to_string (x->vars[i]));
1055         tab_text (table, nvar - 2, 1, TAB_RIGHT | TAT_TITLE,
1056                   var_get_name (x->vars[ROW_VAR]));
1057         for (i = 0; i < n_cols; i++)
1058           table_value_missing (table, nvar + i - 1, 1, TAB_RIGHT, &cols[i],
1059                                x->vars[COL_VAR]);
1060         tab_text (table, nvar + n_cols - 1, 1, TAB_CENTER, _("Total"));
1061       }
1062
1063       tab_hline (table, TAL_1, 0, nvar + n_cols - 1, 2);
1064       tab_vline (table, TAL_1, nvar + n_cols - 1, 0, 1);
1065
1066       /* Title. */
1067       {
1068         char *title = xmalloca (x->nvar * 64 + 128);
1069         char *cp = title;
1070         int i;
1071
1072         if (cmd.pivot == CRS_PIVOT)
1073           for (i = 0; i < nvar; i++)
1074             {
1075               if (i)
1076                 cp = stpcpy (cp, " by ");
1077               cp = stpcpy (cp, var_get_name (x->vars[i]));
1078             }
1079         else
1080           {
1081             cp = spprintf (cp, "%s by %s for",
1082                            var_get_name (x->vars[0]),
1083                            var_get_name (x->vars[1]));
1084             for (i = 2; i < nvar; i++)
1085               {
1086                 char buf[64], *bufp;
1087
1088                 if (i > 2)
1089                   *cp++ = ',';
1090                 *cp++ = ' ';
1091                 cp = stpcpy (cp, var_get_name (x->vars[i]));
1092                 *cp++ = '=';
1093                 format_short (buf, var_get_print_format (x->vars[i]),
1094                               &(*pb)->values[i]);
1095                 for (bufp = buf; isspace ((unsigned char) *bufp); bufp++)
1096                   ;
1097                 cp = stpcpy (cp, bufp);
1098               }
1099           }
1100
1101         cp = stpcpy (cp, " [");
1102         for (i = 0; i < num_cells; i++)
1103           {
1104             struct tuple
1105               {
1106                 int value;
1107                 const char *name;
1108               };
1109
1110             static const struct tuple cell_names[] =
1111               {
1112                 {CRS_CL_COUNT, N_("count")},
1113                 {CRS_CL_ROW, N_("row %")},
1114                 {CRS_CL_COLUMN, N_("column %")},
1115                 {CRS_CL_TOTAL, N_("total %")},
1116                 {CRS_CL_EXPECTED, N_("expected")},
1117                 {CRS_CL_RESIDUAL, N_("residual")},
1118                 {CRS_CL_SRESIDUAL, N_("std. resid.")},
1119                 {CRS_CL_ASRESIDUAL, N_("adj. resid.")},
1120                 {-1, NULL},
1121               };
1122
1123             const struct tuple *t;
1124
1125             for (t = cell_names; t->value != cells[i]; t++)
1126               assert (t->value != -1);
1127             if (i)
1128               cp = stpcpy (cp, ", ");
1129             cp = stpcpy (cp, gettext (t->name));
1130           }
1131         strcpy (cp, "].");
1132
1133         tab_title (table, "%s", title);
1134         freea (title);
1135       }
1136
1137       tab_offset (table, 0, 2);
1138     }
1139   else
1140     table = NULL;
1141
1142   /* Chi-square table initialization. */
1143   if (cmd.a_statistics[CRS_ST_CHISQ])
1144     {
1145       chisq = tab_create (6 + (nvar - 2),
1146                           (pe - pb) / n_cols * 3 / 2 * N_CHISQ + 10, 1);
1147       tab_headers (chisq, 1 + (nvar - 2), 0, 1, 0);
1148
1149       tab_title (chisq, _("Chi-square tests."));
1150
1151       tab_offset (chisq, nvar - 2, 0);
1152       tab_text (chisq, 0, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1153       tab_text (chisq, 1, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1154       tab_text (chisq, 2, 0, TAB_RIGHT | TAT_TITLE, _("df"));
1155       tab_text (chisq, 3, 0, TAB_RIGHT | TAT_TITLE,
1156                 _("Asymp. Sig. (2-sided)"));
1157       tab_text (chisq, 4, 0, TAB_RIGHT | TAT_TITLE,
1158                 _("Exact. Sig. (2-sided)"));
1159       tab_text (chisq, 5, 0, TAB_RIGHT | TAT_TITLE,
1160                 _("Exact. Sig. (1-sided)"));
1161       chisq_fisher = 0;
1162       tab_offset (chisq, 0, 1);
1163     }
1164   else
1165     chisq = NULL;
1166
1167   /* Symmetric measures. */
1168   if (cmd.a_statistics[CRS_ST_PHI] || cmd.a_statistics[CRS_ST_CC]
1169       || cmd.a_statistics[CRS_ST_BTAU] || cmd.a_statistics[CRS_ST_CTAU]
1170       || cmd.a_statistics[CRS_ST_GAMMA] || cmd.a_statistics[CRS_ST_CORR]
1171       || cmd.a_statistics[CRS_ST_KAPPA])
1172     {
1173       sym = tab_create (6 + (nvar - 2), (pe - pb) / n_cols * 7 + 10, 1);
1174       tab_headers (sym, 2 + (nvar - 2), 0, 1, 0);
1175       tab_title (sym, _("Symmetric measures."));
1176
1177       tab_offset (sym, nvar - 2, 0);
1178       tab_text (sym, 0, 0, TAB_LEFT | TAT_TITLE, _("Category"));
1179       tab_text (sym, 1, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1180       tab_text (sym, 2, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1181       tab_text (sym, 3, 0, TAB_RIGHT | TAT_TITLE, _("Asymp. Std. Error"));
1182       tab_text (sym, 4, 0, TAB_RIGHT | TAT_TITLE, _("Approx. T"));
1183       tab_text (sym, 5, 0, TAB_RIGHT | TAT_TITLE, _("Approx. Sig."));
1184       tab_offset (sym, 0, 1);
1185     }
1186   else
1187     sym = NULL;
1188
1189   /* Risk estimate. */
1190   if (cmd.a_statistics[CRS_ST_RISK])
1191     {
1192       risk = tab_create (4 + (nvar - 2), (pe - pb) / n_cols * 4 + 10, 1);
1193       tab_headers (risk, 1 + nvar - 2, 0, 2, 0);
1194       tab_title (risk, _("Risk estimate."));
1195
1196       tab_offset (risk, nvar - 2, 0);
1197       tab_joint_text (risk, 2, 0, 3, 0, TAB_CENTER | TAT_TITLE | TAT_PRINTF,
1198                       _("95%% Confidence Interval"));
1199       tab_text (risk, 0, 1, TAB_LEFT | TAT_TITLE, _("Statistic"));
1200       tab_text (risk, 1, 1, TAB_RIGHT | TAT_TITLE, _("Value"));
1201       tab_text (risk, 2, 1, TAB_RIGHT | TAT_TITLE, _("Lower"));
1202       tab_text (risk, 3, 1, TAB_RIGHT | TAT_TITLE, _("Upper"));
1203       tab_hline (risk, TAL_1, 2, 3, 1);
1204       tab_vline (risk, TAL_1, 2, 0, 1);
1205       tab_offset (risk, 0, 2);
1206     }
1207   else
1208     risk = NULL;
1209
1210   /* Directional measures. */
1211   if (cmd.a_statistics[CRS_ST_LAMBDA] || cmd.a_statistics[CRS_ST_UC]
1212       || cmd.a_statistics[CRS_ST_D] || cmd.a_statistics[CRS_ST_ETA])
1213     {
1214       direct = tab_create (7 + (nvar - 2), (pe - pb) / n_cols * 7 + 10, 1);
1215       tab_headers (direct, 3 + (nvar - 2), 0, 1, 0);
1216       tab_title (direct, _("Directional measures."));
1217
1218       tab_offset (direct, nvar - 2, 0);
1219       tab_text (direct, 0, 0, TAB_LEFT | TAT_TITLE, _("Category"));
1220       tab_text (direct, 1, 0, TAB_LEFT | TAT_TITLE, _("Statistic"));
1221       tab_text (direct, 2, 0, TAB_LEFT | TAT_TITLE, _("Type"));
1222       tab_text (direct, 3, 0, TAB_RIGHT | TAT_TITLE, _("Value"));
1223       tab_text (direct, 4, 0, TAB_RIGHT | TAT_TITLE, _("Asymp. Std. Error"));
1224       tab_text (direct, 5, 0, TAB_RIGHT | TAT_TITLE, _("Approx. T"));
1225       tab_text (direct, 6, 0, TAB_RIGHT | TAT_TITLE, _("Approx. Sig."));
1226       tab_offset (direct, 0, 1);
1227     }
1228   else
1229     direct = NULL;
1230
1231   for (;;)
1232     {
1233       /* Find pivot subtable if applicable. */
1234       te = find_pivot_extent (tb, &tc, 0);
1235       if (te == NULL)
1236         break;
1237
1238       /* Find all the row variable values. */
1239       enum_var_values (tb, te - tb, ROW_VAR, &rows, &n_rows);
1240
1241       /* Allocate memory space for the column and row totals. */
1242       if (n_rows > *maxrows)
1243         {
1244           *row_totp = xnrealloc (*row_totp, n_rows, sizeof **row_totp);
1245           row_tot = *row_totp;
1246           *maxrows = n_rows;
1247         }
1248       if (n_cols > *maxcols)
1249         {
1250           *col_totp = xnrealloc (*col_totp, n_cols, sizeof **col_totp);
1251           col_tot = *col_totp;
1252           *maxcols = n_cols;
1253         }
1254
1255       /* Allocate table space for the matrix. */
1256       if (table && tab_row (table) + (n_rows + 1) * num_cells > tab_nr (table))
1257         tab_realloc (table, -1,
1258                      MAX (tab_nr (table) + (n_rows + 1) * num_cells,
1259                           tab_nr (table) * (pe - pb) / (te - tb)));
1260
1261       if (mode == GENERAL)
1262         {
1263           /* Allocate memory space for the matrix. */
1264           if (n_cols * n_rows > *maxcells)
1265             {
1266               *matp = xnrealloc (*matp, n_cols * n_rows, sizeof **matp);
1267               *maxcells = n_cols * n_rows;
1268             }
1269
1270           mat = *matp;
1271
1272           /* Build the matrix and calculate column totals. */
1273           {
1274             union value *cur_col = cols;
1275             union value *cur_row = rows;
1276             double *mp = mat;
1277             double *cp = col_tot;
1278             struct table_entry **p;
1279
1280             *cp = 0.;
1281             for (p = &tb[0]; p < te; p++)
1282               {
1283                 for (; memcmp (cur_col, &(*p)->values[COL_VAR], sizeof *cur_col);
1284                      cur_row = rows)
1285                   {
1286                     *++cp = 0.;
1287                     for (; cur_row < &rows[n_rows]; cur_row++)
1288                       {
1289                         *mp = 0.;
1290                         mp += n_cols;
1291                       }
1292                     cur_col++;
1293                     mp = &mat[cur_col - cols];
1294                   }
1295
1296                 for (; memcmp (cur_row, &(*p)->values[ROW_VAR], sizeof *cur_row);
1297                      cur_row++)
1298                   {
1299                     *mp = 0.;
1300                     mp += n_cols;
1301                   }
1302
1303                 *cp += *mp = (*p)->u.freq;
1304                 mp += n_cols;
1305                 cur_row++;
1306               }
1307
1308             /* Zero out the rest of the matrix. */
1309             for (; cur_row < &rows[n_rows]; cur_row++)
1310               {
1311                 *mp = 0.;
1312                 mp += n_cols;
1313               }
1314             cur_col++;
1315             if (cur_col < &cols[n_cols])
1316               {
1317                 const int rem_cols = n_cols - (cur_col - cols);
1318                 int c, r;
1319
1320                 for (c = 0; c < rem_cols; c++)
1321                   *++cp = 0.;
1322                 mp = &mat[cur_col - cols];
1323                 for (r = 0; r < n_rows; r++)
1324                   {
1325                     for (c = 0; c < rem_cols; c++)
1326                       *mp++ = 0.;
1327                     mp += n_cols - rem_cols;
1328                   }
1329               }
1330           }
1331         }
1332       else
1333         {
1334           int r, c;
1335           double *tp = col_tot;
1336
1337           assert (mode == INTEGER);
1338           mat = (*tb)->u.data;
1339           ns_cols = n_cols;
1340
1341           /* Calculate column totals. */
1342           for (c = 0; c < n_cols; c++)
1343             {
1344               double cum = 0.;
1345               double *cp = &mat[c];
1346
1347               for (r = 0; r < n_rows; r++)
1348                 cum += cp[r * n_cols];
1349               *tp++ = cum;
1350             }
1351         }
1352
1353       {
1354         double *cp;
1355
1356         for (ns_cols = 0, cp = col_tot; cp < &col_tot[n_cols]; cp++)
1357           ns_cols += *cp != 0.;
1358       }
1359
1360       /* Calculate row totals. */
1361       {
1362         double *mp = mat;
1363         double *rp = row_tot;
1364         int r, c;
1365
1366         for (ns_rows = 0, r = n_rows; r--; )
1367           {
1368             double cum = 0.;
1369             for (c = n_cols; c--; )
1370               cum += *mp++;
1371             *rp++ = cum;
1372             if (cum != 0.)
1373               ns_rows++;
1374           }
1375       }
1376
1377       /* Calculate grand total. */
1378       {
1379         double *tp;
1380         double cum = 0.;
1381         int n;
1382
1383         if (n_rows < n_cols)
1384           tp = row_tot, n = n_rows;
1385         else
1386           tp = col_tot, n = n_cols;
1387         while (n--)
1388           cum += *tp++;
1389         W = cum;
1390       }
1391
1392       /* Find the first variable that differs from the last subtable,
1393          then display the values of the dimensioning variables for
1394          each table that needs it. */
1395       {
1396         int first_difference = nvar - 1;
1397
1398         if (tb != pb)
1399           for (; ; first_difference--)
1400             {
1401               assert (first_difference >= 2);
1402               if (memcmp (&cmp->values[first_difference],
1403                           &(*tb)->values[first_difference],
1404                           sizeof *cmp->values))
1405                 break;
1406             }
1407         cmp = *tb;
1408
1409         if (table)
1410           display_dimensions (table, first_difference, *tb);
1411         if (chisq)
1412           display_dimensions (chisq, first_difference, *tb);
1413         if (sym)
1414           display_dimensions (sym, first_difference, *tb);
1415         if (risk)
1416           display_dimensions (risk, first_difference, *tb);
1417         if (direct)
1418           display_dimensions (direct, first_difference, *tb);
1419       }
1420
1421       if (table)
1422         display_crosstabulation ();
1423       if (cmd.miss == CRS_REPORT)
1424         delete_missing ();
1425       if (chisq)
1426         display_chisq ();
1427       if (sym)
1428         display_symmetric ();
1429       if (risk)
1430         display_risk ();
1431       if (direct)
1432         display_directional ();
1433
1434       tb = te;
1435       free (rows);
1436     }
1437
1438   submit (table);
1439
1440   if (chisq)
1441     {
1442       if (!chisq_fisher)
1443         tab_resize (chisq, 4 + (nvar - 2), -1);
1444       submit (chisq);
1445     }
1446
1447   submit (sym);
1448   submit (risk);
1449   submit (direct);
1450
1451   free (cols);
1452 }
1453
1454 /* Delete missing rows and columns for statistical analysis when
1455    /MISSING=REPORT. */
1456 static void
1457 delete_missing (void)
1458 {
1459   {
1460     int r;
1461
1462     for (r = 0; r < n_rows; r++)
1463       if (var_is_num_missing (x->vars[ROW_VAR], rows[r].f, MV_USER))
1464         {
1465           int c;
1466
1467           for (c = 0; c < n_cols; c++)
1468             mat[c + r * n_cols] = 0.;
1469           ns_rows--;
1470         }
1471   }
1472
1473   {
1474     int c;
1475
1476     for (c = 0; c < n_cols; c++)
1477       if (var_is_num_missing (x->vars[COL_VAR], cols[c].f, MV_USER))
1478         {
1479           int r;
1480
1481           for (r = 0; r < n_rows; r++)
1482             mat[c + r * n_cols] = 0.;
1483           ns_cols--;
1484         }
1485   }
1486 }
1487
1488 /* Prepare table T for submission, and submit it. */
1489 static void
1490 submit (struct tab_table *t)
1491 {
1492   int i;
1493
1494   if (t == NULL)
1495     return;
1496
1497   tab_resize (t, -1, 0);
1498   if (tab_nr (t) == tab_t (t))
1499     {
1500       tab_destroy (t);
1501       return;
1502     }
1503   tab_offset (t, 0, 0);
1504   if (t != table)
1505     for (i = 2; i < nvar; i++)
1506       tab_text (t, nvar - i - 1, 0, TAB_RIGHT | TAT_TITLE,
1507                 var_to_string (x->vars[i]));
1508   tab_box (t, TAL_2, TAL_2, -1, -1, 0, 0, tab_nc (t) - 1, tab_nr (t) - 1);
1509   tab_box (t, -1, -1, -1, TAL_1, tab_l (t), tab_t (t) - 1, tab_nc (t) - 1,
1510            tab_nr (t) - 1);
1511   tab_box (t, -1, -1, -1, TAL_GAP, 0, tab_t (t), tab_l (t) - 1,
1512            tab_nr (t) - 1);
1513   tab_vline (t, TAL_2, tab_l (t), 0, tab_nr (t) - 1);
1514   tab_dim (t, crosstabs_dim);
1515   tab_submit (t);
1516 }
1517
1518 /* Sets the widths of all the columns and heights of all the rows in
1519    table T for driver D. */
1520 static void
1521 crosstabs_dim (struct tab_table *t, struct outp_driver *d)
1522 {
1523   int i;
1524
1525   /* Width of a numerical column. */
1526   int c = outp_string_width (d, "0.000000", OUTP_PROPORTIONAL);
1527   if (cmd.miss == CRS_REPORT)
1528     c += outp_string_width (d, "M", OUTP_PROPORTIONAL);
1529
1530   /* Set width for header columns. */
1531   if (t->l != 0)
1532     {
1533       size_t i;
1534       int w;
1535
1536       w = d->width - c * (t->nc - t->l);
1537       for (i = 0; i <= t->nc; i++)
1538         w -= t->wrv[i];
1539       w /= t->l;
1540
1541       if (w < d->prop_em_width * 8)
1542         w = d->prop_em_width * 8;
1543
1544       if (w > d->prop_em_width * 15)
1545         w = d->prop_em_width * 15;
1546
1547       for (i = 0; i < t->l; i++)
1548         t->w[i] = w;
1549     }
1550
1551   for (i = t->l; i < t->nc; i++)
1552     t->w[i] = c;
1553
1554   for (i = 0; i < t->nr; i++)
1555     t->h[i] = tab_natural_height (t, d, i);
1556 }
1557
1558 static struct table_entry **find_pivot_extent_general (struct table_entry **tp,
1559                                                 int *cnt, int pivot);
1560 static struct table_entry **find_pivot_extent_integer (struct table_entry **tp,
1561                                                 int *cnt, int pivot);
1562
1563 /* Calls find_pivot_extent_general or find_pivot_extent_integer, as
1564    appropriate. */
1565 static struct table_entry **
1566 find_pivot_extent (struct table_entry **tp, int *cnt, int pivot)
1567 {
1568   return (mode == GENERAL
1569           ? find_pivot_extent_general (tp, cnt, pivot)
1570           : find_pivot_extent_integer (tp, cnt, pivot));
1571 }
1572
1573 /* Find the extent of a region in TP that contains one table.  If
1574    PIVOT != 0 that means a set of table entries with identical table
1575    number; otherwise they also have to have the same values for every
1576    dimension after the row and column dimensions.  The table that is
1577    searched starts at TP and has length CNT.  Returns the first entry
1578    after the last one in the table; sets *CNT to the number of
1579    remaining values.  If there are no entries in TP at all, returns
1580    NULL.  A yucky interface, admittedly, but it works. */
1581 static struct table_entry **
1582 find_pivot_extent_general (struct table_entry **tp, int *cnt, int pivot)
1583 {
1584   struct table_entry *fp = *tp;
1585   struct crosstab *x;
1586
1587   if (*cnt == 0)
1588     return NULL;
1589   x = xtab[(*tp)->table];
1590   for (;;)
1591     {
1592       tp++;
1593       if (--*cnt == 0)
1594         break;
1595       assert (*cnt > 0);
1596
1597       if ((*tp)->table != fp->table)
1598         break;
1599       if (pivot)
1600         continue;
1601
1602       if (memcmp (&(*tp)->values[2], &fp->values[2], sizeof (union value) * (x->nvar - 2)))
1603         break;
1604     }
1605
1606   return tp;
1607 }
1608
1609 /* Integer mode correspondent to find_pivot_extent_general().  This
1610    could be optimized somewhat, but I just don't give a crap about
1611    CROSSTABS performance in integer mode, which is just a
1612    CROSSTABS wart as far as I'm concerned.
1613
1614    That said, feel free to send optimization patches to me. */
1615 static struct table_entry **
1616 find_pivot_extent_integer (struct table_entry **tp, int *cnt, int pivot)
1617 {
1618   struct table_entry *fp = *tp;
1619   struct crosstab *x;
1620
1621   if (*cnt == 0)
1622     return NULL;
1623   x = xtab[(*tp)->table];
1624   for (;;)
1625     {
1626       tp++;
1627       if (--*cnt == 0)
1628         break;
1629       assert (*cnt > 0);
1630
1631       if ((*tp)->table != fp->table)
1632         break;
1633       if (pivot)
1634         continue;
1635
1636       if (memcmp (&(*tp)->values[2], &fp->values[2],
1637                   sizeof (union value) * (x->nvar - 2)))
1638         break;
1639     }
1640
1641   return tp;
1642 }
1643
1644 /* Compares `union value's A_ and B_ and returns a strcmp()-like
1645    result.  WIDTH_ points to an int which is either 0 for a
1646    numeric value or a string width for a string value. */
1647 static int
1648 compare_value (const void *a_, const void *b_, const void *width_)
1649 {
1650   const union value *a = a_;
1651   const union value *b = b_;
1652   const int *pwidth = width_;
1653   const int width = *pwidth;
1654
1655   if (width == 0)
1656     return (a->f < b->f) ? -1 : (a->f > b->f);
1657   else
1658     return strncmp (a->s, b->s, width);
1659 }
1660
1661 /* Given an array of ENTRY_CNT table_entry structures starting at
1662    ENTRIES, creates a sorted list of the values that the variable
1663    with index VAR_IDX takes on.  The values are returned as a
1664    malloc()'darray stored in *VALUES, with the number of values
1665    stored in *VALUE_CNT.
1666    */
1667 static void
1668 enum_var_values (struct table_entry **entries, int entry_cnt, int var_idx,
1669                  union value **values, int *value_cnt)
1670 {
1671   const struct variable *v = xtab[(*entries)->table]->vars[var_idx];
1672
1673   if (mode == GENERAL)
1674     {
1675       int width = var_get_width (v);
1676       int i;
1677
1678       *values = xnmalloc (entry_cnt, sizeof **values);
1679       for (i = 0; i < entry_cnt; i++)
1680         (*values)[i] = entries[i]->values[var_idx];
1681       *value_cnt = sort_unique (*values, entry_cnt, sizeof **values,
1682                                 compare_value, &width);
1683     }
1684   else
1685     {
1686       struct var_range *vr = get_var_range (v);
1687       int i;
1688
1689       assert (mode == INTEGER);
1690       *values = xnmalloc (vr->count, sizeof **values);
1691       for (i = 0; i < vr->count; i++)
1692         (*values)[i].f = i + vr->min;
1693       *value_cnt = vr->count;
1694     }
1695 }
1696
1697 /* Sets cell (C,R) in TABLE, with options OPT, to have a value taken
1698    from V, displayed with print format spec from variable VAR.  When
1699    in REPORT missing-value mode, missing values have an M appended. */
1700 static void
1701 table_value_missing (struct tab_table *table, int c, int r, unsigned char opt,
1702                      const union value *v, const struct variable *var)
1703 {
1704   struct substring s;
1705   const struct fmt_spec *print = var_get_print_format (var);
1706
1707   const char *label = var_lookup_value_label (var, v);
1708   if (label)
1709     {
1710       tab_text (table, c, r, TAB_LEFT, label);
1711       return;
1712     }
1713
1714   s.string = tab_alloc (table, print->w);
1715   format_short (s.string, print, v);
1716   s.length = strlen (s.string);
1717   if (cmd.miss == CRS_REPORT && var_is_num_missing (var, v->f, MV_USER))
1718     s.string[s.length++] = 'M';
1719   while (s.length && *s.string == ' ')
1720     {
1721       s.length--;
1722       s.string++;
1723     }
1724   tab_raw (table, c, r, opt, &s);
1725 }
1726
1727 /* Draws a line across TABLE at the current row to indicate the most
1728    major dimension variable with index FIRST_DIFFERENCE out of NVAR
1729    that changed, and puts the values that changed into the table.  TB
1730    and X must be the corresponding table_entry and crosstab,
1731    respectively. */
1732 static void
1733 display_dimensions (struct tab_table *table, int first_difference, struct table_entry *tb)
1734 {
1735   tab_hline (table, TAL_1, nvar - first_difference - 1, tab_nc (table) - 1, 0);
1736
1737   for (; first_difference >= 2; first_difference--)
1738     table_value_missing (table, nvar - first_difference - 1, 0,
1739                          TAB_RIGHT, &tb->values[first_difference],
1740                          x->vars[first_difference]);
1741 }
1742
1743 /* Put VALUE into cell (C,R) of TABLE, suffixed with character
1744    SUFFIX if nonzero.  If MARK_MISSING is true the entry is
1745    additionally suffixed with a letter `M'. */
1746 static void
1747 format_cell_entry (struct tab_table *table, int c, int r, double value,
1748                    char suffix, bool mark_missing)
1749 {
1750   const struct fmt_spec f = {FMT_F, 10, 1};
1751   union value v;
1752   struct substring s;
1753
1754   s.length = 10;
1755   s.string = tab_alloc (table, 16);
1756   v.f = value;
1757   data_out (&v, &f, s.string);
1758   while (*s.string == ' ')
1759     {
1760       s.length--;
1761       s.string++;
1762     }
1763   if (suffix != 0)
1764     s.string[s.length++] = suffix;
1765   if (mark_missing)
1766     s.string[s.length++] = 'M';
1767
1768   tab_raw (table, c, r, TAB_RIGHT, &s);
1769 }
1770
1771 /* Displays the crosstabulation table. */
1772 static void
1773 display_crosstabulation (void)
1774 {
1775   {
1776     int r;
1777
1778     for (r = 0; r < n_rows; r++)
1779       table_value_missing (table, nvar - 2, r * num_cells,
1780                            TAB_RIGHT, &rows[r], x->vars[ROW_VAR]);
1781   }
1782   tab_text (table, nvar - 2, n_rows * num_cells,
1783             TAB_LEFT, _("Total"));
1784
1785   /* Put in the actual cells. */
1786   {
1787     double *mp = mat;
1788     int r, c, i;
1789
1790     tab_offset (table, nvar - 1, -1);
1791     for (r = 0; r < n_rows; r++)
1792       {
1793         if (num_cells > 1)
1794           tab_hline (table, TAL_1, -1, n_cols, 0);
1795         for (c = 0; c < n_cols; c++)
1796           {
1797             bool mark_missing = false;
1798             double expected_value = row_tot[r] * col_tot[c] / W;
1799             if (cmd.miss == CRS_REPORT
1800                 && (var_is_num_missing (x->vars[COL_VAR], cols[c].f, MV_USER)
1801                     || var_is_num_missing (x->vars[ROW_VAR], rows[r].f,
1802                                            MV_USER)))
1803               mark_missing = true;
1804             for (i = 0; i < num_cells; i++)
1805               {
1806                 double v;
1807                 int suffix = 0;
1808
1809                 switch (cells[i])
1810                   {
1811                   case CRS_CL_COUNT:
1812                     v = *mp;
1813                     break;
1814                   case CRS_CL_ROW:
1815                     v = *mp / row_tot[r] * 100.;
1816                     suffix = '%';
1817                     break;
1818                   case CRS_CL_COLUMN:
1819                     v = *mp / col_tot[c] * 100.;
1820                     suffix = '%';
1821                     break;
1822                   case CRS_CL_TOTAL:
1823                     v = *mp / W * 100.;
1824                     suffix = '%';
1825                     break;
1826                   case CRS_CL_EXPECTED:
1827                     v = expected_value;
1828                     break;
1829                   case CRS_CL_RESIDUAL:
1830                     v = *mp - expected_value;
1831                     break;
1832                   case CRS_CL_SRESIDUAL:
1833                     v = (*mp - expected_value) / sqrt (expected_value);
1834                     break;
1835                   case CRS_CL_ASRESIDUAL:
1836                     v = ((*mp - expected_value)
1837                          / sqrt (expected_value
1838                                  * (1. - row_tot[r] / W)
1839                                  * (1. - col_tot[c] / W)));
1840                     break;
1841                   default:
1842                     NOT_REACHED ();
1843                   }
1844
1845                 format_cell_entry (table, c, i, v, suffix, mark_missing);
1846               }
1847
1848             mp++;
1849           }
1850
1851         tab_offset (table, -1, tab_row (table) + num_cells);
1852       }
1853   }
1854
1855   /* Row totals. */
1856   {
1857     int r, i;
1858
1859     tab_offset (table, -1, tab_row (table) - num_cells * n_rows);
1860     for (r = 0; r < n_rows; r++)
1861       {
1862         char suffix = 0;
1863         bool mark_missing = false;
1864
1865         if (cmd.miss == CRS_REPORT
1866             && var_is_num_missing (x->vars[ROW_VAR], rows[r].f, MV_USER))
1867           mark_missing = true;
1868
1869         for (i = 0; i < num_cells; i++)
1870           {
1871             double v;
1872
1873             switch (cells[i])
1874               {
1875               case CRS_CL_COUNT:
1876                 v = row_tot[r];
1877                 break;
1878               case CRS_CL_ROW:
1879                 v = 100.;
1880                 suffix = '%';
1881                 break;
1882               case CRS_CL_COLUMN:
1883                 v = row_tot[r] / W * 100.;
1884                 suffix = '%';
1885                 break;
1886               case CRS_CL_TOTAL:
1887                 v = row_tot[r] / W * 100.;
1888                 suffix = '%';
1889                 break;
1890               case CRS_CL_EXPECTED:
1891               case CRS_CL_RESIDUAL:
1892               case CRS_CL_SRESIDUAL:
1893               case CRS_CL_ASRESIDUAL:
1894                 v = 0.;
1895                 break;
1896               default:
1897                 NOT_REACHED ();
1898               }
1899
1900             format_cell_entry (table, n_cols, 0, v, suffix, mark_missing);
1901             tab_next_row (table);
1902           }
1903       }
1904   }
1905
1906   /* Column totals, grand total. */
1907   {
1908     int c;
1909     int last_row = 0;
1910
1911     if (num_cells > 1)
1912       tab_hline (table, TAL_1, -1, n_cols, 0);
1913     for (c = 0; c <= n_cols; c++)
1914       {
1915         double ct = c < n_cols ? col_tot[c] : W;
1916         bool mark_missing = false;
1917         char suffix = 0;
1918         int i;
1919
1920         if (cmd.miss == CRS_REPORT && c < n_cols
1921             && var_is_num_missing (x->vars[COL_VAR], cols[c].f, MV_USER))
1922           mark_missing = true;
1923
1924         for (i = 0; i < num_cells; i++)
1925           {
1926             double v;
1927
1928             switch (cells[i])
1929               {
1930               case CRS_CL_COUNT:
1931                 v = ct;
1932                 suffix = '%';
1933                 break;
1934               case CRS_CL_ROW:
1935                 v = ct / W * 100.;
1936                 suffix = '%';
1937                 break;
1938               case CRS_CL_COLUMN:
1939                 v = 100.;
1940                 suffix = '%';
1941                 break;
1942               case CRS_CL_TOTAL:
1943                 v = ct / W * 100.;
1944                 suffix = '%';
1945                 break;
1946               case CRS_CL_EXPECTED:
1947               case CRS_CL_RESIDUAL:
1948               case CRS_CL_SRESIDUAL:
1949               case CRS_CL_ASRESIDUAL:
1950                 continue;
1951               default:
1952                 NOT_REACHED ();
1953               }
1954
1955             format_cell_entry (table, c, i, v, suffix, mark_missing);
1956           }
1957         last_row = i;
1958       }
1959
1960     tab_offset (table, -1, tab_row (table) + last_row);
1961   }
1962
1963   tab_offset (table, 0, -1);
1964 }
1965
1966 static void calc_r (double *X, double *Y, double *, double *, double *);
1967 static void calc_chisq (double[N_CHISQ], int[N_CHISQ], double *, double *);
1968
1969 /* Display chi-square statistics. */
1970 static void
1971 display_chisq (void)
1972 {
1973   static const char *chisq_stats[N_CHISQ] =
1974     {
1975       N_("Pearson Chi-Square"),
1976       N_("Likelihood Ratio"),
1977       N_("Fisher's Exact Test"),
1978       N_("Continuity Correction"),
1979       N_("Linear-by-Linear Association"),
1980     };
1981   double chisq_v[N_CHISQ];
1982   double fisher1, fisher2;
1983   int df[N_CHISQ];
1984   int s = 0;
1985
1986   int i;
1987
1988   calc_chisq (chisq_v, df, &fisher1, &fisher2);
1989
1990   tab_offset (chisq, nvar - 2, -1);
1991
1992   for (i = 0; i < N_CHISQ; i++)
1993     {
1994       if ((i != 2 && chisq_v[i] == SYSMIS)
1995           || (i == 2 && fisher1 == SYSMIS))
1996         continue;
1997       s = 1;
1998
1999       tab_text (chisq, 0, 0, TAB_LEFT, gettext (chisq_stats[i]));
2000       if (i != 2)
2001         {
2002           tab_float (chisq, 1, 0, TAB_RIGHT, chisq_v[i], 8, 3);
2003           tab_float (chisq, 2, 0, TAB_RIGHT, df[i], 8, 0);
2004           tab_float (chisq, 3, 0, TAB_RIGHT,
2005                      gsl_cdf_chisq_Q (chisq_v[i], df[i]), 8, 3);
2006         }
2007       else
2008         {
2009           chisq_fisher = 1;
2010           tab_float (chisq, 4, 0, TAB_RIGHT, fisher2, 8, 3);
2011           tab_float (chisq, 5, 0, TAB_RIGHT, fisher1, 8, 3);
2012         }
2013       tab_next_row (chisq);
2014     }
2015
2016   tab_text (chisq, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2017   tab_float (chisq, 1, 0, TAB_RIGHT, W, 8, 0);
2018   tab_next_row (chisq);
2019
2020   tab_offset (chisq, 0, -1);
2021 }
2022
2023 static int calc_symmetric (double[N_SYMMETRIC], double[N_SYMMETRIC],
2024                            double[N_SYMMETRIC]);
2025
2026 /* Display symmetric measures. */
2027 static void
2028 display_symmetric (void)
2029 {
2030   static const char *categories[] =
2031     {
2032       N_("Nominal by Nominal"),
2033       N_("Ordinal by Ordinal"),
2034       N_("Interval by Interval"),
2035       N_("Measure of Agreement"),
2036     };
2037
2038   static const char *stats[N_SYMMETRIC] =
2039     {
2040       N_("Phi"),
2041       N_("Cramer's V"),
2042       N_("Contingency Coefficient"),
2043       N_("Kendall's tau-b"),
2044       N_("Kendall's tau-c"),
2045       N_("Gamma"),
2046       N_("Spearman Correlation"),
2047       N_("Pearson's R"),
2048       N_("Kappa"),
2049     };
2050
2051   static const int stats_categories[N_SYMMETRIC] =
2052     {
2053       0, 0, 0, 1, 1, 1, 1, 2, 3,
2054     };
2055
2056   int last_cat = -1;
2057   double sym_v[N_SYMMETRIC], sym_ase[N_SYMMETRIC], sym_t[N_SYMMETRIC];
2058   int i;
2059
2060   if (!calc_symmetric (sym_v, sym_ase, sym_t))
2061     return;
2062
2063   tab_offset (sym, nvar - 2, -1);
2064
2065   for (i = 0; i < N_SYMMETRIC; i++)
2066     {
2067       if (sym_v[i] == SYSMIS)
2068         continue;
2069
2070       if (stats_categories[i] != last_cat)
2071         {
2072           last_cat = stats_categories[i];
2073           tab_text (sym, 0, 0, TAB_LEFT, gettext (categories[last_cat]));
2074         }
2075
2076       tab_text (sym, 1, 0, TAB_LEFT, gettext (stats[i]));
2077       tab_float (sym, 2, 0, TAB_RIGHT, sym_v[i], 8, 3);
2078       if (sym_ase[i] != SYSMIS)
2079         tab_float (sym, 3, 0, TAB_RIGHT, sym_ase[i], 8, 3);
2080       if (sym_t[i] != SYSMIS)
2081         tab_float (sym, 4, 0, TAB_RIGHT, sym_t[i], 8, 3);
2082       /*tab_float (sym, 5, 0, TAB_RIGHT, normal_sig (sym_v[i]), 8, 3);*/
2083       tab_next_row (sym);
2084     }
2085
2086   tab_text (sym, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2087   tab_float (sym, 2, 0, TAB_RIGHT, W, 8, 0);
2088   tab_next_row (sym);
2089
2090   tab_offset (sym, 0, -1);
2091 }
2092
2093 static int calc_risk (double[], double[], double[], union value *);
2094
2095 /* Display risk estimate. */
2096 static void
2097 display_risk (void)
2098 {
2099   char buf[256];
2100   double risk_v[3], lower[3], upper[3];
2101   union value c[2];
2102   int i;
2103
2104   if (!calc_risk (risk_v, upper, lower, c))
2105     return;
2106
2107   tab_offset (risk, nvar - 2, -1);
2108
2109   for (i = 0; i < 3; i++)
2110     {
2111       if (risk_v[i] == SYSMIS)
2112         continue;
2113
2114       switch (i)
2115         {
2116         case 0:
2117           if (var_is_numeric (x->vars[COL_VAR]))
2118             sprintf (buf, _("Odds Ratio for %s (%g / %g)"),
2119                      var_get_name (x->vars[COL_VAR]), c[0].f, c[1].f);
2120           else
2121             sprintf (buf, _("Odds Ratio for %s (%.*s / %.*s)"),
2122                      var_get_name (x->vars[COL_VAR]),
2123                      var_get_width (x->vars[COL_VAR]), c[0].s,
2124                      var_get_width (x->vars[COL_VAR]), c[1].s);
2125           break;
2126         case 1:
2127         case 2:
2128           if (var_is_numeric (x->vars[ROW_VAR]))
2129             sprintf (buf, _("For cohort %s = %g"),
2130                      var_get_name (x->vars[ROW_VAR]), rows[i - 1].f);
2131           else
2132             sprintf (buf, _("For cohort %s = %.*s"),
2133                      var_get_name (x->vars[ROW_VAR]),
2134                      var_get_width (x->vars[ROW_VAR]), rows[i - 1].s);
2135           break;
2136         }
2137
2138       tab_text (risk, 0, 0, TAB_LEFT, buf);
2139       tab_float (risk, 1, 0, TAB_RIGHT, risk_v[i], 8, 3);
2140       tab_float (risk, 2, 0, TAB_RIGHT, lower[i], 8, 3);
2141       tab_float (risk, 3, 0, TAB_RIGHT, upper[i], 8, 3);
2142       tab_next_row (risk);
2143     }
2144
2145   tab_text (risk, 0, 0, TAB_LEFT, _("N of Valid Cases"));
2146   tab_float (risk, 1, 0, TAB_RIGHT, W, 8, 0);
2147   tab_next_row (risk);
2148
2149   tab_offset (risk, 0, -1);
2150 }
2151
2152 static int calc_directional (double[N_DIRECTIONAL], double[N_DIRECTIONAL],
2153                              double[N_DIRECTIONAL]);
2154
2155 /* Display directional measures. */
2156 static void
2157 display_directional (void)
2158 {
2159   static const char *categories[] =
2160     {
2161       N_("Nominal by Nominal"),
2162       N_("Ordinal by Ordinal"),
2163       N_("Nominal by Interval"),
2164     };
2165
2166   static const char *stats[] =
2167     {
2168       N_("Lambda"),
2169       N_("Goodman and Kruskal tau"),
2170       N_("Uncertainty Coefficient"),
2171       N_("Somers' d"),
2172       N_("Eta"),
2173     };
2174
2175   static const char *types[] =
2176     {
2177       N_("Symmetric"),
2178       N_("%s Dependent"),
2179       N_("%s Dependent"),
2180     };
2181
2182   static const int stats_categories[N_DIRECTIONAL] =
2183     {
2184       0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 2, 2,
2185     };
2186
2187   static const int stats_stats[N_DIRECTIONAL] =
2188     {
2189       0, 0, 0, 1, 1, 2, 2, 2, 3, 3, 3, 4, 4,
2190     };
2191
2192   static const int stats_types[N_DIRECTIONAL] =
2193     {
2194       0, 1, 2, 1, 2, 0, 1, 2, 0, 1, 2, 1, 2,
2195     };
2196
2197   static const int *stats_lookup[] =
2198     {
2199       stats_categories,
2200       stats_stats,
2201       stats_types,
2202     };
2203
2204   static const char **stats_names[] =
2205     {
2206       categories,
2207       stats,
2208       types,
2209     };
2210
2211   int last[3] =
2212     {
2213       -1, -1, -1,
2214     };
2215
2216   double direct_v[N_DIRECTIONAL];
2217   double direct_ase[N_DIRECTIONAL];
2218   double direct_t[N_DIRECTIONAL];
2219
2220   int i;
2221
2222   if (!calc_directional (direct_v, direct_ase, direct_t))
2223     return;
2224
2225   tab_offset (direct, nvar - 2, -1);
2226
2227   for (i = 0; i < N_DIRECTIONAL; i++)
2228     {
2229       if (direct_v[i] == SYSMIS)
2230         continue;
2231
2232       {
2233         int j;
2234
2235         for (j = 0; j < 3; j++)
2236           if (last[j] != stats_lookup[j][i])
2237             {
2238               if (j < 2)
2239                 tab_hline (direct, TAL_1, j, 6, 0);
2240
2241               for (; j < 3; j++)
2242                 {
2243                   const char *string;
2244                   int k = last[j] = stats_lookup[j][i];
2245
2246                   if (k == 0)
2247                     string = NULL;
2248                   else if (k == 1)
2249                     string = var_get_name (x->vars[0]);
2250                   else
2251                     string = var_get_name (x->vars[1]);
2252
2253                   tab_text (direct, j, 0, TAB_LEFT | TAT_PRINTF,
2254                             gettext (stats_names[j][k]), string);
2255                 }
2256             }
2257       }
2258
2259       tab_float (direct, 3, 0, TAB_RIGHT, direct_v[i], 8, 3);
2260       if (direct_ase[i] != SYSMIS)
2261         tab_float (direct, 4, 0, TAB_RIGHT, direct_ase[i], 8, 3);
2262       if (direct_t[i] != SYSMIS)
2263         tab_float (direct, 5, 0, TAB_RIGHT, direct_t[i], 8, 3);
2264       /*tab_float (direct, 6, 0, TAB_RIGHT, normal_sig (direct_v[i]), 8, 3);*/
2265       tab_next_row (direct);
2266     }
2267
2268   tab_offset (direct, 0, -1);
2269 }
2270 \f
2271 /* Statistical calculations. */
2272
2273 /* Returns the value of the gamma (factorial) function for an integer
2274    argument X. */
2275 static double
2276 gamma_int (double x)
2277 {
2278   double r = 1;
2279   int i;
2280
2281   for (i = 2; i < x; i++)
2282     r *= i;
2283   return r;
2284 }
2285
2286 /* Calculate P_r as specified in _SPSS Statistical Algorithms_,
2287    Appendix 5. */
2288 static inline double
2289 Pr (int a, int b, int c, int d)
2290 {
2291   return (gamma_int (a + b + 1.) / gamma_int (a + 1.)
2292           * gamma_int (c + d + 1.) / gamma_int (b + 1.)
2293           * gamma_int (a + c + 1.) / gamma_int (c + 1.)
2294           * gamma_int (b + d + 1.) / gamma_int (d + 1.)
2295           / gamma_int (a + b + c + d + 1.));
2296 }
2297
2298 /* Swap the contents of A and B. */
2299 static inline void
2300 swap (int *a, int *b)
2301 {
2302   int t = *a;
2303   *a = *b;
2304   *b = t;
2305 }
2306
2307 /* Calculate significance for Fisher's exact test as specified in
2308    _SPSS Statistical Algorithms_, Appendix 5. */
2309 static void
2310 calc_fisher (int a, int b, int c, int d, double *fisher1, double *fisher2)
2311 {
2312   int x;
2313
2314   if (MIN (c, d) < MIN (a, b))
2315     swap (&a, &c), swap (&b, &d);
2316   if (MIN (b, d) < MIN (a, c))
2317     swap (&a, &b), swap (&c, &d);
2318   if (b * c < a * d)
2319     {
2320       if (b < c)
2321         swap (&a, &b), swap (&c, &d);
2322       else
2323         swap (&a, &c), swap (&b, &d);
2324     }
2325
2326   *fisher1 = 0.;
2327   for (x = 0; x <= a; x++)
2328     *fisher1 += Pr (a - x, b + x, c + x, d - x);
2329
2330   *fisher2 = *fisher1;
2331   for (x = 1; x <= b; x++)
2332     *fisher2 += Pr (a + x, b - x, c - x, d + x);
2333 }
2334
2335 /* Calculates chi-squares into CHISQ.  MAT is a matrix with N_COLS
2336    columns with values COLS and N_ROWS rows with values ROWS.  Values
2337    in the matrix sum to W. */
2338 static void
2339 calc_chisq (double chisq[N_CHISQ], int df[N_CHISQ],
2340             double *fisher1, double *fisher2)
2341 {
2342   int r, c;
2343
2344   chisq[0] = chisq[1] = 0.;
2345   chisq[2] = chisq[3] = chisq[4] = SYSMIS;
2346   *fisher1 = *fisher2 = SYSMIS;
2347
2348   df[0] = df[1] = (ns_cols - 1) * (ns_rows - 1);
2349
2350   if (ns_rows <= 1 || ns_cols <= 1)
2351     {
2352       chisq[0] = chisq[1] = SYSMIS;
2353       return;
2354     }
2355
2356   for (r = 0; r < n_rows; r++)
2357     for (c = 0; c < n_cols; c++)
2358       {
2359         const double expected = row_tot[r] * col_tot[c] / W;
2360         const double freq = mat[n_cols * r + c];
2361         const double residual = freq - expected;
2362
2363         chisq[0] += residual * residual / expected;
2364         if (freq)
2365           chisq[1] += freq * log (expected / freq);
2366       }
2367
2368   if (chisq[0] == 0.)
2369     chisq[0] = SYSMIS;
2370
2371   if (chisq[1] != 0.)
2372     chisq[1] *= -2.;
2373   else
2374     chisq[1] = SYSMIS;
2375
2376   /* Calculate Yates and Fisher exact test. */
2377   if (ns_cols == 2 && ns_rows == 2)
2378     {
2379       double f11, f12, f21, f22;
2380
2381       {
2382         int nz_cols[2];
2383         int i, j;
2384
2385         for (i = j = 0; i < n_cols; i++)
2386           if (col_tot[i] != 0.)
2387             {
2388               nz_cols[j++] = i;
2389               if (j == 2)
2390                 break;
2391             }
2392
2393         assert (j == 2);
2394
2395         f11 = mat[nz_cols[0]];
2396         f12 = mat[nz_cols[1]];
2397         f21 = mat[nz_cols[0] + n_cols];
2398         f22 = mat[nz_cols[1] + n_cols];
2399       }
2400
2401       /* Yates. */
2402       {
2403         const double x = fabs (f11 * f22 - f12 * f21) - 0.5 * W;
2404
2405         if (x > 0.)
2406           chisq[3] = (W * x * x
2407                       / (f11 + f12) / (f21 + f22)
2408                       / (f11 + f21) / (f12 + f22));
2409         else
2410           chisq[3] = 0.;
2411
2412         df[3] = 1.;
2413       }
2414
2415       /* Fisher. */
2416       if (f11 < 5. || f12 < 5. || f21 < 5. || f22 < 5.)
2417         calc_fisher (f11 + .5, f12 + .5, f21 + .5, f22 + .5, fisher1, fisher2);
2418     }
2419
2420   /* Calculate Mantel-Haenszel. */
2421   if (var_is_numeric (x->vars[ROW_VAR]) && var_is_numeric (x->vars[COL_VAR]))
2422     {
2423       double r, ase_0, ase_1;
2424       calc_r ((double *) rows, (double *) cols, &r, &ase_0, &ase_1);
2425
2426       chisq[4] = (W - 1.) * r * r;
2427       df[4] = 1;
2428     }
2429 }
2430
2431 /* Calculate the value of Pearson's r.  r is stored into R, ase_1 into
2432    ASE_1, and ase_0 into ASE_0.  The row and column values must be
2433    passed in X and Y. */
2434 static void
2435 calc_r (double *X, double *Y, double *r, double *ase_0, double *ase_1)
2436 {
2437   double SX, SY, S, T;
2438   double Xbar, Ybar;
2439   double sum_XYf, sum_X2Y2f;
2440   double sum_Xr, sum_X2r;
2441   double sum_Yc, sum_Y2c;
2442   int i, j;
2443
2444   for (sum_X2Y2f = sum_XYf = 0., i = 0; i < n_rows; i++)
2445     for (j = 0; j < n_cols; j++)
2446       {
2447         double fij = mat[j + i * n_cols];
2448         double product = X[i] * Y[j];
2449         double temp = fij * product;
2450         sum_XYf += temp;
2451         sum_X2Y2f += temp * product;
2452       }
2453
2454   for (sum_Xr = sum_X2r = 0., i = 0; i < n_rows; i++)
2455     {
2456       sum_Xr += X[i] * row_tot[i];
2457       sum_X2r += X[i] * X[i] * row_tot[i];
2458     }
2459   Xbar = sum_Xr / W;
2460
2461   for (sum_Yc = sum_Y2c = 0., i = 0; i < n_cols; i++)
2462     {
2463       sum_Yc += Y[i] * col_tot[i];
2464       sum_Y2c += Y[i] * Y[i] * col_tot[i];
2465     }
2466   Ybar = sum_Yc / W;
2467
2468   S = sum_XYf - sum_Xr * sum_Yc / W;
2469   SX = sum_X2r - sum_Xr * sum_Xr / W;
2470   SY = sum_Y2c - sum_Yc * sum_Yc / W;
2471   T = sqrt (SX * SY);
2472   *r = S / T;
2473   *ase_0 = sqrt ((sum_X2Y2f - (sum_XYf * sum_XYf) / W) / (sum_X2r * sum_Y2c));
2474
2475   {
2476     double s, c, y, t;
2477
2478     for (s = c = 0., i = 0; i < n_rows; i++)
2479       for (j = 0; j < n_cols; j++)
2480         {
2481           double Xresid, Yresid;
2482           double temp;
2483
2484           Xresid = X[i] - Xbar;
2485           Yresid = Y[j] - Ybar;
2486           temp = (T * Xresid * Yresid
2487                   - ((S / (2. * T))
2488                      * (Xresid * Xresid * SY + Yresid * Yresid * SX)));
2489           y = mat[j + i * n_cols] * temp * temp - c;
2490           t = s + y;
2491           c = (t - s) - y;
2492           s = t;
2493         }
2494     *ase_1 = sqrt (s) / (T * T);
2495   }
2496 }
2497
2498 static double somers_d_v[3];
2499 static double somers_d_ase[3];
2500 static double somers_d_t[3];
2501
2502 /* Calculate symmetric statistics and their asymptotic standard
2503    errors.  Returns 0 if none could be calculated. */
2504 static int
2505 calc_symmetric (double v[N_SYMMETRIC], double ase[N_SYMMETRIC],
2506                 double t[N_SYMMETRIC])
2507 {
2508   int q = MIN (ns_rows, ns_cols);
2509
2510   if (q <= 1)
2511     return 0;
2512
2513   {
2514     int i;
2515
2516     if (v)
2517       for (i = 0; i < N_SYMMETRIC; i++)
2518         v[i] = ase[i] = t[i] = SYSMIS;
2519   }
2520
2521   /* Phi, Cramer's V, contingency coefficient. */
2522   if (cmd.a_statistics[CRS_ST_PHI] || cmd.a_statistics[CRS_ST_CC])
2523     {
2524       double Xp = 0.;   /* Pearson chi-square. */
2525
2526       {
2527         int r, c;
2528
2529         for (r = 0; r < n_rows; r++)
2530           for (c = 0; c < n_cols; c++)
2531             {
2532               const double expected = row_tot[r] * col_tot[c] / W;
2533               const double freq = mat[n_cols * r + c];
2534               const double residual = freq - expected;
2535
2536               Xp += residual * residual / expected;
2537             }
2538       }
2539
2540       if (cmd.a_statistics[CRS_ST_PHI])
2541         {
2542           v[0] = sqrt (Xp / W);
2543           v[1] = sqrt (Xp / (W * (q - 1)));
2544         }
2545       if (cmd.a_statistics[CRS_ST_CC])
2546         v[2] = sqrt (Xp / (Xp + W));
2547     }
2548
2549   if (cmd.a_statistics[CRS_ST_BTAU] || cmd.a_statistics[CRS_ST_CTAU]
2550       || cmd.a_statistics[CRS_ST_GAMMA] || cmd.a_statistics[CRS_ST_D])
2551     {
2552       double *cum;
2553       double Dr, Dc;
2554       double P, Q;
2555       double btau_cum, ctau_cum, gamma_cum, d_yx_cum, d_xy_cum;
2556       double btau_var;
2557
2558       {
2559         int r, c;
2560
2561         Dr = Dc = W * W;
2562         for (r = 0; r < n_rows; r++)
2563           Dr -= row_tot[r] * row_tot[r];
2564         for (c = 0; c < n_cols; c++)
2565           Dc -= col_tot[c] * col_tot[c];
2566       }
2567
2568       {
2569         int r, c;
2570
2571         cum = xnmalloc (n_cols * n_rows, sizeof *cum);
2572         for (c = 0; c < n_cols; c++)
2573           {
2574             double ct = 0.;
2575
2576             for (r = 0; r < n_rows; r++)
2577               cum[c + r * n_cols] = ct += mat[c + r * n_cols];
2578           }
2579       }
2580
2581       /* P and Q. */
2582       {
2583         int i, j;
2584         double Cij, Dij;
2585
2586         P = Q = 0.;
2587         for (i = 0; i < n_rows; i++)
2588           {
2589             Cij = Dij = 0.;
2590
2591             for (j = 1; j < n_cols; j++)
2592               Cij += col_tot[j] - cum[j + i * n_cols];
2593
2594             if (i > 0)
2595               for (j = 1; j < n_cols; j++)
2596                 Dij += cum[j + (i - 1) * n_cols];
2597
2598             for (j = 0;;)
2599               {
2600                 double fij = mat[j + i * n_cols];
2601                 P += fij * Cij;
2602                 Q += fij * Dij;
2603
2604                 if (++j == n_cols)
2605                   break;
2606                 assert (j < n_cols);
2607
2608                 Cij -= col_tot[j] - cum[j + i * n_cols];
2609                 Dij += col_tot[j - 1] - cum[j - 1 + i * n_cols];
2610
2611                 if (i > 0)
2612                   {
2613                     Cij += cum[j - 1 + (i - 1) * n_cols];
2614                     Dij -= cum[j + (i - 1) * n_cols];
2615                   }
2616               }
2617           }
2618       }
2619
2620       if (cmd.a_statistics[CRS_ST_BTAU])
2621         v[3] = (P - Q) / sqrt (Dr * Dc);
2622       if (cmd.a_statistics[CRS_ST_CTAU])
2623         v[4] = (q * (P - Q)) / ((W * W) * (q - 1));
2624       if (cmd.a_statistics[CRS_ST_GAMMA])
2625         v[5] = (P - Q) / (P + Q);
2626
2627       /* ASE for tau-b, tau-c, gamma.  Calculations could be
2628          eliminated here, at expense of memory.  */
2629       {
2630         int i, j;
2631         double Cij, Dij;
2632
2633         btau_cum = ctau_cum = gamma_cum = d_yx_cum = d_xy_cum = 0.;
2634         for (i = 0; i < n_rows; i++)
2635           {
2636             Cij = Dij = 0.;
2637
2638             for (j = 1; j < n_cols; j++)
2639               Cij += col_tot[j] - cum[j + i * n_cols];
2640
2641             if (i > 0)
2642               for (j = 1; j < n_cols; j++)
2643                 Dij += cum[j + (i - 1) * n_cols];
2644
2645             for (j = 0;;)
2646               {
2647                 double fij = mat[j + i * n_cols];
2648
2649                 if (cmd.a_statistics[CRS_ST_BTAU])
2650                   {
2651                     const double temp = (2. * sqrt (Dr * Dc) * (Cij - Dij)
2652                                          + v[3] * (row_tot[i] * Dc
2653                                                    + col_tot[j] * Dr));
2654                     btau_cum += fij * temp * temp;
2655                   }
2656
2657                 {
2658                   const double temp = Cij - Dij;
2659                   ctau_cum += fij * temp * temp;
2660                 }
2661
2662                 if (cmd.a_statistics[CRS_ST_GAMMA])
2663                   {
2664                     const double temp = Q * Cij - P * Dij;
2665                     gamma_cum += fij * temp * temp;
2666                   }
2667
2668                 if (cmd.a_statistics[CRS_ST_D])
2669                   {
2670                     d_yx_cum += fij * pow2 (Dr * (Cij - Dij)
2671                                             - (P - Q) * (W - row_tot[i]));
2672                     d_xy_cum += fij * pow2 (Dc * (Dij - Cij)
2673                                             - (Q - P) * (W - col_tot[j]));
2674                   }
2675
2676                 if (++j == n_cols)
2677                   break;
2678                 assert (j < n_cols);
2679
2680                 Cij -= col_tot[j] - cum[j + i * n_cols];
2681                 Dij += col_tot[j - 1] - cum[j - 1 + i * n_cols];
2682
2683                 if (i > 0)
2684                   {
2685                     Cij += cum[j - 1 + (i - 1) * n_cols];
2686                     Dij -= cum[j + (i - 1) * n_cols];
2687                   }
2688               }
2689           }
2690       }
2691
2692       btau_var = ((btau_cum
2693                    - (W * pow2 (W * (P - Q) / sqrt (Dr * Dc) * (Dr + Dc))))
2694                   / pow2 (Dr * Dc));
2695       if (cmd.a_statistics[CRS_ST_BTAU])
2696         {
2697           ase[3] = sqrt (btau_var);
2698           t[3] = v[3] / (2 * sqrt ((ctau_cum - (P - Q) * (P - Q) / W)
2699                                    / (Dr * Dc)));
2700         }
2701       if (cmd.a_statistics[CRS_ST_CTAU])
2702         {
2703           ase[4] = ((2 * q / ((q - 1) * W * W))
2704                     * sqrt (ctau_cum - (P - Q) * (P - Q) / W));
2705           t[4] = v[4] / ase[4];
2706         }
2707       if (cmd.a_statistics[CRS_ST_GAMMA])
2708         {
2709           ase[5] = ((4. / ((P + Q) * (P + Q))) * sqrt (gamma_cum));
2710           t[5] = v[5] / (2. / (P + Q)
2711                          * sqrt (ctau_cum - (P - Q) * (P - Q) / W));
2712         }
2713       if (cmd.a_statistics[CRS_ST_D])
2714         {
2715           somers_d_v[0] = (P - Q) / (.5 * (Dc + Dr));
2716           somers_d_ase[0] = 2. * btau_var / (Dr + Dc) * sqrt (Dr * Dc);
2717           somers_d_t[0] = (somers_d_v[0]
2718                            / (4 / (Dc + Dr)
2719                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2720           somers_d_v[1] = (P - Q) / Dc;
2721           somers_d_ase[1] = 2. / pow2 (Dc) * sqrt (d_xy_cum);
2722           somers_d_t[1] = (somers_d_v[1]
2723                            / (2. / Dc
2724                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2725           somers_d_v[2] = (P - Q) / Dr;
2726           somers_d_ase[2] = 2. / pow2 (Dr) * sqrt (d_yx_cum);
2727           somers_d_t[2] = (somers_d_v[2]
2728                            / (2. / Dr
2729                               * sqrt (ctau_cum - pow2 (P - Q) / W)));
2730         }
2731
2732       free (cum);
2733     }
2734
2735   /* Spearman correlation, Pearson's r. */
2736   if (cmd.a_statistics[CRS_ST_CORR])
2737     {
2738       double *R = xmalloca (sizeof *R * n_rows);
2739       double *C = xmalloca (sizeof *C * n_cols);
2740
2741       {
2742         double y, t, c = 0., s = 0.;
2743         int i = 0;
2744
2745         for (;;)
2746           {
2747             R[i] = s + (row_tot[i] + 1.) / 2.;
2748             y = row_tot[i] - c;
2749             t = s + y;
2750             c = (t - s) - y;
2751             s = t;
2752             if (++i == n_rows)
2753               break;
2754             assert (i < n_rows);
2755           }
2756       }
2757
2758       {
2759         double y, t, c = 0., s = 0.;
2760         int j = 0;
2761
2762         for (;;)
2763           {
2764             C[j] = s + (col_tot[j] + 1.) / 2;
2765             y = col_tot[j] - c;
2766             t = s + y;
2767             c = (t - s) - y;
2768             s = t;
2769             if (++j == n_cols)
2770               break;
2771             assert (j < n_cols);
2772           }
2773       }
2774
2775       calc_r (R, C, &v[6], &t[6], &ase[6]);
2776       t[6] = v[6] / t[6];
2777
2778       freea (R);
2779       freea (C);
2780
2781       calc_r ((double *) rows, (double *) cols, &v[7], &t[7], &ase[7]);
2782       t[7] = v[7] / t[7];
2783     }
2784
2785   /* Cohen's kappa. */
2786   if (cmd.a_statistics[CRS_ST_KAPPA] && ns_rows == ns_cols)
2787     {
2788       double sum_fii, sum_rici, sum_fiiri_ci, sum_fijri_ci2, sum_riciri_ci;
2789       int i, j;
2790
2791       for (sum_fii = sum_rici = sum_fiiri_ci = sum_riciri_ci = 0., i = j = 0;
2792            i < ns_rows; i++, j++)
2793         {
2794           double prod, sum;
2795
2796           while (col_tot[j] == 0.)
2797             j++;
2798
2799           prod = row_tot[i] * col_tot[j];
2800           sum = row_tot[i] + col_tot[j];
2801
2802           sum_fii += mat[j + i * n_cols];
2803           sum_rici += prod;
2804           sum_fiiri_ci += mat[j + i * n_cols] * sum;
2805           sum_riciri_ci += prod * sum;
2806         }
2807       for (sum_fijri_ci2 = 0., i = 0; i < ns_rows; i++)
2808         for (j = 0; j < ns_cols; j++)
2809           {
2810             double sum = row_tot[i] + col_tot[j];
2811             sum_fijri_ci2 += mat[j + i * n_cols] * sum * sum;
2812           }
2813
2814       v[8] = (W * sum_fii - sum_rici) / (W * W - sum_rici);
2815
2816       ase[8] = sqrt ((W * W * sum_rici
2817                       + sum_rici * sum_rici
2818                       - W * sum_riciri_ci)
2819                      / (W * (W * W - sum_rici) * (W * W - sum_rici)));
2820 #if 0
2821       t[8] = v[8] / sqrt (W * (((sum_fii * (W - sum_fii))
2822                                 / pow2 (W * W - sum_rici))
2823                                + ((2. * (W - sum_fii)
2824                                    * (2. * sum_fii * sum_rici
2825                                       - W * sum_fiiri_ci))
2826                                   / cube (W * W - sum_rici))
2827                                + (pow2 (W - sum_fii)
2828                                   * (W * sum_fijri_ci2 - 4.
2829                                      * sum_rici * sum_rici)
2830                                   / pow4 (W * W - sum_rici))));
2831 #else
2832       t[8] = v[8] / ase[8];
2833 #endif
2834     }
2835
2836   return 1;
2837 }
2838
2839 /* Calculate risk estimate. */
2840 static int
2841 calc_risk (double *value, double *upper, double *lower, union value *c)
2842 {
2843   double f11, f12, f21, f22;
2844   double v;
2845
2846   {
2847     int i;
2848
2849     for (i = 0; i < 3; i++)
2850       value[i] = upper[i] = lower[i] = SYSMIS;
2851   }
2852
2853   if (ns_rows != 2 || ns_cols != 2)
2854     return 0;
2855
2856   {
2857     int nz_cols[2];
2858     int i, j;
2859
2860     for (i = j = 0; i < n_cols; i++)
2861       if (col_tot[i] != 0.)
2862         {
2863           nz_cols[j++] = i;
2864           if (j == 2)
2865             break;
2866         }
2867
2868     assert (j == 2);
2869
2870     f11 = mat[nz_cols[0]];
2871     f12 = mat[nz_cols[1]];
2872     f21 = mat[nz_cols[0] + n_cols];
2873     f22 = mat[nz_cols[1] + n_cols];
2874
2875     c[0] = cols[nz_cols[0]];
2876     c[1] = cols[nz_cols[1]];
2877   }
2878
2879   value[0] = (f11 * f22) / (f12 * f21);
2880   v = sqrt (1. / f11 + 1. / f12 + 1. / f21 + 1. / f22);
2881   lower[0] = value[0] * exp (-1.960 * v);
2882   upper[0] = value[0] * exp (1.960 * v);
2883
2884   value[1] = (f11 * (f21 + f22)) / (f21 * (f11 + f12));
2885   v = sqrt ((f12 / (f11 * (f11 + f12)))
2886             + (f22 / (f21 * (f21 + f22))));
2887   lower[1] = value[1] * exp (-1.960 * v);
2888   upper[1] = value[1] * exp (1.960 * v);
2889
2890   value[2] = (f12 * (f21 + f22)) / (f22 * (f11 + f12));
2891   v = sqrt ((f11 / (f12 * (f11 + f12)))
2892             + (f21 / (f22 * (f21 + f22))));
2893   lower[2] = value[2] * exp (-1.960 * v);
2894   upper[2] = value[2] * exp (1.960 * v);
2895
2896   return 1;
2897 }
2898
2899 /* Calculate directional measures. */
2900 static int
2901 calc_directional (double v[N_DIRECTIONAL], double ase[N_DIRECTIONAL],
2902                   double t[N_DIRECTIONAL])
2903 {
2904   {
2905     int i;
2906
2907     for (i = 0; i < N_DIRECTIONAL; i++)
2908       v[i] = ase[i] = t[i] = SYSMIS;
2909   }
2910
2911   /* Lambda. */
2912   if (cmd.a_statistics[CRS_ST_LAMBDA])
2913     {
2914       double *fim = xnmalloc (n_rows, sizeof *fim);
2915       int *fim_index = xnmalloc (n_rows, sizeof *fim_index);
2916       double *fmj = xnmalloc (n_cols, sizeof *fmj);
2917       int *fmj_index = xnmalloc (n_cols, sizeof *fmj_index);
2918       double sum_fim, sum_fmj;
2919       double rm, cm;
2920       int rm_index, cm_index;
2921       int i, j;
2922
2923       /* Find maximum for each row and their sum. */
2924       for (sum_fim = 0., i = 0; i < n_rows; i++)
2925         {
2926           double max = mat[i * n_cols];
2927           int index = 0;
2928
2929           for (j = 1; j < n_cols; j++)
2930             if (mat[j + i * n_cols] > max)
2931               {
2932                 max = mat[j + i * n_cols];
2933                 index = j;
2934               }
2935
2936           sum_fim += fim[i] = max;
2937           fim_index[i] = index;
2938         }
2939
2940       /* Find maximum for each column. */
2941       for (sum_fmj = 0., j = 0; j < n_cols; j++)
2942         {
2943           double max = mat[j];
2944           int index = 0;
2945
2946           for (i = 1; i < n_rows; i++)
2947             if (mat[j + i * n_cols] > max)
2948               {
2949                 max = mat[j + i * n_cols];
2950                 index = i;
2951               }
2952
2953           sum_fmj += fmj[j] = max;
2954           fmj_index[j] = index;
2955         }
2956
2957       /* Find maximum row total. */
2958       rm = row_tot[0];
2959       rm_index = 0;
2960       for (i = 1; i < n_rows; i++)
2961         if (row_tot[i] > rm)
2962           {
2963             rm = row_tot[i];
2964             rm_index = i;
2965           }
2966
2967       /* Find maximum column total. */
2968       cm = col_tot[0];
2969       cm_index = 0;
2970       for (j = 1; j < n_cols; j++)
2971         if (col_tot[j] > cm)
2972           {
2973             cm = col_tot[j];
2974             cm_index = j;
2975           }
2976
2977       v[0] = (sum_fim + sum_fmj - cm - rm) / (2. * W - rm - cm);
2978       v[1] = (sum_fmj - rm) / (W - rm);
2979       v[2] = (sum_fim - cm) / (W - cm);
2980
2981       /* ASE1 for Y given X. */
2982       {
2983         double accum;
2984
2985         for (accum = 0., i = 0; i < n_rows; i++)
2986           for (j = 0; j < n_cols; j++)
2987             {
2988               const int deltaj = j == cm_index;
2989               accum += (mat[j + i * n_cols]
2990                         * pow2 ((j == fim_index[i])
2991                                - deltaj
2992                                + v[0] * deltaj));
2993             }
2994
2995         ase[2] = sqrt (accum - W * v[0]) / (W - cm);
2996       }
2997
2998       /* ASE0 for Y given X. */
2999       {
3000         double accum;
3001
3002         for (accum = 0., i = 0; i < n_rows; i++)
3003           if (cm_index != fim_index[i])
3004             accum += (mat[i * n_cols + fim_index[i]]
3005                       + mat[i * n_cols + cm_index]);
3006         t[2] = v[2] / (sqrt (accum - pow2 (sum_fim - cm) / W) / (W - cm));
3007       }
3008
3009       /* ASE1 for X given Y. */
3010       {
3011         double accum;
3012
3013         for (accum = 0., i = 0; i < n_rows; i++)
3014           for (j = 0; j < n_cols; j++)
3015             {
3016               const int deltaj = i == rm_index;
3017               accum += (mat[j + i * n_cols]
3018                         * pow2 ((i == fmj_index[j])
3019                                - deltaj
3020                                + v[0] * deltaj));
3021             }
3022
3023         ase[1] = sqrt (accum - W * v[0]) / (W - rm);
3024       }
3025
3026       /* ASE0 for X given Y. */
3027       {
3028         double accum;
3029
3030         for (accum = 0., j = 0; j < n_cols; j++)
3031           if (rm_index != fmj_index[j])
3032             accum += (mat[j + n_cols * fmj_index[j]]
3033                       + mat[j + n_cols * rm_index]);
3034         t[1] = v[1] / (sqrt (accum - pow2 (sum_fmj - rm) / W) / (W - rm));
3035       }
3036
3037       /* Symmetric ASE0 and ASE1. */
3038       {
3039         double accum0;
3040         double accum1;
3041
3042         for (accum0 = accum1 = 0., i = 0; i < n_rows; i++)
3043           for (j = 0; j < n_cols; j++)
3044             {
3045               int temp0 = (fmj_index[j] == i) + (fim_index[i] == j);
3046               int temp1 = (i == rm_index) + (j == cm_index);
3047               accum0 += mat[j + i * n_cols] * pow2 (temp0 - temp1);
3048               accum1 += (mat[j + i * n_cols]
3049                          * pow2 (temp0 + (v[0] - 1.) * temp1));
3050             }
3051         ase[0] = sqrt (accum1 - 4. * W * v[0] * v[0]) / (2. * W - rm - cm);
3052         t[0] = v[0] / (sqrt (accum0 - pow2 ((sum_fim + sum_fmj - cm - rm) / W))
3053                        / (2. * W - rm - cm));
3054       }
3055
3056       free (fim);
3057       free (fim_index);
3058       free (fmj);
3059       free (fmj_index);
3060
3061       {
3062         double sum_fij2_ri, sum_fij2_ci;
3063         double sum_ri2, sum_cj2;
3064
3065         for (sum_fij2_ri = sum_fij2_ci = 0., i = 0; i < n_rows; i++)
3066           for (j = 0; j < n_cols; j++)
3067             {
3068               double temp = pow2 (mat[j + i * n_cols]);
3069               sum_fij2_ri += temp / row_tot[i];
3070               sum_fij2_ci += temp / col_tot[j];
3071             }
3072
3073         for (sum_ri2 = 0., i = 0; i < n_rows; i++)
3074           sum_ri2 += row_tot[i] * row_tot[i];
3075
3076         for (sum_cj2 = 0., j = 0; j < n_cols; j++)
3077           sum_cj2 += col_tot[j] * col_tot[j];
3078
3079         v[3] = (W * sum_fij2_ci - sum_ri2) / (W * W - sum_ri2);
3080         v[4] = (W * sum_fij2_ri - sum_cj2) / (W * W - sum_cj2);
3081       }
3082     }
3083
3084   if (cmd.a_statistics[CRS_ST_UC])
3085     {
3086       double UX, UY, UXY, P;
3087       double ase1_yx, ase1_xy, ase1_sym;
3088       int i, j;
3089
3090       for (UX = 0., i = 0; i < n_rows; i++)
3091         if (row_tot[i] > 0.)
3092           UX -= row_tot[i] / W * log (row_tot[i] / W);
3093
3094       for (UY = 0., j = 0; j < n_cols; j++)
3095         if (col_tot[j] > 0.)
3096           UY -= col_tot[j] / W * log (col_tot[j] / W);
3097
3098       for (UXY = P = 0., i = 0; i < n_rows; i++)
3099         for (j = 0; j < n_cols; j++)
3100           {
3101             double entry = mat[j + i * n_cols];
3102
3103             if (entry <= 0.)
3104               continue;
3105
3106             P += entry * pow2 (log (col_tot[j] * row_tot[i] / (W * entry)));
3107             UXY -= entry / W * log (entry / W);
3108           }
3109
3110       for (ase1_yx = ase1_xy = ase1_sym = 0., i = 0; i < n_rows; i++)
3111         for (j = 0; j < n_cols; j++)
3112           {
3113             double entry = mat[j + i * n_cols];
3114
3115             if (entry <= 0.)
3116               continue;
3117
3118             ase1_yx += entry * pow2 (UY * log (entry / row_tot[i])
3119                                     + (UX - UXY) * log (col_tot[j] / W));
3120             ase1_xy += entry * pow2 (UX * log (entry / col_tot[j])
3121                                     + (UY - UXY) * log (row_tot[i] / W));
3122             ase1_sym += entry * pow2 ((UXY
3123                                       * log (row_tot[i] * col_tot[j] / (W * W)))
3124                                      - (UX + UY) * log (entry / W));
3125           }
3126
3127       v[5] = 2. * ((UX + UY - UXY) / (UX + UY));
3128       ase[5] = (2. / (W * pow2 (UX + UY))) * sqrt (ase1_sym);
3129       t[5] = v[5] / ((2. / (W * (UX + UY)))
3130                      * sqrt (P - pow2 (UX + UY - UXY) / W));
3131
3132       v[6] = (UX + UY - UXY) / UX;
3133       ase[6] = sqrt (ase1_xy) / (W * UX * UX);
3134       t[6] = v[6] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UX));
3135
3136       v[7] = (UX + UY - UXY) / UY;
3137       ase[7] = sqrt (ase1_yx) / (W * UY * UY);
3138       t[7] = v[7] / (sqrt (P - W * pow2 (UX + UY - UXY)) / (W * UY));
3139     }
3140
3141   /* Somers' D. */
3142   if (cmd.a_statistics[CRS_ST_D])
3143     {
3144       int i;
3145
3146       if (!sym)
3147         calc_symmetric (NULL, NULL, NULL);
3148       for (i = 0; i < 3; i++)
3149         {
3150           v[8 + i] = somers_d_v[i];
3151           ase[8 + i] = somers_d_ase[i];
3152           t[8 + i] = somers_d_t[i];
3153         }
3154     }
3155
3156   /* Eta. */
3157   if (cmd.a_statistics[CRS_ST_ETA])
3158     {
3159       {
3160         double sum_Xr, sum_X2r;
3161         double SX, SXW;
3162         int i, j;
3163
3164         for (sum_Xr = sum_X2r = 0., i = 0; i < n_rows; i++)
3165           {
3166             sum_Xr += rows[i].f * row_tot[i];
3167             sum_X2r += rows[i].f * rows[i].f * row_tot[i];
3168           }
3169         SX = sum_X2r - sum_Xr * sum_Xr / W;
3170
3171         for (SXW = 0., j = 0; j < n_cols; j++)
3172           {
3173             double cum;
3174
3175             for (cum = 0., i = 0; i < n_rows; i++)
3176               {
3177                 SXW += rows[i].f * rows[i].f * mat[j + i * n_cols];
3178                 cum += rows[i].f * mat[j + i * n_cols];
3179               }
3180
3181             SXW -= cum * cum / col_tot[j];
3182           }
3183         v[11] = sqrt (1. - SXW / SX);
3184       }
3185
3186       {
3187         double sum_Yc, sum_Y2c;
3188         double SY, SYW;
3189         int i, j;
3190
3191         for (sum_Yc = sum_Y2c = 0., i = 0; i < n_cols; i++)
3192           {
3193             sum_Yc += cols[i].f * col_tot[i];
3194             sum_Y2c += cols[i].f * cols[i].f * col_tot[i];
3195           }
3196         SY = sum_Y2c - sum_Yc * sum_Yc / W;
3197
3198         for (SYW = 0., i = 0; i < n_rows; i++)
3199           {
3200             double cum;
3201
3202             for (cum = 0., j = 0; j < n_cols; j++)
3203               {
3204                 SYW += cols[j].f * cols[j].f * mat[j + i * n_cols];
3205                 cum += cols[j].f * mat[j + i * n_cols];
3206               }
3207
3208             SYW -= cum * cum / row_tot[i];
3209           }
3210         v[12] = sqrt (1. - SYW / SY);
3211       }
3212     }
3213
3214   return 1;
3215 }
3216
3217 /* A wrapper around data_out() that limits string output to short
3218    string width and null terminates the result. */
3219 static void
3220 format_short (char *s, const struct fmt_spec *fp, const union value *v)
3221 {
3222   struct fmt_spec fmt_subst;
3223
3224   /* Limit to short string width. */
3225   if (fmt_is_string (fp->type))
3226     {
3227       fmt_subst = *fp;
3228
3229       assert (fmt_subst.type == FMT_A || fmt_subst.type == FMT_AHEX);
3230       if (fmt_subst.type == FMT_A)
3231         fmt_subst.w = MIN (8, fmt_subst.w);
3232       else
3233         fmt_subst.w = MIN (16, fmt_subst.w);
3234
3235       fp = &fmt_subst;
3236     }
3237
3238   /* Format. */
3239   data_out (v, fp, s);
3240
3241   /* Null terminate. */
3242   s[fp->w] = '\0';
3243 }
3244
3245 /*
3246    Local Variables:
3247    mode: c
3248    End:
3249 */