src/crosstabs.q

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