src/language/stats/crosstabs.q

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