src/crosstabs.q

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