src/language/stats/crosstabs.q

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