src/crosstabs.q

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