src/language/stats/crosstabs.q

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