src/language/stats/crosstabs.q

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