src/language/stats/crosstabs.q

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