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