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