src/language/stats/crosstabs.q

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