src/language/stats/crosstabs.q

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