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