src/crosstabs.q

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