src/language/stats/crosstabs.q

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