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