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