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